Fast solve of linear problems

lorenzo.iron · January 30, 2025, 10:47am

Dear GetFEM users,

I would like to fast up as much as possible the execution of a trivial heat conduction problem where the external heat flux is the only changing parameter with time.
In principle, the tangent matrix remains the same whereas the RHS changes at each iteration.

If I add only linear terms they will be assembled only once, which is already a considerable time saving.
On the other hand, I assume that the linsolve procedure is repeated at each iteration, meaning by this the inversion of the tangent matrix, am I right?
For example in case of linsolve_lu, the factorization is repeated each time even thou the tangent matrix is the same?
If so, is there a way (in python) to access the result of the matrix inversion so as to be able to perform just a simple “inverse_matrix*RHS” multiplication at each iteration?

Thank you.

Konstantinos.Poulios · January 30, 2025, 1:26pm

Dear Lorenzo,

Your request makes very much sense. The solution depends very much on which solver you want to use. But if we e.g. talk about mumps, currently if you use

dx = gf.linsolve_mumps(md.tangent_matrix(), md.rhs())

many times for different rhs, the linsolve_mumps command will every time do the same LU factorization before solving. It will also use the unsymmetric mumps version which is slower than the symmetric version.

We should redesign this interface to allow something like

fact = gf.factorize_mumps(md.tangent_matrix(), "sym")
dx = gf.linsolve_mumps(fact, md.rhs())
.... # call gf.linsolve_mumps as many times as necessary
delete(fact) # close connection to MUMPS

the difficulty is that we need to create a new python object for “fact” that will store all information related to an active connection to MUMPS.

The factorize_mumps method should look like this in C++

    id.job = JOB_INIT;
    id.par = 1;
    id.sym = sym ? 2 : 0;
    id.comm_fortran = USE_COMM_WORLD;
    mumps_interf<T>::mumps_c(id);

    id.n = int(gmm::mat_nrows(A));
    id.nz = int(AA.irn.size());
    id.irn = &(AA.irn[0]);
    id.jcn = &(AA.jcn[0]);
    id.a = (MUMPS_T*)(&(AA.a[0]));

    id.ICNTL(1) = -1; // output stream for error messages
    id.ICNTL(2) = -1; // output stream for other messages
    id.ICNTL(3) = -1; // output stream for global information
    id.ICNTL(4) = 0;  // verbosity level
//    id.ICNTL(14) += 80; // small boost to the workspace size
    id.ICNTL(31) = 1;   // only factorization, no solution to follow
    id.ICNTL(33) = 1;   // request determinant calculation

    id.job = 4; // analysis (job=1) + factorization (job=2)
    mumps_interf<T>::mumps_c(id);
    mumps_error_check(id);

Then linsolve mumps should do:

    id.rhs = (MUMPS_T*)(&(rhs[0]));
    id.job = 3;
    mumps_interf<T>::mumps_c(id);
    bool ok = mumps_error_check(id);

And when the object is deleted:

    id.job = JOB_END;
    mumps_interf<T>::mumps_c(id);

lorenzo.iron · January 31, 2025, 2:24pm

Dear Kostas,
thank you for your very comprehensive reply.

Is this interface redesign something that is already in the implementation plans by the GetFEM development team? If it’s not the case, I could try by myself but my experience with C++ is poor for now and it will not be able to succeed in in the short term.

Thank you.

Konstantinos.Poulios · January 31, 2025, 3:01pm

it has been at least in my plans for quite some time. There are 3 python libraries that implement this kind of interface

We should replicate the same logic in the getfem scripting interface. The implementation could be done in the file interface/src/gf_linsolve.cc

But it requires to create a new “mumps context” object, initialized with something like:

  auto mumpsctx = std::make_shared<getfem::some_new_name_for_mumps_context_cpp_object>
                  (arg1, arg2, ...);
  id_type id = store_mumps_ctx_object(shp);
  workspace().set_dependence(id, (const void *)(...));
  out.pop().from_object_id(id, MUMPS_CTX_CLASS_ID);
}

The method store_mumps_ctx_object need to be implemented in interface/src/getfemint.h

and the MUMPS_CTX_CLASS_ID identifier needs to be added to the getfemint_class_id enumeration in interface/src/getfemint.h.

