Fast solve of linear problems

General
21

Dear Lorenzo, without going into the technicalities of your question I would advice you against using the helper black-box add_theta_method_for_first_order method. Why not implementing the theta method yourself, it is super simple and much more transparent. Especially if you want to develop some custom model with turning variables on and off, you want to have full control of the implementation.

Another comment is about using a model or not. You do not need to do everything using the model object. You can assemble a stiffness matrix and rhs on the fly using gf.asm_generic(mim, 2, ...) and gf.asm_generic(mim, 1, ...) respectively. Then you use gf.linsolve_mumps(...) or some other solver for the linear system solution. You do not necessary need to use the model object always.

1 Like
22

Hello Kostas,
thank you for your suggestions.

Concerning the first one, Iā€™ve indeed implemented my theta method without using ā€œadd_theta_method_for_first_orderā€.
The sent example was aimed at showing my attempts to use the model object and his methods, hoping to be able to exploit the efficient already implemented functions of GetFEM and to avoid manipulating the assembled matrices by my self.

About your second suggestion, Ok, for my problem it makes much sense and I can consider to review my model with that new approach, since it will be easier for me than figuring out all the GetFEM internal logics. On the other hand, Iā€™ve not got to know the reasons of the possible ā€œinconsistenciesā€ that I found but I wanted just to share those findings with you for possible improvements.

Thank you for your always valuable help and support.
Have a nice day!
Lorenzo

23

I have now implemented a mumps context object which can be used either from C++ or the supported scripting languages.

For an example on how it can be used, check the unit test check_mumps_ctx.py.

import numpy as np

import getfem as gf

ctx1 = gf.MumpsContext("symmetric")
ctx2 = gf.MumpsContext("unsymmetric")
ctx1.set_ICNTL(4, 0) # silence MUMPS output
ctx2.set_ICNTL(4, 0) # silence MUMPS output

K1 = gf.Spmat('empty', 3, 3)
K2 = gf.Spmat('empty', 3, 3)
rhs = np.ones(3)

K1.add(0, 0, 1)
K1.add(1, 1, 2)
K1.add(2, 2, 2)
K1.add(1, 2, 0.2)
K1.add(2, 1, 0.2)

K2.add(0, 0, 1)
K2.add(1, 1, 2)
K2.add(2, 2, 3)
K2.add(1, 2, 0.2)
K2.add(2, 1, -0.2)

ctx1.set_matrix(K1)
ctx1.set_vector(rhs)

ctx1.analyze()
ctx1.factorize()
x1 = ctx1.solve()

...
24

Thank you Kostas for sharing this implementation.
Very nice, I will try out it soon!

In the meantime, for my urgent purposes I used the Scipy SuperLU factorization instead but Iā€™m going to refactor the code with the MumpsContext with very much pleasure.

Thank you again!
Lorenzo