ensure atmosphere is properly set after ID

Create issue
Issue #1071 open
Frank Löffler created an issue

At the moment GRHydro doesn't ensure properly that the atmosphere is set after initial data. This is no problem as long as the ID thorn takes care of this, the ID doesn't contain atmosphere or the EOS is sufficiently robust to 'recover'. However, this is not always the case. Thus, I propose to set the atmosphere mask in PostInitial, initializing it according to whatever was setup there by some other thorn.

The patch to do that is attached. It is only a change to schedule.ccl as the corresponding function already exists.


Comments (21)

  1. Roland Haas
    • removed comment

    I like the idea of actually marking atmosphere points by a mask and not not only by their density since it will actually freeze their evolution.

    I would like to point out though that this will change current behaviour even for ID thorns that know about atmosphere. Right now the atmosphere mask is really mostly a "do-not-evolve" mask (or set-to-atmosphere) since it is set in GRHydroUpdateAtmosphereMask based on densrhs and taurhs. Also currently the atmosphere mask is cleared after each regrid so if we want to adopt the new logic in the patch of marking all points where rho<rho_min in the atmosphere mask (which would be what the name implies but not what it currently does) then we should also schedule GRHydro_InitAtmosMask in PostRegrid (and all other places where currently SetupAtmosphere is scheduled). Also we schedule AtmosphereReset in MoL_PostStep which clears the atmosphere mask (since wk_atmosphere .eq. 0 is the default) in normal runs.

    With the patch one should be a careful I believe as to what happens in the cells that border the atmosphere.

  2. Frank Löffler reporter
    • removed comment

    Yes, and I planned to propose to schedule GRHydro_InitAtmosMask as well, after testing it a bit more. This would be a more drastic change than only calling GRHydro_InitAtmosMask after Initial. Calling GRHydro_InitAtmosMask during evolution would most likely (I didn't check lately) freeze the surface of the star, which isn't what we want most of the time.

    Not a lot would change with the currently proposed change, although I agree that some data might be a bit different, especially the output for the first timestep. However, that is the whole point of the proposed change. Initial data thorns should not care about an atmosphere. The atmosphere is a 'problem' of the evolution, not of the initial data. Thus, e.g. the TOVSolver doesn't care about it (anymore), does not look at the different parameters in GRHydro and sets rho=press=eps=0 outside of the star.

  3. Roland Haas
    • removed comment

    ok. In that case I would think that rather than setting atmosphere_mask (which is what InitAtmosMask does) it would be better to reset the primitives and conservatives to atmosphere and not set the mask (since the mask will be cleared in the next UpdateMask anyway). I believe there is a routine in GRHydro that does this already: GRHydro_InitialAtmosphereReset which currently runs in CCTK_PostPostInitial (which scheduled is AFTER GRHydro_Rho_Minima_Setup_Final which is where grhydro_rho_min changes from its ID value to its evolution value).

  4. Frank Löffler reporter
    • removed comment

    I would expect this not to work when a hot atmosphere is evolved since calculating the Prim2Con already includes a devision by eps, which at this point might be 0. However, the current approach might also not catch this.

    GRHydro_InitialAtmosphereReset only resets the primitive variables. Scheduling it right before Con2Prim looks suspicious to me. It should instead be called in HydroBase_Prim2ConInitial before Primitive2Conservative(Poly)Cells(M). Problem with that would be that GRHydro_Rho_Minima_Setup_Final didn't run yet at that point, so a relative atmosphere value wouldn't yet be known. Possible solutions would be to run GRHydro_Rho_Minima* earlier, or to do a Prim2Con again after GRHydro_InitialAtmosphereReset.

  5. Frank Löffler reporter
    • removed comment

    Not taking technical difficulties into account, I would expect something like the following schedule:

    - setup ID (prim only) - find relative density minimum (if required, evolution thorn) - apply atmosphere (evolution thorn) - prim2con (evolution thorn)

    At this point, everything should be setup for evolution. GRHydro does an additional con2prim after that, but that shouldn't hurt.

  6. Roland Haas
    • removed comment

    That schedule sounds reasonable to me. The relative density minimum would be computed based on the density maximum on the grid (minimum would not make sense I think since the minimum on the grid at that time might be zero), would expect. Applying atmosphere would then mean to reset all points where rho<grhydro_rho_min to grhydro_rho_min (using the same logic as in AtmosphereReset). I would avoid setting the atmosphere mask since atmosphere_mask is currently (somewhat counter-intuitive) mostly used to indicate evolution failures and not which points have density <= grhydro_rho (which is what I would atmosphere to mean).

    I'd be very happy if the ID schedule became less convoluted.

  7. Frank Löffler reporter
    • removed comment

    After looking at the schedule and in particular how the loops over the refinement levels in Carpet are handled in INITIAL, I see no good way of achieving what we want without circumventing the Cactus/Carpet scheduling by looping over the refinement levels by hand inside a routine scheduled as global-late in INITIAL. Inside that routine we would have to - find the global density maximum - loop over levels - set the atmosphere everywhere - call prim2con everywhere

    Any comments on possible better solutions (that don't involve the new scheduler)?

  8. Erik Schnetter
    • removed comment

    Why does the following not work: - schedule a global-late routine to find the density maximum - schedule a global-late loop-local routine to set the atmosphere and call prim2con

  9. Roland Haas
    • removed comment

    prim2con has to run before FillTimeLevels (in Initialize.cc) runs since it computes the conservatives that are evolved by GRHydro. So prim2con has to run in local/level mode (or global-early, loop-local which is probably also impossible for some reason).

  10. Frank Löffler reporter
    • removed comment

    Yes, at global-early, loop-local the primitives are not set yet. The ID thorns do that only after this and we don't want to change their schedule.

  11. Frank Löffler reporter
    • removed comment

    prim2con could run after FillTimeLevels, but then it would have to manually work on all timelevels. This wouldn't be the only function which would need to understand that, so this isn't really feasable as well.

    We need to fix this, as it is currently really broken. The only way to have a 'nice' (clean, small) 'fix' would be to depreciate the parameter for a relative atmosphere. However, this would depreciate more than just one parameter, and it would change default behavior, making that problematic as well.

    Essentially a lot of functions/thorns are called in Initial/PostInitial and need to work on all levels/timelevels. In between some of them we would need to set the atmosphere, but we cannot know the value for that (if set to relative) before all levels have been set up. I don't see a (nice) way to implement that without changes to Cactus/Carpet atm.

    We could implement a partial fix by scheduling atmosphere setting in case an absolute atmosphere is requested. This is straightforward. However, it would still leave relative atmospheres not working.

  12. Erik Schnetter
    • removed comment

    As I understand, FillTimelevels is the last thing that is run after setting up initial data. Anything that is scheduled after FillTimelevels would be scheduled in an evolution time bin, i.e. in a time bin that is also schedule during evolution, and thus isn't really an "initial data" routine.

    Since FillTimelevels is rather cheap, we could re-run it on all refinement levels instead of just the current level. Would this help? We could do this on those levels that have do_{early,late}_global_mode set.

  13. Frank Löffler reporter
    • removed comment

    Unfortunately, FillTime is only part of the problem. But you might be right that in this particular case it might work. What about this:

    • INITIAL, level mode : initial data setup
    • INITIAL, global-late : density reduction and relative atmosphere value setup
    • INITIAL, global-late, loop-local: everything that needs the primitive variables, like - actual atmosphere setup (has to come first) - multipatch transformation - prim2con

    In particular, this would change the schedule of HydroBase_Prim2ConInitial

    Would this need an additional call to FillTimeLevels or not?

  14. Roland Haas
    • removed comment

    The solution in comment:14 is rather awkward (though it might work...) I find since it requires every user thorn routine in PostInitial (or initial) that wants to use matter variables to be global, loop-local. In particular if there is say a boundary call involved (think MoL_PostStep) this seems dangerous. One thorn affected would be eg. Noise.

  15. Frank Löffler reporter
    • removed comment

    I updated the patch to additionally schedule GRHydro_InitialAtmosphereReset AT CCTK_Initial AFTER HydroBase_Initial BEFORE HydroBase_Prim2ConInitial.

    This doesn't take care of relative atmosphere settings (which are 'sort of' dealt with later, nothing changed there), but it fixes the problem for absolute settings. With this patch all time levels on all refinement levels get an initial atmosphere before any prim2con/con2prim happens (which would/might otherwise fail).

  16. Roland Haas
    • changed status to open
    • removed comment

    This seems to cause GRHydor to use unitialized data unless the ID thorn sets all timelevels. TOVSolver for example does not do so by default. The best solution is likely to have GRHydro::GRHydro_InitialAtmosphereReset only loop over the current timelevel.

    This may require special treatment to get the atmosphere settings re-applied after ID is done and the relative density maximum has been found.

    Unfortunately scheduling around initial data is currently the most complex part of our schedule.

  17. Log in to comment