Merge remote-tracking branch 'upstream/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/capi/index.rst

@@ -84,6 +84,17 @@ Functions
    :return: Return status (negative if an error occurred)
    :rtype: int
 
+.. c:function:: int openmc_cell_get_density(int32_t index, const int32_t* instance, double* density)
+
+   Get the density of a cell
+
+   :param int32_t index: Index in the cells array
+   :param int32_t* instance: Which instance of the cell. If a null pointer is passed, the density
+                             multiplier of the first instance is returned.
+   :param double* density: Density of the cell in [g/cm3]
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
 .. c:function:: int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices)
 
    Set the fill for a cell
@@ -113,8 +124,22 @@ Functions
    :param double T: Temperature in Kelvin
    :param instance: Which instance of the cell. To set the temperature for all
                     instances, pass a null pointer.
-   :param set_contained: If the cell is not filled by a material, whether to set the temperatures
-                         of all filled cells
+   :param bool set_contained: If the cell is not filled by a material, whether
+                              to set the temperatures of all filled cells
+   :type instance: const int32_t*
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_cell_set_density(index index, double density, const int32_t* instance, bool set_contained)
+
+   Set the density of a cell.
+
+   :param int32_t index: Index in the cells array
+   :param double density: Density of the cell in [g/cm3]
+   :param instance: Which instance of the cell. To set the density multiplier for all
+                    instances, pass a null pointer.
+   :param bool set_contained: If the cell is not filled by a material, whether
+                              to set the density multiplier of all filled cells
    :type instance: const int32_t*
    :return: Return status (negative if an error occurred)
    :rtype: int

docs/source/conf.py

@@ -121,9 +121,7 @@
 # -- Options for HTML output ---------------------------------------------------
 
 # The theme to use for HTML and HTML Help pages
-import sphinx_rtd_theme
 html_theme = 'sphinx_rtd_theme'
-html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 html_baseurl = "https://docs.openmc.org/en/stable/"
 
 html_logo = '_images/openmc_logo.png'

docs/source/io_formats/properties.rst

@@ -4,7 +4,7 @@
 Properties File Format
 ======================
 
-The current version of the properties file format is 1.0.
+The current version of the properties file format is 1.1.
 
 **/**
 
@@ -25,6 +25,7 @@ The current version of the properties file format is 1.0.
 **/geometry/cells/cell <uid>/**
 
 :Datasets: - **temperature** (*double[]*) -- Temperature of the cell in [K].
+           - **density** (*double[]*) -- Density of the cell in [g/cm3].
 
 **/materials/**
 

docs/source/io_formats/summary.rst

@@ -4,7 +4,7 @@
 Summary File Format
 ===================
 
-The current version of the summary file format is 6.0.
+The current version of the summary file format is 6.1.
 
 **/**
 
@@ -38,6 +38,7 @@ The current version of the summary file format is 6.0.
              is an array if the cell uses distributed materials, otherwise it is
              a scalar.
            - **temperature** (*double[]*) -- Temperature of the cell in Kelvin.
+           - **density** (*double[]*) -- Density of the cell in [g/cm3].
            - **translation** (*double[3]*) -- Translation applied to the fill
              universe. This dataset is present only if fill_type is set to
              'universe'.

docs/source/methods/neutron_physics.rst

@@ -290,7 +290,10 @@ create and store fission sites for the following generation. First, the average
 number of prompt and delayed neutrons must be determined to decide whether the
 secondary neutrons will be prompt or delayed. This is important because delayed
 neutrons have a markedly different spectrum from prompt neutrons, one that has a
-lower average energy of emission. The total number of neutrons emitted
+lower average energy of emission. Furthermore, in simulations where tracking
+time of neutrons is important, we need to consider the emission time delay of
+the secondary neutrons, which is dependent on the decay constant of the
+delayed neutron precursor. The total number of neutrons emitted
 :math:`\nu_t` is given as a function of incident energy in the ENDF format. Two
 representations exist for :math:`\nu_t`. The first is a polynomial of order
 :math:`N` with coefficients :math:`c_0,c_1,\dots,c_N`. If :math:`\nu_t` has this
