elliptic::dg Namespace Reference

Functionality related to discontinuous Galerkin discretizations of elliptic equations. More...

Namespaces
namespace	Actions
	Actions related to elliptic discontinuous Galerkin schemes.
namespace	OptionTags
	Option tags related to elliptic discontinuous Galerkin schemes.
namespace	subdomain_operator
	Items related to the restriction of the DG operator to an element-centered subdomain.
namespace	Tags
	DataBox tags related to elliptic discontinuous Galerkin schemes.

Classes
struct	InitializeBackground
	Initialize background quantities for the elliptic DG operator, possibly including the metric necessary for normalizing face normals. More...
struct	InitializeFacesAndMortars
	Initialize the geometry on faces and mortars for the elliptic DG operator. More...
struct	InitializeGeometry
	Initialize the background-independent geometry for the elliptic DG operator. More...
struct	MortarData
	Boundary data on both sides of a mortar. More...
struct	ProjectGeometry

Typedefs
template<typename PrimalFields , typename PrimalFluxes >
using	BoundaryData = ::dg::SimpleBoundaryData< tmpl::append< PrimalFields, db::wrap_tags_in<::Tags::NormalDotFlux, PrimalFields > >, tmpl::list<> >
	Data that is projected to mortars and communicated across element boundaries.

Functions
template<typename System , bool Linearized, typename... Args>
void	prepare_mortar_data (Args &&... args)
	Prepare data on mortars so they can be communicated to neighbors. More...
template<typename System , bool Linearized, typename... Args>
void	apply_operator (Args &&... args)
	Apply the elliptic DG operator. More...
template<typename System , typename... Args>
void	impose_inhomogeneous_boundary_conditions_on_source (Args &&... args)
	For linear systems, impose inhomogeneous boundary conditions as contributions to the fixed sources (i.e. the RHS of the equations). More...
DataVector	penalty (const DataVector &element_size, size_t num_points, double penalty_parameter)
	The penalty factor in internal penalty fluxes. More...

Detailed Description

Functionality related to discontinuous Galerkin discretizations of elliptic equations.

The following is a brief overview of the elliptic DG schemes that are implemented here. The scheme is described in detail in [69].

The DG schemes apply to any elliptic PDE that can be formulated in first-order flux-form, as detailed by elliptic::protocols::FirstOrderSystem. The DG discretization of equations in this first-order form amounts to projecting the equations on the set of basis functions that we also use to represent the fields on the computational grid. The currently implemented DG operator uses Lagrange interpolating polynomials w.r.t. Legendre-Gauss or Legendre-Gauss-Lobatto collocation points as basis functions. Skipping all further details here, the discretization results in a linear equation $A(u)=b$ over all grid points and primal variables. Solving the elliptic equations amounts to numerically inverting the DG operator $A$, typically without ever constructing the full matrix but by employing an iterative linear solver that repeatedly applies the DG operator to "test data". Note that the DG operator applies directly to the primal variables. Auxiliary variables are only computed temporarily and don't inflate the size of the operator. This means the DG operator essentially computes second derivatives of the primal variables, modified by the fluxes and sources of the system as well as by DG boundary corrections that couple grid points across element boundaries.

Boundary corrections:

In this implementation we employ the "internal penalty" DG scheme that couples grid points across nearest-neighbor elements through the fluxes:

$$\begin{array}{}\text{(1)}& u^{\ast }& =\frac{1}{2}(u^{\mathrm{i}\mathrm{n}\mathrm{t}}+u^{\mathrm{e}\mathrm{x}\mathrm{t}})\text{(2)}& (n_i{F}^i{)}^{\ast }& =\frac{1}{2}{n}_i(F_{\mathrm{i}\mathrm{n}\mathrm{t}}^i+F_{\mathrm{e}\mathrm{x}\mathrm{t}}^i)-\sigma n_i{F}^i(n_j(u^{\mathrm{i}\mathrm{n}\mathrm{t}}-u^{\mathrm{e}\mathrm{x}\mathrm{t}}))\end{array}$$

Note that $n_i$ denotes the face normal on the "interior" side of the element under consideration. We assume $n_i^{\mathrm{e}\mathrm{x}\mathrm{t}}=-n_i$ in the implementation, i.e. face normals don't depend on the dynamic variables (which may be discontinuous on element faces). This is the case for the problems we are expecting to solve, because those will be on fixed background metrics (e.g. a conformal metric for the XCTS system). Numerically, the face normals on either side of a mortar may nonetheless be different because the two faces adjacent to the mortar may resolve them at different resolutions.