It is not crazy difficult but it requires to understand the logic of the scripting interface in GetFEM, which is not python specific.

lorenzo.iron · February 3, 2025, 9:14am

Thank you Kostas.

I think that for my short term testing purposes I could try out one of the 3 python libraries and input the relevant tangent matrix and RHS vector from GetFEM, just to check the performance improvement for my problem.
I’ve started giving a look to the GetFEM C++ interface to external solvers but I need more time and less work commitments to be able to sort out everything and contribute effectively.

Thank you for your advices.

Konstantinos.Poulios · February 3, 2025, 9:36am

Sounds like a good idea. It should not be difficult to convert the tangent matrix to any format that 3rd party libraries may require. Fore example in the following snippet that I use for eigenmode analysis, getfem.Spmat matrices are converted to scipy.sparse.csc_matrix:

import scipy.sparse.linalg
#import scipy.io

...
II = {}
...
II["u"] = range(*md.interval_of_variable("u").cumsum())
...

    K = gf.Spmat("copy", md.tangent_matrix(), II["u"], II["u"])
    print("K_uu matrix determinant %ge%i" % K.determinant()[0::2])
    K = scipy.sparse.csc_matrix((K.csc_val(),*(K.csc_ind()[::-1])))
    # scipy.io.savemat(f"K_{i}.mat", {'K':K}) # for exporting the matrix to matlab
    # K = scipy.io.loadmat(f"K_{i}.mat")["K"] # for importing the matrix from matlab
    M = gf.asm_generic(mim, 2, "rho(x)*Test2_u.Test_u", -1, md, "select_output", "u", "u")
    M = scipy.sparse.csc_matrix((M.csc_val(),*(M.csc_ind()[::-1])))
    vals,vecs = scipy.sparse.linalg.eigs(K,6,M,sigma=0.)
    vals = vals.real
    vecs = vecs.real.T
    print("EIGENVALUES:", " ".join(map(str,vals)))

    ...
    output = ()
    ...
    for ii,vec in enumerate(vecs):
      output += (mfu, vec, f"eigenmode_{ii}")
    mfout.export_to_vtu("filname.vtu", *output)

lorenzo.iron · February 6, 2025, 12:39pm

Dear Kostas,
I’m trying to introduce a “my_Solver” but I firstly tried to use linsolve_mumps() in order to learn the basics.

I cannot sort out what happens, since the results with the standard Solve() are different respect to my procedure with linsolve_mumps().

I’ve prepared a “minimal” example for test, the following:

import getfem as gf

linear_length = 0.1
linear_elements = 2
X = [round(x*linear_length/linear_elements, 3) for x in range(linear_elements+1)]
mesh = gf.Mesh("cartesian", X, X)

BRG_all = 1
mesh.set_region(BRG_all, mesh.outer_faces())
#print(mesh)

mft = gf.MeshFem(mesh, 1)   # scalar
mft.set_classical_fem(1)    # order 1
mim = gf.MeshIm(mesh, 3)    # degree 3
md = gf.Model("real")
md.add_fem_variable("T", mft) 

Brick_stiffness = md.add_linear_term(mim, '-25*Grad_T.Grad_Test_T')
Brick_natural_convection = md.add_linear_term(mim, '-25*(T-15)*Test_T', BRG_all)

#Brick_radiation = md.add_nonlinear_term(mim, '-0.8*5.670374419E-8*(pow(T+273.15,4)-pow(15+273.15,4))*Test_T', BRG_all)
#md.disable_bricks(Brick_radiation)

Theta_method = md.add_theta_method_for_first_order('T', theta = 1)
Brick_mass = md.add_linear_term(mim, "-7850*580*Dot_T*Test_T")
#md.brick_list()

dt=0.5
md.set_time_step(dt)

# # trial 1 -> calculate solution rate with a linear solver and then update the variables

md.set_variable('T', [900 for x in range(md.nbdof())])                 # K
md.set_variable('Previous_T', [900 for x in range(md.nbdof())])        # K
md.set_variable('Previous_Dot_T', [0 for x in range(md.nbdof())])
print("T = ",md.variable("T"))
print("Dot_T",md.variable("Dot_T"))
print("Previous_T = ",md.variable("Previous_T"))
print("Previous_Dot_T = ",md.variable("Previous_Dot_T"))

