.. _riemann: Riemann solvers =============== The Riemann solver defines the hyperbolic equation that is being solved and does the bulk of the computational work -- it is called at every cell interface every time step and returns the information about waves and speeds that is needed to update the solution. The directory `$CLAW/riemann/src` contains Riemann solvers for many applications, including advection, acoustics, shallow water equations, Euler equations, traffic flow, Burgers' equation, etc. .. _rp1: One-dimensional Riemann solver ------------------------------ Understanding the 1-dimensional solver is a critical first step since in 2 or 3 dimensions the bulk of the work is done by a "normal solver" that solves a 1-dimensional Riemann problem normal to each cell edge. (These are then augmented by transverse solvers as described below.) See :ref:`wp_algorithms` and [LeVeque-FVMHP]_ for more details about how the algorithms in Clawpack use the Riemann solution. The calling sequence for the 1-dimensional Riemann solver subroutine `rp1` is:: subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq) ! Input arguments integer, intent(in) :: maxmx,meqn,mwaves,mbc,mx,maux double precision, intent(in), dimension(meqn, 1-mbc:maxmx+mbc) :: ql,qr double precision, intent(in), dimension(maux, 1-mbc:maxmx+mbc) :: auxl,auxr ! Output arguments double precision, intent(out) :: s(mwaves, 1-mbc:maxmx+mbc) double precision, intent(out) :: fwave(meqn, mwaves, 1-mbc:maxmx+mbc) double precision, intent(out), dimension(meqn, 1-mbc:maxmx+mbc) :: amdq,apdq Note that the integer parameters used in this calling sequence have different names than are now being used in `setrun.py`. The correspondence is: * `meqn = num_eqn`, the number of equations in the system, * `mwaves = num_waves`, the number of waves in each Riemann solution, * `mx = num_cells`, the number of grid cells, * `maux = num_aux`, the number of auxiliary variables, The input data consists of two arrays `ql` and `qr`. The value ql(:,i) is the value :math:`Q_i^L` at the left edge of the i’th cell, while qr(:,i) is the value :math:`Q_i^R` at the right edge of the i’th cell, as indicated in this figure: .. image :: images/qlqr.png :width: 10cm Normally `ql = qr` and both values agree with :math:`Q_i^n` , the cell average at time :math:`t_n`. More flexibility is allowed because in some applications, or in adapting clawpack to implement different algorithms, it is useful to allow different values at each edge. For example, we might want to define a piecewise linear function within the grid cell as illustrated in the figure and then solve the Riemann problems between these values. This approach to high-resolution methods is discussed in Chapter 6 of [LeVeque-FVMHP]_. Note that the Riemann problem at the interface :math:`x_{i−1/2}` between cells :math:`i − 1` and :math:`i` has data: * Left state: :math:`Q_{i-1}^R` = `qr(:,i-1)`, * Right state: :math:`Q_{i}^L =` `ql(:,i)`. This notation is rather confusing since normally we use :math:`q_\ell` to denote the left state and :math:`q_r` to denote the right state in specifying Riemann data. The Riemann solver `rp1` also has input parameters `auxl` and `auxr` that contain values of the auxiliary variables if these are being used (see :ref:`user_setaux`). Normally `auxl = auxr = aux` when the Riemann solver is called from the standard Clawpack routines. The routine `rp1` must solve the Riemann problem for each value of `i`, and return the following (for `i=1-mbc` to `mx+mbc`): * `amdq(1:meqn,i)` = the vector :math:`{\cal A}^-\Delta Q_{i-1/2}`, * `apdq(1:meqn,i)` = the vector :math:`{\cal A}^+\Delta Q_{i-1/2}`, * `wave(1:meqn,i,p)` = the vector :math:`{\cal W}^p_{i-1/2}`, * `s(i,p)` = the wave speed :math:`s^p_{i-1/2}` for each wave. This uses the notation of :ref:`wp_algorithms` and [LeVeque-FVMHP]_. .. _riemann_fwave: f-wave Riemann solvers ---------------------- As described in :ref:`wp_fwave`, for spatially-varying flux functions it is often best to use the f-wave formulation of the wave-propagation algorithms. This can be implemented in Clawpack by providing a suitable Riemann solver that returns f-waves instead of ordinary waves, obtained by decomposing 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. The Riemann solver has the same form and calling sequence as described above, the only difference is that the output argument `wave` should return an array of f-waves. while `amdq` and `apdq` have the same meaning as before. In order to indicate that the Riemann solver returns f-waves, the parameter `runclaw.use_fwaves` in `setrun` should be set to `True`, see :ref:`setrun`. TODO: Continue description -- 2d and 3d, transverse solvers.