GR-AMR incompatibility about athena

[tomidakn]

Not really, I'll take a loot at it. Now I'm back in Princeton and should have more time.

[tomidakn]

Hi @c-white ,

I finally got some time to work on AMR, and I think I can take care of this issue alog with it. Can you tell me how I should restrict/prolongate the variables for GR?

[c-white]

Since the primitives are just used for an initial guess and a fallback plan, I don't think we need to do anything special. Just do the same thing to them as with the conserved values. We should do it for both w and w1 for now to be safe, though in the future just handling w might be enough.

[felker]

I don't expect this will change anything regarding the fix, but I reconfigured the time_integrator.cpp so that it always writes out to prim_out= w1 and then the AthenaArray data pointer is swapped with w:

  if (step <= nsub_steps) {
    // At beginning of this task, phydro->w contains previous substep W(U) output
    // and phydro->w1 is used as a register to store the current substep output.
    // For the second order integrators VL2 and RK2, the prim_old initial guess for the
    // Newton-Raphson solver in GR EOS uses the following abscissae:
    // step=1: W at t^n and
    // step=2: W at t^{n+1/2} (VL2) or t^{n+1} (RK2)
    pmb->peos->ConservedToPrimitive(phydro->u, phydro->w, pfield->b,
                                    phydro->w1, pfield->bcc, pmb->pcoord,
                                    il, iu, jl, ju, kl, ku);
    // swap AthenaArray data pointers so that w now contains the updated w_out
    phydro->w.SwapAthenaArray(phydro->w1);
}
...

[tomidakn]

I think I just need to copy u and w.

Chris, can you test it if I modify the code?

[c-white]

Yes, I'll keep a close eye on what this does to GR.

[felker]

Any update on this, @tomidakn ?

[felker]

I cant remember, was this ever addressed @tomidakn @c-white ? Should I close this regardless, or is it still of interest?

[c-white]

I just checked and new and just-moved blocks still don't have old primitives, and these values are read by the variable inversion function. However, variable inversion has had so many branches added over the years to make it more robust, these values have little to no influence on the outcome.

I guess the problem is solved, though the real test will come when someone is brave enough to do a large GR+AMR production run.

In the long run, there are better options than falling back on old values. One standard in the community is to use a passively advected entropy. That is, recalculate entropy at the beginning of the step, advect it, and, if inversion fails, replace the energy equation with the entropy equation. It would be nice to turn this on with a switch, though it would require a specially-treated passive scalar to be used whether or not the user has enrolled others.

[felker]

Closing until we pick this up again:

• Test a large GR+AMR production run and evaluate stability of variable inversion fallback
• Consider alternative approach of advecting entropy via special passive scalar

@c-white so the idea would be to only allocate this special passive scalar during the timestep in which that switch is flipped? Do you have to calculate and advect the entropy for all cells, just the cell where the inversion failed, or a local neighborhood? What's a good reference for this?

How would the control flow be handled for restarting the entire timestep? Seems like a big change to the task list setup.
Might be a good time to implement --enable-fofc from Athena 4.2 within the scope of that change, although that seems much simpler, since it only affects the cell(s) that produce negative density or pressure in a VL2 timestep.

[c-white]

@felker I should clarify I was envisioning an input parameter that turned on entropy backups, in which case the scalar would exist everywhere for all times. It is simply reset to s ~ p/rho^gamma or whatever every step, and only read by the code in cells where variable inversion fails, requiring no change to the task list. It's just that whenever {D, E, M^i} -> {rho, p, v^i} fails, the eos will try {D, S, M^i} -> {rho, p, v^i} instead. Adiabatic evolution is of course wrong, but arguably less wrong than no evolution or resorting to floors.

All N versions of HARM, as well as probably BHAC, ECHO, Cosmos++, etc. probably do this, though it's seen better in the codes than the code papers.

FOFC is something I would very much like to see too, especially since it's not even wrong, just formally less accurate.

[felker]

Ok that is clear and seems more-or-less straightforward to implement. Although it seems inefficient to preemptively allocate and evolve this extra passive scalar, I suppose the user would likely know whether or not to suspect that their GR problem is susceptible to such failures. Still, it would be unfortunate to start a very expensive calculation without the entropy backup switch active, then find that it is necessary. The user could take a .rst file and restart it mid-simulation with the switch flipped on?

[c-white]

@felker Sure, this could be turned on mid-simulation; it's just a matter of having the code allocate more arrays.

Actually it's made easier by the fact that one would never dump this entropy, either to .athdf or to .rst files, since it can always be calculated from the primitives. So it's not like it would be missing.

Even if GR always turned it on, the marginal cost of 1 passive scalar must be pretty small compared to everything else GR does. But I agree it would be better to make it optional.