Hydrostatic equilibrium

[mfournier01]

Hey,

I'm trying to set up an hydrostatic equilibrium so that its stable on scales of 1e8-1e9 yrs. I made some tries with the already implemented initial conditions provided in the input folder (hse.in), which is by default set with a time limit of 1e6 yrs.

When trying to extend it to a few 1e7 yrs, I quickly observe a square pattern appearing in the box at roughly 50 Myr, which leads to a dramatic collapse of the profile (see file attached). I tried to increase resolution, change the boundary conditions or some profile parameters but did not manage to extend stability much beyond 10 Myr. Can I have any hint on how you manage to sustain equilibrium on larger timescales ?

Thanks by advance !

https://github.com/parthenon-hpc-lab/athenapk/assets/142504260/8b3847e5-5d44-4950-8906-2c3a61cd7582

[pgrete]

pinging @forrestglines and @BenWibking here as they'll like have more insight

[BenWibking]

The square pattern is from an interaction with the boundary. There is something going wrong much before that. You can see a significant change to the density profile already at 10 Myr. I don't know offhand how many free-fall times that is, but that is a useful number to compute.

My guess is that something is wrong with the initial conditions.

[BenWibking]

What is the resolution of this simulation? Generally hydrostatic equilibrium will not hold if it is resolved with fewer than 10-20 cells per pressure scale height.

For marginally-resolved cases, equilibrium can be significantly perturbed by the limiters used in the spatial reconstruction, which requires a modification to the code to fix. But this can also be mitigated by increasing the resolution such that the pressure scale height is well resolved everywhere.

[mfournier01]

The input parameters are the exact ones in the hse.in file, I just extended the tlim to 100 Myr. The box width is 0.2 Mpc, the mesh is set to nx=64, with static refinement down to level 2 at the middle of the box. Boundary conditions are outflow type.

I just gave it an other try with periodic boundary conditions and it looks much better. I still see weird patterns emerging at the very center of the box though, I'll make some more tries by changing smoothing radius and increasing resolution to see if I can improve it.

Thanks for help anyway.

[BenWibking]

I think hse.in was probably created as a computationally inexpensive test case, not with parameters suitable for a production simulation.

However, outflow boundary conditions are generally a bad choice for hydrostatic problems. That should probably be fixed.

[mfournier01]

So to solve the issue definitively, problem was completely solved by 1) switching to periodic boundary conditions and 2) increasing resolution of the box (taking 256^3 cells gives stable profile several 100 Myr at least).

Thanks again !

[forrestglines]

Sorry, I've been busy at a conference last week so didn't respond. Since I wrote the HSE initial conditions (at least for galaxy clusters) I'm probably to blame for some issues that's arisen.

It's also been our experience that more resolution is needed to maintain HSE. This might be a Cartesian grid issue. When I'm done with the spherical coordinates and I've gone through LANL's process to publish those changes, we could revisit the issue ( although we won't have special treatment for the poles or origin).

Another quirk about the HSE initialization is that I integrate a separate curve for each meshblock, then linearly interpolating this curve to cell centers. This has the advantage of integrating a higher resolution curve for higher resolution meshblocks in the initial mesh hierarchy. Consequently, however, we're interpolating from different solutions to HSE in each meshblock, which might lead to discontinuities at meshblock boundaries, not just AMR transitions. I haven't seen this create any noticeable discontinuities in practice, however. We could revisit this by generating one high resolution solution for all meshblock and then using a better interpolation technique (like cubic spline). We can keep an eye on this.

One solution to the boundaries issue might be to fix the boundaries at HSE conditions, which I think would work as long as you don't have initial perturbations.