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
const Matrix &	evolution::dg::subcell::fd::projection_matrix (const Mesh< 1 > &dg_mesh, size_t subcell_extents, const Spectral::Quadrature &subcell_quadrature)
	Computes the projection matrix in 1 dimension going from a DG mesh to a conservative finite difference subcell mesh.
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. More...
const Matrix &	evolution::dg::subcell::fd::projection_matrix (const Mesh< 1 > &dg_mesh, size_t subcell_extents, size_t ghost_zone_size, Side side)
	Computes the projection matrix in 1 dimension going from a DG mesh to a conservative finite difference subcell mesh for only the ghost zones. 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::project_to_faces (const DataVector &dg_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents, const size_t &face_direction)
	Project the variable dg_u onto the subcell grid with extents subcell_extents. More...
template<size_t Dim>
void	evolution::dg::subcell::fd::project_to_faces (gsl::not_null< DataVector * > subcell_u, const DataVector &dg_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents, const size_t &face_direction)
	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_to_faces (const gsl::not_null< Variables< SubcellTagList > * > subcell_u, const Variables< DgTagList > &dg_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents, const size_t &face_direction)
	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_to_faces (const Variables< TagList > &dg_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents, const size_t &face_direction)
	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, ReconstructionMethod reconstruction_method)
	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, ReconstructionMethod reconstruction_method)
	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, const ReconstructionMethod reconstruction_method)
	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, const ReconstructionMethod reconstruction_method)
	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.

Function Documentation

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.

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.

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.

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.

project_to_faces() [1/4]

template<size_t Dim>

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

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.

project_to_faces() [2/4]

template<typename SubcellTagList , typename DgTagList , size_t Dim>

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

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.

project_to_faces() [3/4]

template<typename TagList , size_t Dim>

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

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.

project_to_faces() [4/4]

template<size_t Dim>

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

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.

projection_matrix()

const Matrix & evolution::dg::subcell::fd::projection_matrix	(	const Mesh< 1 > &	dg_mesh,
		size_t	subcell_extents,
		size_t	ghost_zone_size,
		Side	side
	)

Computes the projection matrix in 1 dimension going from a DG mesh to a conservative finite difference subcell mesh for only the ghost zones.

This is used when a neighbor sends DG volume data and we need to switch to FD. In this case we need to project the DG volume data onto the ghost zone cells.

Note

Currently assumes a max ghost zone size of 5 and a minimum ghost zone size of 2.

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,
		ReconstructionMethod	reconstruction_method
	)

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,
		const ReconstructionMethod	reconstruction_method
	)

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,
		const ReconstructionMethod	reconstruction_method
	)

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,
		ReconstructionMethod	reconstruction_method
	)

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.