SpECTRE  v2024.03.19
Finite Difference

Functions needed for (conservative) finite difference methods. More...

Namespaces

namespace  fd
 Functions and classes for finite difference methods.
 
namespace  fd::reconstruction
 Variable and flux vector splitting reconstruction schemes for finite difference methods.
 

Enumerations

enum class  fd::reconstruction::FallbackReconstructorType { Minmod , MonotonisedCentral , None }
 Possible types of reconstruction routines to fall back to from higher-order reconstruction when using adaptive method. More...
 

Functions

template<size_t NonlinearWeightExponent, size_t Dim>
void fd::reconstruction::aoweno_53 (const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables, const double gamma_hi, const double gamma_lo, const double epsilon)
 Performs an adaptive order (AO) WENO reconstruction using a single 5th order stencil and a 3rd-order CWENO scheme. More...
 
template<size_t Dim>
void fd::reconstruction::minmod (const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables)
 Performs minmod reconstruction on the volume_vars in each direction. More...
 
template<size_t Dim>
void fd::reconstruction::monotonicity_preserving_5 (const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables, const double alpha, const double epsilon)
 Performs the fifth order monotonicity-preserving (MP5) reconstruction [170]. More...
 
template<size_t Dim>
void fd::reconstruction::monotonised_central (const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables)
 Performs monotonised central-difference reconstruction on the vars in each direction. More...
 
template<size_t StencilSize>
std::array< std::array< double, StencilSize >, StencilSize > fd::non_uniform_1d_weights (const std::deque< double > &times)
 Returns the weights for a 1D non-uniform finite difference stencil. More...
 
template<Side LowerOrUpperSide, typename Reconstructor , bool UseExteriorCell = true, size_t NumberOfGhostPoints = (Reconstructor::stencil_width() / 2 + 1), size_t Dim, typename... ArgsForReconstructor>
void fd::reconstruction::reconstruct_neighbor (gsl::not_null< DataVector * > face_data, const DataVector &volume_data, const DataVector &neighbor_data, const Index< Dim > &volume_extents, const Index< Dim > &ghost_data_extents, const Direction< Dim > &direction_to_reconstruct, const ArgsForReconstructor &... args_for_reconstructor)
 In a given direction, reconstruct the cells in the neighboring Element (or cluster of cells) nearest to the shared boundary between the current and neighboring Element (or cluster of cells). More...
 
template<size_t Degree, size_t Dim>
void fd::reconstruction::unlimited (const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables)
 Performs unlimited reconstruction on the vars in each direction. More...
 
template<size_t NonlinearWeightExponent, class FallbackReconstructor , size_t Dim>
void fd::reconstruction::wcns5z (const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables, const double epsilon, const size_t max_number_of_extrema)
 Performs fifth order weighted compact nonlinear scheme reconstruction based on the WENO-Z method (WCNS-5Z). Adaptive fallback combined with an auxiliary reconstruction method (e.g. monotonised central) is also supported. More...
 
template<typename LowOrderReconstructor , bool PositivityPreserving, bool Use9thOrder, bool Use7thOrder, size_t Dim>
void fd::reconstruction::positivity_preserving_adaptive_order (const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables, const double four_to_the_alpha_5, const double six_to_the_alpha_7, const double eight_to_the_alpha_9)
 Performs positivity-preserving adaptive-order FD reconstruction. More...
 
template<typename LowOrderReconstructor , bool PositivityPreserving, bool Use9thOrder, bool Use7thOrder, size_t Dim>
void fd::reconstruction::positivity_preserving_adaptive_order (const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::not_null< std::optional< std::array< gsl::span< std::uint8_t >, Dim > > * > reconstruction_order, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables, const double four_to_the_alpha_5, const double six_to_the_alpha_7, const double eight_to_the_alpha_9)
 Performs positivity-preserving adaptive-order FD reconstruction. More...
 

Detailed Description

Functions needed for (conservative) finite difference methods.

Enumeration Type Documentation

◆ FallbackReconstructorType

Possible types of reconstruction routines to fall back to from higher-order reconstruction when using adaptive method.

