evolution::dg::subcell Namespace Reference

Implementation of a generic finite volume/conservative finite difference subcell limiter. More...

Namespaces
namespace	Actions
	Actions for the DG-subcell hybrid solver.
namespace	fd
	Code specific to a conservative finite difference subcell limiter.
namespace	fv
	Code specific to a finite volume subcell limiter.
namespace	OptionTags
	Option tags for the DG-subcell solver.
namespace	Tags
	Tags for the DG-subcell solver

Classes
struct	BackgroundGrVars
	Allocate or assign background general relativity quantities on cell-centered and face-centered FD grid points, for evolution systems run on a curved spacetime without solving Einstein equations (e.g. ValenciaDivclean, ForceFree),. More...
class	GhostData
struct	InitialTciData
	Used to communicate the RDMP and TCI status/decision during initialization. More...
struct	RdmpTciData
	Holds data needed for the relaxed discrete maximum principle troubled-cell indicator. More...
struct	SetInterpolators
	Sets the intrp::IrregularInterpolants for interpolating to ghost zone data at block boundaries. More...
class	SubcellOptions
	Holds the system-agnostic subcell parameters, such as numbers controlling when to switch between DG and subcell. More...

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

Functions
std::ostream &	operator<< (std::ostream &os, ActiveGrid active_grid)
template<typename... CorrectionTags, typename BoundaryCorrection , typename... PackageFieldTags, typename... BoundaryTermsVolumeTags>
void	compute_boundary_terms (const gsl::not_null< Variables< tmpl::list< CorrectionTags... > > * > boundary_corrections_on_face, const BoundaryCorrection &boundary_correction, const Variables< tmpl::list< PackageFieldTags... > > &upper_packaged_data, const Variables< tmpl::list< PackageFieldTags... > > &lower_packaged_data, const db::Access &box, tmpl::list< BoundaryTermsVolumeTags... >)
template<bool OverwriteInternalMortarData, size_t Dim, typename DgPackageFieldTags >
void	correct_package_data (const gsl::not_null< Variables< DgPackageFieldTags > * > lower_packaged_data, const gsl::not_null< Variables< DgPackageFieldTags > * > upper_packaged_data, const size_t logical_dimension_to_operate_in, const Element< Dim > &element, const Mesh< Dim > &subcell_volume_mesh, const std::unordered_map< DirectionalId< Dim >, evolution::dg::MortarData< Dim >, boost::hash< DirectionalId< Dim > > > &mortar_data, const size_t variables_to_offset_in_dg_grid)
	Project the DG package data to the subcells. Data received from a neighboring element doing DG is always projected, while the data we sent to our neighbors before doing a rollback from DG to subcell is only projected if OverwriteInternalMortarData is true. More...
template<typename DgTag , typename SubcellTag , typename DbTagsList >
const DgTag::type &	get_active_tag (const db::DataBox< DbTagsList > &box)
	Retrieve a tag from the active grid.
template<typename DgTag , typename SubcellTag , typename DbTagsList >
const DgTag::type &	get_inactive_tag (const db::DataBox< DbTagsList > &box)
	Retrieve a tag from the inactive grid.
template<typename DbTagsList >
int	get_tci_decision (const db::DataBox< DbTagsList > &box)
	Returns the value of evolution::dg::subcell::Tags::TciDecision.
bool	operator!= (const GhostData &lhs, const GhostData &rhs)
std::ostream &	operator<< (std::ostream &os, const GhostData &ghost_data)
template<bool InsertIntoMap, size_t Dim>
void	insert_or_update_neighbor_volume_data (gsl::not_null< DirectionalIdMap< Dim, GhostData > * > ghost_data_ptr, const DataVector &neighbor_subcell_data, const size_t number_of_rdmp_vars_in_buffer, const DirectionalId< Dim > &directional_element_id, const Mesh< Dim > &neighbor_mesh, const Element< Dim > &element, const Mesh< Dim > &subcell_mesh, size_t number_of_ghost_zones, const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &Neighbor_dg_to_fd_interpolants)
	Check whether neighbor_subcell_data is FD or DG, and either insert or copy into ghost_data_ptr the FD data (projecting if neighbor_subcell_data is DG data). More...
