Reconstruction Schemes

Pydro provides two variants of reconstruction schemes. The first are the standard reconstruction schemes including TVD and WENO-type schemes, while the second are positivity-preserving adaptive-order schemes. The two scheme classes have different interfaces and so must be used slightly differently. The interface for each is documented in the sections below.

Standard reconstruction

The standard reconstruction schemes such as total variation diminishing (TVD) schemes and different flavors of higher-order essentially non-oscillatory or weighted compact schemes are available through a common interface. The Reconstruction.reconstruct() allows reconstructing a series of variables using one of the available schemes (see Reconstruction.Scheme for a list).

class Reconstruction.Scheme

An enum of the various different reconstruction schemes that are supported.

Minmod = 1

Minmod reconstruction.

Minmod reconstruction is performed as

\[\begin{align} \sigma_j=\mathrm{minmod} \left(\frac{q_i-q_{i-1}}{\Delta\xi}, \frac{q_{i+1}-q_i}{\Delta\xi}\right). \end{align}\]

where \(\Delta\xi\) is the grid spacing and \(\mathrm{minmod}(a,b)\) is defined as

\[\begin{split} \begin{align} &\mathrm{minmod}(a,b)= \notag \\ &\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}\end{split}\]

The reconstructed solution at the faces is given by

\[\hat{q}_{i+1/2} = q_i +\frac{\Delta\xi}{2}\sigma_i\]

See, e.g. section 9.3.1 of [2] for a discussion.

Wcns3 = 2

Third order weighted compact nonlinear scheme reconstruction [3].

Third order WCNS3 reconstruction is done by first defining oscillation indicators \(\beta_0\) and \(\beta_1\) as

\[\begin{split}\begin{align} \beta_0 &= (q_i - q_{i-1})^2 \\ \beta_1 &= (q_{i+1} - q_{i})^2 \end{align}\end{split}\]

Then coefficients \(\alpha_k\) are defined as

\[\alpha_k = \frac{c_k}{(\beta_k + \epsilon_k)^2}\]

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

\[\begin{split}\begin{align} \epsilon_0 &= 10^{-17}\left(1 + |q_{i}| + |q_{i-1}|\right) \\ \epsilon_1 &= 10^{-17}\left(1 + |q_{i}| + |q_{i+1}|\right) \end{align}\end{split}\]

and the linear weights are \(c_0=1/4\) and \(c_1=3/4\). Finally, we define the nonlinear weights:

\[\omega_k=\frac{\alpha_k}{\sum_{k=0}^{1}\alpha_k}\]

The reconstruction stencils are given by:

\[\begin{split}\begin{align} q^0_{i+1/2}&=\frac{3}{2}q_i-\frac{1}{2}q_{i-1} \\ q^1_{i+1/2}&=\frac{1}{2}q_i+\frac{1}{2}q_{i+1} \end{align}\end{split}\]

The final reconstructed solution is given by

\[\hat{q}_{i+1/2}=\sum_{k=0}^{1}\omega_k q^k_{i+1/2}\]

Wcns5 = 4

Fifth order weighted compact nonlinear scheme reconstruction [1].

The oscillation indicators are given by

\[\begin{split}\begin{align} \beta_0 &= \frac{1}{4}\left(q_{i-2}-4 q_{i-1}+3 q_{i}\right)^2 + \left(q_{i-2}-2 q_{i-1} + q_{i}\right)^2 \\ \beta_1 &= \frac{1}{4}\left(q_{i-1} - q_{i+1}\right)^2 + \left(q_{i-1} - 2 q_{i+1}\right)^2 \\ \beta_2 &= \frac{1}{4}\left(3 q_{i}-4 q_{i+1}+q_{i+2}\right)^2 + \left(q_{i} - 2 q_{i+1} + q_{i+2}\right)^2 \end{align}\end{split}\]

Then coefficients \(\alpha_k\) are defined as

\[\alpha_k = \frac{c_k}{(\beta_k + \epsilon_k)^2}\]

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

