.. _wp_algorithms:
Wave-propagation algorithms
===========================
.. _wp_1d:
One space dimension
-------------------
Clawpack can be used to solve a system of equations of the form
.. math::
\kappa(x)q_t+f(q)_x = \psi(q,x,t),
where :math:`q=q(x,t)\in{\cal R}^m`. The standard case of a homogeneous
conservation law has :math:`\kappa\equiv 1` and :math:`\psi\equiv 0`,
.. math::
q_t+f(q)_x=0.
:label: cons1d
The flux function :math:`f(q)` can also depend explicitly on :math:`x` and
:math:`t` as well as on :math:`q`.
Hyperbolic systems that are not in conservation form, e.g.,
.. math::
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
:math:`q_{i-1}` and :math:`Q_i`, returns a set of :math:`M_w` waves
:math:`{\cal W}^p_{i-1/2}` and speeds
:math:`s^p_{i-1/2}`
satisfying
.. math::
\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
:math:`{\cal A}^-\Delta Q_{i-1/2}` and
a right-going fluctuation :math:`{\cal A}^+\Delta Q_{i-1/2}`. In
the standard conservative case :eq:`cons1d` these should satisfy
.. math::
{\cal A}^-\Delta Q_{i-1/2}+{\cal A}^+\Delta Q_{i-1/2} = f(Q_i)-f(Q_{i-1})
:label: asum
and the fluctuations then define a "flux-difference splitting".
.. math::
{\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 :math:`s^- = \min(s,0)` and :math:`s^+ = \max(s,0)`. In the
nonconservative case \eqn{claw_1dnoncon}, there is no "flux function"
:math:`f(q)`,
and the constraint :eq:`asum` need not be satisfied.
Godunov's method
----------------
Only the fluctuations are used for the first-order Godunov method, which is
implemented in the form
.. math::
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 :math:`\kappa \equiv 1`.
The Riemann solver must be supplied by the user in the form of a subroutine
`rp1`, as described in :ref:`user_riemann`.
Typically the Riemann solver first computes waves and speeds and then uses
these to compute :math:`{\cal A}^+Q_{i-1/2}` and :math:`{\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
.. math::
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
.. math::
\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 :math:`\tilde{\cal W}_{i-1/2}^p` represents a limited version of the wave
:math:`{\cal W}_{i-1/2}^p`, obtained by comparing :math:`{\cal W}_{i-1/2}^p` to
:math:`{\cal W}_{i-3/2}^p` if :math:`s^p>0` or to :math:`{\cal W}_{i+1/2}^p`
if :math:`s^p<0`.
.. _wp_fwave:
f-wave formulation
-------------------
For equations with spatially-varying flux functions, such as
.. math::
q_t+f(q,x)_x=0.
:label: cons1dvf
it is often convenient to use the f-wave formulation of the algorithm, as
proposed in [BaleLevMitRoss02]_. Rather than decomposing the jump
:math:`Q_i-Q_{i-1}` into waves, the idea of the f-wave approach is to
decompose the flux difference :math:`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 :ref:`riemann_fwave` for more details.
Capacity functions
------------------
When a capacity function :math:`\kappa(x)` is present, the Godunov method becomes
.. math::
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 :math:`q_t=\psi` over a time step.
A fractional step method is used to couple this with the homogeneous
solution, as described in :ref:`user_src`.
Boundary conditions
-------------------
Boundary conditions are imposed by setting values in ghost cells each time
step, as described in :ref:`bc`.
TODO: Continue description -- 2d and 3d, transverse solvers.