# Solver configuration#

Author: Jørgen S. Dokken

In this section, we will go through how to specify what linear algebra solver we would like to use to solve our PDEs, as well as how to verify the implemenation by considering convergence rates.

$-\Delta u = f \text{ in } \Omega$
$u = u_D \text{ on } \partial \Omega.$

Using the manufactured solution $$u_D=\cos(2\pi x)\cos(2\pi y)$$, we obtain $$f=8\pi^2\cos(2\pi x)\cos(2\pi y)$$. We start by creating a generic module for evaluating the analytical solution at any point $$x$$.

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem import dirichletbc, FunctionSpace, Function, locate_dofs_topological
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities_boundary
from ufl import SpatialCoordinate, TestFunction, TrialFunction, div, dx, inner, grad

import numpy
import ufl

def u_ex(mod):
return lambda x: mod.cos(2*mod.pi*x[0])*mod.cos(2*mod.pi*x[1])


Note that the return type of u_ex is a lambda function. Thus, we can create two different lambda functions, one using numpy (which will be used for interpolation) and one using ufl (which will be used for defining the source term)

u_numpy = u_ex(numpy)
u_ufl = u_ex(ufl)


We start by using ufl to define our source term, using ufl.SpatialCoordinate as input to u_ufl.

mesh = create_unit_square(MPI.COMM_WORLD, 30, 30)
x = SpatialCoordinate(mesh)


Next, we define our linear variational problem

V = FunctionSpace(mesh, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
L = f * v * dx
u_bc = Function(V)
u_bc.interpolate(u_numpy)
facets = locate_entities_boundary(mesh, mesh.topology.dim-1, lambda x: numpy.full(x.shape[1], True))
dofs = locate_dofs_topological(V, mesh.topology.dim-1, facets)
bcs = [dirichletbc(u_bc, dofs)]


We start by solving the problem with an LU factorization, a direct solver method (similar to Gaussian elimination).

default_problem = LinearProblem(a, L, bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = default_problem.solve()


We now look at the solver process by inspecting the PETSc-solver. As the view-options in PETSc are not adjusted for notebooks (solver.view() will print output to the terminal if used in a .py file), we write the solver output to file and read it in and print the output.

lu_solver = default_problem.solver
viewer = PETSc.Viewer().createASCII("lu_output.txt")
lu_solver.view(viewer)
solver_output = open("lu_output.txt", "r")
print(line)

KSP Object: (dolfinx_solve_139960817495360) 1 MPI processes

type: preonly

maximum iterations=10000, initial guess is zero

tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.

left preconditioning

using NONE norm type for convergence test

PC Object: (dolfinx_solve_139960817495360) 1 MPI processes

type: lu

out-of-place factorization

tolerance for zero pivot 2.22045e-14

matrix ordering: nd

factor fill ratio given 5., needed 5.08301

Factored matrix follows:

Mat Object: 1 MPI processes

type: seqaij

rows=961, cols=961

package used to perform factorization: petsc

total: nonzeros=32943, allocated nonzeros=32943

not using I-node routines

linear system matrix = precond matrix:

Mat Object: (dolfinx_solve_139960817495360) 1 MPI processes

type: seqaij

rows=961, cols=961

total: nonzeros=6481, allocated nonzeros=6481

total number of mallocs used during MatSetValues calls=0

not using I-node routines


This is a very robust and simple method, and is the recommended method up to a few thousand unknowns and can be efficiently used for many 2D and smaller 3D problems. However, sparse LU decomposition quickly becomes slow, as for a $$N\times N$$-matrix the number of floating point operations scales as $$\sim (2/3)N^3$$.

For large problems, we instead need to use an iterative method which are faster and require less memory.

## Choosing a linear solver and preconditioner#

As the Poisson equation results in a symmetric, positive definite system matrix, the optimal Krylov solver is the conjugate gradient (CG) method. The default preconditioner is the incomplete LU factorization (ILU), which is a popular and robous overall preconditioner. We can change the preconditioner by setting "pc_type" to some of the other preconditioners in petsc, which you can find in the PETSc documentation. You can set any opition in PETSc through the petsc_options input, such as the absolute tolerance ("ksp_atol"), relative tolerance ("ksp_rtol") and maximum number of iterations ("ksp_max_it").

cg_problem = LinearProblem(a, L, bcs=bcs,
petsc_options={"ksp_type": "cg", "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1000})
uh = cg_problem.solve()
cg_solver = cg_problem.solver
viewer = PETSc.Viewer().createASCII("cg_output.txt")
cg_solver.view(viewer)
solver_output = open("cg_output.txt", "r")
print(line)

KSP Object: (dolfinx_solve_139960817498864) 1 MPI processes

type: cg

maximum iterations=1000, initial guess is zero

tolerances:  relative=1e-06, absolute=1e-10, divergence=10000.

left preconditioning

using PRECONDITIONED norm type for convergence test

PC Object: (dolfinx_solve_139960817498864) 1 MPI processes

type: ilu

out-of-place factorization

0 levels of fill

tolerance for zero pivot 2.22045e-14

matrix ordering: natural

factor fill ratio given 1., needed 1.

Factored matrix follows:

Mat Object: 1 MPI processes

type: seqaij

rows=961, cols=961

package used to perform factorization: petsc

total: nonzeros=6481, allocated nonzeros=6481

not using I-node routines

linear system matrix = precond matrix:

Mat Object: (dolfinx_solve_139960817498864) 1 MPI processes

type: seqaij

rows=961, cols=961

total: nonzeros=6481, allocated nonzeros=6481

total number of mallocs used during MatSetValues calls=0

not using I-node routines


For non-symmetrix problems, a Krylov solver for non-symmetrix systems, such as GMRES is a better.

gmres_problem = LinearProblem(a, L, bcs=bcs,
petsc_options={"ksp_type": "gmres", "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1000, "pc_type": "none"})
uh = gmres_problem.solve()
gmres_solver = gmres_problem.solver
viewer = PETSc.Viewer().createASCII("gmres_output.txt")
gmres_solver.view(viewer)
solver_output = open("gmres_output.txt", "r")
print(line)

KSP Object: (dolfinx_solve_139960817495744) 1 MPI processes

type: gmres

restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement

happy breakdown tolerance 1e-30

maximum iterations=1000, initial guess is zero

tolerances:  relative=1e-06, absolute=1e-10, divergence=10000.

left preconditioning

using PRECONDITIONED norm type for convergence test

PC Object: (dolfinx_solve_139960817495744) 1 MPI processes

type: none

linear system matrix = precond matrix:

Mat Object: (dolfinx_solve_139960817495744) 1 MPI processes

type: seqaij

rows=961, cols=961

total: nonzeros=6481, allocated nonzeros=6481

total number of mallocs used during MatSetValues calls=0

not using I-node routines


A remark regarding verification using iterative solvers

When we consider manufactured solutions where we expect the resulting error to be of machine precision, it gets complicated when we use iterative methods. The problem is to keep the error due to the iterative solution smaller than the tolerance used in the iterative test. For linear elements and small meshes, a tolerance of between $$10^{-11}$$ and $$10^{-12}$$ works well in the case of Krylov solvers too.