template<size_t Dim>
void	insert_neighbor_rdmp_and_volume_data (gsl::not_null< RdmpTciData * > rdmp_tci_data_ptr, gsl::not_null< DirectionalIdMap< Dim, GhostData > * > ghost_data_ptr, const DataVector &received_neighbor_subcell_data, size_t number_of_rdmp_vars, const DirectionalId< Dim > &directional_element_id, const Mesh< Dim > &neighbor_mesh, const Element< Dim > &element, const Mesh< Dim > &subcell_mesh, size_t number_of_ghost_zones, const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &neighbor_dg_to_fd_interpolants)
	Check whether the neighbor sent is DG volume or FD ghost data, and orient project DG volume data if necessary. More...
template<size_t VolumeDim, typename DgComputeSubcellNeighborPackagedData >
void	neighbor_reconstructed_face_solution (gsl::not_null< db::Access * > box, gsl::not_null< std::pair< const TimeStepId, DirectionalIdMap< VolumeDim, std::tuple< Mesh< VolumeDim >, Mesh< VolumeDim - 1 >, std::optional< DataVector >, std::optional< DataVector >, ::TimeStepId, int > > > * > received_temporal_id_and_data)
	Invoked in directions where the neighbor is doing subcell, this function computes the neighbor data on the mortar via reconstruction on nearest neighbor subcells. More...
template<size_t Dim>
void	neighbor_tci_decision (gsl::not_null< db::Access * > box, const std::pair< const TimeStepId, DirectionalIdMap< Dim, std::tuple< Mesh< Dim >, Mesh< Dim - 1 >, std::optional< DataVector >, std::optional< DataVector >, ::TimeStepId, int > > > &received_temporal_id_and_data)
	Copies the neighbors' TCI decisions into subcell::Tags::NeighborTciDecisions<Dim>
template<size_t Dim, typename SymmList , typename IndexList >
bool	persson_tci (const Tensor< DataVector, SymmList, IndexList > &tensor, const Mesh< Dim > &dg_mesh, const double alpha, const size_t num_highest_modes)
	Troubled cell indicator using the idea of spectral falloff by [150]. More...
template<typename Metavariables , typename DbTagsList , size_t Dim>
void	prepare_neighbor_data (const gsl::not_null< DirectionMap< Dim, DataVector > * > all_neighbor_data_for_reconstruction, const gsl::not_null< Mesh< Dim > * > ghost_data_mesh, const gsl::not_null< db::DataBox< DbTagsList > * > box, const Variables< db::wrap_tags_in< ::Tags::Flux, typename Metavariables::system::flux_variables, tmpl::size_t< Dim >, Frame::Inertial > > &volume_fluxes)
	Add local data for our and our neighbor's relaxed discrete maximum principle troubled-cell indicator, and for reconstruction. More...
template<typename... EvolvedVarsTags>
int	rdmp_tci (const Variables< tmpl::list< EvolvedVarsTags... > > &active_grid_candidate_evolved_vars, const Variables< tmpl::list< Tags::Inactive< EvolvedVarsTags >... > > &inactive_grid_candidate_evolved_vars, const DataVector &max_of_past_variables, const DataVector &min_of_past_variables, const double rdmp_delta0, const double rdmp_epsilon)
	Troubled cell indicator using a relaxed discrete maximum principle, comparing the candidate solution with the past solution in the element and its neighbors. More...
template<typename... EvolvedVarsTags>
std::pair< DataVector, DataVector >	rdmp_max_min (const Variables< tmpl::list< EvolvedVarsTags... > > &active_grid_evolved_vars, const Variables< tmpl::list< Tags::Inactive< EvolvedVarsTags >... > > &inactive_grid_evolved_vars, const bool include_inactive_grid)
	get the max and min of each component of the active and inactive variables. If include_inactive_grid is false then only the max over the active_grid_evolved_vars for each component is returned.
int	rdmp_tci (const DataVector &max_of_current_variables, const DataVector &min_of_current_variables, const DataVector &max_of_past_variables, const DataVector &min_of_past_variables, double rdmp_delta0, double rdmp_epsilon)
	Check if the current variables satisfy the RDMP. Returns an integer 0 if cell is not troubled and an integer i+1 if the [i]-th element of the input vector is responsible for failing the RDMP. More...