for i in range(3):
    md.assembly('build_all')
    print("\n", md.tangent_matrix())
    print("rhs = ", md.rhs())
    
    rate = gf.linsolve_mumps(md.tangent_matrix(), md.rhs()).flatten()
    print("\nrate = ", rate)
    md.to_variables(md.from_variables()+dt*rate)
    print("T = ",md.variable("T"))
    print("Dot_T = ",md.variable("Dot_T"))
    
    md.shift_variables_for_time_integration()

# # trial 2 -> use standard Solve()

md.set_variable('T', [900 for x in range(md.nbdof())])                 # K
md.set_variable('Previous_T', [900 for x in range(md.nbdof())])        # K
md.set_variable('Previous_Dot_T', [0 for x in range(md.nbdof())])
for i in range(3):
    md.assembly('build_all')
    print("\n", md.tangent_matrix())
    print("rhs = ", md.rhs())
    nit, conv = md.solve('max_iter', 3, 'lsolver', 'mumps')
    print ("nit =", nit, "conv =", conv)
    print("T = ",md.variable("T"))
    print("Dot_T = ",md.variable("Dot_T"))
    md.shift_variables_for_time_integration()

# # trial 3 -> like trial 1 but after running standard Solve()

md.set_variable('T', [900 for x in range(md.nbdof())])                 # K
md.set_variable('Previous_T', [900 for x in range(md.nbdof())])        # K
md.set_variable('Previous_Dot_T', [0 for x in range(md.nbdof())])
print("T = ",md.variable("T"))
print("Dot_T",md.variable("Dot_T"))
print("Previous_T = ",md.variable("Previous_T"))
print("Previous_Dot_T = ",md.variable("Previous_Dot_T"))

for i in range(3):
    md.assembly('build_all')
    print("\n", md.tangent_matrix())
    print("rhs = ", md.rhs())
    rate = gf.linsolve_mumps(md.tangent_matrix(), md.rhs()).flatten()
    print("\nrate = ", rate)
    md.to_variables(md.from_variables()+dt*rate)
    print("T = ",md.variable("T"))
    print("Dot_T = ",md.variable("Dot_T"))
    md.shift_variables_for_time_integration()

In the next post the description of my issues

lorenzo.iron · February 6, 2025, 12:59pm

The example is a 2D thermal FEM 0.1x0.1 square with 4 elements and 9 nodes.
There is the conductivity matrix, a convection one applied on all borders, and a mass term with theta method. There’ s also a radiation term but to begin with it is disabled (commented out). in this condition the problem is linear.
Time step = 0.5 s.
There are 3 Trials with just 3 time steps each.

Trial 1: after assembly(‘build_all’) I see a tangent_matrix with low values and and the RHS too. Since the problem is linear, both terms remain unchanged as expected.

Trial 2: before using of the standard Solve() the matrix and RHS are as per Trial 1 but after having used it the Solve() command the matrix changes to higher values and the RHS too. in this case the temperature result seems to be correct.

Trial 3: I run again in series the same code as per Trial 1 but in this case the matrix is according to Trial 2.

I cannot figure out what the Solve() procedure does and what should I do to correctly set the matrix and RHS with my oun solver.

Thank you.

lorenzo.iron · February 6, 2025, 1:10pm

if I switch on the radiation term (decommenting it) and I deactivate it by disable_bricks() I cannot figure our why in Trial 1 the RHS changes drastically… and the temperature solution too.
In Trial 2 the RHS is still different before running Solve() but afterwards it changes and the temperature result turns out to be correct.
In the following Trial 3 the the matrix is the same as in Trial 2 and the solution seems to be correct.

I don’t understand why the nonlinear Radiation brick affects the assembly even if it is deactivated and how to manage this case.

Thank you

lorenzo.iron · February 10, 2025, 11:50am

