General Concepts (including how to solve a Laplace problem)

In the following we introduce how to construct and use the basic components of a finite element method ending the first part with solving a simple Laplace problem. The required steps are:

  • constructing a tessellation of the computational domain

  • working with functions defined over the grid

  • setting up discrete functions

  • defining the mathematical model to solve

  • solving the (non linear) system arising from the discretization of the model by the Galerkin method

After that we study the more complex problem of solving a time dependent non linear problem.

Setting up the Mesh

After some general import statements we start our tutorial by setting up a simple Cartesian grid over the domain \([0,1]^2\) subdivided into four intervals in each dimension. We can then show the grid using the plot method - note that direct plotting with MatPlotLib is only supported in 2D, other options for visualization, e.g., using ParaView will be discussed later.

Using the dune.grid.structuredGrid function, is the simplest way to construct a grid, more complicated (unstructured) grids can be constructed using a dictionary containing the vertices and the element connectivity or by for example providing a gmsh file. This will be discussed in later section, see for example here for a 3d gmsh file or here for a simple example using a dictonary.

[1]:
import io, numpy, matplotlib
matplotlib.rc( 'image', cmap='jet' )
from matplotlib import pyplot

from dune.grid import structuredGrid as leafGridView
gridView = leafGridView([0, 0], [1, 1], [4, 4])
gridView.plot()
_images/concepts_nb_2_0.png

Grid Functions

These are function that are defined over a given grid and are evaluated by using an element of the grid and local coordinated within that element:

value = gridFunction(element,localCoordinate)

Alternatively one can obtain a LocalFunction from a grid function which can be bound to an element and then evaluate via local coordinate:

localFunction = gridFunction.localFunction()
for e in grid.elements:
    localFunction.bind(e)
    value = localFunction(x)
    localFunction.unbind()

There are multiple ways to construct grid functions. The easiest way it to use UFL expression. Many methods expecting grid functions are argument can also directly handle UFL expression. We can for example integrate a UFL expression over the grid:

[2]:
from ufl import SpatialCoordinate, triangle
x = SpatialCoordinate(triangle)

exact = 1/2*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1

from dune.fem.function import integrate
print( integrate(gridView, exact, order=5) )
1.333333333333333

and plot them using matplotlib or write a vtk file for postprocessing

[3]:
from dune.fem.plotting import plotPointData as plot
plot(exact, grid=gridView)
gridView.writeVTK('exact', pointdata={'exact': exact})
_images/concepts_nb_6_0.png

In some cases it can be necessary to convert a UFL expression into a grid function explicitly - for example to be able to evaluate it over each element in the grid (in a later section) we provide more detail on how to access the grid interface).

[4]:
from dune.fem.function import uflFunction
exact_gf = uflFunction(gridView, name="ufl", order=1, ufl=exact)
mass = 0
for element in gridView.elements:
    mass += exact_gf(element,[0.5,0.5]) * element.geometry.volume
print(mass)
1.328125

Another way to obtain a grid function is to use the gridFunction decorator. This can be obtained from dune.grid but then without UFL support. Using the decorator from dune.fem.function the resulting grid function can be used seamlessly within UFL expressions:

[5]:
from dune.fem.function import gridFunction
@gridFunction(gridView,name="callback",order=1)
def exactLocal(element,xLocal):
    x = element.geometry.toGlobal(xLocal)
    return 1/2.*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1

we can use the same approach but with a function using global coordinates but can then be used like any other grid function:

[6]:
@gridFunction(gridView,name="callback",order=1)
def exactGlobal(x):
    return 1/2.*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1

lf = exactGlobal.localFunction()

mass = 0
for element in gridView.elements:
    lf.bind(element)
    mass += lf([0.5,0.5]) * element.geometry.volume
    lf.unbind()
print(mass)

print( integrate(gridView, [exact,exactLocal,exactGlobal], order=5) )
1.328125
(1.333333, 1.333333, 1.333333)

As pointed out the dune.fem grid function can be used like any other UFL coefficient to form UFL expressions:

[7]:
print( integrate(gridView, abs(exact-exactLocal), order=5) )
gf = uflFunction(gridView, name="ufl", order=1, ufl=exact+exactLocal*exactGlobal)
fig = pyplot.figure(figsize=(20,10))
gf.plot(figure=(fig,121))
exactGlobal.plot(figure=(fig,122))
0.0
_images/concepts_nb_14_1.png

