General Concepts (with a time dependent 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 (hierarchical grid and grid view)
working with functions defined over the grid (grid functions)
setting up discrete spaces and discrete functions
defining the mathematical model to solve (ufl)
solving the (non linear) system arising from the discretization of the model by the Galerkin method (scheme)
Setting up the Grid
Note
The goal of a Dune grid or hierarchical grid is to provide an interface to describe subdivisions (entities or elements for codim 0) of the computational domain \(\Omega\) in a generic way. In the following this is referred to as hierarchicalGrid
. See [BBD+21], Section 3.1, for more details.
After some general import statements we start our tutorial by setting up a simple Cartesian grid over the domain \(\Omega = [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 dictionary.
[1]:
from matplotlib import pyplot as plt
import numpy as np
from dune.grid import structuredGrid as leafGridView
gridView = leafGridView([0, 0], [1, 1], [4, 4])
gridView.plot()

Often the initial grid is a bit coarse - it is very easy to refine it globally. Grids can be refined in a hierarchic manner, meaning that elements are subdivided into several smaller elements. The element to be refined is kept in the grid and remains accessible. We will discuss local refinement in a later section.
Note
A grid view object provides read-only access to the entities of all codimensions in the view, i.e. a subset of the hierarchical grid. In the following this is called gridView
and sometimes abbreviated as grid
. See [BBD+21], Section 3.1.1, for more details.
Note that the gridView
we have constructed can not be changed directly (i.e. read-only access), we need to access the underlying hierarchicalGrid
:
[2]:
from dune.fem import globalRefine
globalRefine(1,gridView.hierarchicalGrid) # the first argument is the number of refinements to perform
gridView.plot()

Grid Functions
Note
A grid functions is a functions that is defined over a given grid view and is evaluated by using an element of this view and local coordinate within that element, i.e. a quadrature point.
For example:
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:
[3]:
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.3333333333333324
and plot them using matplotlib or write a vtk file for postprocessing
[4]:
from dune.fem.plotting import plotPointData as plot
plot(exact, grid=gridView)
gridView.writeVTK('exact', pointdata={'exact': exact})

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).
[5]:
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.33203125
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:
[6]:
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:
[7]:
@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)
1.33203125
As pointed out the dune.fem
grid function can be used like any other UFL coefficient to form UFL expressions:
[8]:
print( integrate(gridView, abs(exact-exactLocal), order=5) )
gf = uflFunction(gridView, name="ufl", order=1, ufl=exact+exactLocal*exactGlobal)
fig = plt.figure(figsize=(20,10))
gf.plot(figure=(fig,121))
exactGlobal.plot(figure=(fig,122))
0.0

Converting UFL expressions to grid functions leads to JIT code generation and is therefore efficient when used in other C++ algorithm (like integrate
). But more complex control structures are hard to express within UFL. 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 as described in the section on
generating C++ Based Grid Functions.
Discrete Spaces and Discrete Functions
The grid functions set up so far did not involve any discretization, they are exactly evaluated at the given point.
Note
A discrete (function) space is a finite dimensional function space defined over a given grid view, typically the Finite Element space involved.
[9]:
from dune.fem.space import lagrange as solutionSpace
space = solutionSpace(gridView, order=2)
A list of available spaces is shown at the bottom of this page.
Note
A special type of grid function is a discrete function living in a discrete (finite dimensional) space.
The easiest way to construct a discrete function is to use the interpolate
method on the discrete function space to create an initialized discrete function or function
to create an uninitialized discrete function.
[10]:
# initialized
u_h = space.interpolate(exact, name='u_h')
# uninitialized
u_h = space.function(name='u_h')
On an existing discrete function the interpolate
method can be used to reinitialize it
[11]:
u_h.interpolate( lambda x: 1/2.*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1 )
u_h.plot()

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 to be carried out.
If a discrete function is already available it is possible to call copy
to obtain further discrete functions:
[12]:
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:
[13]:
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
[14]:
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.33203125
or plot using matplotlib and write a vtk file for postprocessing (using binary data format to reduce size)
[15]:
u_h.plot(gridLines="white")
from dune.grid import OutputType
gridView.writeVTK('uh', pointdata=[u_h])

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
[16]:
print( integrate(gridView, abs(exact-u_h), order=5) )
2.0234046867290766e-05
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
[17]:
for d in u_h.dofVector:
print(d,end=", ")
print()
1.0, 1.0071614583333333, 1.0260416666666667, 1.052734375, 1.0833333333333333, 1.1139322916666667, 1.140625, 1.1595052083333333, 1.1666666666666667, 1.0084635416666667, 1.015625, 1.0345052083333333, 1.0611979166666667, 1.091796875, 1.1223958333333333, 1.1490885416666667, 1.16796875, 1.1751302083333335, 1.0364583333333333, 1.0436197916666667, 1.0625, 1.0891927083333333, 1.1197916666666667, 1.150390625, 1.1770833333333333, 1.1959635416666667, 1.203125, 1.087890625, 1.0950520833333333, 1.1139322916666667, 1.140625, 1.1712239583333333, 1.2018229166666667, 1.228515625, 1.2473958333333333, 1.2545572916666667, 1.1666666666666667, 1.173828125, 1.1927083333333333, 1.2194010416666667, 1.25, 1.2805989583333333, 1.3072916666666667, 1.326171875, 1.3333333333333335, 1.2766927083333333, 1.2838541666666665, 1.302734375, 1.3294270833333333, 1.3600260416666667, 1.390625, 1.4173177083333333, 1.4361979166666667, 1.443359375, 1.421875, 1.4290364583333333, 1.4479166666666665, 1.474609375, 1.5052083333333335, 1.5358072916666665, 1.5625, 1.5813802083333335, 1.5885416666666667, 1.6061197916666665, 1.61328125, 1.6321614583333333, 1.6588541666666665, 1.689453125, 1.7200520833333333, 1.7467447916666665, 1.765625, 1.7727864583333335, 1.8333333333333333, 1.8404947916666665, 1.859375, 1.8860677083333333, 1.9166666666666665, 1.947265625, 1.9739583333333333, 1.9928385416666665, 2.0, 1.0018717447916667, 1.015380859375, 1.0386555989583333, 1.0677897135416667, 1.098876953125, 1.1280110677083333, 1.1512858072916667, 1.164794921875, 1.0103352864583333, 1.0238444010416667, 1.047119140625, 1.0762532552083333, 1.1073404947916667, 1.136474609375, 1.1597493489583333, 1.1732584635416667, 1.038330078125, 1.0518391927083333, 1.0751139322916667, 1.104248046875, 1.1353352864583333, 1.1644694010416667, 1.187744140625, 1.2012532552083335, 1.0897623697916667, 1.103271484375, 1.1265462239583333, 1.1556803385416667, 1.186767578125, 1.2159016927083333, 1.2391764322916667, 1.252685546875, 1.1685384114583333, 1.1820475260416667, 1.205322265625, 1.2344563802083333, 1.2655436197916667, 1.294677734375, 1.3179524739583335, 1.3314615885416667, 1.278564453125, 1.2920735677083333, 1.3153483072916665, 1.344482421875, 1.3755696614583333, 1.4047037760416667, 1.427978515625, 1.4414876302083335, 1.4237467447916665, 1.437255859375, 1.4605305989583333, 1.4896647135416665, 1.520751953125, 1.5498860677083335, 1.5731608072916665, 1.586669921875, 1.6079915364583333, 1.6215006510416665, 1.644775390625, 1.6739095052083333, 1.7049967447916665, 1.734130859375, 1.7574055989583335, 1.7709147135416665, 1.835205078125, 1.8487141927083333, 1.8719889322916665, 1.901123046875, 1.9322102864583333, 1.9613444010416665, 1.984619140625, 1.9981282552083335, 1.0020345052083333, 1.0091959635416667, 1.028076171875, 1.0547688802083333, 1.0853678385416667, 1.115966796875, 1.1426595052083333, 1.1615397135416667, 1.168701171875, 1.019775390625, 1.0269368489583333, 1.0458170572916667, 1.072509765625, 1.1031087239583333, 1.1337076822916667, 1.160400390625, 1.1792805989583333, 1.1864420572916667, 1.0590006510416667, 1.066162109375, 1.0850423177083333, 1.1117350260416667, 1.142333984375, 1.1729329427083333, 1.1996256510416667, 1.218505859375, 1.2256673177083335, 1.1236165364583333, 1.1307779947916667, 1.149658203125, 1.1763509114583333, 1.2069498697916667, 1.237548828125, 1.2642415364583335, 1.2831217447916667, 1.290283203125, 1.217529296875, 1.2246907552083333, 1.2435709635416667, 1.270263671875, 1.3008626302083333, 1.3314615885416667, 1.358154296875, 1.3770345052083335, 1.3841959635416667, 1.3446451822916665, 1.351806640625, 1.3706868489583333, 1.3973795572916665, 1.427978515625, 1.4585774739583333, 1.4852701822916667, 1.504150390625, 1.5113118489583335, 1.5088704427083333, 1.5160319010416665, 1.534912109375, 1.5616048177083333, 1.5922037760416665, 1.622802734375, 1.6494954427083335, 1.6683756510416665, 1.675537109375, 1.714111328125, 1.7212727864583333, 1.7401529947916665, 1.766845703125, 1.7974446614583333, 1.8280436197916665, 1.854736328125, 1.8736165364583335, 1.8807779947916665, 1.00390625, 1.0174153645833333, 1.0406901041666667, 1.06982421875, 1.1009114583333333, 1.1300455729166667, 1.1533203125, 1.1668294270833335, 1.0216471354166667, 1.03515625, 1.0584309895833333, 1.0875651041666667, 1.11865234375, 1.1477864583333333, 1.1710611979166667, 1.1845703125, 1.0608723958333333, 1.0743815104166667, 1.09765625, 1.1267903645833333, 1.1578776041666667, 1.18701171875, 1.2102864583333333, 1.2237955729166667, 1.12548828125, 1.1389973958333333, 1.1622721354166667, 1.19140625, 1.2224934895833333, 1.2516276041666667, 1.27490234375, 1.2884114583333335, 1.2194010416666667, 1.23291015625, 1.2561848958333333, 1.2853190104166667, 1.31640625, 1.3455403645833333, 1.3688151041666667, 1.38232421875, 1.3465169270833333, 1.3600260416666665, 1.38330078125, 1.4124348958333333, 1.4435221354166667, 1.47265625, 1.4959309895833335, 1.5094401041666667, 1.5107421875, 1.5242513020833333, 1.5475260416666665, 1.57666015625, 1.6077473958333335, 1.6368815104166665, 1.66015625, 1.6736653645833335, 1.7159830729166665, 1.7294921875, 1.7527669270833333, 1.7819010416666665, 1.81298828125, 1.8421223958333333, 1.8653971354166665, 1.87890625,
In addition to clear
and assign
we can also add (scalar multiples of) discrete functions together or multiply them by a scalar which will directly change the underlying dofVector - this can only be done in place, so u_h *= 2
changes the dof vector while u_h*2
results in an UFL expression. Here are some examples:
[18]:
u_h_n *= 2
u_h += u_h_n
# this is equivalent to
u_h.axpy(2,u_h_n)
We started with u_h_n=u_h
so after these changes each dof in u_h
should be seven times its original value
[19]:
for d in u_h.dofVector:
print(d,end=", ")
print()
7.0, 7.050130208333333, 7.182291666666667, 7.369140625, 7.583333333333333, 7.797526041666667, 7.984375, 8.116536458333332, 8.166666666666668, 7.059244791666667, 7.109375, 7.241536458333333, 7.428385416666667, 7.642578125, 7.856770833333333, 8.043619791666668, 8.17578125, 8.225911458333334, 7.255208333333333, 7.305338541666667, 7.4375, 7.624348958333333, 7.838541666666667, 8.052734375, 8.239583333333332, 8.371744791666668, 8.421875, 7.615234375, 7.665364583333333, 7.797526041666667, 7.984375, 8.198567708333332, 8.412760416666668, 8.599609375, 8.731770833333332, 8.781901041666668, 8.166666666666668, 8.216796875, 8.348958333333332, 8.535807291666668, 8.75, 8.964192708333332, 9.151041666666668, 9.283203125, 9.333333333333334, 8.936848958333332, 8.986979166666666, 9.119140625, 9.305989583333332, 9.520182291666668, 9.734375, 9.921223958333332, 10.053385416666668, 10.103515625, 9.953125, 10.003255208333332, 10.135416666666666, 10.322265625, 10.536458333333334, 10.750651041666666, 10.9375, 11.069661458333334, 11.119791666666668, 11.242838541666666, 11.29296875, 11.425130208333332, 11.611979166666666, 11.826171875, 12.040364583333332, 12.227213541666666, 12.359375, 12.409505208333334, 12.833333333333332, 12.883463541666666, 13.015625, 13.202473958333332, 13.416666666666666, 13.630859375, 13.817708333333332, 13.949869791666666, 14.0, 7.013102213541667, 7.107666015625, 7.270589192708333, 7.474527994791667, 7.692138671875, 7.896077473958333, 8.059000651041668, 8.153564453125, 7.072347005208333, 7.166910807291667, 7.329833984375, 7.533772786458333, 7.751383463541667, 7.955322265625, 8.118245442708332, 8.212809244791668, 7.268310546875, 7.362874348958333, 7.525797526041667, 7.729736328125, 7.947347005208333, 8.151285807291668, 8.314208984375, 8.408772786458334, 7.628336588541667, 7.722900390625, 7.885823567708333, 8.089762369791668, 8.307373046875, 8.511311848958332, 8.674235026041668, 8.768798828125, 8.179768880208332, 8.274332682291668, 8.437255859375, 8.641194661458332, 8.858805338541668, 9.062744140625, 9.225667317708334, 9.320231119791668, 8.949951171875, 9.044514973958332, 9.207438151041666, 9.411376953125, 9.628987630208332, 9.832926432291668, 9.995849609375, 10.090413411458334, 9.966227213541666, 10.060791015625, 10.223714192708332, 10.427652994791666, 10.645263671875, 10.849202473958334, 11.012125651041666, 11.106689453125, 11.255940755208332, 11.350504557291666, 11.513427734375, 11.717366536458332, 11.934977213541666, 12.138916015625, 12.301839192708334, 12.396402994791666, 12.846435546875, 12.940999348958332, 13.103922526041666, 13.307861328125, 13.525472005208332, 13.729410807291666, 13.892333984375, 13.986897786458334, 7.014241536458333, 7.064371744791667, 7.196533203125, 7.383382161458333, 7.597574869791667, 7.811767578125, 7.998616536458333, 8.130777994791668, 8.180908203125, 7.138427734375, 7.188557942708333, 7.320719401041667, 7.507568359375, 7.721761067708333, 7.935953776041667, 8.122802734375, 8.254964192708332, 8.305094401041668, 7.413004557291667, 7.463134765625, 7.595296223958333, 7.782145182291667, 7.996337890625, 8.210530598958332, 8.397379557291668, 8.529541015625, 8.579671223958334, 7.865315755208333, 7.915445963541667, 8.047607421875, 8.234456380208332, 8.448649088541668, 8.662841796875, 8.849690755208334, 8.981852213541668, 9.031982421875, 8.522705078125, 8.572835286458332, 8.704996744791668, 8.891845703125, 9.106038411458332, 9.320231119791668, 9.507080078125, 9.639241536458334, 9.689371744791668, 9.412516276041666, 9.462646484375, 9.594807942708332, 9.781656901041666, 9.995849609375, 10.210042317708332, 10.396891276041668, 10.529052734375, 10.579182942708334, 10.562093098958332, 10.612223307291666, 10.744384765625, 10.931233723958332, 11.145426432291666, 11.359619140625, 11.546468098958334, 11.678629557291666, 11.728759765625, 11.998779296875, 12.048909505208332, 12.181070963541666, 12.367919921875, 12.582112630208332, 12.796305338541666, 12.983154296875, 13.115315755208334, 13.165445963541666, 7.02734375, 7.121907552083333, 7.284830729166667, 7.48876953125, 7.706380208333333, 7.910319010416667, 8.0732421875, 8.167805989583334, 7.151529947916667, 7.24609375, 7.409016927083333, 7.612955729166667, 7.83056640625, 8.034505208333332, 8.197428385416668, 8.2919921875, 7.426106770833333, 7.520670572916667, 7.68359375, 7.887532552083333, 8.105143229166668, 8.30908203125, 8.472005208333332, 8.566569010416668, 7.87841796875, 7.972981770833333, 8.135904947916668, 8.33984375, 8.557454427083332, 8.761393229166668, 8.92431640625, 9.018880208333334, 8.535807291666668, 8.63037109375, 8.793294270833332, 8.997233072916668, 9.21484375, 9.418782552083332, 9.581705729166668, 9.67626953125, 9.425618489583332, 9.520182291666666, 9.68310546875, 9.887044270833332, 10.104654947916668, 10.30859375, 10.471516927083334, 10.566080729166668, 10.5751953125, 10.669759114583332, 10.832682291666666, 11.03662109375, 11.254231770833334, 11.458170572916666, 11.62109375, 11.715657552083334, 12.011881510416666, 12.1064453125, 12.269368489583332, 12.473307291666666, 12.69091796875, 12.894856770833332, 13.057779947916666, 13.15234375,
In a later section we will see how to extract the underlying dof vector in the form of a numpy
or a petsc
vector.
Handling constant parameters
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 considering
where \(\alpha\) is a real valued parameter. We choose \(f,g_N,g_D\) so that the exact solution for \(\alpha=0\) 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*} which is the grid function exact
we have already defined above.
Although the problem is linear and we could just use dune.fem.assemble
and Scipy to solve, we use a scheme
instead.
Note
A scheme implements the weak form a specific type of equation, where the mathematical model (integrands) is provided as ufl form.
[20]:
from ufl import (TestFunction, TrialFunction, SpatialCoordinate, FacetNormal,
dx, ds, grad, div, inner, dot, sin, cos, pi, conditional)
from dune.ufl import Constant, DirichletBC
from dune.fem.scheme import galerkin as solutionScheme
alpha = Constant(1,name="mass") # start with alpha=1
x = SpatialCoordinate(space)
u = TrialFunction(space)
v = TestFunction(space)
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
a = dot(grad(u), grad(v)) * dx + alpha*u*v * dx
dbc = DirichletBC(space,exact,x[0]<=1e-8)
scheme = solutionScheme([a == b, dbc], solver='cg')
scheme.solve(target = u_h)
h1error = dot(grad(u_h - exact), grad(u_h - exact))
error = np.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: 64 number of dofs: 289 H^1 error: 0.5635700799029998
Note that \(\alpha=1\) so we do not have the correct exact solution. We can print the value of a Constant
with name foo
via scheme.model.foo
and change it’s value using the same attribute:
[21]:
print(scheme.model.mass)
scheme.model.mass = 0 # change to the problem which matches our exact solution
print(scheme.model.mass)
scheme.solve(target = u_h)
h1error = dot(grad(u_h - exact), grad(u_h - exact))
error = np.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: 64 number of dofs: 289 H^1 error: 0.0016470196146885557
A Nonlinear Time Dependent Problem
As a second example we will study the Forchheimer problem [Kie15] which is a scalar, nonlinear parabolic equation \begin{equation} \partial_tu - \nabla\cdot K(\nabla u)\nabla u = f \end{equation} where the diffusion tensor is given by \begin{equation} K(\nabla u) = \frac{2}{1+\sqrt{1+4|\nabla u|}} \end{equation} and \(f=f(x,t)\) is some time dependent forcing term. On the boundary we prescribe Neumann boundary as before and initial conditions \(u=u_0\).
We will solve this problem in variational form and using Crank Nicolson in time
\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*} We solve this problem using a quadratic finite element space. We first setup this space, define the initial data, and a second discrete function to store the solution at the current time step (u_h
will be used to store the next time step):
[22]:
space = solutionSpace(gridView, order=2)
u_h = space.interpolate(0, name='u_h')
u_h_n = u_h.copy(name="previous")
x = SpatialCoordinate(space)
initial = 1/2*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1
Next we setup the form \(a=a(u,v)\) which depends on time \(t\) and the chosen time step \(\Delta t\). Since we want to be able to change these during the simulation (we need to do this for \(t\) and might want to do this for \(\Delta t\) we use the dune.ufl.Constant
which stores a mutable scalar and can be used within a ufl expression:
[23]:
from ufl import FacetNormal, ds, div, dot, sqrt, exp
from dune.ufl import Constant
dt = Constant(0, name="dt") # time step
t = Constant(0, name="t") # current time
u = TrialFunction(space)
v = TestFunction(space)
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
To define the right hand side we use the given function \(u\) to define the forcing \(f\) and the boundary flux \(g\)
[24]:
exact = lambda t: exp(-2*t)*(initial - 1) + 1
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
With the model described as a ufl form, we can again construct a scheme class that provides the solve method which we can use to evolve the solution from one time step to the next. As already mentioned the solve methods always assumes that the problem is non-linear and use a Newton method to solve the problem:
[25]:
scheme = solutionScheme(a == b, solver='cg')
Optionally, we may want to increase the default quadrature orders which are 2 * space.order
for element integrals and 2 * space.order + 1
for surface integrals. Depending on the data this might not be enough. Then we simply set the integration orders by hand like in the following example, by calling the method setQuadratureOrders(interiorOrder, surfaceOrder).
[26]:
#scheme.setQuadratureOrders( 2*space.order, 2*space.order+1 )
Since we have forced the system towards a given solution, we can compute the discretization error. First we define ufl expressions for the \(L^2\) and \(H^1\) norms and will use those to compute the experimental order of convergence of the scheme by computing the time evolution on different grid levels. We use dune.fem.function.uflFunction
to generate code for the errors which can be more efficient than using the ufl expressions directly.
[27]:
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)))
Now we evolve the solution from time \(t=0\) to \(t=T\) in a loop. Since the PDE model has time dependent coefficient (through the forcing term), we need to update the t
constant used to define the model before each step. A second constant we used to define the model was dt
which defines the time step. We keep this constant and set it to \(0.002\) at the beginning of the simulation. This value could be changed in each time step:
[28]:
scheme.model.dt = 0.01
time = 0
u_h.interpolate(initial)
while time < (endTime - 1e-6):
scheme.model.t = time
u_h_n.assign(u_h)
scheme.solve(target=u_h)
time += scheme.model.dt
errors = [np.sqrt(e) for e in integrate(gridView, [l2error,h1error], order=5)]
print('grid size:', gridView.size(0))
print('\t | u_h - u | =', '{:0.5e}'.format(errors[0]))
print('\t | grad(uh - u) | =', '{:0.5e}'.format(errors[1]))
u_h.plot()
gridView.writeVTK('forchheimer', pointdata={'u': u_h, 'l2error': l2error, 'h1error': h1error})
grid size: 64
| u_h - u | = 1.67785e-05
| grad(uh - u) | = 9.99185e-04

