Merge branch 'develop' into virtual_lattice_0.15.2

Author: skywalker_cn

cmake/OpenMCConfig.cmake.in

@@ -1,9 +1,12 @@
 get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY)
 
-find_package(fmt REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../fmt)
-find_package(pugixml REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../pugixml)
-find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl)
-find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor)
+# Compute the install prefix from this file's location
+get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE)
+
+find_package(fmt CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
+find_package(pugixml CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
+find_package(xtl CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
+find_package(xtensor CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
 if(@OPENMC_USE_DAGMC@)
   find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@)
 endif()

docs/source/usersguide/kinetics.rst

@@ -67,6 +67,23 @@ are needed to compute kinetics parameters in OpenMC:
 Obtaining kinetics parameters
 -----------------------------
 
+The ``Model`` class can be used to automatically generate all IFP tallies using
+the Python API with :attr:`openmc.Settings.ifp_n_generation` greater than 0 and
+the :meth:`openmc.Model.add_ifp_kinetics_tallies` method::
+
+    model = openmc.Model(geometry, settings=settings)
+    model.add_kinetics_parameters_tallies(num_groups=6)  # Add 6 precursor groups
+
+Alternatively, each of the tallies can be manually defined using group-wise or
+total :math:`\beta_{\text{eff}}` specified by providing a 6-group
+:class:`openmc.DelayedGroupFilter`::
+
+    beta_tally = openmc.Tally(name="group-beta-score")
+    beta_tally.scores = ["ifp-beta-numerator"]
+
+    # Add DelayedGroupFilter to enable group-wise tallies
+    beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, 7)))]
+
 Here is an example showing how to declare the three available IFP scores in a
 single tally::
 
@@ -95,6 +112,12 @@ for ``ifp-denominator``:
 
     \beta_{\text{eff}} = \frac{S_{\text{ifp-beta-numerator}}}{S_{\text{ifp-denominator}}}
 
+The kinetics parameters can be retrieved directly from a statepoint file using
+the :meth:`openmc.StatePoint.ifp_results` method::
+
+    with openmc.StatePoint(output_path) as sp:
+        generation_time, beta_eff = sp.get_kinetics_parameters()
+
 .. only:: html
 
    .. rubric:: References
@@ -107,4 +130,4 @@ for ``ifp-denominator``:
     of the Iterated Fission Probability Method in OpenMC to Compute Adjoint-Weighted
     Kinetics Parameters", International Conference on Mathematics and Computational
     Methods Applied to Nuclear Science and Engineering (M&C 2025), Denver, April 27-30,
-    2025 (to be presented).
+    2025.

include/openmc/ifp.h

@@ -68,15 +68,14 @@ vector<T> _ifp(const T& value, const vector<T>& data)
 //!
 //! Add the IFP information in the IFP banks using the same index
 //! as the one used to append the fission site to the fission bank.
+//! The information stored are the delayed group number and lifetime
+//! of the neutron that created the fission event.
 //! Multithreading protection is guaranteed by the index returned by the
 //! thread_safe_append call in physics.cpp.
 //!
-//! Needs to be done after the delayed group is found.
-//!
 //! \param[in] p Particle
-//! \param[in] site Fission site
 //! \param[in] idx Bank index from the thread_safe_append call in physics.cpp
-void ifp(const Particle& p, const SourceSite& site, int64_t idx);
+void ifp(const Particle& p, int64_t idx);
 
 //! Resize the IFP banks used in the simulation
 void resize_simulation_ifp_banks();

include/openmc/mesh.h

