Advection routine

Following the notation in [BK19], the quantities of interest that we are advecting are

\[\Psi = (\chi, \chi \vec{u}, \chi w, \chi^\prime)\]

under the advective fluxes \((P \vec{v})\). Say we are doing the advection for the full time-step update, i.e. the quantities \((P \vec{v})^{n+1/2}\) and \(\Psi^n\) are available, where \(n\) indexes the time-level.

Consider (21) in [BK19] for the cells \(i-1\), \(i\) and \(i+1\),

(1)\[\mathcal{A}^{\Delta t / 2}_{x} \mathcal{U}_i = \mathcal{U}_i - \frac{\Delta t}{2 \Delta x} \left[ (Pu)^{n+1/2}_{i+1/2} \Psi_{i+1/2} - (Pu)^{n+1/2}_{i-1/2} \Psi_{i-1/2} \right],\]

where, corresponding to the management.variable.Vars data container,

\[\mathcal{U}_i := \left\{ \rho, \rho \vec{u}, \rho w, P, P \chi^\prime \right\}_i.\]

We see that (1) is a half time-step update,

\[\mathcal{A}^{\Delta t / 2}_{x} \mathcal{U}_i = \mathcal{U}_i - \frac{1}{\Delta x} \left[ \int_{t}^{t+\frac{\Delta t}{2}} \mathcal{F}(\mathcal{U}_i, \mathcal{U_{i+1}}) - \mathcal{F}(\mathcal{U}_{i-1},\mathcal{U}_i) ~ dt \right],\]

where \(\mathcal{F}(\cdot,\cdot)\) is the flux across the cell interface for the respective quantity \(\mathcal{U}\). This is from a second-order MUSCL scheme. For a general advection equation,

\[\mathcal{U}_t + \mathcal{F}_x(\mathcal{U}) = 0,\]

and discretising the above yields

\[\frac{\mathcal{U}^{n+1} - \mathcal{U}^n}{\Delta t} = - \frac{1}{\Delta x} \left[ \mathcal{F}(\mathcal{U}_{i+1/2}) - \mathcal{F}(\mathcal{U}_{i-1/2})\right],\]

where the right-hand side is a central difference scheme evaluated at the cell edges. Rewriting this,

(2)\[\mathcal{U}^{n+1} = \mathcal{U}^n - \frac{\Delta t}{\Delta x} \left[ \mathcal{F}(\mathcal{U}_{i+1/2}) - \mathcal{F}(\mathcal{U})_{i-1/2} \right],\]

we identify the first term in the square bracket on the right with \((P u)^{n+1/2}_{i+1/2} \Psi_{i+1/2}\) and the second with \((P u)^{n+1/2}_{i-1/2} \Psi_{i-1/2}\). Recall that the advective fluxes \((P u)^{n+1/2}\) are available at the onset of the full time-stepping. A cell-to-face averaging is done to obtain these values at \(i+1/2\) and \(i-1/2\). The missing pieces are then \(\Psi^{n+1/2}_{i-1/2}\) and \(\Psi^{n+1/2}_{i+1/2}\).

Todo

Link or insert numerical_fluxes.recompute_advective_fluxes here.

To obtain second-order accuracy in space for the solution, we use a slope limiter to approximate the cell-averages instead of the piecewise constants cell-averages in a first-order Godunov’s method. This sets up a Riemann problem that is solved by a HLLE solver.

Reconstruction

Fig. 2 Recovering the limited slope (in red) in the \(i\)-th cell.

Let’s get the missing pieces, say \(\Psi^{n+1/2}_{i+1/2}\). First, apply the slope limiter. Fig. 2 shows a simple case where the slope in the \(i\)-th cell is the average of the slopes constructed from the cell averages (dotted-lines) of the adjacent cells. This is done in physics.gas_dynamics.recovery.slopes() while physics.gas_dynamics.recovery.limiters() is a switch for the limiter choice specified in the initial coonditions. The limited slope of the \(i\)-th cell is denoted \(S_i\). Now,

(3)\[\Psi(t^{n+1/2}, x_{i+1/2}) = \Psi(t^n, x_{i+1/2} - \frac{\Delta t}{2} u)\]