Also note that the numerical fluxes intentionally don't depend on the auxiliary field values $v$. This property allows us to communicate data for both the primal and auxiliary boundary corrections together, instead of communicating them in two steps. If we were to resort to a two-step communication we could replace the derivatives in $(n_i{F}^i{)}^{\ast }$ with $v$, which would result in a generalized "stabilized central flux" that is slightly less sparse than the internal penalty flux (see e.g. [95], section 7.2). We could also choose to ignore the fluxes in the penalty term, but preliminary tests suggest that this may hurt convergence.

For a Poisson system (see Poisson::FirstOrderSystem) this numerical flux reduces to the standard internal penalty flux (see e.g. [95], section 7.2, or [5]):

$$\begin{array}{}\text{(3)}& u^{\ast }& =\frac{1}{2}(u^{\mathrm{i}\mathrm{n}\mathrm{t}}+u^{\mathrm{e}\mathrm{x}\mathrm{t}})\text{(4)}& (n_i{F}^i{)}^{\ast }& =n_i{v}_i^{\ast }=\frac{1}{2}{n}_i({\partial }_i{u}^{\mathrm{i}\mathrm{n}\mathrm{t}}+{\partial }_i{u}^{\mathrm{e}\mathrm{x}\mathrm{t}})-\sigma (u^{\mathrm{i}\mathrm{n}\mathrm{t}}-u^{\mathrm{e}\mathrm{x}\mathrm{t}})\end{array}$$

where a sum over repeated indices is assumed, since the equation is formulated on a Euclidean geometry.

The penalty factor $\sigma$ is responsible for removing zero eigenmodes and impacts the conditioning of the linear operator to be solved. See elliptic::dg::penalty for details. For the element size that goes into computing the penalty we choose $h=\frac{J_{\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}}}{J_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}}$, i.e. the ratio of Jacobi determinants from logical to inertial coordinates in the element volume and on the element face, both evaluated on the face (see [200]). Since both $N_{\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}}$ and $h$ can be different on either side of the element boundary we take the maximum of $N_{\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}}$ and the pointwise minimum of $h$ across the element boundary as is done in [200]. Note that we use the number of points $N_{\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}}$ where [200] uses the polynomial degree $N_{\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}}-1$ because we found unstable configurations on curved meshes when using the polynomial degree. Optimizing the penalty on curved meshes is subject to further investigation.

Discontinuous fluxes:

The DG operator also supports systems with potentially discontinuous fluxes, such as elasticity with layered materials. The way to handle the discontinuous fluxes in the DG scheme is described in [202]. Essentially, we evaluate the penalty term in Eq. $\text{2}$ on both sides of an element boundary and take the average. The other terms in the numerical flux remain unchanged.

Function Documentation

apply_operator()

template<typename System , bool Linearized, typename... Args>

void elliptic::dg::apply_operator

(

Args &&...

args

)

Apply the elliptic DG operator.

This function applies the elliptic DG operator on an element, assuming all data on mortars is already available. Use the elliptic::dg::prepare_mortar_data function to prepare mortar data on neighboring elements, then communicate the data and insert them on the "remote" side of the mortars before calling this function.

impose_inhomogeneous_boundary_conditions_on_source()

template<typename System , typename... Args>

void elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source

(

Args &&...

args

)

For linear systems, impose inhomogeneous boundary conditions as contributions to the fixed sources (i.e. the RHS of the equations).

This function exists because the DG operator must typically be linear, but even for linear elliptic equations we typically apply boundary conditions with a constant, and therefore nonlinear, contribution. Standard examples are inhomogeneous (i.e. non-zero) Dirichlet or Neumann boundary conditions. This nonlinear contribution can be added to the fixed sources, leaving only the linearized boundary conditions in the DG operator. For standard constant Dirichlet or Neumann boundary conditions the linearization is of course just zero.

This function essentially feeds zero variables through the nonlinear operator and subtracts the result from the fixed sources: b -= A(x=0).

penalty()

DataVector elliptic::dg::penalty	(	const DataVector &	element_size,
		size_t	num_points,
		double	penalty_parameter
	)

The penalty factor in internal penalty fluxes.

The penalty factor is computed as

$$\begin{array}{}\text{(5)}& \sigma =C\frac{N_{\text{points}}^2}{h}\end{array}$$

where $N_{\text{points}}=1+N_p$ is the number of points (or one plus the polynomial degree) and $h$ is a measure of the element size. Both quantities are taken perpendicular to the face of the DG element that the penalty is being computed on. $C$ is the "penalty parameter". For details see section 7.2 in [95].

prepare_mortar_data()

template<typename System , bool Linearized, typename... Args>

void elliptic::dg::prepare_mortar_data

(

Args &&...

args

)

Prepare data on mortars so they can be communicated to neighbors.

Call this function on all elements and communicate the mortar data, then call elliptic::dg::apply_operator.