SpECTRE
v2022.05.05

Implementation of a generic finite volume/conservative finite difference subcell limiter. More...
Namespaces  
namespace  Actions 
Actions for the DGsubcell 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 DGsubcell solver.  
namespace  Tags 
Tags for the DGsubcell solver  
Classes  
struct  RdmpTciData 
Holds data needed for the relaxed discrete maximum principle troubledcell indicator. More...  
class  SubcellOptions 
Holds the systemagnostic 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 DGsubcell evolution.  
Functions  
std::ostream &  operator<< (std::ostream &os, ActiveGrid active_grid) 
template<typename... CorrectionTags, typename BoundaryCorrection , typename... PackageFieldTags>  
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) 
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< std::pair< Direction< Dim >, ElementId< Dim > >, evolution::dg::MortarData< Dim >, boost::hash< std::pair< Direction< Dim >, ElementId< Dim > > > > &mortar_data) 
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 Metavariables , typename DbTagsList >  
void  neighbor_reconstructed_face_solution (const gsl::not_null< db::DataBox< DbTagsList > * > box, const gsl::not_null< std::pair< const TimeStepId, FixedHashMap< maximum_number_of_neighbors(Metavariables::volume_dim), std::pair< Direction< Metavariables::volume_dim >, ElementId< Metavariables::volume_dim > >, std::tuple< Mesh< Metavariables::volume_dim  1 >, std::optional< std::vector< double > >, std::optional< std::vector< double > >, ::TimeStepId >, boost::hash< std::pair< Direction< Metavariables::volume_dim >, ElementId< Metavariables::volume_dim > > > > > * > 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, typename SymmList , typename IndexList >  
bool  persson_tci (const Tensor< DataVector, SymmList, IndexList > &tensor, const Mesh< Dim > &dg_mesh, const double alpha) 
Troubled cell indicator using spectral falloff of [99]. More...  
template<typename Metavariables , typename DbTagsList >  
DirectionMap< Metavariables::volume_dim, std::vector< double > >  prepare_neighbor_data (const gsl::not_null< db::DataBox< DbTagsList > * > box) 
Add local data for our and our neighbor's relaxed discrete maximum principle troubledcell indicator, and compute and slice data needed for the neighbor cells. More...  
template<typename... EvolvedVarsTags>  
bool  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 std::vector< double > &max_of_past_variables, const std::vector< double > &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< std::vector< double >, std::vector< double > >  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.  
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 TagList >  
DirectionMap< Dim, std::vector< double > >  slice_data (const Variables< TagList > &volume_subcell_vars, const Index< Dim > &subcell_extents, const size_t number_of_ghost_points, const DirectionMap< Dim, bool > &directions_to_slice) 
Slice the subcell variables needed for neighbors to perform reconstruction. More...  
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) 
bool  operator== (const SubcellOptions &lhs, const SubcellOptions &rhs) 
bool  operator!= (const SubcellOptions &lhs, const SubcellOptions &rhs) 
template<typename... EvolvedVarsTags>  
bool  two_mesh_rdmp_tci (const Variables< tmpl::list< EvolvedVarsTags... > > &dg_evolved_vars, const Variables< tmpl::list< Tags::Inactive< EvolvedVarsTags >... > > &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 2ndorder flux divergence on a Cartesian mesh to the cellcentered 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 2ndorder flux divergence on a Cartesian mesh to the cellcentered 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 2ndorder flux divergence on a Cartesian mesh to the cellcentered time derivatives.  
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) 
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) 
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) 
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) 
Slice a Variables object on subcell mesh for a given direction and slicing depth (number of ghost points).  
template<size_t Dim>  
void  tci_status (gsl::not_null< Scalar< DataVector > * > status, const Mesh< Dim > &dg_mesh, const Mesh< Dim > &subcell_mesh, subcell::ActiveGrid active_grid, const std::deque< subcell::ActiveGrid > &tci_history) 
Set status to 1 if the cell is marked as needing subcell, otherwise set status to 0 . status has grid points on the currently active mesh. More...  
template<size_t Dim>  
Scalar< DataVector >  tci_status (const Mesh< Dim > &dg_mesh, const Mesh< Dim > &subcell_mesh, subcell::ActiveGrid active_grid, const std::deque< subcell::ActiveGrid > &tci_history) 
Set status to 1 if the cell is marked as needing subcell, otherwise set status to 0 . status has grid points on the currently active mesh. More...  
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 [44]. Other implementations of a subcell limiter exist, e.g. [113] [23] [66]. Our implementation and that of [44] are a generalization of the Multidimensional Optimal Order Detection (MOOD) algorithm [29] [41] [42] [81].
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< std::pair< Direction< Dim >, ElementId< Dim > >, evolution::dg::MortarData< Dim >, boost::hash< std::pair< Direction< Dim >, ElementId< Dim > > > > &  mortar_data  
) 
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 DGFD/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 DGsubcell 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 blockadaptive finite difference AMR code [27]
void evolution::dg::subcell::neighbor_reconstructed_face_solution  (  const gsl::not_null< db::DataBox< DbTagsList > * >  box, 
const gsl::not_null< std::pair< const TimeStepId, FixedHashMap< maximum_number_of_neighbors(Metavariables::volume_dim), std::pair< Direction< Metavariables::volume_dim >, ElementId< Metavariables::volume_dim > >, std::tuple< Mesh< Metavariables::volume_dim  1 >, std::optional< std::vector< double > >, std::optional< std::vector< double > >, ::TimeStepId >, boost::hash< std::pair< Direction< Metavariables::volume_dim >, ElementId< Metavariables::volume_dim > > > > > * >  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::NeighborDataForReconstruction
. Additionally, the max/min of the evolved variables from neighboring elements that is used for the relaxed discrete maximum principle troubledcell 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::DgOuterCorrectionPackageData::apply
, which must return a
which holds the reconstructed dg_packaged_data
on the face (stored in the std::vector<double>
) for the boundary correction. A std::vector<std::pair<Direction<volume_dim>, ElementId<volume_dim>>>
holding the list of mortars that need to be reconstructed to is passed in as the last argument to Metavariables::SubcellOptions::DgOuterCorrectionPackageData::apply
.
bool evolution::dg::subcell::persson_tci  (  const Tensor< DataVector, SymmList, IndexList > &  tensor, 
const Mesh< Dim > &  dg_mesh,  
const double  alpha  
) 
Troubled cell indicator using spectral falloff of [99].
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 higherdimensional tensor product bases is done dimensionbydimension):
\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)=c_N P_N(x). \end{align*}
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. We define the discontinuity indicator \(s^\Omega\) as
\begin{align*} s^\Omega=\log_{10}\left(\frac{(\hat{U}, \hat{U})}{(U, U)}\right), \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). A cell is troubled if \(s^\Omega > \alpha \log_{10}(N)\). Typically, \(\alpha=4\) is a good choice.
The parameter zero_cutoff
is used to avoid division and logarithms of small numbers, which can be wildly fluctuating because of roundoff errors. We do not check the TCI for tensor components when \(L_2(\hat{U}) \leq \epsilon L_2(U)\), where \(\epsilon\) is the zero_cutoff
. If all components are skipped the TCI returns false
, i.e. the cell is not troubled.
DirectionMap< Metavariables::volume_dim, std::vector< double > > evolution::dg::subcell::prepare_neighbor_data  (  const gsl::not_null< db::DataBox< DbTagsList > * >  box  ) 
Add local data for our and our neighbor's relaxed discrete maximum principle troubledcell indicator, and compute and slice data needed for the neighbor cells.
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::GhostDataToSlice
, which returns a Variables
of the tensors to slice and send to the neighbors. Note that the Tags::Inactive<variables_tag>
are already projected onto the subcells before this function is called. This is because the projection needs to be done anyway for the a posteriori TCI. That is, the evolved variables are projected at the end of the time step rather than the beginning, while this function is called at the beginning. The main reason for having the mutator GhostDataToSlice
is to allow sending primitive or characteristic variables for reconstruction.
bool 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 std::vector< double > &  max_of_past_variables,  
const std::vector< double > &  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}\).
DirectionMap< Dim, std::vector< double > > evolution::dg::subcell::slice_data  (  const Variables< TagList > &  volume_subcell_vars, 
const Index< Dim > &  subcell_extents,  
const size_t  number_of_ghost_points,  
const DirectionMap< Dim, bool > &  directions_to_slice  
) 
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 xvariesfastest).
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  
) 
Slice a single volume tensor for a given direction and slicing depth (number of ghost points).
Note that the last argument has the type Direction<Dim>
, not a DirectionMap (cf. evolution::dg::subcell::slice_data
)
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  
) 
Slice a single volume tensor for a given direction and slicing depth (number of ghost points).
Note that the last argument has the type Direction<Dim>
, not a DirectionMap (cf. evolution::dg::subcell::slice_data
)
Scalar< DataVector > evolution::dg::subcell::tci_status  (  const Mesh< Dim > &  dg_mesh, 
const Mesh< Dim > &  subcell_mesh,  
subcell::ActiveGrid  active_grid,  
const std::deque< subcell::ActiveGrid > &  tci_history  
) 
Set status
to 1
if the cell is marked as needing subcell, otherwise set status
to 0
. status
has grid points on the currently active mesh.
status
of 0
, indicating a smooth solution, even though the active mesh is the subcell mesh. This can occur when using a multistep time integration method like Adams Bashforth, where the FD/FV scheme is used as long as any of the timestepper history is flagged for FD/FV subcell evolution. An example of this would be just after a shock leaves the cell: the solution is now smooth (status
takes value 0
), but the multistep method continues to take FD timesteps (status
is defined on the subcell mesh) until the entire history is flagged as smooth. Thus, tci_history.front()
(which corresponds to the most recent TCI status) is used to set the status when using multistep time integrators, while the active_grid
is used when using other time steppers. void evolution::dg::subcell::tci_status  (  gsl::not_null< Scalar< DataVector > * >  status, 
const Mesh< Dim > &  dg_mesh,  
const Mesh< Dim > &  subcell_mesh,  
subcell::ActiveGrid  active_grid,  
const std::deque< subcell::ActiveGrid > &  tci_history  
) 
Set status
to 1
if the cell is marked as needing subcell, otherwise set status
to 0
. status
has grid points on the currently active mesh.
status
of 0
, indicating a smooth solution, even though the active mesh is the subcell mesh. This can occur when using a multistep time integration method like Adams Bashforth, where the FD/FV scheme is used as long as any of the timestepper history is flagged for FD/FV subcell evolution. An example of this would be just after a shock leaves the cell: the solution is now smooth (status
takes value 0
), but the multistep method continues to take FD timesteps (status
is defined on the subcell mesh) until the entire history is flagged as smooth. Thus, tci_history.front()
(which corresponds to the most recent TCI status) is used to set the status when using multistep time integrators, while the active_grid
is used when using other time steppers. bool evolution::dg::subcell::two_mesh_rdmp_tci  (  const Variables< tmpl::list< EvolvedVarsTags... > > &  dg_evolved_vars, 
const Variables< tmpl::list< Tags::Inactive< EvolvedVarsTags >... > > &  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 variablespecific manner.