fixed dual-energy formalism synchronization.

[mabruzzo]

Background

A brief summary of the problem (for future generations), the dual-energy formalism (see the end of section 4.1.1 of the Enzo-paper for a summary) is parameterized by quantities that $\eta_1$ and $\eta_2$

• about $\eta_1$: The whole point of tracking the internal-energy field is to have extra precision for the thermal energy when the kinetic energy is a very large fraction of the total energy). Essentially, $\eta_1$ governs when this this extra precision becomes necessary. We don't don't bother using the internal-energy field unless $((Etot-KE) / Etot) < \eta_1$, where $\eta_1$ is usually 0.001.

• about $\eta_2$: you can think of this as a synchronization parameter that tells us when we should updating the internal energy-field field to use the value computed from the total energy field. This is most relevant for including the effects of shock-heating.

It's worth mentioning that in most descriptions of the dual-energy formalism, only $\eta_2$ is discussed in the context of synchronization. However, $\eta_1$ becomes relevant in Cholla during synchronization of the internal energy because immediately after we synchronize the internal energy field, we use recompute the total energy from the kinetic energy and the total energy (If we didn't recompute the total energy - we would not need to consider $\eta_1$ at this point).

What I changed

The logic I needed to change was choosing which energy is assigned to the Internal Energy field by the Select_Internal_Energy_3D kernel (& in the 2D and 1D versions) since the subsequent kernel that gets called synchronizes total energy with gas energy.

• Essentially, I just needed to change the condition checked by this if-statement:
  • Previously the line read: if (U_total / Emax > eta_2) {
  • Now, the line reads: if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) {

[bvillasen]

I expect @evaneschneider is going to want my opinion about this... here it is.

As mentioned, the dual energy condition is important in preventing numerical errors when the kinetic energy is much greater than the thermal energy, but it can also suppress shock heating if it is not calibrated correctly.

Currently, Cholla applies the eta_1 condition when computing the pressure from the conserved variables and the eta_2 condition when synchronizing the advected thermal energy with (E_total - E_kin).

Note that the current values of eta_1 and eta_2 (defined here) are:

#define DE_ETA_1 0.001 
#define DE_ETA_2 0.035 

The value of eta_2 has been calibrated to produce gas temperatures in cosmological simulations that match the expectation from combining the IGM temperature and the temperature of the virialized gas that collapses into dark matter halos. If you are interested, a discussion is presented in sections 2.3 and 3.4 of Villasenor et al. 2021. I found from varying eta_2 that if eta_2 is too small, then the gas gets spuriously heated to high temperatures very early after the simulation starts (maybe z~ 50 to z~ 20), which is incorrect. If eta_2 is too high, shock heating is overly suppressed, and the gas infalling into DM halos fails to heat up to its virial temperature.

Since the (U_total / E > eta_1) condition is much less restrictive than the (U_total / Emax > eta_2) because E <= E_max and more importantly eta_1 < eta_2. If we go ahead with the proposed change and allow (U_total / E > eta_1) to be sufficient criteria to choose (E_total - E_kin) over the advected internal energy, then I can almost guarantee that we will again see the early spurious heating in cosmological simulations.

If you decide to proceed with this change, I recommend leaving the current condition (only the eta_2 condition) when running cosmological simulations.

@mabruzzo, is there a particular reason for you to suggest this change? Have you observed unexpected behavior somewhere? If I recall correctly, from the Enzo paper, the eta_1 condition was not used for synchronization of the thermal energy but only used for the pressure calculation (It could be that I'm wrong about that; it's been many years since I worked on this)

[evaneschneider]

@bvillasen Thanks for your comments! A little more background: I asked @mabruzzo to make this change because we've been running some galaxy sims with supernova lately and seeing incorrect results inside the supernova bubbles -- in particular, with dual energy turned on and the current version of Select_Internal_Energy, the gas inside the bubbles ends up being cold, rather than hot. I found that when I turned dual energy off, the problem was fixed, and realized that the issue is the way that the total energy synchronization after the selection is currently applied.

I believe that in the original formulation of DE in Cholla, both eta criteria had to be met in order for the value of the total energy to be updated with internal energy value at the end of the function. As currently written, ONLY the eta_2 check is applied, which means that in any circumstance other than when the eta_2 check is met, the value of the total energy is updated to contain the advected internal energy instead of the conservatively calculated value. This I think is definitely wrong, since it applies all over the place if there are no nearby shocks, and has nothing to do with whether you actually need the additional precision (as Matthew noted above). It also gave rise to the cold gas in bubbles we saw above.

To summarize, I actually think the "or" statement in this corrected version should be an "and" -- you only update the total energy with the advected internal energy if both criteria are met, and you use the conservatively calculated value of the internal energy in all other cases. That is the version that I have in my branch (where things seem to be behaving correctly now), but I haven't done extensive tests, or looked at the cosmological case, which I imagine could be different.

That said, as Matthew points out, the discussion of dual energy in Enzo doesn't actually talk about synchronization in the context of eta_2. So it's possible that a better solution is to not synchronize the total energy with the internal energy at the end of the step -- basically let the advected internal energy do its thing at all times unless the eta_2 criterion is met, and just use eta_1 to figure out when to use that quantity elsewhere in the code, as @bvillasen currently has it written. My understanding from @mabruzzo is that is more similar to what Enzo actually does, and it would avoid the problem we are currently having with the total energy being updated to contain the wrong value of the internal energy.

Happy to hear more thoughts on this if anyone has them!

[mabruzzo]

Hi Bruno, thanks for reaching out! We definitely value your opinions on this!

I've actually got some experience with the dual-energy formalism (Several years back, I needed to implement it in the VL+CT integrator I had implemented in Enzo-E[^1] and I've spoken at length about the topic with Greg Bryan). When I came across your paper a little while back I found it very insightful! (Calibrating the value of $\eta_2$ off of the virial temperature was a super-clever idea!).