Note
FD reconstructions with input arguments or options (e.g. AoWeno, Wcns5z, ..) are currently not supported as fallback types. We may need it at some point in the future but it would require using a different option parsing method (e.g. factory creation).

Function Documentation

◆ aoweno_53()

template<size_t NonlinearWeightExponent, size_t Dim>
void fd::reconstruction::aoweno_53 ( const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_upper_side_of_face_vars,
const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_lower_side_of_face_vars,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Index< Dim > &  volume_extents,
const size_t  number_of_variables,
const double  gamma_hi,
const double  gamma_lo,
const double  epsilon 
)

Performs an adaptive order (AO) WENO reconstruction using a single 5th order stencil and a 3rd-order CWENO scheme.

The AO-WENO(5,3) scheme is based on the scheme presented in [9] but adjusted to do reconstruction on variables instead of fluxes. The Legendre basis functions on the domain \(\xi\in[-1/2,1/2]\) are given by:

\begin{align*} L_0(\xi) &= 1 \\ L_1(\xi) &= \xi \\ L_2(\xi) &= \xi^2-\frac{1}{12} \\ L_3(\xi) &= \xi^3 - \frac{3}{20} \xi \\ L_4(\xi) &= \xi^4 - \frac{3}{14} \xi^2 + \frac{3}{560} \end{align*}

The oscillation indicators are given by

\begin{align*} \beta_l = \sum_{k=1}^{p}\sum_{m=1}^{p}\sigma_{km} u_{k,l} u_{m,l} \end{align*}

where \(p\) is the maximum degree of the basis function used for the stencil \(l\), and

\begin{align*} \sigma_{km}=\sum_{i=1}^{p}\int_{-1/2}^{1/2} \frac{d^i L_k(\xi)}{d\xi^i}\frac{d^i L_m(\xi)}{d\xi^i}d\xi \end{align*}

We write the 3rd-order reconstructed polynomial \(P^{r3}_l(\xi)\) associated with the stencil \(S^{r3}_l\) as

\begin{align*} P^{r3}_l(\xi) = u_0 +u_{\xi} L_1(\xi) + u_{\xi 2} L_2(\xi) \end{align*}

For the stencil \(S^{r3}_1\) we get

\begin{align*} u_0 &= \frac{25}{24} u_{j} -\frac{1}{12} u_{j-1} + \frac{1}{24} u_{j-2} \\ u_{\xi} &= \frac{1}{2}u_{j-2} - 2 u_{j-1} + \frac{3}{2} u_j \\ u_{\xi2} &= \frac{1}{2} u_{j-2} - u_{j-1} + \frac{1}{2} u_j \end{align*}

For the stencil \(S^{r3}_2\) we get

\begin{align*} u_0 &= \frac{1}{24} u_{j-1} + \frac{11}{12} u_{j} + \frac{1}{24} u_{j+1} \\ u_{\xi} &= \frac{1}{2}(u_{j+1} - u_{j-1}) \\ u_{\xi2} &= \frac{1}{2} u_{j-1} - u_{j} + \frac{1}{2} u_{j+1} \end{align*}

For the stencil \(S^{r3}_3\) we get

\begin{align*} u_0 &= \frac{25}{24} u_{j} -\frac{1}{12} u_{j+1} + \frac{1}{24} u_{j+2} \\ u_{\xi} &= -\frac{1}{2}u_{j+2} + 2 u_{j+1} - \frac{3}{2} u_j \\ u_{\xi2} &= \frac{1}{2} u_{j+2} - u_{j+1} + \frac{1}{2} u_j \end{align*}

The oscillation indicator for the 3rd-order stencils is given by

\begin{align*} \beta^{r3}_l = \left(u_{\xi}\right)^2+\frac{13}{3}\left(u_{\xi2}\right)^2 \end{align*}

We write the 5th-order reconstructed polynomial \(P^{r5}(\xi)\) as

\begin{align*} P^{r5}_l(\xi) = u_0 +u_{\xi} L_1(\xi) + u_{\xi 2} L_2(\xi) + u_{\xi3} L_3(\xi) + u_{\xi4} L_4(\xi) \end{align*}

