|
SpECTRE
v2024.08.03
|
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 [178]. 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 > ×) |
| 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... | |
Functions needed for (conservative) finite difference methods.
|
strong |
Possible types of reconstruction routines to fall back to from higher-order reconstruction when using adaptive method.
| 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]\).
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.
if constexpr to avoid conditionals inside tight loops. | 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 [161] 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.
| 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 [178].
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}\).
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}
| 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 [161] 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.
| 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".
| times | A 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.
| 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}
| 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}
| 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.
| 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.
| 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 [95]
\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 [140] (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.