CEP1

Status	Active
Author	Kyle Mandli <kyle.mandli@gmail.com>
Created	Nov. 4, 2014
Updated	March 19, 2015
Discussion	Riemann #87
Implementation	Riemann #91

Motivation

Currently, a Riemann solver routine must include a loop over a 1D slice of the grid. This was done for efficiency reasons, to avoid the overhead of having the Riemann solver function call inside the lowest-level loop. It seems likely that modern compilers can avoid this overhead, allowing us to move the Riemann solver function call inside the loop.

Thus, a Riemann solver would take as input just the two states involved in a specific Riemann problem. It is more natural to write such a Riemann solver, since it corresponds directly to the mathematical concept of a Riemann problem.

In the long-term, we envision a transition to using pointwise Riemann solvers exclusively. However, as a first step we will develop pointwise solvers that can be drop-in replacements for existing vectorized Riemann solvers. This will be facilitated by a new layer of code that implements the loop over a 1D slice.

Specification

The 1D pointwise solver definition:

subroutine rp1_ptwise(num_eqn, num_aux, num_waves, q_l, q_r, aux_l, aux_r, wave, s, amdq, apdq)

    implicit none

    integer, intent(in) :: num_eqn, num_aux, num_waves
    real(kind=8), intent(in) :: q_l(num_eqn), q_r(num_eqn)
    real(kind=8), intent(in) :: aux_l(num_aux), aux_r(num_aux)
    real(kind=8), intent(in out) :: wave(num_eqn, num_waves)
    real(kind=8), intent(in out) :: s(num_waves)
    real(kind=8), intent(in out) :: apdq(num_eqn), amdq(num_eqn)

end subroutine rp1_ptwise

The 2D normal pointwise solver definition is the same, except for the function name and the addition of an ixy argument:

subroutine rpn2_ptwise(ixy,num_eqn, num_aux, num_waves, q_l, q_r, aux_l, aux_r, wave, s, amdq, apdq)

    implicit none

    integer, intent(in) :: num_eqn, num_aux, num_waves, ixy
    real(kind=8), intent(in) :: q_l(num_eqn), q_r(num_eqn)
    real(kind=8), intent(in) :: aux_l(num_aux), aux_r(num_aux)
    real(kind=8), intent(in out) :: wave(num_eqn, num_waves)
    real(kind=8), intent(in out) :: s(num_waves)
    real(kind=8), intent(in out) :: apdq(num_eqn), amdq(num_eqn)

end subroutine rpn2_ptwise

The 2D transverse solver definition is similar. Note that we pass in three left and right states and aux values, corresponding to the values in adjacent cells:

subroutine rpt2_ptwise(ixy, imp, ,num_eqn, num_aux, num_waves, q_l, q_r, aux_l, aux_r, asdq, bmasdq, bpasdq)

    implicit none

    integer, intent(in) :: num_eqn, num_aux, num_waves, ixy, imp
    real(kind=8), intent(in) :: q_l(num_eqn,3), q_r(num_eqn,3)
    real(kind=8), intent(in) :: aux_l(num_aux,3), aux_r(num_aux,3)
    real(kind=8), intent(in out) :: asdq(num_eqn)
    real(kind=8), intent(in out) :: bmasdq(num_eqn), bpasdq(num_eqn)

end subroutine rpt2_ptwise

The extension to 3D is left as an exercise for the reader. ;-)

Reference implementation

There is now an example implementation working and tested in https://github.com/mandli/riemann/tree/ptwise which should work with the Fortran and PyClaw packages. Some changes to PyClaw were helpful and can be found at https://github.com/mandli/pyclaw/tree/ptwise-support. There's also an example for classic usage at http://github.com/mandli/classic/tree/ptwise-example.

Backwards compatibility issues

This will be fully backward-compatible.