Merge pull request yt-project#4782 from jzuhone/gamer_sr_update

Updates for GAMER SR fields

doc/source/examining/loading_data.rst

@@ -1868,12 +1868,20 @@ are supported, and the following fields are defined:
 * ``("gas", "four_velocity_[txyz]")``: Four-velocity fields :math:`U_t, U_x, U_y, U_z`
 * ``("gas", "lorentz_factor")``: Lorentz factor :math:`\gamma = \sqrt{1+U_iU^i/c^2}`
   (where :math:`i` runs over the spatial indices)
+* ``("gas", "specific_reduced_enthalpy")``: Specific reduced enthalpy :math:`\tilde{h} = \epsilon + p/\rho`
+* ``("gas", "specific_enthalpy")``: Specific enthalpy :math:`h = c^2 + \epsilon + p/\rho`
 
 These, and other fields following them (3-velocity, energy densities, etc.) are
 computed in the same manner as in the
 `GAMER-SR paper <https://ui.adsabs.harvard.edu/abs/2021MNRAS.504.3298T/abstract>`_
 to avoid catastrophic cancellations.
 
+All of the special relativistic fields will only be available if the ``Temp`` and
+``Enth`` fields are present in the dataset, which can be ensured if the runtime
+options ``OPT__OUTPUT_TEMP = 1`` and ``OPT__OUTPUT_ENTHALPY  = 1`` are set in the
+``Input__Parameter`` file when running the simulation. This greatly speeds up
+calculations of the above derived fields in yt.
+
 .. rubric:: Caveats
 
 * GAMER data in raw binary format (i.e., ``OPT__OUTPUT_TOTAL = "C-binary"``) is not

tests/tests.yaml

@@ -76,7 +76,7 @@ answer_tests:
     - yt/frontends/gadget/tests/test_outputs.py:test_bigendian_field_access
     - yt/frontends/gadget/tests/test_outputs.py:test_magneticum
 
-  local_gamer_011:  # PR 3971
+  local_gamer_012:  # PR 4782
     - yt/frontends/gamer/tests/test_outputs.py:test_jet
     - yt/frontends/gamer/tests/test_outputs.py:test_psiDM
     - yt/frontends/gamer/tests/test_outputs.py:test_plummer

yt/frontends/gamer/cfields.pyx

@@ -7,102 +7,49 @@ cimport numpy as np
 import numpy as np
 
 
-cdef np.float64_t h_eos_tb(np.float64_t kT, np.float64_t g) noexcept nogil:
-    cdef np.float64_t x
-    x = 2.25 * kT * kT
-    return 2.5 * kT + x / (1.0 + math.sqrt(x + 1.0))
-
-cdef np.float64_t h_eos(np.float64_t kT, np.float64_t g) noexcept nogil:
-    return g * kT / (g - 1.0)
-
 cdef np.float64_t gamma_eos(np.float64_t kT, np.float64_t g) noexcept nogil:
     return g
 
 cdef np.float64_t gamma_eos_tb(np.float64_t kT, np.float64_t g) noexcept nogil:
     cdef np.float64_t x, c_p, c_v
-    x = 2.25 * kT / math.sqrt(2.25 * kT * kT + 1.0)
+    x = 2.25 * kT / ( math.sqrt(2.25 * kT * kT + 1.0) + 1.0 )
     c_p = 2.5 + x
     c_v = 1.5 + x
     return c_p / c_v
 
-cdef np.float64_t cs_eos_tb(np.float64_t kT, np.float64_t c, np.float64_t g) noexcept nogil:
+cdef np.float64_t cs_eos_tb(np.float64_t kT, np.float64_t h, np.float64_t g) noexcept nogil:
     cdef np.float64_t hp, cs2
-    hp = h_eos_tb(kT, 0.0) + 1.0
+    hp = h + 1.0
     cs2 = kT / (3.0 * hp)
     cs2 *= (5.0 * hp - 8.0 * kT) / (hp - kT)
-    return c * math.sqrt(cs2)
+    return math.sqrt(cs2)
 
-cdef np.float64_t cs_eos(np.float64_t kT, np.float64_t c, np.float64_t g) noexcept nogil:
+cdef np.float64_t cs_eos(np.float64_t kT, np.float64_t h, np.float64_t g) noexcept nogil:
     cdef np.float64_t hp, cs2
