Run this tutorial

Click here to run this tutorial on try on
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
%load_ext wurlitzer
%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

import inspect
from discretize_elliptic_ipdg import 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)

        The grid, given as a GridProvider from dune.xt.grid.
        Diffusion function with values in R^{d x d}, anything that dune.xt.functions.GridFunction
        can handle.
        Right hand side source function with values in R, anything that
        dune.xt.functions.GridFunction can handle.
        Usually one of -1, 0, 1, determines the IPDG scheme (see above).
        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`.
        Usually 1 or diffusion, determines the IPDG scheme (see above).

        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(
    a_h += LocalElementIntegralBilinearForm(LocalLaplaceIntegrand(diffusion))
    a_h += (
            LocalLaplaceIPDGInnerCouplingIntegrand(symmetry_factor, diffusion, weight)
            + LocalIPDGInnerPenaltyIntegrand(penalty_parameter, weight)
    a_h += (
            LocalIPDGBoundaryPenaltyIntegrand(penalty_parameter, weight)
            + LocalLaplaceIPDGDirichletCouplingIntegrand(symmetry_factor, diffusion)
        ApplyOnCustomBoundaryIntersections(grid, boundary_info, DirichletBoundary()),

    walker = Walker(grid)

    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)

Download the code: example__ipdg_stationary_heat_equation.ipynb