Merge pull request #4817 from V-Nathir/fix_ramses_temperature

[BUG] Fix missing mu in temperature calculation RAMSES

yt/frontends/ramses/fields.py

@@ -1,10 +1,12 @@
 import os
+import warnings
 from functools import partial
 
 import numpy as np
 
 from yt import units
 from yt._typing import KnownFieldsT
+from yt.fields.field_detector import FieldDetector
 from yt.fields.field_info_container import FieldInfoContainer
 from yt.frontends.ramses.io import convert_ramses_conformal_time_to_physical_age
 from yt.utilities.cython_fortran_utils import FortranFile
@@ -196,18 +198,45 @@ def star_age(field, data):
         )
 
     def setup_fluid_fields(self):
-        def _temperature(field, data):
+        def _temperature_over_mu(field, data):
             rv = data["gas", "pressure"] / data["gas", "density"]
             rv *= mass_hydrogen_cgs / boltzmann_constant_cgs
             return rv
 
+        self.add_field(
+            ("gas", "temperature_over_mu"),
+            sampling_type="cell",
+            function=_temperature_over_mu,
+            units=self.ds.unit_system["temperature"],
+        )
+        found_cooling_fields = self.create_cooling_fields()
+
+        if found_cooling_fields:
+
+            def _temperature(field, data):
+                return data["gas", "temperature_over_mu"] * data["gas", "mu"]
+
+        else:
+
+            def _temperature(field, data):
+                if not isinstance(data, FieldDetector):
+                    warnings.warn(
+                        "Trying to calculate temperature but the cooling tables "
+                        "couldn't be found or read. yt will return T/µ instead of "
+                        "T — this is equivalent to assuming µ=1.0. To suppress this, "
+                        "derive the temperature from temperature_over_mu with "
+                        "some values for mu.",
+                        category=RuntimeWarning,
+                        stacklevel=1,
+                    )
+                return data["gas", "temperature_over_mu"]
+
         self.add_field(
             ("gas", "temperature"),
             sampling_type="cell",
             function=_temperature,
             units=self.ds.unit_system["temperature"],
         )
-        self.create_cooling_fields()
 
         self.species_names = [
             known_species_names[fn]
@@ -373,7 +402,8 @@ def _photon_flux(field, data):
                     units=flux_unit,
                 )
 
-    def create_cooling_fields(self):
+    def create_cooling_fields(self) -> bool:
+        "Create cooling fields from the cooling files. Return True if successful."
         num = os.path.basename(self.ds.parameter_filename).split(".")[0].split("_")[1]
         filename = "%s/cooling_%05i.out" % (
             os.path.dirname(self.ds.parameter_filename),
@@ -382,15 +412,15 @@ def create_cooling_fields(self):
 
         if not os.path.exists(filename):
             mylog.warning("This output has no cooling fields")
-            return
+            return False
 
         # Function to create the cooling fields
         def _create_field(name, interp_object, unit):
             def _func(field, data):
-                shape = data["gas", "temperature"].shape
+                shape = data["gas", "temperature_over_mu"].shape
                 d = {
                     "lognH": np.log10(_X * data["gas", "density"] / mh).ravel(),
-                    "logT": np.log10(data["gas", "temperature"]).ravel(),
+                    "logT": np.log10(data["gas", "temperature_over_mu"]).ravel(),
                 }
                 rv = interp_object(d).reshape(shape)
                 if name[-1] != "mu":
@@ -425,7 +455,7 @@ def _func(field, data):
                         "This cooling file format is no longer supported. "
                         "Cooling field loading skipped."
                     )
-                    return
+                    return False
                 if var.size == n1 * n2:
                     tvals[tname] = {
                         "data": var.reshape((n1, n2), order="F"),
@@ -446,7 +476,7 @@ def _func(field, data):
             ["lognH", "logT"],
             truncate=True,
         )
-        _create_field(("gas", "mu"), interp, tvals["mu"]["unit"])
+        _create_field(("gas", "mu"), interp, "dimensionless")
 
         # Add the number density field, based on mu
         def _number_density(field, data):
@@ -504,3 +534,5 @@ def _net_cool(field, data):
             function=_net_cool,
             units=cooling_function_units,
         )
+
+        return True

yt/frontends/ramses/tests/test_outputs_pytest.py

@@ -55,3 +55,28 @@ def test_field_config_2(custom_ramses_fields_conf):
         assert ("ramses", f) in ds.field_list
     for f in custom_grav:
         assert ("gravity", f) in ds.field_list
+
+
+@requires_file(output_00080)
+@requires_file(ramses_new_format)
+def test_warning_T2():
+    ds1 = yt.load(output_00080)
+    ds2 = yt.load(ramses_new_format)
+
+    # Should not raise warnings
+    ds1.r["gas", "temperature_over_mu"]
+    ds2.r["gas", "temperature_over_mu"]
+
+    # We cannot read the cooling tables of output_00080
+    # so this should raise a warning
+    with pytest.warns(
+        RuntimeWarning,
+        match=(
+            "Trying to calculate temperature but the cooling tables couldn't be "
+            r"found or read\. yt will return T/µ instead of T.*"
+        ),
+    ):
+        ds1.r["gas", "temperature"]
+
+    # But this one should not
+    ds2.r["gas", "temperature"]

yt/testing.py

@@ -6,11 +6,10 @@
 import sys
 import tempfile
 import unittest
-from collections.abc import Mapping
 from functools import wraps
 from importlib.util import find_spec
 from shutil import which
-from typing import Callable, Union
+from typing import TYPE_CHECKING, Callable, Union
 from unittest import SkipTest
 
 import matplotlib
@@ -20,13 +19,17 @@
 from unyt.exceptions import UnitOperationError
 
 from yt._maintenance.deprecation import issue_deprecation_warning
-from yt._typing import AnyFieldKey
 from yt.config import ytcfg
 from yt.frontends.stream.data_structures import StreamParticlesDataset
 from yt.funcs import is_sequence
 from yt.loaders import load, load_particles
 from yt.units.yt_array import YTArray, YTQuantity
 
+if TYPE_CHECKING:
+    from collections.abc import Mapping
+
+    from yt._typing import AnyFieldKey
+
 if sys.version_info < (3, 10):
     from yt._maintenance.backports import zip