\[\begin{split}\begin{align} \epsilon_0 &= 2\times10^{-16}\left(1 + |q_{i}| + |q_{i-1}| + |q_{i-2}|\right) \\ \epsilon_1 &= 2\times10^{-16}\left(1 + |q_{i}| + |q_{i+1}| + |q_{i-1}|\right) \\ \epsilon_2 &= 2\times10^{-16}\left(1 + |q_{i}| + |q_{i+1}| + |q_{i+2}|\right) \end{align}\end{split}\]

and the linear weights are \(c_0=1/16, c_1=10/16\), and \(c_2=5/16\). Finally, we define the nonlinear weights:

\[\omega_k=\frac{\alpha_k}{\sum_{k=0}^{2}\alpha_k}\]

The reconstruction stencils are given by:

\[\begin{split}\begin{align} q^0_{i+1/2}&=\frac{3}{8}q_{i-2} - \frac{5}{4}q_{i-1} + \frac{15}{8}q_i, \\ q^1_{i+1/2}&=-\frac{1}{8}q_{i-1} + \frac{3}{4}q_i + \frac{3}{8}q_{i+1}, \\ q^2_{i+1/2}&=\frac{3}{8}q_i + \frac{3}{4}q_{i+1} - \frac{1}{8}q_{i+2}. \end{align}\end{split}\]

The final reconstructed solution is given by

\[\hat{q}_{i+1/2}=\sum_{k=0}^{2}\omega_k q^k_{i+1/2}\]

Wcns5Weno = 6

Fifth order weighted compact nonlinear scheme reconstruction with the Jiang and Shu [4] weights.

Follows the procedure of Wcns5() except using the oscillation indicators given by

\[\begin{split}\begin{align} \beta_0 &=\frac{1}{4}\left(q_{i-2}-4q_{i-1}+3q_i\right)^2 +\frac{13}{12}\left(q_{i-2}-2q_{i-1}+q_i\right)^2 \\ \beta_1 &=\frac{1}{4}\left(q_{i-1}-q_{i+1}\right)^2 +\frac{13}{12}\left(q_{i-1} - 2 q_{i+1}\right)^2 \\ \beta_2 &=\frac{1}{4}\left(-3q_i+4q_{i+1}-q_{i+2}\right)^2 +\frac{13}{12}\left(q_i-2q_{i+1}+q_{i+2}\right)^2. \end{align}\end{split}\]

Wcns5z = 5

Fifth order weighted compact nonlinear scheme reconstruction with the \(Z\) oscillation indicator.

Follows the procedure of Wcns5() except using the oscillation indicators given by

\[\beta_k^Z =\frac{\beta_k+\epsilon_k}{\beta_k + \tau_5 + \epsilon_k}\]

where

\[\tau_5 = |\beta_2 - \beta_0|\]

and the oscillation indicators are the ones from Jiang and Shu [4], as described in Wcns5Weno().

Weno3 = 3

Third order weighted essentially non-oscillarity reconstruction.

The same as the Wcns3() reconstruction except with \(c_0=1/3\) and \(c_1=2/3\).

Reconstruction.reconstruct(vars_to_reconstruct, scheme, order_used)

Reconstructs all variables using the requested scheme.

Parameters

• vars_to_reconstruct (list of list of double) – The variables at the cell centers.

• scheme (Reconstruction.Scheme) – The reconstruction scheme to use.

• order_used (list of int) – Filled by the function and is used to return the order of the reconstruction used.

Returns

(list of list of double) The face reconstructed variables. Each variable is of length 2 * number_of_cells

Positivity-preserving reconstruction

The PPAO schemes need to be wrapped in a function specific for each evolution system so that the appropriate variables can have their positivity preserved.

ReconstructionPpao.adaptive_order_1(q, i, j, recons)

First-order reconstruction.

First-order reconstruction is given by

\[\hat{q}_{i + 1/2} = q_i\]

ReconstructionPpao.adaptive_order_3(q, i, j, recons, keep_positive, alpha=3.0, eps=1e-36)

Uses a third-order centered stencil for reconstruction

