# Alternative Linear Solvers (Scipy and Petsc)

Here we look at different ways of solving PDEs using external packages and python functionality. Different linear algebra backends can be accessed by changing setting the storage parameter during construction of the discrete space. All discrete functions and operators/schemes based on this space will then use this backend. Available backends are numpy,istl,petsc. The default is numpy which uses simple data structures and linear solvers implemented in the dune-fem package. The simplicity of the data structure makes it possible to use the buffer protocol to seamlessly move between C++ and Numpy/Scipy data structures on the python side. A degrees of freedom vector (dof vector) can be retrieved from a discrete function over the numpy space by using the as_numpy method. Similar methods are available for the other storages, i.e., as_istl,as_petsc. The same methods are also available to retrieve the underlying matrix structures of linear operators.

We will revisit the nonlinear time dependent problem studied in the introduction which after discretizing in time had the variational formulation $$\begin{split} \int_{\Omega} \frac{u^{n+1}-u^n}{\Delta t} \varphi + \frac{1}{2}K(\nabla u^{n+1}) \nabla u^{n+1} \cdot \nabla \varphi \ + \frac{1}{2}K(\nabla u^n) \nabla u^n \cdot \nabla \varphi v\ dx \\ - \int_{\Omega} \frac{1}{2}(f(x,t^n)+f(x,t^n+\Delta t) \varphi\ dx - \int_{\partial \Omega} \frac{1}{2}(g(x,t^n)+g(x,t^n+\Delta t)) v\ ds = 0. \end{split}$$ on a domain $$\Omega=[0,1]^2$$. We choose $$f,g$$ so that the exact solution is given by \begin{align*} u(x,t) = e^{-2t}\left(\frac{1}{2}(x^2 + y^2) - \frac{1}{3}(x^3 - y^3)\right) + 1 \end{align*} The following code was described in the introduction:

[1]:

import matplotlib
matplotlib.rc( 'image', cmap='jet' )
import numpy, sys, io

from dune.grid import structuredGrid as leafGridView
from dune.fem.space import lagrange as solutionSpace
from dune.fem.scheme import galerkin as solutionScheme
from dune.fem.function import gridFunction, integrate, uflFunction
from dune.ufl import Constant
from ufl import TestFunction, TrialFunction, SpatialCoordinate, FacetNormal, \
dx, ds, div, grad, dot, inner, sqrt, exp, sin,\
conditional

gridView = leafGridView([0, 0], [1, 1], [4, 4])
space = solutionSpace(gridView, order=2)

x = SpatialCoordinate(space)
initial = 1/2*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1
exact   = lambda t: exp(-2*t)*(initial - 1) + 1

u_h   = space.interpolate(initial, name='u_h')
u_h_n = u_h.copy(name="previous")

u = TrialFunction(space)
v = TestFunction(space)
dt = Constant(0, name="dt")    # time step
t  = Constant(0, name="t")     # current time

abs_du = lambda u: sqrt(inner(grad(u), grad(u)))
K = lambda u: 2/(1 + sqrt(1 + 4*abs_du(u)))
a = ( dot((u - u_h_n)/dt, v) \
+ 0.5*dot(K(u)*grad(u), grad(v)) \
+ 0.5*dot(K(u_h_n)*grad(u_h_n), grad(v)) ) * dx

f = lambda s: -2*exp(-2*s)*(initial - 1) - div( K(exact(s))*grad(exact(s)) )
g = lambda s: K(exact(s))*grad(exact(s))
n = FacetNormal(space)
b = 0.5*(f(t)+f(t+dt))*v*dx + 0.5*dot(g(t)+g(t+dt),n)*v*ds


When creating a scheme, it is possible to set the linear solver as well as parameters for the internal Newton solver and the linear solver and preconditioning. See a list of available solvers and preconditioning methods here.

[2]:

scheme = solutionScheme(a == b, solver='cg')

endTime    = 0.25
exact_end  = exact(endTime)
l2error = uflFunction(gridView, name="l2error", order=u_h.space.order, ufl=dot(u_h - exact_end, u_h - exact_end))
h1error = uflFunction(gridView, name="h1error", order=u_h.space.order, ufl=dot(grad(u_h - exact_end), grad(u_h - exact_end)))


We define a function to evolve the solution from time 0 to the end time. The first argument is a class with a solve method that moves the solution from one time level to the next - i.e., solves for $$u^{n+1}$$ given $$u^n$$:

[3]:

def evolve(scheme, u_h, u_h_n, endTime):
time = 0
while time < (endTime - 1e-6):
scheme.model.t = time
u_h_n.assign(u_h)
scheme.solve(target=u_h)
time += scheme.model.dt


We can simply use the scheme in this function to produce the solution at the final time. We combine this with a loop to compute the error over two grids and estimate the convergence rate:

[4]:

scheme.model.dt = 0.005

errors = 0,0
loops = 2
for eocLoop in range(loops):
u_h.interpolate(initial)
evolve(scheme, u_h, u_h_n, endTime)
errors_old = errors
errors = [sqrt(e) for e in integrate(gridView, [l2error,h1error], order=5)]
if eocLoop == 0:
eocs = ['-','-']
else:
eocs = [ round(numpy.log(e/e_old)/numpy.log(0.5),2) \
for e,e_old in zip(errors,errors_old) ]
print('Forchheimer: step:', eocLoop, ', size:', gridView.size(0))
print('\t | u_h - u | =', '{:0.5e}'.format(errors[0]), ', eoc =', eocs[0])
print('\t | grad(uh - u) | =', '{:0.5e}'.format(errors[1]), ', eoc =', eocs[1])
u_h.plot()
if eocLoop < loops-1:
gridView.hierarchicalGrid.globalRefine(1)
scheme.model.dt /= 2

Forchheimer: step: 0 , size: 16
| u_h - u | = 1.30293e-04 , eoc = -
| grad(uh - u) | = 3.99906e-03 , eoc = -

Forchheimer: step: 1 , size: 64
| u_h - u | = 1.61439e-05 , eoc = 3.01
| grad(uh - u) | = 9.99164e-04 , eoc = 2.0


## Using Scipy

We implement a simple Newton Krylov solver using a linear solver from Scipy. We can use the as_numpy method to access the degrees of freedom as Numpy vector based on the python buffer protocol. So no data is copied and changes to the dofs made on the python side are automatically carried over to the C++ side. from Scipy.

[5]:

from dune.fem.operator import linear as linearOperator
import numpy as np
from scipy.sparse.linalg import spsolve as solver
class Scheme:
def __init__(self, scheme):
self.model = scheme.model
self.jacobian = linearOperator(scheme)

def solve(self, target):
# create a copy of target for the residual
res = target.copy(name="residual")

# extract numpy vectors from target and res
sol_coeff = target.as_numpy
res_coeff = res.as_numpy

n = 0
while True:
scheme(target, res)
absF = numpy.sqrt( np.dot(res_coeff,res_coeff) )
if absF < 1e-10:
break
scheme.jacobian(target,self.jacobian)
sol_coeff -= solver(self.jacobian.as_numpy, res_coeff)
n += 1

scheme_cls = Scheme(scheme)

u_h.interpolate(initial)                # reset u_h to initial
evolve(scheme_cls, u_h, u_h_n, endTime)
error = u_h - exact_end
print("Forchheimer(numpy) size: ", gridView.size(0), "L^2, H^1 error:",'{:0.5e}, {:0.5e}'.format(
*[ sqrt(e) for e in integrate(gridView,[error**2,inner(grad(error),grad(error))], order=5) ]))

Forchheimer(numpy) size:  64 L^2, H^1 error: 1.61699e-05, 9.99162e-04


Using a non linear solver from the Scipy package

[6]:

from scipy.optimize import newton_krylov
from scipy.sparse.linalg import LinearOperator
from scipy.sparse.linalg import cg as solver

class Scheme2:
def __init__(self, scheme):
self.scheme = scheme
self.model = scheme.model
self.res = u_h.copy(name="residual")

# non linear function
def f(self, x_coeff):
# the following converts a given numpy array
# into a discrete function over the given space
x = space.function("tmp", dofVector=x_coeff)
scheme(x, self.res)
return self.res.as_numpy

# class for the derivative DS of S
class Df(LinearOperator):
def __init__(self, x_coeff):
self.shape = (x_coeff.shape[0], x_coeff.shape[0])
self.dtype = x_coeff.dtype
x = space.function("tmp", dofVector=x_coeff)
self.jacobian = linearOperator(scheme, ubar=x)
# reassemble the matrix DF(u) given a DoF vector for u
def update(self, x_coeff, f):
x = space.function("tmp", dofVector=x_coeff)
scheme.jacobian(x, self.jacobian)
# compute DS(u)^{-1}x for a given DoF vector x
def _matvec(self, x_coeff):
return solver(self.jacobian.as_numpy, x_coeff, tol=1e-10)[0]

def solve(self, target):
sol_coeff = target.as_numpy
# call the newton krylov solver from scipy
sol_coeff[:] = newton_krylov(self.f, sol_coeff,
verbose=0, f_tol=1e-8,
inner_M=self.Df(sol_coeff))

scheme2_cls = Scheme2(scheme)
u_h.interpolate(initial)
evolve(scheme2_cls, u_h, u_h_n, endTime)
error = u_h - exact_end
print("Forchheimer(scipy) size: ", gridView.size(0), "L^2, H^1 error:",'{:0.5e}, {:0.5e}'.format(
*[ sqrt(e) for e in integrate(gridView,[error**2,inner(grad(error),grad(error))], order=5) ]))

Forchheimer(scipy) size:  64 L^2, H^1 error: 1.61699e-05, 9.99162e-04


## Handling Dirichlet boundary conditions

We look at a simple Poisson problem with Dirichlet BCs to show how to use external solvers like the cg method from Scipy in this case. We solve $$-\triangle u=10\chi_\omega$$ where $$\chi_\omega$$ is a characteristic function with $$\omega=\{x\colon |x|^2<0.6\}$$. For the boundary we prescribe trivial Neuman at the top and bottom boundaries and Dirichlet values $$u=-1$$ and $$u=1$$ at the left and right boundaries, respectively. We will use the CG solver from scipy.sparse.linalg:

Note: since we are not needing to invert the operator we will use the dune.fem.operator.galerkin class to setup the problem. This is similar to dune.fem.scheme.galerkin we have been using so far but can be used to model operators between different spaces.

[7]:

from dune.ufl import DirichletBC
from dune.fem.operator import galerkin
from scipy.sparse.linalg import cg as solver
model  = ( inner(grad(u), grad(v)) -
conditional(dot(x,x)<0.6,10.,0.) * v ) * dx
dbcs   = [ DirichletBC(space,-1,1),
DirichletBC(space, 1,2) ]
op     = galerkin([model, *dbcs], space)
A = linearOperator(op).as_numpy
sol = space.interpolate(0, name="u_h")
rhs = sol.copy()
op(sol, rhs)
rhs.as_numpy[:] *= -1


So far everything is as before. Dirichlet boundary conditions are handled in the matrix through changing all rows associated with boundary degrees of freedom to unit rows - associated columns are not changed so the matrix will not be symmetric anymore. For solving the system we need to modify the right hand side and the initial guess for the iterative solver to include the boundary values (to counter the missing symmetry). We can use the first of the three versions of the setConstraints methods on the scheme class discussed previously.

[8]:

op.setConstraints(rhs)
op.setConstraints(sol)
rk = sol.copy("residual")
def cb(xk): # a callback to print the residual norm in each step
x_h = space.function("iterate", dofVector=xk)
op(x_h,rk)
print(rk.as_numpy[:].dot(rk.as_numpy[:]), flush=True, end='\t')
sol.as_numpy[:], _ = solver(A, rhs.as_numpy, x0=sol.as_numpy,
callback=cb, tol=1e-10)
sol.plot()

21.55296640308155       6.325538337303643       4.953832681894503       2.517384423575954       2.5399150209107493      1.2556383326063998      1.5569218960849052      1.2786201081471515      0.9125357940577181      0.3567921408073909      0.2599475839416715      0.14086624759863325     0.1575839322086894      0.2104364147699787      0.2209842133400043      0.36115057029326153     0.22948222295623288     0.08966742907847959     0.06571321824701241     0.027183843287078593    0.009725308020726454    0.008752431767516089    0.007409479346256096    0.012270132675777288    0.0059332308589063265   0.0022248065677152105   0.003570461229314596    0.004905376772491767    0.003104404010258965    0.0018845766417821344   0.00042391511594883096  0.00019993006312547264  8.035559607735921e-05   3.939121668647674e-05   5.3498421755238035e-05  3.8926838321958035e-05  1.0753629248137449e-05  5.197442708772677e-06   3.315808852882241e-06   3.4676669782103684e-06  1.949682814935511e-06   7.432979169339984e-07   3.2758055946840397e-07  1.151364056864077e-07   7.510339459465767e-08   3.830090417816041e-08   1.3063486034709437e-08  5.799521579470168e-09   2.3524139712879676e-09  1.1619915790661837e-09  6.2865080163072e-10     1.5449859639395997e-10  5.350843041337233e-11   4.571537742091988e-11   1.840069336259542e-11   1.0276595193853602e-11  5.328734121365251e-12   4.313640880303699e-12   3.3080552059697887e-12  8.007039260332559e-13   4.816491317214374e-13   5.265750545323643e-13   3.8623256183123157e-13  1.6186073147441346e-13  6.662699914788252e-14   2.3766576368020617e-14  1.1203094816854793e-14  3.751271698934124e-15   1.159769006787277e-15   4.253765267616772e-16   1.130018017169068e-16   4.4139679901149106e-17  6.4813975686812045e-18  2.4436655787619954e-18  8.767380574109113e-19   5.894905001239205e-19   1.382737290474406e-19


## Using Petsc and Petsc4Py

Switching to a storage based on the PETSc solver package and solving the system using the dune-fem bindings

[9]:

from dune.generator import ConfigurationError
try:
import petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc
spacePetsc = solutionSpace(gridView, order=2, storage='petsc')
except ModuleNotFoundError:
print("petsc4py not found: skipping example")
print("petsc4py module not found so skipping example - ignored")
petsc4py = None
except ConfigurationError:
print("petsc4py found but petsc was not found during configuration of dune")
petsc4py = None
pass

if petsc4py is not None:
# first we will use the petsc solver available in the dune-fem package
# (using the sor preconditioner)
schemePetsc = solutionScheme(a == b, space=spacePetsc,
parameters={"linear.preconditioning.method":"sor"})
schemePetsc.model.dt = scheme.model.dt
u_h = spacePetsc.interpolate(initial, name='u_h')
u_h_n = u_h.copy(name="previous")
evolve(schemePetsc, u_h, u_h_n, endTime)
error = u_h - exact_end
print("Forchheimer(petsc) size: ", gridView.size(0), "L^2, H^1 error:",'{:0.5e}, {:0.5e}'.format(
*[ sqrt(e) for e in integrate(gridView,[error**2,inner(grad(error),grad(error))], order=5) ]))
else:
print("petsc module not found so skipping example - ignored")
spacePetsc = None
pass

Forchheimer(petsc) size:  64 L^2, H^1 error: 1.61684e-05, 9.99164e-04


Implementing a Newton Krylov solver using the binding provided by petsc4py

[10]:

if petsc4py is not None and spacePetsc is not None:
class Scheme3:
def __init__(self, scheme):
self.model = scheme.model
self.jacobian = linearOperator(scheme)
self.ksp = PETSc.KSP()
self.ksp.create(PETSc.COMM_WORLD)
# use conjugate gradients method
self.ksp.setType("cg")
# and incomplete Cholesky
self.ksp.getPC().setType("icc")
self.ksp.setOperators(self.jacobian.as_petsc)
self.ksp.setFromOptions()
def solve(self, target):
res = target.copy(name="residual")
sol_coeff = target.as_petsc
res_coeff = res.as_petsc
n = 0
while True:
schemePetsc(target, res)
absF = numpy.sqrt( res_coeff.dot(res_coeff) )
if absF < 1e-10:
break
schemePetsc.jacobian(target, self.jacobian)
self.ksp.solve(res_coeff, res_coeff)
sol_coeff -= res_coeff
n += 1

u_h.interpolate(initial)
scheme3_cls = Scheme3(schemePetsc)
evolve(scheme3_cls, u_h, u_h_n, endTime)
error = u_h - exact_end
print("Forchheimer(petsc) size: ", gridView.size(0), "L^2, H^1 error:",'{:0.5e}, {:0.5e}'.format(
*[ sqrt(e) for e in integrate(gridView,[error**2,inner(grad(error),grad(error))], order=5) ]))

Forchheimer(petsc) size:  64 L^2, H^1 error: 1.61699e-05, 9.99162e-04


Using the petsc4py bindings for the non linear KSP solvers from PETSc

[11]:

if petsc4py is not None and spacePetsc is not None:
class Scheme4:
def __init__(self, scheme):
self.model = scheme.model
self.res = scheme.space.interpolate([0],name="residual")
self.scheme = scheme
self.jacobian = linearOperator(self.scheme)
self.snes = PETSc.SNES().create()
self.snes.setFunction(self.f, self.res.as_petsc.duplicate())
self.snes.setUseMF(False)
self.snes.setJacobian(self.Df, self.jacobian.as_petsc, self.jacobian.as_petsc)
self.snes.getKSP().setType("cg")
self.snes.setFromOptions()

def f(self, snes, x, f):
# setup discrete function using the provide petsc vectors
inDF = self.scheme.space.function("tmp",dofVector=x)
outDF = self.scheme.space.function("tmp",dofVector=f)
self.scheme(inDF,outDF)

def Df(self, snes, x, m, b):
inDF = self.scheme.space.function("tmp",dofVector=x)
self.scheme.jacobian(inDF, self.jacobian)
return PETSc.Mat.Structure.SAME_NONZERO_PATTERN

def solve(self, target):
sol_coeff = target.as_petsc
self.res.clear()
self.snes.solve(self.res.as_petsc, sol_coeff)

u_h.interpolate(initial)
scheme4_cls = Scheme4(schemePetsc)
evolve(scheme4_cls, u_h, u_h_n, endTime)
error = u_h - exact_end
print("Forchheimer(petsc4py) size: ", gridView.size(0), "L^2, H^1 error:",'{:0.5e}, {:0.5e}'.format(
*[ sqrt(e) for e in integrate(gridView,[error**2,inner(grad(error),grad(error))], order=5) ]))

Forchheimer(petsc4py) size:  64 L^2, H^1 error: 1.61699e-05, 9.99162e-04


## Available solvers and parameters

Upon creation of a discrete function space one also have to specifies the storage which is tied to the solver backend. As mentioned, different linear algebra backends can be accessed by changing setting the storage parameter during construction of the discrete space. All discrete functions and operators/schemes based on this space will then use this backend. Available backends are numpy,istl,petsc. Note that not all methods which are available in dune-istl or PETSc have been forwarded to be used with dune-fem.

[12]:

space = solutionSpace(gridView, order=2, storage='numpy')


Switching is as simple as passing storage='istl' or storage='petsc'. Here is a summary of the available backends

Solver

Description

numpy

the storage is based on a raw C pointer which can be

directly accessed as a numpy.array using the Python buffer protocol

To change the underlying vector of a discrete function ‘u_h’ use

‘uh.as_numpy[:]’.

As shown in the examples, linear operators return a

scipy.sparse_matrix through the ‘as_numpy’ property.

istl

data is stored in a block vector/matrix from the dune.istl package

Access through ‘as_istl’

petsc

data is stored in a petsc vector/matrix which can also be used with

petsc4py on the python side using ‘as_petsc’

When creating a scheme, there is the possibility to select a linear solver for the internal Newton method. In addition the behavior of the solver can be customized through a parameter dictionary. This allows to set tolerances, verbosity, but also which preconditioner to use.

For details see the help available for a scheme:

[13]:

help(scheme)

Help on Scheme in module dune.generated.femscheme_1830e0b499d58f20cb34ce353560f035 object:

class Scheme(pybind11_builtins.pybind11_object)
|  A scheme finds a solution u=ufl.TrialFunction for a given variational equation.
|  The main method is solve which takes a discrete functions as target argument to
|  store the solution. The method always uses a Newton method to solve the problem.
|  The linear solver used in each iteration of the Newton method can be chosen
|  using the solver parameter in the constructor of the scheme. Available solvers are:
|  ------------------------------------------
|  |  Solver  |         Storage             |
|  |   name   |  numpy  |  istl   |  petsc  |
|  |----------|---------|---------|---------|
|  | bicg     |   ---   |   ---   |    x    |
|  | bicgstab |    x    |    x    |    x    |
|  | cg       |    x    |    x    |    x    |
|  | gmres    |    x    |    x    |    x    |
|  | gradient |   ---   |    x    |   ---   |
|  | loop     |   ---   |    x    |   ---   |
|  | minres   |   ---   |    x    |    x    |
|  | preonly  |   ---   |   ---   |    x    |
|  | superlu  |   ---   |    x    |   ---   |
|  ------------------------------------------
|
|  In addition the direct solvers from the suitesparse package can be used with the
|  numpy storage. In this case provide a tuple as solver argument with "suitesparse" as
|  first argument and the solver to use as second, e.g.,
|  'solver=("suitesparse","umfpack")'.
|
|  The detailed behavior of the schemes can be customized by providing a
|  parameters dictionary to the scheme constructor, e.g.,
|     {"newton.tolerance": 1e-3, # tolerance for newton solver
|      "newton.verbose": False,  # toggle iteration output
|      "newton.linear.tolerance": 1e-5, # tolerance for linear solver
|      "newton.linear.errormeasure": "absolute", # or "relative" or "residualreduction"
|      "newton.linear.preconditioning.method": "jacobi", # (see table below)
|      "newton.linear.preconditioning.hypre.method": "boomeramg", #  "pilu-t" "parasails"
|      "newton.linear.preconditioning.iteration": 3, # iterations for preconditioner
|      "newton.linear.preconditioning.relaxation": 1.0, # omega for SOR and ILU
|      "newton.linear.maxiterations":1000, # max number of linear iterations
|      "newton.linear.verbose": False,     # toggle linear iteration output
|      "newton.linear.preconditioning.level": 0} # fill-in level for ILU preconditioning
|  -----------------------------------------------
|  |  Precondition | (x = parallel | s = serial) |
|  |  method       |  numpy  |   istl  |  petsc  |
|  |---------------|---------|---------|---------|
|  | amg-ilu       |   ---   |    x    |   ---   |
|  | amg-jacobi    |   ---   |    x    |   ---   |
|  | gauss-seidel  |    x    |    x    |    x    |
|  | hypre         |   ---   |   ---   |    x    |
|  | icc           |   ---   |   ---   |    x    |
|  | ildl          |   ---   |    x    |   ---   |
|  | ilu           |   ---   |    s    |    s    |
|  | jacobi        |    x    |    x    |    x    |
|  | kspoptions    |   ---   |   ---   |    x    |
|  | lu            |   ---   |   ---   |    s    |
|  | ml            |   ---   |   ---   |    x    |
|  | none          |    x    |    x    |    x    |
|  | oas           |   ---   |   ---   |    x    |
|  | pcgamg        |   ---   |   ---   |    x    |
|  | sor           |    x    |    x    |    x    |
|  | ssor          |    x    |    x    |    x    |
|  -----------------------------------------------
|
|  The functionality of some of the preconditioners listed for petsc will
|  depend on the petsc installation.
|
|  Method resolution order:
|      Scheme
|      pybind11_builtins.pybind11_object
|      builtins.object
|
|  Methods defined here:
|
|  __call__(...)
|
|  __init__(...)
|
|  inverseLinearOperator(...)
|
|  jacobian(...)
|
|  setErrorMeasure(...)
|
|  setQuadratureOrders(...)
|
|  solve(scheme, target, rhs=None)
|
|  ----------------------------------------------------------------------
|  Readonly properties defined here:
|
|  cppIncludes
|
|  cppTypeName
|
|  dimRange
|
|  domainSpace
|
|  parameterHelp
|
|  rangeSpace
|
|  space
|
|  ----------------------------------------------------------------------
|  Data descriptors defined here:
|
|  __dict__
|
|  ----------------------------------------------------------------------
|  Data and other attributes defined here:
|
|  LinearInverseOperator = <class 'dune.generated.femscheme_1830e0b499d58...
|
|
|  ----------------------------------------------------------------------
|  Static methods inherited from pybind11_builtins.pybind11_object:
|
|  __new__(*args, **kwargs) from pybind11_builtins.pybind11_type
|      Create and return a new object.  See help(type) for accurate signature.



This page was generated from the notebook solvers_nb.ipynb and is part of the tutorial for the dune-fem python bindings