with

\begin{align*} u_0 &= -\frac{17}{5760} u_{j-2} + \frac{77}{1440} u_{j-1} + \frac{863}{960} u_{j} + \frac{77}{1440} u_{j+1} - \frac{17}{5760} u_{j+2} \\ u_{\xi} &= \frac{17}{240} u_{j-2} - \frac{77}{120} u_{j-1} + \frac{77}{120} u_{j+1} - \frac{17}{240} u_{j+2} \\ u_{\xi2} &= -\frac{11}{336} u_{j-2} + \frac{53}{84} u_{j-1} - \frac{67}{56} u_{j} + \frac{53}{84} u_{j+1} - \frac{11}{336} u_{j+2} \\ u_{\xi3} &= -\frac{1}{12} u_{j-2} + \frac{1}{6} u_{j-1} -\frac{1}{6} u_{j+1} + \frac{1}{12} u_{j+2} \\ u_{\xi4} &= \frac{1}{24} u_{j-2} -\frac{1}{6} u_{j-1} +\frac{1}{4} u_{j} - \frac{1}{6}u_{j+1}+\frac{1}{24}u_{j+2} \end{align*}

The oscillation indicator is given by

\begin{align*} \beta^{r5}&=\left(u_{\xi}+\frac{1}{10}u_{\xi3}\right)^2 \\ &+\frac{13}{3}\left(u_{\xi2}+\frac{123}{455}u_{\xi4}\right)^2\\ &+\frac{781}{20}(u_{\xi3})^2+\frac{1421461}{2275}(u_{\xi4})^2 \end{align*}

There are two linear weights \(\gamma_{\mathrm{hi}}\) and \(\gamma_{\mathrm{lo}}\). \(\gamma_{\mathrm{hi}}\) controls how much of the 5th-order stencil is used in smooth regions, while \(\gamma_{\mathrm{lo}}\) controls the linear weight of the central 3rd-order stencil. For larger \(\gamma_{\mathrm{lo}}\), the 3rd-order method is more centrally weighted. The linear weights for the stencils are given by

\begin{align*} \gamma^{r5}&=\gamma_{\mathrm{hi}} \\ \gamma^{r3}_1& = (1-\gamma_{\mathrm{hi}})(1-\gamma_{\mathrm{lo}})/2 \\ \gamma^{r3}_2& = (1-\gamma_{\mathrm{hi}})\gamma_{\mathrm{lo}} \\ \gamma^{r3}_3& = (1-\gamma_{\mathrm{hi}})(1-\gamma_{\mathrm{lo}})/2 \end{align*}

We use the standard nonlinear weights instead of the "Z" weights of [24]

\begin{align*} w^{r5}&=\frac{\gamma^{r5}}{(\beta^{r5}+\epsilon)^q} \\ w^{r3}_l&=\frac{\gamma^{r3}_l}{(\beta^{r3}_l+\epsilon)^q} \end{align*}

where \(\epsilon\) is a small number used to avoid division by zero. The normalized nonlinear weights are denoted by \(\bar{w}^{r5}\) and \(\bar{w}_l^{r3}\). The final reconstructed polynomial \(P(\xi)\) is given by

\begin{align*} P(\xi)&=\frac{\bar{w}^{r5}}{\gamma^{r5}} \left(P^{r5}(\xi)-\gamma^{r3}_1P^{r3}_1(\xi) -\gamma^{r3}_2P^{r3}_2(\xi)-\gamma^{r3}_3P^{r3}_3(\xi)\right) \\ &+\bar{w}^{r3}_1P^{r3}_1(\xi)+\bar{w}^{r3}_2P^{r3}_2(\xi) +\bar{w}^{r3}_3P^{r3}_3(\xi) \end{align*}

The weights \(\gamma_{\mathrm{hi}}\) and \(\gamma_{\mathrm{lo}}\) are typically chosen to be in the range \([0.85,0.95]\).

First alternative oscillation indicators

