# Access to the matrix of the linear problem

**URL:** <https://fenicsproject.discourse.group/t/access-to-the-matrix-of-the-linear-problem/16456>\
**Category:** dolfinx\
**Created:** [December 5, 2024, 3:23pm UTC](https://fenicsproject.discourse.group/t/access-to-the-matrix-of-the-linear-problem/16456 "2024-12-05T15:23:39Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![domi](https://avatars.discourse-cdn.com/v4/letter/d/d9b06d/32.png) [@domi](https://fenicsproject.discourse.group/u/domi)\
**Post date:** [December 5, 2024, 3:23pm UTC](https://fenicsproject.discourse.group/t/access-to-the-matrix-of-the-linear-problem/16456/1 "2024-12-05T15:23:39Z")

</div>

Hello,

when solving a nonlinear problem (e.g. nonlinear Poisson Equation) within fenicsx, one can use the following code from the tutorial:

```auto
import ufl
import numpy as np

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import mesh, fem
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

def q(u):
     return 1 + u**2

domain = mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)
x = ufl.SpatialCoordinate(domain)
u_ufl = 1 + x[0] + 2 * x[1]
f = - ufl.div(q(u_ufl) * ufl.grad(u_ufl))

V = fem.functionspace(domain, ("Lagrange", 1))

def u_exact(x):
     return eval(str(u_ufl))

u_D = fem.Function(V)
u_D.interpolate(u_exact)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets))
uh = fem.Function(V)
v = ufl.TestFunction(V)
F = q(uh) * ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - f * v * ufl.dx

problem = NonlinearProblem(F, uh, bcs=[bc])
solver = NewtonSolver(MPI.COMM_WORLD, problem)

# direct solver
ksp = PETSc.KSP().create()
option_prefix = ksp.getOptionsPrefix()
opts = PETSc.Options()
opts['ksp_type'] = "preonly"
opts['pc_type'] = "lu"
ksp.setFromOptions()

n, converged = solver.solve(uh)
print(f"Number of interations: {n:d}")

```

I wonder how it is possible to access the linearized systems that are solved within the NewtonSolver (e.g. to visualize A and RHS b). I found that in dolfinx/fem/petsc.py there is the function

```auto
    def J(self, x: PETSc.Vec, A: PETSc.Mat) -> None:
        """Assemble the Jacobian matrix.

        Args:
            x: The vector containing the latest solution

        """
        A.zeroEntries()
        assemble_matrix_mat(A, self._a, self.bcs)
        A.assemble()

```

which is called in every iteration within the NewtonSolver. However, I suspect that directly editing code here (e.g. adding A.view()) is not the proper way to do this. I tried using KSPMonitor but this did not work. Maybe one has to overwrite the NewtonSolver class in order to “expose” the linearized systems? Any help is greatly apprechiated.

Thank you.

cheers

domi

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 5, 2024, 3:54pm UTC](https://fenicsproject.discourse.group/t/access-to-the-matrix-of-the-linear-problem/16456/2 "2024-12-05T15:54:56Z")

</div>

The Newton solver stores J under `solver.A`: [dolfinx/python/dolfinx/nls/petsc.py at main · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/nls/petsc.py#L55-L58)  
which you can access post solve to look at the elements.  
There are also various additional call-back functions you can overload:  
(set\_form), or `set_update`: [dolfinx/cpp/dolfinx/nls/NewtonSolver.h at main · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/main/cpp/dolfinx/nls/NewtonSolver.h#L114-L124)  
See for instance: [dolfinx/python/test/unit/nls/test\_newton.py at 56394e9e2fa0ccaa1847b39841615b62b6493c7a · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/56394e9e2fa0ccaa1847b39841615b62b6493c7a/python/test/unit/nls/test_newton.py#L133-L142)

---

<div class="post-metadata">

**Author:** ![domi](https://avatars.discourse-cdn.com/v4/letter/d/d9b06d/32.png) [@domi](https://fenicsproject.discourse.group/u/domi)\
**Post date:** [December 5, 2024, 7:51pm UTC](https://fenicsproject.discourse.group/t/access-to-the-matrix-of-the-linear-problem/16456/3 "2024-12-05T19:51:00Z")

</div>

Hi dokken,

since solver.A returns a matrix containing numbers I assume that it represents the Jacobian evaluated at the the last linearization step, right? If so, I fear that is not sufficient for my purpose.

Regarding the overloaded call-back functions I do not know how this would work in my case. If it is not too much work, could you show me how to implement it in my specific example, since I don’t get it from the test\_newton.py.

Thank you.

Cheers

domi

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 5, 2024, 8:06pm UTC](https://fenicsproject.discourse.group/t/access-to-the-matrix-of-the-linear-problem/16456/4 "2024-12-05T20:06:56Z")

</div>

> [@domi](#):
>
> since solver.A returns a matrix containing numbers I assume that it represents the Jacobian evaluated at the the last linearization step, right? If so, I fear that is not sufficient for my purpose.

Could you specify what you want to do with `A.view()`?

You could do something like:

```python
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)

def extended_J(x, A)->None:
     problem.J(x, A)
     A.view()

solver.setJ(extended_J, solver.A)

```

which means that you do not interfere much with the solver itself.

---

<div class="post-metadata">

**Author:** ![domi](https://avatars.discourse-cdn.com/v4/letter/d/d9b06d/32.png) [@domi](https://fenicsproject.discourse.group/u/domi)\
**Post date:** [December 5, 2024, 8:51pm UTC](https://fenicsproject.discourse.group/t/access-to-the-matrix-of-the-linear-problem/16456/5 "2024-12-05T20:51:51Z")

</div>

Hi dokken,

this works as expected, thanks a lot.  
I have a nonlinear problem with multiple coupled dependent variables which can be solved using a Newton method and a direct solver for inverting the linearized problems. However, when I try to use an iterative solver (gmres) combined with a preconditioner, the residual norm stagnates and I don’t understand the reason for this. My hope is now that by inspecting each linearized problem (and relating the entries to the individual depdendent variables) I’ll be able to find the cause for this behavior.

regards

domi