Separating entropy fixes from Riemann solvers

[David Ketcheson]

Gorune Ohannessian, who is developing CUDAClaw, brought to me a suggestion today regarding the implementation of entropy fixes.  The suggestion has both aesthetic and performance consequences.

The computation of fluctuations (amdq, apdq) can be done in a generic way independent of details of the Riemann solver.  Thus it makes sense to do this outside of the Riemann solver, and not write duplicate code in every Riemann solver for computing fluctuations.  In CUDAClaw, this also has the performance advantage of reducing the amount of shared memory that needs to be used, because the Riemann solver only needs to return waves and speeds.

An exception to the above is that an entropy fix modifies the fluctuations directly, and is performed inside the Riemann solver.  The entropy fix usually does not modify the waves themselves.  It seems that the entropy fix needs to go inside the Riemann solver, and so the fluctuations also must be computed there.

However, most Clawpack solvers use a Harten-Hyman entropy fix (pp. 324-325 of Randy's FV book).  As far as I can tell, that entropy fix only requires knowledge of the waves and speeds, and can be performed in a totally black-box way for any system of equations.  Thus, we could take the entropy fix computation out of the Riemann solver and put it in as a separate subroutine called afterward.  This would also have the advantage of automatically providing an entropy fix for any system.  The user could choose whether to turn on the entropy fix, in order to avoid performance loss in systems where it is not needed.

I think this proposal would remove code that is duplicated across many Riemann solvers (i.e., the computation of fluctuations), would improve performance on GPUs, and would enable anyone to have a simple entropy fix for new systems.  Does this seem like a good idea to you?

-David

[Randall J LeVeque]

This sounds attractive, but unfortunately the Harten-Hyman fix
generally needs to know more about the system than just the waves and
speeds that come out of something like a Roe solver. The point of the
fix is to effectively spread a jump propagating at a single speed s
into corrections that affect the cells on both sides of the interface
when the wave should really be a transonic rarefaction. But to detect
that and use the H-H formulas, we need an estimate of the wave speed
for this $k$th family on each side of the rarefaction wave, i.e. an
estimate of the characteristic speed in each intermediate state
bounding this wave (the $\lambda^k_{l,r}$ in formula (15.49)). The
intermediate states are easily compute from the waves, but then
formulas are needed for the characteristic speed in these states.
These are easily computed for systems like Euler or shallow water but
of course are different than the speed coming out of the Roe solver.

Perhaps there's still a nice way to turn it into more of a black box
if it's assumed there is also another function available that returns
the eigenvalues of the Jacobian f'(q) for any intermediate state q,
for example.

- Randy

> --
> You received this message because you are subscribed to the Google Groups
> "claw-dev" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to claw-dev+u...@googlegroups.com.
> To post to this group, send email to claw...@googlegroups.com.
> Visit this group at http://groups.google.com/group/claw-dev?hl=en.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[David Ketcheson]

On Fri, Apr 26, 2013 at 9:52 PM, Randall J LeVeque <r...@uw.edu> wrote:

This sounds attractive, but unfortunately the Harten-Hyman fix
generally needs to know more about the system than just the waves and
speeds that come out of something like a Roe solver.  The point of the
fix is to effectively spread a jump propagating at a single speed s
into corrections that affect the cells on both sides of the interface
when the wave should really be a transonic rarefaction.  But to detect
that and use the H-H formulas, we need an estimate of the wave speed
for this $k$th family on each side of the rarefaction wave, i.e. an
estimate of the characteristic speed in each intermediate state
bounding this wave (the $\lambda^k_{l,r}$ in formula (15.49)).  The
intermediate states are easily compute from the waves, but then
formulas are needed for the characteristic speed in these states.
These are easily computed for systems like Euler or shallow water but
of course are different than the speed coming out of the Roe solver.

Perhaps there's still a nice way to turn it into more of a black box
if it's assumed there is also another function available that returns
the eigenvalues of the Jacobian f'(q) for any intermediate state q,
for example.

You're right, and in fact Gorune realized this and proposed to me what you've said in your second paragraph (require a routine that computes eigenvalues of the flux Jacobian).  I muddled it when I transmitted his suggestion here.

I suppose it wouldn't be difficult to at least support this as an option?  It would involve changes similar to those required to support f-wave solvers.

-David

[Randall J LeVeque]

Such a routine could have other uses, e.g. setting the initial dt in
an AMR calculation. In GeoClaw we added something like this.

Of course the exact eigenvalues might not be available for some
systems but as long as it's optional for cases where you want to use
it, this seems reasonable.

- Randy

[David Ketcheson]

Okay, I think Gorune will go ahead and implement this in CUDAClaw as a first step.

[Donna]

I did something like this once - separated out the HH entropy fix from the Riemann problem.   As others have noted, it requires knowledge of the eigenvalues of the system.   

I've attached the Riemann solver and entropy fix that I once wrote for a detonation problem (Euler + one more variable), and have included in comments where one can call something like "calculate_speeds(q)".    The updates to apdq and amdq are also fairly general and based on modified speeds that come out of the entropy fix.  

(I tried posting this before, but never saw my post show up;  my apologizes if this now shows up multiple times)

[Calhoun Donna]

I once did something like this for Euler - that is, split off the HH entropy fix  from rpn2.  It wasn't black-box, because, as Randy said, these routines depend on knowing the eigenvalues of the system.  But I've indicated in comments where one might substitute in something like "calculate_speed(q)" in place of where the eigenvalues are explicitly computed.   I also incorporated this fix into the update of apdq and amdq in a way that I think is fairly general. 

I've attached a code that I once wrote to implement this idea for a detonation problem (Euler + an extra variable that tracks the detonation front). 

Please send any comments or questions, 

Donna