\[\hat{q}_{i+1/2}=\frac{1}{8}q_{i-1} + \frac{3}{4}q_{i} + \frac{3}{8}q_{i+1}\]

How oscillatory the resulting polynomial is can be determined by comparing

\[s^j_N = \frac{1}{2}\log_{10}\left(\frac{\bar{\kappa}_N} {\kappa_N + \epsilon}\right),\]

where

\[\begin{split}\begin{align} \bar{\kappa}_3 &= \frac{2}{5}\left(\frac{3}{4}q_{i+1} - \frac{3}{2}q_{i} + \frac{3}{4}q_{i-1}\right)^{2}, \\ \kappa_3 &= \left(\frac{3}{8}q_{i+1} + \frac{1}{4}q_{i} + \frac{3}{8}q_{i-1}\right) \left(\frac{31}{20}q_{i+1} - \frac{1}{10}q_{i} + \frac{11}{20}q_{i-1}\right), \end{align}\end{split}\]

to

\[-\alpha\log_{10}(4)\]

Typically \(\alpha\sim4\) so that the coefficients decay as \(1/3^4\).

Parameters

• q (list of double) – The variable values at the cell centers.

• i (int) – The index into the reconstructed array

• j (int) – The index of the cell whose faces are being reconstructed in q

• recons (list of double) – The array of the reconstructed variable.

• keep_positive (bool) – If True then returns False if the reconstructed solution is not positive.

• alpha (double) – The expected decay of increasing coefficients in the method.

• eps (double) – The epsilon parameter to ignore small values and impose an absolute tolerance.

Returns

(bool) True if the reconstruction was successful, otherwise False

ReconstructionPpao.adaptive_order_5(q, i, j, recons, keep_positive, alpha=5.0, eps=1e-36)

Uses a fifth-order centered stencil for reconstruction

\[\hat{q}_{i+1/2}=\frac{3}{128}q_{i-2} - \frac{5}{32}q_{i-1} + \frac{45}{64}q_{i} + \frac{15}{32}q_{i+1} - \frac{5}{128}q_{i+2}\]

Parameters

• q (list of double) – The variable values at the cell centers.

• i (int) – The index into the reconstructed array

• j (int) – The index of the cell whose faces are being reconstructed in q

• recons (list of double) – The array of the reconstructed variable.

• keep_positive (bool) – If True then returns False if the reconstructed solution is not positive.

• alpha (double) – The expected decay of increasing coefficients in the method.

• eps (double) – The epsilon parameter to ignore small values and impose an absolute tolerance.

Returns

(bool) True if the reconstruction was successful, otherwise False

ReconstructionPpao.adaptive_order_7(q, i, j, recons, alpha=5.0, eps=1e-36)

Uses a seventh-order centered stencil for reconstruction

\[\begin{split}\begin{align} \hat{q}_{i+1/2}&=-\frac{5}{1024}q_{i-3} + \frac{21}{512}q_{i-2} - \frac{175}{1024}q_{i-1} + \frac{175}{256}q_{i} \\ &+ \frac{525}{1024}q_{i+1} - \frac{35}{512}q_{i+2} + \frac{7}{1024}q_{i+3} \end{align}\end{split}\]

Parameters

• q (list of double) – The variable values at the cell centers.

• i (int) – The index into the reconstructed array

• j (int) – The index of the cell whose faces are being reconstructed in q

• recons (list of double) – The array of the reconstructed variable.

• keep_positive (bool) – If True then returns False if the reconstructed solution is not positive.

• alpha (double) – The expected decay of increasing coefficients in the method.

• eps (double) – The epsilon parameter to ignore small values and impose an absolute tolerance.

Returns

(bool) True if the reconstruction was successful, otherwise False

ReconstructionPpao.adaptive_order_9(q, i, j, recons, alpha=5.0, eps=1e-36)

Uses a ninth-order centered stencil for reconstruction

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

Parameters

• q (list of double) – The variable values at the cell centers.

• i (int) – The index into the reconstructed array

