SpECTRE  v2022.01.03
DG-Subcell

Functions and classes specific to the discontinuous Galerkin method supplemented with a finite volume or finite difference subcell limiter. Can also be thought of as a DG-FD hybrid method. More...

## Namespaces

namespace  evolution::dg::subcell
Implementation of a generic finite volume/conservative finite difference subcell limiter.

namespace  evolution::dg::subcell::fv
Code specific to a finite volume subcell limiter.

namespace  evolution::dg::subcell::fd
Code specific to a conservative finite difference subcell limiter.

namespace  evolution::dg::subcell::fd::Actions
Actions specific to using a finite-difference subcell method.

namespace  evolution::dg::subcell::fd::Tags
Tags for the DG-subcell finite difference solver

namespace  evolution::dg::subcell::Tags
Tags for the DG-subcell solver

namespace  evolution::dg::subcell::OptionTags
Option tags for the DG-subcell solver.

namespace  evolution::dg::subcell::Actions
Actions for the DG-subcell hybrid solver.

## Enumerations

enum class  evolution::dg::subcell::ActiveGrid { Dg , Subcell }
The grid that is currently being used for the DG-subcell evolution.

## Functions

template<size_t Dim>
const Matrixevolution::dg::subcell::fd::projection_matrix (const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
Computes the projection matrix in Dim dimensions going from a DG mesh to a conservative finite difference subcell mesh.

template<size_t Dim>
const Matrixevolution::dg::subcell::fd::reconstruction_matrix (const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
Computes the matrix needed for reconstructing the DG solution from the subcell solution. More...

template<size_t Dim>
DataVector evolution::dg::subcell::fd::project (const DataVector &dg_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
Project the variable dg_u onto the subcell grid with extents subcell_extents. More...

template<size_t Dim>
void evolution::dg::subcell::fd::project (gsl::not_null< DataVector * > subcell_u, const DataVector &dg_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
Project the variable dg_u onto the subcell grid with extents subcell_extents. More...

template<typename SubcellTagList , typename DgTagList , size_t Dim>
void evolution::dg::subcell::fd::project (const gsl::not_null< Variables< SubcellTagList > * > subcell_u, const Variables< DgTagList > &dg_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
Project the variable dg_u onto the subcell grid with extents subcell_extents. More...

template<typename TagList , size_t Dim>
Variables< TagList > evolution::dg::subcell::fd::project (const Variables< TagList > &dg_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
Project the variable dg_u onto the subcell grid with extents subcell_extents. More...

template<size_t Dim>
DataVector evolution::dg::subcell::fd::reconstruct (const DataVector &subcell_u_times_projected_det_jac, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
reconstruct the variable subcell_u_times_projected_det_jac onto the DG grid dg_mesh. More...

template<size_t Dim>
void evolution::dg::subcell::fd::reconstruct (gsl::not_null< DataVector * > dg_u, const DataVector &subcell_u_times_projected_det_jac, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
reconstruct the variable subcell_u_times_projected_det_jac onto the DG grid dg_mesh. More...

template<typename SubcellTagList , typename DgTagList , size_t Dim>
void evolution::dg::subcell::fd::reconstruct (const gsl::not_null< Variables< DgTagList > * > dg_u, const Variables< SubcellTagList > &subcell_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
reconstruct the variable subcell_u_times_projected_det_jac onto the DG grid dg_mesh. More...

template<typename TagList , size_t Dim>
Variables< TagList > evolution::dg::subcell::fd::reconstruct (const Variables< TagList > &subcell_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents)
reconstruct the variable subcell_u_times_projected_det_jac onto the DG grid dg_mesh. More...

## Detailed Description

Functions and classes specific to the discontinuous Galerkin method supplemented with a finite volume or finite difference subcell limiter. Can also be thought of as a DG-FD hybrid method.

## ◆ project() [1/4]

template<size_t Dim>
 DataVector evolution::dg::subcell::fd::project ( const DataVector & dg_u, const Mesh< Dim > & dg_mesh, const Index< Dim > & subcell_extents )

Project the variable dg_u onto the subcell grid with extents subcell_extents.

Note
In the return-by-gsl::not_null with Variables interface, the SubcellTagList and the DgtagList must be the same when all tag prefixes are removed. Typically the Tags::Inactive prefix will be used.

## ◆ project() [2/4]

template<typename SubcellTagList , typename DgTagList , size_t Dim>
 void evolution::dg::subcell::fd::project ( const gsl::not_null< Variables< SubcellTagList > * > subcell_u, const Variables< DgTagList > & dg_u, const Mesh< Dim > & dg_mesh, const Index< Dim > & subcell_extents )

Project the variable dg_u onto the subcell grid with extents subcell_extents.

Note
In the return-by-gsl::not_null with Variables interface, the SubcellTagList and the DgtagList must be the same when all tag prefixes are removed. Typically the Tags::Inactive prefix will be used.

## ◆ project() [3/4]

template<typename TagList , size_t Dim>
 Variables< TagList > evolution::dg::subcell::fd::project ( const Variables< TagList > & dg_u, const Mesh< Dim > & dg_mesh, const Index< Dim > & subcell_extents )

Project the variable dg_u onto the subcell grid with extents subcell_extents.

Note
In the return-by-gsl::not_null with Variables interface, the SubcellTagList and the DgtagList must be the same when all tag prefixes are removed. Typically the Tags::Inactive prefix will be used.

## ◆ project() [4/4]

template<size_t Dim>
 void evolution::dg::subcell::fd::project ( gsl::not_null< DataVector * > subcell_u, const DataVector & dg_u, const Mesh< Dim > & dg_mesh, const Index< Dim > & subcell_extents )

Project the variable dg_u onto the subcell grid with extents subcell_extents.

Note
In the return-by-gsl::not_null with Variables interface, the SubcellTagList and the DgtagList must be the same when all tag prefixes are removed. Typically the Tags::Inactive prefix will be used.

## ◆ reconstruct() [1/4]

template<size_t Dim>
 DataVector evolution::dg::subcell::fd::reconstruct ( const DataVector & subcell_u_times_projected_det_jac, const Mesh< Dim > & dg_mesh, const Index< Dim > & subcell_extents )

reconstruct the variable subcell_u_times_projected_det_jac onto the DG grid dg_mesh.

In general we wish that the reconstruction operator is the pseudo-inverse of the projection operator. On curved meshes this means we either need to compute a (time-dependent) reconstruction and projection matrix on each DG element, or we expand the determinant of the Jacobian on the basis, accepting the aliasing errors from that. We accept the aliasing errors in favor of the significantly reduced computational overhead. This means that the projection and reconstruction operators are only inverses of each other if both operate on $$u J$$ where $$u$$ is the variable being projected and $$J$$ is the determinant of the Jacobian. That is, the matrices are guaranteed to satisfy $$\mathcal{R}(\mathcal{P}(u J))=u J$$. If the mesh is regular Cartesian, then this isn't an issue. Furthermore, if we reconstruct $$uJ/\mathcal{P}(J)$$ we again recover the exact DG solution. Doing the latter has the advantage that, in general, we are ideally projecting to the subcells much more often than reconstructing from them (a statement that we would rather use DG more than the subcells).

## ◆ reconstruct() [2/4]

template<typename SubcellTagList , typename DgTagList , size_t Dim>
 void evolution::dg::subcell::fd::reconstruct ( const gsl::not_null< Variables< DgTagList > * > dg_u, const Variables< SubcellTagList > & subcell_u, const Mesh< Dim > & dg_mesh, const Index< Dim > & subcell_extents )

reconstruct the variable subcell_u_times_projected_det_jac onto the DG grid dg_mesh.

In general we wish that the reconstruction operator is the pseudo-inverse of the projection operator. On curved meshes this means we either need to compute a (time-dependent) reconstruction and projection matrix on each DG element, or we expand the determinant of the Jacobian on the basis, accepting the aliasing errors from that. We accept the aliasing errors in favor of the significantly reduced computational overhead. This means that the projection and reconstruction operators are only inverses of each other if both operate on $$u J$$ where $$u$$ is the variable being projected and $$J$$ is the determinant of the Jacobian. That is, the matrices are guaranteed to satisfy $$\mathcal{R}(\mathcal{P}(u J))=u J$$. If the mesh is regular Cartesian, then this isn't an issue. Furthermore, if we reconstruct $$uJ/\mathcal{P}(J)$$ we again recover the exact DG solution. Doing the latter has the advantage that, in general, we are ideally projecting to the subcells much more often than reconstructing from them (a statement that we would rather use DG more than the subcells).

## ◆ reconstruct() [3/4]

template<typename TagList , size_t Dim>
 Variables< TagList > evolution::dg::subcell::fd::reconstruct ( const Variables< TagList > & subcell_u, const Mesh< Dim > & dg_mesh, const Index< Dim > & subcell_extents )

reconstruct the variable subcell_u_times_projected_det_jac onto the DG grid dg_mesh.

In general we wish that the reconstruction operator is the pseudo-inverse of the projection operator. On curved meshes this means we either need to compute a (time-dependent) reconstruction and projection matrix on each DG element, or we expand the determinant of the Jacobian on the basis, accepting the aliasing errors from that. We accept the aliasing errors in favor of the significantly reduced computational overhead. This means that the projection and reconstruction operators are only inverses of each other if both operate on $$u J$$ where $$u$$ is the variable being projected and $$J$$ is the determinant of the Jacobian. That is, the matrices are guaranteed to satisfy $$\mathcal{R}(\mathcal{P}(u J))=u J$$. If the mesh is regular Cartesian, then this isn't an issue. Furthermore, if we reconstruct $$uJ/\mathcal{P}(J)$$ we again recover the exact DG solution. Doing the latter has the advantage that, in general, we are ideally projecting to the subcells much more often than reconstructing from them (a statement that we would rather use DG more than the subcells).

## ◆ reconstruct() [4/4]

template<size_t Dim>
 void evolution::dg::subcell::fd::reconstruct ( gsl::not_null< DataVector * > dg_u, const DataVector & subcell_u_times_projected_det_jac, const Mesh< Dim > & dg_mesh, const Index< Dim > & subcell_extents )

reconstruct the variable subcell_u_times_projected_det_jac onto the DG grid dg_mesh.

In general we wish that the reconstruction operator is the pseudo-inverse of the projection operator. On curved meshes this means we either need to compute a (time-dependent) reconstruction and projection matrix on each DG element, or we expand the determinant of the Jacobian on the basis, accepting the aliasing errors from that. We accept the aliasing errors in favor of the significantly reduced computational overhead. This means that the projection and reconstruction operators are only inverses of each other if both operate on $$u J$$ where $$u$$ is the variable being projected and $$J$$ is the determinant of the Jacobian. That is, the matrices are guaranteed to satisfy $$\mathcal{R}(\mathcal{P}(u J))=u J$$. If the mesh is regular Cartesian, then this isn't an issue. Furthermore, if we reconstruct $$uJ/\mathcal{P}(J)$$ we again recover the exact DG solution. Doing the latter has the advantage that, in general, we are ideally projecting to the subcells much more often than reconstructing from them (a statement that we would rather use DG more than the subcells).

## ◆ reconstruction_matrix()

template<size_t Dim>
 const Matrix & evolution::dg::subcell::fd::reconstruction_matrix ( const Mesh< Dim > & dg_mesh, const Index< Dim > & subcell_extents )

Computes the matrix needed for reconstructing the DG solution from the subcell solution.

Reconstructing the DG solution from the FD solution is a bit more involved than projecting the DG solution to the FD subcells. Denoting the projection operator by $$\mathcal{P}$$ and the reconstruction operator by $$\mathcal{R}$$, we desire the property

\begin{align*} \mathcal{R}(\mathcal{P}(u_{\breve{\imath}} J_{\breve{\imath}}))=u_{\breve{\imath}} J_{\breve{\imath}}, \end{align*}

where $$\breve{\imath}$$ denotes a grid point on the DG grid, $$u$$ is the solution on the DG grid, and $$J$$ is the determinant of the Jacobian on the DG grid. We also require that the integral of the conserved variables over the subcells is equal to the integral over the DG element. That is,

\begin{align*} \int_{\Omega}u \,d^3x =\int_{\Omega} \underline{u} \,d^3x \Longrightarrow \int_{\Omega}u J \,d^3\xi=\int_{\Omega} \underline{u} J \,d^3\xi, \end{align*}

where $$\underline{u}$$ is the solution on the subcells. Because the number of subcell points is larger than the number of DG points, we need to solve a constrained linear least squares problem to reconstruct the DG solution from the subcells.

The final reconstruction matrix is given by

\begin{align*} R_{\breve{\jmath}\underline{i}} &=\left\{(2 \mathcal{P}\otimes\mathcal{P})^{-1}2\mathcal{P} - (2 \mathcal{P}\otimes\mathcal{P})^{-1}\vec{w}\left[\mathbf{w}(2 \mathcal{P}\otimes\mathcal{P})^{-1}\vec{w}\right]^{-1}\mathbf{w}(2 \mathcal{P}\otimes\mathcal{P})^{-1}2\mathcal{P} + (2 \mathcal{P}\otimes\mathcal{P})^{-1}\vec{w}\left[\mathbf{w}(2 \mathcal{P}\otimes\mathcal{P})^{-1}\vec{w}\right]^{-1}\vec{\underline{w}} \right\}_{\breve{\jmath}\underline{i}}, \end{align*}

where $$\vec{w}$$ is the vector of integration weights on the DG element, $$\mathbf{w}=w_{\breve{l}}\delta_{\breve{l}\breve{\jmath}}$$, and $$\vec{\underline{w}}$$ is the vector of integration weights over the subcells. The integration weights $$\vec{\underline{w}}$$ on the subcells are those for 6th-order integration on a uniform mesh.