We can refine the grid once and recompute the solution on the finer grid.
[29]:
globalRefine(1, gridView.hierarchicalGrid)
scheme.model.dt /= 2
u_h.interpolate(initial)
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
errorsFine = [np.sqrt(e) for e in integrate(gridView, [l2error,h1error], order=5)]
print('grid size:', gridView.size(0))
print('\t | u_h - u | =', '{:0.5e}'.format(errorsFine[0]))
print('\t | grad(uh - u) | =', '{:0.5e}'.format(errorsFine[1]))
u_h.plot()
grid size: 256
| u_h - u | = 2.35195e-06
| grad(uh - u) | = 2.49777e-04

To check that everything is working as expected let’s compute the experimental order of convergence (convergence should be cubic for \(L^2\)-error and quadratic for the \(H^1\)-error:
[30]:
eocs = [ round(np.log(fine/coarse)/np.log(0.5),2)
for fine,coarse in zip(errorsFine,errors) ]
print("EOCs:",eocs)
EOCs: [2.83, 2.0]
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:
[31]:
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:
[32]:
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:
[33]:
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:
[34]:
components("scheme")
available entries for this category are:
entry function module
-----------------------------------------
bbdg bbdgScheme dune.vem.vem
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
vem vemScheme dune.vem.vem
-----------------------------------------