From ad985d1320b99ded28e880c605605f416e24e30b Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Tue, 5 Mar 2024 13:42:22 +0100 Subject: [PATCH 01/11] Add support for self-shielding in heating/cooling rates calculation --- yt/frontends/ramses/data_structures.py | 57 ++++++++++++++++++-------- yt/frontends/ramses/fields.py | 7 +++- 2 files changed, 46 insertions(+), 18 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 5b811cf07b2..1b308938cbb 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -795,6 +795,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, @@ -809,6 +815,7 @@ def __init__( max_level=None, max_level_convention=None, default_species_fields=None, + self_shielding=None, ): # Here we want to initiate a traceback, if the reader is not built. if isinstance(fields, str): @@ -880,6 +887,7 @@ def __init__( self.fluid_types += (FH.ftype,) self.storage_filename = storage_filename + self.self_shielding = self_shielding @staticmethod def _sanitize_max_level(max_level, max_level_convention): @@ -1111,27 +1119,42 @@ def caster(val): self.group_size = rheader["ncpu"] // self.num_groups # Read namelist.txt file (if any) - self.read_namelist() + has_namelist = self.read_namelist() + + if self.self_shielding is None and has_namelist: + nml = self.ds.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.self_shielding = physics_params.get("self_shielding", False) + self.self_shielding |= cooling_params.get("self_shielding", False) - 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 e: + nml = f"An error occurred when reading the namelist: {str(e)}" + 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 4c871ec08b8..3808b1590b5 100644 --- a/yt/frontends/ramses/fields.py +++ b/yt/frontends/ramses/fields.py @@ -418,8 +418,13 @@ def create_cooling_fields(self) -> bool: def _create_field(name, interp_object, unit): def _func(field, data): shape = data["gas", "temperature_over_mu"].shape + nH = _X * data["gas", "density"].to("mp/cm**3") / mh + 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) From a365c271a02045aa322e5f2d8d31a27c5a8a4ab3 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 19 Sep 2024 09:18:08 +0200 Subject: [PATCH 02/11] Correct densities by metallicity MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Hydrogen abundances are defined as (1 - Y) × (1 - Z), where Y is the hard-coded primordial Helium fraction and Z is metallicity. Co-authored-by: Víctor Rufo Pastor <67118212+V-Nathir@users.noreply.github.com> --- yt/frontends/ramses/data_structures.py | 2 +- yt/frontends/ramses/fields.py | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 1b308938cbb..8acf546cd31 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -799,7 +799,7 @@ class RAMSESDataset(Dataset): force_cosmological: bool | None _force_max_level: tuple[int, str] _bbox: list[list[float]] | None - self_shielding: bool | None = None + _self_shielding: bool | None = None def __init__( self, diff --git a/yt/frontends/ramses/fields.py b/yt/frontends/ramses/fields.py index 3808b1590b5..8a8940ed8cb 100644 --- a/yt/frontends/ramses/fields.py +++ b/yt/frontends/ramses/fields.py @@ -418,7 +418,13 @@ def create_cooling_fields(self) -> bool: def _create_field(name, interp_object, unit): def _func(field, data): shape = data["gas", "temperature_over_mu"].shape - nH = _X * data["gas", "density"].to("mp/cm**3") / mh + # Ramses assumes a fraction X of Hydrogen within the non-metal gas. + # It has to be corrected by metallicity. + if ("gas", "metallicity") in data.ds.derived_field_list: + Z = data["gas", "metallicity"] + else: + Z = 0 + nH = (1 - _Y) * (1 - Z) * data["gas", "density"] / mh if data.ds.self_shielding: boost = np.maximum(np.exp(-nH / 0.01), 1e-20) else: From 9b9e2e45f81a30b74545e8cce2114e07b0aca080 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 12 Apr 2024 09:44:58 +0200 Subject: [PATCH 03/11] Fix failing tests due to self_shielding attribute Co-authored-by: Chris Havlin --- yt/frontends/ramses/data_structures.py | 2 +- yt/frontends/ramses/fields.py | 5 +---- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 8acf546cd31..99805e2eb68 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -1122,7 +1122,7 @@ def caster(val): has_namelist = self.read_namelist() if self.self_shielding is None and has_namelist: - nml = self.ds.parameters["namelist"] + nml = self.parameters["namelist"] # "self_shielding" is stored in physics_params in older versions of the code physics_params = nml.get("physics_params", default={}) diff --git a/yt/frontends/ramses/fields.py b/yt/frontends/ramses/fields.py index 8a8940ed8cb..9bfff673c98 100644 --- a/yt/frontends/ramses/fields.py +++ b/yt/frontends/ramses/fields.py @@ -420,10 +420,7 @@ 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. - if ("gas", "metallicity") in data.ds.derived_field_list: - Z = data["gas", "metallicity"] - else: - Z = 0 + Z = data["gas", "metallicity"] nH = (1 - _Y) * (1 - Z) * data["gas", "density"] / mh if data.ds.self_shielding: boost = np.maximum(np.exp(-nH / 0.01), 1e-20) From 2ed7c9f59594adc9bc3a1083104d3edddeae2109 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 19 Sep 2024 09:42:02 +0200 Subject: [PATCH 04/11] Read self_shielding as a property --- yt/frontends/ramses/data_structures.py | 44 ++++++++++++++++++-------- 1 file changed, 31 insertions(+), 13 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 99805e2eb68..668b6502b26 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -887,8 +887,38 @@ def __init__( self.fluid_types += (FH.ftype,) 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 @@ -1118,19 +1148,7 @@ def caster(val): if self.num_groups > 0: self.group_size = rheader["ncpu"] // self.num_groups - # Read namelist.txt file (if any) - has_namelist = self.read_namelist() - - if self.self_shielding is None and has_namelist: - 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.self_shielding = physics_params.get("self_shielding", False) - self.self_shielding |= cooling_params.get("self_shielding", False) + self.read_namelist() def read_namelist(self) -> bool: """Read the namelist.txt file in the output folder, if present""" From dcdf38efef38f3614d8499be57b73072efc383ed Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 19 Sep 2024 09:42:16 +0200 Subject: [PATCH 05/11] More useful warning if f90nml is missing --- yt/frontends/ramses/data_structures.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 668b6502b26..d4af1290b7e 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -1159,8 +1159,11 @@ def read_namelist(self) -> bool: 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 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 From 106e22d3424ea370e1373093421621234debbade Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 19 Sep 2024 11:47:21 +0200 Subject: [PATCH 06/11] Test self-shielding --- yt/frontends/ramses/tests/test_outputs.py | 102 ++++++++++++++++++++++ 1 file changed, 102 insertions(+) diff --git a/yt/frontends/ramses/tests/test_outputs.py b/yt/frontends/ramses/tests/test_outputs.py index 064c57b78f8..ebeed7b7a0c 100644 --- a/yt/frontends/ramses/tests/test_outputs.py +++ b/yt/frontends/ramses/tests/test_outputs.py @@ -1,4 +1,5 @@ import os +from contextlib import contextmanager import numpy as np from numpy.testing import assert_equal, assert_raises @@ -683,3 +684,104 @@ 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(fname: str, section: str, val: bool): + with open(fname) as f: + old_content = f.read() + f.seek(0) + 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 + + with open(fname, "w") as f: + f.write(old_content) + + +@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 + namelist = os.path.join(ds.root_folder, "namelist.txt") + # Can be set in any of these two + for section in ("physics_params", "cooling_params"): + for val in (True, False): + with override_self_shielding(namelist, section, val): + # Autodetection should work + ds = yt.load(ramses_new_format) + assert ds.parameters["namelist"] is not None + assert ds.self_shielding is val + + # Manually set should ignore the namelist + 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 + + ################################################## + # 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() From 4774ae7981f38d426da17219ab635a510ab7a850 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 4 Oct 2024 09:15:52 +0200 Subject: [PATCH 07/11] Add docstring --- yt/frontends/ramses/data_structures.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index d4af1290b7e..d2325a6e7bf 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -832,6 +832,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 From 3d1ecb03249f25bc65ed61f64013c6f9e95d5f74 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 4 Oct 2024 09:16:34 +0200 Subject: [PATCH 08/11] Convert nH to 1/mp**3 --- yt/frontends/ramses/fields.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/frontends/ramses/fields.py b/yt/frontends/ramses/fields.py index 9bfff673c98..1b6c974bec6 100644 --- a/yt/frontends/ramses/fields.py +++ b/yt/frontends/ramses/fields.py @@ -421,7 +421,7 @@ def _func(field, data): # 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 + 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: From abda72c382d59356df651cbf8512a7427cb14fb7 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 4 Oct 2024 09:22:29 +0200 Subject: [PATCH 09/11] Work in temporary directory --- yt/frontends/ramses/tests/test_outputs.py | 46 +++++++++++++---------- 1 file changed, 26 insertions(+), 20 deletions(-) diff --git a/yt/frontends/ramses/tests/test_outputs.py b/yt/frontends/ramses/tests/test_outputs.py index ebeed7b7a0c..f2a552d750a 100644 --- a/yt/frontends/ramses/tests/test_outputs.py +++ b/yt/frontends/ramses/tests/test_outputs.py @@ -1,5 +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 @@ -689,23 +692,28 @@ def _dummy_field(field, data): # Simple context manager that overrides the value of # the self_shielding flag in the namelist.txt file @contextmanager -def override_self_shielding(fname: str, section: str, val: bool): - with open(fname) as f: - old_content = f.read() - f.seek(0) - namelist_data = f90nml.read(f).todict() - sec = namelist_data.get(section, {}) - sec["self_shielding"] = val - namelist_data[section] = sec +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() - with open(fname, "w") as f: - new_nml = f90nml.Namelist(namelist_data) - new_nml.write(f) + # Copy content of `root` in `tmpdir` recursively + copytree(root, tmpdir) - yield + fname = Path(tmpdir) / "namelist.txt" - with open(fname, "w") as f: - f.write(old_content) + 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) @@ -720,21 +728,19 @@ def test_self_shielding_logic(): ################################################## # Modify the namelist in-situ, reload and check - namelist = os.path.join(ds.root_folder, "namelist.txt") - # Can be set in any of these two for section in ("physics_params", "cooling_params"): for val in (True, False): - with override_self_shielding(namelist, section, val): + with override_self_shielding(ds.root_folder, section, val) as tmp_output: # Autodetection should work - ds = yt.load(ramses_new_format) + 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(ramses_new_format, self_shielding=True) + ds = yt.load(tmp_output, self_shielding=True) assert ds.self_shielding is True - ds = yt.load(ramses_new_format, self_shielding=False) + ds = yt.load(tmp_output, self_shielding=False) assert ds.self_shielding is False ################################################## From d8b87a7354ea4acf3cf05faa422288cabef5baa3 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Sat, 5 Oct 2024 11:42:17 +0200 Subject: [PATCH 10/11] Allow same folder to be created twice Co-authored-by: Chris Havlin --- yt/frontends/ramses/tests/test_outputs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/frontends/ramses/tests/test_outputs.py b/yt/frontends/ramses/tests/test_outputs.py index f2a552d750a..aabd04ae8dd 100644 --- a/yt/frontends/ramses/tests/test_outputs.py +++ b/yt/frontends/ramses/tests/test_outputs.py @@ -699,7 +699,7 @@ def override_self_shielding(root: Path, section: str, val: bool): tmpdir.mkdir() # Copy content of `root` in `tmpdir` recursively - copytree(root, tmpdir) + copytree(root, tmpdir, dirs_exist_ok=True) fname = Path(tmpdir) / "namelist.txt" From 5976261ad345f8760c90fea071a3cf11d7b138de Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 9 Oct 2024 14:40:49 +0200 Subject: [PATCH 11/11] Get root folder from original ds --- yt/frontends/ramses/tests/test_outputs.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/yt/frontends/ramses/tests/test_outputs.py b/yt/frontends/ramses/tests/test_outputs.py index aabd04ae8dd..a68724b42c9 100644 --- a/yt/frontends/ramses/tests/test_outputs.py +++ b/yt/frontends/ramses/tests/test_outputs.py @@ -728,9 +728,10 @@ def test_self_shielding_logic(): ################################################## # 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(ds.root_folder, section, val) as tmp_output: + 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