Namespaces | Functions
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.
 

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) noexcept
 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) noexcept
 Performs minmod reconstruction on the volume_vars in each direction. More...
 
template<size_t Dim>
void fd::reconstruction::monotised_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) noexcept
 Performs monotised 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) noexcept
 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...
 

Detailed Description

Functions needed for (conservative) finite difference methods.

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

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 [16]

\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 
)
noexcept

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 [95] 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.

◆ monotised_central()

template<size_t Dim>
void fd::reconstruction::monotised_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 
)
noexcept

Performs monotised 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 monotised central-difference reconstruction [95] 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 
)
noexcept

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