SpECTRE  v2022.06.14
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::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<Side LowerOrUpperSide, typename Reconstructor , 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, 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)
 Performs fifth order weighted compact nonlinear scheme reconstruction based on the WENO-Z method. More...
 
template<size_t NonlinearWeightExponent, size_t Dim>
void fd::reconstruction::wcns5z_mc (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 size_t max_number_of_extrema, const double epsilon)
 Performs adaptive reconstruction using WCNS-5Z and MC schemes. 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 [5] 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 [17]

\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 [107] 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.

◆ 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 [107] 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.

◆ reconstruct_neighbor()

template<Side LowerOrUpperSide, typename Reconstructor , 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}

◆ wcns5z()

template<size_t NonlinearWeightExponent, 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 
)

Performs fifth order weighted compact nonlinear scheme reconstruction based on the WENO-Z method.

Using the WENO oscillation indicators given by [68]

\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 [17])

\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 [95] (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 [20])

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

◆ wcns5z_mc()

template<size_t NonlinearWeightExponent, size_t Dim>
void fd::reconstruction::wcns5z_mc ( 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 size_t  max_number_of_extrema,
const double  epsilon 
)

Performs adaptive reconstruction using WCNS-5Z and MC schemes.

For each finite difference stencils first check how many extrema are in a given stencil. If the number of extrema is less than or equal to a non-negative integer max_number_of_extrema which is given as an input parameter, perform fifth order weighted compact nonlinear reconstruction with Z oscillation indicator (WCNS-5Z); otherwise, switch to the monotonised central (MC) reconstruction.

See the documentations of wcns5z() and monotonised_central() for descriptions on each reconstruction schemes.