I made some progress.
Using md.matrix_term() instead of md.tangent_matrix() I could retrieve all the matrix components
C = mass matrix
K = conductivity matrix + convection matrix
and together with the RHS term from the model
RHS = md.rhs()
I was able to construct my rhs named Q
Q^n =RHS^n-K \,(T^{n-1}-dt\,(1-\theta)\,\dot T^{n-1})
and build my \theta-method for linear problems, that has given the same results as the standard md.Solve(), by solving the following quation for \dot T^n:
[C+K*dt*\theta] * \dot T^n = Q^n
and then calculating
T^n = T^{n-1}+dt*(\theta\, \dot T^n+(1-\theta)\,\dot T^{n-1})
I implemented also the following method, giving directly T as a result without needing \dot T, valid also for explicit problems.
Q^n = RHS^n + (C/dt-(1-\theta)K)T^{n-1}
[C/dt+K\,\theta] T^n = Q^n

What I don’t still understand, and I hope you can help me in this connection, is why the RHS is assembled differently if I run md.Solve() before running my own solver, as shown running the following code snippet.

import getfem as gf

linear_length = 0.1
linear_elements = 1
X = [round(x*linear_length/linear_elements, 3) for x in range(linear_elements+1)]
mesh = gf.Mesh("cartesian", X, X)
BRG_all = 1
mesh.set_region(BRG_all, mesh.outer_faces())
#print(mesh)

mft = gf.MeshFem(mesh, 1)     # scalar
mft.set_classical_fem(1)      # order 1
mim = gf.MeshIm(mesh, 3)      # degree 3
md = gf.Model("real")
md.add_fem_variable("T", mft)
#md.variable_list()

conductivity = md.add_linear_term(mim, '-25*Grad_T.Grad_Test_T')
convection = md.add_linear_term(mim, '-25*(T-15)*Test_T', BRG_all)
theta = 1
md.add_theta_method_for_first_order('T', theta)
mass = md.add_linear_term(mim, "-7850*580*Dot_T*Test_T")

dt=0.5
md.set_time_step(dt)


# implicit theta method similar to the standard Solve()

print("\n\n RUNNING MY SOLVER")

md.set_variable('T', [900 for x in range(md.nbdof())])
md.set_variable('Previous_T', [900 for x in range(md.nbdof())])
md.set_variable('Dot_T', [0 for x in range(md.nbdof())])
md.set_variable('Previous_Dot_T', [0 for x in range(md.nbdof())])

md.assembly('build_all')

C = md.matrix_term(mass, 0)
K = md.matrix_term(conductivity, 0)+md.matrix_term(convection, 0)
RHS = md.rhs()

for i in range(3):
    md.assembly('build_all')
    M=C+theta*K*dt
    Q=RHS-K*(md.variable("Previous_T")-dt*(1-theta)*md.variable("Previous_Dot_T"))
    md.set_variable('Dot_T', gf.linsolve_mumps(M, Q).flatten())
    md.set_variable("T", md.variable("Previous_T")+dt*(theta*md.variable("Dot_T")+(1-theta)*md.variable("Previous_Dot_T")))
    print("T = ",md.variable("T"))
    md.shift_variables_for_time_integration()


# implicit or explitic theta method

print("\n\n RUNNING MY SECOND VERSION OF THE SOLVER")

md.set_variable('T', [900 for x in range(md.nbdof())])
md.set_variable('Previous_T', [900 for x in range(md.nbdof())])
md.set_variable('Dot_T', [0 for x in range(md.nbdof())])
md.set_variable('Previous_Dot_T', [0 for x in range(md.nbdof())])

md.assembly('build_all')

C = md.matrix_term(mass, 0)*(1/dt)
K = md.matrix_term(conductivity, 0)+md.matrix_term(convection, 0)
RHS = md.rhs()

theta=0.5

for i in range(3):
    md.assembly('build_all')
    M=C+theta*K
    Q=RHS+(C+(theta-1)*K)*md.variable("T")
    md.set_variable('T', gf.linsolve_mumps(M, Q).flatten())
    print("T = ",md.variable("T"))


# md.Solve()

print("\n\n RUNNING STANDARD SOLVE")

md.set_variable('T', [900 for x in range(md.nbdof())])
md.set_variable('Previous_T', [900 for x in range(md.nbdof())])
md.set_variable('Dot_T', [0 for x in range(md.nbdof())])
md.set_variable('Previous_Dot_T', [0 for x in range(md.nbdof())])