@@ -990,25 +990,26 @@ class LibMesh : public UnstructuredMesh {
 
   libMesh::MeshBase* mesh_ptr() const { return m_; };
 
-private:
-  void initialize() override;
-  void set_mesh_pointer_from_filename(const std::string& filename);
-  void build_eqn_sys();
-
+protected:
   // Methods
 
   //! Translate a bin value to an element reference
-  const libMesh::Elem& get_element_from_bin(int bin) const;
+  virtual const libMesh::Elem& get_element_from_bin(int bin) const;
 
   //! Translate an element pointer to a bin index
-  int get_bin_from_element(const libMesh::Elem* elem) const;
+  virtual int get_bin_from_element(const libMesh::Elem* elem) const;
+
+  libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set
+                         //!< during intialization
+private:
+  void initialize() override;
+  void set_mesh_pointer_from_filename(const std::string& filename);
+  void build_eqn_sys();
 
   // Data members
   unique_ptr<libMesh::MeshBase> unique_m_ =
     nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is
              //!< created inside OpenMC
-  libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set
-                         //!< during intialization
   vector<unique_ptr<libMesh::PointLocatorBase>>
     pl_; //!< per-thread point locators
   unique_ptr<libMesh::EquationSystems>
@@ -1022,8 +1023,34 @@ class LibMesh : public UnstructuredMesh {
   libMesh::BoundingBox bbox_; //!< bounding box of the mesh
   libMesh::dof_id_type
     first_element_id_; //!< id of the first element in the mesh
+};
+
+class AdaptiveLibMesh : public LibMesh {
+public:
+  // Constructor
+  AdaptiveLibMesh(
+    libMesh::MeshBase& input_mesh, double length_multiplier = 1.0);
+
+  // Overridden methods
+  int n_bins() const override;
+
+  void add_score(const std::string& var_name) override;
+
+  void set_score_data(const std::string& var_name, const vector<double>& values,
+    const vector<double>& std_dev) override;
+
+  void write(const std::string& filename) const override;
+
+protected:
+  // Overridden methods
+  int get_bin_from_element(const libMesh::Elem* elem) const override;
+
+  const libMesh::Elem& get_element_from_bin(int bin) const override;
+
+private:
+  // Data members
+  const libMesh::dof_id_type num_active_; //!< cached number of active elements
 
-  const bool adaptive_; //!< whether this mesh has adaptivity enabled or not
   std::vector<libMesh::dof_id_type>
     bin_to_elem_map_; //!< mapping bin indices to dof indices for active
                       //!< elements

include/openmc/particle_data.h

@@ -662,6 +662,7 @@ class ParticleData : public GeometryState {
   int& event_mt() { return event_mt_; } // MT number of collision
   const int& event_mt() const { return event_mt_; }
   int& delayed_group() { return delayed_group_; } // delayed group
+  const int& delayed_group() const { return delayed_group_; }
   const int& parent_nuclide() const { return parent_nuclide_; }
   int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide
 

openmc/deplete/transfer_rates.py

@@ -366,7 +366,7 @@ def set_external_source_rate(
         rate : float
             External source rate in units of mass per time. A positive or
             negative value corresponds to a feed or removal rate, respectively.
-        units : {'g/s', 'g/min', 'g/h', 'g/d', 'g/a'}
+        rate_units : {'g/s', 'g/min', 'g/h', 'g/d', 'g/a'}
             Units for values specified in the `rate` argument. 's' for seconds,
             'min' for minutes, 'h' for hours, 'a' for Julian years.
         timesteps : list of int, optional

openmc/mesh.py

@@ -287,6 +287,7 @@ def get_homogenized_materials(
             model: openmc.Model,
             n_samples: int | tuple[int, int, int] = 10_000,
             include_void: bool = True,
+            material_volumes: MeshMaterialVolumes | None = None,
             **kwargs
     ) -> list[openmc.Material]:
         """Generate homogenized materials over each element in a mesh.
@@ -305,16 +306,23 @@ def get_homogenized_materials(
             the x, y, and z dimensions.
         include_void : bool, optional
             Whether homogenization should include voids.
+        material_volumes : MeshMaterialVolumes, optional
+            Previously computed mesh material volumes to use for homogenization.
+            If not provided, they will be computed by calling
+            :meth:`material_volumes`.
         **kwargs
-            Keyword-arguments passed to :meth:`MeshBase.material_volumes`.
+            Keyword-arguments passed to :meth:`material_volumes`.
 
         Returns
         -------
         list of openmc.Material
             Homogenized material in each mesh element
 
         """
-        vols = self.material_volumes(model, n_samples, **kwargs)
+        if material_volumes is None:
+            vols = self.material_volumes(model, n_samples, **kwargs)
+        else:
+            vols = material_volumes
         mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)]
 
         # Create homogenized material for each element
@@ -424,7 +432,6 @@ def material_volumes(
 
         # Restore original tallies
         model.tallies = original_tallies
-
         return volumes
 
 

openmc/model/model.py

@@ -1,7 +1,7 @@
 from __future__ import annotations
 from collections.abc import Iterable, Sequence
 import copy
-from functools import lru_cache
+from functools import cache
 from pathlib import Path
 import math
 from numbers import Integral, Real
@@ -160,7 +160,7 @@ def is_initialized(self) -> bool:
             return False
 
     @property
-    @lru_cache(maxsize=None)
+    @cache
     def _materials_by_id(self) -> dict:
         """Dictionary mapping material ID --> material"""
         if self.materials:
@@ -170,14 +170,14 @@ def _materials_by_id(self) -> dict:
         return {mat.id: mat for mat in mats}
 
     @property
-    @lru_cache(maxsize=None)
+    @cache
     def _cells_by_id(self) -> dict:
         """Dictionary mapping cell ID --> cell"""
         cells = self.geometry.get_all_cells()
         return {cell.id: cell for cell in cells.values()}
 
     @property
-    @lru_cache(maxsize=None)
+    @cache
     def _cells_by_name(self) -> dict[int, openmc.Cell]:
         # Get the names maps, but since names are not unique, store a set for
         # each name key. In this way when the user requests a change by a name,
@@ -190,7 +190,7 @@ def _cells_by_name(self) -> dict[int, openmc.Cell]:
         return result
 
     @property
-    @lru_cache(maxsize=None)
+    @cache
     def _materials_by_name(self) -> dict[int, openmc.Material]:
         if self.materials is None:
             mats = self.geometry.get_all_materials().values()
@@ -203,6 +203,37 @@ def _materials_by_name(self) -> dict[int, openmc.Material]:
             result[mat.name].add(mat)
         return result
 
+    def add_kinetics_parameters_tallies(self, num_groups: int | None = None):
+        """Add tallies for calculating kinetics parameters using the IFP method.
+
+        This method adds tallies to the model for calculating two kinetics
+        parameters, the generation time and the effective delayed neutron
+        fraction (beta effective). After a model is run, these parameters can be
+        determined through the :meth:`openmc.StatePoint.ifp_results` method.
+
+        Parameters
+        ----------
+        num_groups : int, optional
+            Number of precursor groups to filter the delayed neutron fraction.
+            If None, only the total effective delayed neutron fraction is
+            tallied.
+
+        """
+        if not any('ifp-time-numerator' in t.scores for t in self.tallies):
+            gen_time_tally = openmc.Tally(name='IFP time numerator')
+            gen_time_tally.scores = ['ifp-time-numerator']
+            self.tallies.append(gen_time_tally)
+        if not any('ifp-beta-numerator' in t.scores for t in self.tallies):
+            beta_tally = openmc.Tally(name='IFP beta numerator')
+            beta_tally.scores = ['ifp-beta-numerator']
+            if num_groups is not None:
+                beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))]
+            self.tallies.append(beta_tally)
+        if not any('ifp-denominator' in t.scores for t in self.tallies):
+            denom_tally = openmc.Tally(name='IFP denominator')
+            denom_tally.scores = ['ifp-denominator']
+            self.tallies.append(denom_tally)
+
     @classmethod
     def from_xml(
         cls,

openmc/settings.py

@@ -1714,12 +1714,18 @@ def _create_random_ray_subelement(self, root, mesh_memo=None):
                             domain_elem = ET.SubElement(mesh_elem, 'domain')
                             domain_elem.set('id', str(domain.id))
                             domain_elem.set('type', domain.__class__.__name__.lower())
-                        if mesh_memo is not None and mesh.id not in mesh_memo:
+                        # See if a <mesh> element already exists -- if not, add it
+                        path = f"./mesh[@id='{mesh.id}']"
+                        if root.find(path) is None:
                             root.append(mesh.to_xml_element())
-                            mesh_memo.add(mesh.id)
+                            if mesh_memo is not None:
+                                mesh_memo.add(mesh.id)
                 else:
                     subelement = ET.SubElement(element, key)
-                    subelement.text = str(value)
+                    if isinstance(value, bool):
+                        subelement.text = str(value).lower()
+                    else:
+                        subelement.text = str(value)
 
     def _create_source_rejection_fraction_subelement(self, root):
         if self._source_rejection_fraction is not None:
@@ -2103,7 +2109,7 @@ def _max_tracks_from_xml_element(self, root):
         if text is not None:
             self.max_tracks = int(text)
 
-    def _random_ray_from_xml_element(self, root):
+    def _random_ray_from_xml_element(self, root, meshes=None):
         elem = root.find('random_ray')
         if elem is not None:
             self.random_ray = {}
@@ -2130,7 +2136,11 @@ def _random_ray_from_xml_element(self, root):
                 elif child.tag == 'source_region_meshes':
                     self.random_ray['source_region_meshes'] = []
                     for mesh_elem in child.findall('mesh'):
-                        mesh = MeshBase.from_xml_element(mesh_elem)
+                        mesh_id = int(get_text(mesh_elem, 'id'))
+                        if meshes and mesh_id in meshes:
+                            mesh = meshes[mesh_id]
+                        else:
+                            mesh = MeshBase.from_xml_element(mesh_elem)
                         domains = []
                         for domain_elem in mesh_elem.findall('domain'):
                             domain_id = int(get_text(domain_elem, "id"))
@@ -2329,7 +2339,7 @@ def from_xml_element(cls, elem, meshes=None):
         settings._max_history_splits_from_xml_element(elem)
         settings._max_tracks_from_xml_element(elem)
         settings._max_secondaries_from_xml_element(elem)
-        settings._random_ray_from_xml_element(elem)
+        settings._random_ray_from_xml_element(elem, meshes)
         settings._use_decay_photons_from_xml_element(elem)
         settings._source_rejection_fraction_from_xml_element(elem)