cholla-hydro / cholla

A GPU-based hydro code
https://github.com/cholla-hydro/cholla/wiki
MIT License
66 stars 32 forks source link

`Calc_dt_[123]D` ignores dual-energy #416

Open mabruzzo opened 1 month ago

mabruzzo commented 1 month ago

I realized yesterday that Calc_dt_1D, Calc_dt_2D, and Calc_dt_3D all ignore the dual-energy formalism.

It would probably be more correct to use the internal energy in the regime where kinetic energy strongly dominates over the total energy. I'm not sure whether actually translates to any issues in existing simulations, but it should be addressed at some point.

bvillasen commented 1 month ago

Isn't calc_dt called after the energies have been synchronized? maybe the first timestep is the exception

mabruzzo commented 1 month ago

Isn't calc_dt called after the energies have been synchronized?

Yes, but there are definitely some hypersonic regimes (primarily in cosmological simulations) where it would be more correct to use the internal energy (if the mach number is really high, then you are going to lose a lot of precision in the pressure that you use for this calculation).

(I just slightly rephrased what I said in my original comment. It was perhaps a little too strongly worded)

bvillasen commented 1 month ago

I don't follow; after the energies have been synchronized, the thermal energy you get from both cases is the same, and it's being chosen based on the eta_2 condition. So then, why would it matter which internal energy one you use for calc_dt ... Unless something is significantly changing the energies between the energy_syncronization call and the calc_dt call, I don't see what is the issue

bvillasen commented 1 month ago

Ok I should restate. From a general correctness point of view, I agree that a dual energy condition (eta_1 condition, I think) should probably be applied during calc_dt. What I was trying to say is that, at least for cosmological simulations where eta_2>eta_1, and the energies are synchronized before calc_dt, this change would not make a difference. But yeah, after some thought, it's probably a good idea to add the DE check to calc_dt, even if it's just for consistency.

mabruzzo commented 1 month ago

In a purely mathematical sense, sure that's correct. But, there is a numerical difference. If $E$ and $E{\rm kin}$ are both numbers of very comparable magnitude that far exceed the difference between the values, then computing pressure from $E-E{\rm kin}$ is not very numerically robust.

The second paragraph of section 4.1.1 of Bryan+ 2014, refers to this exact scenario as "a numerically disastrous situation".

mabruzzo commented 1 month ago

To clarify, the main point I am making here is that using could be throwing away a lot of the precision in the pressure value.

If $E{kin} >> e{int}$ (lets say it is several orders of magnitude larger) and we previously set $E=E{kin} + e{int}$, then $E-E{kin}$ won't be exactly the same as $e{int}$. $E-E{kin}$ fundamentally can't have as much precision as $e{int}$ in this scenario. I believe the official term for this phenomenon is catastrophic cancellation.

(As I said at the start, I'm not sure whether this actually translates to any real issues in existing simulations)

bvillasen commented 1 month ago

Ok, yeah I see your point now. Agree that it's a good idea to do DE check in the calc_dt In the other hand, if E_kin >> e_int then the gas velocity would be much higher than the sound speed, so delta_t wouldn't care much about the pressure, right? as long as the numerical pressure is positive, getting a more precise pressure value wouldn't make a significant difference in delta_t, or am I missing something?

mabruzzo commented 1 month ago

Yeah, that's a fair point -- it's probably not very important.

Although if the gas is relatively cool (namely pressure is on the lower side and density is on the higher side), that could be slightly more problematic. But I suspect that could be a pathological scenario