Converting UFL expressions to grid functions leads to JIT code generation and is therefore efficient when used in other C++ algorithm (like integrate). On the other hand using the gridFunction decorator leads to a callback into Python for each evaluation and is therefore much less efficient. An alternative approach is based on writing small C++ snippets implementing the grid function:

[8]:
from dune.fem.function import cppFunction
code="""
#include <cmath>
#include <dune/common/fvector.hh>
template <class GridView>
auto aTimesExact(double a) {
  // the return value needs to be a lambda which either returns a `double`
  // or a `Dune::FieldVector<double,R>` with `R>=1`
  return [a](const auto& en,const auto& xLocal) -> auto {
    auto x = en.geometry().global(xLocal);
    return a*(1./2.*(std::pow(x[0],2)+std::pow(x[1],2)) - 1./3.*(std::pow(x[0],3) - std::pow(x[1],3)) + 1.);
  };
}
"""
exactCpp = cppFunction(gridView, name="exactCpp", order=2,
                       fctName="aTimesExact",includes=io.StringIO(code),
                       args=[2.])
print( integrate(gridView, abs(2*exact-exactCpp), order=5) )
exactCpp.plot()
0.0
_images/concepts_nb_16_1.png

As the above example shows it is easy to pass in parameters to the C++ implementation - here a double 2. Note that it is not possible to obtain a reference to this double so to make sure a change to the constant on the Python side carries over to the C++ side an option is to use a dune.common.FieldVector instead. These parameters can also be more complex e.g. other grid # function - note that a UFL expression can not be directly passed in - it first needs to be converted into a grid function using uflFunction.

[9]:
code2="""
#include <cmath>
#include <dune/common/fvector.hh>
template <class GridView, class GF>
auto aTimesExact(const GF &gf,Dune::FieldVector<double,1> &a) {
  return [lgf=localFunction(gf),&a] (const auto& en,const auto& xLocal) mutable -> auto {
    lgf.bind(en); // lambda must be mutable so that the non const function can be called
    return a[0]*lgf(xLocal);
  };
}
"""
from dune.common import FieldVector
a = FieldVector([2])
exactCpp2 = cppFunction(gridView, name="exactCpp", order=2,
                        fctName="aTimesExact",includes=io.StringIO(code2),
                        args=[exact_gf,a])
print( integrate(gridView, abs(exactCpp-exactCpp2), order=5) )
a[0] = 4
print( integrate(gridView, abs(2*exactCpp-exactCpp2), order=5) )
0.0
0.0

The above is just one of a few ways C++ code snippets can be used together with Python code to improve efficiency or extend the existing binding to Dune. In the above example there is no advantage of using the C++ code over the code generated based on the UFL expression. For more complicated functions e.g. with many if statements or based on more information from the given element like it’s neighbors the expressibility of UFL might not be sufficient or lead to hard to read code. In these cases directly providing C++ code (or Python code) can be a reasonable alternative.

Discrete Functions

Note that the grid functions set up so far did not involve any discretization, they are exactly evaluated at the given point. A special type of grid functions are discrete functions living in a discrete (finite dimensional) space. .. index:: Functions; discrete functions

[10]:
from dune.fem.space import lagrange as solutionSpace
space = solutionSpace(gridView, order=2)

The easiest way to construct a discrete function is to use the interpolate method on the discrete function space.

[11]:
u_h = space.interpolate(exact, name='u_h')

On an existing discrete function the interpolate method can be used to reinitialize it

[12]:
u_h.interpolate( cppFunction(gridView, name="exactCpp", order=2,
                 fctName="aTimesExact",includes=io.StringIO(code),
                 args=[2.]) )
u_h.interpolate( lambda x: 1/2.*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1 )

Note that in the last example the Python lambda is used as a callback automatically using the same concept used in the gridFunction decorator. As pointed out above there are some methods where these conversions are implicit and no explicit generation of a grid function has be carried out.

If a discrete function is already available it is possible to call copy to obtain further discrete functions:

[13]:
u_h_n = u_h.copy(name="previous")

Finally, clear can be called on a discrete function which sets all coefficient to zero and assign can be used to copy all coefficients between two discrete function over the same space:

[14]:
u_h_n.clear()
u_h_n.assign( u_h )

All the things we did above with grid functions can be done with discrete functions, e.g., evaluate locally

