coarse_cons_ in ApplyPhysicalBoundariesOnCoarseLevel has zero densities, triggering density floor

[sanghyukmoon]

Prerequisite checklist

Place an X in between the brackets on each line as you complete these checks:

• [ X ] Did you check that the issue hasn't already been reported?
• [ X ] Did you check the documentation in the Wiki for an answer?
• [ X ] Are you running the latest version of Athena++?

Summary of issue

ConservedToPrimitive is called in various places, including Cons2Prim step in the time integrator, Mesh::Initialize, and ApplyPhysicalBoundariesOnCoarseLevel. In addition to converting a conserved variables into primitive variables, it also applies a density floor https://github.com/PrincetonUniversity/athena/blob/e34cd25ab4fda6e2026a4039d25e365bf4b496c3/src/eos/adiabatic_hydro.cpp#L62 to prevent too small or negative density that could arise in regions of strong rarefaction.

While using AMR, I found when ConservedToPrimitive is called by ApplyPhysicalBoundariesOnCoarseLevel, https://github.com/PrincetonUniversity/athena/blob/e34cd25ab4fda6e2026a4039d25e365bf4b496c3/src/bvals/bvals_refine.cpp#L388-L391 the density array in coarse_cons_ contains zeros, triggering the density floor. This zero density has nothing to do with strong rarefaction and seems to be simply not initialized.

As I am not familiar with the underlying code structure related to the AMR boundary conditions, I am not sure if this is a problem at all, but decided to report here just in case. The practical problem I am experiencing is that, when I do some bookkeeping of how much flooring has occurred, the above issue produces too many flooring event that seems to be unrelated to the actual "bad cells". I solicit your comment on this issue, whether this is a problem or not.

Steps to reproduce

I created a minimal working example, by modifying the linear_wave problem. The refinement condition is such that the AMR is triggered at ncycle=1, in the Meshblock owned by rank=0. https://github.com/sanghyukmoon/athena/tree/refinement-physical-bc

./configure.py --prob linear_wave -mpi
mpirun -np 2 bin/athena -i inputs/hydro/athinput.linear_wave3d 1>out.txt 2>err.txt

The stderr dump contains some diagnostic info when ConservedToPrimitive is called and when the flooring is triggerred.

ncycle = 0 [Mesh::Initialize] Calling ProlongateBoundaries!
ncycle = 0 [Mesh::Initialize] Calling ConservedToPrimitive!
ncycle = 0 [TimeIntegratorTaskList::Prolongation] Calling ProlongateBoundaries!
ncycle = 0 [TimeIntegratorTaskList::Primitives] Calling ConservedToPrimitive!
ncycle = 0 [TimeIntegratorTaskList::Prolongation] Calling ProlongateBoundaries!
ncycle = 0 [TimeIntegratorTaskList::Primitives] Calling ConservedToPrimitive!
ncycle = 1 [TimeIntegratorTaskList::Prolongation] Calling ProlongateBoundaries!
ncycle = 1 [TimeIntegratorTaskList::Primitives] Calling ConservedToPrimitive!
ncycle = 1 [TimeIntegratorTaskList::Prolongation] Calling ProlongateBoundaries!
ncycle = 1 [TimeIntegratorTaskList::Primitives] Calling ConservedToPrimitive!
ncycle = 2 [Mesh::Initialize] Calling ProlongateBoundaries!
ncycle = 2 [BoundaryValues::ProlongateBoundaries] Calling ApplyPhysicalBoundariesOnCoarseLevel!
ncycle = 2 [BoundaryValues::ApplyPhysicalBoundariesOnCoarseLevel] Calling ConservedToPrimitive!
[ConservedToPrimitive] rank 0: density floor applied. k=1 j=5 i=2 z=-1.8750000000000000e-01 y=1.3125000000000000e+00 x=1.8750000000000000e-01 old=0.0000000000000000e+00
[ConservedToPrimitive] rank 0: density floor applied. k=2 j=5 i=2 z=1.8750000000000000e-01 y=1.3125000000000000e+00 x=1.8750000000000000e-01 old=0.0000000000000000e+00
[ConservedToPrimitive] rank 0: density floor applied. k=3 j=5 i=2 z=5.6250000000000000e-01 y=1.3125000000000000e+00 x=1.8750000000000000e-01 old=0.0000000000000000e+00
[ConservedToPrimitive] rank 0: density floor applied. k=4 j=5 i=2 z=9.3750000000000000e-01 y=1.3125000000000000e+00 x=1.8750000000000000e-01 old=0.0000000000000000e+00
[ConservedToPrimitive] rank 0: density floor applied. k=5 j=1 i=2 z=1.3125000000000000e+00 y=-1.8750000000000000e-01 x=1.8750000000000000e-01 old=0.0000000000000000e+00
[ConservedToPrimitive] rank 0: density floor applied. k=5 j=2 i=2 z=1.3125000000000000e+00 y=1.8750000000000000e-01 x=1.8750000000000000e-01 old=0.0000000000000000e+00
[ConservedToPrimitive] rank 0: density floor applied. k=5 j=3 i=2 z=1.3125000000000000e+00 y=5.6250000000000000e-01 x=1.8750000000000000e-01 old=0.0000000000000000e+00
[ConservedToPrimitive] rank 0: density floor applied. k=5 j=4 i=2 z=1.3125000000000000e+00 y=9.3750000000000000e-01 x=1.8750000000000000e-01 old=0.0000000000000000e+00

[tomidakn]

Thank you for the nice Christmas present. I'll look into it.

[sanghyukmoon]

Thank you for the nice Christmas present. I'll look into it.

Now I am guilty :joy: Happy holidays everyone (including myself)!

[tomidakn]

Hi @sanghyukmoon,

I think this is fine - let me explain. Suppose a situation like this figure. We are now processing the boundaries of MeshBlock A. Neighbor MeshBlocks on the same level and coarser level send NGHOST (=2) cells to MeshBlock A. The data from the coarser MeshBlocks are directly stored in the Coarse Buffer (open symbols), while the data from the MeshBlocks on the same level are restricted into the Coarse Buffer (filled symbols). Then, to prolongate the cell marked with a red triangle, we only need cells in the cyan stencil. However, to make the region rectangular, we intentionally expand this region and apply ConservativeToPrimitive to the yellow cells. As some cells (which belong to the same level as MeshBlock A. In this particular example, i=2 and j=5, matching your example) are not filled, the density floor is triggered. But those values should not affect the results.

I am sorry to trigger your floor detector wrongly. The easiest cure to this is initializing the Coarse Buffer data with some positive value. Can you check if this is the case? If so, we do not have to fix this issue.

[sanghyukmoon]

Hi @tomidakn,

Thank you very much for your detailed explanation with the diagram! Now I better understand what is going on under the hood thanks to you. Indeed, the issue goes away when the coarse buffer is initialized with a value greater than density_floor.