@@ -306,8 +309,8 @@ interpolation law. The number of prompt neutrons released per fission event
 :math:`\nu_p` is also given as a function of incident energy and can be
 specified in a polynomial or tabular format. The number of delayed neutrons
 released per fission event :math:`\nu_d` can only be specified in a tabular
-format. In practice, we only need to determine :math:`nu_t` and
-:math:`nu_d`. Once these have been determined, we can calculated the delayed
+format. In practice, we only need to determine :math:`\nu_t` and
+:math:`\nu_d`. Once these have been determined, we can calculate the delayed
 neutron fraction
 
 .. math::
@@ -335,8 +338,14 @@ neutrons. Otherwise, we produce :math:`\lfloor \nu \rfloor + 1` neutrons. Then,
 for each fission site produced, we sample the outgoing angle and energy
 according to the algorithms given in :ref:`sample-angle` and
 :ref:`sample-energy` respectively. If the neutron is to be born delayed, then
-there is an extra step of sampling a delayed neutron precursor group since they
-each have an associated secondary energy distribution.
+there is an extra step of sampling a delayed neutron precursor group to get the
+associated secondary energy distribution and the decay constant
+:math:`\lambda`, which is needed to sample the emission delay time :math:`t_d`:
+
+.. math::
+    :label: sample-delay-time
+
+    t_d = -\frac{\ln \xi}{\lambda}.
 
 The sampled outgoing angle and energy of fission neutrons along with the
 position of the collision site are stored in an array called the fission

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.

examples/pincell_pulsed/run_pulse.py

@@ -0,0 +1,101 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import openmc
+
+###############################################################################
+# Create materials for the problem
+
+uo2 = openmc.Material(name="UO2 fuel at 2.4% wt enrichment")
+uo2.set_density("g/cm3", 10.29769)
+uo2.add_element("U", 1.0, enrichment=2.4)
+uo2.add_element("O", 2.0)
+
+helium = openmc.Material(name="Helium for gap")
+helium.set_density("g/cm3", 0.001598)
+helium.add_element("He", 2.4044e-4)
+
+zircaloy = openmc.Material(name="Zircaloy 4")
+zircaloy.set_density("g/cm3", 6.55)
+zircaloy.add_element("Sn", 0.014, "wo")
+zircaloy.add_element("Fe", 0.00165, "wo")
+zircaloy.add_element("Cr", 0.001, "wo")
+zircaloy.add_element("Zr", 0.98335, "wo")
+
+borated_water = openmc.Material(name="Borated water")
+borated_water.set_density("g/cm3", 0.740582)
+borated_water.add_element("B", 2.0e-4)  # 3x the original pincell
+borated_water.add_element("H", 5.0e-2)
+borated_water.add_element("O", 2.4e-2)
+borated_water.add_s_alpha_beta("c_H_in_H2O")
+
+###############################################################################
+# Define problem geometry
+
+# Create cylindrical surfaces
+fuel_or = openmc.ZCylinder(r=0.39218, name="Fuel OR")
+clad_ir = openmc.ZCylinder(r=0.40005, name="Clad IR")
+clad_or = openmc.ZCylinder(r=0.45720, name="Clad OR")
+
+# Create a region represented as the inside of a rectangular prism
+pitch = 1.25984
+box = openmc.model.RectangularPrism(pitch, pitch, boundary_type="reflective")
+
+# Create cells, mapping materials to regions
+fuel = openmc.Cell(fill=uo2, region=-fuel_or)
+gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
+clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
+water = openmc.Cell(fill=borated_water, region=+clad_or & -box)
+
+# Create a model and assign geometry
+model = openmc.Model()
+model.geometry = openmc.Geometry([fuel, gap, clad, water])
+
+###############################################################################
+# Define problem settings
+
+# Set the mode
+model.settings.run_mode = "fixed source"
+
+# Indicate how many batches and particles to run
+model.settings.batches = 10
+model.settings.particles = 10000
+
+# Set time cutoff (we only care about t < 100 seconds, see tally below)
+model.settings.cutoff = {"time_neutron": 100}
+
+# Create the neutron pulse source (by default, isotropic direction, t=0)
+space = openmc.stats.Point()  # At the origin (0, 0, 0)
+energy = openmc.stats.delta_function(14.1e6)  # At 14.1 MeV
+model.settings.source = openmc.IndependentSource(space=space, energy=energy)
+
+###############################################################################
+# Define tallies
+
+# Create time filter
+t_grid = np.insert(np.logspace(-6, 2, 100), 0, 0.0)
+time_filter = openmc.TimeFilter(t_grid)
+
+# Tally for total neutron density in time
+density_tally = openmc.Tally(name="Density")
+density_tally.filters = [time_filter]
+density_tally.scores = ["inverse-velocity"]
+
+# Add tallies to model
+model.tallies = openmc.Tallies([density_tally])
+
+
+# Run the model
+model.run(apply_tally_results=True)
+
+# Bin-averaged result
+density_mean = density_tally.mean.ravel() / np.diff(t_grid)
+
+# Plot particle density versus time
+fig, ax = plt.subplots()
+ax.stairs(density_mean, t_grid)
+ax.set_xscale("log")
+ax.set_yscale("log")
+ax.set_xlabel("Time [s]")
+ax.set_ylabel("Total density")
+ax.grid()
+plt.show()