-    hp = h_eos(kT, g) + 1.0
+    hp = h + 1.0
     cs2 = g / hp * kT
-    return c * math.sqrt(cs2)
+    return math.sqrt(cs2)
+
 
 ctypedef np.float64_t (*f2_type)(np.float64_t, np.float64_t) noexcept nogil
 ctypedef np.float64_t (*f3_type)(np.float64_t, np.float64_t, np.float64_t) noexcept nogil
 
+
 cdef class SRHDFields:
-    cdef f2_type h
     cdef f2_type gamma
     cdef f3_type cs
-    cdef np.float64_t _gamma, _c, _c2
+    cdef np.float64_t _gamma
 
-    def __init__(self, int eos, np.float64_t gamma, np.float64_t clight):
+    def __init__(self, int eos, np.float64_t gamma):
         self._gamma = gamma
-        self._c = clight
-        self._c2 = clight * clight
         # Select aux functions based on eos no.
         if eos == 1:
-            self.h = h_eos
             self.gamma = gamma_eos
             self.cs = cs_eos
         else:
-            self.h = h_eos_tb
             self.gamma = gamma_eos_tb
             self.cs = cs_eos_tb
 
-    cdef np.float64_t _lorentz_factor(
-            self,
-            np.float64_t rho,
-            np.float64_t mx,
-            np.float64_t my,
-            np.float64_t mz,
-            np.float64_t kT,
-        ) noexcept nogil:
-        cdef np.float64_t u2
-        cdef np.float64_t fac
-
-        fac = (1.0 / (rho * (self.h(kT, self._gamma) + 1.0))) ** 2
-        u2 = (mx * mx + my * my + mz * mz) * fac
-        return math.sqrt(1.0 + u2 / self._c2)
-
-    cdef np.float64_t _four_vel(
-        self,
-        np.float64_t rho,
-        np.float64_t kT,
-        np.float64_t mi
-    ) noexcept nogil:
-        return mi / (rho * (self.h(kT, self._gamma) + 1.0))
-
-    @cython.boundscheck(False)
-    @cython.wraparound(False)
-    @cython.cdivision(True)
-    def lorentz_factor(self, dens, momx, momy, momz, temp):
-        cdef np.float64_t[:] rho = dens.ravel()
-
-        cdef np.float64_t[:] mx = momx.ravel()
-        cdef np.float64_t[:] my = momy.ravel()
-        cdef np.float64_t[:] mz = momz.ravel()
-        cdef np.float64_t[:] kT = temp.ravel()
-
-        out = np.empty_like(dens)
-        cdef np.float64_t[:] outp = out.ravel()
-        cdef int i
-
-        for i in range(outp.shape[0]):
-            outp[i] = self._lorentz_factor(rho[i], mx[i], my[i], mz[i], kT[i])
-        return out
-
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
@@ -117,156 +64,124 @@ cdef class SRHDFields:
             outp[i] = self.gamma(kT[i], self._gamma)
         return out
 
-    @cython.boundscheck(False)
-    @cython.wraparound(False)
-    @cython.cdivision(True)
-    def sound_speed(self, temp):
-        cdef np.float64_t[:] kT = temp.ravel()
-        out = np.empty_like(kT)
-        cdef np.float64_t[:] outp = out.ravel()
+    cdef np.float64_t _lorentz_factor(
+            self,
+            np.float64_t rho,
+            np.float64_t mx,
+            np.float64_t my,
+            np.float64_t mz,
+            np.float64_t h,
+        ) noexcept nogil:
+        cdef np.float64_t u2
+        cdef np.float64_t fac
 
-        cdef int i
-        for i in range(outp.shape[0]):
-            outp[i] = self.cs(kT[i], self._c, self._gamma)
-        return out
+        fac = (1.0 / (rho * (h + 1.0))) ** 2
+        u2 = (mx * mx + my * my + mz * mz) * fac
+        return math.sqrt(1.0 + u2)
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
-    def four_velocity_xyz(self, mom, dens, temp):
+    def lorentz_factor(self, dens, momx, momy, momz, enth):
         cdef np.float64_t[:] rho = dens.ravel()