for i in range(3):
    nit, conv = md.solve()
    print("T = ",md.variable("T"))
    md.shift_variables_for_time_integration()


# test after having run md.Solve()

print("\n\n RUNNING AGAIN MY SOLVER")

md.set_variable('T', [900 for x in range(md.nbdof())])
md.set_variable('Previous_T', [900 for x in range(md.nbdof())])
md.set_variable('Dot_T', [0 for x in range(md.nbdof())])
md.set_variable('Previous_Dot_T', [0 for x in range(md.nbdof())])

md.assembly('build_all')

C = md.matrix_term(mass, 0)
K = md.matrix_term(conductivity, 0)+md.matrix_term(convection, 0)
RHS = md.rhs()

for i in range(3):
    md.assembly('build_all')
    M=C+theta*K*dt
    Q=RHS-K*(md.variable("Previous_T")-dt*(1-theta)*md.variable("Previous_Dot_T"))
    md.set_variable('Dot_T', gf.linsolve_mumps(M, Q).flatten())
    md.set_variable("T", md.variable("Previous_T")+dt*(theta*md.variable("Dot_T")+(1-theta)*md.variable("Previous_Dot_T")))
    print("T = ",md.variable("T"))
    md.shift_variables_for_time_integration()

print("\n\n RUNNING AGAIN MY SECOND SOLVER")

md.set_variable('T', [900 for x in range(md.nbdof())])
md.set_variable('Previous_T', [900 for x in range(md.nbdof())])
md.set_variable('Dot_T', [0 for x in range(md.nbdof())])
md.set_variable('Previous_Dot_T', [0 for x in range(md.nbdof())])

md.assembly('build_all')

C = md.matrix_term(mass, 0)*(1/dt)
K = md.matrix_term(conductivity, 0)+md.matrix_term(convection, 0)
RHS = md.rhs()

theta=0.5

for i in range(3):
    md.assembly('build_all')
    M=C+theta*K
    Q=RHS+(C+(theta-1)*K)*md.variable("T")
    md.set_variable('T', gf.linsolve_mumps(M, Q).flatten())
    print("T = ",md.variable("T"))

It seems to me that the md.Solve() procedure replaces the initial RHS with Q^n =RHS^n-K \,(T^{n-1}-dt\,(1-\theta)\,\dot T^{n-1})
and as a consequence I’m no more able to access the original RHS after running md.Assembly().

Is there a way to reset the RHS calculation procedure or have you got any suggestion how should I proceed?
Thank you

lorenzo.iron · February 10, 2025, 1:47pm

Finally I found a way out.
My convection term is

id = md.add_linear_term(mim, '-25*(T-25)*Test_T', BRG_all)

Initially, i couldn’t use the method

md.brick_term_rhs(id,...)

because it returned an array of zeroes. I realized only this morning that it was due to “md.add_linear_term” function that returned 2, that correspond to the matrix term, whereas the corresponding source term Brick number is 1, but I could not manage it with the returned id.

So, with these tools understood I should be able step forward now.
Thank you

Konstantinos.Poulios · February 10, 2025, 7:17pm

Hi Lorenzo, good that you figured things out. I would like to warn you about some inconsistency. Linear and nonlinear problems are treated differently. For nonlinear problems rhs is the residual vector (which would be the equivalent of b-K.x if the problem was linear), for linear problems rhs is the vector b in the system K.x=b.

You can try this by running the following code before and after exchanging add_linear_term with add_nonlinear_term:

import getfem as gf
m = gf.Mesh("cartesian",[10,20,30],[110,120,130])
mf = gf.MeshFem(m,1)
mf.set_classical_fem(1)
md = gf.Model("real")
mim = gf.MeshIm(m,3)
md.add_fem_variable("u", mf)
md.add_linear_term(mim,"(u-1000)*Test_u")
#md.add_nonlinear_term(mim,"(u-1000)*Test_u")
md.tangent_matrix()
md.rhs()
md.assembly("build_all")
#md.tangent_matrix().full()
print(md.rhs())
md.set_variable("u", md.interpolation("1000",mf))
md.assembly("build_all")
print(md.rhs())