[15]:
localUh = u_h.localFunction()
mass = 0
for element in gridView.elements:
    localUh.bind(element) # using u_h(element,[0.5,0.5]) also works
    mass += localUh([0.5,0.5]) * element.geometry.volume
    localUh.unbind()
print(mass)
1.328125

or plot using matplotlib and write a vtk file for postprocessing (using binary data format to reduce size)

[16]:
u_h.plot(gridLines="white")
from dune.grid import OutputType
gridView.writeVTK('uh', pointdata=[u_h], outputType=OutputType.appendedraw)
_images/concepts_nb_33_0.png

Note: the discrete function u_h already has a name attribute given in the interpolate call. This is used by default in the vtk file. An alternative name can be given by using a dictionary as shown previously.

Of course a discrete function can also be used as a coefficient in a UFL expression

[17]:
print( integrate(gridView, abs(exact-u_h), order=5) )
0.00016187237493796603

The main difference between a grid function and a discrete function is that the latter has a dofVector attached to it. We can iterate over that easily

[18]:
for d in u_h.dofVector:
    print(d,end=", ")
print()
1.0, 1.0260416666666667, 1.0833333333333333, 1.140625, 1.1666666666666667, 1.0364583333333333, 1.0625, 1.1197916666666667, 1.1770833333333333, 1.203125, 1.1666666666666667, 1.1927083333333333, 1.25, 1.3072916666666667, 1.3333333333333335, 1.421875, 1.4479166666666665, 1.5052083333333335, 1.5625, 1.5885416666666667, 1.8333333333333333, 1.859375, 1.9166666666666665, 1.9739583333333333, 2.0, 1.0071614583333333, 1.052734375, 1.1139322916666667, 1.1595052083333333, 1.0436197916666667, 1.0891927083333333, 1.150390625, 1.1959635416666667, 1.173828125, 1.2194010416666667, 1.2805989583333333, 1.326171875, 1.4290364583333333, 1.474609375, 1.5358072916666665, 1.5813802083333335, 1.8404947916666665, 1.8860677083333333, 1.947265625, 1.9928385416666665, 1.0084635416666667, 1.0345052083333333, 1.091796875, 1.1490885416666667, 1.1751302083333335, 1.087890625, 1.1139322916666667, 1.1712239583333333, 1.228515625, 1.2545572916666667, 1.2766927083333333, 1.302734375, 1.3600260416666667, 1.4173177083333333, 1.443359375, 1.6061197916666665, 1.6321614583333333, 1.689453125, 1.7467447916666665, 1.7727864583333335, 1.015625, 1.0611979166666667, 1.1223958333333333, 1.16796875, 1.0950520833333333, 1.140625, 1.2018229166666667, 1.2473958333333333, 1.2838541666666665, 1.3294270833333333, 1.390625, 1.4361979166666667, 1.61328125, 1.6588541666666665, 1.7200520833333333, 1.765625,

In a later section we will see how to extract the underlying dof vector in the form of a numpy or a petsc vector.

Models and Schemes

We consider a scalar boundary value problem \begin{align*} -\triangle u &= f & \text{in}\;\Omega:=(0,1)^2 \\ \nabla u\cdot n &= g_N & \text{on}\;\Gamma_N \\ u &= g_D & \text{on}\;\Gamma_D \end{align*} and \(f=f(x)\) is some forcing term. For the boundary conditions we set \(\Gamma_D={0}\times[0,1]\) and take \(\Gamma_N\) to be the remaining boundary of \(\Omega\).

We will solve this problem in variational form \begin{align*} \int \nabla u \cdot \nabla \varphi \ - \int_{\Omega} f(x) \varphi\ dx - \int_{\Gamma_N} g_N(x) v\ ds = 0. \end{align*} We choose \(f,g_N,g_D\) so that the exact solution is given by \begin{align*} u(x) = \left(\frac{1}{2}(x^2 + y^2) - \frac{1}{3}(x^3 - y^3)\right) + 1 \end{align*} Note: in a later section we discuss more general boundary conditions.

[19]:
from ufl import TestFunction, TrialFunction
from dune.ufl import DirichletBC
u = TrialFunction(space)
v = TestFunction(space)

from ufl import dx, grad, div, grad, dot, inner, sqrt, conditional, FacetNormal, ds
a = dot(grad(u), grad(v)) * dx