-        cdef np.float64_t[:] mi = mom.ravel()
-        cdef np.float64_t[:] kT = temp.ravel()
-
-        out = np.empty_like(dens)
-        cdef np.float64_t[:] outp = out.ravel()
-        cdef int i
-
-        for i in range(outp.shape[0]):
-            outp[i] = self._four_vel(rho[i], kT[i], mi[i])
-        return out
 
-    @cython.boundscheck(False)
-    @cython.wraparound(False)
-    @cython.cdivision(True)
-    def velocity_xyz(self, dens, momx, momy, momz, temp, momi):
-        cdef np.float64_t[:] rho = dens.ravel()
         cdef np.float64_t[:] mx = momx.ravel()
         cdef np.float64_t[:] my = momy.ravel()
         cdef np.float64_t[:] mz = momz.ravel()
-        cdef np.float64_t[:] kT = temp.ravel()
-        cdef np.float64_t[:] mi = momi.ravel()
+        cdef np.float64_t[:] h = enth.ravel()
 
         out = np.empty_like(dens)
         cdef np.float64_t[:] outp = out.ravel()
-        cdef np.float64_t lf
         cdef int i
 
         for i in range(outp.shape[0]):
-            lf = self._lorentz_factor(rho[i], mx[i], my[i], mz[i], kT[i])
-            outp[i] = self._four_vel(rho[i], kT[i], mi[i]) / lf
+            outp[i] = self._lorentz_factor(rho[i], mx[i], my[i], mz[i], h[i])
         return out
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
-    def density(self, dens, momx, momy, momz, temp):
-        cdef np.float64_t[:] rho = dens.ravel()
-
-        cdef np.float64_t[:] mx = momx.ravel()
-        cdef np.float64_t[:] my = momy.ravel()
-        cdef np.float64_t[:] mz = momz.ravel()
+    def sound_speed(self, temp, enth):
         cdef np.float64_t[:] kT = temp.ravel()
-
-        out = np.empty_like(dens)
+        cdef np.float64_t[:] h = enth.ravel()
+        out = np.empty_like(kT)
         cdef np.float64_t[:] outp = out.ravel()
-        cdef np.float64_t lf
-        cdef int i
 
-        for i in range(outp.shape[0]):
-            lf = self._lorentz_factor(rho[i], mx[i], my[i], mz[i], kT[i])
-            outp[i] = rho[i] / lf
-        return out
-
-    @cython.boundscheck(False)
-    @cython.wraparound(False)
-    @cython.cdivision(True)
-    def pressure(self, dens, momx, momy, momz, temp):
-        cdef np.float64_t[:] rho = dens.ravel()
-        cdef np.float64_t[:] mx = momx.ravel()
-        cdef np.float64_t[:] my = momy.ravel()
-        cdef np.float64_t[:] mz = momz.ravel()
-        cdef np.float64_t[:] kT = temp.ravel()
-
-        out = np.empty_like(dens)
-        cdef np.float64_t[:] outp = out.ravel()
-        cdef np.float64_t lf
         cdef int i
-
         for i in range(outp.shape[0]):
-            lf = self._lorentz_factor(rho[i], mx[i], my[i], mz[i], kT[i])
-            outp[i] = rho[i] / lf * self._c2 * kT[i]
+            outp[i] = self.cs(kT[i], h[i], self._gamma)
         return out
 
+    cdef np.float64_t _four_vel(
+            self,
+            np.float64_t rho,
+            np.float64_t mi,
+            np.float64_t h,
+        ) noexcept nogil:
+        return mi / (rho * (h + 1.0))
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
-    def specific_thermal_energy(self, dens, momx, momy, momz, temp):
-        cdef np.float64_t[:] kT = temp.ravel()
+    def four_velocity_xyz(self, dens, mom, enth):
+        cdef np.float64_t[:] rho = dens.ravel()
+        cdef np.float64_t[:] mi = mom.ravel()
+        cdef np.float64_t[:] h = enth.ravel()
 
         out = np.empty_like(dens)
         cdef np.float64_t[:] outp = out.ravel()
         cdef int i
 
         for i in range(outp.shape[0]):