void	pup (PUP::er &p, RdmpTciData &rdmp_tci_data)
void	operator| (PUP::er &p, RdmpTciData &rdmp_tci_data)
bool	operator== (const RdmpTciData &lhs, const RdmpTciData &rhs)
bool	operator!= (const RdmpTciData &lhs, const RdmpTciData &rhs)
std::ostream &	operator<< (std::ostream &os, const RdmpTciData &t)
template<size_t Dim, typename DbTagsList >
void	store_reconstruction_order_in_databox (const gsl::not_null< db::DataBox< DbTagsList > * > box, const std::optional< std::array< gsl::span< std::uint8_t >, Dim > > &reconstruction_order)
	Copies the reconstruction_order into the DataBox tag evolution::dg::subcell::Tags::ReconstructionOrder if reconstruction_order has a value.
template<size_t Dim, typename VectorType >
void	slice_tensor_for_subcell (const gsl::not_null< Tensor< VectorType, Symmetry<>, index_list<> > * > sliced_scalar, const Tensor< VectorType, Symmetry<>, index_list<> > &volume_scalar, const Index< Dim > &subcell_extents, size_t number_of_ghost_points, const Direction< Dim > &direction, const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &fd_to_neighbor_fd_interpolants)
bool	operator!= (const SubcellOptions &lhs, const SubcellOptions &rhs)
template<typename... DgEvolvedVarsTags, typename... SubcellEvolvedVarsTags>
int	two_mesh_rdmp_tci (const Variables< tmpl::list< DgEvolvedVarsTags... > > &dg_evolved_vars, const Variables< tmpl::list< SubcellEvolvedVarsTags... > > &subcell_evolved_vars, const double rdmp_delta0, const double rdmp_epsilon)
	Troubled cell indicator using a relaxed discrete maximum principle, comparing the solution on two grids at the same point in time. More...
void	add_cartesian_flux_divergence (gsl::not_null< DataVector * > dt_var, double one_over_delta, const DataVector &inv_jacobian, const DataVector &boundary_correction, const Index< 1 > &subcell_extents, size_t dimension)
	Compute and add the 2nd-order flux divergence on a Cartesian mesh to the cell-centered time derivatives.
void	add_cartesian_flux_divergence (gsl::not_null< DataVector * > dt_var, double one_over_delta, const DataVector &inv_jacobian, const DataVector &boundary_correction, const Index< 2 > &subcell_extents, size_t dimension)
	Compute and add the 2nd-order flux divergence on a Cartesian mesh to the cell-centered time derivatives.
void	add_cartesian_flux_divergence (gsl::not_null< DataVector * > dt_var, double one_over_delta, const DataVector &inv_jacobian, const DataVector &boundary_correction, const Index< 3 > &subcell_extents, size_t dimension)
	Compute and add the 2nd-order flux divergence on a Cartesian mesh to the cell-centered time derivatives.
template<size_t Dim>
DirectionMap< Dim, DataVector >	slice_data (const DataVector &volume_subcell_vars, const Index< Dim > &subcell_extents, const size_t number_of_ghost_points, const std::unordered_set< Direction< Dim > > &directions_to_slice, const size_t additional_buffer, const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &fd_to_neighbor_fd_interpolants)
	Slice the subcell variables needed for neighbors to perform reconstruction. More...
template<size_t Dim, typename TagList >
DirectionMap< Dim, DataVector >	slice_data (const Variables< TagList > &volume_subcell_vars, const Index< Dim > &subcell_extents, const size_t number_of_ghost_points, const std::unordered_set< Direction< Dim > > &directions_to_slice, const size_t additional_buffer, const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &fd_to_neighbor_fd_interpolants)
	Slice the subcell variables needed for neighbors to perform reconstruction. More...
template<size_t Dim, typename VectorType , typename... Structure>
void	slice_tensor_for_subcell (const gsl::not_null< Tensor< VectorType, Structure... > * > sliced_tensor, const Tensor< VectorType, Structure... > &volume_tensor, const Index< Dim > &subcell_extents, size_t number_of_ghost_points, const Direction< Dim > &direction, const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &fd_to_neighbor_fd_interpolants)
	Slice a single volume tensor for a given direction and slicing depth (number of ghost points). More...
template<size_t Dim, typename VectorType , typename... Structure>
Tensor< VectorType, Structure... >	slice_tensor_for_subcell (const Tensor< VectorType, Structure... > &volume_tensor, const Index< Dim > &subcell_extents, size_t number_of_ghost_points, const Direction< Dim > &direction, const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &fd_to_neighbor_fd_interpolants)
	Slice a single volume tensor for a given direction and slicing depth (number of ghost points). More...