Instead of integrating over just the cell, we can instead integrate the basis functions over the entire fit interval, \([-5/2,5/2]\). Using this interval for the \(S^{r3}_l\) and the \(S^{r5}\) stencils we get

\begin{align*} \beta^{r3}_l&=(u_{\xi})^2 + \frac{37}{3} (u_{\xi2})^2 \\ \beta^{r5} &=\left(u_{\xi}+\frac{61}{10}u_{\xi3}\right)^2 + \frac{569}{4}u_{\xi3}^2 \\ &+ \frac{1}{8190742}\left(5383 u_{\xi2} + 167158 u_{\xi4}\right)^2 +\frac{4410763}{501474}(u_{\xi2})^2 \\ &=u_{\xi}^2 + \frac{61}{5}u_{\xi}u_{\xi3} + \frac{37}{3}u_{\xi2}^2 + \frac{1538}{7}u_{\xi2}u_{\xi4} \\ &+ \frac{8973}{50}u_{\xi3}^2 + \frac{167158}{49}u_{\xi4}^2 \end{align*}

Note that the indicator is manifestly non-negative, a required property of oscillation indicators. These indicators weight high modes more, which means the scheme is more sensitive to high-frequency features in the solution.

Note
currently it is the alternative indicators that are used. However, an option to control which are used can readily be added, probably best done as a template parameter with if constexpr to avoid conditionals inside tight loops.

◆ minmod()

template<size_t Dim>
void fd::reconstruction::minmod ( const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_upper_side_of_face_vars,
const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_lower_side_of_face_vars,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Index< Dim > &  volume_extents,
const size_t  number_of_variables 
)

Performs minmod reconstruction on the volume_vars in each direction.

On a 1d mesh with spacing \(\Delta x\) we denote the solution at the \(j\)th point by \(u_j\). The minmod reconstruction [154] first computes:

\begin{align} \sigma_j\equiv \textrm{minmod}\left(\frac{u_j-u_{j-1}}{\Delta x}, \frac{u_{j+1}-u_j}{\Delta x}\right) \end{align}

where

