Run this tutorial
Click here to run this tutorial on mybinder.org:
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 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)
Download the code:
example__ipdg_stationary_heat_equation.md
example__ipdg_stationary_heat_equation.ipynb