.. _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.