\begin{align} \mathrm{minmod}(a,b)= \left\{ \begin{array}{ll} \mathrm{sgn}(a)\min(\lvert a\rvert, \lvert b\rvert) & \mathrm{if} \; \mathrm{sgn}(a)=\mathrm{sgn}(b) \\ 0 & \mathrm{otherwise} \end{array}\right. \end{align}

The reconstructed solution \(u_j(x)\) in the \(j\)th cell is given by

\begin{align} u_j(x)=u_j + \sigma_j (x-x_j), \end{align}

where \(x_j\) is the coordinate \(x\) at the center of the cell.

◆ monotonicity_preserving_5()

template<size_t Dim>
void fd::reconstruction::monotonicity_preserving_5 ( const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_upper_side_of_face_vars,
const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_lower_side_of_face_vars,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Index< Dim > &  volume_extents,
const size_t  number_of_variables,
const double  alpha,
const double  epsilon 
)

Performs the fifth order monotonicity-preserving (MP5) reconstruction [170].

First, calculate the original interface value \(q_{j+1/2}^\text{OR}\) with the (unlimited) fifth order finite difference reconstruction

\begin{align} q_{j+1/2}^\text{OR} = \frac{1}{128}( 3 q_{j-2} - 20 q_{j-1} + 90 q_{j} + 60 q_{j+1} - 5 q_{j+2}) . \end{align}

and compute

\begin{align} q^\text{MP} = q_j + \text{minmod}(q_{j+1} - q_j, \alpha(q_j - q_{j-1})) \end{align}

for a given value of \(\alpha\), which is usually set to \(\alpha=4\).

If \( (q_{j+1/2}^\text{OR} - q_j)(q_{j+1/2}^\text{OR} - q^\text{MP}) \leq \epsilon\) where \(\epsilon\) is a small tolerance value, use \(q_{1/2} = q_{j+1/2}^\text{OR}\).

Note
A proper value of \(\epsilon\) may depend on the scale of the quantity \(q\) to be reconstructed; reference [170] suggests 1e-10. For hydro simulations with atmosphere treatment, setting \(\epsilon=0.0\) would be safe.

Otherwise, calculate

\begin{align} d_{j+1} & = q_{j+2} - 2q_{j+1} + q_{j} \\ d_j & = q_{j+1} - 2q_j + q_{j-1} \\ d_{j-1} & = q_{j} - 2q_{j-1} + q_{j-2} , \end{align}

\begin{align} d^\text{M4}_{j+1/2} = \text{minmod}(4d_j - d_{j+1}, 4d_{j+1}-d_j, d_j, d_{j+1}) \\ d^\text{M4}_{j-1/2} = \text{minmod}(4d_j - d_{j-1}, 4d_{j-1}-d_j, d_j, d_{j-1}), \end{align}

\begin{align} q^\text{UL} & = q_j + \alpha(q_j - q_{j-1}) \\ q^\text{AV} & = (q_j + q_{j+1}) / 2 \\ q^\text{MD} & = q^\text{AV} - d^\text{M4}_{j+1/2} / 2 \\ q^\text{LC} & = q_j + (q_j - q_{j-1}) / 2 + (4/3) d^\text{M4}_{j-1/2} \end{align}

and

\begin{align} q^\text{min} & = \text{max}[ \text{min}(q_j, q_{j+1}, q^\text{MD}), \text{min}(q_j, q^\text{UL}, q^\text{LC}) ] \\ q^\text{max} & = \text{min}[ \text{max}(q_j, q_{j+1}, q^\text{MD}), \text{max}(q_j, q^\text{UL}, q^\text{LC}) ]. \end{align}

The reconstructed value is given as

\begin{align} q_{i+1/2} = \text{median}(q_{j+1/2}^\text{OR}, q^\text{min}, q^\text{max}) \end{align}

where the median function can be expressed as

\begin{align} \text{median}(x,y,z) = x + \text{minmod}(y-x, z-x). \end{align}

◆ monotonised_central()

template<size_t Dim>
void fd::reconstruction::monotonised_central ( const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_upper_side_of_face_vars,
const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_lower_side_of_face_vars,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Index< Dim > &  volume_extents,
const size_t  number_of_variables 
)

Performs monotonised central-difference reconstruction on the vars in each direction.

On a 1d mesh with spacing \(\Delta x\) we denote the solution at the \(j\)th point by \(u_j\). The monotonised central-difference reconstruction [154] first computes:

\begin{align} \sigma_j\equiv \textrm{minmod} \left(\frac{u_{j+1} - u_{j-1}}{2\Delta x}, 2\frac{u_j-u_{j-1}}{\Delta x}, 2\frac{u_{j+1}-u_j}{\Delta x}\right) \end{align}

where

\begin{align} \mathrm{minmod}(a,b,c)= \left\{ \begin{array}{ll} \mathrm{sgn}(a)\min(\lvert a\rvert, \lvert b\rvert, \lvert c\rvert) & \mathrm{if} \; \mathrm{sgn}(a)=\mathrm{sgn}(b)=\mathrm{sgn}(c) \\ 0 & \mathrm{otherwise} \end{array}\right. \end{align}

The reconstructed solution \(u_j(x)\) in the \(j\)th cell is given by

\begin{align} u_j(x)=u_j + \sigma_j (x-x_j) \end{align}

where \(x_j\) is the coordinate \(x\) at the center of the cell.

◆ non_uniform_1d_weights()

template<size_t StencilSize>
std::array< std::array< double, StencilSize >, StencilSize > fd::non_uniform_1d_weights ( const std::deque< double > &  times)

Returns the weights for a 1D non-uniform finite difference stencil.

These weights are for the Lagrange interpolation polynomial and its derivatives evaluated at times[0]. If the number of times is not the same as StencilSize, then an error will occur.

The reason that times was chosen to be monotonically decreasing is because the intended use of this function is with data that is stored in a std::deque where the zeroth element is the "most recent" in time and then later elements are further in the "past".

Parameters
timesA monotonically decreasing sequence of times

Returns: A 2D array of all finite difference weights that correspond to the input times. The outer dimension is the Nth derivative, and the inner dimension loops over the weights for each of the times.

Note
Only stencil sizes 2, 3, and 4 are implemented right now. If you need more, you'll either need to add it yourself or generalize the algorithm.

◆ positivity_preserving_adaptive_order() [1/2]

template<typename LowOrderReconstructor , bool PositivityPreserving, bool Use9thOrder, bool Use7thOrder, size_t Dim>
void fd::reconstruction::positivity_preserving_adaptive_order ( const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_upper_side_of_face_vars,
const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_lower_side_of_face_vars,
const gsl::not_null< std::optional< std::array< gsl::span< std::uint8_t >, Dim > > * >  reconstruction_order,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Index< Dim > &  volume_extents,
const size_t  number_of_variables,
const double  four_to_the_alpha_5,
const double  six_to_the_alpha_7,
const double  eight_to_the_alpha_9 
)

Performs positivity-preserving adaptive-order FD reconstruction.

Performs a fifth-order unlimited reconstruction. If the reconstructed values at the interfaces aren't positive (when PositivityPreserving is true) or when the Persson TCI condition:

\begin{align} 4^\alpha \int_{x_{i-5/2}}^{x_{i+5/2}} \hat{u}^2(x) dx > \int_{x_{i-5/2}}^{x_{i+5/2}} u^2(x) dx \end{align}

is satisfied, where \(\hat{u}\) is the polynomial with only the largest modal coefficient non-zero, then the LowOrderReconstructor is used.

If PositivityPreserving is true then if the low-order reconstructed solution isn't positive we use first-order (constant in space) reconstruction.

If Use9thOrder is true then first a ninth-order reconstruction is used, followed by fifth-order. If Use7thOrder is true then seventh-order reconstruction is used before fifth-order (but after ninth-order if Use9thOrder is also true). This allows using the highest possible order locally for reconstruction.

Fifth order unlimited reconstruction is:

\begin{align} \hat{u}_{i+1/2}=\frac{3}{128}u_{i-2} - \frac{5}{32}u_{i-1} + \frac{45}{64}u_{i} + \frac{15}{32}u_{i+1} - \frac{5}{128}u_{i+2} \end{align}

Seventh order unlimited reconstruction is:

\begin{align} \hat{u}_{i+1/2}&=-\frac{5}{1024}u_{i-3} + \frac{21}{512}u_{i-2} - \frac{175}{1024}u_{i-1} + \frac{175}{256}u_{i} \\ &+ \frac{525}{1024}u_{i+1} - \frac{35}{512}u_{i+2} + \frac{7}{1024}u_{i+3} \end{align}

Ninth order unlimited reconstruction is:

\begin{align} \hat{u}_{i+1/2}&=\frac{35}{32768}u_{i-4} - \frac{45}{4096}u_{i-3} + \frac{441}{8291}u_{i-2} - \frac{735}{4096}u_{i-1} \\ &+ \frac{11025}{16384}u_{i} + \frac{2205}{4096}u_{i+1} - \frac{735}{8192}u_{i+2} + \frac{63}{4096}u_{i+3} \\ &- \frac{45}{32768}u_{i+4} \end{align}

◆ positivity_preserving_adaptive_order() [2/2]

template<typename LowOrderReconstructor , bool PositivityPreserving, bool Use9thOrder, bool Use7thOrder, size_t Dim>
void fd::reconstruction::positivity_preserving_adaptive_order ( const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_upper_side_of_face_vars,
const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_lower_side_of_face_vars,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Index< Dim > &  volume_extents,
const size_t  number_of_variables,
const double  four_to_the_alpha_5,
const double  six_to_the_alpha_7,
const double  eight_to_the_alpha_9 
)

Performs positivity-preserving adaptive-order FD reconstruction.

Performs a fifth-order unlimited reconstruction. If the reconstructed values at the interfaces aren't positive (when PositivityPreserving is true) or when the Persson TCI condition:

\begin{align} 4^\alpha \int_{x_{i-5/2}}^{x_{i+5/2}} \hat{u}^2(x) dx > \int_{x_{i-5/2}}^{x_{i+5/2}} u^2(x) dx \end{align}

is satisfied, where \(\hat{u}\) is the polynomial with only the largest modal coefficient non-zero, then the LowOrderReconstructor is used.

If PositivityPreserving is true then if the low-order reconstructed solution isn't positive we use first-order (constant in space) reconstruction.

If Use9thOrder is true then first a ninth-order reconstruction is used, followed by fifth-order. If Use7thOrder is true then seventh-order reconstruction is used before fifth-order (but after ninth-order if Use9thOrder is also true). This allows using the highest possible order locally for reconstruction.

Fifth order unlimited reconstruction is:

\begin{align} \hat{u}_{i+1/2}=\frac{3}{128}u_{i-2} - \frac{5}{32}u_{i-1} + \frac{45}{64}u_{i} + \frac{15}{32}u_{i+1} - \frac{5}{128}u_{i+2} \end{align}

Seventh order unlimited reconstruction is:

\begin{align} \hat{u}_{i+1/2}&=-\frac{5}{1024}u_{i-3} + \frac{21}{512}u_{i-2} - \frac{175}{1024}u_{i-1} + \frac{175}{256}u_{i} \\ &+ \frac{525}{1024}u_{i+1} - \frac{35}{512}u_{i+2} + \frac{7}{1024}u_{i+3} \end{align}

Ninth order unlimited reconstruction is:

\begin{align} \hat{u}_{i+1/2}&=\frac{35}{32768}u_{i-4} - \frac{45}{4096}u_{i-3} + \frac{441}{8291}u_{i-2} - \frac{735}{4096}u_{i-1} \\ &+ \frac{11025}{16384}u_{i} + \frac{2205}{4096}u_{i+1} - \frac{735}{8192}u_{i+2} + \frac{63}{4096}u_{i+3} \\ &- \frac{45}{32768}u_{i+4} \end{align}

◆ reconstruct_neighbor()

template<Side LowerOrUpperSide, typename Reconstructor , bool UseExteriorCell = true, size_t NumberOfGhostPoints = (Reconstructor::stencil_width() / 2 + 1), size_t Dim, typename... ArgsForReconstructor>
void fd::reconstruction::reconstruct_neighbor ( gsl::not_null< DataVector * >  face_data,
const DataVector volume_data,
const DataVector neighbor_data,
const Index< Dim > &  volume_extents,
const Index< Dim > &  ghost_data_extents,
const Direction< Dim > &  direction_to_reconstruct,
const ArgsForReconstructor &...  args_for_reconstructor 
)

In a given direction, reconstruct the cells in the neighboring Element (or cluster of cells) nearest to the shared boundary between the current and neighboring Element (or cluster of cells).

This is needed if one is sending reconstruction and flux data separately, or if one is using DG-FD hybrid schemes. Below is an ASCII diagram of what is reconstructed.

Self | Neighbor
x x x x x | o o o
^+
Reconstruct to right/+ side of the interface

◆ unlimited()

template<size_t Degree, size_t Dim>
void fd::reconstruction::unlimited ( const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_upper_side_of_face_vars,
const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_lower_side_of_face_vars,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Index< Dim > &  volume_extents,
const size_t  number_of_variables 
)

Performs unlimited reconstruction on the vars in each direction.

On a 1d mesh with spacing \(\Delta x\) we denote the solution at the \(j\)th point by \(u_j\). The degree 2 reconstruction is

\begin{align} u_{j+1/2} &= -\frac{1}{8} u_{j-1} + \frac{3}{4} u_j + \frac{3}{8} u_{j+1}, \\ u_{j-1/2} &= \frac{3}{8} u_{j-1} + \frac{3}{4} u_j - \frac{1}{8} u_{j+1}. \end{align}

The degree 4 reconstruction is

\begin{align} u_{j-1/2} &= - \frac{5}{128}u_{j-2} + \frac{15}{32}u_{j-1} + \frac{45}{64}u_{j} -\frac{5}{32}u_{j+1} + \frac{3}{128}u_{j+2}, \\ u_{j+1/2} &= \frac{3}{128}u_{j-2} - \frac{5}{32}u_{j-1} + \frac{45}{64}u_{j} + \frac{15}{32}u_{j+1} - \frac{5}{128}u_{j+2}. \end{align}

Degree 6 and 8 reconstruction is also supported.

◆ wcns5z()

template<size_t NonlinearWeightExponent, class FallbackReconstructor , size_t Dim>
void fd::reconstruction::wcns5z ( const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_upper_side_of_face_vars,
const gsl::not_null< std::array< gsl::span< double >, Dim > * >  reconstructed_lower_side_of_face_vars,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Index< Dim > &  volume_extents,
const size_t  number_of_variables,
const double  epsilon,
const size_t  max_number_of_extrema 
)

Performs fifth order weighted compact nonlinear scheme reconstruction based on the WENO-Z method (WCNS-5Z). Adaptive fallback combined with an auxiliary reconstruction method (e.g. monotonised central) is also supported.

Using the WENO oscillation indicators given by [91]

\begin{align} \beta_0 & = \frac{13}{12} (q_{i-2} - 2q_{i-1} + q_{i})^2 + \frac{1}{4} (q_{i-2} - 4q_{i-1} + 3q_{i})^2 \\ \beta_1 & = \frac{13}{12} (q_{i-1} - 2q_{i} + q_{i+1})^2 + \frac{1}{4} (q_{i+1} - q_{i-1})^2 \\ \beta_2 & = \frac{13}{12} (q_{i} - 2q_{i+1} + q_{i+2})^2 + \frac{1}{4} (3q_{i} - 4q_{i+1} + q_{i+2})^2 , \end{align}

compute the modified nonlinear weights (cf. WENO-Z scheme: see Eq. 42 of [24])

\begin{align} \omega_k^z = \frac{\alpha_k^z}{\sum_{l=0}^z \alpha_l^z}, \quad \alpha_k^z = c_k \left(1 + \left( \frac{\tau_5}{\beta_k + \epsilon_k} \right)^q \right), \quad k = 0, 1, 2 \end{align}

where \(\epsilon_k\) is a factor used to avoid division by zero and is set to

\begin{align} \epsilon_k = \varepsilon (1 + |q_{i+k}| + |q_{i+k-1}| + |q_{i+k-2}|), \quad k = 0, 1, 2, \end{align}

linear weights \(c_k\) are adopted from [134] (see Table 17 of it)

\begin{align} c_0 = 1/16, \quad c_1 = 10/16, \quad c_2 = 5/16, \end{align}

and \(\tau_5 \equiv |\beta_2-\beta_0|\).

Reconstruction stencils are given by Lagrange interpolations (e.g. see Eq. 29d - 29f of [28])

\begin{align} q_{i+1/2}^0 &= \frac{3}{8}q_{i-2} -\frac{5}{4}q_{i-1} +\frac{15}{8}q_{i} \\ q_{i+1/2}^1 &= -\frac{1}{8}q_{i-1} +\frac{3}{4}q_{i} +\frac{3}{8}q_{i+1} \\ q_{i+1/2}^2 &= \frac{3}{8}q_{i} +\frac{3}{4}q_{i+1} -\frac{1}{8}q_{i+2} \end{align}

and the final reconstructed solution is given by

\begin{align} q_{i+1/2} = \sum_{k=0}^2 \omega_k q_{i+1/2}^k . \end{align}

The nonlinear weights exponent \(q (=1 \text{ or } 2)\) and the small factor \(\varepsilon\) can be chosen via an input file.

If the template parameter FallbackReconstructor is set to one of the FD reconstructor structs of which names are listed in fd::reconstruction::FallbackReconstructorType, adaptive reconstruction is performed as follows. For each finite difference stencils, first check how many extrema are in the stencil. If the number of local extrema is less than or equal to a non-negative integer max_number_of_extrema which is given as an input parameter, perform the WCNS-5Z reconstruction; otherwise switch to the given FallbackReconstructor for performing reconstruction. If FallbackReconstructor is set to void, the adaptive method is disabled and WCNS-5Z is always used.