hybrid lambda term

[asnodin]

This implements an additional -\overline{\lambda}\rho g term in the momentum equation that should lead to better energy conservation when hydrostatic equilibrium is enforced (the form of the term(s) may be slightly different in the spherical case, that is not considered here). \overline{\lambda} is stored as BaseState<amrex::Real> lambdabar and included in MakeVelForce() by subtracting \rho\overline{\lambda} from rhopert. lambdabar is updated only at the end of the timestep in the nodal projection when Pi, from which it can be obtained, is also updated. \lambda (i.e. the non-horizontally averaged quantity) also appears in a correction to the temperature that is not included here (that will be done with Pi instead). (The thermodynamically consistent temperature correction depends on the term implemented here, but can in principle be used without it, with unknown accuracy.) The term is enabled using the use_lambdabar_term switch. Although this has been implemented and tested with planar geometry in mind, other than perhaps the special case phi degeneracy problem discussed below, I don't see any obvious reason why this should not (almost) already work with the spherical case.

lambdabar is obtained by averaging the cell-centered Pi and dividing by the time-centered p0 * gamma1bar:

auto pg = 0.25 * (p0_new + p0_old) * (gamma1bar_new + gamma1bar_old);
Average(snew, lambdabar, Pi);
lambdabar /= pg;

However, the update of Pi requires careful consideration.

In the case where the upper and lower boundaries are impenetrable, as with a closed box, the phiobtained from the nodal projection is only determined up to a constant (with the existing velocity constraint). It is important to determine this constant because it affects the value of Pi and hence lambdabar (and also the temperature correction). Whenever the new zero_average_lambda switch is set, the code adds a constant onto phi such that the volume average of \lambda is zero (an energy conservation constraint). The present implementation assumes an equidistantly-spaced radial grid and simply sums up the horizontally-averaged values -- I guess that could easily be generalized if required.

A related issue is that the existing velocity constraint is not really the right constraint when this term is included. An additional term involving \overline{\phi} should be present, but that would involve some work (perhaps this could be implemented by adding a horizontal average operator in the linear solvers?) Instead, three options have been added to work with the existing constraint via different values of lambda_update_method:

• lambda_update_method=0 => just use the standard constraint and accept that the \overline{\phi} quantity is in some sense lagged. (This is essentially a "do nothing" option and involves no additional code.)
• lambda_update_method=1 => similar to 0, but subtract/add the lambdabar term in the same fashion as with the gpi term before/after the projection
• lambdabar_update_method=N >=2. => iterate the nodal projection, updating only the lambdabar term in unew that is used to construct Vproj each time.

It is so far not clear from tests we have done if any of these options is advantageous over others (because the effect of including the term overall has been quite small). The iterative version converges in tests we have performed in a closed box, but it seems that very little, if anything, is gained by performing more than one additional iteration, i.e. N=2. (This may change after further study and perhaps one then might set a default tolerance for convergence, but it may be simply that one never needs to iterate in most applications.) Rather than changing the code in MaestroNodalProj.cpp significantly, the iteration is presently implemented from within AdvanceTimeStepAverage().