Access to the matrix of the linear problem

dolfinx
1

Hello,

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

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

    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

2

The Newton solver stores J under solver.A: dolfinx/python/dolfinx/nls/petsc.py at main · FEniCS/dolfinx · GitHub
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
See for instance: dolfinx/python/test/unit/nls/test_newton.py at 56394e9e2fa0ccaa1847b39841615b62b6493c7a · FEniCS/dolfinx · GitHub

3

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

4

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

You could do something like:

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.

1 Like
5

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