Run this tutorial

Click here to run this tutorial on mybinder.org: try on mybinder.org
Please note that starting the notebook server may take a couple of minutes.

Example: prolongations and computing products and norms

This is work in progress (WIP), still missing:

  • documentation and explanation

  • assembled matrix as product

# wurlitzer: display dune's output in the notebook
%load_ext wurlitzer


import numpy as np

1: prolongation onto a finer grid

Let’s set up a 2d grid first, as seen in other tutorials and examples.

from dune.xt.grid import Dim, Cube, Simplex, make_cube_grid
from dune.xt.functions import ExpressionFunction, GridFunction as GF

d = 2
omega = ([0, 0], [1, 1])
grid = make_cube_grid(Dim(d), Simplex(), lower_left=omega[0], upper_right=omega[1], num_elements=[2, 2])
grid.global_refine(1)

print(f'grid has {grid.size(0)} elements, {grid.size(d - 1)} edges and {grid.size(d)} vertices')

grid has 16 elements, 28 edges and 13 vertices

from dune.gdt import default_interpolation, prolong, DiscontinuousLagrangeSpace

V_H = DiscontinuousLagrangeSpace(grid, order=1)

f = ExpressionFunction(dim_domain=Dim(d), variable='x', order=10, expression='exp(x[0]*x[1])', name='f')
f_H = default_interpolation(GF(grid, f), V_H)

print(f'f_H has {len(f_H.dofs.vector)} DoFs')

f_H has 48 DoFs

reference_grid = make_cube_grid(Dim(d), Simplex(), lower_left=omega[0], upper_right=omega[1], num_elements=[16, 16])
reference_grid.global_refine(1)

print(f'grid has {reference_grid.size(0)} elements')

grid has 1024 elements

from dune.gdt import DiscreteFunction

V_h = DiscontinuousLagrangeSpace(reference_grid, order=1)

f_h = prolong(f_H, V_h)
print(f'f_h has {len(f_h.dofs.vector)} DoFs')

f_h has 3072 DoFs

2: computing products and norms

u = GF(grid,
       ExpressionFunction(dim_domain=Dim(d), variable='x', order=10, expression='exp(x[0]*x[1])',
                          gradient_expressions=['x[1]*exp(x[0]*x[1])', 'x[0]*exp(x[0]*x[1])'], name='h'))
\[ (u, v)_{H^1} = (u, v)_{L^2} + (u, v)_{H^1_0} \]
\[ ||u||_{H^1} = \sqrt{||u||_{L^2}^2 + |u|_{H^1_0}^2} \]

We collect the local integrands of each product in a BilinearForm (via +=) and evaluate a(u, u) with apply2(u, u), which walks the grid and returns the scalar:

from dune.gdt import (
    BilinearForm,
    LocalElementProductIntegrand,
    LocalElementIntegralBilinearForm,
    LocalLaplaceIntegrand,
)

L2_product = BilinearForm(grid)
L2_product += LocalElementIntegralBilinearForm(LocalElementProductIntegrand(GF(grid, 1)))

H1_semi_product = BilinearForm(grid)
H1_semi_product += LocalElementIntegralBilinearForm(LocalLaplaceIntegrand(GF(grid, 1, dim_range=(Dim(d), Dim(d)))))

H1_product = BilinearForm(grid)
H1_product += LocalElementIntegralBilinearForm(LocalElementProductIntegrand(GF(grid, 1)))
H1_product += LocalElementIntegralBilinearForm(LocalLaplaceIntegrand(GF(grid, 1, dim_range=(Dim(d), Dim(d)))))

l2 = L2_product.apply2(u, u)
h1_semi = H1_semi_product.apply2(u, u)
h1 = H1_product.apply2(u, u)

print(f'|| u ||_L2 = {np.sqrt(l2)}')
print(f' | u |_H1  = {np.sqrt(h1_semi)}')
print(f'|| u ||_H1 = {np.sqrt(h1)}')

|| u ||_L2 = 1.357179337917508
 | u |_H1  = 1.2638291121558571
|| u ||_H1 = 1.8545079616984306

np.allclose(l2 + h1_semi, h1)

True

Equivalently, norm returns the induced norm directly (norm(u) == sqrt(apply2(u, u))):

print(f'|| u ||_L2 = {L2_product.norm(u)}')
print(f' | u |_H1  = {H1_semi_product.norm(u)}')
print(f'|| u ||_H1 = {H1_product.norm(u)}')

|| u ||_L2 = 1.357179337917508
 | u |_H1  = 1.2638291121558571
|| u ||_H1 = 1.8545079616984306

3: computing errors w.r.t. a known reference solution

# suppose we have a DiscreteFunction on a coarse grid, usually the solution of a discretization
# (we use the interpolation here for simplicity)

u_H = default_interpolation(u, V_H)

# prolong it onto the reference grid
u_h = prolong(u_H, V_h)

# define error function using the analytically known solution u
error = GF(grid, u - u_h)

# compute norms on reference grid
L2_product = BilinearForm(reference_grid)
L2_product += LocalElementIntegralBilinearForm(LocalElementProductIntegrand(GF(grid, 1)))
print(f'|| error ||_L2 = {L2_product.norm(error)}')

|| error ||_L2 = 0.02215546188074599

Download the code: example__prolongations_products_and_norms.md example__prolongations_products_and_norms.ipynb