Temperature feedback support in the random ray solver.

docs/source/usersguide/geometry.rst

@@ -248,6 +248,28 @@ The classes :class:`Halfspace`, :class:`Intersection`, :class:`Union`, and
 :class:`Complement` and all instances of :class:`openmc.Region` and can be
 assigned to the :attr:`Cell.region` attribute.
 
+Cells also contain :attr:`Cell.temperature` and :attr:`Cell.density`
+attributes which override the temperature and density of the fill. These can
+be quite useful when temperatures and densities are spatially varying, as the
+alternative would be to add a unique :class:`Material` for each permutation of
+temperature, density, and composition. You can set the temperature (K) and
+density (g/cc) of a cell like so::
+
+  fuel.temperature = 800.0
+  fuel.density = 10.0
+
+The real utility of cell temperatures and densities occurs when a cell is
+replicated across the geometry, such as when a cell is the root geometric element
+in a replicated :ref:`universe<usersguide_universes>` or :ref:`lattice
+<usersguide_lattices>`. In those cases, you can provide a list of temperatures
+and densities to apply a temperature/density field to all of the distributed cells::
+
+  fuel.temperature = [800.0, 900.0, 800.0, 900.0]
+  fuel.density = [10.0, 9.0, 10.0, 9.0]
+
+In this example, the fuel cell is distributed four times in the geometry. Each
+distributed instance then receives its own temperature and density.
+
 .. _usersguide_universes:
 
 ---------

docs/source/usersguide/random_ray.rst

@@ -646,7 +646,9 @@ model to use these multigroup cross sections. An example is given below::
       overwrite_mgxs_library=False,
       mgxs_path="mgxs.h5",
       correction=None,
-      source_energy=None
+      source_energy=None,
+      temperatures=None,
+      temperature_settings=None
   )
 
 The most important parameter to set is the ``method`` parameter, which can be
@@ -733,6 +735,20 @@ distribution for MGXS generation as::
 
   source_energy = openmc.stats.delta_function(2.45e6)
 
+The ``temperatures`` parameter can be provided if temperature-dependent
+multi-group cross sections are desired for multi-physics simulations. An
+individual cross section generation calculation is run for each temperature
+provided, where the materials in the model are set to the temperature. The
+temperature settings used during cross section generation can be specified with the
+``temperature_settings`` parameter. If no ``temperature_settings`` are provided,
+the settings contained in the model will be used. The valid keys and values in the
+``temperature_settings`` dictionary are identical to
+:attr:`openmc.Settings.temperature_settings`; more information can be found in
+:class:`openmc.Settings` . This approach yields isothermal cross section interpolation
+tables, which can be inaccurate for systems with large differences between temperatures
+in each material (often the case in fission reactors). If a more sophisticated
+temperature-dependence is required, we recommend generating cross sections manually.
+
 Ultimately, the methods described above are all just approximations.
 Approximations in the generated MGXS data will fundamentally limit the potential
 accuracy of the random ray solver. However, the methods described above are all

include/openmc/mgxs.h

@@ -113,6 +113,11 @@ class Mgxs {
     const vector<Mgxs*>& micros, const vector<double>& atom_densities,
     int num_group, int num_delay);
 
+  //! \brief Get the number of temperature data points.
+  //!
+  //! @return The number of temperature data points for this MGXS
+  inline int n_temperature_points() { return kTs.size(); }
+
   //! \brief Provides a cross section value given certain parameters
   //!
   //! @param xstype Type of cross section requested, according to the

include/openmc/random_ray/flat_source_domain.h