f   = -div( grad(exact) )
g_N = grad(exact)
n   = FacetNormal(space)
b   = f*v*dx + dot(g_N,n)*conditional(x[0]>=1e-8,1,0)*v*ds
dbc = DirichletBC(space,exact,x[0]<=1e-8)

With the model described as a ufl form, we can construct a scheme class that provides the solve method which we can use to compute the solution:

[20]:
from dune.fem.scheme import galerkin as solutionScheme
scheme = solutionScheme([a == b, dbc], solver='cg')
scheme.solve(target = u_h)
[20]:
{'converged': True, 'iterations': 1, 'linear_iterations': 49}

We can compute the error between the exact and the discrete solution by using the integrate function described above:

[21]:
h1error = dot(grad(u_h - exact), grad(u_h - exact))
error = sqrt(integrate(gridView, h1error, order=5))
print("Number of elements:",gridView.size(0),
      "number of dofs:",space.size,"H^1 error:", error)
Number of elements: 16 number of dofs: 81 H^1 error: 0.00658807845868413

To verify that the discrete scheme is converging to the exact solution we can compute the experimental order of convergence (EOC): \begin{align*} {\rm eoc} = \frac{\log{e_h/e_H}}{\log{h/H}} \end{align*} where \(h,H\) refer to the spacing of two grids and \(e_h,e_H\) are measures for the accuracy of the discrete solution, i.e., \(e_h=\|u-u_h\|,e_H=\|u-u_h\) using a suitable norm (here \(h_1\) as used before):

[22]:
import dune.fem
loops = 2
for eocLoop in range(loops):
    error_old = error
    # gridView.hierarchicalGrid.globalRefine(1)
    dune.fem.globalRefine(1,gridView.hierarchicalGrid)
    u_h.interpolate(0)
    scheme.solve(target = u_h)
    error = sqrt(integrate(gridView, h1error, order=5))
    eoc = round(numpy.log(error/error_old)/numpy.log(0.5),2)
    print("EOC:",eoc,
          "Number of elements:", gridView.size(0),
          "number of dofs:", space.size,"H^1 error:", error)
EOC: -0.0 Number of elements: 16 number of dofs: 81 H^1 error: 0.0065880784587151
EOC: -0.0 Number of elements: 16 number of dofs: 81 H^1 error: 0.0065880784587151

We have already seen how grid function can be used within UFL expressions and forms - in the next section we will give another example for this in the context of a time dependent problems. In addition we can also use the Constant class to add constants to UFL expressions/forms which then can be changed without requiring any recompilation of the model. An example would again be in a time dependent problem a time varying coefficient. Being able to change the value of the time in the model without recompilation is crucial for an efficient code. We will demonstrate this here by adding a mass term with a constant mass \(m\). We will not change the right hand side so the exact solution we used so far will only be valid for \(m=0\).

[23]:
from dune.ufl import Constant

m = Constant(1,name="mass")   # start with m=1
a = dot(grad(u), grad(v)) * dx + m*u*v * dx
scheme = solutionScheme([a == b, dbc], solver='cg')

scheme.solve(target = u_h)

h1error = dot(grad(u_h - exact), grad(u_h - exact))
error = sqrt(integrate(gridView, h1error, order=5))
print("Number of elements:",gridView.size(0),
      "number of dofs:",space.size,"H^1 error:", error)
Number of elements: 16 number of dofs: 81 H^1 error: 0.5635954310234561

We can print the value of a Constant with name foo via scheme.model.foo and change it’s value using the same attribute:

[24]:
print(scheme.model.mass)
scheme.model.mass = 0      # go back to original problem
print(scheme.model.mass)

scheme.solve(target = u_h)

h1error = dot(grad(u_h - exact), grad(u_h - exact))
error = sqrt(integrate(gridView, h1error, order=5))
print("Number of elements:",gridView.size(0),
      "number of dofs:",space.size,"H^1 error:", error)
1.0
0.0
Number of elements: 16 number of dofs: 81 H^1 error: 0.0065880784586842085

Parallelization

Parallelization is available using either distributed memory based on MPI or multithreading using OpenMP.

MPI

It is straightforward to use MPI for parallelization. It requires a parallel grid, in fact most of the DUNE grids work in parallel except albertaGrid or polyGrid. Most iterative solvers in DUNE work for parallel runs. Some of the preconditioning methods also work in parallel. Running a parallel job can be done by using mpirun

