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

\[-\triangle u + \alpha u = f\]

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
-----------------------------------------

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