lorenzo.iron · February 11, 2025, 12:06pm

Hello Kostas,
I’ve another question.

If I add to your model an im_data with add_linear_term, the brick returns -1.

Since I’m adding several bricks, it is of big help for me to have the reference to them, both for activating/deactivating them and also for “my_Solver” logics.
The issue disappears using add_nonlinear_term.

Why?

Here below your modified snippet showing the issue:

import getfem as gf
m = gf.Mesh("cartesian",[10,20,30],[110,120,130])
mf = gf.MeshFem(m,1)
mf.set_classical_fem(1)
md = gf.Model("real")
mim = gf.MeshIm(m,3)
md.add_fem_variable("u", mf)
md.add_linear_term(mim,"(u-1000)*Test_u")
#md.add_nonlinear_term(mim,"(u-1000)*Test_u")

mimd = gf.MeshImData(mim, -1, 1)
md.add_im_data("q_im_data", mimd)
q_term_id = md.add_linear_term(mim,'q_im_data*Test_u')
#q_term_id = md.add_nonlinear_term(mim,'q_im_data*Test_u')
print("q_term_id =", q_term_id)

md.tangent_matrix()
md.rhs()
md.assembly("build_all")
#md.tangent_matrix().full()
print(md.rhs())
md.set_variable("u", md.interpolation("1000",mf))
md.assembly("build_all")
print(md.rhs())

Thank you

lorenzo.iron · February 11, 2025, 12:13pm

If i check the Bricks with

md.brick_list()

I can see an additional source term, but the numbering is not returned by the brick and I cannot manage it with easy.

lorenzo.iron · February 11, 2025, 1:30pm

Ok, for my easy task the solution is using md.add_source_term(mim, ‘…’) instead of md.add_linear_term(). This way the Brick returns a positive and sequential ID that I can refer to with confidence during my automatic routines for model adaptations.
Anyhow, I would like to understand what “-1” means when using md.add_linear_term(). In case the dependence on a fem variable is missing, wouldn’t be better to return the ID of the corresponding source term?
Thank you.

Konstantinos.Poulios · February 12, 2025, 12:52pm

I could fix this behavior like this

but there are cases that trigger both add_source_term_ and md.add_brick.

Then the function will return only the index of the last one.

For example, calling add_linear_term with (u-1000)*Test_u will add two bricks, one for -1000*Test_u and one for u*Test_u. This is done only for linear bricks, as a kind of optimization. So it will still be confusing for the user, because the returned id will be only for the term u*Test_u.

In that sense, in the current implementation the -1 result can be interpreted as if only a rhs has been added to a linear model, and this cannot be turned on/off upon demand.

lorenzo.iron · February 12, 2025, 6:56pm

As you told, in the current implementation when running add_linear_term with (u-1000)*Test_u, two bricks are generated but the method returns the ID_u=2 of the fem variable term whereas the ID_source=1 of the source term is to be inferred like ID_source=ID_u-1.
If you use add_linear_term with 1000*test_u it returns ID=-1 but in this case I cannot automatically infer the ID_source. This is my main issue.
On the other hand, if you use add_source_term it returns an integer ID. I’d have expected that also add_linear_term would return the same ID, since practically the effect on the liner system is the same.
If not possible, maybe it could return an array like (1,2) for the first case and (-1, 1) for the second case, giving access to the source term ID.
Otherwise, is there a way to cycle the list the the available bricks and recognize their nature, like if they are linear, nonlinear, source terms, activated, etc?

Konstantinos.Poulios · February 12, 2025, 10:03pm

Yes, that is exactly what the patch in my previous screenshot would fix. I will apply this patch.

Regarding changing the return type to a list of integers, no I would not do it, it is too radical change for the c++ interface.

Konstantinos.Poulios · February 13, 2025, 8:37am

It should be fixed now with this commit

Thanks for reporting this issue.

lorenzo.iron · February 16, 2025, 9:25pm

Dear Kostas,
thank you for your commitment.

I would like to show a couple of issues I found during my previous test, still related to the way GetFEM solves linear problems, that I deem deserves attention and a fix.

