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.
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)
f = -div(grad(u_ufl(x)))
Next, we define our linear variational problem
V = FunctionSpace(mesh, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
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")
for line in solver_output.readlines():
print(line)
KSP Object: (dolfinx_solve_140368388634464) 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_140368388634464) 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_140368388634464) 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")
for line in solver_output.readlines():
print(line)
KSP Object: (dolfinx_solve_140368388824832) 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_140368388824832) 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_140368388824832) 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")
for line in solver_output.readlines():
print(line)
KSP Object: (dolfinx_solve_140368388817776) 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_140368388817776) 1 MPI processes
type: none
linear system matrix = precond matrix:
Mat Object: (dolfinx_solve_140368388817776) 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.