Wave-propagation algorithms

One space dimension

Clawpack can be used to solve a system of equations of the form

\[\kappa(x)q_t+f(q)_x = \psi(q,x,t),\]

where \(q=q(x,t)\in{\cal R}^m\). The standard case of a homogeneous conservation law has \(\kappa\equiv 1\) and \(\psi\equiv 0\),

(1)\[q_t+f(q)_x=0.\]

The flux function \(f(q)\) can also depend explicitly on \(x\) and \(t\) as well as on \(q\). Hyperbolic systems that are not in conservation form, e.g.,

\[q_t+A(x,t)q_x=0,\]

can also be solved. See [LeVeque-FVMHP] for more details about the wave-propagation algorithms that are briefly described here.

The basic requirement on the homogeneous system is that it be hyperbolic in the sense that a Riemann solver can be specified that, for any two states \(q_{i-1}\) and \(Q_i\), returns a set of \(M_w\) waves \({\cal W}^p_{i-1/2}\) and speeds \(s^p_{i-1/2}\) satisfying

\[\sum_{p=1}^{M_w} {\cal W}^p_{i-1/2} = Q_i-Q_{i-1} \equiv \Delta Q_{i-1/2}.\]

The Riemann solver must also return a left-going fluctuation \({\cal A}^-\Delta Q_{i-1/2}\) and a right-going fluctuation \({\cal A}^+\Delta Q_{i-1/2}\). In the standard conservative case (1) these should satisfy

(2)\[{\cal A}^-\Delta Q_{i-1/2}+{\cal A}^+\Delta Q_{i-1/2} = f(Q_i)-f(Q_{i-1})\]

and the fluctuations then define a “flux-difference splitting”.

\[{\cal A}^-\Delta Q_{i-1/2} = \sum_p (s^p_{i-1/2})^- {\cal W}^p_{i-1/2}, \qquad {\cal A}^+\Delta Q_{i-1/2} = \sum_p (s^p_{i-1/2})^+ {\cal W}^p_{i-1/2},\]

where \(s^- = \min(s,0)\) and \(s^+ = \max(s,0)\). In the nonconservative case eqn{claw_1dnoncon}, there is no “flux function” \(f(q)\), and the constraint (2) need not be satisfied.

Godunov’s method

Only the fluctuations are used for the first-order Godunov method, which is implemented in the form

\[Q_i^{n+1} = Q_i^n - \frac{\Delta t}{\Delta x}\left[{\cal A}^+\Delta Q_{i-1/2} + {\cal A}^-\Delta Q_{i+1/2}\right],\]

assuming \(\kappa \equiv 1\).

The Riemann solver must be supplied by the user in the form of a subroutine rp1, as described in Specifying the Riemann solver.

Typically the Riemann solver first computes waves and speeds and then uses these to compute \({\cal A}^+Q_{i-1/2}\) and \({\cal A}^-Q_{i-1/2}\) internally in the Riemann solver.

High-resolution methods

The waves and speeds must also be returned by the Riemann solver in order to use the high-resolution methods described in [LeVeque-FVMHP], which reduce to a second-order accurate Lax-Wendroff type method when limiters are not used. By introducing wave limiters, non-physical oscillations near discontinuities or steep gradients in the solution can be avoided. The limiters are based on the theory of TVD methods for nonlinear scalar equations and extended in a natural way to systems of equations.

These methods take the form

\[Q_i^{n+1} = Q_i^n - \frac{\Delta t}{\Delta x}\left[{\cal A}^+Q_{i-1/2} + {\cal A}^-Q_{i+1/2}\right] - \frac{\Delta t}{\Delta x}(\tilde F_{i+1/2} - \tilde F_{i-1/2})\]

where

\[\tilde F_{i-1/2} = \frac 1 2 \sum_{p=1}^{M_w} |s^p_{i-1/2}| \left( 1-\frac{\Delta t}{\Delta x} |s^p_{i-1/2}|\right) \tilde{\cal W}_{i-1/2}^p.\]

Here \(\tilde{\cal W}_{i-1/2}^p\) represents a limited version of the wave \({\cal W}_{i-1/2}^p\), obtained by comparing \({\cal W}_{i-1/2}^p\) to \({\cal W}_{i-3/2}^p\) if \(s^p>0\) or to \({\cal W}_{i+1/2}^p\) if \(s^p<0\).

f-wave formulation

For equations with spatially-varying flux functions, such as

(3)\[q_t+f(q,x)_x=0.\]

it is often convenient to use the f-wave formulation of the algorithm, as proposed in [BaleLevMitRoss02]. Rather than decomposing the jump \(Q_i-Q_{i-1}\) into waves, the idea of the f-wave approach is to decompose the flux difference \(f(Q_i,x_i) - f(Q_{i-1},x_{i-1})\) into f-waves using appropriate eigenvectors of the Jacobian matrices to either side of the interface. See f-wave Riemann solvers for more details.

Capacity functions

When a capacity function \(\kappa(x)\) is present, the Godunov method becomes

\[Q_i^{n+1} = Q_i^n - \frac{\Delta t}{\kappa_i \Delta x} \left[{\cal A}^+Q_{i-1/2} + {\cal A}^-Q_{i+1/2}\right],\]

See [LeVeque-FVMHP] for discussion of this algorithm and its extension to the high-resolution method. Capacity functions are useful in particular for solving equations on a mapped grid.

Source terms

If the equation has a source term, a routine src1 must also be supplied that solves the source term equation \(q_t=\psi\) over a time step. A fractional step method is used to couple this with the homogeneous solution, as described in Using src for source terms.

Boundary conditions

Boundary conditions are imposed by setting values in ghost cells each time step, as described in Boundary conditions.

TODO: Continue description – 2d and 3d, transverse solvers.