Merge pull request #5159 from cphyc/bugfix/order-of-operation-matters

[BUGFIX] RAMSES: Use one registry per dataset

Author: chrishavlin

yt/frontends/ramses/field_handlers.py

@@ -42,6 +42,7 @@ def setup_handler(self, domain):
         self.ds = ds = domain.ds
         self.domain = domain
         self.domain_id = domain.domain_id
+
         basename = os.path.abspath(ds.root_folder)
         iout = int(os.path.basename(ds.parameter_filename).split(".")[0].split("_")[1])
 
@@ -166,6 +167,9 @@ def __init_subclass__(cls, *args, **kwargs):
             register_field_handler(cls)
 
         cls._unique_registry = {}
+        cls.parameters = {}
+        cls.rt_parameters = {}
+        cls._detected_field_list = {}
         return cls
 
     def __init__(self, domain):
@@ -231,6 +235,10 @@ def level_count(self):
 
         return self._level_count
 
+    @property
+    def field_list(self):
+        return self._detected_field_list[self.ds.unique_identifier]
+
     @cached_property
     def offset(self):
         """
@@ -242,7 +250,7 @@ def offset(self):
         It should be generic enough for most of the cases, but if the
         *structure* of your fluid file is non-canonical, change this.
         """
-        nvars = len(self.field_list)
+        nvars = len(self._detected_field_list[self.ds.unique_identifier])
         with FortranFile(self.fname) as fd:
             # Skip headers
             nskip = len(self.attrs)
@@ -265,7 +273,7 @@ def offset(self):
                 fd,
                 min_level,
                 self.domain.domain_id,
-                self.parameters["nvar"],
+                self.parameters[self.ds.unique_identifier]["nvar"],
                 self.domain.amr_header,
                 Nskip=nvars * 8,
             )
@@ -314,7 +322,7 @@ def detect_fields(cls, ds):
         attrs = cls.attrs
         with FortranFile(fname) as fd:
             hvals = fd.read_attrs(attrs)
-        cls.parameters = hvals
+        cls.parameters[ds.unique_identifier] = hvals
 
         # Store some metadata
         ds.gamma = hvals["gamma"]
@@ -445,7 +453,9 @@ def detect_fields(cls, ds):
             count_extra += 1
         if count_extra > 0:
             mylog.debug("Detected %s extra fluid fields.", count_extra)
-        cls.field_list = [(cls.ftype, e) for e in fields]
+        cls._detected_field_list[ds.unique_identifier] = [
+            (cls.ftype, e) for e in fields
+        ]
 
         cls.set_detected_fields(ds, fields)
 
@@ -476,9 +486,9 @@ def detect_fields(cls, ds):
         basedir = os.path.split(ds.parameter_filename)[0]
         fname = os.path.join(basedir, cls.fname.format(iout=iout, icpu=1))
         with FortranFile(fname) as fd:
-            cls.parameters = fd.read_attrs(cls.attrs)
+            cls.parameters[ds.unique_identifier] = fd.read_attrs(cls.attrs)
 
-        nvar = cls.parameters["nvar"]
+        nvar = cls.parameters[ds.unique_identifier]["nvar"]
         ndim = ds.dimensionality
 
         fields = cls.load_fields_from_yt_config()
@@ -497,7 +507,9 @@ def detect_fields(cls, ds):
                 for i in range(nvar - ndetected):
                     fields.append(f"var{i}")
 
-        cls.field_list = [(cls.ftype, e) for e in fields]
+        cls._detected_field_list[ds.unique_identifier] = [
+            (cls.ftype, e) for e in fields
+        ]
 
         cls.set_detected_fields(ds, fields)
 
@@ -572,7 +584,7 @@ def read_rhs(cast):
             # Touchy part, we have to read the photon group properties
             mylog.debug("Not reading photon group properties")
 
-            cls.rt_parameters = rheader
+            cls.rt_parameters[ds.unique_identifier] = rheader
 
         ngroups = rheader["nGroups"]
 
@@ -581,7 +593,7 @@ def read_rhs(cast):
         fname = os.path.join(basedir, cls.fname.format(iout=iout, icpu=1))
         fname_desc = os.path.join(basedir, cls.file_descriptor)
         with FortranFile(fname) as fd:
-            cls.parameters = fd.read_attrs(cls.attrs)
+            cls.parameters[ds.unique_identifier] = fd.read_attrs(cls.attrs)
 
         ok = False
 
@@ -615,16 +627,18 @@ def read_rhs(cast):
             for ng in range(ngroups):
                 fields.extend([t % (ng + 1) for t in tmp])
 
-        cls.field_list = [(cls.ftype, e) for e in fields]
+        cls._detected_field_list[ds.unique_identifier] = [
+            (cls.ftype, e) for e in fields
+        ]
 
         cls.set_detected_fields(ds, fields)
         return fields
 
     @classmethod
     def get_rt_parameters(cls, ds):
-        if cls.rt_parameters:
-            return cls.rt_parameters
+        if cls.rt_parameters[ds.unique_identifier]:
+            return cls.rt_parameters[ds.unique_identifier]
 
         # Call detect fields to get the rt_parameters
         cls.detect_fields(ds)
-        return cls.rt_parameters
+        return cls.rt_parameters[ds.unique_identifier]

yt/frontends/ramses/tests/test_outputs.py

@@ -390,8 +390,8 @@ def test_ramses_field_detection():
     fields_1 = set(DETECTED_FIELDS[ds1.unique_identifier]["ramses"])
 
     # Check the right number of variables has been loaded
-    assert P1["nvar"] == 10
-    assert len(fields_1) == P1["nvar"]
+    assert P1[ds1.unique_identifier]["nvar"] == 10
+    assert len(fields_1) == P1[ds1.unique_identifier]["nvar"]
 
     # Now load another dataset
     ds2 = yt.load(output_00080)
@@ -400,8 +400,8 @@ def test_ramses_field_detection():
     fields_2 = set(DETECTED_FIELDS[ds2.unique_identifier]["ramses"])
 
     # Check the right number of variables has been loaded
-    assert P2["nvar"] == 6
-    assert len(fields_2) == P2["nvar"]
+    assert P2[ds2.unique_identifier]["nvar"] == 6
+    assert len(fields_2) == P2[ds2.unique_identifier]["nvar"]
 
 
 @requires_file(ramses_new_format)
@@ -794,3 +794,21 @@ def test_self_shielding_loading():
 
     # Also make sure the difference is large for some cells
     assert (np.abs(diff) > 0.1).any()
+
+
+@requires_file(output_00080)
+@requires_file(ramses_mhd_128)
+def test_order_does_not_matter():
+    for order in (1, 2):
+        ds0 = yt.load(output_00080)
+        ds1 = yt.load(ramses_mhd_128)
+
+        # This should not raise any exception
+        if order == 1:
+            _sp1 = ds1.all_data()
+            sp0 = ds0.all_data()
+        else:
+            sp0 = ds0.all_data()
+            _sp1 = ds1.all_data()
+
+        sp0["gas", "velocity_x"].max().to("km/s")