template<size_t Dim, typename TagList >
void	slice_variable (const gsl::not_null< Variables< TagList > * > &sliced_subcell_vars, const Variables< TagList > &volume_subcell_vars, const Index< Dim > &subcell_extents, const size_t ghost_zone_size, const Direction< Dim > &direction, const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &fd_to_neighbor_fd_interpolants)
	Slice a Variables object on subcell mesh for a given direction and slicing depth (number of ghost points).
template<size_t Dim, typename TagList >
Variables< TagList >	slice_variable (const Variables< TagList > &volume_subcell_vars, const Index< Dim > &subcell_extents, const size_t ghost_zone_size, const Direction< Dim > &direction, const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &fd_to_neighbor_fd_interpolants)
	Slice a Variables object on subcell mesh for a given direction and slicing depth (number of ghost points).

Detailed Description

Implementation of a generic finite volume/conservative finite difference subcell limiter.

Our implementation of a finite volume (FV) or finite difference (FD) subcell limiter (SCL) follows [59]. Other implementations of a subcell limiter exist, e.g. [169] [33] [90]. Our implementation and that of [59] are a generalization of the Multidimensional Optimal Order Detection (MOOD) algorithm [39] [55] [56] [120].

Function Documentation

correct_package_data()

template<bool OverwriteInternalMortarData, size_t Dim, typename DgPackageFieldTags >

void evolution::dg::subcell::correct_package_data	(	const gsl::not_null< Variables< DgPackageFieldTags > * >	lower_packaged_data,
		const gsl::not_null< Variables< DgPackageFieldTags > * >	upper_packaged_data,
		const size_t	logical_dimension_to_operate_in,
		const Element< Dim > &	element,
		const Mesh< Dim > &	subcell_volume_mesh,
		const std::unordered_map< DirectionalId< Dim >, evolution::dg::MortarData< Dim >, boost::hash< DirectionalId< Dim > > > &	mortar_data,
		const size_t	variables_to_offset_in_dg_grid
	)

Project the DG package data to the subcells. Data received from a neighboring element doing DG is always projected, while the data we sent to our neighbors before doing a rollback from DG to subcell is only projected if OverwriteInternalMortarData is true.

In order for the hybrid DG-FD/FV scheme to be conservative between elements using DG and elements using subcell, the boundary terms must be the same on both elements. In practice this means the boundary corrections \(G+D\) must be computed on the same grid. Consider the element doing subcell which receives data from an element doing DG. In this case the DG element's ingredients going into \(G+D\) are projected to the subcells and then \(G+D\) are computed on the subcells. Similarly, for strict conservation the element doing DG must first project the data it sent to the neighbor to the subcells, then compute \(G+D\) on the subcells, and finally reconstrct \(G+D\) back to the DG grid before lifting \(G+D\) to the volume.

This function updates the packaged_data (ingredients into \(G+D\)) received by an element doing subcell by projecting the neighbor's DG data onto the subcells. Note that this is only half of what is required for strict conservation, the DG element must also compute \(G+D\) on the subcells. Note that we currently do not perform the other half of the correction needed to be strictly conservative.

If we are retaking a time step after the DG step failed then maintaining conservation requires additional care. If OverwriteInternalMortarData is true then the local (the element switching from DG to subcell) ingredients into \(G+D\) are projected and overwrite the data computed from the FD reconstruction to the interface. However, even this is insufficient to guarantee conservation. To guarantee conservation (which we do not currently do) the correction \(G+D\) must be computed on the DG grid and then projected to the subcells.

Note that our practical experience shows that since the DG-subcell hybrid scheme switches to the subcell solver before the local solution contains discontinuities, strict conservation is not necessary between DG and FD/FV regions. This was also observed with a block-adaptive finite difference AMR code [37]

The variable variables_to_offset_in_dg_grid is used in cases like the combined generalized harmonic and GRMHD system where the DG grid uses boundary corrections for both GH and GRMHD, but the subcell grid only does GRMHD.

insert_neighbor_rdmp_and_volume_data()

template<size_t Dim>

void evolution::dg::subcell::insert_neighbor_rdmp_and_volume_data	(	gsl::not_null< RdmpTciData * >	rdmp_tci_data_ptr,
		gsl::not_null< DirectionalIdMap< Dim, GhostData > * >	ghost_data_ptr,
		const DataVector &	received_neighbor_subcell_data,
		size_t	number_of_rdmp_vars,
		const DirectionalId< Dim > &	directional_element_id,
		const Mesh< Dim > &	neighbor_mesh,
		const Element< Dim > &	element,
		const Mesh< Dim > &	subcell_mesh,
		size_t	number_of_ghost_zones,
		const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &	neighbor_dg_to_fd_interpolants
	)