include/openmc/capi.h

@@ -17,6 +17,8 @@ int openmc_cell_get_fill(
 int openmc_cell_get_id(int32_t index, int32_t* id);
 int openmc_cell_get_temperature(
   int32_t index, const int32_t* instance, double* T);
+int openmc_cell_get_density(
+  int32_t index, const int32_t* instance, double* rho);
 int openmc_cell_get_translation(int32_t index, double xyz[]);
 int openmc_cell_get_rotation(int32_t index, double rot[], size_t* n);
 int openmc_cell_get_name(int32_t index, const char** name);
@@ -27,6 +29,8 @@ int openmc_cell_set_fill(
 int openmc_cell_set_id(int32_t index, int32_t id);
 int openmc_cell_set_temperature(
   int32_t index, double T, const int32_t* instance, bool set_contained = false);
+int openmc_cell_set_density(int32_t index, double rho, const int32_t* instance,
+  bool set_contained = false);
 int openmc_cell_set_translation(int32_t index, const double xyz[]);
 int openmc_cell_set_rotation(int32_t index, const double rot[], size_t rot_len);
 int openmc_dagmc_universe_get_cell_ids(

include/openmc/cell.h

@@ -216,6 +216,18 @@ class Cell {
   //! \return Temperature in [K]
   double temperature(int32_t instance = -1) const;
 
+  //! Get the density multiplier of a cell instance
+  //! \param[in] instance Instance index. If -1 is given, the density multiplier
+  //! for the first instance is returned.
+  //! \return Density multiplier
+  double density_mult(int32_t instance = -1) const;
+
+  //! Get the density of a cell instance in g/cm3
+  //! \param[in] instance Instance index. If -1 is given, the density
+  //! for the first instance is returned.
+  //! \return Density in [g/cm3]
+  double density(int32_t instance = -1) const;
+
   //! Set the temperature of a cell instance
   //! \param[in] T Temperature in [K]
   //! \param[in] instance Instance index. If -1 is given, the temperature for
@@ -226,6 +238,16 @@ class Cell {
   void set_temperature(
     double T, int32_t instance = -1, bool set_contained = false);
 
+  //! Set the density of a cell instance
+  //! \param[in] density Density [g/cm3]
+  //! \param[in] instance Instance index. If -1 is given, the density
+  //!   for all instances is set.
+  //! \param[in] set_contained If this cell is not filled with a material,
+  //!   collect all contained cells with material fills and set their
+  //!   densities.
+  void set_density(
+    double density, int32_t instance = -1, bool set_contained = false);
+
   int32_t n_instances() const;
 
   //! Set the rotation matrix of a cell instance
@@ -341,6 +363,9 @@ class Cell {
   //! T. The units are sqrt(eV).
   vector<double> sqrtkT_;
 
+  //! \brief Unitless density multiplier(s) within this cell.
+  vector<double> density_mult_;
+
   //! \brief Neighboring cells in the same universe.
   NeighborList neighbors_;