From 788241b3e03462fc5df4f6cba1dbfa6e3544342e Mon Sep 17 00:00:00 2001 From: John ZuHone Date: Sun, 12 Nov 2023 23:53:49 -0500 Subject: [PATCH 1/6] First definitions of cosmic ray fields --- yt/frontends/gamer/data_structures.py | 1 + yt/frontends/gamer/fields.py | 18 ++++++++++++++++++ 2 files changed, 19 insertions(+) diff --git a/yt/frontends/gamer/data_structures.py b/yt/frontends/gamer/data_structures.py index 6100bdf091a..728f9dd08e5 100644 --- a/yt/frontends/gamer/data_structures.py +++ b/yt/frontends/gamer/data_structures.py @@ -356,6 +356,7 @@ def _parse_parameter_file(self): # make aliases to some frequently used variables if parameters["Model"] == "Hydro": self.gamma = parameters["Gamma"] + self.gamma_cr = self.parameters.get("CR_Gamma", None) self.eos = parameters.get("EoS", 1) # Assume gamma-law by default # default to 0.6 for old data format self.mu = parameters.get( diff --git a/yt/frontends/gamer/fields.py b/yt/frontends/gamer/fields.py index 1c794677a57..fe182886ceb 100644 --- a/yt/frontends/gamer/fields.py +++ b/yt/frontends/gamer/fields.py @@ -23,6 +23,7 @@ class GAMERFieldInfo(FieldInfoContainer): ("MomY", (mom_units, ["momentum_density_y"], None)), ("MomZ", (mom_units, ["momentum_density_z"], None)), ("Engy", (erg_units, ["total_energy_density"], None)), + ("CRay", (erg_units, ["cosmic_ray_energy_density"], None)), ("Pote", (pot_units, ["gravitational_potential"], None)), # MHD fields on disk (CC=cell-centered) ("CCMagX", (b_units, [], "B_x")), @@ -289,6 +290,9 @@ def et(data): if self.ds.mhd: # magnetic_energy is a yt internal field Et -= data["gas", "magnetic_energy_density"] + if self.ds.gamma_cr is not None: + # cosmic rays are included in this dataset + Et -= data["gas", "cosmic_ray_energy_density"] return Et # thermal energy per mass (i.e., specific) @@ -334,6 +338,20 @@ def _thermal_energy_density(field, data): units=unit_system["pressure"], ) + if self.ds.gamma_cr is not None: + + def _cr_pressure(field, data): + return (data.ds.gamma_cr - 1.0) * data[ + "gas", "cosmic_ray_energy_density" + ] + + self.add_field( + ("gas", "cosmic_ray_pressure"), + _cr_pressure, + sampling_type="cell", + units=self.ds.unit_system["pressure"], + ) + # mean molecular weight if hasattr(self.ds, "mu"): From 3f96629c4925cbcf446bfbff3e160ac9e1eab992 Mon Sep 17 00:00:00 2001 From: John ZuHone Date: Wed, 15 Nov 2023 11:24:17 -0500 Subject: [PATCH 2/6] Fix problems with EOS detection --- yt/frontends/gamer/cfields.pyx | 16 ++++++++-------- yt/frontends/gamer/fields.py | 8 +++++--- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/yt/frontends/gamer/cfields.pyx b/yt/frontends/gamer/cfields.pyx index 76eb5c75d4a..dcdab90c0fd 100644 --- a/yt/frontends/gamer/cfields.pyx +++ b/yt/frontends/gamer/cfields.pyx @@ -7,7 +7,7 @@ cimport numpy as np import numpy as np -cdef np.float64_t h_eos4(np.float64_t kT, np.float64_t g) noexcept nogil: +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)) @@ -18,14 +18,14 @@ cdef np.float64_t h_eos(np.float64_t kT, np.float64_t g) noexcept nogil: cdef np.float64_t gamma_eos(np.float64_t kT, np.float64_t g) noexcept nogil: return g -cdef np.float64_t gamma_eos4(np.float64_t kT, np.float64_t g) noexcept nogil: +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) c_p = 2.5 + x c_v = 1.5 + x return c_p / c_v -cdef np.float64_t cs_eos4(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 c, np.float64_t g) noexcept nogil: cdef np.float64_t hp, cs2 hp = h_eos4(kT, 0.0) + 1.0 cs2 = kT / (3.0 * hp) @@ -52,14 +52,14 @@ cdef class SRHDFields: self._c = clight self._c2 = clight * clight # Select aux functions based on eos no. - if (eos == 4): - self.h = h_eos4 - self.gamma = gamma_eos4 - self.cs = cs_eos4 - else: + 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, diff --git a/yt/frontends/gamer/fields.py b/yt/frontends/gamer/fields.py index fe182886ceb..b5e4ac4a06b 100644 --- a/yt/frontends/gamer/fields.py +++ b/yt/frontends/gamer/fields.py @@ -62,10 +62,12 @@ def setup_fluid_fields(self): if self.ds.srhd: c2 = pc.clight * pc.clight c = pc.clight.in_units("code_length / code_time") - if self.ds.eos == 4: - fgen = SRHDFields(self.ds.eos, 0.0, c.d) - else: + if self.ds.eos == 1: + # gamma-law EOS fgen = SRHDFields(self.ds.eos, self.ds.gamma, c.d) + else: + # Taub-Mathews EOS + fgen = SRHDFields(self.ds.eos, 0.0, c.d) def _sound_speed(field, data): out = fgen.sound_speed(data["gamer", "Temp"].d) From 01bc2caf237bf301351ba5863b8019cba435f109 Mon Sep 17 00:00:00 2001 From: John ZuHone Date: Wed, 15 Nov 2023 11:36:35 -0500 Subject: [PATCH 3/6] Bugfix --- yt/frontends/gamer/cfields.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/frontends/gamer/cfields.pyx b/yt/frontends/gamer/cfields.pyx index dcdab90c0fd..3e87a15ae8e 100644 --- a/yt/frontends/gamer/cfields.pyx +++ b/yt/frontends/gamer/cfields.pyx @@ -27,7 +27,7 @@ cdef np.float64_t gamma_eos_tb(np.float64_t kT, np.float64_t g) noexcept nogil: 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 hp, cs2 - hp = h_eos4(kT, 0.0) + 1.0 + hp = h_eos_tb(kT, 0.0) + 1.0 cs2 = kT / (3.0 * hp) cs2 *= (5.0 * hp - 8.0 * kT) / (hp - kT) return c * math.sqrt(cs2) From a8f60b6dd56d888542103d64f16b2b9b53fb32a6 Mon Sep 17 00:00:00 2001 From: John ZuHone Date: Tue, 21 Nov 2023 09:41:29 -0500 Subject: [PATCH 4/6] Bugfix --- yt/frontends/gamer/fields.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/frontends/gamer/fields.py b/yt/frontends/gamer/fields.py index b5e4ac4a06b..b3eabdcb27a 100644 --- a/yt/frontends/gamer/fields.py +++ b/yt/frontends/gamer/fields.py @@ -292,7 +292,7 @@ def et(data): if self.ds.mhd: # magnetic_energy is a yt internal field Et -= data["gas", "magnetic_energy_density"] - if self.ds.gamma_cr is not None: + if getattr(self.ds, "gamma_cr", None): # cosmic rays are included in this dataset Et -= data["gas", "cosmic_ray_energy_density"] return Et @@ -340,7 +340,7 @@ def _thermal_energy_density(field, data): units=unit_system["pressure"], ) - if self.ds.gamma_cr is not None: + if getattr(self.ds, "gamma_cr", None): def _cr_pressure(field, data): return (data.ds.gamma_cr - 1.0) * data[ From 43cb8016d99f0e1050849813050bc93e332b3d15 Mon Sep 17 00:00:00 2001 From: John ZuHone Date: Tue, 21 Nov 2023 14:14:04 -0500 Subject: [PATCH 5/6] Simple test for cosmic rays --- yt/frontends/gamer/tests/test_outputs.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/yt/frontends/gamer/tests/test_outputs.py b/yt/frontends/gamer/tests/test_outputs.py index c7aae4a8ac9..87c70cb019b 100644 --- a/yt/frontends/gamer/tests/test_outputs.py +++ b/yt/frontends/gamer/tests/test_outputs.py @@ -125,3 +125,20 @@ def test_stress_energy(): else: Tmunu += p assert_array_almost_equal(sp[f"T{mu}{nu}"], Tmunu) + + +cr_shock = "CRShockTube/Data_000005" + + +@requires_ds(cr_shock) +def test_cosmic_rays(): + ds = data_dir_load(cr_shock) + assert_array_almost_equal(ds.gamma_cr, 4.0 / 3.0) + ad = ds.all_data() + p_cr = ad["gas", "cosmic_ray_pressure"] + e_cr = ad["gas", "cosmic_ray_energy_density"] + assert_array_almost_equal(p_cr, e_cr / 3.0) + e_kin = ad["gas", "kinetic_energy_density"] + e_int = ad["gas", "thermal_energy_density"] + e_tot = ad["gas", "total_energy_density"] + assert_array_almost_equal(e_tot, e_kin + e_int + e_cr) From 95b434a017cce77c768b0102806d0a589de36bf1 Mon Sep 17 00:00:00 2001 From: John ZuHone Date: Wed, 29 Nov 2023 23:06:53 -0500 Subject: [PATCH 6/6] Catch Lorentz factor and velocity fields that end up in the file if the correct options are turned on --- yt/frontends/gamer/fields.py | 56 ++++++++++++++++++++++++------------ 1 file changed, 37 insertions(+), 19 deletions(-) diff --git a/yt/frontends/gamer/fields.py b/yt/frontends/gamer/fields.py index b3eabdcb27a..6cbeda17eb1 100644 --- a/yt/frontends/gamer/fields.py +++ b/yt/frontends/gamer/fields.py @@ -109,15 +109,22 @@ def _four_velocity(field, data): ) # lorentz factor - def _lorentz_factor(field, data): - out = fgen.lorentz_factor( - data["gamer", "Dens"].d, - data["gamer", "MomX"].d, - data["gamer", "MomY"].d, - data["gamer", "MomZ"].d, - data["gamer", "Temp"].d, - ) - return data.ds.arr(out, "dimensionless") + if ("gamer", "Lrtz") in self.field_list: + + def _lorentz_factor(field, data): + return data["gamer", "Lrtz"] + + else: + + def _lorentz_factor(field, data): + out = fgen.lorentz_factor( + data["gamer", "Dens"].d, + data["gamer", "MomX"].d, + data["gamer", "MomY"].d, + data["gamer", "MomZ"].d, + data["gamer", "Temp"].d, + ) + return data.ds.arr(out, "dimensionless") self.add_field( ("gas", "lorentz_factor"), @@ -128,16 +135,27 @@ def _lorentz_factor(field, data): # velocity def velocity_xyz(v): - def _velocity(field, data): - out = fgen.velocity_xyz( - data["gamer", "Dens"].d, - data["gamer", "MomX"].d, - data["gamer", "MomY"].d, - data["gamer", "MomZ"].d, - data["gamer", "Temp"].d, - data["gamer", f"Mom{v.upper()}"].d, - ) - return data.ds.arr(out, "code_velocity").to(unit_system["velocity"]) + if ("gamer", f"Vel{v.upper()}") in self.field_list: + + def _velocity(field, data): + return data.ds.arr(data["gamer", f"Vel{v.upper()}"].d, "c").to( + unit_system["velocity"] + ) + + else: + + def _velocity(field, data): + out = fgen.velocity_xyz( + data["gamer", "Dens"].d, + data["gamer", "MomX"].d, + data["gamer", "MomY"].d, + data["gamer", "MomZ"].d, + data["gamer", "Temp"].d, + data["gamer", f"Mom{v.upper()}"].d, + ) + return data.ds.arr(out, "code_velocity").to( + unit_system["velocity"] + ) return _velocity