by tracing the characteristic of \(\Psi\) backwards in time. Fig. 3 illustrates this.

Fig. 3 Backward tracing of the characteristic from \((t^{n+1/2}, x_{i+1/2})\).

\(u\) is the velocity at the cell face \(i+1/2\), with

\[u_{i+1/2} = \frac{(P u)^{n+1/2}_{i+1/2}}{(P_i + P_{i+1})/2},\]

as \(u\) is a constant along the characteristic at \(x_{i+1/2}\).

As we now have the limited slope in the cell, the value for \(\Psi\) anywhere in the cell can be obtained by linear interpolation,

(4)\[\Psi_i(t^n, x) = \Psi_i^n + (x - x^n_i) S_i^n.\]

Plugging (4) into (3),

(5)\[\begin{split}\Psi(t^{n+1/2}, x_{i+1/2}) &= \Psi^n_i + \left( x_{i+1/2}^n - \frac{\Delta t}{2} u - x^n_i \right) S^n_i \\ &= \Psi^n_i + \frac{\Delta x}{2} \left( 1 - \frac{\Delta t}{\Delta x} u \right) S^n_i \\ &=: \Psi^-_{i+1/2},\end{split}\]

where \(x_{i+1/2} - x_{i} = \Delta x /2\) is used. \(\Psi^-_{i+1/2}\) is the solution of \(\Psi\) from the left for the cell interface at \(i+1/2\). See Fig. 4 for more details.

The Riemann problem

Fig. 4 Discontinuity at the cell interface \(i+1/2\) arising from the recovery of \(\Psi_{i+1/2}^{n+1/2}\).

(5) supposes that the characteristic for \(\Psi^{n+1/2}_{i+1/2}\) is traced from the \(i\)-th cell. Another possibility is from the \((i+1)\)-th cell,

(6)\[\begin{split}\Psi (t^{n+1/2}, x_{i+1/2}) &= \Psi^n_{i+1} - \frac{\Delta x}{2} \left( 1 + \frac{\Delta t}{\Delta x} u \right) S^n_{i+1} \\ &=: \Psi^{+}_{i+1/2}.\end{split}\]

(2), (5) and (6) constitute the Riemann problem. (5) and (6) are computed in physics.gas_dynamics.recovery.recovery().

HLLE Solver

The Riemann problem is solved by selecting the choice of (5) or (6) based on the direction of the flux,

(7)\[\Psi_{i+1/2} = \sigma \Psi^{-}_{i+1/2} + (1-\sigma) \Psi^+_{i+1/2},\]

where

(8)\[\sigma = \text{sign}\left[ (Pu)^{n+1/2}_{i+1/2} \right].\]

(7) and (8) are computed in physics.gas_dynamics.numerical_flux.hll_solver(). With (7), (2) is updated, yielding the solution to (1).

Strang splitting

(1) demonstrates the solution of a substep obtained from the Strang splitting in the advection scheme. Writing the advection in the full time-stepping of \(\mathcal{U}^*\) as a general operator,

\[\mathcal{U}^{**} = \mathcal{A}_{\text{full}}^{\Delta t} ~ \mathcal{U}^*,\]

where \(\mathcal{A}_{\text{full}}^{\Delta t}\) is the advection operator. Splitting the operator dimension-wise,

\[\mathcal{U}^{**} = \mathcal{A}_{x}^{\Delta t/2} \mathcal{A}_{y}^{\Delta t/2} \mathcal{A}_{z}^{\Delta t/2} \mathcal{A}_{z}^{\Delta t/2} \mathcal{A}_{y}^{\Delta t/2} \mathcal{A}_{x}^{\Delta t/2} \mathcal{U}^*,\]

recalling that the Strang splitting is a second-order operator splitting method. The Strang-splitting is computed in physics.gas_dynamics.explicit.advect() which calls physics.gas_dynamics.explicit.explicit_step_and_flux() for each substep.

More details on the methods discussed here can be found in [LeVeque92].

References

[LeVeque92]

LeVeque, R.J., 1992: Numerical methods for conservation laws. Basel: Birkhäuser.