Check whether the neighbor sent is DG volume or FD ghost data, and orient project DG volume data if necessary.

This is intended to be used by the ReceiveDataForReconstruction action.

insert_or_update_neighbor_volume_data()

template<bool InsertIntoMap, size_t Dim>

void evolution::dg::subcell::insert_or_update_neighbor_volume_data	(	gsl::not_null< DirectionalIdMap< Dim, GhostData > * >	ghost_data_ptr,
		const DataVector &	neighbor_subcell_data,
		const size_t	number_of_rdmp_vars_in_buffer,
		const DirectionalId< Dim > &	directional_element_id,
		const Mesh< Dim > &	neighbor_mesh,
		const Element< Dim > &	element,
		const Mesh< Dim > &	subcell_mesh,
		size_t	number_of_ghost_zones,
		const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &	Neighbor_dg_to_fd_interpolants
	)

Check whether neighbor_subcell_data is FD or DG, and either insert or copy into ghost_data_ptr the FD data (projecting if neighbor_subcell_data is DG data).

This is intended to be used during a rollback from DG to make sure neighbor data is projected to the FD grid.

neighbor_reconstructed_face_solution()

template<size_t VolumeDim, typename DgComputeSubcellNeighborPackagedData >

void evolution::dg::subcell::neighbor_reconstructed_face_solution	(	gsl::not_null< db::Access * >	box,
		gsl::not_null< std::pair< const TimeStepId, DirectionalIdMap< VolumeDim, std::tuple< Mesh< VolumeDim >, Mesh< VolumeDim - 1 >, std::optional< DataVector >, std::optional< DataVector >, ::TimeStepId, int > > > * >	received_temporal_id_and_data
	)

Invoked in directions where the neighbor is doing subcell, this function computes the neighbor data on the mortar via reconstruction on nearest neighbor subcells.

The data needed for reconstruction is copied over into subcell::Tags::GhostDataForReconstruction. Additionally, the max/min of the evolved variables from neighboring elements that is used for the relaxed discrete maximum principle troubled-cell indicator is combined with the data from the local element and stored in subcell::Tags::DataForRdmpTci. We handle the RDMP data now because it is sent in the same buffer as the data for reconstruction.

A list of all the directions that are doing subcell is created and then passed to the mutator Metavariables::SubcellOptions::DgComputeSubcellNeighborPackagedData::apply, which must return a

DirectionalIdMap<volume_dim, DataVector>

which holds the reconstructed dg_packaged_data on the face (stored in the DataVector) for the boundary correction. A std::vector<DirectionalId<volume_dim>> holding the list of mortars that need to be reconstructed to is passed in as the last argument to Metavariables::SubcellOptions::DgComputeSubcellNeighborPackagedData::apply.

persson_tci()

template<size_t Dim, typename SymmList , typename IndexList >

bool evolution::dg::subcell::persson_tci	(	const Tensor< DataVector, SymmList, IndexList > &	tensor,
		const Mesh< Dim > &	dg_mesh,
		const double	alpha,
		const size_t	num_highest_modes
	)

Troubled cell indicator using the idea of spectral falloff by [150].

Consider a discontinuity sensing quantity \(U\), which is typically a scalar but could be a tensor of any rank. Let \(U\) have the 1d spectral decomposition (generalization to higher-dimensional tensor product bases is done dimension-by-dimension):

\begin{align*} U(x)=\sum_{i=0}^{N}c_i P_i(x), \end{align*}

where \(P_i(x)\) are the basis functions, in our case the Legendre polynomials, and \(c_i\) are the spectral coefficients. We then define a filtered solution \(\hat{U}\) as

\begin{align*} \hat{U}(x)=\sum_{i=N+1-M}^{N} c_i P_i(x). \end{align*}

where \(M\) is the number of highest modes to include in the filtered solution.

Note that when an exponential filter is being used to deal with aliasing, lower modes can be included in \(\hat{U}\). The main goal of \(\hat{U}\) is to measure how much power is in the highest modes, which are the modes responsible for Gibbs phenomena.

A cell is troubled if

