diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 6658ee91c0..34cd6607c3 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -790,6 +790,12 @@ class RAMSESDataset(Dataset): _field_info_class = RAMSESFieldInfo gamma = 1.4 # This will get replaced on hydro_fn open + # RAMSES-specific parameters + force_cosmological: bool | None + _force_max_level: tuple[int, str] + _bbox: list[list[float]] | None + _self_shielding: bool | None = None + def __init__( self, filename, @@ -804,6 +810,7 @@ def __init__( max_level=None, max_level_convention=None, default_species_fields=None, + self_shielding=None, use_conformal_time=None, ): # Here we want to initiate a traceback, if the reader is not built. @@ -821,6 +828,10 @@ def __init__( cosmological: If set to None, automatically detect cosmological simulation. If a boolean, force its value. + + self_shielding: + If set to True, assume gas is self-shielded above 0.01 mp/cm^3. + This affects the fields related to cooling and the mean molecular weight. """ self._fields_in_file = fields @@ -887,6 +898,37 @@ def __init__( self.storage_filename = storage_filename + self.self_shielding = self_shielding + + @property + def self_shielding(self) -> bool: + if self._self_shielding is not None: + return self._self_shielding + + # Read namelist.txt file (if any) + has_namelist = self.read_namelist() + + if not has_namelist: + self._self_shielding = False + return self._self_shielding + + nml = self.parameters["namelist"] + + # "self_shielding" is stored in physics_params in older versions of the code + physics_params = nml.get("physics_params", default={}) + # and in "cooling_params" in more recent ones + cooling_params = nml.get("cooling_params", default={}) + + self_shielding = physics_params.get("self_shielding", False) + self_shielding |= cooling_params.get("self_shielding", False) + + self._self_shielding = self_shielding + return self_shielding + + @self_shielding.setter + def self_shielding(self, value): + self._self_shielding = value + @staticmethod def _sanitize_max_level(max_level, max_level_convention): # NOTE: the initialisation of the dataset class sets @@ -1091,28 +1133,34 @@ def caster(val): if self.num_groups > 0: self.group_size = rheader["ncpu"] // self.num_groups - # Read namelist.txt file (if any) self.read_namelist() - def read_namelist(self): + def read_namelist(self) -> bool: """Read the namelist.txt file in the output folder, if present""" namelist_file = os.path.join(self.root_folder, "namelist.txt") - if os.path.exists(namelist_file): - try: - with open(namelist_file) as f: - nml = f90nml.read(f) - except ImportError as e: - nml = f"An error occurred when reading the namelist: {str(e)}" - except (ValueError, StopIteration, AssertionError) as err: - # Note: f90nml may raise a StopIteration, a ValueError or an AssertionError if - # the namelist is not valid. - mylog.warning( - "Could not parse `namelist.txt` file as it was malformed:", - exc_info=err, - ) - return + if not os.path.exists(namelist_file): + return False + + try: + with open(namelist_file) as f: + nml = f90nml.read(f) + except ImportError as err: + mylog.warning( + "`namelist.txt` file found but missing package f90nml to read it:", + exc_info=err, + ) + return False + except (ValueError, StopIteration, AssertionError) as err: + # Note: f90nml may raise a StopIteration, a ValueError or an AssertionError if + # the namelist is not valid. + mylog.warning( + "Could not parse `namelist.txt` file as it was malformed:", + exc_info=err, + ) + return False - self.parameters["namelist"] = nml + self.parameters["namelist"] = nml + return True @classmethod def _is_valid(cls, filename: str, *args, **kwargs) -> bool: diff --git a/yt/frontends/ramses/fields.py b/yt/frontends/ramses/fields.py index 4f8ddea926..e1da99d73e 100644 --- a/yt/frontends/ramses/fields.py +++ b/yt/frontends/ramses/fields.py @@ -436,8 +436,16 @@ def create_cooling_fields(self) -> bool: def _create_field(name, interp_object, unit): def _func(field, data): shape = data["gas", "temperature_over_mu"].shape + # Ramses assumes a fraction X of Hydrogen within the non-metal gas. + # It has to be corrected by metallicity. + Z = data["gas", "metallicity"] + nH = ((1 - _Y) * (1 - Z) * data["gas", "density"] / mh).to("cm**-3") + if data.ds.self_shielding: + boost = np.maximum(np.exp(-nH / 0.01), 1e-20) + else: + boost = 1 d = { - "lognH": np.log10(_X * data["gas", "density"] / mh).ravel(), + "lognH": np.log10(nH / boost).ravel(), "logT": np.log10(data["gas", "temperature_over_mu"]).ravel(), } rv = interp_object(d).reshape(shape) diff --git a/yt/frontends/ramses/tests/test_outputs.py b/yt/frontends/ramses/tests/test_outputs.py index bc76f27d7e..e38a505a5d 100644 --- a/yt/frontends/ramses/tests/test_outputs.py +++ b/yt/frontends/ramses/tests/test_outputs.py @@ -1,4 +1,8 @@ import os +from contextlib import contextmanager +from pathlib import Path +from shutil import copytree +from tempfile import TemporaryDirectory import numpy as np from numpy.testing import assert_equal, assert_raises @@ -685,3 +689,108 @@ def _dummy_field(field, data): v1 = ad["gas", "test"] np.testing.assert_allclose(v0, v1) + + +# Simple context manager that overrides the value of +# the self_shielding flag in the namelist.txt file +@contextmanager +def override_self_shielding(root: Path, section: str, val: bool): + # Copy content of root in a temporary folder + with TemporaryDirectory() as tmp: + tmpdir = Path(tmp) / root.name + tmpdir.mkdir() + + # Copy content of `root` in `tmpdir` recursively + copytree(root, tmpdir, dirs_exist_ok=True) + + fname = Path(tmpdir) / "namelist.txt" + + with open(fname) as f: + namelist_data = f90nml.read(f).todict() + sec = namelist_data.get(section, {}) + sec["self_shielding"] = val + namelist_data[section] = sec + + with open(fname, "w") as f: + new_nml = f90nml.Namelist(namelist_data) + new_nml.write(f) + + yield tmpdir + + +@requires_file(ramses_new_format) +@requires_module("f90nml") +def test_self_shielding_logic(): + ################################################## + # Check that the self_shielding flag is correctly read + ds = yt.load(ramses_new_format) + + assert ds.parameters["namelist"] is not None + assert ds.self_shielding is False + + ################################################## + # Modify the namelist in-situ, reload and check + root_folder = Path(ds.root_folder) + for section in ("physics_params", "cooling_params"): + for val in (True, False): + with override_self_shielding(root_folder, section, val) as tmp_output: + # Autodetection should work + ds = yt.load(tmp_output) + assert ds.parameters["namelist"] is not None + assert ds.self_shielding is val + + # Manually set should ignore the namelist + ds = yt.load(tmp_output, self_shielding=True) + assert ds.self_shielding is True + + ds = yt.load(tmp_output, self_shielding=False) + assert ds.self_shielding is False + + ################################################## + # Manually set should work, even if namelist does not contain the flag + ds = yt.load(ramses_new_format, self_shielding=True) + assert ds.self_shielding is True + + ds = yt.load(ramses_new_format, self_shielding=False) + assert ds.self_shielding is False + + +ramses_star_formation = "ramses_star_formation/output_00013" + + +@requires_file(ramses_star_formation) +@requires_module("f90nml") +def test_self_shielding_loading(): + ds_off = yt.load(ramses_star_formation, self_shielding=False) + ds_on = yt.load(ramses_star_formation, self_shielding=True) + + # We do not expect any significant change at densities << 1e-2 + reg_on = ds_on.all_data().cut_region( + ["obj['gas', 'density'].to('mp/cm**3') < 1e-6"] + ) + reg_off = ds_off.all_data().cut_region( + ["obj['gas', 'density'].to('mp/cm**3') < 1e-6"] + ) + + np.testing.assert_allclose( + reg_on["gas", "cooling_total"], + reg_off["gas", "cooling_total"], + rtol=1e-5, + ) + + # We expect large differences from 1e-2 mp/cc + reg_on = ds_on.all_data().cut_region( + ["obj['gas', 'density'].to('mp/cm**3') > 1e-2"] + ) + reg_off = ds_off.all_data().cut_region( + ["obj['gas', 'density'].to('mp/cm**3') > 1e-2"] + ) + # Primordial cooling decreases with density (except at really high temperature) + # so we expect the boosted version to cool *less* efficiently + diff = ( + reg_on["gas", "cooling_primordial"] - reg_off["gas", "cooling_primordial"] + ) / reg_off["gas", "cooling_primordial"] + assert (np.isclose(diff, 0, atol=1e-5) | (diff < 0)).all() + + # Also make sure the difference is large for some cells + assert (np.abs(diff) > 0.1).any()