Unexpected behavior

Before saying anything else, let me address:

@mabruzzo, is there a particular reason for you to suggest this change? Have you observed unexpected behavior somewhere?

@evanschneider and I have been trying to debug issues with supernovae bubbles in galaxy-scale simulations. When running these simulations we were finding that the cavities formed by supernovae were filled with lots of cool gas. A little while back Evan realized that the dual-energy formalism was making the problem worse. More recently Evan realized that the application of $\eta_1$ during internal_energy-synchronization was in the original version of CGOLs and adding it back in seems to improve things.

Highlighting subtle differences

If I recall correctly, from the Enzo paper, the eta_1 condition was not used for synchronization of the thermal energy but only used for the pressure calculation (It could be that I'm wrong about that; it's been many years since I worked on this).

So your recollection is exactly right! As we've both established, the description in the Enzo-paper talks about tracking the the separately-advected internal energy, $e$, and the total energy, $E$. And indeed, the paper only talks about using $\eta_1$ in the context of the pressure calculation.

However, there is a subtle, important difference between the description in that paper and what Cholla actually does (and for that matter what Enzo does). The difference comes down to "synchronization".

In the Enzo-paper, the only discussion of "synchronization" is relates to updating the value of $e$ based on the comparison of $\rho_i (Ei - {\bf v}^2/2) / \max{i-1,i,i+1} (\rho E)$ with $\eta_2$.

• Their choice to use the word "synchronize" is perhaps a poor choice or words. To me it implies a bi-directional relationship. But they never talk about altering $E$ based on the value of $e$.
• In a pure-hydro sim, the only way $e$ can affect $E$ is by affecting the pressure that is used in flux calculations. Thus, the dual-energy formalism has no effect on a pure-hydro simulation when $\eta_1=0$ (they actually make this point in the paper).

Now, in practice, there is a strong temptation to somehow update $E$ during this "synchronization"-step (especially if we are applying floors elsewhere) to reflect the value of $e$. I'm now going to highlight 3 different variations of the modified "synchronization":

1. What Cholla currently does: Cholla's dev branch currently does the "$\eta_2$-neighbor-comparison" to update the value of $e$. Then we always recompute $E$ from $e$ and the kinetic energy.
   • In this case, even when $\eta_1$ is zero, a pure hydro simulation can still be impacted by the dual-energy formalism.
   • In other words, the value of $e$ will always affect the computed value of $p$. $\eta_1$ just controls whether you compute $p$ directly from $e$ or if you compute $p$ from $(e + 0.5 \rho v^2) - 0.5 \rho v^2$. This is almost equivalent to running simulations with $\eta_1=1$ (in fact, under the current system it would always be preferable to use $\eta_1=1$)
2. Proposed change in this PR: This PR currently proposes that when $(E-E_{\rm kinetic})/E > \eta1$ that we have adequate precision in the thermal-energy computed from $E$; so we assign $e$ to be $E-E{\rm kinetic}$. When that condition is not-satisfied, we fall back to using the "$\eta_2$-neighbor-comparison" to update the value of $e$.
   • In this case, when $\eta_1$ is zero, a pure hydro simulation won't be impacted by the dual-energy formalism.
3. What Enzo's ppm-solver traditionally [^2] did: the ppm solver updated $e$ based on $\eta_1$. Then immediately afterwards, they used the $\eta_1$-condition (that's normally used for computing pressure) to determine whether they should overwrite $E$ with $(e + 0.5 \rho v^2)$.
   • This approach has the nice feature of updating $E$ to reflect $e$ while still ensuring that $\eta_1$ still has the same meaning intended in the original description of the dual energy -formalism