• j (int) – The index of the cell whose faces are being reconstructed in q

• recons (list of double) – The array of the reconstructed variable.

• keep_positive (bool) – If True then returns False if the reconstructed solution is not positive.

• alpha (double) – The expected decay of increasing coefficients in the method.

• eps (double) – The epsilon parameter to ignore small values and impose an absolute tolerance.

Returns

(bool) True if the reconstruction was successful, otherwise False

ReconstructionPpao.adaptive_order_wcns3(q, i, j, recons, keep_positive, eps=1e-17, c0=0.25, c1=0.75, liu_indicators=True, exponent=2)

A general third order weight compact nonlinear scheme.

Third order WCNS3 reconstruction is done by first defining oscillation indicators \(\beta_0\) and \(\beta_1\) as

\[\begin{split}\begin{align} \beta_0 &= (q_i - q_{i-1})^2 \\ \beta_1 &= (q_{i+1} - q_{i})^2 \end{align}\end{split}\]

We refer to these as the standard oscillation indicators, but also provide the improved oscillation indicators of Liu [5]:

\[\begin{split}\begin{align} \beta_0 &= \frac{1}{4}\left(\lvert q_{i+1} - q_{i-1}\rvert - \lvert4 q_i - 3 q_{i-1} - q_{i+1}\rvert\right)^2, \\ \beta_1 &= \frac{1}{4}\left(\lvert q_{i+1} - q_{i-1}\rvert - \lvert3 q_{i+1} + q_{i-1} - 4 q_{i}\rvert\right)^2. \end{align}\end{split}\]

Then coefficients \(\alpha_k\) are defined as

\[\alpha_k = \frac{c_k}{(\beta_k + \epsilon_k)^p}\]

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

\[\begin{split}\begin{align} \epsilon_0 &= \epsilon\left(1 + |q_{i}| + |q_{i-1}|\right) \\ \epsilon_1 &= \epsilon\left(1 + |q_{i}| + |q_{i+1}|\right) \end{align}\end{split}\]

and the linear weights are \(c_0=1/4\) and \(c_1=3/4\). Finally, we define the nonlinear weights:

\[\omega_k=\frac{\alpha_k}{\sum_{k=0}^{1}\alpha_k}\]

The reconstruction stencils are given by:

\[\begin{split}\begin{align} q^0_{i+1/2}&= \frac{3}{2}q_{i} - \frac{1}{2}q_{i-1}, \\ q^1_{i+1/2}&=\frac{1}{2}q_i + \frac{1}{2}q_{i+1}, \end{align}\end{split}\]

The final reconstructed solution is given by

\[\hat{q}_{i+1/2}=\sum_{k=0}^{1}\omega_k q^k_{i+1/2}\]

Parameters

• q (list of double) – The variable values at the cell centers.

• i (int) – The index into the reconstructed array

• j (int) – The index of the cell whose faces are being reconstructed in q

• recons (list of double) – The array of the reconstructed variable.

• keep_positive (bool) – If True then returns False if the reconstructed solution is not positive.

• eps (double) – The epsilon parameter to avoid division by zero.

• c0 (double) – The optimal linear weight \(c_{0}\). For 3rd order use \(1/4\). For 2nd order but increased robustness use \(1/2\).

• c1 (double) – The optimal linear weight \(c_{1}\). For 3rd order use \(3/4\). For 2nd order but increased robustness use \(1/2\).

• liu_indicators (bool) – If True use the oscillation indicators of [5]

• exponent (int) – The exponent \(p\) in denominator of the \(\alpha_k\)

Returns

(bool) True if the reconstruction was successful, otherwise False

ReconstructionPpao.adaptive_order_wcns3z(q, i, j, recons, keep_positive, eps=1e-17, c0=0.25, c1=0.75, liu_indicators=True)

A general third order weight compact nonlinear scheme using the Z weights.

The same as adaptive_order_wcns3() except that the Z weights are used. First we define

\[\tau_3 = |\beta_1 - \beta_0|\]

Then the new \(\alpha_k\) are given by