\begin{align*} \frac{(\hat{U}, \hat{U})}{(U, U)} > (N + 1 - M)^{-\alpha} \end{align*}

where \((\cdot,\cdot)\) is an inner product, which we take to be the Euclidean \(L_2\) norm (i.e. we do not divide by the number of grid points since that cancels out anyway) computed as

\begin{align*} (U, U) \equiv \sqrt{\sum_{i=0}^N U_i^2} \end{align*}

where \(U_i\) are nodal values of the quantity \(U\) at grid points.

Typically, \(\alpha=4.0\) and \(M=1\) is a good choice.

prepare_neighbor_data()

template<typename Metavariables , typename DbTagsList , size_t Dim>

void evolution::dg::subcell::prepare_neighbor_data	(	const gsl::not_null< DirectionMap< Dim, DataVector > * >	all_neighbor_data_for_reconstruction,
		const gsl::not_null< Mesh< Dim > * >	ghost_data_mesh,
		const gsl::not_null< db::DataBox< DbTagsList > * >	box,
		const Variables< db::wrap_tags_in< ::Tags::Flux, typename Metavariables::system::flux_variables, tmpl::size_t< Dim >, Frame::Inertial > > &	volume_fluxes
	)

Add local data for our and our neighbor's relaxed discrete maximum principle troubled-cell indicator, and for reconstruction.

The local maximum and minimum of the evolved variables is added to Tags::DataForRdmpTci for the RDMP TCI. Then the data needed by neighbor elements to do reconstruction on the FD grid is sent. The data to be sent is computed in the mutator Metavariables::SubcellOptions::GhostVariables, which returns a Variables of the tensors to send to the neighbors. The main reason for having the mutator GhostVariables is to allow sending primitive or characteristic variables for reconstruction.

Note

If all neighbors are using DG then we send our DG volume data without orienting it. This elides the expense of projection and slicing. If any neighbors are doing FD, we project and slice to all neighbors. A future optimization would be to measure the cost of slicing data, and figure out how many neighbors need to be doing FD before it's worth projecting and slicing vs. just projecting to the ghost cells. Another optimization is to always send DG volume data to neighbors that we think are doing DG, so that we can elide any slicing or projection cost in that direction.

rdmp_tci() [1/2]

int evolution::dg::subcell::rdmp_tci	(	const DataVector &	max_of_current_variables,
		const DataVector &	min_of_current_variables,
		const DataVector &	max_of_past_variables,
		const DataVector &	min_of_past_variables,
		double	rdmp_delta0,
		double	rdmp_epsilon
	)

Check if the current variables satisfy the RDMP. Returns an integer 0 if cell is not troubled and an integer i+1 if the [i]-th element of the input vector is responsible for failing the RDMP.

Let the candidate solution be denoted by \(u^\star_{\alpha}(t^{n+1})\). Then the RDMP requires that

\begin{align*} \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right) - \delta_\alpha \le u^\star_{\alpha}(t^{n+1}) \le \max_{\forall\mathcal{N}} \left(u_{\alpha}(t^n)\right) + \delta_\alpha \end{align*}

where \(\mathcal{N}\) are either the Neumann or Voronoi neighbors and the element itself, and \(\delta_\alpha\) is a parameter defined below that relaxes the discrete maximum principle (DMP). When computing \(\max(u_\alpha)\) and \(\min(u_\alpha)\) over a DG element that is not using subcells we first project the DG solution to the subcells and then compute the maximum and minimum over both the DG grid and the subcell grid. However, when a DG element is using subcells we compute the maximum and minimum of \(u_\alpha(t^n)\) over the subcells only. Note that the maximum and minimum values of \(u^\star_\alpha\) are always computed over both the DG and the subcell grids, even when using the RDMP to check if the reconstructed DG solution would be admissible.

The parameter \(\delta_\alpha\) is given by:

\begin{align*} \delta_\alpha = \max\left(\delta_{0},\epsilon \left(\max_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right) - \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)\right) \right), \end{align*}

where we typically take \(\delta_{0}=10^{-4}\) and \(\epsilon=10^{-3}\).

If all checks are passed and cell is not troubled, returns an integer 0. Otherwise returns an 1-based index of the element in the input DataVector that fails the check.

e.g. Suppose we have three variables to check RDMP so that max_of_current_variables.size() == 3. If RDMP TCI flags max_of_current_variables[1], min_of_current_variables[1], .. (and so on) as troubled, returned integer value is 2.