@@ -100,17 +100,18 @@ class FlatSourceDomain {
   // in model::cells
   vector<int64_t> source_region_offsets_;
 
-  // 2D arrays stored in 1D representing values for all materials x energy
-  // groups
+  // 3D arrays stored in 1D representing values for all materials x temperature
+  // points x energy groups
   int n_materials_;
+  int ntemperature_;
   vector<double> sigma_t_;
   vector<double> nu_sigma_f_;
   vector<double> sigma_f_;
   vector<double> chi_;
   vector<double> kappa_fission_;
 
-  // 3D arrays stored in 1D representing values for all materials x energy
-  // groups x energy groups
+  // 4D arrays stored in 1D representing values for all materials x temperature
+  // points x energy groups x energy groups
   vector<double> sigma_s_;
 
   // The abstract container holding all source region-specific data

include/openmc/random_ray/random_ray.h

@@ -65,6 +65,7 @@ class RandomRay : public Particle {
   vector<double> mesh_fractional_lengths_;
 
   int negroups_;
+  int ntemperature_;
   FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source
                                        // data needed for ray transport
   double distance_travelled_ {0};

include/openmc/random_ray/source_region.h

@@ -146,6 +146,7 @@ class SourceRegionHandle {
 
   // Scalar fields
   int* material_;
+  int* temperature_idx_;
   double* density_mult_;
   int* is_small_;
   int* n_hits_;
@@ -199,6 +200,9 @@ class SourceRegionHandle {
   double& density_mult() { return *density_mult_; }
   const double density_mult() const { return *density_mult_; }
 
+  int& temperature_idx() { return *temperature_idx_; }
+  const int temperature_idx() const { return *temperature_idx_; }
+
   int& is_small() { return *is_small_; }
   const int is_small() const { return *is_small_; }
 
@@ -319,8 +323,9 @@ class SourceRegion {
 
   //---------------------------------------
   // Scalar fields
-
-  int material_ {0};          //!< Index in openmc::model::materials array
+  int material_ {0}; //!< Index in openmc::model::materials array
+  int temperature_idx_ {
+    0}; //!< Index into the MGXS array representing temperature
   double density_mult_ {1.0}; //!< A density multiplier queried from the cell
                               //!< corresponding to the source region.
   OpenMPMutex lock_;
@@ -400,6 +405,9 @@ class SourceRegionContainer {
   int& material(int64_t sr) { return material_[sr]; }
   const int material(int64_t sr) const { return material_[sr]; }
 
+  int& temperature_idx(int64_t sr) { return temperature_idx_[sr]; }
+  const int temperature_idx(int64_t sr) const { return temperature_idx_[sr]; }
+
   double& density_mult(int64_t sr) { return density_mult_[sr]; }
   const double density_mult(int64_t sr) const { return density_mult_[sr]; }
 
@@ -634,6 +642,7 @@ class SourceRegionContainer {
 
   // SoA storage for scalar fields (one item per source region)
   vector<int> material_;
+  vector<int> temperature_idx_;
   vector<double> density_mult_;
   vector<int> is_small_;
   vector<int> n_hits_;

openmc/examples.py

@@ -656,10 +656,19 @@ def slab_mg(num_regions=1, mat_names=None, mgxslib_name='2g.h5') -> openmc.Model
     return model
 
 
-def _generate_c5g7_materials() -> openmc.Materials:
+def _generate_c5g7_materials(second_temp = False) -> openmc.Materials:
     """Generate materials utilizing multi-group cross sections based on the
     the C5G7 Benchmark.
 
+    Parameters
+    ----------
+    second_temp : bool, optional
+        Whether or not the cross sections should contain two temperature datapoints.
+        The first data point is the C5G7 cross sections, which corresponds to a temperature
+        of 294 K. The second data point is the C5G7 cross sections multiplied by 1/2,
+        which corresponds to a temperature of 394 K. This temperature dependence is
+        ficticious; it is used for testing temperature feedback in the random ray solver.
+
     Returns
     -------
     materials : openmc.Materials
@@ -672,9 +681,45 @@ def _generate_c5g7_materials() -> openmc.Materials:
     assembly transport calculations without spatial homogenization"
     """
     # Instantiate the energy group data
+    # MGXS for the UO2 pins.
     group_edges = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6]
     groups = openmc.mgxs.EnergyGroups(group_edges)
 
+    uo2_total = np.array([0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678,
+                          0.5644058])
+    uo2_abs = np.array([8.0248e-03, 3.7174e-03, 2.6769e-02, 9.6236e-02, 3.0020e-02,
+                        1.1126e-01, 2.8278e-01])
+    uo2_scatter_matrix = np.array(
+        [[[0.1275370, 0.0423780, 0.0000094, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
+          [0.0000000, 0.3244560, 0.0016314, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
+          [0.0000000, 0.0000000, 0.4509400, 0.0026792, 0.0000000, 0.0000000, 0.0000000],
+          [0.0000000, 0.0000000, 0.0000000, 0.4525650, 0.0055664, 0.0000000, 0.0000000],
+          [0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000],
+          [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090],
+          [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800]]])
+    uo2_scatter_matrix = np.rollaxis(uo2_scatter_matrix, 0, 3)
+    uo2_fission = np.array([7.21206e-03, 8.19301e-04, 6.45320e-03, 1.85648e-02, 1.78084e-02,
+                            8.30348e-02, 2.16004e-01])
+    uo2_nu_fission = np.array([2.005998e-02, 2.027303e-03, 1.570599e-02, 4.518301e-02,
+                               4.334208e-02, 2.020901e-01, 5.257105e-01])
+    uo2_chi = np.array([5.8791e-01, 4.1176e-01, 3.3906e-04, 1.1761e-07, 0.0000e+00,
+                        0.0000e+00, 0.0000e+00])
+
+    # MGXS for the H2O moderator.
+    h2o_total = np.array([0.15920605, 0.412969593, 0.59030986, 0.58435, 0.718, 1.2544497,
+                          2.650379])
+    h2o_abs = np.array([6.0105e-04, 1.5793e-05, 3.3716e-04, 1.9406e-03, 5.7416e-03,
+                        1.5001e-02, 3.7239e-02])
+    h2o_scatter_matrix = np.array(
+        [[[0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000],
+          [0.0000000, 0.2823340, 0.1299400, 0.0006234, 0.0000480, 0.0000074, 0.0000010],
+          [0.0000000, 0.0000000, 0.3452560, 0.2245700, 0.0169990, 0.0026443, 0.0005034],
+          [0.0000000, 0.0000000, 0.0000000, 0.0910284, 0.4155100, 0.0637320, 0.0121390],
+          [0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290],
+          [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
+          [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]])
+    h2o_scatter_matrix = np.rollaxis(h2o_scatter_matrix, 0, 3)
+
     # Instantiate the 7-group (C5G7) cross section data
     uo2_xsdata = openmc.XSdata('UO2', groups)
     uo2_xsdata.order = 0
@@ -707,29 +752,33 @@ def _generate_c5g7_materials() -> openmc.Materials:
     uo2_xsdata.set_nu_fission(nu_fission)
     uo2_xsdata.set_chi([5.8791e-01, 4.1176e-01, 3.3906e-04, 1.1761e-07, 0.0000e+00,
                         0.0000e+00, 0.0000e+00])
+    uo2_xsdata.set_total(uo2_total, temperature=294.0)
+    uo2_xsdata.set_absorption(uo2_abs, temperature=294.0)
+    uo2_xsdata.set_scatter_matrix(uo2_scatter_matrix, temperature=294.0)
+    uo2_xsdata.set_fission(uo2_fission, temperature=294.0)
+    uo2_xsdata.set_nu_fission(uo2_nu_fission, temperature=294.0)
+    uo2_xsdata.set_chi(uo2_chi, temperature=294.0)
 
     h2o_xsdata = openmc.XSdata('LWTR', groups)
     h2o_xsdata.order = 0
-    h2o_xsdata.set_total([0.15920605, 0.412969593, 0.59030986, 0.58435,
-                          0.718, 1.2544497, 2.650379])
-    h2o_xsdata.set_absorption([6.0105e-04, 1.5793e-05, 3.3716e-04,
-                               1.9406e-03, 5.7416e-03, 1.5001e-02,
-                               3.7239e-02])
-    scatter_matrix = np.array(
-        [[[0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000],
-          [0.0000000, 0.2823340, 0.1299400, 0.0006234,
-              0.0000480, 0.0000074, 0.0000010],
-          [0.0000000, 0.0000000, 0.3452560, 0.2245700,
-              0.0169990, 0.0026443, 0.0005034],
-          [0.0000000, 0.0000000, 0.0000000, 0.0910284,
-              0.4155100, 0.0637320, 0.0121390],
-          [0.0000000, 0.0000000, 0.0000000, 0.0000714,
-              0.1391380, 0.5118200, 0.0612290],
-          [0.0000000, 0.0000000, 0.0000000, 0.0000000,
-              0.0022157, 0.6999130, 0.5373200],
-          [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]])
-    scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
-    h2o_xsdata.set_scatter_matrix(scatter_matrix)
+    h2o_xsdata.set_total(h2o_total, temperature=294.0)
+    h2o_xsdata.set_absorption(h2o_abs, temperature=294.0)
+    h2o_xsdata.set_scatter_matrix(h2o_scatter_matrix, temperature=294.0)
+
+    # Add the second temperature data point if requested.
+    if second_temp:
+        uo2_xsdata.add_temperature(394.0)
+        uo2_xsdata.set_total(0.5 * uo2_total, temperature=394.0)
+        uo2_xsdata.set_absorption(0.5 * uo2_abs, temperature=394.0)
+        uo2_xsdata.set_scatter_matrix(0.5 * uo2_scatter_matrix, temperature=394.0)
+        uo2_xsdata.set_fission(0.5 * uo2_fission, temperature=394.0)
+        uo2_xsdata.set_nu_fission(0.5 * uo2_nu_fission, temperature=394.0)
+        uo2_xsdata.set_chi(uo2_chi, temperature=394.0)
+
+        h2o_xsdata.add_temperature(394.0)
+        h2o_xsdata.set_total(0.5 * h2o_total, temperature=394.0)
+        h2o_xsdata.set_absorption(0.5 * h2o_abs, temperature=394.0)
+        h2o_xsdata.set_scatter_matrix(0.5 * h2o_scatter_matrix, temperature=394.0)
 
     mg_cross_sections = openmc.MGXSLibrary(groups)
     mg_cross_sections.add_xsdatas([uo2_xsdata, h2o_xsdata])
@@ -823,10 +872,19 @@ def _generate_subdivided_pin_cell(uo2, water) -> openmc.Universe:
     return pincell
 
 
-def random_ray_pin_cell() -> openmc.Model:
+def random_ray_pin_cell(second_temp = False) -> openmc.Model:
     """Create a PWR pin cell example using C5G7 cross section data.
     cross section data.
 
+    Parameters
+    ----------
+    second_temp : bool, optional
+        Whether or not the cross sections should contain two temperature datapoints.
+        The first data point is the C5G7 cross sections, which corresponds to a temperature
+        of 294 K. The second data point is the C5G7 cross sections multiplied by 1/2,
+        which corresponds to a temperature of 3934 K. This temperature dependence is
+        ficticious; it is used for testing temperature feedback in the random ray solver.
+
     Returns
     -------
     model : openmc.Model
@@ -837,7 +895,7 @@ def random_ray_pin_cell() -> openmc.Model:
 
     ###########################################################################
     # Create Materials for the problem
-    materials = _generate_c5g7_materials()
+    materials = _generate_c5g7_materials(second_temp)
     uo2 = materials[0]
     water = materials[1]
 
@@ -897,13 +955,22 @@ def random_ray_pin_cell() -> openmc.Model:
     return model
 
 
-def random_ray_lattice() -> openmc.Model:
+def random_ray_lattice(second_temp = False) -> openmc.Model:
     """Create a 2x2 PWR pin cell asymmetrical lattice example.
 
     This model is a 2x2 reflective lattice of fuel pins with one of the lattice
     locations having just moderator instead of a fuel pin. It uses C5G7
     cross section data.
 
+    Parameters
+    ----------
+    second_temp : bool, optional
+        Whether or not the cross sections should contain two temperature datapoints.
+        The first data point is the C5G7 cross sections, which corresponds to a temperature
+        of 294 K. The second data point is the C5G7 cross sections multiplied by 1/2,
+        which corresponds to a temperature of 3934 K. This temperature dependence is
+        ficticious; it is used for testing temperature feedback in the random ray solver.
+
     Returns
     -------
     model : openmc.Model
@@ -914,7 +981,7 @@ def random_ray_lattice() -> openmc.Model:
 
     ###########################################################################
     # Create Materials for the problem
-    materials = _generate_c5g7_materials()
+    materials = _generate_c5g7_materials(second_temp)
     uo2 = materials[0]
     water = materials[1]