Sudden Density Decrease Before Cons-Prim Conversion

[hwhitehead]

I am experiencing a halting problem in my simulation where the timestep suddenly plummets, sometimes by 15 orders of magnitude in a single timestep. I have tracked this down to an excessively high velocity that appears within a single timestep. This is accompanied by a sudden drop in density.

My current belief is that for some reason a density drop is occuring before the cons-prim conversion which for a set momenta leads to incredibly high velocities. However none of my custom source terms decrement the density, they only update the momenta of gas.

I have checked the functionality of my source terms and they do not (to the best of my knowledge) suddenly boost the momenta of cells.I have turned off viscosity and am running without magnetic fields/relativity.

Obviously this is quite a specialised problem that is likely directly tied to my own custom additions. I just wanted to see if anyone else had experienced similar problems.

[hwhitehead]

Bumping this thread as the issue persists.

Further investigation into this issue has revealed that these plumetting densities are proceeded by a NEGATIVE density. This change does not occur during the user enrolled source functions, so I am at a loss as to what is causing it.

I have attached an image of a tracked failure. rho_min1 and v_max1 are the extremal (conserved) densities and velocity magnitudes before user source functions are applied. rho_min2 and v_max2 are the extremal values after user souce functions are applied.

You can see that the source functions do NOT change the conserved density explicitly, yet in the next time step the density is negative. They do change the conserved momenta, and therefore the velocity, but this should not lead to negative densities.

If anything like this has been experienced by other code users, advice would be much appreciated.

[c-white]

Negative densities and other violations of physical constraints are a fact of life for multidimensional simulations. Mass might simply be advected out of a cell in several directions in one step, such that the change in mass exceeds the original amount. This can happen even while obeying stability constraints like the CFL condition.

In one dimension it's easy enough to employ algorithms that avoid negative densities. However, with multiple dimensions, this requires particular care, with things like multidimensional Riemann solvers (Athena, like most codes, does not have these) and possibly quite small timesteps.

The code knows that a negative density will lead to nans, so it replaces anything below the floor with the floor (hydro/dfloor). Given how large a drop happened, I'm guessing there is no explicit floor set, so the code is using a very small default value. I think the first thing to try is raising the floor. A very rough rule of thumb is to have the floors be 1e-6 to 1e-8 times the typical values in the simulation. With a higher floor, a single cell getting a bad value won't cause the simulation to get too slow, and in a reasonable time some combination of physics, explicit diffusion, and numerical diffusion will fix the cell.

[hwhitehead]

Reopening this thread as unfortunately density floors have been unable to improve the situation. The code seems unable to relax from its incredibly short timestep, the density remains at the floor and the velocity remains incredibly high. This is with a VL2 integrator and xorder = 2. I have tried setting a hard ceiling for the timestep and while this stops some of the errors, I am unsure how robust this is.

I have considered some potential fixes but I don't know how reliable/feasible they would be in Athena++.

1. Rollbacks: is it possible to repeat a cycle at a smaller timetep if an negative density occurs? I know there are registers for the hydro values, so perhaps it would be possible to redo an integration at a smaller timestep to avoid such errors. I don't know that this would prevent negative densities however.
2. Smoothing: while clearly unphysical, you could presumably remove errors by replacing erroneous cell values with the average of their neighbours. I'm considering this the worst case fix for the moment
3. Smarter timestepping: at the moment I set my minimum timestep by the maximal acceleration that the source terms can exert, such that $dt < 0.1 \sqrt{L/a}$. While I think this is stable on its own, it does not seem to ensure stability with diffusive physics. I there exists a robust methodology for selecting a timestep with general source terms I would be interested to know more.

I don't wish to take up too much of the moderators' time as I am aware that this has become an increasingly specific problem, but if there is further advice for combating this problem I would greatly appreciate it as it has become a significant problem in my research.

[c-white]

Rollbacks in some form can help some simulations, but they haven't been implemented in this version of the code. There is a fair amount of careful bookkeeping that has to be done in order to make this compatible with MPI.

Smoothing isn't a switch that can be turned on, but some people have implemented it for their own problems, so that might be a way forward.

I can't promise results, but it might help to have a better idea of what exactly the source terms are, and what physics is being included. Also, it would help to see the state of the simulation right around when problems develop.

[hwhitehead]

The two key source terms that may introduce errors are gravitational attraction and accretion onto moving particles within the simulation. These particles are tracked by storing their properites in the ruser_mesh_data arrays and integrating them with the simulation. This is done consistently across MPI ranks. The particles' gravitational well is softened and so this source is unable to accelerate gas to massive values. Similarly while accretion can decrement a cells density, it does so geometrically and will never reduce the density through zero on its own.

The problems first appear early in the simulation but can appear seemingly stochastically later as well. When the negative densities do occur, they do so in a single cell. I have had some luck with applying smoothing to cells that fail, lifting their density from the dfloor to the average of their neighbours. I am worried about what errors this might introduce into the code however.

The simulation is hydrodynamic with viscosity and an isothermal equation of state.

[leeyp2009]

Have you resolved this problem? I think a similar problem happens to me. The time step suddenly becomes very large at a later time with the fluid density being very low (close to the density floor of 1e-9). Actually, I also implement very similar source terms as yours, e.g., a gravitational body in the disk and accretion onto it. Thanks.

[hwhitehead]

Unfortunately I have been unable to resolve this problem as of yet. I have tried introducing smoothing in density and monentum space local to the error to reduce the impact of these errors when they do occur, but this routine seems to create issues of its own. I'd be interested to hear if you have managed to find a fix.