Added postprocess() and run() methods to Benchmark class (#58)

* Benchmark build methods -> private

* + _build_model() abstractmethod

* + _build_model abstractmethod

* + _build_model in OpenmcBenchmark

* + self.model in OpenmcBenchmark constructor

* Benchmark build methods -> private

* + _build_model() abstractmethod

* + _build_model abstractmethod

* + _build_model in OpenmcBenchmark

* + self.model in OpenmcBenchmark constructor

* fix tests

* fix test

* + abstractmethod postprocess

* okt specs tally surf id -> 7

* + self._build_source() fix in build settings

* _build_geometry check if mesh.h5m is already present

* + warning if mesh.h5m already present

* update todo.md

* + xarray as dependency

* + tally postprocess, df -> xarray

* updated todo file

* + store results in benchmark_results.h5

* update todo

* + run() abstractmethod

* + run() method in OpenmcBenchmark

* + throw error if unsupported run mode

* update todo

* + pydagmc dependency

* - pydagmc dependency

* + fix polar angle -> cos(angle) source

* + pydagmc as dependency

* fix pydagmc dependency

* update tests

* + fix tests

* typo in surfaces_by_id

* updated todo

* drop 'surface' col in result df

* mode='a' in to_netcdf result method

* add vectorial logo

* change logo

* added 3rd dimension to result xarray

* + small fix on xarray definition

* Experimental configuration of OKTAVIAN experiments

* OKTAVIAN general information

* references

* add al results

* + point to rignt lfs submodule branch

* specify branch in submodule

* + lfs change

* test cad visualization

* test for cad visualization

* update name for oktavian al benchmark

* deleted wrong file

* change src

* change references

* add .glb file for cad visualization

* add html page

* add new reference style

* add interactive view

* + abstractmethod postprocess

* okt specs tally surf id -> 7

* + self._build_source() fix in build settings

* _build_geometry check if mesh.h5m is already present

* + warning if mesh.h5m already present

* update todo.md

* + xarray as dependency

* + tally postprocess, df -> xarray

* updated todo file

* + store results in benchmark_results.h5

* update todo

* + run() abstractmethod

* + run() method in OpenmcBenchmark

* + throw error if unsupported run mode

* update todo

* + pydagmc dependency

* - pydagmc dependency

* + fix polar angle -> cos(angle) source

* + pydagmc as dependency

* fix pydagmc dependency

* update tests

* + fix tests

* typo in surfaces_by_id

* updated todo

* drop 'surface' col in result df

* mode='a' in to_netcdf result method

* added 3rd dimension to result xarray

* + small fix on xarray definition

---------

Co-authored-by: davidepettinari <[email protected]>
Co-authored-by: Davide Pettinari <[email protected]>
Co-authored-by: Ethan Peterson <[email protected]>

Author: 3 people

TODO.md

@@ -1,16 +1,31 @@
 # To do
 ### Manage CAD geometry
 
-### Benchmark modules
+
+## Benchmarks API
+- Add capability of merging multiple surface results for oktavian (should happen only with surface tallies)
+- Add expected tally results shape in `schema` and `specifications` to deal with the point here above
+- Try to manage results as numpy arrays instead of pandas dataframes in `OpenmcBenchmark.postprocess()`
+- Should we use `h5netcdf` engine for saving/opening datasets in `h5` files or stick with the default `netcdf4`?
+- Do we wanna use `makefun` decorator for wrapped methods (e.g. `OpenmcBenchmark.run()`)?
 
 ## Tests
 - Add tests to everything
 
 ## Notebooks
+- Use `specifications` validation against `schema`  function and script
+- Open and inspect a `specifications` file (just `metadata`)
+- Instantiate an `OpenmcBenchmark` object, build model, run simulation
+- Postprocessing and visualization
 
-## clean postprocessing
+## Documentation
+- introduction and motivation
+- V&V explanation
+- Definition of benchmark `specifications`
+- Benchmarks
+    - Oktavian description
+    - Oktavian results
 
 ## Miscellanea
 
-
 ## Questions

pyproject.toml

@@ -22,8 +22,10 @@ dependencies = [
     "openapi_schema_validator",
     "pydantic",
     "numpy",
+    "xarray",
     "tables",
     "cad-to-dagmc",
+    "pydagmc @ git+https://github.com/svalinn/pydagmc.git"
 ]
 
 [tool.setuptools.package-data]

src/openmc_fusion_benchmarks/benchmark.py

@@ -1,10 +1,14 @@
 import yaml
 from pathlib import Path
+import warnings
 from abc import ABC, abstractmethod
 import numpy as np
+import xarray as xr
+import h5py
 from .validate import validate_benchmark
 
 import openmc
+import pydagmc
 from cad_to_dagmc import CadToDagmc
 
 
@@ -54,6 +58,16 @@ def _build_model(self):
         """Build the whole model for the benchmark."""
         pass
 
+    @abstractmethod
+    def postprocess(self):
+        """Post-process the model after running."""
+        pass
+
+    @abstractmethod
+    def run(self):
+        """Run the benchmark simulation."""
+        pass
+
     def _read_metadata(self):
         """Read metadata from the benchmark specification."""
         metadata = self._benchmark_spec['metadata']
@@ -175,9 +189,13 @@ def build_mesh(cad_file: str, material_tags, set_size: dict, global_mesh_size_mi
         # Get the STEP file
         cad_file = LFS_DIR / "benchmarks" / f"{geometry_data['cad_file']}"
 
-        # Generate the mesh
-        build_mesh(cad_file=cad_file, material_tags=material_tags, set_size=set_size,
-                   global_mesh_size_min=global_mesh_size_min, global_mesh_size_max=global_mesh_size_max)
+        # Generate the mesh if mesh.h5m not already present
+        if Path("mesh.h5m").exists():
+            warnings.warn(
+                f"Mesh file already exists. Skipping mesh generation.")
+        else:
+            build_mesh(cad_file=cad_file, material_tags=material_tags, set_size=set_size,
+                       global_mesh_size_min=global_mesh_size_min, global_mesh_size_max=global_mesh_size_max)
 
         # download the h5m file
         # download_from_drive(benchmark_name=self.name, file_format='h5m')
@@ -251,7 +269,7 @@ def angular_conversion(values, units):
             angular_sources = []
             # Openmc needs one source per angle bin:
             angles = source['angular_energy_distribution']
-            abins = np.array(angles['angle']['bins'])
+            abins = np.cos(angles['angle']['bins'])
             for i in range(len(abins)-1):
                 lb = angular_conversion(
                     abins[i], angles['angle']['units'])
@@ -344,10 +362,13 @@ def _build_settings(self):
         settings_data = self._benchmark_spec['settings']
 
         settings = openmc.Settings()
-        if settings_data['run_mode'] == 'fixed source':
+        if settings_data['run_mode'] == 'fixed_source':
             settings.run_mode = 'fixed source'
         elif settings_data['run_mode'] == 'k-eigenvalue':
             settings.run_mode = 'eigenvalue'
+        else:
+            raise ValueError(
+                f"Unsupported run mode: {settings_data['run_mode']}")
         settings.batches = int(settings_data['batches'])
         settings.particles = int(settings_data['particles_per_batch'])
         settings.photon_transport = settings_data['photon_transport']
@@ -356,7 +377,7 @@ def _build_settings(self):
         # electron treatment
         settings.output = {'tallies': False}
 
-        source = self.build_source()
+        source = self._build_source()
         settings.source = source
 
         return settings
@@ -373,3 +394,78 @@ def _build_model(self):
             tallies=tallies
         )
         return model
+
+    def postprocess(self, statepoint: str = 'statepoint.100.h5', mesh: str = 'mesh.h5m'):
+        """Post-process the model after running."""
+        # Retrieve tallies data from specifications
+        tallies_data = self._benchmark_spec['tallies']
+        # Read openmc statepoint file
+        sp = openmc.StatePoint(statepoint)
+        # Read mesh file
+        mesh = pydagmc.DAGModel(mesh)
+
+        # Cycle tallies in specifications
+        for spec_t in tallies_data:
+            # Get corresponding tally from statepoint
+            df = sp.get_tally(name=spec_t['name']).get_pandas_dataframe()
+
+            # Preparing tally dataframe
+            df = df.drop(columns=['surface', 'cell', 'particle', 'nuclide',
+                                  'score', 'energyfunction'], errors='ignore')
+            # Cyle tally filters
+            norm = 1
+            for f in spec_t['filters']:
+                if f['type'] == 'cell':
+                    # Get cell volumes for normalization
+                    norm = [mesh.volumes_by_id[v].area for v in f['values']]
+                elif f['type'] == 'surface':
+                    # Get surface areas for normalization
+                    norm = [mesh.surfaces_by_id[v].area for v in f['values']]
+                elif f['type'] == 'material':
+                    raise NotImplementedError(
+                        'Material filter not implemented in postprocess yet.')
+
+                # Normalize the tally data
+                df['mean'] = df['mean'] / norm
+                df['std. dev.'] = df['std. dev.'] / norm
+
+            # Convert to xarray and add dimensions
+            t = xr.DataArray(
+                df.values[np.newaxis, :, :],  # shape: (1, r, c)
+                dims=["case", "row", "column"],
+                coords={
+                    "case": ["0"],
+                    "column": df.columns,
+                    "row": np.arange(df.shape[0])
+                },
+                name=spec_t['name']
+            )
+
+            # Save the tally data to a netCDF file
+            t.to_netcdf(f"benchmark_results.h5",
+                        group=f"{spec_t['name']}", mode='a')
+
+        # Add some metadata attributes
+        with h5py.File(f"benchmark_results.h5", "a") as f:
+            f.attrs['benchmark_name'] = self.name
+            # f.attrs['benchmark_version'] = self._benchmark_spec['metadata'].get(
+            #     'version', 'N/A')
+            # f.attrs['description'] = self._benchmark_spec['metadata'].get(
+            #     'description', 'N/A')
+            # f.attrs['title'] = self._benchmark_spec['metadata'].get(
+            #     'title', 'N/A')
+            # f.attrs['literature'] = self._benchmark_spec['metadata'].get(
+            #     'references', [])
+            f.attrs['code'] = f"openmc {sp.version[0]}.{sp.version[1]}.{sp.version[2]}"
+
+        return
+
+    def run(self, *args, **kwargs):
+        """Run the benchmark simulation."""
+        # Run the OpenMC model
+        statepoint = self.model.run(*args, **kwargs)
+
+        # Post-process the results
+        self.postprocess(statepoint=statepoint)
+
+        return

src/openmc_fusion_benchmarks/benchmarks/oktavian_al/specifications.yaml

@@ -146,7 +146,7 @@ tallies:
     scores: [current]
     filters:
       - type: surface
-        values: [6]
+        values: [7]
       - type: energy
         values: [97122., 101090., 105210., 109500., 113970., 118620.,
                 123470., 128500., 133750., 139210., 144890., 150800.,
@@ -179,7 +179,7 @@ tallies:
     scores: [current]
     filters:
       - type: surface
-        values: [6]
+        values: [7]
       - type: energy
         values: [500000.,   600000.,   700000.,   800000.,   900000.,  1000000.,
                 1100000.,  1200000.,  1300000.,  1400000.,  1500000.,  1600000.,

test/test_benchmark.py

@@ -14,6 +14,8 @@ def _build_source(self): return "source"
     def _build_settings(self): return "settings"
     def _build_tallies(self): return "tallies"
     def _build_model(self): return "model"
+    def postprocess(self): return "postprocess"
+    def run(self): return "run"
 
 
 @pytest.fixture
@@ -23,7 +25,7 @@ def valid_yaml():
         "materials": [],
         "geometry": {},
         "sources": [],
-        "tallies": []
+        "tallies": [],
     })