Namespaces | Classes | Enumerations | Functions
evolution::dg::subcell Namespace Reference

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

Namespaces

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

Classes

struct  NeighborData
 Holds neighbor data needed for the TCI and reconstruction. More...
 
class  SubcellOptions
 Holds the system-agnostic subcell parameters, such as numbers controlling when to switch between DG and subcell. More...
 

Enumerations

enum  ActiveGrid { Dg, Subcell }
 

Functions

std::ostreamoperator<< (std::ostream &os, ActiveGrid active_grid) noexcept
 
void pup (PUP::er &p, NeighborData &nhbr_data) noexcept
 
void operator| (PUP::er &p, NeighborData &nhbr_data) noexcept
 
bool operator== (const NeighborData &lhs, const NeighborData &rhs) noexcept
 
bool operator!= (const NeighborData &lhs, const NeighborData &rhs) noexcept
 
std::ostreamoperator<< (std::ostream &os, const NeighborData &t) noexcept
 
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 double zero_cutoff) noexcept
 Troubled cell indicator using spectral falloff of [71]. 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) noexcept
 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) noexcept
 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.
 
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) noexcept
 Slice the subcell variables needed for neighbors to perform reconstruction. More...
 
bool operator== (const SubcellOptions &lhs, const SubcellOptions &rhs) noexcept
 
bool operator!= (const SubcellOptions &lhs, const SubcellOptions &rhs) noexcept
 
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) noexcept
 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< DataVectortci_status (const Mesh< Dim > &dg_mesh, const Mesh< Dim > &subcell_mesh, subcell::ActiveGrid active_grid, const std::deque< subcell::ActiveGrid > &tci_history) noexcept
 
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) noexcept
 Troubled cell indicator using a relaxed discrete maximum principle, comparing the solution on two grids at the same point in time. More...
 

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 [30]. Other implementations of a subcell limiter exist, e.g. [81] [17] [47]. Our implementation and that of [30] are a generalization of the Multidimensional Optimal Order Detection (MOOD) algorithm [19] [27] [28] [60].

Function Documentation

◆ 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 double  zero_cutoff 
)
noexcept

Troubled cell indicator using spectral falloff of [71].

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)=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.

◆ rdmp_tci()

template<typename... EvolvedVarsTags>
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 
)
noexcept

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}\).

◆ slice_data()

template<size_t Dim, typename TagList >
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 
)
noexcept

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).

◆ tci_status()

template<size_t Dim>
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 
)
noexcept

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.

Note
It is possible to encounter a 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.

◆ two_mesh_rdmp_tci()

template<typename... EvolvedVarsTags>
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 
)
noexcept

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.