\[\alpha_k = c_k\left(1 + \frac{\tau_3} {\beta_k + \epsilon_k}\right),\]

Parameters

• q (list of double) – The variable values at the cell centers.

• i (int) – The index into the reconstructed array

• j (int) – The index of the cell whose faces are being reconstructed in q

• recons (list of double) – The array of the reconstructed variable.

• keep_positive (bool) – If True then returns False if the reconstructed solution is not positive.

• eps (double) – The epsilon parameter to avoid division by zero.

• c0 (double) – The optimal linear weight \(c_{0}\). For 3rd order use \(1/4\). For 2nd order but increased robustness use \(1/2\).

• c1 (double) – The optimal linear weight \(c_{1}\). For 3rd order use \(3/4\). For 2nd order but increased robustness use \(1/2\).

• liu_indicators (bool) – If True use the oscillation indicators of [5]

Returns

(bool) True if the reconstruction was successful, otherwise False

ReconstructionPpao.adaptive_order_weno3_robust(q, i, j, recons, keep_positive, eps=1e-17, c1=1.0, c2=1000.0, c3=1.0, exponent=4, wenoz=False)

A robust WENO3 reconstruction using 5 points.

The individual polynomials stencils for the reconstruction are written as

\[\begin{align} u(\xi) = u_0 + u_\xi P_1(\xi) + u_{\xi\xi} P_2(\xi) \end{align}\]

The left-, central-, and right-biased stencils for the one-dimensional coefficients are:

\[\begin{split}\begin{align} u_{\xi}^{(L)}&=\frac{1}{2} u_{-2} - 2u_{-1} + \frac{3}{2} u_0 \\ u_{\xi\xi}^{(L)}&=\frac{u_{-2} - 2 u_{-1} + u_0}{2} \\ u_{\xi}^{(C)}&=\frac{1}{2}(u_1 - u_{-1}) \\ u_{\xi\xi}^{(C)}&=\frac{u_{-1} - 2 u_0 + u_1}{2} \\ u_{\xi}^{(R)}&=-\frac{3}{2}u_0 + 2 u_1 - \frac{1}{2} u_2 \\ u_{\xi\xi}^{(R)}&=\frac{u_{0} - 2 u_{1} + u_2}{2}. \end{align}\end{split}\]

The oscillation indicators are given by

\[\beta_{(i)} = \left(u_\xi^{(i)}\right)^2 + \frac{13}{3}\left(u_{\xi\xi}^{(i)}\right)^2,\]

where \(i\in\{L,C,R\}\). The nonlinear weights are:

\[\begin{split}\begin{align} \omega_k &= \frac{\alpha_k}{\sum_{l=0}^{2}\alpha_l} \\ \alpha_k &= \frac{\lambda_k}{(\beta_k + \epsilon_k)^p} \end{align}\end{split}\]

where \(p\) is usually chosen to be 4 or 8, and \(\lambda_0=1\), \(\lambda_1=10^5\), and \(\lambda_2=1\).

To obtain the WENOZ weights use \(p=1\) and with the new oscillation indicators

\[\beta_k^Z=\frac{\beta_k}{\beta_k + \tau_5 + \epsilon_k}\]

where

\[\tau_5 = |\beta_3 - \beta_1|.\]

Parameters

• q (list of double) – The variable values at the cell centers.

• i (int) – The index into the reconstructed array

• j (int) – The index of the cell whose faces are being reconstructed in q

• recons (list of double) – The array of the reconstructed variable.

• keep_positive (bool) – If True then returns False if the reconstructed solution is not positive.

• eps (double) – The epsilon parameter to avoid division by zero.

• c0 (double) – The linear weight \(\lambda_{0}\).

• c1 (double) – The linear weight \(\lambda_{1}\).

• c2 (double) – The linear weight \(\lambda_{2}\).

• exponent (double) – The exponent \(p\) in denominator of the \(\alpha_k\).

• wenoz (bool) – If True then use the WENOZ weights.

Returns

(bool) True if the reconstruction was successful, otherwise False