Discussion of your concerns

Since the (U_total / E > eta_1) condition is much less restrictive than the (U_total / Emax > eta_2) because E <= E_max and more importantly eta_1 < eta_2. If we go ahead with the proposed change and allow (U_total / E > eta_1) to be sufficient criteria to choose (E_total - E_kin) over the advected internal energy, then I can almost guarantee that we will again see the early spurious heating in cosmological simulations.

This is a very good point! You make a compelling case that this may probably change the temperature structure of the cosmological simulation. I certainly don't want to do that.

With that said, this highlights a potential issue with how well the current choice of $\eta_2$ generalizes that should be explored in the future.

• In effect, $\eta_1$ can be thought of as a parameter controlling the precision where we think truncation error in $E$ will significantly alter the outcome of the simulations.
• $\eta_2$ primarily exists to balance spurious heating from the contamination of errors introduced during the advection of $E$ (due to truncation errors) while still accounting for shock-heating. However, if a cell and its neighbors (technically the proposed change doesn't check the neighbors) all satisfy the $\eta_1$ condition, then the dual-energy formalism shouldn't really be necessary[^3].

Thus, if the proposed change makes a big difference in the temperature of halos, this suggests that our calibrated $\eta_2$ suppresses other sources of "spurious" heating (beyond the truncation error).

• There's nothing inherently wrong with this - and it will work well as long as the simulation conditions match the conditions you used to calibrate $\eta_2$. But if this is true, I do think caution needs to be exercised when running simulations with different reconstruction-methods, different integrator-schemes (VL vs simple), different resolution.

Where to go from here

So I think we have 3 options: A. We do what you suggested and only apply the change in non-cosmological simulations

• if we go this route, in the long term we should confirm that this is equivalent to running cosmological simulations with the proposed changes and configuring DE_ETA_1 to be 1 (or some really large numberer). If that is indeed the case, we should do that.

B. We can go with the approach historically used in Enzo's ppm solver (variation number 3 from above). This doesn't let $\eta_1$ affect the value stored in $e$, but it would restore meaning to $\eta_1$ when computing pressure. (Overall it's more consistent with the)

C. We could go with the approach described in the Enzo-paper (i.e. we never overwrite the value of $E$ based on $e$)

I'm a little agnostic here. On the one hand, I like choices 2 or 3 a lot in principle (and they are somewhat battle-tested). On the other hand, the fact that Enzo-Classic doesn't really use choices 2 or 3 and hasn't for the last 12 years mutes my feelings on this

[^1]: As a complete aside: It's actually on my todo-list to go back and "fix" this implementation in Enzo-E -- the current version is more similar to the one described by Teyssier (2015) and we definitely see suppression of shock heating in our Zeldovich Pancake test problem.

[^2]: Specifically, I'm talking about the formulation used by Enzo-E's ppm-solver - which you find here. While the name suggests this is used to compute thermal pressure, this is explicitly called after updating the hydro quantities. You can find the identical file here in Enzo-Classic. The reason I'm linking to the Enzo-E version is it uses a slightly older implementation of the ppm-solver, which uses an implementation of the dual-energy formalism that is more consistent with the description in the Enzo-paper. For the last 12 years, Enzo-Classic does something different -- it's nearly equivalent to running the simulation with $\eta_1=1$. There are some other quirks in this implementation that make me question how intentional the changes in the implementation are (so I actually opened an issue on that repo to ask about it).

[^3]: It could arguably be more accurate to ignore the dual-energy formalism in these cases. After all, there are inaccuracies in the advection of $e$. It's not a conservative update and we don't use the time-centered pressure in the source-terms. Furthermore, it may be possible to improve the velocity divergence calculation -- rather than using the cell-centered values from neighboring cells OR the reconstructed values, we could use the interface values implicitly computed within the Riemann solver (from the reconstructed values)

EDIT: Changed enumeration of options from numbers to letters.

[evaneschneider]

I'll update my comment based on @mabruzzo's more extensive and better articulated explanation: The reason I prefer "and" to "or" in this check is that my thinking on dual energy in general is that it should only be necessary in a scenario where we lack precision in the internal energy calculation due to truncation error. So, if the eta_1 criterion isn't met, we should not need it (and definitely should not change the value of the total energy based on the non-conservative update). I guess that is most similar in philosophy to @mabruzzo's option 3, since it essentially imposes quite a strict set of requirements before you are allowed to overwrite the value in E. That said, if @bvillasen thinks that abandoning the synchronization at the end of the function (and therefore, only applying the eta_2 check to adjust the value of e) would make this fix more compatible with the cosmological sims, I'm happy for us to go that route.

[brantr]

IMO, we should do whatever is the most physically motivated. I agree that “and” seems preferable to “or”, but we need to check the behavior of the cosmological simulation test before we know whether this should be used generally or if the behavior should be optional.On Dec 11, 2023, at 6:47 AM, Evan Schneider @.***> wrote: I'll update my comment based on @mabruzzo more extensive and better articulated explanation: The reason I prefer "and" to "or" in this check is that my thinking on dual energy in general is that it should only be necessary in a scenario where we lack precision in the internal energy calculation due to truncation error. So, if the eta_1 criterion isn't met, we should not need it (and definitely should not change the value of the total energy based on the non-conservative update). I guess that is most similar in philosophy to @mabruzzo's option 3, since it essentially imposes quite a strict set of requirements before you are allowed to overwrite the value in E. That said, if @bvillasen thinks that abandoning the synchronization at the end of the function (and therefore, only applying the eta_2 check to adjust the value of e) would make this fix more compatible with the cosmological sims, I'm happy for us to go that route.

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you are subscribed to this thread.Message ID: @.***>

[mabruzzo]

About the current implementation: @evaneschneider I think you and I are conceptually on the same page, but I still think "or" is correct (it's possible I am misunderstanding the code).

In the currently proposed change:

if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) {
      U = U_total;
    } else {
      U = U_advected;
    }

We are saying here, that we overwrite internal_energy with the value computed from total_energy if EITHER:

• $E - {\bf v}^2/ / E > \eta_1$ -- the total energy has enough precision to make the shock-heating unnecessary
• $\rho (E - {\bf v}^2/2) / \max (\rho E) > \eta_2$ -- there's possibly nearby shock-heating that must be included

If we made that an "and" statement, then its possible that we won't overwrite the internal_energy in cases where there is strong shock heating (in that scenario the precision of the internal-energy computed from total_energy would be low)

[mabruzzo]

Just a note - I realized that the enumerated options don't match up with the enumerated variations of the dual-energy formalism.

Therefore I just changed the enumerated options to use letters instead of numbers (to try to avoid confusion)

[evaneschneider]

About the current implementation: @evaneschneider I think you and I are conceptually on the same page, but I still think "or" is correct (it's possible I am misunderstanding the code).

In the currently proposed change:

if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) {
      U = U_total;
    } else {
      U = U_advected;
    }

We are saying here, that we overwrite internal_energy with the value computed from total_energy if EITHER:

• E−v2//E>η1 -- the total energy has enough precision to make the shock-heating unnecessary
• ρ(E−v2/2)/max(ρE)>η2 -- there's possibly nearby shock-heating that must be included

If we made that an "and" statement, then its possible that we won't overwrite the internal_energy in cases where there is strong shock heating (in that scenario the precision of the internal-energy computed from total_energy would be low)

Ah yes, I see. I agree!

[brantr]

Could we have a meeting to discuss this issue?On Dec 11, 2023, at 7:33 AM, Evan Schneider @.***> wrote:

About the current implementation: @evaneschneider I think you and I are conceptually on the same page, but I still think "or" is correct (it's possible I am misunderstanding the code). In the currently proposed change: if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) { U = U_total; } else { U = U_advected; } We are saying here, that we overwrite internal_energy with the value computed from total_energy if EITHER:

E−v2//E>η1 -- the total energy has enough precision to make the shock-heating unnecessary ρ(E−v2/2)/max(ρE)>η2 -- there's possibly nearby shock-heating that must be included

If we made that an "and" statement, then its possible that we won't overwrite the internal_energy in cases where there is strong shock heating (in that scenario the precision of the internal-energy computed from total_energy would be low)

Ah yes, I see. I agree!

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you commented.Message ID: @.***>

[mabruzzo]

Sorry for bombarding you with messages --

Just to clarify:

• If we go with option B:
  • We would revert the implementation of Select_Internal_Energy_3D back to the current implementation found on dev.
  • We would modify Sync_Energies_3D to only overwrite total_energy with the sum of internal_energy and the kinetic energy when $(E - 0.5 {\bf v}^2)/ E < \eta_1$ (Yes, the < operation is intention).
• If we go with option C:
  • We would revert the implementation of Select_Internal_Energy_3D back to the current implementation found on dev.
  • We would delete Sync_Energies_3D (despite its name, that function only modifies total_energy)

Again, I'm fine with anything.

Maybe scheduling a meeting would be good.

[evaneschneider]

HI @mabruzzo, can we get this one finalized today so I can merge it in?