Functionality related to discontinuous Galerkin discretizations of elliptic equations. More...
Namespaces | |
| Actions | |
| Actions related to elliptic discontinuous Galerkin schemes. | |
| OptionTags | |
| Option tags related to elliptic discontinuous Galerkin schemes. | |
| subdomain_operator | |
| Items related to the restriction of the DG operator to an element-centered subdomain. | |
| 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 | NormalizeFaceNormal |
Normalize face normals, possibly using a background metric. Set InvMetricTag to void to normalize face normals with the Euclidean magnitude. More... | |
Typedefs | |
| template<typename PrimalFields , typename PrimalFluxes > | |
| using | BoundaryData = ::dg::SimpleBoundaryData< tmpl::append< db::wrap_tags_in<::Tags::NormalDotFlux, PrimalFields >, db::wrap_tags_in<::Tags::NormalDotFlux, PrimalFluxes >, db::wrap_tags_in< Tags::NormalDotFluxForJump, PrimalFields >, tmpl::list< Tags::ElementSize > >, tmpl::list< Tags::PerpendicularNumPoints > > |
| Data that is projected to mortars and communicated across element boundaries. | |
Functions | |
| template<typename PrimalMortarFields , typename PrimalMortarFluxes , size_t Dim> | |
| BoundaryData< PrimalMortarFields, PrimalMortarFluxes > | zero_boundary_data_on_mortar (const Direction< Dim > &direction, const Mesh< Dim > &mesh, const Scalar< DataVector > &face_normal_magnitude, const Mesh< Dim - 1 > &mortar_mesh, const ::dg::MortarSize< Dim - 1 > &mortar_size) noexcept |
Construct elliptic::dg::BoundaryData assuming the variable data on the element is zero, and project it to the mortar. | |
| template<typename System , bool Linearized, typename... Args> | |
| void | prepare_mortar_data (Args &&... args) noexcept |
| Prepare data on mortars so they can be communicated to neighbors. More... | |
| template<typename System , bool Linearized, typename... Args> | |
| void | apply_operator (Args &&... args) noexcept |
| Apply the elliptic DG operator. More... | |
| template<typename System , typename... Args> | |
| void | impose_inhomogeneous_boundary_conditions_on_source (Args &&... args) noexcept |
| 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) noexcept |
| The penalty factor in internal penalty fluxes. More... | |
Functionality related to discontinuous Galerkin discretizations of elliptic equations.
The following is a brief overview of the elliptic DG schemes that are implemented here. A publication that describes the schemes in detail is in preparation and will be referenced here.
The DG schemes apply to any elliptic PDE that can be formulated in first-order flux-form
\begin{equation} -\partial_i F^i_\alpha + S_\alpha = f_\alpha(x) \end{equation}
in terms of fluxes \(F^i_\alpha\), sources \(S_\alpha\) and fixed-sources \(f_\alpha(x)\). The fluxes and sources are functionals of the "primal" system variables \(u_A(x)\) and their corresponding "auxiliary" variables \(v_A(x)\). The uppercase index enumerates the variables such that \(v_A\) is the auxiliary variable corresponding to \(u_A\). Greek indices enumerate all variables. In documentation related to particular elliptic systems we generally use the canonical system-specific symbols for the fields in place of these indices. See the Poisson::FirstOrderSystem and the Elasticity::FirstOrderSystem for examples.
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-Lobatto collocation points as basis functions. Support for Legendre-Gauss collocation points can be added if needed. 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.
\begin{align} (n_i F^i_v)^* &= \frac{1}{2} n_i \left( F^i_v(u^\mathrm{int}) + F^i_v(u^\mathrm{ext})\right) \\ (n_i F^i_u)^* &= \frac{1}{2} n_i \left( F^i_u(\partial_j F^j_v(u^\mathrm{int}) - S_v(u^\mathrm{int})) + F^i_u(\partial_j F^j_v(u^\mathrm{ext}) - S_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) \end{align}
Note that \(n_i\) denotes the face normal on the "interior" side of the element under consideration. We assume \(n^\mathrm{ext}_i=-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 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. [55], 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. [55], section 7.2, or [2]):
\begin{align} (n_i F^i_{v_j})^* &= n_j u^* = \frac{1}{2} n_j \left( u^\mathrm{int} + u^\mathrm{ext}\right) \\ (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) \end{align}
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{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 [106]). 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 [106]. Note that we use the number of points \(N_\mathrm{points}\) where [106] 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.
|
noexcept |
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.
|
noexcept |
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).
|
noexcept |
The penalty factor in internal penalty fluxes.
The penalty factor is computed as
\begin{equation} \sigma = C \frac{N_\text{points}^2}{h} \end{equation}
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 [55].
|
noexcept |
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.