From eacc9b75d00e5abd49d11c868acf8299ed9930e3 Mon Sep 17 00:00:00 2001
From: Rongjun Huang <u6569836@anu.edu.au>
Date: Wed, 20 Nov 2024 18:24:48 +1100
Subject: [PATCH 01/26] Add Quokka dataset and hierarchy support with metadata
 and header parsing (#3)

Updated the following classes to read header file and `metadata.yaml` in
Quokka simulations.

```
class QuokkaHierarchy(BoxlibHierarchy)
class QuokkaDataset(AMReXDataset)
```

---------

Co-authored-by: Rongjun-ANU <Rongjun.Huang@alumni.anu.edu.au>
---
 yt/frontends/amrex/data_structures.py | 157 ++++++++++++++++++++++++++
 1 file changed, 157 insertions(+)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index ddf0052e812..44199e11bcc 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1,3 +1,4 @@
+import yaml
 import glob
 import os
 import re
@@ -1327,11 +1328,167 @@ def _set_code_unit_attributes(self):
         setdefaultattr(self, "velocity_unit", self.length_unit / self.time_unit)
 
 
+class QuokkaHierarchy(BoxlibHierarchy):
+    def __init__(self, ds, dataset_type="quokka_native"):
+        super().__init__(ds, dataset_type)
+
+        if "particles" in self.ds.parameters:
+            # extra beyond the base real fields that all Boxlib
+            # particles have, i.e. the xyz positions
+            quokka_extra_real_fields = [
+                "particle_velocity_x",
+                "particle_velocity_y",
+                "particle_velocity_z",
+            ]
+
+            is_checkpoint = True
+
+            self._read_particles(
+                "What", # RJ: Quokka uses What for particles, Tracer for Castro, DM for Nyx
+                is_checkpoint,
+                quokka_extra_real_fields[0 : self.ds.dimensionality],
+            )
+
+
 class QuokkaDataset(AMReXDataset):
     # match any plotfiles that have a metadata.yaml file in the root
     _subtype_keyword = ""
     _default_cparam_filename = "metadata.yaml"
 
+    def __init__(
+        self,
+        output_dir,
+        cparam_filename=None,
+        fparam_filename=None,
+        dataset_type="boxlib_native",
+        storage_filename=None,
+        units_override=None,
+        unit_system="cgs",
+        default_species_fields=None,
+    ):
+        # Initialize cparam_filename to the default if not provided
+        if cparam_filename is None:
+            cparam_filename = self._default_cparam_filename
+
+        super().__init__(
+            output_dir,
+            cparam_filename,
+            fparam_filename,
+            dataset_type,
+            storage_filename,
+            units_override,
+            unit_system,
+            default_species_fields=default_species_fields,
+        )
+
+        # Parse the metadata from metadata.yaml after initialization
+        self._parse_metadata_file()
+
+    def _parse_parameter_file(self):
+        # Call parent method to initialize core setup by yt
+        super()._parse_parameter_file()
+
+        # Define the path to the Header file
+        header_filename = os.path.join(self.output_dir, "Header")
+        if not os.path.exists(header_filename):
+            raise FileNotFoundError(f"Header file not found: {header_filename}")
+
+        with open(header_filename, 'r') as f:
+            # Parse header version
+            self.parameters['header_version'] = f.readline().strip()
+
+            # Number of fields
+            num_fields = int(f.readline().strip())
+            self.parameters['fields'] = [f.readline().strip() for _ in range(num_fields)]
+
+            # Dimensionality
+            self.parameters['dimensionality'] = int(f.readline().strip())
+
+            # Simulation time
+            self.parameters['current_time'] = float(f.readline().strip())
+
+            # Refinement levels
+            self.parameters['refinement_level'] = int(f.readline().strip())
+
+            # Domain edges
+            self.parameters['domain_left_edge'] = list(map(float, f.readline().strip().split()))
+            self.parameters['domain_right_edge'] = list(map(float, f.readline().strip().split()))
+
+            # skip empty line
+            f.readline()
+
+            # Grid info
+            self.parameters['grid_info'] = f.readline().strip()
+
+            # Timestamp
+            self.parameters['timestamp'] = list(map(int, f.readline().strip().split()))
+
+            # Grid sizes for all refinement levels
+            self.parameters['grid_sizes'] = []
+            for _ in range(self.parameters['refinement_level'] + 1):
+                self.parameters['grid_sizes'].append(list(map(float, f.readline().strip().split())))
+
+            # Skip placeholders
+            f.readline()  # Placeholder line 1
+            f.readline()  # Placeholder line 2
+
+            # Parse data for all refinement levels
+            self.parameters['refinement_levels'] = []
+
+            for level in range(self.parameters['refinement_level'] + 1):
+                level_data = {}
+
+                # Parse metadata line
+                metadata = f.readline().strip().split()
+                level_data["level"] = int(metadata[0])
+                level_data["num_groups"] = int(metadata[1])
+                level_data["current_time"] = float(metadata[2])
+
+                # Parse timestamp
+                level_data["timestamp"] = int(f.readline().strip())
+
+                # Parse group information
+                level_data["groups"] = []
+                for _ in range(level_data["num_groups"]):
+                    group = {}
+                    for axis in range(self.parameters['dimensionality']):
+                        left_edge, right_edge = map(float, f.readline().strip().split())
+                        group[f"axis_{axis}"] = {"left_edge": left_edge, "right_edge": right_edge}
+                    level_data["groups"].append(group)
+
+                # Append level data to refinement levels
+                self.parameters['refinement_levels'].append(level_data)
+
+                # Skip level marker (e.g., Level_0/Cell)
+                marker = f.readline().strip()
+                if not marker.startswith(f"Level_{level}/Cell"):
+                    raise ValueError(f"Unexpected marker '{marker}' for level {level}")
+
+
+        # hydro method is set by the base class -- override it here
+        self.parameters["HydroMethod"] = "Quokka"
+
+        # Print parsed information for verification
+        print("Parsed header parameters:", self.parameters)
+
+    def _parse_metadata_file(self):
+        # Construct the full path to the metadata file
+        metadata_filename = os.path.join(self.output_dir, self.cparam_filename)
+        try:
+            with open(metadata_filename, 'r') as f:
+                # Load metadata using yaml
+                metadata = yaml.safe_load(f)
+                # Update dataset parameters with metadata if it exists
+                if metadata:
+                    self.parameters.update(metadata)
+                    print("Metadata loaded successfully:", metadata)
+                else:
+                    print("Warning: Metadata file is empty.")
+        except FileNotFoundError:
+            print(f"Error: Metadata file '{metadata_filename}' not found.")
+        except yaml.YAMLError as e:
+            print(f"Error parsing metadata file: {e}")
+
 
 def _guess_pcast(vals):
     # Now we guess some things about the parameter and its type

From cd0dab6433fce453a7d10b08d619b194c88772cf Mon Sep 17 00:00:00 2001
From: Rongjun Huang <u6569836@anu.edu.au>
Date: Tue, 10 Dec 2024 16:36:20 +1100
Subject: [PATCH 02/26] Rongjun update support for temperature, scalars, rad,
 velocities, BFields fields (#5)

This pull request introduces support for the Quokka simulation framework
by adding new dataset and hierarchy classes, and updating field
information. The most important changes include the addition of
`QuokkaDataset` and `QuokkaHierarchy` classes, the parsing of metadata
and parameter files, and the setup of fluid and radiation fields.

### Support for Quokka Framework:

* Added `QuokkaDataset` and `QuokkaHierarchy` classes to handle
Quokka-specific data structures and hierarchy.
(`yt/frontends/amrex/data_structures.py`,
[yt/frontends/amrex/data_structures.pyR1332-R1571](diffhunk://#diff-7aa59ffaeb7262d9580ae294f19155633ce68e2dda607cc4e2ddbf1bc20b66b3R1332-R1571))
* Implemented `_parse_parameter_file` and `_parse_metadata_file` methods
in `QuokkaDataset` to read and parse metadata from `metadata.yaml` and
other parameter files. (`yt/frontends/amrex/data_structures.py`,
[yt/frontends/amrex/data_structures.pyR1332-R1571](diffhunk://#diff-7aa59ffaeb7262d9580ae294f19155633ce68e2dda607cc4e2ddbf1bc20b66b3R1332-R1571))

### Field Information Updates:

* Added `QuokkaFieldInfo` class to define known fields and set up fluid
and radiation fields dynamically. (`yt/frontends/amrex/fields.py`,
[yt/frontends/amrex/fields.pyR506-R597](diffhunk://#diff-565dd12ba00274995d4c360bd7cb3f7d43068a0f2c7f93fe2d33978c1e90c98fR506-R597))

### API Updates:

* Included `QuokkaDataset`, `QuokkaHierarchy`, and `QuokkaFieldInfo` in
the API imports to make them accessible. (`yt/frontends/amrex/api.py`,
[[1]](diffhunk://#diff-92252639991f47b99474012776b030614e644789c56c7a9c39661ee5a301f87bR14-R15)
[[2]](diffhunk://#diff-92252639991f47b99474012776b030614e644789c56c7a9c39661ee5a301f87bR24)

### Miscellaneous:

* Added `yaml` import for parsing metadata files.
(`yt/frontends/amrex/data_structures.py`,
[yt/frontends/amrex/data_structures.pyR1](diffhunk://#diff-7aa59ffaeb7262d9580ae294f19155633ce68e2dda607cc4e2ddbf1bc20b66b3R1))

---------

Co-authored-by: Rongjun-ANU <Rongjun.Huang@alumni.anu.edu.au>
---
 yt/frontends/amrex/api.py             |   3 +
 yt/frontends/amrex/data_structures.py | 152 +++++++++++++++++++++-----
 yt/frontends/amrex/fields.py          | 108 ++++++++++++++++++
 3 files changed, 234 insertions(+), 29 deletions(-)

diff --git a/yt/frontends/amrex/api.py b/yt/frontends/amrex/api.py
index 334cc8bb63f..e7e3727188e 100644
--- a/yt/frontends/amrex/api.py
+++ b/yt/frontends/amrex/api.py
@@ -11,6 +11,8 @@
     NyxHierarchy,
     OrionDataset,
     OrionHierarchy,
+    QuokkaDataset,
+    QuokkaHierarchy,
     WarpXDataset,
     WarpXHierarchy,
 )
@@ -19,6 +21,7 @@
     CastroFieldInfo,
     MaestroFieldInfo,
     NyxFieldInfo,
+    QuokkaFieldInfo,    
     WarpXFieldInfo,
 )
 from .io import IOHandlerBoxlib
diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 44199e11bcc..2a91785becb 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -23,6 +23,7 @@
     CastroFieldInfo,
     MaestroFieldInfo,
     NyxFieldInfo,
+    QuokkaFieldInfo,
     WarpXFieldInfo,
 )
 
@@ -1332,9 +1333,17 @@ class QuokkaHierarchy(BoxlibHierarchy):
     def __init__(self, ds, dataset_type="quokka_native"):
         super().__init__(ds, dataset_type)
 
-        if "particles" in self.ds.parameters:
-            # extra beyond the base real fields that all Boxlib
-            # particles have, i.e. the xyz positions
+        possible_particle_types = ["Rad_particles", "CIC_particles", "tracer_particles"]
+        found_particle_type = None
+
+        # Iterate through possible particle types and check if they exist
+        for particle_type in possible_particle_types:
+            if particle_type in self.ds.parameters.get("particles", []):
+                found_particle_type = particle_type
+                break
+
+        if found_particle_type:
+            # Extra fields beyond base real fields (e.g., xyz positions)
             quokka_extra_real_fields = [
                 "particle_velocity_x",
                 "particle_velocity_y",
@@ -1344,14 +1353,21 @@ def __init__(self, ds, dataset_type="quokka_native"):
             is_checkpoint = True
 
             self._read_particles(
-                "What", # RJ: Quokka uses What for particles, Tracer for Castro, DM for Nyx
+                found_particle_type,  # Use the detected particle type
                 is_checkpoint,
                 quokka_extra_real_fields[0 : self.ds.dimensionality],
             )
+        else:
+            mylog.warning(
+                "No recognized particle types found in the dataset. Checked for: %s",
+                ", ".join(possible_particle_types),
+            )
 
 
 class QuokkaDataset(AMReXDataset):
     # match any plotfiles that have a metadata.yaml file in the root
+    _index_class = QuokkaHierarchy
+    _field_info_class = QuokkaFieldInfo
     _subtype_keyword = ""
     _default_cparam_filename = "metadata.yaml"
 
@@ -1384,6 +1400,10 @@ def __init__(
         # Parse the metadata from metadata.yaml after initialization
         self._parse_metadata_file()
 
+        # Add radiation fields in fluid_types only if detected
+        if self.parameters['radiation_field_groups'] > 0: 
+            self.fluid_types += ("rad",)
+
     def _parse_parameter_file(self):
         # Call parent method to initialize core setup by yt
         super()._parse_parameter_file()
@@ -1399,22 +1419,88 @@ def _parse_parameter_file(self):
 
             # Number of fields
             num_fields = int(f.readline().strip())
-            self.parameters['fields'] = [f.readline().strip() for _ in range(num_fields)]
-
-            # Dimensionality
-            self.parameters['dimensionality'] = int(f.readline().strip())
+            self.parameters['fields'] = []
+
+            # Parse fixed gas fields (6 mandatory fields)
+            gas_fields = [
+                "gasDensity",
+                "x-GasMomentum",
+                "y-GasMomentum",
+                "z-GasMomentum",
+                "gasEnergy",
+                "gasInternalEnergy",
+            ]
+            for field in gas_fields:
+                self.parameters['fields'].append(f.readline().strip())
+
+            # Metadata flags
+            temperature_present = False
+            velocity_present = False
+            bfield_present = False
+            scalar_count = 0
+            rad_group_count = 0
+
+            # Check for optional temperature field
+            next_field = f.readline().strip()
+            if next_field == "gasTemperature":
+                self.parameters['fields'].append(next_field)
+                temperature_present = True
+            else:
+                # Rewind the file pointer if temperature is absent
+                f.seek(f.tell() - len(next_field) - 1)
+
+            # Dynamically parse all remaining fields
+            remaining_fields = num_fields - len(self.parameters['fields'])
+
+            while remaining_fields > 0:
+                current_field = f.readline().strip()
+                remaining_fields -= 1
+
+                if current_field.startswith("scalar_"):
+                    self.parameters['fields'].append(current_field)
+                    scalar_count += 1
+                elif current_field.startswith("radEnergy-Group"):
+                    # Parse radEnergy and corresponding flux fields for a group
+                    self.parameters['fields'].append(current_field)  # radEnergy
+                    self.parameters['fields'].append(f.readline().strip())  # x-RadFlux
+                    self.parameters['fields'].append(f.readline().strip())  # y-RadFlux
+                    self.parameters['fields'].append(f.readline().strip())  # z-RadFlux
+                    remaining_fields -= 3  # Account for the flux fields
+                    rad_group_count += 1
+                elif current_field == "x-velocity":
+                    # Parse velocity fields (x, y, z)
+                    self.parameters['fields'].append(current_field)
+                    self.parameters['fields'].append(f.readline().strip())
+                    self.parameters['fields'].append(f.readline().strip())
+                    remaining_fields -= 2
+                    velocity_present = True
+                elif current_field == "x-BField":
+                    # Parse BField fields (x, y, z)
+                    self.parameters['fields'].append(current_field)
+                    self.parameters['fields'].append(f.readline().strip())
+                    self.parameters['fields'].append(f.readline().strip())
+                    remaining_fields -= 2
+                    bfield_present = True
+                else:
+                    raise ValueError(f"Unexpected field: {current_field}")
 
-            # Simulation time
-            self.parameters['current_time'] = float(f.readline().strip())
+            # Add metadata for radiation groups, scalars, and field existence flags
+            self.parameters['radiation_field_groups'] = rad_group_count
+            self.parameters['scalar_field_count'] = scalar_count
+            self.parameters['temperature_field'] = temperature_present
+            self.parameters['velocity_fields'] = velocity_present
+            self.parameters['Bfields'] = bfield_present
 
-            # Refinement levels
-            self.parameters['refinement_level'] = int(f.readline().strip())
+            # Parse remaining metadata
+            self.parameters['dimensionality'] = int(f.readline().strip())  # Dimensionality
+            self.parameters['current_time'] = float(f.readline().strip())  # Simulation time
+            self.parameters['refinement_level'] = int(f.readline().strip())  # Refinement levels
 
             # Domain edges
             self.parameters['domain_left_edge'] = list(map(float, f.readline().strip().split()))
             self.parameters['domain_right_edge'] = list(map(float, f.readline().strip().split()))
 
-            # skip empty line
+            # Skip empty line
             f.readline()
 
             # Grid info
@@ -1433,7 +1519,8 @@ def _parse_parameter_file(self):
             f.readline()  # Placeholder line 2
 
             # Parse data for all refinement levels
-            self.parameters['refinement_levels'] = []
+            self.parameters['refinement_details'] = []
+            max_refinement_level = -1
 
             for level in range(self.parameters['refinement_level'] + 1):
                 level_data = {}
@@ -1441,35 +1528,42 @@ def _parse_parameter_file(self):
                 # Parse metadata line
                 metadata = f.readline().strip().split()
                 level_data["level"] = int(metadata[0])
-                level_data["num_groups"] = int(metadata[1])
+                level_data["num_boxes"] = int(metadata[1])  # Number of boxes in current level
                 level_data["current_time"] = float(metadata[2])
 
+                # Update max refinement level
+                max_refinement_level = max(max_refinement_level, level_data["level"])
+
                 # Parse timestamp
                 level_data["timestamp"] = int(f.readline().strip())
 
-                # Parse group information
-                level_data["groups"] = []
-                for _ in range(level_data["num_groups"]):
-                    group = {}
+                # Parse box information
+                level_data["boxes"] = []
+                for _ in range(level_data["num_boxes"]):
+                    box = {}
                     for axis in range(self.parameters['dimensionality']):
                         left_edge, right_edge = map(float, f.readline().strip().split())
-                        group[f"axis_{axis}"] = {"left_edge": left_edge, "right_edge": right_edge}
-                    level_data["groups"].append(group)
+                        box[f"axis_{axis}"] = {"left_edge": left_edge, "right_edge": right_edge}
+                    level_data["boxes"].append(box)
 
                 # Append level data to refinement levels
-                self.parameters['refinement_levels'].append(level_data)
+                self.parameters['refinement_details'].append(level_data)
 
                 # Skip level marker (e.g., Level_0/Cell)
                 marker = f.readline().strip()
                 if not marker.startswith(f"Level_{level}/Cell"):
                     raise ValueError(f"Unexpected marker '{marker}' for level {level}")
 
+            # Update parameters for refinement
+            self.parameters["max_refinement_level"] = max_refinement_level
+            self.parameters["refinement"] = max_refinement_level > 0
 
-        # hydro method is set by the base class -- override it here
+        # Set hydro method explicitly
         self.parameters["HydroMethod"] = "Quokka"
 
-        # Print parsed information for verification
-        print("Parsed header parameters:", self.parameters)
+        # Debug output for verification
+        mylog.debug("Header loaded successfully")
+        mylog.debug("Parsed header parameters: %s", self.parameters)
 
     def _parse_metadata_file(self):
         # Construct the full path to the metadata file
@@ -1481,13 +1575,13 @@ def _parse_metadata_file(self):
                 # Update dataset parameters with metadata if it exists
                 if metadata:
                     self.parameters.update(metadata)
-                    print("Metadata loaded successfully:", metadata)
+                    mylog.debug("Metadata loaded successfully")
                 else:
-                    print("Warning: Metadata file is empty.")
+                    mylog.debug("Warning: Metadata file is empty.")
         except FileNotFoundError:
-            print(f"Error: Metadata file '{metadata_filename}' not found.")
+            mylog.debug(f"Error: Metadata file '{metadata_filename}' not found.")
         except yaml.YAMLError as e:
-            print(f"Error parsing metadata file: {e}")
+            mylog.debug(f"Error parsing metadata file: {e}")
 
 
 def _guess_pcast(vals):
diff --git a/yt/frontends/amrex/fields.py b/yt/frontends/amrex/fields.py
index 6354b7d65e0..a5e7071f7e4 100644
--- a/yt/frontends/amrex/fields.py
+++ b/yt/frontends/amrex/fields.py
@@ -503,6 +503,114 @@ def setup_fluid_fields(self):
                     units=unit_system["frequency"],
                 )
 
+class QuokkaFieldInfo(FieldInfoContainer):
+    known_other_fields: KnownFieldsT = (
+        ("gasDensity", ("code_mass / code_length**3", ["density"], r"\rho")),
+        ("gasEnergy", ("code_mass / (code_length * code_time**2)", ["total_energy_density"], r"\rho E")),
+        ("gasInternalEnergy", ("code_mass / (code_length * code_time**2)", ["internal_energy_density"], r"\rho e")),
+        ("x-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_x"], r"\rho u")),
+        ("y-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_y"], r"\rho v")),
+        ("z-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_z"], r"\rho w")),
+        ("scalar_0", ("", ["scalar_0"], "Scalar 0")),
+        ("scalar_1", ("", ["scalar_1"], "Scalar 1")),
+        ("scalar_2", ("", ["scalar_2"], "Scalar 2")),
+    )
+
+    known_particle_fields: KnownFieldsT = (
+        ("particle_mass", ("code_mass", [], None)),
+        ("particle_position_x", ("code_length", [], None)),
+        ("particle_position_y", ("code_length", [], None)),
+        ("particle_position_z", ("code_length", [], None)),
+        ("particle_momentum_x", ("code_mass*code_length/code_time", [], None)),
+        ("particle_momentum_y", ("code_mass*code_length/code_time", [], None)),
+        ("particle_momentum_z", ("code_mass*code_length/code_time", [], None)),
+        # Note that these are *internal* agmomen
+        ("particle_angmomen_x", ("code_length**2/code_time", [], None)),
+        ("particle_angmomen_y", ("code_length**2/code_time", [], None)),
+        ("particle_angmomen_z", ("code_length**2/code_time", [], None)),
+        ("particle_id", ("", ["particle_index"], None)),
+        ("particle_mdot", ("code_mass/code_time", [], None)),
+    )
+
+    def setup_fluid_fields(self):
+        # Define momentum-to-velocity conversion
+        def _get_cell_velocity(axis):
+            def velocity(field, data):
+                # Divide momentum by density for cell-centered velocity
+                return data["boxlib", f"{axis}-GasMomentum"] / data["boxlib", "gasDensity"]
+
+            return velocity
+
+        # Add cell-centered velocity fields dynamically for each axis
+        for axis in "xyz":
+            field_name = f"velocity_{axis}"
+            momentum_name = f"{axis}-GasMomentum"
+
+            if ("boxlib", momentum_name) in self.field_list:
+                self.add_field(
+                    ("gas", field_name),
+                    sampling_type="cell",
+                    function=_get_cell_velocity(axis),
+                    units=self.ds.unit_system["length"] / self.ds.unit_system["time"],
+                    display_name=f"v_{axis} (cell-centered)",
+                )
+
+        # Add face-centered velocities dynamically for each axis
+        for axis in "xyz":
+            face_velocity_name = f"velocity_face_{axis}"
+
+            if ("boxlib", f"{axis}-velocity") in self.field_list:
+                self.add_field(
+                    ("gas", face_velocity_name),
+                    sampling_type="cell",
+                    function=lambda field, data, axis=axis: data["boxlib", f"{axis}-velocity"] * self.ds.unit_system["length"] / self.ds.unit_system["time"],
+                    units=self.ds.unit_system["length"] / self.ds.unit_system["time"],
+                    display_name=f"v_{axis} (face-centered)",
+                )
+
+        # Add magnetic fields dynamically for each axis (Placeholder for now)
+        for axis in "xyz":
+            bfield_name = f"BField_{axis}"
+            boxlib_bfield = f"{axis}-BField"
+
+            if ("boxlib", boxlib_bfield) in self.field_list:
+                self.add_field(
+                    ("gas", bfield_name),
+                    sampling_type="cell",
+                    function=lambda field, data, axis=axis: data["boxlib", boxlib_bfield],
+                    units="",  # Leave units empty for now
+                    display_name=f"B_{axis} (magnetic field)",
+                )
+
+        # Call the new radiation fields setup
+        self.setup_radiation_fields()
+
+    def setup_radiation_fields(self):
+        # Dynamically add radiation fields
+        num_groups = self.ds.parameters.get("radiation_field_groups", 0)
+        for group in range(num_groups):
+            # Add energy density
+            energy_field = f"radEnergy-Group{group}"
+            if ("boxlib", energy_field) in self.field_list:
+                self.add_field(
+                    ("rad", f"energy_density_{group}"),
+                    sampling_type="cell",
+                    function=lambda field, data: data["boxlib", energy_field] * self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
+                    units=self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
+                    display_name=f"Radiation Energy Density Group {group}",
+                )
+
+            # Add flux density for each axis
+            for axis in "xyz":
+                flux_field = f"{axis}-RadFlux-Group{group}"
+                if ("boxlib", flux_field) in self.field_list:
+                    self.add_field(
+                        ("rad", f"flux_density_{axis}_{group}"),
+                        sampling_type="cell",
+                        function=lambda field, data, flux_field=flux_field: data["boxlib", flux_field] * self.ds.unit_system["energy"] / self.ds.unit_system["length"]**2 / self.ds.unit_system["time"],
+                        units=self.ds.unit_system["energy"] / self.ds.unit_system["length"]**2 / self.ds.unit_system["time"],
+                        display_name=f"Radiation Flux Density {axis.upper()} Group {group}",
+                    )
 
 substance_expr_re = re.compile(r"\(([a-zA-Z][a-zA-Z0-9]*)\)")
 substance_elements_re = re.compile(r"(?P<element>[a-zA-Z]+)(?P<digits>\d*)")

From f5acf45ef1f936bbab5149678136457d29197ee8 Mon Sep 17 00:00:00 2001
From: Rongjun Huang <u6569836@anu.edu.au>
Date: Tue, 10 Dec 2024 23:51:24 +1100
Subject: [PATCH 03/26] add support of multiple particle types (#6)

slight update on `class QuokkaHierarchy(BoxlibHierarchy)` to fix the
bug, and also add support on multi particle types.

---------

Co-authored-by: Rongjun-ANU <Rongjun.Huang@alumni.anu.edu.au>
---
 yt/frontends/amrex/data_structures.py | 45 ++++++++++++++++-----------
 1 file changed, 26 insertions(+), 19 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 2a91785becb..fb67665762d 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1333,17 +1333,17 @@ class QuokkaHierarchy(BoxlibHierarchy):
     def __init__(self, ds, dataset_type="quokka_native"):
         super().__init__(ds, dataset_type)
 
-        possible_particle_types = ["Rad_particles", "CIC_particles", "tracer_particles"]
-        found_particle_type = None
+        # Dynamically detect particle types by searching for "*_particles" directories
+        particle_dirs = glob.glob(os.path.join(self.ds.output_dir, "*_particles"))
+        detected_particle_types = []
 
-        # Iterate through possible particle types and check if they exist
-        for particle_type in possible_particle_types:
-            if particle_type in self.ds.parameters.get("particles", []):
-                found_particle_type = particle_type
-                break
+        for pdir in particle_dirs:
+            ptype = os.path.basename(pdir)
+            header_file = os.path.join(pdir, "Header")
+            if os.path.exists(header_file):
+                detected_particle_types.append(ptype)
 
-        if found_particle_type:
-            # Extra fields beyond base real fields (e.g., xyz positions)
+        if detected_particle_types:
             quokka_extra_real_fields = [
                 "particle_velocity_x",
                 "particle_velocity_y",
@@ -1352,16 +1352,23 @@ def __init__(self, ds, dataset_type="quokka_native"):
 
             is_checkpoint = True
 
-            self._read_particles(
-                found_particle_type,  # Use the detected particle type
-                is_checkpoint,
-                quokka_extra_real_fields[0 : self.ds.dimensionality],
-            )
+            for ptype in detected_particle_types:
+                self._read_particles(
+                    ptype,  # Use the detected particle type
+                    is_checkpoint,
+                    quokka_extra_real_fields[0 : self.ds.dimensionality],
+                )
+
+            # Update the dataset's particle types
+            self.ds.parameters["particles"] = len(detected_particle_types)
+            self.ds.particle_types = tuple(detected_particle_types)
+            self.ds.particle_types_raw = self.ds.particle_types
         else:
-            mylog.warning(
-                "No recognized particle types found in the dataset. Checked for: %s",
-                ", ".join(possible_particle_types),
-            )
+            # No particle types detected
+            self.ds.parameters["particles"] = 0
+            self.ds.particle_types = ()
+            self.ds.particle_types_raw = ()
+            mylog.debug("No recognized particle types found in the dataset.")
 
 
 class QuokkaDataset(AMReXDataset):
@@ -1849,4 +1856,4 @@ def _set_code_unit_attributes(self):
         setdefaultattr(self, "mass_unit", self.quan(1.0, "kg"))
         setdefaultattr(self, "time_unit", self.quan(1.0, "s"))
         setdefaultattr(self, "velocity_unit", self.quan(1.0, "m/s"))
-        setdefaultattr(self, "magnetic_unit", self.quan(1.0, "T"))
+        setdefaultattr(self, "magnetic_unit", self.quan(1.0, "T"))
\ No newline at end of file

From ae0dd250fcd113016fbad30d17d42aabe6cb65a3 Mon Sep 17 00:00:00 2001
From: Rongjun Huang <u6569836@anu.edu.au>
Date: Tue, 7 Jan 2025 11:23:08 +1100
Subject: [PATCH 04/26] further support for particles and their units (#8)

add support for `Fields.json` in particle folders. also readout the unit of particles.

---------

Co-authored-by: Rongjun-ANU <Rongjun.Huang@alumni.anu.edu.au>
---
 yt/frontends/amrex/data_structures.py | 153 ++++++++++++++++++++------
 yt/frontends/amrex/fields.py          |  76 +++++++++++++
 2 files changed, 198 insertions(+), 31 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 7f0fd4b63b9..2daffb325f1 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1,4 +1,5 @@
 import yaml
+import json
 import glob
 import os
 import re
@@ -1336,44 +1337,39 @@ def __init__(self, ds, dataset_type="quokka_native"):
 
         # Dynamically detect particle types by searching for "*_particles" directories
         particle_dirs = glob.glob(os.path.join(self.ds.output_dir, "*_particles"))
-        detected_particle_types = []
-
+        
         for pdir in particle_dirs:
             ptype = os.path.basename(pdir)
             header_file = os.path.join(pdir, "Header")
-            if os.path.exists(header_file):
-                detected_particle_types.append(ptype)
 
-        if detected_particle_types:
-            quokka_extra_real_fields = [
-                "particle_velocity_x",
-                "particle_velocity_y",
-                "particle_velocity_z",
-            ]
+            if os.path.exists(header_file):
+                # Read the Header to determine the number of fields
+                with open(header_file, 'r') as f:
+                    _ = f.readline().strip()  # Skip version info
+                    _ = f.readline().strip()  # Skip number of particles
+                    num_fields = int(f.readline().strip())  # Number of fields
+
+                    # Read the field names
+                    fields = [f.readline().strip() for _ in range(num_fields)]
+
+                # Define extra real fields specific to Quokka
+                quokka_extra_real_fields = [
+                    "particle_velocity_x",
+                    "particle_velocity_y",
+                    "particle_velocity_z",
+                ]
 
-            is_checkpoint = True
+                is_checkpoint = False
 
-            for ptype in detected_particle_types:
-                self._read_particles(
-                    ptype,  # Use the detected particle type
-                    is_checkpoint,
-                    quokka_extra_real_fields[0 : self.ds.dimensionality],
-                )
+                # Pass the fields to _read_particles
+                extra_fields = quokka_extra_real_fields[:ds.dimensionality] + fields
 
-            # Update the dataset's particle types
-            self.ds.parameters["particles"] = len(detected_particle_types)
-            self.ds.particle_types = tuple(detected_particle_types)
-            self.ds.particle_types_raw = self.ds.particle_types
-        else:
-            # No particle types detected
-            self.ds.parameters["particles"] = 0
-            self.ds.particle_types = ()
-            self.ds.particle_types_raw = ()
-            mylog.debug("No recognized particle types found in the dataset.")
+                # Read particle data for this type
+                self._read_particles(ptype, is_checkpoint, extra_fields)
 
 
 class QuokkaDataset(AMReXDataset):
-    # match any plotfiles that have a metadata.yaml file in the root
+    # Match any plotfiles that have a metadata.yaml file in the root
     _index_class = QuokkaHierarchy
     _field_info_class = QuokkaFieldInfo
     _subtype_keyword = ""
@@ -1390,7 +1386,6 @@ def __init__(
         unit_system="cgs",
         default_species_fields=None,
     ):
-        # Initialize cparam_filename to the default if not provided
         if cparam_filename is None:
             cparam_filename = self._default_cparam_filename
 
@@ -1409,7 +1404,7 @@ def __init__(
         self._parse_metadata_file()
 
         # Add radiation fields in fluid_types only if detected
-        if self.parameters['radiation_field_groups'] > 0: 
+        if self.parameters.get('radiation_field_groups', 0) > 0:
             self.fluid_types += ("rad",)
 
     def _parse_parameter_file(self):
@@ -1423,7 +1418,7 @@ def _parse_parameter_file(self):
 
         with open(header_filename, 'r') as f:
             # Parse header version
-            self.parameters['header_version'] = f.readline().strip()
+            self.parameters['plot_file_type'] = f.readline().strip()
 
             # Number of fields
             num_fields = int(f.readline().strip())
@@ -1566,6 +1561,102 @@ def _parse_parameter_file(self):
             self.parameters["max_refinement_level"] = max_refinement_level
             self.parameters["refinement"] = max_refinement_level > 0
 
+        # Handle particle information from "*_particles/Header"
+        particle_dirs = glob.glob(os.path.join(self.output_dir, "*_particles"))
+        detected_particle_types = []
+        particle_info = {}
+
+        for pdir in particle_dirs:
+            particle_type = os.path.basename(pdir)
+            header_file = os.path.join(pdir, "Header")
+            fields_json_file = os.path.join(pdir, "Fields.json")
+
+            if os.path.exists(header_file):
+                detected_particle_types.append(particle_type)
+
+                # Parse the Header
+                with open(header_file, 'r') as f:
+                    f.readline().strip()  # Skip version line
+                    num_particles = int(f.readline().strip())  # Second line
+                    num_fields = int(f.readline().strip())  # Third line
+                    fields = [f.readline().strip() for _ in range(num_fields)]  # Remaining lines
+
+                # Default field names and units
+                field_names = fields[:]
+                field_units = {field: "dimensionless" for field in fields}
+
+                # Initialize variables
+                json_field_names = None
+
+                # Check and parse Fields.json
+                if os.path.exists(fields_json_file):
+                    try:
+                        with open(fields_json_file, 'r') as f:
+                            field_data = json.load(f)
+                            json_field_names = list(field_data.get("fields", []))
+                            raw_units = field_data.get("units", {})
+
+                            # Validate field names count
+                            if len(json_field_names) == len(fields):
+                                # Replace the default fields (real_comp*) with those from Fields.json
+                                field_names = json_field_names
+                                field_units = {}
+
+                                # Translate raw unit dimensions into readable strings
+                                base_units = ["M", "L", "T", "Θ"]  # Mass, Length, Time, Temperature
+                                for field, unit_dims in raw_units.items():
+                                    field_units[field] = " ".join(
+                                        f"{base_units[i]}^{exp}" if exp != 0 else ""
+                                        for i, exp in enumerate(unit_dims)
+                                    ).strip() or "dimensionless"
+
+                                # Remove any `real_comp*` entries from the field_units
+                                for idx in range(len(fields)):
+                                    real_comp_field = f"real_comp{idx}"
+                                    if real_comp_field in field_units:
+                                        del field_units[real_comp_field]
+
+                            else:
+                                mylog.debug(
+                                    "Field names in Fields.json do not match the number of fields in Header for '%s'. "
+                                    "Using names from Header.",
+                                    particle_type,
+                                )
+                    except Exception as e:
+                        mylog.debug(
+                            "Failed to parse Fields.json for particle type '%s': %s",
+                            particle_type,
+                            e,
+                        )
+                else:
+                    # Ensure all fields have units (handle missing units gracefully)
+                    field_units = {field: "dimensionless" for field in fields}
+
+                # Explicitly remove real_comp* entries if Fields.json replaces the fields
+                if field_names == json_field_names:
+                    for idx in range(len(fields)):
+                        real_comp_field = f"real_comp{idx}"
+                        if real_comp_field in field_units:
+                            del field_units[real_comp_field]
+
+                # Add particle info
+                particle_info[particle_type] = {
+                    "num_particles": num_particles,
+                    "num_fields": num_fields,
+                    "fields": field_names,
+                    "units": field_units,
+                }
+
+        # Update parameters with particle info
+        self.parameters.update({
+            "particles": len(detected_particle_types),
+            "particle_types": tuple(detected_particle_types),
+            "particle_info": particle_info,
+        })
+
+        # Debugging: Log parameters for verification
+        mylog.debug("Updated ds.parameters with particle info: %s", self.parameters)
+
         # Set hydro method explicitly
         self.parameters["HydroMethod"] = "Quokka"
 
diff --git a/yt/frontends/amrex/fields.py b/yt/frontends/amrex/fields.py
index a5e7071f7e4..87d8b1dc8b2 100644
--- a/yt/frontends/amrex/fields.py
+++ b/yt/frontends/amrex/fields.py
@@ -611,6 +611,82 @@ def setup_radiation_fields(self):
                         units=self.ds.unit_system["energy"] / self.ds.unit_system["length"]**2 / self.ds.unit_system["time"],
                         display_name=f"Radiation Flux Density {axis.upper()} Group {group}",
                     )
+        
+        # Dynamically set up custom particle fields (convert `real_comp` to physics) for all particle types
+        for ptype in self.ds.particle_types:
+            self.setup_custom_particle_fields(ptype)
+
+    def setup_custom_particle_fields(self, ptype):
+        """
+        Dynamically set up custom particle fields (real_comp) for the given particle type, 
+        interpreting dimensional expressions and applying appropriate unit conversions.
+
+        Parameters:
+        -----------
+        ptype : str
+            The particle type (e.g., 'Rad_particles', 'CIC_particles', etc.) for which fields will be dynamically added.
+        """
+        particle_info = self.ds.parameters.get("particle_info", {}).get(ptype, {})
+        field_map = particle_info.get("fields", [])
+        unit_map = particle_info.get("units", {})
+
+        def interpret_units(dim_expr):
+            """
+            Interpret a dimensional expression like 'M^a L^b T^c Θ^d' and return the corresponding unit factor.
+
+            Parameters:
+            -----------
+            dim_expr : str
+                Dimensional expression (e.g., 'M^1 L^2 T^-3' for luminosity).
+
+            Returns:
+            --------
+            conversion_factor : unyt_quantity
+                The unit conversion factor constructed from the dimensional expression.
+            """
+            if not dim_expr or dim_expr == "dimensionless":
+                return 1.0  # No conversion needed for dimensionless fields
+
+            # Parse the dimensional expression
+            dimensions = {"M": 0, "L": 0, "T": 0, "Θ": 0}  # Default to zero exponents
+            for term in dim_expr.split():
+                base, exponent = term.split("^")
+                dimensions[base] = int(exponent)
+
+            # Construct the conversion factor using the unit system. Use code_mass, code_length, code_time format for units
+            conversion_factor = (
+                f"code_mass**{dimensions['M']} * "
+                f"code_length**{dimensions['L']} * "
+                f"code_time**{dimensions['T']} * "
+                f"code_temperature**{dimensions['Θ']}"
+            )
+            # Convert the string expression to actual unit quantity
+            conversion_factor = self.ds.quan(1.0, conversion_factor)
+            return conversion_factor
+
+        for idx, field_name in enumerate(field_map):
+            # Construct the `real_comp` field name
+            real_comp_field = f"particle_real_comp{idx}"
+
+            # Check if the `real_comp` field exists in the dataset
+            if (ptype, real_comp_field) in self.field_list:
+                # Retrieve the dimensional expression for the new field
+                dim_expr = unit_map.get(field_name, "dimensionless")
+
+                # Interpret the dimensional expression to get the conversion factor
+                conversion_factor = interpret_units(dim_expr)
+
+                # Add the field with the appropriate units
+                self.add_field(
+                    (ptype, field_name),
+                    sampling_type="particle",
+                    function=lambda field, data, real_comp_field=real_comp_field, conv=conversion_factor: (
+                        data[ptype, real_comp_field] * conv
+                    ),
+                    units=conversion_factor.units if hasattr(conversion_factor, "units") else "",
+                    display_name=field_name.replace("_", " ").capitalize(),
+                )
+
 
 substance_expr_re = re.compile(r"\(([a-zA-Z][a-zA-Z0-9]*)\)")
 substance_elements_re = re.compile(r"(?P<element>[a-zA-Z]+)(?P<digits>\d*)")

From cdb199775e5e9f70fb3fe944b850b8b438ab53eb Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Tue, 14 Jan 2025 16:30:19 +1100
Subject: [PATCH 05/26] clean up

---
 yt/frontends/amrex/data_structures.py | 47 +++++----------------------
 1 file changed, 8 insertions(+), 39 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 2daffb325f1..f7a0064ef7a 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1440,58 +1440,27 @@ def _parse_parameter_file(self):
             temperature_present = False
             velocity_present = False
             bfield_present = False
-            scalar_count = 0
             rad_group_count = 0
 
-            # Check for optional temperature field
-            next_field = f.readline().strip()
-            if next_field == "gasTemperature":
-                self.parameters['fields'].append(next_field)
-                temperature_present = True
-            else:
-                # Rewind the file pointer if temperature is absent
-                f.seek(f.tell() - len(next_field) - 1)
-
             # Dynamically parse all remaining fields
-            remaining_fields = num_fields - len(self.parameters['fields'])
+            remaining_fields = num_fields - len(gas_fields)
+
+            for i in range(remaining_fields):
 
-            while remaining_fields > 0:
                 current_field = f.readline().strip()
-                remaining_fields -= 1
-
-                if current_field.startswith("scalar_"):
-                    self.parameters['fields'].append(current_field)
-                    scalar_count += 1
-                elif current_field.startswith("radEnergy-Group"):
-                    # Parse radEnergy and corresponding flux fields for a group
-                    self.parameters['fields'].append(current_field)  # radEnergy
-                    self.parameters['fields'].append(f.readline().strip())  # x-RadFlux
-                    self.parameters['fields'].append(f.readline().strip())  # y-RadFlux
-                    self.parameters['fields'].append(f.readline().strip())  # z-RadFlux
-                    remaining_fields -= 3  # Account for the flux fields
+                if current_field.startswith("radEnergy-Group"):
                     rad_group_count += 1
                 elif current_field == "x-velocity":
-                    # Parse velocity fields (x, y, z)
-                    self.parameters['fields'].append(current_field)
-                    self.parameters['fields'].append(f.readline().strip())
-                    self.parameters['fields'].append(f.readline().strip())
-                    remaining_fields -= 2
                     velocity_present = True
                 elif current_field == "x-BField":
-                    # Parse BField fields (x, y, z)
-                    self.parameters['fields'].append(current_field)
-                    self.parameters['fields'].append(f.readline().strip())
-                    self.parameters['fields'].append(f.readline().strip())
-                    remaining_fields -= 2
                     bfield_present = True
-                else:
-                    raise ValueError(f"Unexpected field: {current_field}")
+                elif current_field == "gasTemperature":
+                    temperature_present = True
+
+                self.parameters['fields'].append(current_field)
 
             # Add metadata for radiation groups, scalars, and field existence flags
             self.parameters['radiation_field_groups'] = rad_group_count
-            self.parameters['scalar_field_count'] = scalar_count
-            self.parameters['temperature_field'] = temperature_present
-            self.parameters['velocity_fields'] = velocity_present
             self.parameters['Bfields'] = bfield_present
 
             # Parse remaining metadata

From 6b455426a959fa967a70fb254759f2f9c05f38c1 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Tue, 14 Jan 2025 16:40:12 +1100
Subject: [PATCH 06/26] more cleanup

---
 yt/frontends/amrex/data_structures.py | 26 ++------------------------
 1 file changed, 2 insertions(+), 24 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index f7a0064ef7a..0cf9e2cd54c 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1420,43 +1420,21 @@ def _parse_parameter_file(self):
             # Parse header version
             self.parameters['plot_file_type'] = f.readline().strip()
 
+            # expecting mandatory fields: gasDensity, x-GasMomentum, y-GasMomentum, z-GasMomentum, gasEnergy, gasInternalEnergy
             # Number of fields
             num_fields = int(f.readline().strip())
             self.parameters['fields'] = []
 
-            # Parse fixed gas fields (6 mandatory fields)
-            gas_fields = [
-                "gasDensity",
-                "x-GasMomentum",
-                "y-GasMomentum",
-                "z-GasMomentum",
-                "gasEnergy",
-                "gasInternalEnergy",
-            ]
-            for field in gas_fields:
-                self.parameters['fields'].append(f.readline().strip())
-
             # Metadata flags
-            temperature_present = False
-            velocity_present = False
             bfield_present = False
             rad_group_count = 0
 
-            # Dynamically parse all remaining fields
-            remaining_fields = num_fields - len(gas_fields)
-
-            for i in range(remaining_fields):
-
+            for i in range(num_fields):
                 current_field = f.readline().strip()
                 if current_field.startswith("radEnergy-Group"):
                     rad_group_count += 1
-                elif current_field == "x-velocity":
-                    velocity_present = True
                 elif current_field == "x-BField":
                     bfield_present = True
-                elif current_field == "gasTemperature":
-                    temperature_present = True
-
                 self.parameters['fields'].append(current_field)
 
             # Add metadata for radiation groups, scalars, and field existence flags

From c7ad3e60f65e62b3ba689a7156f43f9b4dc5db12 Mon Sep 17 00:00:00 2001
From: Rongjun Huang <u6569836@anu.edu.au>
Date: Wed, 15 Jan 2025 01:34:57 +1100
Subject: [PATCH 07/26] Rongjun read quokka metadata header (#10)

update for Bfields

---------

Co-authored-by: Rongjun-ANU <Rongjun.Huang@alumni.anu.edu.au>
---
 yt/frontends/amrex/data_structures.py |  3 ---
 yt/frontends/amrex/fields.py          | 25 ++++++++++++++++++-------
 2 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 0cf9e2cd54c..a83450443a0 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1433,13 +1433,10 @@ def _parse_parameter_file(self):
                 current_field = f.readline().strip()
                 if current_field.startswith("radEnergy-Group"):
                     rad_group_count += 1
-                elif current_field == "x-BField":
-                    bfield_present = True
                 self.parameters['fields'].append(current_field)
 
             # Add metadata for radiation groups, scalars, and field existence flags
             self.parameters['radiation_field_groups'] = rad_group_count
-            self.parameters['Bfields'] = bfield_present
 
             # Parse remaining metadata
             self.parameters['dimensionality'] = int(f.readline().strip())  # Dimensionality
diff --git a/yt/frontends/amrex/fields.py b/yt/frontends/amrex/fields.py
index 87d8b1dc8b2..5fedeb2844b 100644
--- a/yt/frontends/amrex/fields.py
+++ b/yt/frontends/amrex/fields.py
@@ -511,6 +511,9 @@ class QuokkaFieldInfo(FieldInfoContainer):
         ("x-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_x"], r"\rho u")),
         ("y-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_y"], r"\rho v")),
         ("z-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_z"], r"\rho w")),
+        # Temperature field is not present in early Quokka datasets
+        ("gasTemperature", ("K", ["temperature"], r"T")), 
+        # Scalar fields are not always present
         ("scalar_0", ("", ["scalar_0"], "Scalar 0")),
         ("scalar_1", ("", ["scalar_1"], "Scalar 1")),
         ("scalar_2", ("", ["scalar_2"], "Scalar 2")),
@@ -568,23 +571,31 @@ def velocity(field, data):
                     display_name=f"v_{axis} (face-centered)",
                 )
 
-        # Add magnetic fields dynamically for each axis (Placeholder for now)
+        # Call the Bfields and radiation fields setup
+        self.setup_Bfields() # magnetic fields are still placeholder here
+        self.setup_radiation_fields()
+
+    def setup_Bfields(self):
+        """
+        Dynamically add magnetic fields based on presence of Bfield fields in ds.parameters['fields']
+        """
+        # Check if any field name contains 'Bfield'
+        if not any('Bfield' in field for field in self.ds.parameters.get('fields', [])):
+            return
+
         for axis in "xyz":
             bfield_name = f"BField_{axis}"
             boxlib_bfield = f"{axis}-BField"
 
             if ("boxlib", boxlib_bfield) in self.field_list:
                 self.add_field(
-                    ("gas", bfield_name),
+                    ("mag", f"{axis}-field"),
                     sampling_type="cell",
-                    function=lambda field, data, axis=axis: data["boxlib", boxlib_bfield],
-                    units="",  # Leave units empty for now
+                    function=lambda field, data, axis=axis: data["boxlib", f"{axis}-BField"] * self.ds.unit_system["magnetic_field_strength"],
+                    units=self.ds.unit_system["magnetic_field_strength"],
                     display_name=f"B_{axis} (magnetic field)",
                 )
 
-        # Call the new radiation fields setup
-        self.setup_radiation_fields()
-
     def setup_radiation_fields(self):
         # Dynamically add radiation fields
         num_groups = self.ds.parameters.get("radiation_field_groups", 0)

From 0fefa614e932579140283f346f8ea59bea3e2762 Mon Sep 17 00:00:00 2001
From: Rongjun Huang <u6569836@anu.edu.au>
Date: Thu, 20 Mar 2025 15:10:59 +1100
Subject: [PATCH 08/26] rename Fields.json to Fields.yaml

rename Fields.json to Fields.yaml

---------

Co-authored-by: Rongjun-ANU <Rongjun.Huang@alumni.anu.edu.au>
---
 yt/frontends/amrex/data_structures.py | 35 ++++++++++++++-------------
 1 file changed, 18 insertions(+), 17 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index a83450443a0..d72f66207bc 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1,5 +1,4 @@
 import yaml
-import json
 import glob
 import os
 import re
@@ -1406,6 +1405,9 @@ def __init__(
         # Add radiation fields in fluid_types only if detected
         if self.parameters.get('radiation_field_groups', 0) > 0:
             self.fluid_types += ("rad",)
+        # Add magnetic fields in fluid_types only if detected
+        if "Bfield" in self.parameters['fields']:
+            self.fluid_types += ("mag",)
 
     def _parse_parameter_file(self):
         # Call parent method to initialize core setup by yt
@@ -1426,7 +1428,6 @@ def _parse_parameter_file(self):
             self.parameters['fields'] = []
 
             # Metadata flags
-            bfield_present = False
             rad_group_count = 0
 
             for i in range(num_fields):
@@ -1513,7 +1514,7 @@ def _parse_parameter_file(self):
         for pdir in particle_dirs:
             particle_type = os.path.basename(pdir)
             header_file = os.path.join(pdir, "Header")
-            fields_json_file = os.path.join(pdir, "Fields.json")
+            fields_yaml_file = os.path.join(pdir, "Fields.yaml")
 
             if os.path.exists(header_file):
                 detected_particle_types.append(particle_type)
@@ -1530,20 +1531,20 @@ def _parse_parameter_file(self):
                 field_units = {field: "dimensionless" for field in fields}
 
                 # Initialize variables
-                json_field_names = None
+                yaml_field_names = None
 
-                # Check and parse Fields.json
-                if os.path.exists(fields_json_file):
+                # Check and parse Fields.yaml
+                if os.path.exists(fields_yaml_file):
                     try:
-                        with open(fields_json_file, 'r') as f:
-                            field_data = json.load(f)
-                            json_field_names = list(field_data.get("fields", []))
-                            raw_units = field_data.get("units", {})
+                        with open(fields_yaml_file, 'r') as f:
+                            field_data = yaml.safe_load(f)
+                            yaml_field_names = list(field_data.keys())  # Extract field names
+                            raw_units = field_data
 
                             # Validate field names count
-                            if len(json_field_names) == len(fields):
-                                # Replace the default fields (real_comp*) with those from Fields.json
-                                field_names = json_field_names
+                            if len(yaml_field_names) == len(fields):
+                                # Replace the default fields (real_comp*) with those from Fields.yaml
+                                field_names = yaml_field_names
                                 field_units = {}
 
                                 # Translate raw unit dimensions into readable strings
@@ -1562,13 +1563,13 @@ def _parse_parameter_file(self):
 
                             else:
                                 mylog.debug(
-                                    "Field names in Fields.json do not match the number of fields in Header for '%s'. "
+                                    "Field names in Fields.yaml do not match the number of fields in Header for '%s'. "
                                     "Using names from Header.",
                                     particle_type,
                                 )
                     except Exception as e:
                         mylog.debug(
-                            "Failed to parse Fields.json for particle type '%s': %s",
+                            "Failed to parse Fields.yaml for particle type '%s': %s",
                             particle_type,
                             e,
                         )
@@ -1576,8 +1577,8 @@ def _parse_parameter_file(self):
                     # Ensure all fields have units (handle missing units gracefully)
                     field_units = {field: "dimensionless" for field in fields}
 
-                # Explicitly remove real_comp* entries if Fields.json replaces the fields
-                if field_names == json_field_names:
+                # Explicitly remove real_comp* entries if Fields.yaml replaces the fields
+                if field_names == yaml_field_names:
                     for idx in range(len(fields)):
                         real_comp_field = f"real_comp{idx}"
                         if real_comp_field in field_units:

From 5988844bc8654621066e8a92de473e0b2c0adc9e Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Thu, 20 Mar 2025 19:08:04 +1100
Subject: [PATCH 09/26] add quokka = ["PyYAML>=6.0.1"] (#13)

add quokka = ["PyYAML>=6.0.1"] (#13)
---
 pyproject.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/pyproject.toml b/pyproject.toml
index 937f4a0afd4..b904aa14b20 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -121,6 +121,7 @@ open-pmd = ["yt[HDF5]"]
 owls = ["yt[HDF5]"]
 owls-subfind = ["yt[HDF5]"]
 parthenon = ["yt[HDF5]"]
+quokka = ["PyYAML>=6.0.1"]
 ramses = ["yt[Fortran]", "scipy"]
 rockstar = []
 sdf = ["requests>=2.20.0"]
@@ -181,6 +182,7 @@ full = [
     "yt[owls]",
     "yt[owls_subfind]",
     "yt[parthenon]",
+    "yt[quokka]",
     "yt[ramses]",
     "yt[rockstar]",
     "yt[sdf]",

From 8aad9dcf81270841440b3fc838a8d2b1923d8676 Mon Sep 17 00:00:00 2001
From: Rongjun Huang <u6569836@anu.edu.au>
Date: Thu, 20 Mar 2025 22:39:04 +1100
Subject: [PATCH 10/26] Remove rebundancy (#14)

Co-authored-by: Rongjun-ANU <Rongjun.Huang@alumni.anu.edu.au>
---
 yt/frontends/amrex/data_structures.py | 10 ----------
 1 file changed, 10 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index d72f66207bc..cd7520efe20 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1573,16 +1573,6 @@ def _parse_parameter_file(self):
                             particle_type,
                             e,
                         )
-                else:
-                    # Ensure all fields have units (handle missing units gracefully)
-                    field_units = {field: "dimensionless" for field in fields}
-
-                # Explicitly remove real_comp* entries if Fields.yaml replaces the fields
-                if field_names == yaml_field_names:
-                    for idx in range(len(fields)):
-                        real_comp_field = f"real_comp{idx}"
-                        if real_comp_field in field_units:
-                            del field_units[real_comp_field]
 
                 # Add particle info
                 particle_info[particle_type] = {

From 9c1103d66a0a82ae36310691affd69f80b4bcc41 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Fri, 21 Mar 2025 22:44:46 +1100
Subject: [PATCH 11/26] supported structured metadata.yaml (#15)

Recently I have updated the format of metadata.yaml from Quokka
simulations, and also include `quokka_version` variable in
metadata.yaml. See https://github.com/quokka-astro/quokka/pull/935 .
This PR updates the Quokka frontend to support this new, nested yaml
metadata.

## PR Checklist

<!-- Note that some of these check boxes may not apply to all pull
requests -->
- [ ] New features are documented, with docstrings and narrative docs
- [ ] Adds a test for any bugs fixed. Adds tests for new features.

<!--We understand that PRs can sometimes be overwhelming, especially as
the
reviews start coming in. Please let us know if the reviews are unclear
or the
recommended next step seems overly demanding , or if you would like help
in
addressing a reviewer's comments. And please ping us if you've been
waiting
too long to hear back on your PR.-->
---
 yt/frontends/amrex/data_structures.py | 39 ++++++++++++++++++++++++---
 1 file changed, 35 insertions(+), 4 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index cd7520efe20..73e36c66832 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1606,16 +1606,47 @@ def _parse_metadata_file(self):
             with open(metadata_filename, 'r') as f:
                 # Load metadata using yaml
                 metadata = yaml.safe_load(f)
-                # Update dataset parameters with metadata if it exists
-                if metadata:
+                
+                if not metadata:
+                    mylog.debug("Warning: Metadata file is empty.")
+                    return
+                    
+                # Get quokka version, default to 0.0 if not present
+                quokka_version = metadata.get('quokka_version', '0.0')
+                
+                # Define version comparison function
+                is_newer_or_equal = lambda ver, compare_to: (
+                    tuple(map(int, str(ver).split('.'))) >= 
+                    tuple(map(int, str(compare_to).split('.')))
+                )
+                
+                # Compare with version 25.03 (year 2025, month 3)
+                if not is_newer_or_equal(quokka_version, '25.03'):
+                    # For older versions
                     self.parameters.update(metadata)
-                    mylog.debug("Metadata loaded successfully")
+                    mylog.debug("Metadata loaded successfully (older versions)")
                 else:
-                    mylog.debug("Warning: Metadata file is empty.")
+                    # For newer versions, process each section according to rules
+                    for key, value in metadata.items():
+                        if key == 'units':
+                            # For units section, read everything underneath
+                            for unit_name, unit_value in value.items():
+                                self.parameters[unit_name] = unit_value
+                        elif key == 'constants':
+                            # Ignore constants section
+                            continue
+                        else:
+                            # For everything else, read in directly
+                            self.parameters[key] = value
+                    
+                    mylog.debug("Metadata loaded successfully (version 25.03 or newer)")
+                
         except FileNotFoundError:
             mylog.debug(f"Error: Metadata file '{metadata_filename}' not found.")
         except yaml.YAMLError as e:
             mylog.debug(f"Error parsing metadata file: {e}")
+        except ValueError as e:
+            mylog.debug(f"Error parsing version numbers: {e}")
 
 
 def _guess_pcast(vals):

From 7755781916d9be511b4c50fe4e26842b68d0fea2 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Mon, 24 Mar 2025 11:34:55 +1100
Subject: [PATCH 12/26] add PyYAML to amrex (#16)

## PR Summary

<!--Please provide at least 1-2 sentences describing the pull request in
detail.  Why is this change required?  What problem does it solve?-->

<!--If it fixes an open issue, please link to the issue here.-->

## PR Checklist

<!-- Note that some of these check boxes may not apply to all pull
requests -->
- [ ] New features are documented, with docstrings and narrative docs
- [ ] Adds a test for any bugs fixed. Adds tests for new features.

<!--We understand that PRs can sometimes be overwhelming, especially as
the
reviews start coming in. Please let us know if the reviews are unclear
or the
recommended next step seems overly demanding , or if you would like help
in
addressing a reviewer's comments. And please ping us if you've been
waiting
too long to hear back on your PR.-->
---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index b904aa14b20..8c80640823e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -89,7 +89,7 @@ Fortran = ["f90nml>=1.1"]
 # We also normalize all target names to lower case for consistency.
 adaptahop = []
 ahf = []
-amrex = []
+amrex = ["PyYAML>=6.0.1"]
 amrvac = ["yt[Fortran]"]
 art = []
 arepo = ["yt[HDF5]"]

From 075769d3774ff47adb7ba571a18be0ace3d790d0 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Mon, 24 Mar 2025 15:04:15 +1100
Subject: [PATCH 13/26] Chong/update-loading_data-rst (#17)

update YT docs for Quokka frontend
---
 doc/source/examining/loading_data.rst | 155 ++++++++++++++++++++++++++
 1 file changed, 155 insertions(+)

diff --git a/doc/source/examining/loading_data.rst b/doc/source/examining/loading_data.rst
index 55b209efa1e..629b388377f 100644
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -2133,6 +2133,161 @@ zero.
 * Particles may be difficult to integrate.
 * Data must already reside in memory.
 
+
+.. _loading-quokka-data:
+
+QUOKKA Data
+-----------
+
+`QUOKKA <https://github.com/quokka-astro/quokka>`_ data is supported and cared
+for by ChongChong He and Rongjun Huang. QUOKKA is based on the AMReX framework.
+Uniform-grid and AMR datasets in cartesian coordinates and particle data are
+fully supported.
+
+A standard QUOKKA dataset includes the following components:
+
+.. code-block:: bash
+
+    plt00007/
+    ├── Level_0/         # QUOKKA data
+    ├── Header           # Dataset header information
+    ├── metadata.yaml    # Configuration file
+    ├── XXX_particles/     # Optional, particle data (e.g., Rad or Sink)
+    │   ├── Fields.yaml    # Particle field names and units
+    │   ├── Header
+    │   └── Level_0/
+
+Loading QUOKKA data is straightforward with the ``load`` command:
+
+.. code-block:: python
+
+    import yt
+    ds = yt.load("sample/RadBeam/plt00007")
+
+Boxlib Fields and Units
+^^^^^^^^^^^^^^^^^^^^^^^
+
+The QUOKKA frontend is built upon the AMReX framework, which follows the BoxLib
+data format. When a QUOKKA dataset is loaded, yt automatically:
+
+- Detects and loads all available fields in the dataset from the ``Header`` file
+- Maps the native field names to (field_type, field_name) tuples
+- Assigns physical units based on metadata information
+
+QUOKKA datasets include six mandatory gas fields:
+
+.. code-block:: python
+
+    [('boxlib', 'gasDensity'),
+     ('boxlib', 'gasEnergy'),
+     ('boxlib', 'gasInternalEnergy'),
+     ('boxlib', 'x-GasMomentum'),
+     ('boxlib', 'y-GasMomentum'),
+     ('boxlib', 'z-GasMomentum')]
+
+For simulations with radiative transfer, the following fields are also available:
+
+.. code-block:: python
+
+    [('boxlib', 'radEnergy-Group0'),
+     ('boxlib', 'x-RadFlux-Group0'),
+     ('boxlib', 'y-RadFlux-Group0'),
+     ('boxlib', 'z-RadFlux-Group0'),
+     ('boxlib', 'radEnergy-Group1'),
+     ('boxlib', 'x-RadFlux-Group1'),
+     ('boxlib', 'y-RadFlux-Group1'),
+     ('boxlib', 'z-RadFlux-Group1'),
+     ...
+     ]
+
+To see all available native fields in your dataset, do ``print(ds.field_list)``.
+
+Derived Fields
+^^^^^^^^^^^^^
+
+When a QUOKKA dataset is loaded by yt, the native fields are mapped to derived
+fields with appropriate physical units. The frontend handles:
+
+- Converting dimensionless native fields to physically meaningful units
+- Creating convenient derived fields for analysis and visualization
+
+You can view all available derived fields with ``print(ds.derived_field_list)``.
+For example, the native ('boxlib', 'gasDensity') field is mapped to ('gas',
+'density') with proper units:
+
+.. code-block:: python
+
+    # Accessing density with units
+    print(ds.r[('boxlib', 'gasDensity')])
+    print(ds.r[('gas', 'density')])
+
+.. code-block:: python
+
+    unyt_array([1., 1., 1., ..., 1., 1., 1.], 'code_mass/code_length**3')
+    unyt_array([1., 1., 1., ..., 1., 1., 1.], 'g/cm**3')
+
+Particle Data
+^^^^^^^^^^^^
+
+QUOKKA supports various types of particle data. When particles are present,
+they're stored in dedicated directories:
+
+.. code-block:: bash
+
+    dataset_folder/
+    ├── Level_0/
+    ├── Header
+    ├── metadata.yaml
+    ├── XXX_particles/     # Particle fields (e.g., Rad or Sink or CIC)
+    │   ├── Fields.yaml    # Particle field names and units
+    │   ├── Header
+    │   └── Level_0/
+
+The particle types and fields are automatically detected and loaded. The
+``Fields.yaml`` file defines particle field names and their units, which YT will
+read and automatically convert to physical units.
+
+To access particle data in your analysis:
+
+.. code-block:: python
+
+    # Access particle data
+    ad = ds.all_data()
+
+    # print particle info
+    print(ds.parameters['particle_info'])
+    
+    # List particle fields available for the selected particle type
+    print(ad["Rad_particles", "particle_position_x"])
+    print(ad["Rad_particles", "luminosity"])
+    
+    # Convert to physical units
+    print(ad["Rad_particles", "luminosity"].to("erg/s"))
+
+To annotate particles in a plot, use the ``annotate_particles`` method:
+
+.. code-block:: python
+
+    yt.SlicePlot(ds, 'z', ('gas', 'density')).annotate_particles(
+      1, p_size=400., col='blue', marker='*', ptype='Rad_particles')
+
+Support for Legacy Datasets
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The ``metadata.yaml`` file is the key to determining if a dataset is a QUOKKA
+data object. For data generated by early versions of QUOKKA that might be
+missing this file, you can bypass metadata parsing:
+
+.. code-block:: python
+
+    from yt.frontends.amrex.data_structures import QuokkaDataset
+
+    class OldQuokkaDataset(QuokkaDataset):
+        def _parse_metadata_file(self):
+            pass
+
+    ds_old = OldQuokkaDataset("sample/RadBeam/plt007")
+
 .. _loading-semi-structured-mesh-data:
 
 Semi-Structured Grid Data

From 932b7c46d8c545caad9ec30ea45a2e9ee2b0454e Mon Sep 17 00:00:00 2001
From: Rongjun Huang <u6569836@anu.edu.au>
Date: Sun, 27 Apr 2025 16:08:56 +1000
Subject: [PATCH 14/26] test_quokka (#19)

update the `yt/frontends/amrex/tests/test_outputs.py` to run pytest for
quokka dataset

---------

Co-authored-by: Rongjun-ANU <Rongjun.Huang@alumni.anu.edu.au>
Co-authored-by: ChongChong He <chongchong.he@anu.edu.au>
---
 yt/frontends/amrex/tests/test_outputs.py | 64 ++++++++++++++++++++++++
 1 file changed, 64 insertions(+)

diff --git a/yt/frontends/amrex/tests/test_outputs.py b/yt/frontends/amrex/tests/test_outputs.py
index 14e1031243b..d537777c957 100644
--- a/yt/frontends/amrex/tests/test_outputs.py
+++ b/yt/frontends/amrex/tests/test_outputs.py
@@ -1,3 +1,4 @@
+import os
 import numpy as np
 from numpy.testing import assert_allclose, assert_equal
 
@@ -7,6 +8,7 @@
     MaestroDataset,
     NyxDataset,
     OrionDataset,
+    QuokkaDataset,
     WarpXDataset,
 )
 from yt.loaders import load
@@ -443,3 +445,65 @@ def test_maestro_parameters():
     # Check an int parameter
     assert ds.parameters["s0_interp_type"] == 3
     assert type(ds.parameters["s0_interp_type"]) is int  # noqa: E721
+
+
+quokka = "RadiatingParticles_plt026"
+
+
+@requires_file(quokka)
+def test_quokka():
+    # Test QuokkaDataset instance type
+    ds = data_dir_load(quokka)
+    assert isinstance(ds, QuokkaDataset)
+        
+    # Test basic parameters 
+    assert ds.parameters["HydroMethod"] == "Quokka"
+    assert ds.parameters["quokka_version"] == 25.03
+    assert ds.parameters["dimensionality"] == 2
+    assert_equal(ds.parameters["fields"], [
+        'gasDensity',
+        'x-GasMomentum',
+        'y-GasMomentum',
+        'z-GasMomentum',
+        'gasEnergy',
+        'gasInternalEnergy',
+        'radEnergy-Group0',
+        'x-RadFlux-Group0',
+        'y-RadFlux-Group0',
+        'z-RadFlux-Group0'
+    ])
+        
+    # Check domain dimensions and boundaries
+    assert_equal(ds.domain_dimensions, [64, 64, 1])
+    assert_equal(ds.domain_left_edge, [0.0, 0.0, 0.0])
+    assert_equal(ds.domain_right_edge, [1.0, 1.0, 1.0])
+        
+    # Test derived fields and unit conversions
+    ad = ds.all_data()
+    assert ("gas", "density") in ds.derived_field_list
+    density = ad["gas", "density"]
+    assert str(density.units) == "g/cm**3"
+        
+    # Check momentum/velocity relationship which confirms velocity is correctly derived
+    momentum_x = ad["boxlib", "x-GasMomentum"]
+    gas_density = ad["boxlib", "gasDensity"]
+    velocity_x = ad["gas", "velocity_x"]
+    assert_allclose(momentum_x / gas_density, velocity_x)
+        
+    # Test radiation fields specifically
+    assert ds.parameters["radiation_field_groups"] == 1
+    assert "rad" in ds.fluid_types
+    assert ("rad", "energy_density_0") in ds.derived_field_list
+        
+    # Test particle IO with specific expectations
+    assert ds.parameters["particle_types"] == ("Rad_particles",)
+        
+    # Verify particle fields and units
+    particle_info = ds.parameters["particle_info"]["Rad_particles"]
+    assert particle_info["num_particles"] == 2
+    assert particle_info["num_fields"] == 3
+    assert set(particle_info["fields"]) == {"luminosity", "birth_time", "end_time"}
+    assert particle_info["units"]["luminosity"] == "M^1 L^2 T^-3"
+    assert particle_info["units"]["birth_time"] == "T^1"
+    assert particle_info["units"]["end_time"] == "T^1"
+        
\ No newline at end of file

From 5ad1e7828fb5eaa54eeb9582e43ce2549ec4cb0d Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Thu, 15 May 2025 16:05:27 +1000
Subject: [PATCH 15/26] Add face-centered dataset support (#20)

Add support for reading face-centered dataset. See changes in
loading_data.rst for details.

## PR Summary

<!--Please provide at least 1-2 sentences describing the pull request in
detail.  Why is this change required?  What problem does it solve?-->

<!--If it fixes an open issue, please link to the issue here.-->

## PR Checklist

<!-- Note that some of these check boxes may not apply to all pull
requests -->
- [x] New features are documented, with docstrings and narrative docs
- [ ] Adds a test for any bugs fixed. Adds tests for new features.

<!--We understand that PRs can sometimes be overwhelming, especially as
the
reviews start coming in. Please let us know if the reviews are unclear
or the
recommended next step seems overly demanding , or if you would like help
in
addressing a reviewer's comments. And please ping us if you've been
waiting
too long to hear back on your PR.-->
---
 doc/source/examining/loading_data.rst |  5 +++
 yt/frontends/amrex/data_structures.py | 48 +++++++++++++++++++++++++++
 2 files changed, 53 insertions(+)

diff --git a/doc/source/examining/loading_data.rst b/doc/source/examining/loading_data.rst
index 629b388377f..7ffa4c6e67f 100644
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -2271,6 +2271,11 @@ To annotate particles in a plot, use the ``annotate_particles`` method:
     yt.SlicePlot(ds, 'z', ('gas', 'density')).annotate_particles(
       1, p_size=400., col='blue', marker='*', ptype='Rad_particles')
 
+Face-Centered Variables
+^^^^^^^^^^^^^^^^^^^^^^^
+
+QUOKKA datasets can also include face-centered variables. These are stored in the `fc_vars` directory which is inside the plotfile directory. When loading a Quokka dataset, YT automatically detects and loads any face-centered variables from the `fc_vars/x*****`, `fc_vars/y*****`, and `fc_vars/z*****` subdirectories. These are made available as attributes of the main dataset (`ds.ds_fc_x`, `ds.ds_fc_y`, `ds.ds_fc_z`). Each face-centered dataset is a full YT dataset object with its own fields and mesh structure, and includes a reference back to the parent dataset, `ds.ds_fc_x.parent_ds`. This allows for easy access to both cell-centered and face-centered data in the same analysis workflow.
+
 Support for Legacy Datasets
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 73e36c66832..acd5ae77353 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -17,6 +17,7 @@
 from yt.utilities.io_handler import io_registry
 from yt.utilities.lib.misc_utilities import get_box_grids_level
 from yt.utilities.parallel_tools.parallel_analysis_interface import parallel_root_only
+from yt.loaders import load
 
 from .fields import (
     BoxlibFieldInfo,
@@ -1408,6 +1409,53 @@ def __init__(
         # Add magnetic fields in fluid_types only if detected
         if "Bfield" in self.parameters['fields']:
             self.fluid_types += ("mag",)
+            
+        # Check for face-centered variables directories
+        self._load_face_centered_datasets()
+
+    def _load_face_centered_datasets(self):
+        """
+        Check for face-centered variables directories (fc_vars/x, fc_vars/y, fc_vars/z)
+        and load them as additional datasets if they exist.
+        """
+
+        fc_dir = os.path.join(self.output_dir, "fc_vars")
+        if not os.path.isdir(fc_dir):
+            mylog.debug(f"No face-centered variables directory found at {fc_dir}. Skipping face-centered datasets.")
+            return
+
+        try:
+            out_id = int(os.path.basename(self.output_dir)[-5:])
+        except ValueError:
+            raise ValueError(f"Output directory {self.output_dir} does not end with at least 5 digits. Cannot infer output index. Cannot load face-centered datasets.")
+        
+        # Define the possible face-centered directories
+        fc_dirs = {
+            'x': os.path.join(fc_dir, f"x{out_id:05d}"),
+            'y': os.path.join(fc_dir, f"y{out_id:05d}"),
+            'z': os.path.join(fc_dir, f"z{out_id:05d}"),
+        }
+        
+        # For each direction, check if the directory exists and load it
+        for direction, fc_dir in fc_dirs.items():
+            if os.path.isdir(fc_dir):
+                try:
+                    # Load the face-centered dataset
+                    mylog.info(f"Loading face-centered {direction} dataset from {fc_dir}")
+                    fc_ds = load(fc_dir)
+                    
+                    # Store the dataset as an attribute
+                    setattr(self, f"ds_fc_{direction}", fc_ds)
+                    
+                    # Add a reference back to the parent dataset
+                    fc_ds.parent_ds = self
+                    
+                    # Add information about this being a face-centered dataset
+                    fc_ds.fc_direction = direction
+                except Exception as e:
+                    mylog.warning(f"Failed to load face-centered {direction} dataset: {e}")
+            else:
+                mylog.debug(f"No face-centered {direction} dataset found at {fc_dir}")
 
     def _parse_parameter_file(self):
         # Call parent method to initialize core setup by yt

From 225041047b2b07fb442716e6b26fad0185fc2927 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Thu, 15 May 2025 16:35:58 +1000
Subject: [PATCH 16/26] format

---
 doc/source/examining/loading_data.rst | 15 ++++++++++++---
 1 file changed, 12 insertions(+), 3 deletions(-)

diff --git a/doc/source/examining/loading_data.rst b/doc/source/examining/loading_data.rst
index 7ffa4c6e67f..00ee079673e 100644
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -2256,11 +2256,11 @@ To access particle data in your analysis:
 
     # print particle info
     print(ds.parameters['particle_info'])
-    
+
     # List particle fields available for the selected particle type
     print(ad["Rad_particles", "particle_position_x"])
     print(ad["Rad_particles", "luminosity"])
-    
+
     # Convert to physical units
     print(ad["Rad_particles", "luminosity"].to("erg/s"))
 
@@ -2274,7 +2274,16 @@ To annotate particles in a plot, use the ``annotate_particles`` method:
 Face-Centered Variables
 ^^^^^^^^^^^^^^^^^^^^^^^
 
-QUOKKA datasets can also include face-centered variables. These are stored in the `fc_vars` directory which is inside the plotfile directory. When loading a Quokka dataset, YT automatically detects and loads any face-centered variables from the `fc_vars/x*****`, `fc_vars/y*****`, and `fc_vars/z*****` subdirectories. These are made available as attributes of the main dataset (`ds.ds_fc_x`, `ds.ds_fc_y`, `ds.ds_fc_z`). Each face-centered dataset is a full YT dataset object with its own fields and mesh structure, and includes a reference back to the parent dataset, `ds.ds_fc_x.parent_ds`. This allows for easy access to both cell-centered and face-centered data in the same analysis workflow.
+QUOKKA datasets can also include face-centered variables. These are stored in
+the `fc_vars` directory which is inside the plotfile directory. When loading a
+Quokka dataset, YT automatically detects and loads any face-centered variables
+from the `fc_vars/x*****`, `fc_vars/y*****`, and `fc_vars/z*****`
+subdirectories. These are made available as attributes of the main dataset
+(`ds.ds_fc_x`, `ds.ds_fc_y`, `ds.ds_fc_z`). Each face-centered dataset is a full
+YT dataset object with its own fields and mesh structure, and includes a
+reference back to the parent dataset, `ds.ds_fc_x.parent_ds`. This allows for
+easy access to both cell-centered and face-centered data in the same analysis
+workflow.
 
 Support for Legacy Datasets
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^

From fa2289a9e9a81fe516236a137f0d5b8d74ebe7e8 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Thu, 15 May 2025 17:06:39 +1000
Subject: [PATCH 17/26] fix pre-commit.ci errors

---
 yt/frontends/amrex/data_structures.py | 12 +++++-------
 yt/frontends/amrex/fields.py          |  3 +--
 2 files changed, 6 insertions(+), 9 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 29e69487e54..0b7b87f95b0 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1425,8 +1425,8 @@ def _load_face_centered_datasets(self):
 
         try:
             out_id = int(os.path.basename(self.output_dir)[-5:])
-        except ValueError:
-            raise ValueError(f"Output directory {self.output_dir} does not end with at least 5 digits. Cannot infer output index. Cannot load face-centered datasets.")
+        except ValueError as err:
+            raise ValueError(f"Output directory {self.output_dir} does not end with at least 5 digits. Cannot infer output index. Cannot load face-centered datasets.") from err
         
         # Define the possible face-centered directories
         fc_dirs = {
@@ -1477,7 +1477,7 @@ def _parse_parameter_file(self):
             # Metadata flags
             rad_group_count = 0
 
-            for i in range(num_fields):
+            for _i in range(num_fields):
                 current_field = f.readline().strip()
                 if current_field.startswith("radEnergy-Group"):
                     rad_group_count += 1
@@ -1662,10 +1662,8 @@ def _parse_metadata_file(self):
                 quokka_version = metadata.get('quokka_version', '0.0')
                 
                 # Define version comparison function
-                is_newer_or_equal = lambda ver, compare_to: (
-                    tuple(map(int, str(ver).split('.'))) >= 
-                    tuple(map(int, str(compare_to).split('.')))
-                )
+                def is_newer_or_equal(ver, compare_to):
+                    return tuple(map(int, str(ver).split('.'))) >= tuple(map(int, str(compare_to).split('.')))
                 
                 # Compare with version 25.03 (year 2025, month 3)
                 if not is_newer_or_equal(quokka_version, '25.03'):
diff --git a/yt/frontends/amrex/fields.py b/yt/frontends/amrex/fields.py
index 1b518df929a..06cb8f839b9 100644
--- a/yt/frontends/amrex/fields.py
+++ b/yt/frontends/amrex/fields.py
@@ -584,7 +584,6 @@ def setup_Bfields(self):
             return
 
         for axis in "xyz":
-            bfield_name = f"BField_{axis}"
             boxlib_bfield = f"{axis}-BField"
 
             if ("boxlib", boxlib_bfield) in self.field_list:
@@ -606,7 +605,7 @@ def setup_radiation_fields(self):
                 self.add_field(
                     ("rad", f"energy_density_{group}"),
                     sampling_type="cell",
-                    function=lambda field, data: data["boxlib", energy_field] * self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
+                    function=lambda _, data: data["boxlib", energy_field] * self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
                     units=self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
                     display_name=f"Radiation Energy Density Group {group}",
                 )

From 8935ab795e9314a2fc62b8a20a0d0a227a13c571 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Thu, 15 May 2025 17:12:49 +1000
Subject: [PATCH 18/26] fix lambda error

---
 yt/frontends/amrex/fields.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/yt/frontends/amrex/fields.py b/yt/frontends/amrex/fields.py
index 06cb8f839b9..092a8325081 100644
--- a/yt/frontends/amrex/fields.py
+++ b/yt/frontends/amrex/fields.py
@@ -605,7 +605,7 @@ def setup_radiation_fields(self):
                 self.add_field(
                     ("rad", f"energy_density_{group}"),
                     sampling_type="cell",
-                    function=lambda _, data: data["boxlib", energy_field] * self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
+                    function=lambda _, data, ef=energy_field: data["boxlib", ef] * self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
                     units=self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
                     display_name=f"Radiation Energy Density Group {group}",
                 )

From bf56f910a6459fce0edcce79cb54297711b97dc7 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Thu, 15 May 2025 17:16:24 +1000
Subject: [PATCH 19/26] remove trailing white space

---
 yt/frontends/amrex/data_structures.py | 26 +++++++++++++-------------
 yt/frontends/amrex/fields.py          |  6 +++---
 2 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 0b7b87f95b0..832f5a38759 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1336,7 +1336,7 @@ def __init__(self, ds, dataset_type="quokka_native"):
 
         # Dynamically detect particle types by searching for "*_particles" directories
         particle_dirs = glob.glob(os.path.join(self.ds.output_dir, "*_particles"))
-        
+
         for pdir in particle_dirs:
             ptype = os.path.basename(pdir)
             header_file = os.path.join(pdir, "Header")
@@ -1408,7 +1408,7 @@ def __init__(
         # Add magnetic fields in fluid_types only if detected
         if "Bfield" in self.parameters['fields']:
             self.fluid_types += ("mag",)
-            
+
         # Check for face-centered variables directories
         self._load_face_centered_datasets()
 
@@ -1427,14 +1427,14 @@ def _load_face_centered_datasets(self):
             out_id = int(os.path.basename(self.output_dir)[-5:])
         except ValueError as err:
             raise ValueError(f"Output directory {self.output_dir} does not end with at least 5 digits. Cannot infer output index. Cannot load face-centered datasets.") from err
-        
+
         # Define the possible face-centered directories
         fc_dirs = {
             'x': os.path.join(fc_dir, f"x{out_id:05d}"),
             'y': os.path.join(fc_dir, f"y{out_id:05d}"),
             'z': os.path.join(fc_dir, f"z{out_id:05d}"),
         }
-        
+
         # For each direction, check if the directory exists and load it
         for direction, fc_dir in fc_dirs.items():
             if os.path.isdir(fc_dir):
@@ -1442,13 +1442,13 @@ def _load_face_centered_datasets(self):
                     # Load the face-centered dataset
                     mylog.info(f"Loading face-centered {direction} dataset from {fc_dir}")
                     fc_ds = load(fc_dir)
-                    
+
                     # Store the dataset as an attribute
                     setattr(self, f"ds_fc_{direction}", fc_ds)
-                    
+
                     # Add a reference back to the parent dataset
                     fc_ds.parent_ds = self
-                    
+
                     # Add information about this being a face-centered dataset
                     fc_ds.fc_direction = direction
                 except Exception as e:
@@ -1653,18 +1653,18 @@ def _parse_metadata_file(self):
             with open(metadata_filename, 'r') as f:
                 # Load metadata using yaml
                 metadata = yaml.safe_load(f)
-                
+
                 if not metadata:
                     mylog.debug("Warning: Metadata file is empty.")
                     return
-                    
+
                 # Get quokka version, default to 0.0 if not present
                 quokka_version = metadata.get('quokka_version', '0.0')
-                
+
                 # Define version comparison function
                 def is_newer_or_equal(ver, compare_to):
                     return tuple(map(int, str(ver).split('.'))) >= tuple(map(int, str(compare_to).split('.')))
-                
+
                 # Compare with version 25.03 (year 2025, month 3)
                 if not is_newer_or_equal(quokka_version, '25.03'):
                     # For older versions
@@ -1683,9 +1683,9 @@ def is_newer_or_equal(ver, compare_to):
                         else:
                             # For everything else, read in directly
                             self.parameters[key] = value
-                    
+
                     mylog.debug("Metadata loaded successfully (version 25.03 or newer)")
-                
+
         except FileNotFoundError:
             mylog.debug(f"Error: Metadata file '{metadata_filename}' not found.")
         except yaml.YAMLError as e:
diff --git a/yt/frontends/amrex/fields.py b/yt/frontends/amrex/fields.py
index 092a8325081..163b4238d9d 100644
--- a/yt/frontends/amrex/fields.py
+++ b/yt/frontends/amrex/fields.py
@@ -512,7 +512,7 @@ class QuokkaFieldInfo(FieldInfoContainer):
         ("y-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_y"], r"\rho v")),
         ("z-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_z"], r"\rho w")),
         # Temperature field is not present in early Quokka datasets
-        ("gasTemperature", ("K", ["temperature"], r"T")), 
+        ("gasTemperature", ("K", ["temperature"], r"T")),
         # Scalar fields are not always present
         ("scalar_0", ("", ["scalar_0"], "Scalar 0")),
         ("scalar_1", ("", ["scalar_1"], "Scalar 1")),
@@ -621,14 +621,14 @@ def setup_radiation_fields(self):
                         units=self.ds.unit_system["energy"] / self.ds.unit_system["length"]**2 / self.ds.unit_system["time"],
                         display_name=f"Radiation Flux Density {axis.upper()} Group {group}",
                     )
-        
+
         # Dynamically set up custom particle fields (convert `real_comp` to physics) for all particle types
         for ptype in self.ds.particle_types:
             self.setup_custom_particle_fields(ptype)
 
     def setup_custom_particle_fields(self, ptype):
         """
-        Dynamically set up custom particle fields (real_comp) for the given particle type, 
+        Dynamically set up custom particle fields (real_comp) for the given particle type,
         interpreting dimensional expressions and applying appropriate unit conversions.
 
         Parameters:

From 9997397b58f9ec3be83b37e67e8aadd32ae39c2b Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Fri, 16 May 2025 12:39:07 +1000
Subject: [PATCH 20/26] fix style issue

---
 doc/source/examining/loading_data.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/source/examining/loading_data.rst b/doc/source/examining/loading_data.rst
index 00ee079673e..cf6f5e35742 100644
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -2275,13 +2275,13 @@ Face-Centered Variables
 ^^^^^^^^^^^^^^^^^^^^^^^
 
 QUOKKA datasets can also include face-centered variables. These are stored in
-the `fc_vars` directory which is inside the plotfile directory. When loading a
+the ``fc_vars`` directory which is inside the plotfile directory. When loading a
 Quokka dataset, YT automatically detects and loads any face-centered variables
-from the `fc_vars/x*****`, `fc_vars/y*****`, and `fc_vars/z*****`
+from the ``fc_vars/x*****``, ``fc_vars/y*****``, and ``fc_vars/z*****``
 subdirectories. These are made available as attributes of the main dataset
 (`ds.ds_fc_x`, `ds.ds_fc_y`, `ds.ds_fc_z`). Each face-centered dataset is a full
 YT dataset object with its own fields and mesh structure, and includes a
-reference back to the parent dataset, `ds.ds_fc_x.parent_ds`. This allows for
+reference back to the parent dataset, ``ds.ds_fc_x.parent_ds``. This allows for
 easy access to both cell-centered and face-centered data in the same analysis
 workflow.
 

From 145d0eac4be89e6ad64d3b722113936f1923bcdd Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Fri, 16 May 2025 12:42:07 +1000
Subject: [PATCH 21/26] fix style

---
 doc/source/examining/loading_data.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/source/examining/loading_data.rst b/doc/source/examining/loading_data.rst
index cf6f5e35742..042d1d0b192 100644
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -2279,7 +2279,7 @@ the ``fc_vars`` directory which is inside the plotfile directory. When loading a
 Quokka dataset, YT automatically detects and loads any face-centered variables
 from the ``fc_vars/x*****``, ``fc_vars/y*****``, and ``fc_vars/z*****``
 subdirectories. These are made available as attributes of the main dataset
-(`ds.ds_fc_x`, `ds.ds_fc_y`, `ds.ds_fc_z`). Each face-centered dataset is a full
+(``ds.ds_fc_x``, ``ds.ds_fc_y``, ``ds.ds_fc_z``). Each face-centered dataset is a full
 YT dataset object with its own fields and mesh structure, and includes a
 reference back to the parent dataset, ``ds.ds_fc_x.parent_ds``. This allows for
 easy access to both cell-centered and face-centered data in the same analysis

From 1542a026a3d2d30256cf10a4e081e704069c5535 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Fri, 16 May 2025 12:47:24 +1000
Subject: [PATCH 22/26] pre-commit run all

---
 doc/source/examining/loading_data.rst    |  53 ++++----
 yt/frontends/amrex/api.py                |   2 +-
 yt/frontends/amrex/data_structures.py    | 153 ++++++++++++++---------
 yt/frontends/amrex/fields.py             |  94 +++++++++++---
 yt/frontends/amrex/tests/test_outputs.py |  44 +++----
 5 files changed, 227 insertions(+), 119 deletions(-)

diff --git a/doc/source/examining/loading_data.rst b/doc/source/examining/loading_data.rst
index 042d1d0b192..c02dc2ae78d 100644
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -2162,6 +2162,7 @@ Loading QUOKKA data is straightforward with the ``load`` command:
 .. code-block:: python
 
     import yt
+
     ds = yt.load("sample/RadBeam/plt00007")
 
 Boxlib Fields and Units
@@ -2178,27 +2179,30 @@ QUOKKA datasets include six mandatory gas fields:
 
 .. code-block:: python
 
-    [('boxlib', 'gasDensity'),
-     ('boxlib', 'gasEnergy'),
-     ('boxlib', 'gasInternalEnergy'),
-     ('boxlib', 'x-GasMomentum'),
-     ('boxlib', 'y-GasMomentum'),
-     ('boxlib', 'z-GasMomentum')]
+    [
+        ("boxlib", "gasDensity"),
+        ("boxlib", "gasEnergy"),
+        ("boxlib", "gasInternalEnergy"),
+        ("boxlib", "x-GasMomentum"),
+        ("boxlib", "y-GasMomentum"),
+        ("boxlib", "z-GasMomentum"),
+    ]
 
 For simulations with radiative transfer, the following fields are also available:
 
 .. code-block:: python
 
-    [('boxlib', 'radEnergy-Group0'),
-     ('boxlib', 'x-RadFlux-Group0'),
-     ('boxlib', 'y-RadFlux-Group0'),
-     ('boxlib', 'z-RadFlux-Group0'),
-     ('boxlib', 'radEnergy-Group1'),
-     ('boxlib', 'x-RadFlux-Group1'),
-     ('boxlib', 'y-RadFlux-Group1'),
-     ('boxlib', 'z-RadFlux-Group1'),
-     ...
-     ]
+    [
+        ("boxlib", "radEnergy-Group0"),
+        ("boxlib", "x-RadFlux-Group0"),
+        ("boxlib", "y-RadFlux-Group0"),
+        ("boxlib", "z-RadFlux-Group0"),
+        ("boxlib", "radEnergy-Group1"),
+        ("boxlib", "x-RadFlux-Group1"),
+        ("boxlib", "y-RadFlux-Group1"),
+        ("boxlib", "z-RadFlux-Group1"),
+        ...,
+    ]
 
 To see all available native fields in your dataset, do ``print(ds.field_list)``.
 
@@ -2218,13 +2222,13 @@ For example, the native ('boxlib', 'gasDensity') field is mapped to ('gas',
 .. code-block:: python
 
     # Accessing density with units
-    print(ds.r[('boxlib', 'gasDensity')])
-    print(ds.r[('gas', 'density')])
+    print(ds.r[("boxlib", "gasDensity")])
+    print(ds.r[("gas", "density")])
 
 .. code-block:: python
 
-    unyt_array([1., 1., 1., ..., 1., 1., 1.], 'code_mass/code_length**3')
-    unyt_array([1., 1., 1., ..., 1., 1., 1.], 'g/cm**3')
+    unyt_array([1.0, 1.0, 1.0, ..., 1.0, 1.0, 1.0], "code_mass/code_length**3")
+    unyt_array([1.0, 1.0, 1.0, ..., 1.0, 1.0, 1.0], "g/cm**3")
 
 Particle Data
 ^^^^^^^^^^^^
@@ -2255,7 +2259,7 @@ To access particle data in your analysis:
     ad = ds.all_data()
 
     # print particle info
-    print(ds.parameters['particle_info'])
+    print(ds.parameters["particle_info"])
 
     # List particle fields available for the selected particle type
     print(ad["Rad_particles", "particle_position_x"])
@@ -2268,8 +2272,9 @@ To annotate particles in a plot, use the ``annotate_particles`` method:
 
 .. code-block:: python
 
-    yt.SlicePlot(ds, 'z', ('gas', 'density')).annotate_particles(
-      1, p_size=400., col='blue', marker='*', ptype='Rad_particles')
+    yt.SlicePlot(ds, "z", ("gas", "density")).annotate_particles(
+        1, p_size=400.0, col="blue", marker="*", ptype="Rad_particles"
+    )
 
 Face-Centered Variables
 ^^^^^^^^^^^^^^^^^^^^^^^
@@ -2296,10 +2301,12 @@ missing this file, you can bypass metadata parsing:
 
     from yt.frontends.amrex.data_structures import QuokkaDataset
 
+
     class OldQuokkaDataset(QuokkaDataset):
         def _parse_metadata_file(self):
             pass
 
+
     ds_old = OldQuokkaDataset("sample/RadBeam/plt007")
 
 .. _loading-semi-structured-mesh-data:
diff --git a/yt/frontends/amrex/api.py b/yt/frontends/amrex/api.py
index e7e3727188e..c1e0161772b 100644
--- a/yt/frontends/amrex/api.py
+++ b/yt/frontends/amrex/api.py
@@ -21,7 +21,7 @@
     CastroFieldInfo,
     MaestroFieldInfo,
     NyxFieldInfo,
-    QuokkaFieldInfo,    
+    QuokkaFieldInfo,
     WarpXFieldInfo,
 )
 from .io import IOHandlerBoxlib
diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 832f5a38759..7a27f4efbf4 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -1,4 +1,3 @@
-import yaml
 import glob
 import os
 import re
@@ -7,6 +6,7 @@
 from stat import ST_CTIME
 
 import numpy as np
+import yaml
 
 from yt.data_objects.index_subobjects.grid_patch import AMRGridPatch
 from yt.data_objects.static_output import Dataset
@@ -14,10 +14,10 @@
 from yt.funcs import mylog, setdefaultattr
 from yt.geometry.api import Geometry
 from yt.geometry.grid_geometry_handler import GridIndex
+from yt.loaders import load
 from yt.utilities.io_handler import io_registry
 from yt.utilities.lib.misc_utilities import get_box_grids_level
 from yt.utilities.parallel_tools.parallel_analysis_interface import parallel_root_only
-from yt.loaders import load
 
 from .fields import (
     BoxlibFieldInfo,
@@ -1343,7 +1343,7 @@ def __init__(self, ds, dataset_type="quokka_native"):
 
             if os.path.exists(header_file):
                 # Read the Header to determine the number of fields
-                with open(header_file, 'r') as f:
+                with open(header_file) as f:
                     _ = f.readline().strip()  # Skip version info
                     _ = f.readline().strip()  # Skip number of particles
                     num_fields = int(f.readline().strip())  # Number of fields
@@ -1361,7 +1361,7 @@ def __init__(self, ds, dataset_type="quokka_native"):
                 is_checkpoint = False
 
                 # Pass the fields to _read_particles
-                extra_fields = quokka_extra_real_fields[:ds.dimensionality] + fields
+                extra_fields = quokka_extra_real_fields[: ds.dimensionality] + fields
 
                 # Read particle data for this type
                 self._read_particles(ptype, is_checkpoint, extra_fields)
@@ -1403,10 +1403,10 @@ def __init__(
         self._parse_metadata_file()
 
         # Add radiation fields in fluid_types only if detected
-        if self.parameters.get('radiation_field_groups', 0) > 0:
+        if self.parameters.get("radiation_field_groups", 0) > 0:
             self.fluid_types += ("rad",)
         # Add magnetic fields in fluid_types only if detected
-        if "Bfield" in self.parameters['fields']:
+        if "Bfield" in self.parameters["fields"]:
             self.fluid_types += ("mag",)
 
         # Check for face-centered variables directories
@@ -1420,19 +1420,23 @@ def _load_face_centered_datasets(self):
 
         fc_dir = os.path.join(self.output_dir, "fc_vars")
         if not os.path.isdir(fc_dir):
-            mylog.debug(f"No face-centered variables directory found at {fc_dir}. Skipping face-centered datasets.")
+            mylog.debug(
+                f"No face-centered variables directory found at {fc_dir}. Skipping face-centered datasets."
+            )
             return
 
         try:
             out_id = int(os.path.basename(self.output_dir)[-5:])
         except ValueError as err:
-            raise ValueError(f"Output directory {self.output_dir} does not end with at least 5 digits. Cannot infer output index. Cannot load face-centered datasets.") from err
+            raise ValueError(
+                f"Output directory {self.output_dir} does not end with at least 5 digits. Cannot infer output index. Cannot load face-centered datasets."
+            ) from err
 
         # Define the possible face-centered directories
         fc_dirs = {
-            'x': os.path.join(fc_dir, f"x{out_id:05d}"),
-            'y': os.path.join(fc_dir, f"y{out_id:05d}"),
-            'z': os.path.join(fc_dir, f"z{out_id:05d}"),
+            "x": os.path.join(fc_dir, f"x{out_id:05d}"),
+            "y": os.path.join(fc_dir, f"y{out_id:05d}"),
+            "z": os.path.join(fc_dir, f"z{out_id:05d}"),
         }
 
         # For each direction, check if the directory exists and load it
@@ -1440,7 +1444,9 @@ def _load_face_centered_datasets(self):
             if os.path.isdir(fc_dir):
                 try:
                     # Load the face-centered dataset
-                    mylog.info(f"Loading face-centered {direction} dataset from {fc_dir}")
+                    mylog.info(
+                        f"Loading face-centered {direction} dataset from {fc_dir}"
+                    )
                     fc_ds = load(fc_dir)
 
                     # Store the dataset as an attribute
@@ -1452,7 +1458,9 @@ def _load_face_centered_datasets(self):
                     # Add information about this being a face-centered dataset
                     fc_ds.fc_direction = direction
                 except Exception as e:
-                    mylog.warning(f"Failed to load face-centered {direction} dataset: {e}")
+                    mylog.warning(
+                        f"Failed to load face-centered {direction} dataset: {e}"
+                    )
             else:
                 mylog.debug(f"No face-centered {direction} dataset found at {fc_dir}")
 
@@ -1465,14 +1473,14 @@ def _parse_parameter_file(self):
         if not os.path.exists(header_filename):
             raise FileNotFoundError(f"Header file not found: {header_filename}")
 
-        with open(header_filename, 'r') as f:
+        with open(header_filename) as f:
             # Parse header version
-            self.parameters['plot_file_type'] = f.readline().strip()
+            self.parameters["plot_file_type"] = f.readline().strip()
 
             # expecting mandatory fields: gasDensity, x-GasMomentum, y-GasMomentum, z-GasMomentum, gasEnergy, gasInternalEnergy
             # Number of fields
             num_fields = int(f.readline().strip())
-            self.parameters['fields'] = []
+            self.parameters["fields"] = []
 
             # Metadata flags
             rad_group_count = 0
@@ -1481,49 +1489,63 @@ def _parse_parameter_file(self):
                 current_field = f.readline().strip()
                 if current_field.startswith("radEnergy-Group"):
                     rad_group_count += 1
-                self.parameters['fields'].append(current_field)
+                self.parameters["fields"].append(current_field)
 
             # Add metadata for radiation groups, scalars, and field existence flags
-            self.parameters['radiation_field_groups'] = rad_group_count
+            self.parameters["radiation_field_groups"] = rad_group_count
 
             # Parse remaining metadata
-            self.parameters['dimensionality'] = int(f.readline().strip())  # Dimensionality
-            self.parameters['current_time'] = float(f.readline().strip())  # Simulation time
-            self.parameters['refinement_level'] = int(f.readline().strip())  # Refinement levels
+            self.parameters["dimensionality"] = int(
+                f.readline().strip()
+            )  # Dimensionality
+            self.parameters["current_time"] = float(
+                f.readline().strip()
+            )  # Simulation time
+            self.parameters["refinement_level"] = int(
+                f.readline().strip()
+            )  # Refinement levels
 
             # Domain edges
-            self.parameters['domain_left_edge'] = list(map(float, f.readline().strip().split()))
-            self.parameters['domain_right_edge'] = list(map(float, f.readline().strip().split()))
+            self.parameters["domain_left_edge"] = list(
+                map(float, f.readline().strip().split())
+            )
+            self.parameters["domain_right_edge"] = list(
+                map(float, f.readline().strip().split())
+            )
 
             # Skip empty line
             f.readline()
 
             # Grid info
-            self.parameters['grid_info'] = f.readline().strip()
+            self.parameters["grid_info"] = f.readline().strip()
 
             # Timestamp
-            self.parameters['timestamp'] = list(map(int, f.readline().strip().split()))
+            self.parameters["timestamp"] = list(map(int, f.readline().strip().split()))
 
             # Grid sizes for all refinement levels
-            self.parameters['grid_sizes'] = []
-            for _ in range(self.parameters['refinement_level'] + 1):
-                self.parameters['grid_sizes'].append(list(map(float, f.readline().strip().split())))
+            self.parameters["grid_sizes"] = []
+            for _ in range(self.parameters["refinement_level"] + 1):
+                self.parameters["grid_sizes"].append(
+                    list(map(float, f.readline().strip().split()))
+                )
 
             # Skip placeholders
             f.readline()  # Placeholder line 1
             f.readline()  # Placeholder line 2
 
             # Parse data for all refinement levels
-            self.parameters['refinement_details'] = []
+            self.parameters["refinement_details"] = []
             max_refinement_level = -1
 
-            for level in range(self.parameters['refinement_level'] + 1):
+            for level in range(self.parameters["refinement_level"] + 1):
                 level_data = {}
 
                 # Parse metadata line
                 metadata = f.readline().strip().split()
                 level_data["level"] = int(metadata[0])
-                level_data["num_boxes"] = int(metadata[1])  # Number of boxes in current level
+                level_data["num_boxes"] = int(
+                    metadata[1]
+                )  # Number of boxes in current level
                 level_data["current_time"] = float(metadata[2])
 
                 # Update max refinement level
@@ -1536,13 +1558,16 @@ def _parse_parameter_file(self):
                 level_data["boxes"] = []
                 for _ in range(level_data["num_boxes"]):
                     box = {}
-                    for axis in range(self.parameters['dimensionality']):
+                    for axis in range(self.parameters["dimensionality"]):
                         left_edge, right_edge = map(float, f.readline().strip().split())
-                        box[f"axis_{axis}"] = {"left_edge": left_edge, "right_edge": right_edge}
+                        box[f"axis_{axis}"] = {
+                            "left_edge": left_edge,
+                            "right_edge": right_edge,
+                        }
                     level_data["boxes"].append(box)
 
                 # Append level data to refinement levels
-                self.parameters['refinement_details'].append(level_data)
+                self.parameters["refinement_details"].append(level_data)
 
                 # Skip level marker (e.g., Level_0/Cell)
                 marker = f.readline().strip()
@@ -1567,15 +1592,17 @@ def _parse_parameter_file(self):
                 detected_particle_types.append(particle_type)
 
                 # Parse the Header
-                with open(header_file, 'r') as f:
+                with open(header_file) as f:
                     f.readline().strip()  # Skip version line
                     num_particles = int(f.readline().strip())  # Second line
                     num_fields = int(f.readline().strip())  # Third line
-                    fields = [f.readline().strip() for _ in range(num_fields)]  # Remaining lines
+                    fields = [
+                        f.readline().strip() for _ in range(num_fields)
+                    ]  # Remaining lines
 
                 # Default field names and units
                 field_names = fields[:]
-                field_units = {field: "dimensionless" for field in fields}
+                field_units = dict.fromkeys(fields, "dimensionless")
 
                 # Initialize variables
                 yaml_field_names = None
@@ -1583,9 +1610,11 @@ def _parse_parameter_file(self):
                 # Check and parse Fields.yaml
                 if os.path.exists(fields_yaml_file):
                     try:
-                        with open(fields_yaml_file, 'r') as f:
+                        with open(fields_yaml_file) as f:
                             field_data = yaml.safe_load(f)
-                            yaml_field_names = list(field_data.keys())  # Extract field names
+                            yaml_field_names = list(
+                                field_data.keys()
+                            )  # Extract field names
                             raw_units = field_data
 
                             # Validate field names count
@@ -1595,12 +1624,20 @@ def _parse_parameter_file(self):
                                 field_units = {}
 
                                 # Translate raw unit dimensions into readable strings
-                                base_units = ["M", "L", "T", "Θ"]  # Mass, Length, Time, Temperature
+                                base_units = [
+                                    "M",
+                                    "L",
+                                    "T",
+                                    "Θ",
+                                ]  # Mass, Length, Time, Temperature
                                 for field, unit_dims in raw_units.items():
-                                    field_units[field] = " ".join(
-                                        f"{base_units[i]}^{exp}" if exp != 0 else ""
-                                        for i, exp in enumerate(unit_dims)
-                                    ).strip() or "dimensionless"
+                                    field_units[field] = (
+                                        " ".join(
+                                            f"{base_units[i]}^{exp}" if exp != 0 else ""
+                                            for i, exp in enumerate(unit_dims)
+                                        ).strip()
+                                        or "dimensionless"
+                                    )
 
                                 # Remove any `real_comp*` entries from the field_units
                                 for idx in range(len(fields)):
@@ -1630,11 +1667,13 @@ def _parse_parameter_file(self):
                 }
 
         # Update parameters with particle info
-        self.parameters.update({
-            "particles": len(detected_particle_types),
-            "particle_types": tuple(detected_particle_types),
-            "particle_info": particle_info,
-        })
+        self.parameters.update(
+            {
+                "particles": len(detected_particle_types),
+                "particle_types": tuple(detected_particle_types),
+                "particle_info": particle_info,
+            }
+        )
 
         # Debugging: Log parameters for verification
         mylog.debug("Updated ds.parameters with particle info: %s", self.parameters)
@@ -1650,7 +1689,7 @@ def _parse_metadata_file(self):
         # Construct the full path to the metadata file
         metadata_filename = os.path.join(self.output_dir, self.cparam_filename)
         try:
-            with open(metadata_filename, 'r') as f:
+            with open(metadata_filename) as f:
                 # Load metadata using yaml
                 metadata = yaml.safe_load(f)
 
@@ -1659,25 +1698,27 @@ def _parse_metadata_file(self):
                     return
 
                 # Get quokka version, default to 0.0 if not present
-                quokka_version = metadata.get('quokka_version', '0.0')
+                quokka_version = metadata.get("quokka_version", "0.0")
 
                 # Define version comparison function
                 def is_newer_or_equal(ver, compare_to):
-                    return tuple(map(int, str(ver).split('.'))) >= tuple(map(int, str(compare_to).split('.')))
+                    return tuple(map(int, str(ver).split("."))) >= tuple(
+                        map(int, str(compare_to).split("."))
+                    )
 
                 # Compare with version 25.03 (year 2025, month 3)
-                if not is_newer_or_equal(quokka_version, '25.03'):
+                if not is_newer_or_equal(quokka_version, "25.03"):
                     # For older versions
                     self.parameters.update(metadata)
                     mylog.debug("Metadata loaded successfully (older versions)")
                 else:
                     # For newer versions, process each section according to rules
                     for key, value in metadata.items():
-                        if key == 'units':
+                        if key == "units":
                             # For units section, read everything underneath
                             for unit_name, unit_value in value.items():
                                 self.parameters[unit_name] = unit_value
-                        elif key == 'constants':
+                        elif key == "constants":
                             # Ignore constants section
                             continue
                         else:
@@ -1959,4 +2000,4 @@ def _set_code_unit_attributes(self):
         setdefaultattr(self, "mass_unit", self.quan(1.0, "kg"))
         setdefaultattr(self, "time_unit", self.quan(1.0, "s"))
         setdefaultattr(self, "velocity_unit", self.quan(1.0, "m/s"))
-        setdefaultattr(self, "magnetic_unit", self.quan(1.0, "T"))
\ No newline at end of file
+        setdefaultattr(self, "magnetic_unit", self.quan(1.0, "T"))
diff --git a/yt/frontends/amrex/fields.py b/yt/frontends/amrex/fields.py
index 163b4238d9d..6a891cc346f 100644
--- a/yt/frontends/amrex/fields.py
+++ b/yt/frontends/amrex/fields.py
@@ -503,14 +503,50 @@ def setup_fluid_fields(self):
                     units=unit_system["frequency"],
                 )
 
+
 class QuokkaFieldInfo(FieldInfoContainer):
     known_other_fields: KnownFieldsT = (
         ("gasDensity", ("code_mass / code_length**3", ["density"], r"\rho")),
-        ("gasEnergy", ("code_mass / (code_length * code_time**2)", ["total_energy_density"], r"\rho E")),
-        ("gasInternalEnergy", ("code_mass / (code_length * code_time**2)", ["internal_energy_density"], r"\rho e")),
-        ("x-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_x"], r"\rho u")),
-        ("y-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_y"], r"\rho v")),
-        ("z-GasMomentum", ("code_mass / (code_length**2 * code_time)", ["momentum_density_z"], r"\rho w")),
+        (
+            "gasEnergy",
+            (
+                "code_mass / (code_length * code_time**2)",
+                ["total_energy_density"],
+                r"\rho E",
+            ),
+        ),
+        (
+            "gasInternalEnergy",
+            (
+                "code_mass / (code_length * code_time**2)",
+                ["internal_energy_density"],
+                r"\rho e",
+            ),
+        ),
+        (
+            "x-GasMomentum",
+            (
+                "code_mass / (code_length**2 * code_time)",
+                ["momentum_density_x"],
+                r"\rho u",
+            ),
+        ),
+        (
+            "y-GasMomentum",
+            (
+                "code_mass / (code_length**2 * code_time)",
+                ["momentum_density_y"],
+                r"\rho v",
+            ),
+        ),
+        (
+            "z-GasMomentum",
+            (
+                "code_mass / (code_length**2 * code_time)",
+                ["momentum_density_z"],
+                r"\rho w",
+            ),
+        ),
         # Temperature field is not present in early Quokka datasets
         ("gasTemperature", ("K", ["temperature"], r"T")),
         # Scalar fields are not always present
@@ -540,7 +576,9 @@ def setup_fluid_fields(self):
         def _get_cell_velocity(axis):
             def velocity(field, data):
                 # Divide momentum by density for cell-centered velocity
-                return data["boxlib", f"{axis}-GasMomentum"] / data["boxlib", "gasDensity"]
+                return (
+                    data["boxlib", f"{axis}-GasMomentum"] / data["boxlib", "gasDensity"]
+                )
 
             return velocity
 
@@ -566,13 +604,17 @@ def velocity(field, data):
                 self.add_field(
                     ("gas", face_velocity_name),
                     sampling_type="cell",
-                    function=lambda field, data, axis=axis: data["boxlib", f"{axis}-velocity"] * self.ds.unit_system["length"] / self.ds.unit_system["time"],
+                    function=lambda field, data, axis=axis: data[
+                        "boxlib", f"{axis}-velocity"
+                    ]
+                    * self.ds.unit_system["length"]
+                    / self.ds.unit_system["time"],
                     units=self.ds.unit_system["length"] / self.ds.unit_system["time"],
                     display_name=f"v_{axis} (face-centered)",
                 )
 
         # Call the Bfields and radiation fields setup
-        self.setup_Bfields() # magnetic fields are still placeholder here
+        self.setup_Bfields()  # magnetic fields are still placeholder here
         self.setup_radiation_fields()
 
     def setup_Bfields(self):
@@ -580,7 +622,7 @@ def setup_Bfields(self):
         Dynamically add magnetic fields based on presence of Bfield fields in ds.parameters['fields']
         """
         # Check if any field name contains 'Bfield'
-        if not any('Bfield' in field for field in self.ds.parameters.get('fields', [])):
+        if not any("Bfield" in field for field in self.ds.parameters.get("fields", [])):
             return
 
         for axis in "xyz":
@@ -590,7 +632,10 @@ def setup_Bfields(self):
                 self.add_field(
                     ("mag", f"{axis}-field"),
                     sampling_type="cell",
-                    function=lambda field, data, axis=axis: data["boxlib", f"{axis}-BField"] * self.ds.unit_system["magnetic_field_strength"],
+                    function=lambda field, data, axis=axis: data[
+                        "boxlib", f"{axis}-BField"
+                    ]
+                    * self.ds.unit_system["magnetic_field_strength"],
                     units=self.ds.unit_system["magnetic_field_strength"],
                     display_name=f"B_{axis} (magnetic field)",
                 )
@@ -605,8 +650,11 @@ def setup_radiation_fields(self):
                 self.add_field(
                     ("rad", f"energy_density_{group}"),
                     sampling_type="cell",
-                    function=lambda _, data, ef=energy_field: data["boxlib", ef] * self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
-                    units=self.ds.unit_system["energy"] / self.ds.unit_system["length"]**3,
+                    function=lambda _, data, ef=energy_field: data["boxlib", ef]
+                    * self.ds.unit_system["energy"]
+                    / self.ds.unit_system["length"] ** 3,
+                    units=self.ds.unit_system["energy"]
+                    / self.ds.unit_system["length"] ** 3,
                     display_name=f"Radiation Energy Density Group {group}",
                 )
 
@@ -617,8 +665,15 @@ def setup_radiation_fields(self):
                     self.add_field(
                         ("rad", f"flux_density_{axis}_{group}"),
                         sampling_type="cell",
-                        function=lambda field, data, flux_field=flux_field: data["boxlib", flux_field] * self.ds.unit_system["energy"] / self.ds.unit_system["length"]**2 / self.ds.unit_system["time"],
-                        units=self.ds.unit_system["energy"] / self.ds.unit_system["length"]**2 / self.ds.unit_system["time"],
+                        function=lambda field, data, flux_field=flux_field: data[
+                            "boxlib", flux_field
+                        ]
+                        * self.ds.unit_system["energy"]
+                        / self.ds.unit_system["length"] ** 2
+                        / self.ds.unit_system["time"],
+                        units=self.ds.unit_system["energy"]
+                        / self.ds.unit_system["length"] ** 2
+                        / self.ds.unit_system["time"],
                         display_name=f"Radiation Flux Density {axis.upper()} Group {group}",
                     )
 
@@ -690,10 +745,13 @@ def interpret_units(dim_expr):
                 self.add_field(
                     (ptype, field_name),
                     sampling_type="particle",
-                    function=lambda field, data, real_comp_field=real_comp_field, conv=conversion_factor: (
-                        data[ptype, real_comp_field] * conv
-                    ),
-                    units=conversion_factor.units if hasattr(conversion_factor, "units") else "",
+                    function=lambda field,
+                    data,
+                    real_comp_field=real_comp_field,
+                    conv=conversion_factor: (data[ptype, real_comp_field] * conv),
+                    units=conversion_factor.units
+                    if hasattr(conversion_factor, "units")
+                    else "",
                     display_name=field_name.replace("_", " ").capitalize(),
                 )
 
diff --git a/yt/frontends/amrex/tests/test_outputs.py b/yt/frontends/amrex/tests/test_outputs.py
index 457cba2b030..681b6025c61 100644
--- a/yt/frontends/amrex/tests/test_outputs.py
+++ b/yt/frontends/amrex/tests/test_outputs.py
@@ -1,4 +1,3 @@
-import os
 import numpy as np
 from numpy.testing import assert_allclose, assert_equal
 
@@ -455,49 +454,52 @@ def test_quokka():
     # Test QuokkaDataset instance type
     ds = data_dir_load(quokka)
     assert isinstance(ds, QuokkaDataset)
-        
-    # Test basic parameters 
+
+    # Test basic parameters
     assert ds.parameters["HydroMethod"] == "Quokka"
     assert ds.parameters["quokka_version"] == 25.03
     assert ds.parameters["dimensionality"] == 2
-    assert_equal(ds.parameters["fields"], [
-        'gasDensity',
-        'x-GasMomentum',
-        'y-GasMomentum',
-        'z-GasMomentum',
-        'gasEnergy',
-        'gasInternalEnergy',
-        'radEnergy-Group0',
-        'x-RadFlux-Group0',
-        'y-RadFlux-Group0',
-        'z-RadFlux-Group0'
-    ])
-        
+    assert_equal(
+        ds.parameters["fields"],
+        [
+            "gasDensity",
+            "x-GasMomentum",
+            "y-GasMomentum",
+            "z-GasMomentum",
+            "gasEnergy",
+            "gasInternalEnergy",
+            "radEnergy-Group0",
+            "x-RadFlux-Group0",
+            "y-RadFlux-Group0",
+            "z-RadFlux-Group0",
+        ],
+    )
+
     # Check domain dimensions and boundaries
     assert_equal(ds.domain_dimensions, [64, 64, 1])
     assert_equal(ds.domain_left_edge, [0.0, 0.0, 0.0])
     assert_equal(ds.domain_right_edge, [1.0, 1.0, 1.0])
-        
+
     # Test derived fields and unit conversions
     ad = ds.all_data()
     assert ("gas", "density") in ds.derived_field_list
     density = ad["gas", "density"]
     assert str(density.units) == "g/cm**3"
-        
+
     # Check momentum/velocity relationship which confirms velocity is correctly derived
     momentum_x = ad["boxlib", "x-GasMomentum"]
     gas_density = ad["boxlib", "gasDensity"]
     velocity_x = ad["gas", "velocity_x"]
     assert_allclose(momentum_x / gas_density, velocity_x)
-        
+
     # Test radiation fields specifically
     assert ds.parameters["radiation_field_groups"] == 1
     assert "rad" in ds.fluid_types
     assert ("rad", "energy_density_0") in ds.derived_field_list
-        
+
     # Test particle IO with specific expectations
     assert ds.parameters["particle_types"] == ("Rad_particles",)
-        
+
     # Verify particle fields and units
     particle_info = ds.parameters["particle_info"]["Rad_particles"]
     assert particle_info["num_particles"] == 2

From 75bfe3f7f0375aa8f4aa0a7e66946e0ef1f2b4e4 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Sun, 18 May 2025 11:21:47 +1000
Subject: [PATCH 23/26] add tests of face-centered data

---
 yt/frontends/amrex/tests/test_outputs.py | 26 ++++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/yt/frontends/amrex/tests/test_outputs.py b/yt/frontends/amrex/tests/test_outputs.py
index 681b6025c61..1f78ef002a1 100644
--- a/yt/frontends/amrex/tests/test_outputs.py
+++ b/yt/frontends/amrex/tests/test_outputs.py
@@ -510,6 +510,32 @@ def test_quokka():
     assert particle_info["units"]["end_time"] == "T^1"
 
 
+quokka_face_centered = "HydroWave/plt00004"
+
+
+@requires_file(quokka_face_centered)
+def test_quokka_face_centered():
+    ds = data_dir_load(quokka_face_centered)
+
+    # Check if face-centered datasets were loaded
+    assert hasattr(ds, "ds_fc_x")
+    assert hasattr(ds, "ds_fc_y")
+    assert hasattr(ds, "ds_fc_z")
+
+    # Check that the face-centered datasets have the correct attributes
+    assert ds.ds_fc_x.fc_direction == "x"
+    assert ds.ds_fc_y.fc_direction == "y"
+    assert ds.ds_fc_z.fc_direction == "z"
+    
+    # Check that the face-centered datasets have a reference to the parent dataset
+    assert ds.ds_fc_x.parent_ds is ds
+    assert ds.ds_fc_y.parent_ds is ds
+    assert ds.ds_fc_z.parent_ds is ds
+    
+    assert isinstance(ds, QuokkaDataset)
+    assert isinstance(ds.ds_fc_x, QuokkaDataset)
+
+
 # test loading non-Cartesian coordinate systems in different dimensionalities
 
 

From dfc69a821a6f2dcfa65ec7fcfecee2d4d28ab55f Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Sun, 18 May 2025 12:14:50 +1000
Subject: [PATCH 24/26] rename dataset name

---
 yt/frontends/amrex/tests/test_outputs.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/yt/frontends/amrex/tests/test_outputs.py b/yt/frontends/amrex/tests/test_outputs.py
index 1f78ef002a1..fde0702e9c7 100644
--- a/yt/frontends/amrex/tests/test_outputs.py
+++ b/yt/frontends/amrex/tests/test_outputs.py
@@ -446,7 +446,7 @@ def test_maestro_parameters():
     assert type(ds.parameters["s0_interp_type"]) is int  # noqa: E721
 
 
-quokka = "RadiatingParticles_plt026"
+quokka = "quokka_RadiatingParticles_plt026"
 
 
 @requires_file(quokka)
@@ -510,7 +510,7 @@ def test_quokka():
     assert particle_info["units"]["end_time"] == "T^1"
 
 
-quokka_face_centered = "HydroWave/plt00004"
+quokka_face_centered = "quokka_HydroWave_plt00004"
 
 
 @requires_file(quokka_face_centered)

From d6799ce6fa02dcf45998a43bb3199964436ca2a2 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Sun, 18 May 2025 12:18:16 +1000
Subject: [PATCH 25/26] whitespace

---
 yt/frontends/amrex/tests/test_outputs.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/yt/frontends/amrex/tests/test_outputs.py b/yt/frontends/amrex/tests/test_outputs.py
index fde0702e9c7..79247295925 100644
--- a/yt/frontends/amrex/tests/test_outputs.py
+++ b/yt/frontends/amrex/tests/test_outputs.py
@@ -526,12 +526,12 @@ def test_quokka_face_centered():
     assert ds.ds_fc_x.fc_direction == "x"
     assert ds.ds_fc_y.fc_direction == "y"
     assert ds.ds_fc_z.fc_direction == "z"
-    
+
     # Check that the face-centered datasets have a reference to the parent dataset
     assert ds.ds_fc_x.parent_ds is ds
     assert ds.ds_fc_y.parent_ds is ds
     assert ds.ds_fc_z.parent_ds is ds
-    
+
     assert isinstance(ds, QuokkaDataset)
     assert isinstance(ds.ds_fc_x, QuokkaDataset)
 

From 581ad28a563e4789ac53b6138232c2307351deb0 Mon Sep 17 00:00:00 2001
From: ChongChong He <chongchong.he@anu.edu.au>
Date: Wed, 11 Jun 2025 00:26:01 +1000
Subject: [PATCH 26/26] make PyYAML a dependency of Quokka only

---
 pyproject.toml                        | 4 ++--
 yt/frontends/amrex/data_structures.py | 6 +++++-
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 03300c8360f..16a5f9b41ee 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -89,7 +89,7 @@ Fortran = ["f90nml>=1.1"]
 # We also normalize all target names to lower case for consistency.
 adaptahop = []
 ahf = []
-amrex = ["PyYAML>=6.0.1"]
+amrex = []
 amrvac = ["yt[Fortran]"]
 art = []
 arepo = ["yt[HDF5]"]
@@ -121,7 +121,7 @@ open-pmd = ["yt[HDF5]"]
 owls = ["yt[HDF5]"]
 owls-subfind = ["yt[HDF5]"]
 parthenon = ["yt[HDF5]"]
-quokka = ["PyYAML>=6.0.1"]
+quokka = ["yt[amrex]", "PyYAML>=6.0.1"]
 ramses = ["yt[Fortran]", "scipy"]
 rockstar = []
 sdf = ["requests>=2.20.0"]
diff --git a/yt/frontends/amrex/data_structures.py b/yt/frontends/amrex/data_structures.py
index 7a27f4efbf4..4e584de6b11 100644
--- a/yt/frontends/amrex/data_structures.py
+++ b/yt/frontends/amrex/data_structures.py
@@ -6,7 +6,6 @@
 from stat import ST_CTIME
 
 import numpy as np
-import yaml
 
 from yt.data_objects.index_subobjects.grid_patch import AMRGridPatch
 from yt.data_objects.static_output import Dataset
@@ -1373,6 +1372,7 @@ class QuokkaDataset(AMReXDataset):
     _field_info_class = QuokkaFieldInfo
     _subtype_keyword = ""
     _default_cparam_filename = "metadata.yaml"
+    _load_requirements = ["yaml"]
 
     def __init__(
         self,
@@ -1465,6 +1465,8 @@ def _load_face_centered_datasets(self):
                 mylog.debug(f"No face-centered {direction} dataset found at {fc_dir}")
 
     def _parse_parameter_file(self):
+        import yaml
+
         # Call parent method to initialize core setup by yt
         super()._parse_parameter_file()
 
@@ -1686,6 +1688,8 @@ def _parse_parameter_file(self):
         mylog.debug("Parsed header parameters: %s", self.parameters)
 
     def _parse_metadata_file(self):
+        import yaml
+
         # Construct the full path to the metadata file
         metadata_filename = os.path.join(self.output_dir, self.cparam_filename)
         try: