# Fast solve of linear problems

**URL:** <https://getfem.discourse.group/t/fast-solve-of-linear-problems/61>\
**Category:** General\
**Created:** [January 30, 2025, 10:47am UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61 "2025-01-30T10:47:37Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [January 30, 2025, 10:47am UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/1 "2025-01-30T10:47:37Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [January 30, 2025, 1:26pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/2 "2025-01-30T13:26:48Z")

</div>

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

```python
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

```python
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++

```cpp
    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:

```auto
    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:

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

```

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [January 31, 2025, 2:24pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/3 "2025-01-31T14:24:04Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [January 31, 2025, 3:01pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/4 "2025-01-31T15:01:37Z")

</div>

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

> **[GitHub - PyMumps/pymumps: PyMUMPS: A parallel sparse direct solver](https://github.com/pymumps/pymumps)**
>
> PyMUMPS: A parallel sparse direct solver

> **[GitHub - akhmerov/python-mumps: Python interface for the MUMPS sparse linear...](https://github.com/akhmerov/python-mumps)**
>
> Python interface for the MUMPS sparse linear solver (mirror)

> **[GitHub - basnijholt/mumpy: Python bindings for the MUMPS package MUMPS: a...](https://github.com/basnijholt/mumpy)**
>
> Python bindings for the MUMPS package MUMPS: a parallel sparse direct solver

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:

```cpp
  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.

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 3, 2025, 9:14am UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/5 "2025-02-03T09:14:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [February 3, 2025, 9:36am UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/6 "2025-02-03T09:36:59Z")

</div>

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`:

```python
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)

```

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 6, 2025, 12:39pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/7 "2025-02-06T12:39:51Z")

</div>

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:

```auto
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

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 6, 2025, 12:59pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/8 "2025-02-06T12:59:47Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 6, 2025, 1:10pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/9 "2025-02-06T13:10:44Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 10, 2025, 11:50am UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/10 "2025-02-10T11:50:27Z")

</div>

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.

```auto
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

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 10, 2025, 1:47pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/11 "2025-02-10T13:47:17Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [February 10, 2025, 7:17pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/12 "2025-02-10T19:17:20Z")

</div>

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:

```auto
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())

```

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 11, 2025, 12:06pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/13 "2025-02-11T12:06:09Z")

</div>

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:

```auto
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

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 11, 2025, 12:13pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/14 "2025-02-11T12:13:16Z")

</div>

If i check the Bricks with

```auto
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.

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 11, 2025, 1:30pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/15 "2025-02-11T13:30:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [February 12, 2025, 12:52pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/16 "2025-02-12T12:52:42Z")

</div>

I could fix this behavior like this

 ![image](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/73f2fddd109d3ba45720f7e1fb1792dd23716e7e.png)  
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.

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 12, 2025, 6:56pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/17 "2025-02-12T18:56:14Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [February 12, 2025, 10:03pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/18 "2025-02-12T22:03:35Z")

</div>

> [@lorenzo.iron](#):
>
> 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.

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.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [February 13, 2025, 8:37am UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/19 "2025-02-13T08:37:37Z")

</div>

It should be fixed now with [this commit](https://git.savannah.nongnu.org/cgit/getfem.git/commit/?id=4375c20325c203fabdb535d5b28ee1d9f4e1f92e)

Thanks for reporting this issue.

---

<div class="post-metadata">

**Author:** ![lorenzo.iron](https://avatars.discourse-cdn.com/v4/letter/l/b5ac83/32.png) [@lorenzo.iron](https://getfem.discourse.group/u/lorenzo.iron)\
**Post date:** [February 16, 2025, 9:25pm UTC](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61/20 "2025-02-16T21:25:11Z")

</div>

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.

1. 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')` the`md.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.

```auto
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()

```

1. 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

```auto
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

[Next page](https://getfem.discourse.group/t/fast-solve-of-linear-problems/61.md?page=2)