Hello Ian,
My memory is that the idea is related to the smoothness of the metric/gauge terms: these should be C^2, so are sufficiently differentiable to be "best" approximated (at the second order that Whisky/GRHydro worked to) by the midpoint rule.
Ok, that is partially what we are doing (ie we do not reconstruct lapse explicitly (it is not stored as a volume average anyway) and compute it as just the next neighbour average. The confusion for us was that essentially GRHydro / Whisky do something like this:
num_flux_(i-1/2) = FOO_(i-1/2)(prims, gxx, beta/lapse) num_flux_(i+1/2) = FOO_(i+1/2)(prims, gxx, beta/lapse)
where num_flux is the result of the Rieman solver and reconstruction etc.
*Then* when it computes the actual right hand side terms for MoL it uses:
(1)
RHS_i = 0.5*(lapse_(i) + lapse_(i+1)) num_flux(i+1/2) - 0.5*(lapse_(i) + lapse_(i-1)) num_flux(i-1/2)
instead of the more natural (given the papers and also GRHydro docs):
(2)
RHS_i = num_flux(i+1/2) - num_flux(i-1/2)
(the two num_flux in method (1) and method (2) would differ in that one absorbed the extra lapse term while the other did not). There are always some metric terms in the fluxes (eg shift is always present and so of course is a component of the 3-metric itself).
However, if included within the fluxes then they would implicitly go through the reconstruction/Riemann solve stage, so would - in places - be limited to lower accuracy.
I would not do this now - keeping the fully densitized terms together is cleaner and a better way to extend to high order.
Ok, so from this I understand that we are ok doing things "different" from GRHydro and use method (2), at least until we have evidence for one method being clearly superior to the other. At least there is no clear "you must never do this" no go claim.
Yours, Roland