elliptic::dg::NumericalFluxes::FirstOrderInternalPenalty< Dim, FluxesComputerTag, FieldTagsList, AuxiliaryFieldTagsList > Struct Template Reference

The internal penalty flux for first-order elliptic equations. More...

#include <FirstOrderInternalPenalty.hpp>

## Detailed Description

### template<size_t Dim, typename FluxesComputerTag, typename FieldTagsList, typename AuxiliaryFieldTagsList> struct elliptic::dg::NumericalFluxes::FirstOrderInternalPenalty< Dim, FluxesComputerTag, FieldTagsList, AuxiliaryFieldTagsList >

The internal penalty flux for first-order elliptic equations.

### Details

Computes the internal penalty numerical flux (see e.g. [34], section 7.2) dotted with the interface unit normal.

We implement here a suggested generalization of the internal penalty flux for any set of elliptic PDEs up to second order in derivatives. It is designed for fluxes (i.e. principal parts of the PDEs) that may depend on the dynamic fields, but do so at most linearly. This is the case for the velocity potential equation of binary neutron stars, for example. The linearity is only necessary to impose inhomogeneous boundary conditions as contributions to the fixed source (see elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource).

We reduce the second-order elliptic PDEs to first-order form by introducing an auxiliary variable $$v$$ for each primal variable $$u$$ (i.e. for each variable whose second derivative appears in the equations). Then, the equations take the flux-form

\begin{align} -\partial_i F^i_A + S_A = f_A \end{align}

where the fluxes $$F^i_A(u,v)$$ and sources $$S_A(u,v)$$ are indexed (with capital letters) by the variables and we have defined $$f_A(x)$$ as those sources that are independent of the variables. Note that the fluxes are functions of the primal and auxiliary variables; e.g. for a Poisson system $$\{u,v_i\}$$ on a spatial metric $$\gamma_{ij}$$ they are simply $$F^i_u(v)=\sqrt{\gamma}\gamma^{ij} v_j$$ and $$F^i_{v_j}(u)=u\delta^i_j$$ (see Poisson::FirstOrderSystem) or for an Elasticity system $$\{\xi^i,S_{ij}\}$$ with Young's tensor $$Y^{ijkl}$$ they are $$F^i_{\xi^j}(S)=Y^{ijkl}S_{kl}$$ and $$F^i_{S_{jk}}(\xi)=\delta^i_{(j}\xi_{k)}$$. Since each primal flux is commonly only a function of the corresponding auxiliary variable and vice-versa, we omit the other function arguments here to lighten the notational load.

We now choose the internal penalty numerical fluxes $$(n_i F^i_A)^*$$ as follows for each primal variable $$u$$ and their corresponding auxiliary variable $$v$$:

\begin{align} (n_i F^i_u)^* &= \frac{1}{2} n_i \left( F^i_u(\partial_j F^j_v(u^\mathrm{int})) + F^i_u(\partial_j F^j_v(u^\mathrm{ext})) \right) - \sigma n_i \left(F^i_u(n_j F^j_v(u^\mathrm{int})) - F^i_u( n_j F^j_v(u^\mathrm{ext}))\right) \\ (n_i F^i_v)^* &= \frac{1}{2} n_i \left(F^i_v(u^\mathrm{int}) + F^i_v(u^\mathrm{ext})\right) \end{align}

Note that we have assumed $$n^\mathrm{ext}_i=-n_i$$ here, 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).

Also note that the numerical fluxes intentionally don't depend on the auxiliary field values $$v$$. This property allows us to use the numerical fluxes also for the second-order (or primal) DG formulation, where we remove the need for an auxiliary variable. For the first-order system we could replace the divergence in $$(n_i F^i_u)^*$$ with $$v$$, which would result in a generalized "stabilized central flux" that is slightly less sparse than the internal penalty flux (see e.g. [34], 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 above) this numeric flux reduces to the standard internal penalty flux (see e.g. [34], section 7.2, or [2])

\begin{align} (n_i F^i_u)^* &= n_i v_i^* = \frac{1}{2} n_i \left(\partial_i u^\mathrm{int} + \partial_i u^\mathrm{ext}\right) - \sigma \left(u^\mathrm{int} - u^\mathrm{ext}\right) \\ (n_i F^i_{v_j})^* &= n_j u^* = \frac{1}{2} n_j \left(u^\mathrm{int} + u^\mathrm{ext}\right) \end{align}

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

This generalization of the internal penalty flux is based on unpublished work by Nils L. Fischer (nils..nosp@m.fisc.nosp@m.her@a.nosp@m.ei.m.nosp@m.pg.de).

The penalty factor $$\sigma$$ is responsible for removing zero eigenmodes and impacts the conditioning of the linear operator to be solved. It can be chosen as $$\sigma=C\frac{N_\mathrm{points}^2}{h}$$ where $$N_\mathrm{points}$$ is the number of collocation points (i.e. the polynomial degree plus 1), $$h$$ is a measure of the element size in inertial coordinates (both measured perpendicular to the interface under consideration) and $$C\geq 1$$ is a free parameter (see e.g. [34], section 7.2). For the element size we choose $$h=\frac{J_\mathrm{volume}}{J_\mathrm{face}}$$, 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 [69]). Since both $$N_\mathrm{points}$$ and $$h$$ can be different on either side of the element boundary we take the maximum of $$N_\mathrm{points}$$ and the pointwise minimum of $$h$$ across the element boundary as is done in [69]. Note that we use the number of points $$N_\mathrm{points}$$ where [69] uses the polynomial degree $$N_\mathrm{points} - 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.

The documentation for this struct was generated from the following file: