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 mostly revisit the nonlinear time dependent problem studied at the end of the concepts section which after discretizing in time had the variational formulation \begin{equation} \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} \end{equation} 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 concepts section

[1]:
import numpy, sys, io
import matplotlib.pyplot as plt

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
from dune.fem import integrate
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 at the end of this section.

[2]:
scheme = solutionScheme(a == b, solver='cg')

endTime    = 0.25
exact_end  = exact(endTime)
l2error = gridFunction(name="l2error", expr=dot(u_h - exact_end, u_h - exact_end))
h1error = gridFunction(name="h1error", expr=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 the non-linear problem for \(u^{n+1}\) given \(u^n\):

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

We can simply use the galerkinScheme instance 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]:
dt.value = 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([l2error,h1error])]
    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)
        dt.value /= 2
Forchheimer: step: 0 , size: 16
         | u_h - u | = 1.55256e-04 , eoc = -
         | grad(uh - u) | = 3.99906e-03 , eoc = -
_images/solvers_nb_7_1.jpg
Forchheimer: step: 1 , size: 64
         | u_h - u | = 1.92874e-05 , eoc = 3.01
         | grad(uh - u) | = 9.99164e-04 , eoc = 2.0
_images/solvers_nb_7_3.jpg

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.

The most important step is accessing the data structures setup on the C++ side in Python. In this case we would like to use the underlying dof vector from a discrete function as numpy arrays and system matrices assembled by the schemes and operators as scipy sparse matrices. In the introduction we already discussed the as_numpy method. So

[5]:
vecu_h = u_h.as_numpy

provides access to the underlying dof vector without copying. So changes to the numpy array vecu_h carries over to the discrete function u_h. Just remember to make changes using vecu_h[:] to change the actual memory buffer.

A scheme describing an operator L provides a method linear which returns an object that stores a sparse matrix structure. The object describes the operator linearized around zero. To linearize around a different value use the jacobian method on the scheme that linearized the the operator L around a given grid function ubar and fills the linear operator structure passed in as second argument. It is also possible to pass assemble=False to the linear method to avoid an the linearization around zero to reduce computational cost:

[6]:
linOp = scheme.linear()                  # linearized around 0
linOp = scheme.linear(assemble=False)    # empty (non valid) linear operator
scheme.jacobian(space.zero, linOp)       # linearized around zero

Here we linearize around zero. But that argument could be any grid function. A second version of this method will return an addition discrete function rhs which equals -L[ubar] such that DL[ubar](u-ubar) - rhs are the first terms in the Taylor expansion of the operator L:

[7]:
rhs = u_h.copy()
scheme.jacobian(u_h, linOp, rhs)

One can now easily access the underlying sparse matrix by again using as_numpy (and again the underlying data buffers are not copied):

[8]:
A = linOp.as_numpy
print(type(A))
# plt.spy(A, precision=1e-8, markersize=1)
<class 'scipy.sparse._csr.csr_matrix'>

Now we have all the ingredients to write a simple Newton solver to solve our non-linear time dependent PDE.

[9]:
import numpy as np
from scipy.sparse.linalg import spsolve as solver
class Scheme:
  def __init__(self, scheme):
      self.model = scheme.model
      self.jacobian = scheme.linear()

  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([error**2,inner(grad(error),grad(error))]) ]))
Forchheimer(numpy) size:  64 L^2, H^1 error: 1.93091e-05, 9.99162e-04

We can also use a non linear solver from the Scipy package

[10]:
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 = scheme.linear()
            self.update(x_coeff,None)
        # 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, atol=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([error**2,inner(grad(error),grad(error))]) ]))
Forchheimer(scipy) size:  64 L^2, H^1 error: 1.93091e-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.

Tip

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. See here for a summary of the concepts and API for operators and schemes.

[11]:
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)
sol = space.interpolate(0, name="u_h")
rhs = sol.copy()
lin = op.linear()

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 in the section on more general boundary conditions.

[12]:
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(lin.as_numpy, rhs.as_numpy, x0=sol.as_numpy,
                            callback=cb, tol=1e-10, atol=1e-10)
