|
SpECTRE
v2023.01.13
|
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... | |
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) |
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) |
| 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... | |
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, 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-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. [72], 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. [72], section 7.2, or [4]):
\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 [146]). 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 [146]. Note that we use the number of points \(N_\mathrm{points}\) where [146] 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.
| 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.
| 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).
| 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{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 [72].
| 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.