Dual Energy Formalism

Cholla has the capability to employ the dual-energy formalism in order to simulate high Mach number flows.

Brief Primer

To ensure the conservation of mass, momentum, and energy, the hydro solvers provided in Cholla evolve the conserved quantities $(\rho, \rho v_x, \rho v_y, \rho v_z, E)$. The bulk of the work during each timestep is to calculate the fluxes in each conserved quantites between every cell. The flux calculations require knowledge of the thermal pressure, which is given by $p = (\gamma - 1) (E-E_{\rm kinetic})$.

Problems arise in simulations with large Mach numbers (such as cosmological simulations). When the Mach number is very large, then $(E - E_{\rm kinetic})$ is very small. Numerical issue arise, since the thermal pressure linearly scales with this small difference between two very large numbers.

About the Dual Energy Formalism

The solution to this numerical problem is to use the "dual-energy formalism" (more details are provided in Bryan+ 2014). The core idea is to track an extra separately-advected "thermal energy" field at each cell-location, in addition to the total energy field and use this "thermal energy" field in cases where $(E - E_{\rm thermal})$ provides insufficient precision.

The dual-energy formalism is parameterized by 2 parameters $\eta_1$ and $\eta_2$. It's easiest to understand their meaning by discussing how they are used. The Bryan+ 2014 paper describes 2 main steps:

1. during a given timestep, when we want to compute thermal pressure, we compare quotient of the "thermal energy" field divided by $E$ to $\eta_1$.
   • When the ratio is smaller than $\eta_1$ we use the "thermal energy" field. When it exceeds $\eta_1$, we use $(E-E_{\rm kinetic})$
   • In effect, $\eta_1$ directly parameterizes the precision where the dual-energy formalism kicks in.
   • It's worth mentioning that running a simulation with $\eta_1=0$ is equivalent to running a simulation without the dual energy formalism.
2. Near the end of each timestep (after updating all the fields with fluxes & and any source terms), the "thermal energy" energy is optionally overwritten with the value taken from $(E-E_{\rm kinetic})$
   • to motivate this step, it's important to understand that when we separately advecting the "thermal energy" and add the $-p(\nabla \cdot {\bf v})\Delta t/ \rho$ source term, we are effectively ignoring the effects of shock heating
   • Consequently, we might want to overwrite the "thermal energy" if we want to capture the effects of shock heating.
   • The precise condition that dictates when we overwrite the "thermal energy" field involves a comparison of $\eta_2$ and the values in neighboring cells. When $\eta_2$ is too high, we would effectively exclude shock-heating from weaker shocks. When $eta_2$ is too low we may include spurious heating that is introduced by the truncation error of $(E-E_{\rm kinetic})$.
   • NOTE: Bryan+ 2014 call this step "synchronization" - we find that name somewhat confusing since it may imply a bidirectional update (updating both "thermal energy" and $E$)

In practice Cholla does something slightly different.

1. It implements step 1 exactly as it's described above
2. Step 2 is pretty much the same. However, we also overwrite the "thermal energy" field when the quotient of "thermal energy" divided by $E$ exceeds $\eta_1$.
3. We add an additional step, after step 2, where we overwrite $E$ with the sum of the "thermal energy" field and $E_{\rm kinetic}$.
   • this is useful for bookkeeping reasons throughout the rest of the codebase
   • Something like this step is commonly implemented in other simulation codes (for example, Enzo effectively does this in any configuration involving radiative cooling physics).

Configuring the Dual Energy Formalism

To use the dual-energy formalism, you just need to define the DE macro at compile-time

The DE_ETA_1 and DE_ETA_2 macros are automatically defined within Cholla. At the time of writing this page, Cholla sets DE_ETA_1 to different values based on whether it is configured in cosmology. In cosmological simulations, DE_ETA_1 is always set to 1. This means that the separately advected "thermal energy" field is always prioritized over $E$ when computing the thermal pressure.

At this time, the dual-energy formalism is NOT compatible with MHD