Hi everyone, I believe the routines in GRHydro_SlopeLimiter.F90 are all either mislabelled or simply the wrong formulae. As a quick reminder, for the upwind slope ("dupw") and downwind slope ("dloc"), we expect to get the following results:
If dupw * dloc <= 0, delta=0 [slope limiters go flat near extrema, by design]
else
if (minmod)
delta = sign(dupw) * min( abs(dupw), abs(dloc))
else if (MC)
delta = sign(dupw) * min( 2*abs(dupw), 2*abs(dloc), abs(dupw+dloc)/2 )
else if (superbee)
delta = sign(dupw) * max[ min(2*abs(dupw), dloc) , min(2*abs(dloc),dupw) ]
On the rather nice wikipedia page, http://en.wikipedia.org/wiki/Flux_limiter
they take variables r = dupw/dloc and phi = delta / dloc, but the results are the same.
Attached are a spreadsheet where I tried to reproduce what the current routine actually calculates for various inputs, along with the recommended code fix for GRHydro_TVDReconstruct.F90.
In summmary: the current "MC2" and "minmod" options work like they should. "MC1" is actually minmod with a floor for very small values. minmod2, minmod3, and superbee are all incorrect expressions.
Cheers, Josh