sol.plot()
22.364801469753214      6.5503071305549865      5.5020389790076845      2.651703574864195       2.847846537860448       1.3560658457805435      1.6365176519674605      1.1996200851668588      0.849170828428025       0.4116851221505151      0.3318991774651997      0.2518006403242052      0.23615171539855195     0.22530579266807363     0.2214195904289021      0.21926871065555859     0.21829153537638452     0.21653684764907868     0.21675782107045938     0.2158110638718003      0.21559315250503033     0.21544432425537108     0.21544150394440717     0.2153861339325574      0.21538859378926503     0.21538505227526789     0.21538418947503998     0.21538241666417576     0.21538282256904034     0.21538287372050155     0.21538276135478887     0.21538277558134716     0.21538281353457012     0.21538283870317626     0.2153828390521391      0.21538283539716077     0.21538283650631837     0.21538283693719715     0.2153828375348075      0.21538283758833912     0.21538283758714327     0.21538283755554377
_images/solvers_nb_24_1.jpg

Using PETSc and Petsc4Py

The following requires that a PETSc installation was found during the configuration of dune. Furthermore some examples make use of the Python package `petsc4py.

[13]:
from dune.common.checkconfiguration import assertCMakeHave, ConfigurationError
try:
    assertCMakeHave("HAVE_PETSC")
    petsc = True
except ConfigurationError:
    print("Dune not configured with petsc - skipping example")
    petsc = False
try:
    import petsc4py
    petsc4py.init(sys.argv)
    from petsc4py import PETSc
except ModuleNotFoundError:
    print("petsc4py module not found -- skipping example")
    petsc4py = None

Switching to a storage based on the PETSc solver package and solving the system using the dune-fem bindings can be achieved by using the storage argument to the space constructor

[14]:
if petsc:
    spacePetsc = solutionSpace(gridView, order=2, storage='petsc')
    # 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"})
    dt.value = 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([error**2,inner(grad(error),grad(error))]) ]))
Forchheimer(petsc) size:  64 L^2, H^1 error: 1.93079e-05, 9.99164e-04

Implementing a Newton Krylov solver using the binding provided by petsc4py

[15]:
if petsc4py is not None and petsc:
    class Scheme3:
      def __init__(self, scheme):
          self.model = scheme.model
          self.jacobian = scheme.linear()
          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([error**2,inner(grad(error),grad(error))]) ]))
Forchheimer(petsc) size:  64 L^2, H^1 error: 1.93091e-05, 9.99162e-04

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

[16]:
if petsc4py is not None and petsc 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 = scheme.linear()
            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([error**2,inner(grad(error),grad(error))]) ]))
Forchheimer(petsc4py) size:  64 L^2, H^1 error: 1.93091e-05, 9.99162e-04

Accessing and reusing values of parameters

Sometimes it is necessary to extract which parameters were read and which values were used, e.g., for debugging purposes like finding spelling in the parameters provided to a scheme. Note that this information can only be reliably obtained after usage of the scheme, e.g., after calling solve as shown in the example below. To add logging to a set of parameters passed to a scheme one simply needs to add a logging key to the parameter dictionary provided to the scheme with a tag (string) that is used in the output.

As an example we will solve the simple Laplace equation from the introduction but pass some preconditioning parameters to the scheme.

[17]:
import dune.fem
from dune.grid import structuredGrid
from dune.fem.space import lagrange
from dune.fem.scheme import galerkin
from ufl import (TestFunction, TrialFunction, SpatialCoordinate,
                 dx, grad, inner, dot, sin, cos, pi )
gridView = structuredGrid([0, 0], [1, 1], [200, 200])
space = lagrange(gridView, order=1, storage="istl")
u_h   = space.interpolate(0, name='u_h')
x = SpatialCoordinate(space)
u = TrialFunction(space)
v = TestFunction(space)

f = (8*pi**2+1) * cos(2*pi*x[0])*cos(2*pi*x[1])
a = ( inner(grad(u),grad(v)) + u*v ) * dx
l = f*v * dx
scheme = galerkin( a==l, solver="cg", parameters=
                   {"newton.linear.tolerance": 1e-12,
                    "newton.linear.verbose": True,
                    "newton.linear.preconditioning.method": "amg-ilu",
                    "fem.solver.newton.linear.errormeasure": "relative",
                    "logging": "precon-amg"
                   } )
info = scheme.solve(target=u_h)
=== Dune::CGSolver
 Iter          Defect            Rate
    0          0.19886
    1         0.138041         0.694164
    2        0.0533951         0.386805
    3        0.0335279         0.627921
    4       0.00936836          0.27942
    5        0.0158677          1.69376
    6       0.00553746         0.348977
    7       0.00407438         0.735784
    8       0.00238342         0.584978
    9       0.00097925         0.410859
   10       0.00030568         0.312157
   11      0.000177396         0.580333
   12      4.54516e-05         0.256215
   13      2.72129e-05         0.598721
   14        2.549e-05         0.936688
   15      1.32301e-05          0.51903
   16      7.51289e-06         0.567865
   17      2.13798e-06         0.284575
   18      1.08267e-06         0.506398
   19      7.67007e-07         0.708442
   20      2.79923e-07         0.364955
   21      5.71861e-08         0.204292
   22       2.5419e-08         0.444496
   23      1.38267e-08         0.543951
   24      6.60951e-09         0.478026
   25      1.93829e-09         0.293257
   26      6.68847e-10         0.345072
   27      3.43101e-10         0.512974
   28      1.06161e-10         0.309416
   29       2.6079e-11         0.245655
   30      4.90648e-12         0.188139
   31      6.17261e-12          1.25805
   32      2.10807e-12         0.341519
   33      9.47557e-13         0.449491
=== rate=0.453848, T=0.228761, TIT=0.00693215, IT=33

We use the pprint (pretty print) module if available to get nicer output.

[18]:
try:
    from pprint import pprint as _pprint
    pprint = lambda *args,**kwargs: _pprint(*args,**kwargs,width=200,compact=False)
except ImportError:
    pprint = print

pprint(dune.fem.parameter.log())
{'default': {('fem.dofmanager.memoryfactor', '1.1'), ('fem.threads.communicationthread', 'false'), ('fem.dofmanager.clearresizedarrays', 'true')},
 'precon-amg': {('fem.solver.newton.lineSearch', 'none'),
                ('fem.solver.newton.linear.errormeasure', 'absolute'),
                ('fem.solver.newton.linear.matrix.overflowfraction', '1'),
                ('fem.solver.newton.linear.maxiterations', '2147483647'),
                ('fem.solver.newton.linear.preconditioning.iterations', '1'),
                ('fem.solver.newton.linear.preconditioning.relaxation', '1.1'),
                ('fem.solver.newton.linear.threading', 'true'),
                ('fem.solver.newton.linear.tolerance.strategy', 'none'),
                ('fem.solver.newton.maxiterations', '2147483647'),
                ('fem.solver.newton.maxlinesearchiterations', '2147483647'),
                ('fem.solver.newton.tolerance', '1e-06'),
                ('fem.solver.newton.verbose', 'false'),
                ('newton.linear.method', 'cg'),
                ('newton.linear.preconditioning.method', 'amg-ilu'),
                ('newton.linear.tolerance', '1e-12'),
                ('newton.linear.verbose', 'True')},
 'program code': {('fem.adaptation.method', 'callback')}}

Note above that all parameters are printed including some default ones used in other parts of the code. If multiple schemes with different logging parameter strings are used, all would be shown using the log method as shown above. To access only the parameters used in the scheme simply use either dune.fem.parameter.log()["tag"]) or access the parameter log through the scheme:

[19]:
pprint(scheme.parameterLog())
{('fem.solver.newton.lineSearch', 'none'),
 ('fem.solver.newton.linear.errormeasure', 'absolute'),
 ('fem.solver.newton.linear.matrix.overflowfraction', '1'),
 ('fem.solver.newton.linear.maxiterations', '2147483647'),
 ('fem.solver.newton.linear.preconditioning.iterations', '1'),
 ('fem.solver.newton.linear.preconditioning.relaxation', '1.1'),
 ('fem.solver.newton.linear.threading', 'true'),
 ('fem.solver.newton.linear.tolerance.strategy', 'none'),
 ('fem.solver.newton.maxiterations', '2147483647'),
 ('fem.solver.newton.maxlinesearchiterations', '2147483647'),
 ('fem.solver.newton.tolerance', '1e-06'),
 ('fem.solver.newton.verbose', 'false'),
 ('newton.linear.method', 'cg'),
 ('newton.linear.preconditioning.method', 'amg-ilu'),
 ('newton.linear.tolerance', '1e-12'),
 ('newton.linear.verbose', 'True')}

One can easily reuse these parameters to construct another scheme by converting the result of the above call to a dictionary. As an example change the above problem to a PDE with Dirichlet conditions but turn of verbose output of the solver.

Note

the logging parameter has to be set if we want to use the parameterLog method on the scheme.

[20]:
param = dict(scheme.parameterLog()) # this method returns a set of pairs which we can convert to a dictionary
param["logging"] = "Dirichlet" # only needed to use the `parameterLog` method
param["newton.linear.verbose"] = False
scheme2 = galerkin( [a==l,DirichletBC(space,0)], parameters=param )
u_h.clear()
info = scheme2.solve(target=u_h)
pprint(scheme2.parameterLog())
{('fem.solver.newton.lineSearch', 'none'),
 ('fem.solver.newton.linear.errormeasure', 'absolute'),
 ('fem.solver.newton.linear.matrix.overflowfraction', '1'),
 ('fem.solver.newton.linear.maxiterations', '2147483647'),
 ('fem.solver.newton.linear.preconditioning.iterations', '1'),
 ('fem.solver.newton.linear.preconditioning.relaxation', '1.1'),
 ('fem.solver.newton.linear.threading', 'true'),
 ('fem.solver.newton.linear.tolerance.strategy', 'none'),
 ('fem.solver.newton.maxiterations', '2147483647'),
 ('fem.solver.newton.maxlinesearchiterations', '2147483647'),
 ('fem.solver.newton.tolerance', '1e-06'),
 ('fem.solver.newton.verbose', 'false'),
 ('newton.linear.method', 'cg'),
 ('newton.linear.preconditioning.method', 'amg-ilu'),
 ('newton.linear.tolerance', '1e-12'),
 ('newton.linear.verbose', 'False')}

Parameter hints

Tip

To get information about available values for some parameters (those with string arguments) a possible approach is to provide a non valid string, e.g., "help".

[21]:
scheme = galerkin( a==l, solver="cg", parameters=
                   {"newton.linear.tolerance": 1e-12,
                    "newton.linear.verbose": True,
                    "newton.linear.preconditioning.method": "help",
                    "fem.solver.newton.linear.errormeasure": "relative",
                    "logging": "precon-amg"
                   } )
try:
    scheme.solve(target=u_h)
except RuntimeError as rte:
    print(rte)
ParameterInvalid [getEnumeration:/home/dedner/DUNE_PYTHON/dune-fem/dune/fem/io/parameter/reader.hh:301]: Parameter 'fem.solver.newton.linear.preconditioning.method' invalid.

Parameter 'fem.solver.newton.linear.preconditioning.method' invalid.
Valid values are: none, ssor, sor, ilu, gauss-seidel, jacobi, amg-ilu, amg-jacobi, ildl

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.

[22]:
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.csr_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:

[23]:
help(scheme)
Help on Scheme in module dune.generated.femscheme_9b7005253a7d0b23b665ce4b07198dde 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__(...)
 |
 |  dirichletIndices = _opDirichletIndices(self, id=None)
 |
 |  inverseLinearOperator(...)
 |
 |  jacobian(...)
 |
 |  linear = _schemeLinear(self, assemble=True, parameters=None)
 |
 |  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_9b7005253a7d0...
 |
 |
 |  ----------------------------------------------------------------------
 |  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 DOI