mpirun -np 4 python script.py

in order to use 4 MPI processes. Example scripts that run in parallel are, for example, the Re-entrant Corner Problem.

OpenMP

It is straightforward to enable some multithreading support. Note that this will speedup assembly and evaluation of the spatial operators but not in general the linear algebra backend so that the overall speedup of the code might not be as high as expected. Since we rely mostly on external packages for the linear algebra, speedup in this step will depend on the multithreading support available in the chosen linear algebra backend. We discuss in a later section how to switch between different implementation for solving linear systems of equations, e.g., how to use the thread parallel solvers from scipy.

By default only a single thread is used. To enable multithreading simply add

[25]:
from dune.fem import threading
print("Using",threading.use,"threads")
threading.use = 4 # use 4 threads
print("Using",threading.use,"threads")
Using 1 threads
Using 4 threads

At startup a maximum number of threads is selected based on the hardware concurrency. This number can be changed by setting the environment variable DUNE_NUM_THREADS which sets both the maximum and set the number of threads to use. To get this number or set the number of threads to use to the maximum use

[26]:
print("Maximum number of threads available:",threading.max)
threading.useMax()
print("Using",threading.use,"threads")
Maximum number of threads available: 12
Using 12 threads

Listing Available Dune Components

The available realization of a given interface, i.e., the available grid implementations, depends on the modules found during configuration. Getting access to all available components is straightforward:

[27]:
from dune.utility import components
# to get a list of all available components:
components()
# to get for example all available grid implementations:
components("grid")
available categories are:
 discretefunction,function,grid,model,operator,scheme,solver,space,view
available entries for this category are:
entry         function          module
------------------------------------------
agglomerate   polyGrid          dune.vem.vem
alberta       albertaGrid       dune.grid
alu           aluGrid           dune.alugrid
aluconform    aluConformGrid    dune.alugrid
alucube       aluCubeGrid       dune.alugrid
alusimplex    aluSimplexGrid    dune.alugrid
oned          onedGrid          dune.grid
sp            spIsotropicGrid   dune.spgrid
spanisotropic spAnisotropicGrid dune.spgrid
spbisection   spBisectionGrid   dune.spgrid
spisotropic   spIsotropicGrid   dune.spgrid
ug            ugGrid            dune.grid
yasp          yaspGrid          dune.grid
------------------------------------------
[ ]:

Available discrete function spaces are:

[28]:
components("space")
available entries for this category are:
entry          function       module
------------------------------------------
bbdg           bbdgSpace      dune.vem.vem
bdm            bdm            dune.fem.space
combined       combined       dune.fem.space
dglagrange     dglagrange     dune.fem.space
dglegendre     dglegendre     dune.fem.space
dglegendrehp   dglegendrehp   dune.fem.space
dgonb          dgonb          dune.fem.space
dgonbhp        dgonbhp        dune.fem.space
finitevolume   finiteVolume   dune.fem.space
lagrange       lagrange       dune.fem.space
lagrangehp     lagrangehp     dune.fem.space
p1bubble       p1Bubble       dune.fem.space
product        product        dune.fem.space
rannacherturek rannacherTurek dune.fem.space
raviartthomas  raviartThomas  dune.fem.space
vem            vemSpace       dune.vem.vem
------------------------------------------

Available discrete functions are:

[29]:
components("function")
available entries for this category are:
entry      function          module
--------------------------------------------
cpp        cppFunction       dune.fem.function
discrete   discreteFunction  dune.fem.function
global     globalFunction    dune.fem.function
levels     levelFunction     dune.fem.function
local      localFunction     dune.fem.function
partitions partitionFunction dune.fem.function
ufl        uflFunction       dune.fem.function
--------------------------------------------

Available schemes:

[30]:
components("scheme")
available entries for this category are:
entry      function         module
-----------------------------------------
bbdg       bbdgScheme       dune.vem.vem
burgers    burgers          dune.fem.scheme
dg         dg               dune.fem.scheme
dggalerkin dgGalerkin       dune.fem.scheme
galerkin   galerkin         dune.fem.scheme
h1         h1               dune.fem.scheme
h1galerkin h1Galerkin       dune.fem.scheme
linearized linearized       dune.fem.scheme
rungekutta rungeKuttaSolver dune.femdg
stokes     stokes           dune.fem.scheme
vem        vemScheme        dune.vem.vem
-----------------------------------------

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