-            outp[i] = self._c2 * (self.h(kT[i], self._gamma)- kT[i])
+            outp[i] = self._four_vel(rho[i], mi[i], h[i])
         return out
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
-    def kinetic_energy_density(self, dens, momx, momy, momz, temp):
+    def kinetic_energy_density(self, dens, momx, momy, momz, temp, enth):
         cdef np.float64_t[:] rho = dens.ravel()
         cdef np.float64_t[:] mx = momx.ravel()
         cdef np.float64_t[:] my = momy.ravel()
         cdef np.float64_t[:] mz = momz.ravel()
         cdef np.float64_t[:] kT = temp.ravel()
+        cdef np.float64_t[:] h = enth.ravel()
 
         out = np.empty_like(dens)
         cdef np.float64_t[:] outp = out.ravel()
-        cdef np.float64_t lf, u2, hp, ux, uy, uz
+        cdef np.float64_t lf, u2, ux, uy, uz
         cdef int i
 
         for i in range(outp.shape[0]):
-            hp = self.h(kT[i], self._gamma) + 1.0
-            ux = self._four_vel(rho[i], kT[i], mx[i])
-            uy = self._four_vel(rho[i], kT[i], my[i])
-            uz = self._four_vel(rho[i], kT[i], mz[i])
+            ux = self._four_vel(rho[i], mx[i], h[i])
+            uy = self._four_vel(rho[i], my[i], h[i])
+            uz = self._four_vel(rho[i], mz[i], h[i])
             u2 = ux**2 + uy**2 + uz**2
-            lf = self._lorentz_factor(rho[i], mx[i], my[i], mz[i], kT[i])
-            gm1 = u2 / self._c2 / (lf + 1.0)
-            p = rho[i] / lf * self._c2 * kT[i]
-            outp[i] = gm1 * (rho[i] * hp * self._c2 + p)
+            lf = self._lorentz_factor(rho[i], mx[i], my[i], mz[i], h[i])
+            gm1 = u2 / (lf + 1.0)
+            p = rho[i] / lf * kT[i]
+            outp[i] = gm1 * (rho[i] * (h[i] + 1.0) + p)
         return out
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
-    def mach_number(self, dens, momx, momy, momz, temp):
+    def mach_number(self, dens, momx, momy, momz, temp, enth):
         cdef np.float64_t[:] rho = dens.ravel()
         cdef np.float64_t[:] mx = momx.ravel()
         cdef np.float64_t[:] my = momy.ravel()
         cdef np.float64_t[:] mz = momz.ravel()
         cdef np.float64_t[:] kT = temp.ravel()
+        cdef np.float64_t[:] h = enth.ravel()
 
         out = np.empty_like(dens)
         cdef np.float64_t[:] outp = out.ravel()
-        cdef np.float64_t cs
+        cdef np.float64_t cs, us, u
         cdef int i
 
         for i in range(outp.shape[0]):
-            cs = self.cs(kT[i], self._c, self._gamma)
-            us = cs / math.sqrt(1.0 - cs**2 / self._c2)
-            u = math.sqrt(mx[i]**2 + my[i]**2 + mz[i]**2) / (rho[i] * (self.h(kT[i], self._gamma) + 1.0))
+            cs = self.cs(kT[i], h[i], self._gamma)
+            us = cs / math.sqrt(1.0 - cs**2)
+            u = math.sqrt(mx[i]**2 + my[i]**2 + mz[i]**2) / (rho[i] * (h[i] + 1.0))
             outp[i] = u / us
         return out

yt/frontends/gamer/data_structures.py

@@ -255,7 +255,7 @@ def __init__(
         self.storage_filename = storage_filename
 
     def _set_code_unit_attributes(self):
-        if self.parameters["Opt__Unit"]:
+        if self.opt_unit:
             # GAMER units are always in CGS
             setdefaultattr(
                 self, "length_unit", self.quan(self.parameters["Unit_L"], "cm")
@@ -367,10 +367,11 @@ def _parse_parameter_file(self):
         else:
             self.mhd = 0
             self.srhd = 0
+            # set dummy value of mu here to avoid more complicated workarounds later
+            self.mu = 1.0
 
         # old data format (version < 2210) did not contain any information of code units
-        self.parameters.setdefault("Opt__Unit", 0)
-
+        self.opt_unit = self.parameters.get("Opt__Unit", 0)
         self.geometry = Geometry(geometry_parameters[parameters.get("Coordinate", 1)])
 
     @classmethod