Run this tutorial

Please note that starting the notebook server may take a couple of minutes.

# Tutorial: discontinuous IPDG for the stationary heat equation#

This tutorial shows how to solve the stationary heat equation with homogeneous Dirichlet boundary conditions using interior penalty (IP) discontinuous Galerkin (DG) Finite Elements with dune-gdt.

This is work in progress [WIP], still missing:

• mathematical theory on IPDG methods

• explanation of the IPDG implementation

• non-homonegenous Dirichlet boundary values

• Neumann boundary values

• Robin boundary values

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

import numpy as np
np.warnings.filterwarnings('ignore') # silence numpys warnings

from dune.xt.grid import Dim
from dune.xt.functions import ConstantFunction, ExpressionFunction

d = 2
omega = ([0, 0], [1, 1])

kappa = ConstantFunction(dim_domain=Dim(d), dim_range=Dim(1), value=[1.], name='kappa')
# note that we need to prescribe the approximation order, which determines the quadrature on each element
f = ExpressionFunction(dim_domain=Dim(d), variable='x', expression='exp(x[0]*x[1])', order=3, name='f')

from dune.xt.grid import Simplex, make_cube_grid, visualize_grid

grid = make_cube_grid(Dim(d), Simplex(), lower_left=omega[0], upper_right=omega[1], num_elements=[2, 2])
grid.global_refine(1) # we need to refine once to obtain a symmetric grid

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

_ = visualize_grid(grid)

grid has 16 elements, 28 edges and 13 vertices

GridParameterBlock: Parameter 'refinementedge' not specified, defaulting to 'ARBITRARY'.


## 1.9: everything in a single function#

For a better overview, the above discretization code is also available in a single function in the file discretize_elliptic_ipdg.py.

import inspect
from discretize_elliptic_ipdg import discretize_elliptic_ipdg_dirichlet_zero

print(inspect.getsource(discretize_elliptic_ipdg_dirichlet_zero))

def discretize_elliptic_ipdg_dirichlet_zero(
grid, diffusion, source, symmetry_factor=1, penalty_parameter=None, weight=1
):

"""
Discretizes the stationary heat equation with homogeneous Dirichlet boundary values everywhere
with dune-gdt using an interior penalty (IP) discontinuous Galerkin (DG) method based on first
order Lagrange finite elements.

The type of IPDG scheme is determined by symmetry_factor and weight:

* with weight==1 we obtain

- symmetry_factor==-1: non-symmetric interior penalty scheme (NIPDG)
- symmetry_factor==0: incomplete interior penalty scheme (IIPDG)
- symmetry_factor==1: symmetric interior penalty scheme (SIPDG)

* with weight=diffusion, we obtain

- symmetry_factor==1: symmetric weighted interior penalty scheme (SWIPDG)

Parameters
----------
grid
The grid, given as a GridProvider from dune.xt.grid.
diffusion
Diffusion function with values in R^{d x d}, anything that dune.xt.functions.GridFunction
can handle.
source
Right hand side source function with values in R, anything that
dune.xt.functions.GridFunction can handle.
symmetry_factor
Usually one of -1, 0, 1, determines the IPDG scheme (see above).
penalty_parameter
Positive number to ensure coercivity of the resulting diffusion bilinear form,
e.g. 16 for simplicial non-degenerate grids in 2d and diffusion==1.
Determined automatically if None.
weight
Usually 1 or diffusion, determines the IPDG scheme (see above).

Returns
-------
u_h
The computed solution as a dune.gdt.DiscreteFunction for postprocessing.
"""

d = grid.dimension
diffusion = GF(grid, diffusion, dim_range=(Dim(d), Dim(d)))
source = GF(grid, source)
weight = GF(grid, weight, dim_range=(Dim(d), Dim(d)))

boundary_info = AllDirichletBoundaryInfo(grid)

V_h = DiscontinuousLagrangeSpace(grid, order=1)
if not penalty_parameter:
# TODO: check if we need to include diffusion for the coercivity here!
# TODO: each is a grid walk, compute this in one grid walk with the sparsity pattern
C_G = estimate_element_to_intersection_equivalence_constant(grid)
C_M_times_1_plus_C_T = estimate_combined_inverse_trace_inequality_constant(V_h)
penalty_parameter = C_G * C_M_times_1_plus_C_T
if symmetry_factor == 1:
penalty_parameter *= 4
assert isinstance(penalty_parameter, Number)
assert penalty_parameter > 0

l_h = VectorFunctional(grid, source_space=V_h)
l_h += LocalElementIntegralFunctional(
LocalElementProductIntegrand(GF(grid, 1)).with_ansatz(source)
)

a_h = MatrixOperator(
grid,
source_space=V_h,
range_space=V_h,
sparsity_pattern=make_element_and_intersection_sparsity_pattern(V_h),
)
a_h += LocalElementIntegralBilinearForm(LocalLaplaceIntegrand(diffusion))
a_h += (
LocalCouplingIntersectionIntegralBilinearForm(
LocalLaplaceIPDGInnerCouplingIntegrand(symmetry_factor, diffusion, weight)
+ LocalIPDGInnerPenaltyIntegrand(penalty_parameter, weight)
),
{},
ApplyOnInnerIntersectionsOnce(grid),
)
a_h += (
LocalIntersectionIntegralBilinearForm(
LocalIPDGBoundaryPenaltyIntegrand(penalty_parameter, weight)
+ LocalLaplaceIPDGDirichletCouplingIntegrand(symmetry_factor, diffusion)
),
{},
ApplyOnCustomBoundaryIntersections(grid, boundary_info, DirichletBoundary()),
)

walker = Walker(grid)
walker.append(a_h)
walker.append(l_h)
walker.walk()

u_h = DiscreteFunction(V_h, name="u_h")
a_h.apply_inverse(l_h.vector, u_h.dofs.vector)

return u_h

from dune.gdt import visualize_function

u_h = discretize_elliptic_ipdg_dirichlet_zero(
grid, kappa, f,
symmetry_factor=1, penalty_parameter=16, weight=1) # SIPDG scheme

_ = visualize_function(u_h)