I was not able to figure out why md.tangent_matrix() returns different marixes when called after md.assembly('build_matrix') or after md.solve() when the model contains a theta method. I realized that GetFEM uses the following formula for the theta method when the problem is linear: \left( \frac{C}{\theta \, dt} + K \right) T^n = Q+C\left(\frac {T^{n-1}}{\theta \, dt}+\frac{1-\theta}{\theta}\dot T^{n-1}\right) where C is the mass matrix, K is the stiffness matrix and Q is the source term vector. After md.solve() the md.tangent_matrix() method returns exactly \left( \frac{C}{\theta \, dt} + K \right) but before running md.solve() and just after md.assembly('build_matrix') themd.tangent_matrix() method returns just \left( C + K \right) that is not the right tangent matrix! This way, when we want to use our solvers we are obliged to compose the matrices by hand or optionally we have to run md.solve() at least one time to retrieve the correct tangent matrix. In the following snippet you can see the issue.

import getfem as gf

linear_length = 0.01
linear_elements = 2
X = [round(x*linear_length/linear_elements, 3) for x in range(linear_elements+1)]
mesh = gf.Mesh("cartesian", X, X)
BRG_all = 1
mesh.set_region(BRG_all, mesh.outer_faces())

mft = gf.MeshFem(mesh, 1) 
mft.set_classical_fem(1)
mim = gf.MeshIm(mesh, 3)
md = gf.Model("real")
md.add_fem_variable("T", mft)

conductivity = md.add_linear_term(mim, '-25*Grad_T.Grad_Test_T')
convection = md.add_linear_term(mim, '-25*(T-15)*Test_T', BRG_all)
source_term = convection-1

theta = 0.5
md.add_theta_method_for_first_order('T', theta)
mass = md.add_linear_term(mim, "-7850*580*Dot_T*Test_T")

gf.util_trace_level(0)
gf.util_warning_level(0)

dt=0.5
md.set_time_step(dt)

md.set_variable('T', [900 for x in range(md.nbdof())])
md.set_variable('Previous_T', md.variable('T'))
md.set_variable('Dot_T', [0 for x in range(md.nbdof())])
md.set_variable('Previous_Dot_T', md.variable('Dot_T'))

md.set_time(0)

md.assembly('build_matrix')
C = md.matrix_term(mass, 0)
K = md.matrix_term(conductivity, 0)+md.matrix_term(convection, 0)
M = gf.Spmat('add',C*(1/dt/theta),K)

print("my wrong tangent matrix C+K", (C+K).full())
print("md.tangent_matrix()", md.tangent_matrix().full())

nit, conv = md.solve() 
print("my tangent matrix M", M.full())
print("md.tangent_matrix()", md.tangent_matrix().full())
md.shift_variables_for_time_integration()

It would be useful for me to be able to change a model from linear to nonlinear by activating/deactivating some linear/nonlinear bricks. As I mentioned in a previous post above, if I add to the model a nonlinear Brick and I deactivate it, the md.solve() method starts returning the md.rhs() as a residual vector as like as like it was a nonlinear model, but it is linear since the nonlinear Brick has been deactivated! GetFEM doesn’t recognize this and returns a correct RHS before running md.solve() but returns a residual afterwards. This inconsistency makes it difficult to manage the matrixes with custom solvers and obliges to generate the RHS by hand, that I believe it is a pity and not in the philosophy of GetFEM.
You can check by adding to my previous snippet the following one

radiation = md.add_nonlinear_term(mim, '-4.5363e-08*(pow(T,4)-pow(15,4))*Test_T', BRG_all)
md.disable_bricks(radiation)

Please, let me know what you think about the matter.
Thank you

Topic		Replies	Views
Problem solving transient brick with MUMPS Questions and issues	2	65	January 21, 2025
Problem during boundary condition assembly Questions and issues	133	440	June 23, 2025
Parallel mumps solver and custom solver General	9	77	July 15, 2025
Issue with demo_laplacian.py Questions and issues	59	257	January 24, 2025
Simple 2D isotropic linear elasticity problem on the unit square with a mesh of 2 triangles General	9	61	August 6, 2025

Fast solve of linear problems

Related topics