Once cell is marked as troubled, checks for the remaining part of the input std::vectors are skipped. In the example above, for instance if [1]-th component of inputs is flagged as troubled, checking the remaining index [2] is skipped.

rdmp_tci() [2/2]

template<typename... EvolvedVarsTags>

int evolution::dg::subcell::rdmp_tci	(	const Variables< tmpl::list< EvolvedVarsTags... > > &	active_grid_candidate_evolved_vars,
		const Variables< tmpl::list< Tags::Inactive< EvolvedVarsTags >... > > &	inactive_grid_candidate_evolved_vars,
		const DataVector &	max_of_past_variables,
		const DataVector &	min_of_past_variables,
		const double	rdmp_delta0,
		const double	rdmp_epsilon
	)

Troubled cell indicator using a relaxed discrete maximum principle, comparing the candidate solution with the past solution in the element and its neighbors.

Let the candidate solution be denoted by \(u^\star_{\alpha}(t^{n+1})\). Then the RDMP requires that

\begin{align*} \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right) - \delta_\alpha \le u^\star_{\alpha}(t^{n+1}) \le \max_{\forall\mathcal{N}} \left(u_{\alpha}(t^n)\right) + \delta_\alpha \end{align*}

where \(\mathcal{N}\) are either the Neumann or Voronoi neighbors and the element itself, and \(\delta_\alpha\) is a parameter defined below that relaxes the discrete maximum principle (DMP). When computing \(\max(u_\alpha)\) and \(\min(u_\alpha)\) over a DG element that is not using subcells we first project the DG solution to the subcells and then compute the maximum and minimum over both the DG grid and the subcell grid. However, when a DG element is using subcells we compute the maximum and minimum of \(u_\alpha(t^n)\) over the subcells only. Note that the maximum and minimum values of \(u^\star_\alpha\) are always computed over both the DG and the subcell grids, even when using the RDMP to check if the reconstructed DG solution would be admissible.

The parameter \(\delta_\alpha\) is given by:

\begin{align*} \delta_\alpha = \max\left(\delta_{0},\epsilon \left(\max_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right) - \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)\right) \right), \end{align*}

where we typically take \(\delta_{0}=10^{-4}\) and \(\epsilon=10^{-3}\).

If all checks are passed and cell is not troubled, returns an integer 0. Otherwise returns 1-based index of the tag in the input Variables that fails the check. For instance, if we have following two Variables objects as candidate solutions on active and inactive grids

• Variables<tmpl::list<DgVar1, DgVar2, DgVar3>>
• Variables<tmpl::list<SubVar1, SubVar2, SubVar3>>

and TCI flags the second pair DgVar2 and SubVar2 not satisfying two-mesh RDMP criteria, returned value is 2 since the second pair of tags failed the check.

Note

Once a single pair of tags fails to satisfy the check, checks for the remaining part of the input variables are skipped. In the example above, for instance if the second pair (DgVar2,SubVar2) is flagged as troubled, the third pair (DgVar3,SubVar3) is ignored and not checked.

slice_data() [1/2]

template<size_t Dim>

DirectionMap< Dim, DataVector > evolution::dg::subcell::slice_data	(	const DataVector &	volume_subcell_vars,
		const Index< Dim > &	subcell_extents,
		const size_t	number_of_ghost_points,
		const std::unordered_set< Direction< Dim > > &	directions_to_slice,
		const size_t	additional_buffer,
		const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &	fd_to_neighbor_fd_interpolants
	)

Slice the subcell variables needed for neighbors to perform reconstruction.

Note that we slice to a grid that is against the boundary of the element but is several ghost points deep. This is in contrast to the slicing used in the DG method which is to the boundary of the element only.

The number_of_ghost_points will depend on the number of neighboring points the reconstruction method needs that is used on the subcell. The directions_to_slice determines in which directions data is sliced. Generally this will be the directions in which the element has neighbors.

The data always has the same ordering as the volume data (tags have the same ordering, grid points are x-varies-fastest).

The additional_buffer argument is used to add extra padding to the result storage to be used for example for sending the RDMP TCI data. This eliminates expensive data copying.

slice_data() [2/2]

template<size_t Dim, typename TagList >

