[Cluster] Negative pressure

[mfournier01]

Hey,

I'm currently starting my PhD project and am setting up a simulation starting from the cluster.in environment. I've slightly modified the default set up so that I start with a perturbed density field (basically a clumpy atmosphere). When I let the simulation run, everything goes fine (the gas is cooling down, and start collapsing due to gravity).

However at some point the pressure reach negative values in some cells for a reason I could not identify. I added some couts in the code to check for the values of the variables involved in the computation of the pressure, and it look like the energy density is balancing the kinetic term so that they cancel out (w_p being defined in the code as gm1*(u_e - e_k)) and the result is slightly negative.

Any idea where this could come from ? I work with periodic boundary conditions by the way.

Thanks a lot !

[BenWibking]

This is a standard problem with kinetic energy dominated fluid flows. Since the equations are solved in conservation law form using the total energy, the numerical solution of the equations does not guarantee that the internal energy (and therefore the pressure) is positive.

There are a few solutions. A common one is to fall back to a positivity-preserving method for the cells that end up with a negative pressure, e.g., a positivity-preserving Riemann solver with piecewise-constant reconstruction and forward Euler time integration. This is what the "first-order flux correction option" does, except that it does not change the time integration. You can enable that and that should help. The code for that is here: https://github.com/parthenon-hpc-lab/athenapk/blob/ea838bf6e876ff6314394cb69f45708985a6bc2e/src/hydro/hydro.cpp#L962