DirectionMap< Dim, DataVector > evolution::dg::subcell::slice_data	(	const Variables< TagList > &	volume_subcell_vars,
		const Index< Dim > &	subcell_extents,
		const size_t	number_of_ghost_points,
		const std::unordered_set< Direction< Dim > > &	directions_to_slice,
		const size_t	additional_buffer,
		const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &	fd_to_neighbor_fd_interpolants
	)

Slice the subcell variables needed for neighbors to perform reconstruction.

Note that we slice to a grid that is against the boundary of the element but is several ghost points deep. This is in contrast to the slicing used in the DG method which is to the boundary of the element only.

The number_of_ghost_points will depend on the number of neighboring points the reconstruction method needs that is used on the subcell. The directions_to_slice determines in which directions data is sliced. Generally this will be the directions in which the element has neighbors.

The data always has the same ordering as the volume data (tags have the same ordering, grid points are x-varies-fastest).

The additional_buffer argument is used to add extra padding to the result storage to be used for example for sending the RDMP TCI data. This eliminates expensive data copying.

slice_tensor_for_subcell() [1/2]

template<size_t Dim, typename VectorType , typename... Structure>

void evolution::dg::subcell::slice_tensor_for_subcell	(	const gsl::not_null< Tensor< VectorType, Structure... > * >	sliced_tensor,
		const Tensor< VectorType, Structure... > &	volume_tensor,
		const Index< Dim > &	subcell_extents,
		size_t	number_of_ghost_points,
		const Direction< Dim > &	direction,
		const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &	fd_to_neighbor_fd_interpolants
	)

Slice a single volume tensor for a given direction and slicing depth (number of ghost points).

Note that the second to last argument has the type Direction<Dim>, not a DirectionMap (cf. evolution::dg::subcell::slice_data)

slice_tensor_for_subcell() [2/2]

template<size_t Dim, typename VectorType , typename... Structure>

Tensor< VectorType, Structure... > evolution::dg::subcell::slice_tensor_for_subcell	(	const Tensor< VectorType, Structure... > &	volume_tensor,
		const Index< Dim > &	subcell_extents,
		size_t	number_of_ghost_points,
		const Direction< Dim > &	direction,
		const DirectionalIdMap< Dim, std::optional< intrp::Irregular< Dim > > > &	fd_to_neighbor_fd_interpolants
	)

Slice a single volume tensor for a given direction and slicing depth (number of ghost points).

Note that the second to last argument has the type Direction<Dim>, not a DirectionMap (cf. evolution::dg::subcell::slice_data)

two_mesh_rdmp_tci()

template<typename... DgEvolvedVarsTags, typename... SubcellEvolvedVarsTags>

int evolution::dg::subcell::two_mesh_rdmp_tci	(	const Variables< tmpl::list< DgEvolvedVarsTags... > > &	dg_evolved_vars,
		const Variables< tmpl::list< SubcellEvolvedVarsTags... > > &	subcell_evolved_vars,
		const double	rdmp_delta0,
		const double	rdmp_epsilon
	)

Troubled cell indicator using a relaxed discrete maximum principle, comparing the solution on two grids at the same point in time.

Checks that the subcell solution \(\underline{u}\) and the DG solution \(u\) satisfy

\begin{align*} \min(u)-\delta \le \underline{u} \le \max(u)+\delta \end{align*}

where

\begin{align*} \delta = \max\left[\delta_0, \epsilon(\max(u) - \min(u))\right] \end{align*}

where \(\delta_0\) and \(\epsilon\) are constants controlling the maximum absolute and relative change allowed when projecting the DG solution to the subcell grid. We currently specify one value of \(\delta_0\) and \(\epsilon\) for all variables, but this could be generalized to choosing the allowed variation in a variable-specific manner.

If all checks are passed and cell is not troubled, returns an integer 0. Otherwise returns 1-based index of the tag in the input Variables that fails the check. For instance, if we have

• Variables<tmpl::list<DgVar1, DgVar2, DgVar3>> for dg_evolved_vars
• Variables<tmpl::list<SubVar1, SubVar2, SubVar3>> for subcell_evolved_vars

as inputs and TCI flags the second pair DgVar2 and SubVar2 not satisfying two-mesh RDMP criteria, returned value is 2 since the second pair of tags failed the check.

Note

Once a single pair of tags fails to satisfy the check, checks for the remaining part of the input variables are skipped. In the example above, for instance if the second pair (DgVar2,SubVar2) is flagged, the third pair (DgVar3,SubVar3) is ignored and not checked.