Temperature feedback support in the random ray solver. (#3737)

This commit is contained in:
Kevin Sawatzky 2026-02-09 16:54:59 -06:00 committed by GitHub
parent d0346e94ac
commit 4f4930b633
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
29 changed files with 1515 additions and 255 deletions

View file

@ -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:
---------

View file

@ -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

View file

@ -113,6 +113,11 @@ public:
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

View file

@ -100,17 +100,18 @@ public:
// 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

View file

@ -65,6 +65,7 @@ private:
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};

View file

@ -146,6 +146,7 @@ public:
// Scalar fields
int* material_;
int* temperature_idx_;
double* density_mult_;
int* is_small_;
int* n_hits_;
@ -199,6 +200,9 @@ public:
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 @@ public:
//---------------------------------------
// 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 @@ public:
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 @@ private:
// 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_;

View file

@ -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
fictitious; 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
fictitious; 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
fictitious; 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]

View file

@ -13,6 +13,7 @@ import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from ..tallies import ESTIMATOR_TYPES
ROOM_TEMPERATURE_KELVIN = 294.0
class Library:
"""A multi-energy-group and multi-delayed-group cross section library for
@ -954,7 +955,7 @@ class Library:
return pickle.load(f)
def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro',
subdomain=None, apply_domain_chi=False):
subdomain=None, apply_domain_chi=False, temperature=ROOM_TEMPERATURE_KELVIN):
"""Generates an openmc.XSdata object describing a multi-group cross section
dataset for writing to an openmc.MGXSLibrary object.
@ -990,6 +991,9 @@ class Library:
downstream multigroup solvers that precompute a material-specific
chi before the transport solve provides group-wise fluxes. Defaults
to False.
temperature : float, optional
The temperature to set in the XSdata object. Defaults to 294 K
(room temperature).
Returns
-------
@ -1036,6 +1040,7 @@ class Library:
else:
representation = 'isotropic'
xsdata = openmc.XSdata(name, self.energy_groups,
temperatures=[temperature],
representation=representation)
xsdata.num_delayed_groups = self.num_delayed_groups
if self.num_polar > 1 or self.num_azimuthal > 1:
@ -1053,45 +1058,61 @@ class Library:
# Now get xs data itself
if 'nu-transport' in self.mgxs_types and self.correction == 'P0':
mymgxs = self.get_mgxs(domain, 'nu-transport')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
xsdata.set_total_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
elif 'transport' in self.mgxs_types and self.correction == 'P0':
mymgxs = self.get_mgxs(domain, 'transport')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
xsdata.set_total_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
elif 'total' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'total')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
xsdata.set_total_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'absorption' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'absorption')
xsdata.set_absorption_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_absorption_mgxs(mymgxs,
temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'fission')
xsdata.set_fission_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide], subdomain=subdomain)
xsdata.set_fission_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'kappa-fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'kappa-fission')
xsdata.set_kappa_fission_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_kappa_fission_mgxs(mymgxs,
temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'inverse-velocity' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'inverse-velocity')
xsdata.set_inverse_velocity_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_inverse_velocity_mgxs(mymgxs,
temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'nu-fission matrix' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'nu-fission matrix')
xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_nu_fission_mgxs(mymgxs,
temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
@ -1101,7 +1122,9 @@ class Library:
nuc = "sum"
else:
nuc = nuclide
xsdata.set_chi_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuc],
xsdata.set_chi_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuc],
subdomain=subdomain)
if 'chi-prompt' in self.mgxs_types:
@ -1110,8 +1133,10 @@ class Library:
nuc = "sum"
else:
nuc = nuclide
xsdata.set_chi_prompt_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuc], subdomain=subdomain)
xsdata.set_chi_prompt_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuc],
subdomain=subdomain)
if 'chi-delayed' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'chi-delayed')
@ -1119,53 +1144,61 @@ class Library:
nuc = "sum"
else:
nuc = nuclide
xsdata.set_chi_delayed_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuc], subdomain=subdomain)
xsdata.set_chi_delayed_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuc],
subdomain=subdomain)
if 'nu-fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'nu-fission')
xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_nu_fission_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'prompt-nu-fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'prompt-nu-fission')
xsdata.set_prompt_nu_fission_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_prompt_nu_fission_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'prompt-nu-fission matrix' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'prompt-nu-fission matrix')
xsdata.set_prompt_nu_fission_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_prompt_nu_fission_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'delayed-nu-fission' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'delayed-nu-fission')
xsdata.set_delayed_nu_fission_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_delayed_nu_fission_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'delayed-nu-fission matrix' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'delayed-nu-fission matrix')
xsdata.set_delayed_nu_fission_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_delayed_nu_fission_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
if 'beta' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'beta')
xsdata.set_beta_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomain=subdomain)
xsdata.set_beta_mgxs(mymgxs, temperature=temperature, xs_type=xs_type,
nuclide=[nuclide], subdomain=subdomain)
if 'decay-rate' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'decay-rate')
xsdata.set_decay_rate_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomain=subdomain)
xsdata.set_decay_rate_mgxs(mymgxs, temperature=temperature, xs_type=xs_type,
nuclide=[nuclide], subdomain=subdomain)
# If multiplicity matrix is available, prefer that
if 'multiplicity matrix' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'multiplicity matrix')
xsdata.set_multiplicity_matrix_mgxs(mymgxs, xs_type=xs_type,
xsdata.set_multiplicity_matrix_mgxs(mymgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
using_multiplicity = True
@ -1176,6 +1209,7 @@ class Library:
scatt_mgxs = self.get_mgxs(domain, 'scatter matrix')
nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix')
xsdata.set_multiplicity_matrix_mgxs(nuscatt_mgxs, scatt_mgxs,
temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
@ -1188,6 +1222,7 @@ class Library:
nuscatt_mgxs = \
self.get_mgxs(domain, 'consistent nu-scatter matrix')
xsdata.set_multiplicity_matrix_mgxs(nuscatt_mgxs, scatt_mgxs,
temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
@ -1202,7 +1237,8 @@ class Library:
else:
nuscatt_mgxs = \
self.get_mgxs(domain, 'consistent nu-scatter matrix')
xsdata.set_scatter_matrix_mgxs(nuscatt_mgxs, xs_type=xs_type,
xsdata.set_scatter_matrix_mgxs(nuscatt_mgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
else:
@ -1213,7 +1249,8 @@ class Library:
else:
nuscatt_mgxs = \
self.get_mgxs(domain, 'consistent nu-scatter matrix')
xsdata.set_scatter_matrix_mgxs(nuscatt_mgxs, xs_type=xs_type,
xsdata.set_scatter_matrix_mgxs(nuscatt_mgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
@ -1253,14 +1290,15 @@ class Library:
'are ignored since multiplicity or nu-scatter matrices '\
'were not tallied for ' + xsdata_name
warn(msg, RuntimeWarning)
xsdata.set_scatter_matrix_mgxs(scatt_mgxs, xs_type=xs_type,
xsdata.set_scatter_matrix_mgxs(scatt_mgxs, temperature=temperature,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
return xsdata
def create_mg_library(self, xs_type='macro', xsdata_names=None,
apply_domain_chi=False):
apply_domain_chi=False, temperature=ROOM_TEMPERATURE_KELVIN):
"""Creates an openmc.MGXSLibrary object to contain the MGXS data for the
Multi-Group mode of OpenMC.
@ -1286,6 +1324,9 @@ class Library:
downstream multigroup solvers that precompute a material-specific
chi before the transport solve provides group-wise fluxes. Defaults
to False.
temperature : float, optional
The temperature to set in the MGXSLibrary object. Defaults to 294 K
(room temperature).
Returns
-------

View file

@ -515,6 +515,75 @@ class XSdata:
if np.sum(array) > 0:
self._fissionable = True
def add_temperature_data(self, other):
"""This method adds temperature-dependent cross section
values from another XSdata object to this XSdata object.
Note: if a temperature datapoint from 'other' already exists in this
object, it will be overridden.
Parameters
----------
other: openmc.XSdata
The other XSdata object to fetch data from
"""
# Sanity check to make sure they have the same name, energy group structure,
# and delayed group structure
check_value('name', other.name, self.name)
check_value('energy_groups', other.energy_groups, [self.energy_groups])
check_value('delayed_groups', other.num_delayed_groups, [self.num_delayed_groups])
# Add the temperature data.
for temp in other.temperatures:
if temp not in self.temperatures:
self.add_temperature(temp)
if np.all(other.absorption[other._temperature_index(temp)] != None):
self.set_absorption(other.absorption[other._temperature_index(temp)], temp)
if np.all(other.beta[other._temperature_index(temp)] != None):
self.set_beta(other.beta[other._temperature_index(temp)], temp)
if np.all(other.chi[other._temperature_index(temp)] != None):
self.set_chi(other.chi[other._temperature_index(temp)], temp)
if np.all(other.chi_delayed[other._temperature_index(temp)] != None):
self.set_chi_delayed(other.chi_delayed[other._temperature_index(temp)], temp)
if np.all(other.chi_prompt[other._temperature_index(temp)] != None):
self.set_chi_prompt(other.chi_prompt[other._temperature_index(temp)], temp)
if np.all(other.decay_rate[other._temperature_index(temp)] != None):
self.set_decay_rate(other.decay_rate[other._temperature_index(temp)], temp)
if np.all(other.delayed_nu_fission[other._temperature_index(temp)] != None):
self.set_delayed_nu_fission(other.delayed_nu_fission[other._temperature_index(temp)], temp)
if np.all(other.fission[other._temperature_index(temp)] != None):
self.set_fission(other.fission[other._temperature_index(temp)], temp)
if np.all(other.inverse_velocity[other._temperature_index(temp)] != None):
self.set_inverse_velocity(other.inverse_velocity[other._temperature_index(temp)], temp)
if np.all(other.kappa_fission[other._temperature_index(temp)] != None):
self.set_kappa_fission(other.kappa_fission[other._temperature_index(temp)], temp)
if np.all(other.multiplicity_matrix[other._temperature_index(temp)] != None):
self.set_multiplicity_matrix(other.multiplicity_matrix[other._temperature_index(temp)], temp)
if np.all(other.nu_fission[other._temperature_index(temp)] != None):
self.set_nu_fission(other.nu_fission[other._temperature_index(temp)], temp)
if np.all(other.prompt_nu_fission[other._temperature_index(temp)] != None):
self.set_prompt_nu_fission(other.prompt_nu_fission[other._temperature_index(temp)], temp)
if np.all(other.scatter_matrix[other._temperature_index(temp)] != None):
self.set_scatter_matrix(other.scatter_matrix[other._temperature_index(temp)], temp)
if np.all(other.fission[other._temperature_index(temp)] != None):
self.set_total(other.total[other._temperature_index(temp)], temp)
def set_total(self, total, temperature=ROOM_TEMPERATURE_KELVIN):
"""This method sets the cross section for this XSdata object at the
provided temperature.

View file

@ -1695,8 +1695,8 @@ class Model:
self.geometry.get_all_materials().values()
)
@staticmethod
def _auto_generate_mgxs_lib(
self,
model: openmc.model.model,
groups: openmc.mgxs.EnergyGroups,
correction: str | none,
@ -1861,6 +1861,85 @@ class Model:
return sources
@staticmethod
def _isothermal_infinite_media_mgxs(
material: openmc.Material,
groups: openmc.mgxs.EnergyGroups,
nparticles: int,
correction: str | None,
directory: PathLike,
source: openmc.IndependentSource,
temperature_settings: dict,
temperature: float | None = None,
) -> openmc.XSdata:
"""Generate a single MGXS set for one material, where the geometry is an
infinite medium composed of that material at an isothermal temperature value.
Parameters
----------
material : openmc.Material
The material to generate MGXS for
groups : openmc.mgxs.EnergyGroups
Energy group structure for the MGXS.
nparticles : int
Number of particles to simulate per batch when generating MGXS.
correction : str
Transport correction to apply to the MGXS. Options are None and
"P0".
directory : str
Directory to run the simulation in, so as to contain XML files.
source : openmc.IndependentSource
Source to use when generating MGXS.
temperature_settings : dict
A dictionary of temperature settings to use when generating MGXS.
Valid entries for temperature_settings are the same as the valid
entries in openmc.Settings.temperature_settings.
temperature : float, optional
The isothermal temperature value to apply to the material. If not specified,
defaults to the temperature in the material.
Returns
-------
data : openmc.XSdata
The material MGXS for the given temperature isotherm.
"""
model = openmc.Model()
# Set materials on the model
model.materials = [material]
if temperature != None:
model.materials[-1].temperature = temperature
# Settings
model.settings.batches = 100
model.settings.particles = nparticles
model.settings.source = source
model.settings.run_mode = 'fixed source'
model.settings.create_fission_neutrons = False
model.settings.output = {'summary': True, 'tallies': False}
model.settings.temperature = temperature_settings
# Geometry
box = openmc.model.RectangularPrism(
100000.0, 100000.0, boundary_type='reflective')
name = material.name
infinite_cell = openmc.Cell(name=name, fill=model.materials[-1], region=-box)
infinite_universe = openmc.Universe(name=name, cells=[infinite_cell])
model.geometry.root_universe = infinite_universe
# Generate MGXS
mgxs_lib = Model._auto_generate_mgxs_lib(
model, groups, correction, directory)
if temperature != None:
return mgxs_lib.get_xsdata(domain=material, xsdata_name=name,
temperature=temperature)
else:
return mgxs_lib.get_xsdata(domain=material, xsdata_name=name)
def _generate_infinite_medium_mgxs(
self,
groups: openmc.mgxs.EnergyGroups,
@ -1869,13 +1948,17 @@ class Model:
correction: str | None,
directory: PathLike,
source_energy: openmc.stats.Univariate | None = None,
):
temperatures: Sequence[float] | None = None,
temperature_settings: dict | None = None,
) -> None:
"""Generate a MGXS library by running multiple OpenMC simulations, each
representing an infinite medium simulation of a single isolated
material. A discrete source is used to sample particles, with an equal
strength spread across each of the energy groups. This is a highly naive
method that ignores all spatial self shielding effects and all resonance
shielding effects between materials.
shielding effects between materials. If temperature data points are provided,
isothermal cross sections are generated at each temperature point for
each material to build a temperature interpolation table.
Note that in all cases, a discrete source that is uniform over all
energy groups is created (strength = 0.01) to ensure that total cross
@ -1907,50 +1990,78 @@ class Model:
source_energy : openmc.stats.Univariate, optional
Energy distribution to use when generating MGXS data, replacing any
existing sources in the model.
temperatures : Sequence[float], optional
A list of temperatures to generate MGXS at. Each infinite material region
is isothermal at a given temperature data point for cross
section generation.
temperature_settings : dict, optional
A dictionary of temperature settings to use when generating MGXS.
Valid entries for temperature_settings are the same as the valid
entries in openmc.Settings.temperature_settings.
"""
mgxs_sets = []
for material in self.materials:
model = openmc.Model()
# Set materials on the model
model.materials = [material]
src = self._create_mgxs_sources(
groups,
spatial_dist=openmc.stats.Point(),
source_energy=source_energy
)
# Settings
model.settings.batches = 100
model.settings.particles = nparticles
temp_settings = {}
if temperature_settings == None:
temp_settings = self.settings.temperature
else:
temp_settings = temperature_settings
model.settings.source = self._create_mgxs_sources(
groups,
spatial_dist=openmc.stats.Point(),
source_energy=source_energy
)
if temperatures == None:
mgxs_sets = []
for material in self.materials:
xs_data = Model._isothermal_infinite_media_mgxs(
material,
groups,
nparticles,
correction,
directory,
src,
temp_settings
)
mgxs_sets.append(xs_data)
model.settings.run_mode = 'fixed source'
model.settings.create_fission_neutrons = False
# Write the file to disk.
mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
for mgxs_set in mgxs_sets:
mgxs_file.add_xsdata(mgxs_set)
mgxs_file.export_to_hdf5(mgxs_path)
else:
# Build a series of XSData objects, one for each isothermal temperature value.
raw_mgxs_sets = {}
for temperature in temperatures:
raw_mgxs_sets[temperature] = []
for material in self.materials:
xs_data = Model._isothermal_infinite_media_mgxs(
material,
groups,
nparticles,
correction,
directory,
src,
temp_settings,
temperature
)
raw_mgxs_sets[temperature].append(xs_data)
model.settings.output = {'summary': True, 'tallies': False}
# Unpack the isothermal XSData objects and build a single XSData object per material.
mgxs_sets = []
for m in range(len(self.materials)):
mgxs_sets.append(openmc.XSdata(self.materials[m].name, groups))
mgxs_sets[-1].order = 0
for temperature in temperatures:
mgxs_sets[-1].add_temperature_data(raw_mgxs_sets[temperature][m])
# Geometry
box = openmc.model.RectangularPrism(
100000.0, 100000.0, boundary_type='reflective')
name = material.name
infinite_cell = openmc.Cell(name=name, fill=material, region=-box)
infinite_universe = openmc.Universe(name=name, cells=[infinite_cell])
model.geometry.root_universe = infinite_universe
# Add MGXS Tallies
mgxs_lib = self._auto_generate_mgxs_lib(
model, groups, correction, directory)
# Create a MGXS File which can then be written to disk
mgxs_set = mgxs_lib.get_xsdata(domain=material, xsdata_name=name)
mgxs_sets.append(mgxs_set)
# Write the file to disk
mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
for mgxs_set in mgxs_sets:
mgxs_file.add_xsdata(mgxs_set)
mgxs_file.export_to_hdf5(mgxs_path)
# Write the file to disk.
mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
for mgxs_set in mgxs_sets:
mgxs_file.add_xsdata(mgxs_set)
mgxs_file.export_to_hdf5(mgxs_path)
@staticmethod
def _create_stochastic_slab_geometry(
@ -2026,6 +2137,89 @@ class Model:
return geometry, box
@staticmethod
def _isothermal_stochastic_slab_mgxs(
stoch_geom: openmc.Geometry,
groups: openmc.mgxs.EnergyGroups,
nparticles: int,
correction: str | None,
directory: PathLike,
source: openmc.IndependentSource,
temperature_settings: dict,
temperature: float | None = None,
) -> dict[str, openmc.XSdata]:
"""Generate MGXS assuming a stochastic "sandwich" of materials in a layered
slab geometry. If a temperature is specified, all materials in the slab have
their temperatures set to be isothermal at this temperature.
Parameters
----------
stoch_geom : openmc.Geometry
The stochastic slab geometry.
groups : openmc.mgxs.EnergyGroups
Energy group structure for the MGXS.
nparticles : int
Number of particles to simulate per batch when generating MGXS.
correction : str
Transport correction to apply to the MGXS. Options are None and
"P0".
directory : str
Directory to run the simulation in, so as to contain XML files.
source : openmc.IndependentSource
Source to use when generating MGXS.
temperature_settings : dict
A dictionary of temperature settings to use when generating MGXS.
Valid entries for temperature_settings are the same as the valid
entries in openmc.Settings.temperature_settings.
temperature : float, optional
The isothermal temperature value to apply to the materials in the
slab. If not specified, defaults to the temperature in the materials.
Returns
-------
data : dict[str, openmc.XSdata]
A dictionary where the key is the name of the material and the value is the isothermal MGXS.
"""
model = openmc.Model()
model.geometry = stoch_geom
if temperature != None:
for material in model.geometry.get_all_materials().values():
material.temperature = temperature
# Settings
model.settings.batches = 200
model.settings.inactive = 100
model.settings.particles = nparticles
model.settings.output = {'summary': True, 'tallies': False}
model.settings.temperature = temperature_settings
# Define the sources
model.settings.source = source
model.settings.run_mode = 'fixed source'
model.settings.create_fission_neutrons = False
model.settings.output = {'summary': True, 'tallies': False}
# Generate MGXS
mgxs_lib = Model._auto_generate_mgxs_lib(
model, groups, correction, directory)
# Fetch all of the isothermal results.
if temperature != None:
return {
mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name,
temperature=temperature)
for mat in mgxs_lib.domains
}
else:
return {
mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name)
for mat in mgxs_lib.domains
}
def _generate_stochastic_slab_mgxs(
self,
groups: openmc.mgxs.EnergyGroups,
@ -2034,6 +2228,8 @@ class Model:
correction: str | None,
directory: PathLike,
source_energy: openmc.stats.Univariate | None = None,
temperatures: Sequence[float] | None = None,
temperature_settings: dict | None = None,
) -> None:
"""Generate MGXS assuming a stochastic "sandwich" of materials in a layered
slab geometry. While geometry-specific spatial shielding effects are not
@ -2043,7 +2239,9 @@ class Model:
will generate cross sections for all materials in the problem regardless
of type. If this is a fixed source problem, a discrete source is used to
sample particles, with an equal strength spread across each of the
energy groups.
energy groups. If temperature data points are provided,
isothermal cross sections are generated at each temperature point for
the stochastic slab to build a temperature interpolation table.
Parameters
----------
@ -2075,41 +2273,152 @@ class Model:
no sources are defined on the model and the run mode is
'eigenvalue', then a default Watt spectrum source (strength = 0.99)
is added.
temperatures : Sequence[float], optional
A list of temperatures to generate MGXS at. Each infinite material region
is isothermal at a given temperature data point for cross
section generation.
temperature_settings : dict, optional
A dictionary of temperature settings to use when generating MGXS.
Valid entries for temperature_settings are the same as the valid
entries in openmc.Settings.temperature_settings.
"""
model = openmc.Model()
model.materials = self.materials
# Stochastic slab geometry
geo, spatial_distribution = Model._create_stochastic_slab_geometry(
self.materials)
src = self._create_mgxs_sources(
groups,
spatial_dist=spatial_distribution,
source_energy=source_energy
)
temp_settings = {}
if temperature_settings == None:
temp_settings = self.settings.temperature
else:
temp_settings = temperature_settings
if temperatures == None:
mgxs_sets = Model._isothermal_stochastic_slab_mgxs(
geo,
groups,
nparticles,
correction,
directory,
src,
temp_settings
).values()
# Write the file to disk.
mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
for mgxs_set in mgxs_sets:
mgxs_file.add_xsdata(mgxs_set)
mgxs_file.export_to_hdf5(mgxs_path)
else:
# Build a series of XSData objects, one for each isothermal temperature value.
raw_mgxs_sets = {}
for temperature in temperatures:
raw_mgxs_sets[temperature] = Model._isothermal_stochastic_slab_mgxs(
geo,
groups,
nparticles,
correction,
directory,
src,
temp_settings,
temperature
)
# Unpack the isothermal XSData objects and build a single XSData object per material.
mgxs_sets = []
for mat in self.materials:
mgxs_sets.append(openmc.XSdata(mat.name, groups))
mgxs_sets[-1].order = 0
for temperature in temperatures:
mgxs_sets[-1].add_temperature_data(raw_mgxs_sets[temperature][mat.name])
# Write the file to disk.
mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
for mgxs_set in mgxs_sets:
mgxs_file.add_xsdata(mgxs_set)
mgxs_file.export_to_hdf5(mgxs_path)
@staticmethod
def _isothermal_materialwise_mgxs(
input_model: openmc.Model,
groups: openmc.mgxs.EnergyGroups,
nparticles: int,
correction: str | None,
directory: PathLike,
temperature_settings: dict,
temperature: float | None = None,
) -> dict[str, openmc.XSdata]:
"""Generate a material-wise MGXS library for the model by running the
original continuous energy OpenMC simulation. If a temperature is
specified, each material in the input model is set to that temperature.
Otherwise, the original material temperatures are used. If temperature
data points are provided, isothermal cross sections are generated at
each temperature point for the whole model to build a temperature
interpolation table.
Parameters
----------
input_model : openmc.Model
The model to use when computing material-wise MGXS.
groups : openmc.mgxs.EnergyGroups
Energy group structure for the MGXS.
nparticles : int
Number of particles to simulate per batch when generating MGXS.
correction : str
Transport correction to apply to the MGXS. Options are None and
"P0".
directory : str
Directory to run the simulation in, so as to contain XML files.
temperature_settings : dict
A dictionary of temperature settings to use when generating MGXS.
Valid entries for temperature_settings are the same as the valid
entries in openmc.Settings.temperature_settings.
temperature : float, optional
The isothermal temperature value to apply to the materials in the
input model. If not specified, defaults to the temperatures in the
materials.
Returns
-------
data : dict[str, openmc.XSdata]
A dictionary where the key is the name of the material and the value is the isothermal MGXS.
"""
model = copy.deepcopy(input_model)
model.tallies = openmc.Tallies()
if temperature != None:
for material in model.geometry.get_all_materials().values():
material.temperature = temperature
# Settings
model.settings.batches = 200
model.settings.inactive = 100
model.settings.particles = nparticles
model.settings.output = {'summary': True, 'tallies': False}
model.settings.temperature = temperature_settings
# Stochastic slab geometry
model.geometry, spatial_distribution = Model._create_stochastic_slab_geometry(
model.materials)
# Define the sources
model.settings.source = self._create_mgxs_sources(
groups,
spatial_dist=spatial_distribution,
source_energy=source_energy
)
model.settings.run_mode = 'fixed source'
model.settings.create_fission_neutrons = False
model.settings.output = {'summary': True, 'tallies': False}
# Add MGXS Tallies
mgxs_lib = self._auto_generate_mgxs_lib(
# Generate MGXS
mgxs_lib = Model._auto_generate_mgxs_lib(
model, groups, correction, directory)
names = [mat.name for mat in mgxs_lib.domains]
# Create a MGXS File which can then be written to disk
mgxs_file = mgxs_lib.create_mg_library(xs_type='macro', xsdata_names=names)
mgxs_file.export_to_hdf5(mgxs_path)
# Fetch all of the isothermal results.
if temperature != None:
return {
mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name,
temperature=temperature)
for mat in mgxs_lib.domains
}
else:
return {
mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name)
for mat in mgxs_lib.domains
}
def _generate_material_wise_mgxs(
self,
@ -2118,6 +2427,8 @@ class Model:
mgxs_path: PathLike,
correction: str | None,
directory: PathLike,
temperatures: Sequence[float] | None = None,
temperature_settings: dict | None = None,
) -> None:
"""Generate a material-wise MGXS library for the model by running the
original continuous energy OpenMC simulation of the full material
@ -2142,26 +2453,63 @@ class Model:
"P0".
directory : PathLike
Directory to run the simulation in, so as to contain XML files.
temperatures : Sequence[float], optional
A list of temperatures to generate MGXS at. Each infinite material region
is isothermal at a given temperature data point for cross
section generation.
temperature_settings : dict, optional
A dictionary of temperature settings to use when generating MGXS.
Valid entries for temperature_settings are the same as the valid
entries in openmc.Settings.temperature_settings.
"""
model = copy.deepcopy(self)
model.tallies = openmc.Tallies()
temp_settings = {}
if temperature_settings == None:
temp_settings = self.settings.temperature
else:
temp_settings = temperature_settings
# Settings
model.settings.batches = 200
model.settings.inactive = 100
model.settings.particles = nparticles
model.settings.output = {'summary': True, 'tallies': False}
if temperatures == None:
mgxs_sets = Model._isothermal_materialwise_mgxs(
self,
groups,
nparticles,
correction,
directory,
temp_settings
).values()
# Add MGXS Tallies
mgxs_lib = self._auto_generate_mgxs_lib(
model, groups, correction, directory)
# Write the file to disk.
mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
for mgxs_set in mgxs_sets:
mgxs_file.add_xsdata(mgxs_set)
mgxs_file.export_to_hdf5(mgxs_path)
else:
# Build a series of XSData objects, one for each isothermal temperature value.
raw_mgxs_sets = {}
for temperature in temperatures:
raw_mgxs_sets[temperature] = Model._isothermal_materialwise_mgxs(
self,
groups,
nparticles,
correction,
directory,
temp_settings,
temperature
)
names = [mat.name for mat in mgxs_lib.domains]
# Unpack the isothermal XSData objects and build a single XSData object per material.
mgxs_sets = []
for mat in self.materials:
mgxs_sets.append(openmc.XSdata(mat.name, groups))
mgxs_sets[-1].order = 0
for temperature in temperatures:
mgxs_sets[-1].add_temperature_data(raw_mgxs_sets[temperature][mat.name])
# Create a MGXS File which can then be written to disk
mgxs_file = mgxs_lib.create_mg_library(
xs_type='macro', xsdata_names=names)
mgxs_file.export_to_hdf5(mgxs_path)
# Write the file to disk.
mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
for mgxs_set in mgxs_sets:
mgxs_file.add_xsdata(mgxs_set)
mgxs_file.export_to_hdf5(mgxs_path)
def convert_to_multigroup(
self,
@ -2172,6 +2520,8 @@ class Model:
mgxs_path: PathLike = "mgxs.h5",
correction: str | None = None,
source_energy: openmc.stats.Univariate | None = None,
temperatures: Sequence[float] | None = None,
temperature_settings: dict | None = None,
):
"""Convert all materials from continuous energy to multigroup.
@ -2212,6 +2562,14 @@ class Model:
'eigenvalue', then a default Watt spectrum source (strength = 0.99)
is added. Note that this argument is only used when using the
"stochastic_slab" or "infinite_medium" MGXS generation methods.
temperatures : Sequence[float], optional
A list of temperatures to generate MGXS at. Each infinite material region
is isothermal at a given temperature data point for cross
section generation.
temperature_settings : dict, optional
A dictionary of temperature settings to use when generating MGXS.
Valid entries for temperature_settings are the same as the valid
entries in openmc.Settings.temperature_settings.
"""
if isinstance(groups, str):
groups = openmc.mgxs.EnergyGroups(groups)
@ -2241,13 +2599,16 @@ class Model:
if not Path(mgxs_path).is_file() or overwrite_mgxs_library:
if method == "infinite_medium":
self._generate_infinite_medium_mgxs(
groups, nparticles, mgxs_path, correction, tmpdir, source_energy)
groups, nparticles, mgxs_path, correction, tmpdir, source_energy,
temperatures, temperature_settings)
elif method == "material_wise":
self._generate_material_wise_mgxs(
groups, nparticles, mgxs_path, correction, tmpdir)
groups, nparticles, mgxs_path, correction, tmpdir,
temperatures, temperature_settings)
elif method == "stochastic_slab":
self._generate_stochastic_slab_mgxs(
groups, nparticles, mgxs_path, correction, tmpdir, source_energy)
groups, nparticles, mgxs_path, correction, tmpdir, source_energy,
temperatures, temperature_settings)
else:
raise ValueError(
f'MGXS generation method "{method}" not recognized')

View file

@ -169,13 +169,13 @@ vector<vector<double>> MgxsInterface::get_mat_kTs()
continue;
// Get temperature of cell (rounding to nearest integer)
double sqrtkT =
cell->sqrtkT_.size() == 1 ? cell->sqrtkT_[j] : cell->sqrtkT_[0];
double kT = sqrtkT * sqrtkT;
for (int k = 0; k < cell->sqrtkT_.size(); ++k) {
double kT = cell->sqrtkT_[k] * cell->sqrtkT_[k];
// Add temperature if it hasn't already been added
if (!contains(kTs[i_material], kT)) {
kTs[i_material].push_back(kT);
// Add temperature if it hasn't already been added
if (!contains(kTs[i_material], kT)) {
kTs[i_material].push_back(kT);
}
}
}
}

View file

@ -109,11 +109,13 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
// Add scattering + fission source
int material = srh.material();
int temp = srh.temperature_idx();
double density_mult = srh.density_mult();
if (material != MATERIAL_VOID) {
double inverse_k_eff = 1.0 / k_eff_;
int material_offset = material * negroups_;
int scatter_offset = material * negroups_ * negroups_;
const int material_offset = (material * ntemperature_ + temp) * negroups_;
const int scatter_offset =
(material * ntemperature_ + temp) * negroups_ * negroups_;
for (int g_out = 0; g_out < negroups_; g_out++) {
double sigma_t = sigma_t_[material_offset + g_out] * density_mult;
double scatter_source = 0.0;
@ -193,6 +195,7 @@ void FlatSourceDomain::set_flux_to_flux_plus_source(
int64_t sr, double volume, int g)
{
int material = source_regions_.material(sr);
int temp = source_regions_.temperature_idx(sr);
if (material == MATERIAL_VOID) {
source_regions_.scalar_flux_new(sr, g) /= volume;
if (settings::run_mode == RunMode::FIXED_SOURCE) {
@ -201,8 +204,9 @@ void FlatSourceDomain::set_flux_to_flux_plus_source(
source_regions_.volume_sq(sr);
}
} else {
double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g] *
source_regions_.density_mult(sr);
double sigma_t =
sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] *
source_regions_.density_mult(sr);
source_regions_.scalar_flux_new(sr, g) /= (sigma_t * volume);
source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);
}
@ -328,6 +332,7 @@ void FlatSourceDomain::compute_k_eff()
}
int material = source_regions_.material(sr);
int temp = source_regions_.temperature_idx(sr);
if (material == MATERIAL_VOID) {
continue;
}
@ -336,8 +341,9 @@ void FlatSourceDomain::compute_k_eff()
double sr_fission_source_new = 0;
for (int g = 0; g < negroups_; g++) {
double nu_sigma_f = nu_sigma_f_[material * negroups_ + g] *
source_regions_.density_mult(sr);
double nu_sigma_f =
nu_sigma_f_[(material * ntemperature_ + temp) * negroups_ + g] *
source_regions_.density_mult(sr);
sr_fission_source_old +=
nu_sigma_f * source_regions_.scalar_flux_old(sr, g);
sr_fission_source_new +=
@ -560,6 +566,7 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const
#pragma omp parallel for reduction(+ : simulation_external_source_strength)
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
int temp = source_regions_.temperature_idx(sr);
double volume = source_regions_.volume(sr) * simulation_volume_;
for (int g = 0; g < negroups_; g++) {
// For non-void regions, we store the external source pre-divided by
@ -567,8 +574,8 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const
// to get the total source strength in the expected units.
double sigma_t = 1.0;
if (material != MATERIAL_VOID) {
sigma_t =
sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr);
sigma_t = sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] *
source_regions_.density_mult(sr);
}
simulation_external_source_strength +=
source_regions_.external_source(sr, g) * sigma_t * volume;
@ -624,9 +631,9 @@ void FlatSourceDomain::random_ray_tally()
// source strength.
double volume = source_regions_.volume(sr) * simulation_volume_;
int material = source_regions_.material(sr);
double material = source_regions_.material(sr);
int temp = source_regions_.temperature_idx(sr);
double density_mult = source_regions_.density_mult(sr);
for (int g = 0; g < negroups_; g++) {
double flux =
source_regions_.scalar_flux_new(sr, g) * source_normalization_factor;
@ -643,21 +650,27 @@ void FlatSourceDomain::random_ray_tally()
case SCORE_TOTAL:
if (material != MATERIAL_VOID) {
score =
flux * volume * sigma_t_[material * negroups_ + g] * density_mult;
flux * volume *
sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] *
density_mult;
}
break;
case SCORE_FISSION:
if (material != MATERIAL_VOID) {
score =
flux * volume * sigma_f_[material * negroups_ + g] * density_mult;
flux * volume *
sigma_f_[(material * ntemperature_ + temp) * negroups_ + g] *
density_mult;
}
break;
case SCORE_NU_FISSION:
if (material != MATERIAL_VOID) {
score = flux * volume * nu_sigma_f_[material * negroups_ + g] *
density_mult;
score =
flux * volume *
nu_sigma_f_[(material * ntemperature_ + temp) * negroups_ + g] *
density_mult;
}
break;
@ -666,8 +679,10 @@ void FlatSourceDomain::random_ray_tally()
break;
case SCORE_KAPPA_FISSION:
score = flux * volume * kappa_fission_[material * negroups_ + g] *
density_mult;
score =
flux * volume *
kappa_fission_[(material * ntemperature_ + temp) * negroups_ + g] *
density_mult;
break;
default:
@ -925,12 +940,14 @@ void FlatSourceDomain::output_to_vtk() const
float total_fission = 0.0;
if (fsr >= 0) {
int mat = source_regions_.material(fsr);
int temp = source_regions_.temperature_idx(fsr);
if (mat != MATERIAL_VOID) {
for (int g = 0; g < negroups_; g++) {
int64_t source_element = fsr * negroups_ + g;
float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g);
double sigma_f = sigma_f_[mat * negroups_ + g] *
source_regions_.density_mult(fsr);
double sigma_f =
sigma_f_[(mat * ntemperature_ + temp) * negroups_ + g] *
source_regions_.density_mult(fsr);
total_fission += sigma_f * flux;
}
}
@ -944,6 +961,7 @@ void FlatSourceDomain::output_to_vtk() const
for (int i = 0; i < Nx * Ny * Nz; i++) {
int64_t fsr = voxel_indices[i];
int mat = source_regions_.material(fsr);
int temp = source_regions_.temperature_idx(fsr);
float total_external = 0.0f;
if (fsr >= 0) {
for (int g = 0; g < negroups_; g++) {
@ -951,7 +969,7 @@ void FlatSourceDomain::output_to_vtk() const
// multiply it back to get the true external source.
double sigma_t = 1.0;
if (mat != MATERIAL_VOID) {
sigma_t = sigma_t_[mat * negroups_ + g] *
sigma_t = sigma_t_[(mat * ntemperature_ + temp) * negroups_ + g] *
source_regions_.density_mult(fsr);
}
total_external += source_regions_.external_source(fsr, g) * sigma_t;
@ -1131,67 +1149,74 @@ void FlatSourceDomain::flatten_xs()
{
// Temperature and angle indices, if using multiple temperature
// data sets and/or anisotropic data sets.
// TODO: Currently assumes we are only using single temp/single angle data.
const int t = 0;
// TODO: Currently assumes we are only using single angle data.
const int a = 0;
n_materials_ = data::mg.macro_xs_.size();
ntemperature_ = 1;
for (int i = 0; i < n_materials_; i++) {
ntemperature_ =
std::max(ntemperature_, data::mg.macro_xs_[i].n_temperature_points());
}
for (int i = 0; i < n_materials_; i++) {
auto& m = data::mg.macro_xs_[i];
for (int g_out = 0; g_out < negroups_; g_out++) {
if (m.exists_in_model) {
double sigma_t =
m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a);
sigma_t_.push_back(sigma_t);
for (int t = 0; t < ntemperature_; t++) {
for (int g_out = 0; g_out < negroups_; g_out++) {
if (m.exists_in_model && t < m.n_temperature_points()) {
double sigma_t =
m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a);
sigma_t_.push_back(sigma_t);
if (sigma_t < MINIMUM_MACRO_XS) {
Material* mat = model::materials[i].get();
warning(fmt::format(
"Material \"{}\" (id: {}) has a group {} total cross section "
"({:.3e}) below the minimum threshold "
"({:.3e}). Material will be treated as pure void.",
mat->name(), mat->id(), g_out, sigma_t, MINIMUM_MACRO_XS));
}
if (sigma_t < MINIMUM_MACRO_XS) {
Material* mat = model::materials[i].get();
warning(fmt::format(
"Material \"{}\" (id: {}) has a group {} total cross section "
"({:.3e}) below the minimum threshold "
"({:.3e}). Material will be treated as pure void.",
mat->name(), mat->id(), g_out, sigma_t, MINIMUM_MACRO_XS));
}
double nu_sigma_f =
m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a);
nu_sigma_f_.push_back(nu_sigma_f);
double nu_sigma_f =
m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a);
nu_sigma_f_.push_back(nu_sigma_f);
double sigma_f =
m.get_xs(MgxsType::FISSION, g_out, NULL, NULL, NULL, t, a);
sigma_f_.push_back(sigma_f);
double sigma_f =
m.get_xs(MgxsType::FISSION, g_out, NULL, NULL, NULL, t, a);
sigma_f_.push_back(sigma_f);
double chi =
m.get_xs(MgxsType::CHI_PROMPT, g_out, &g_out, NULL, NULL, t, a);
if (!std::isfinite(chi)) {
// MGXS interface may return NaN in some cases, such as when material
// is fissionable but has very small sigma_f.
chi = 0.0;
}
chi_.push_back(chi);
double chi =
m.get_xs(MgxsType::CHI_PROMPT, g_out, &g_out, NULL, NULL, t, a);
if (!std::isfinite(chi)) {
// MGXS interface may return NaN in some cases, such as when
// material is fissionable but has very small sigma_f.
chi = 0.0;
}
chi_.push_back(chi);
double kappa_fission =
m.get_xs(MgxsType::KAPPA_FISSION, g_out, NULL, NULL, NULL, t, a);
kappa_fission_.push_back(kappa_fission);
double kappa_fission =
m.get_xs(MgxsType::KAPPA_FISSION, g_out, NULL, NULL, NULL, t, a);
kappa_fission_.push_back(kappa_fission);
for (int g_in = 0; g_in < negroups_; g_in++) {
double sigma_s =
m.get_xs(MgxsType::NU_SCATTER, g_in, &g_out, NULL, NULL, t, a);
sigma_s_.push_back(sigma_s);
// For transport corrected XS data, diagonal elements may be negative.
// In this case, set a flag to enable transport stabilization for the
// simulation.
if (g_out == g_in && sigma_s < 0.0)
is_transport_stabilization_needed_ = true;
}
} else {
sigma_t_.push_back(0);
nu_sigma_f_.push_back(0);
sigma_f_.push_back(0);
chi_.push_back(0);
kappa_fission_.push_back(0);
for (int g_in = 0; g_in < negroups_; g_in++) {
sigma_s_.push_back(0);
for (int g_in = 0; g_in < negroups_; g_in++) {
double sigma_s =
m.get_xs(MgxsType::NU_SCATTER, g_in, &g_out, NULL, NULL, t, a);
sigma_s_.push_back(sigma_s);
// For transport corrected XS data, diagonal elements may be
// negative. In this case, set a flag to enable transport
// stabilization for the simulation.
if (g_out == g_in && sigma_s < 0.0)
is_transport_stabilization_needed_ = true;
}
} else {
sigma_t_.push_back(0);
nu_sigma_f_.push_back(0);
sigma_f_.push_back(0);
chi_.push_back(0);
kappa_fission_.push_back(0);
for (int g_in = 0; g_in < negroups_; g_in++) {
sigma_s_.push_back(0);
}
}
}
}
@ -1263,12 +1288,14 @@ void FlatSourceDomain::set_adjoint_sources()
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
int temp = source_regions_.temperature_idx(sr);
if (material == MATERIAL_VOID) {
continue;
}
for (int g = 0; g < negroups_; g++) {
double sigma_t =
sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr);
sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] *
source_regions_.density_mult(sr);
source_regions_.external_source(sr, g) /= sigma_t;
}
}
@ -1279,15 +1306,17 @@ void FlatSourceDomain::transpose_scattering_matrix()
// Transpose the inner two dimensions for each material
#pragma omp parallel for
for (int m = 0; m < n_materials_; ++m) {
int material_offset = m * negroups_ * negroups_;
for (int i = 0; i < negroups_; ++i) {
for (int j = i + 1; j < negroups_; ++j) {
// Calculate indices of the elements to swap
int idx1 = material_offset + i * negroups_ + j;
int idx2 = material_offset + j * negroups_ + i;
for (int t = 0; t < ntemperature_; t++) {
int material_offset = (m * ntemperature_ + t) * negroups_ * negroups_;
for (int i = 0; i < negroups_; ++i) {
for (int j = i + 1; j < negroups_; ++j) {
// Calculate indices of the elements to swap
int idx1 = material_offset + i * negroups_ + j;
int idx2 = material_offset + j * negroups_ + i;
// Swap the elements to transpose the matrix
std::swap(sigma_s_[idx1], sigma_s_[idx2]);
// Swap the elements to transpose the matrix
std::swap(sigma_s_[idx1], sigma_s_[idx2]);
}
}
}
}
@ -1507,18 +1536,26 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
int gs_i_cell = gs.lowest_coord().cell();
Cell& cell = *model::cells[gs_i_cell];
int material = cell.material(gs.cell_instance());
int temp = 0;
// If material total XS is extremely low, just set it to void to avoid
// problems with 1/Sigma_t
for (int g = 0; g < negroups_; g++) {
double sigma_t = sigma_t_[material * negroups_ + g];
if (sigma_t < MINIMUM_MACRO_XS) {
material = MATERIAL_VOID;
break;
if (material != MATERIAL_VOID) {
temp = data::mg.macro_xs_[material].get_temperature_index(
cell.sqrtkT(gs.cell_instance()));
for (int g = 0; g < negroups_; g++) {
double sigma_t =
sigma_t_[(material * ntemperature_ + temp) * negroups_ + g];
if (sigma_t < MINIMUM_MACRO_XS) {
material = MATERIAL_VOID;
temp = 0;
break;
}
}
}
handle.material() = material;
handle.temperature_idx() = temp;
handle.density_mult() = cell.density_mult(gs.cell_instance());
@ -1551,7 +1588,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
if (material != C_NONE) {
for (int g = 0; g < negroups_; g++) {
double sigma_t =
sigma_t_[material * negroups_ + g] * handle.density_mult();
sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] *
handle.density_mult();
handle.external_source(g) /= sigma_t;
}
}
@ -1626,6 +1664,7 @@ void FlatSourceDomain::apply_transport_stabilization()
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
int temp = source_regions_.temperature_idx(sr);
double density_mult = source_regions_.density_mult(sr);
if (material == MATERIAL_VOID) {
continue;
@ -1634,10 +1673,14 @@ void FlatSourceDomain::apply_transport_stabilization()
// Only apply stabilization if the diagonal (in-group) scattering XS is
// negative
double sigma_s =
sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g] *
sigma_s_[((material * ntemperature_ + temp) * negroups_ + g) *
negroups_ +
g] *
density_mult;
if (sigma_s < 0.0) {
double sigma_t = sigma_t_[material * negroups_ + g] * density_mult;
double sigma_t =
sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] *
density_mult;
double phi_new = source_regions_.scalar_flux_new(sr, g);
double phi_old = source_regions_.scalar_flux_old(sr, g);

View file

@ -43,13 +43,16 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
// Add scattering + fission source
int material = srh.material();
int temp = srh.temperature_idx();
double density_mult = srh.density_mult();
if (material != MATERIAL_VOID) {
double inverse_k_eff = 1.0 / k_eff_;
MomentMatrix invM = srh.mom_matrix().inverse();
for (int g_out = 0; g_out < negroups_; g_out++) {
double sigma_t = sigma_t_[material * negroups_ + g_out] * density_mult;
double sigma_t =
sigma_t_[(material * ntemperature_ + temp) * negroups_ + g_out] *
density_mult;
double scatter_flat = 0.0f;
double fission_flat = 0.0f;
@ -62,12 +65,16 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
MomentArray flux_linear = srh.flux_moments_old(g_in);
// Handles for cross sections
double sigma_s = sigma_s_[material * negroups_ * negroups_ +
g_out * negroups_ + g_in] *
density_mult;
double sigma_s =
sigma_s_[((material * ntemperature_ + temp) * negroups_ + g_out) *
negroups_ +
g_in] *
density_mult;
double nu_sigma_f =
nu_sigma_f_[material * negroups_ + g_in] * density_mult;
double chi = chi_[material * negroups_ + g_out];
nu_sigma_f_[(material * ntemperature_ + temp) * negroups_ + g_in] *
density_mult;
double chi =
chi_[(material * ntemperature_ + temp) * negroups_ + g_out];
// Compute source terms for flat and linear components of the flux
scatter_flat += sigma_s * flux_flat;

View file

@ -432,11 +432,13 @@ void RandomRay::attenuate_flux_flat_source(
// Get material
int material = srh.material();
int temp = srh.temperature_idx();
// MOC incoming flux attenuation + source contribution/attenuation equation
for (int g = 0; g < negroups_; g++) {
float sigma_t =
domain_->sigma_t_[material * negroups_ + g] * srh.density_mult();
domain_->sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] *
srh.density_mult();
float tau = sigma_t * distance;
float exponential = cjosey_exponential(tau); // exponential = 1 - exp(-tau)
float new_delta_psi = (angular_flux_[g] - srh.source(g)) * exponential;
@ -531,6 +533,7 @@ void RandomRay::attenuate_flux_linear_source(
n_event()++;
int material = srh.material();
int temp = srh.temperature_idx();
Position& centroid = srh.centroid();
Position midpoint = r + u() * (distance / 2.0);
@ -560,7 +563,8 @@ void RandomRay::attenuate_flux_linear_source(
// Compute tau, the optical thickness of the ray segment
float sigma_t =
domain_->sigma_t_[material * negroups_ + g] * srh.density_mult();
domain_->sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] *
srh.density_mult();
float tau = sigma_t * distance;
// If tau is very small, set it to zero to avoid numerical issues.
@ -765,6 +769,7 @@ void RandomRay::attenuate_flux_linear_source_void(
void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
{
domain_ = domain;
ntemperature_ = domain->ntemperature_;
// Reset particle event counter
n_event() = 0;

View file

@ -79,10 +79,6 @@ void validate_random_ray_inputs()
fatal_error("Anisotropic MGXS detected. Only isotropic XS data sets "
"supported in random ray mode.");
}
if (material.get_xsdata().size() > 1) {
warning("Non-isothermal MGXS detected. Only isothermal XS data sets "
"supported in random ray mode. Using lowest temperature.");
}
for (int g = 0; g < data::mg.num_energy_groups_; g++) {
if (material.exists_in_model) {
// Temperature and angle indices, if using multiple temperature

View file

@ -11,11 +11,11 @@ namespace openmc {
//==============================================================================
SourceRegionHandle::SourceRegionHandle(SourceRegion& sr)
: negroups_(sr.scalar_flux_old_.size()), material_(&sr.material_),
density_mult_(&sr.density_mult_), is_small_(&sr.is_small_),
n_hits_(&sr.n_hits_), is_linear_(sr.source_gradients_.size() > 0),
lock_(&sr.lock_), volume_(&sr.volume_), volume_t_(&sr.volume_t_),
volume_sq_(&sr.volume_sq_), volume_sq_t_(&sr.volume_sq_t_),
volume_naive_(&sr.volume_naive_),
temperature_idx_(&sr.temperature_idx_), density_mult_(&sr.density_mult_),
is_small_(&sr.is_small_), n_hits_(&sr.n_hits_),
is_linear_(sr.source_gradients_.size() > 0), lock_(&sr.lock_),
volume_(&sr.volume_), volume_t_(&sr.volume_t_), volume_sq_(&sr.volume_sq_),
volume_sq_t_(&sr.volume_sq_t_), volume_naive_(&sr.volume_naive_),
position_recorded_(&sr.position_recorded_),
external_source_present_(&sr.external_source_present_),
position_(&sr.position_), centroid_(&sr.centroid_),
@ -71,6 +71,7 @@ void SourceRegionContainer::push_back(const SourceRegion& sr)
// Scalar fields
material_.push_back(sr.material_);
temperature_idx_.push_back(sr.temperature_idx_);
density_mult_.push_back(sr.density_mult_);
is_small_.push_back(sr.is_small_);
n_hits_.push_back(sr.n_hits_);
@ -125,6 +126,7 @@ void SourceRegionContainer::assign(
// Clear existing data
n_source_regions_ = 0;
material_.clear();
temperature_idx_.clear();
density_mult_.clear();
is_small_.clear();
n_hits_.clear();
@ -183,6 +185,7 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr)
SourceRegionHandle handle;
handle.negroups_ = negroups();
handle.material_ = &material(sr);
handle.temperature_idx_ = &temperature_idx(sr);
handle.density_mult_ = &density_mult(sr);
handle.is_small_ = &is_small(sr);
handle.n_hits_ = &n_hits(sr);

View file

@ -0,0 +1,70 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="UO2__2_4__" depletable="true" temperature="294.0">
<density value="1.0" units="macro"/>
<macroscopic name="UO2__2_4__"/>
</material>
<material id="2" name="Zircaloy" temperature="294.0">
<density value="1.0" units="macro"/>
<macroscopic name="Zircaloy"/>
</material>
<material id="3" name="Hot_borated_water" temperature="294.0">
<density value="1.0" units="macro"/>
<macroscopic name="Hot_borated_water"/>
</material>
</materials>
<geometry>
<cell id="1" name="Fuel" material="1" region="-1" universe="0">
<temperature>395.0</temperature>
</cell>
<cell id="2" name="Cladding" material="2" region="1 -2" universe="0"/>
<cell id="3" name="Water" material="3" region="2 3 -4 5 -6" universe="0"/>
<surface id="1" name="Fuel OR" type="z-cylinder" coeffs="0 0 0.39218"/>
<surface id="2" name="Clad OR" type="z-cylinder" coeffs="0 0 0.4572"/>
<surface id="3" name="left" type="x-plane" boundary="reflective" coeffs="-0.63"/>
<surface id="4" name="right" type="x-plane" boundary="reflective" coeffs="0.63"/>
<surface id="5" name="bottom" type="y-plane" boundary="reflective" coeffs="-0.63"/>
<surface id="6" name="top" type="y-plane" boundary="reflective" coeffs="0.63"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
</space>
<constraints>
<fissionable>true</fissionable>
</constraints>
</source>
<energy_mode>multi-group</energy_mode>
<temperature_method>nearest</temperature_method>
<temperature_multipole>true</temperature_multipole>
<temperature_range>200.0 400.0</temperature_range>
<temperature_tolerance>200.0</temperature_tolerance>
<random_ray>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
</space>
</source>
<distance_inactive>30.0</distance_inactive>
<distance_active>150.0</distance_active>
<source_region_meshes>
<mesh id="1">
<domain id="0" type="universe"/>
</mesh>
</source_region_meshes>
<source_shape>linear</source_shape>
</random_ray>
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-0.63 -0.63</lower_left>
<upper_right>0.63 0.63</upper_right>
</mesh>
</settings>
</model>

View file

@ -0,0 +1,2 @@
k-combined:
7.499800E-01 1.615317E-02

View file

@ -0,0 +1,70 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="UO2__2_4__" depletable="true" temperature="294.0">
<density value="1.0" units="macro"/>
<macroscopic name="UO2__2_4__"/>
</material>
<material id="2" name="Zircaloy" temperature="294.0">
<density value="1.0" units="macro"/>
<macroscopic name="Zircaloy"/>
</material>
<material id="3" name="Hot_borated_water" temperature="294.0">
<density value="1.0" units="macro"/>
<macroscopic name="Hot_borated_water"/>
</material>
</materials>
<geometry>
<cell id="1" name="Fuel" material="1" region="-1" universe="0">
<temperature>395.0</temperature>
</cell>
<cell id="2" name="Cladding" material="2" region="1 -2" universe="0"/>
<cell id="3" name="Water" material="3" region="2 3 -4 5 -6" universe="0"/>
<surface id="1" name="Fuel OR" type="z-cylinder" coeffs="0 0 0.39218"/>
<surface id="2" name="Clad OR" type="z-cylinder" coeffs="0 0 0.4572"/>
<surface id="3" name="left" type="x-plane" boundary="reflective" coeffs="-0.63"/>
<surface id="4" name="right" type="x-plane" boundary="reflective" coeffs="0.63"/>
<surface id="5" name="bottom" type="y-plane" boundary="reflective" coeffs="-0.63"/>
<surface id="6" name="top" type="y-plane" boundary="reflective" coeffs="0.63"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
</space>
<constraints>
<fissionable>true</fissionable>
</constraints>
</source>
<energy_mode>multi-group</energy_mode>
<temperature_method>nearest</temperature_method>
<temperature_multipole>true</temperature_multipole>
<temperature_range>200.0 400.0</temperature_range>
<temperature_tolerance>200.0</temperature_tolerance>
<random_ray>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
</space>
</source>
<distance_inactive>30.0</distance_inactive>
<distance_active>150.0</distance_active>
<source_region_meshes>
<mesh id="1">
<domain id="0" type="universe"/>
</mesh>
</source_region_meshes>
<source_shape>linear</source_shape>
</random_ray>
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-0.63 -0.63</lower_left>
<upper_right>0.63 0.63</upper_right>
</mesh>
</settings>
</model>

View file

@ -0,0 +1,2 @@
k-combined:
7.367927E-01 6.850805E-03

View file

@ -0,0 +1,70 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="UO2__2_4__" depletable="true" temperature="294.0">
<density value="1.0" units="macro"/>
<macroscopic name="UO2__2_4__"/>
</material>
<material id="2" name="Zircaloy" temperature="294.0">
<density value="1.0" units="macro"/>
<macroscopic name="Zircaloy"/>
</material>
<material id="3" name="Hot_borated_water" temperature="294.0">
<density value="1.0" units="macro"/>
<macroscopic name="Hot_borated_water"/>
</material>
</materials>
<geometry>
<cell id="1" name="Fuel" material="1" region="-1" universe="0">
<temperature>395.0</temperature>
</cell>
<cell id="2" name="Cladding" material="2" region="1 -2" universe="0"/>
<cell id="3" name="Water" material="3" region="2 3 -4 5 -6" universe="0"/>
<surface id="1" name="Fuel OR" type="z-cylinder" coeffs="0 0 0.39218"/>
<surface id="2" name="Clad OR" type="z-cylinder" coeffs="0 0 0.4572"/>
<surface id="3" name="left" type="x-plane" boundary="reflective" coeffs="-0.63"/>
<surface id="4" name="right" type="x-plane" boundary="reflective" coeffs="0.63"/>
<surface id="5" name="bottom" type="y-plane" boundary="reflective" coeffs="-0.63"/>
<surface id="6" name="top" type="y-plane" boundary="reflective" coeffs="0.63"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
</space>
<constraints>
<fissionable>true</fissionable>
</constraints>
</source>
<energy_mode>multi-group</energy_mode>
<temperature_method>nearest</temperature_method>
<temperature_multipole>true</temperature_multipole>
<temperature_range>200.0 400.0</temperature_range>
<temperature_tolerance>200.0</temperature_tolerance>
<random_ray>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
</space>
</source>
<distance_inactive>30.0</distance_inactive>
<distance_active>150.0</distance_active>
<source_region_meshes>
<mesh id="1">
<domain id="0" type="universe"/>
</mesh>
</source_region_meshes>
<source_shape>linear</source_shape>
</random_ray>
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-0.63 -0.63</lower_left>
<upper_right>0.63 0.63</upper_right>
</mesh>
</settings>
</model>

View file

@ -0,0 +1,2 @@
k-combined:
6.431774E-01 2.076589E-02

View file

@ -0,0 +1,73 @@
import os
import openmc
from openmc.examples import pwr_pin_cell
from openmc import RegularMesh
from openmc.utility_funcs import change_directory
import pytest
from tests.testing_harness import TolerantPyAPITestHarness
class MGXSTestHarness(TolerantPyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = 'mgxs.h5'
if os.path.exists(f):
os.remove(f)
@pytest.mark.parametrize("method", ["material_wise", "stochastic_slab", "infinite_medium"])
def test_random_ray_auto_convert(method):
with change_directory(method):
openmc.reset_auto_ids()
# Start with a normal continuous energy model
model = pwr_pin_cell()
temp_settings = {
'method' : 'nearest',
'tolerance' : 200.0,
'range' : (200.0, 400.0),
'multipole' : True
}
# Convert to a multi-group model
model.convert_to_multigroup(
method=method, groups='CASMO-2', nparticles=100,
overwrite_mgxs_library=False, mgxs_path="mgxs.h5",
temperatures=[294.0, 394.0], temperature_settings=temp_settings
)
# Convert to a random ray model
model.convert_to_random_ray()
model.settings.temperature = temp_settings
# Set all material temperatures to room temperature
for mat in model.geometry.get_all_materials().values():
mat.temperature = 294.0
# Set the cell temperature of the fuel such that it moves up to the next
# temperature bin.
for cell in model.geometry.get_all_cells().values():
if cell.name == "Fuel":
cell.temperature = [395.0]
# Set the number of particles
model.settings.particles = 100
# Overlay a basic 2x2 mesh
n = 2
mesh = RegularMesh()
mesh.dimension = (n, n)
bbox = model.geometry.bounding_box
mesh.lower_left = (bbox.lower_left[0], bbox.lower_left[1])
mesh.upper_right = (bbox.upper_right[0], bbox.upper_right[1])
model.settings.random_ray['source_region_meshes'] = [
(mesh, [model.geometry.root_universe])]
# Set the source shape to linear
model.settings.random_ray['source_shape'] = 'linear'
harness = MGXSTestHarness('statepoint.10.h5', model)
harness.main()

View file

@ -22,7 +22,7 @@ def test_random_ray_basic(run_mode):
if run_mode == "eigen":
openmc.reset_auto_ids()
model = random_ray_lattice()
# Double the densities of the lower-left fuel pin -> cell instances [0, 9).
# Double the densities of the lower-left fuel pin -> cell instances [0, 8).
for id, cell in model.geometry.get_all_cells().items():
if cell.fill.name == "UO2 fuel":
cell.density = [((i < 8) + 1.0) for i in range(24)]

View file

@ -0,0 +1,118 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="UO2 fuel">
<density value="1.0" units="macro"/>
<macroscopic name="UO2"/>
</material>
<material id="2" name="Water">
<density value="1.0" units="macro"/>
<macroscopic name="LWTR"/>
</material>
</materials>
<geometry>
<cell id="1" name="fuel inner a" material="1" region="-2" universe="1">
<temperature>395.0 395.0 395.0 395.0 395.0 395.0 395.0 395.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0</temperature>
</cell>
<cell id="2" name="fuel inner b" material="1" region="2 -3" universe="1">
<temperature>395.0 395.0 395.0 395.0 395.0 395.0 395.0 395.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0</temperature>
</cell>
<cell id="3" name="fuel inner c" material="1" region="3 -1" universe="1">
<temperature>395.0 395.0 395.0 395.0 395.0 395.0 395.0 395.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0 295.0</temperature>
</cell>
<cell id="4" name="moderator inner a" material="2" region="1 -4" universe="1"/>
<cell id="5" name="moderator outer b" material="2" region="4 -5" universe="1"/>
<cell id="6" name="moderator outer c" material="2" region="5" universe="1"/>
<cell id="7" name="azimuthal_cell_0" fill="1" region="6 -7" universe="2"/>
<cell id="8" name="azimuthal_cell_1" fill="1" region="7 -8" universe="2"/>
<cell id="9" name="azimuthal_cell_2" fill="1" region="8 -9" universe="2"/>
<cell id="10" name="azimuthal_cell_3" fill="1" region="9 -10" universe="2"/>
<cell id="11" name="azimuthal_cell_4" fill="1" region="10 -11" universe="2"/>
<cell id="12" name="azimuthal_cell_5" fill="1" region="11 -12" universe="2"/>
<cell id="13" name="azimuthal_cell_6" fill="1" region="12 -13" universe="2"/>
<cell id="14" name="azimuthal_cell_7" fill="1" region="13 -6" universe="2"/>
<cell id="15" name="moderator infinite" material="2" universe="3"/>
<cell id="16" fill="4" universe="5"/>
<cell id="17" name="assembly" fill="6" region="14 -15 16 -17" universe="7"/>
<lattice id="4">
<pitch>0.126 0.126</pitch>
<dimension>10 10</dimension>
<lower_left>-0.63 -0.63</lower_left>
<universes>
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 </universes>
</lattice>
<lattice id="6">
<pitch>1.26 1.26</pitch>
<dimension>2 2</dimension>
<lower_left>-1.26 -1.26</lower_left>
<universes>
2 2
2 5 </universes>
</lattice>
<surface id="1" name="Fuel OR" type="z-cylinder" coeffs="0.0 0.0 0.54"/>
<surface id="2" name="inner ring a" type="z-cylinder" coeffs="0.0 0.0 0.33"/>
<surface id="3" name="inner ring b" type="z-cylinder" coeffs="0.0 0.0 0.45"/>
<surface id="4" name="outer ring a" type="z-cylinder" coeffs="0.0 0.0 0.6"/>
<surface id="5" name="outer ring b" type="z-cylinder" coeffs="0.0 0.0 0.69"/>
<surface id="6" type="plane" coeffs="-0.0 1.0 0 0"/>
<surface id="7" type="plane" coeffs="-0.7071067811865475 0.7071067811865476 0 0"/>
<surface id="8" type="plane" coeffs="-1.0 6.123233995736766e-17 0 0"/>
<surface id="9" type="plane" coeffs="-0.7071067811865476 -0.7071067811865475 0 0"/>
<surface id="10" type="plane" coeffs="-1.2246467991473532e-16 -1.0 0 0"/>
<surface id="11" type="plane" coeffs="0.7071067811865475 -0.7071067811865477 0 0"/>
<surface id="12" type="plane" coeffs="1.0 -1.8369701987210297e-16 0 0"/>
<surface id="13" type="plane" coeffs="0.7071067811865477 0.7071067811865474 0 0"/>
<surface id="14" name="minimum x" type="x-plane" boundary="reflective" coeffs="-1.26"/>
<surface id="15" name="maximum x" type="x-plane" boundary="reflective" coeffs="1.26"/>
<surface id="16" name="minimum y" type="y-plane" boundary="reflective" coeffs="-1.26"/>
<surface id="17" name="maximum y" type="y-plane" boundary="reflective" coeffs="1.26"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<energy_mode>multi-group</energy_mode>
<temperature_method>nearest</temperature_method>
<temperature_range>200.0 400.0</temperature_range>
<temperature_tolerance>10.0</temperature_tolerance>
<random_ray>
<distance_active>100.0</distance_active>
<distance_inactive>20.0</distance_inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
</random_ray>
</settings>
<tallies>
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-1.26 -1.26</lower_left>
<upper_right>1.26 1.26</upper_right>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="2" type="energy">
<bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
</filter>
<tally id="1" name="Mesh tally">
<filters>1 2</filters>
<scores>flux fission nu-fission</scores>
<estimator>analog</estimator>
</tally>
</tallies>
</model>

View file

@ -0,0 +1,171 @@
k-combined:
8.721099E-01 6.686066E-03
tally 1:
1.530055E+00
4.684582E-01
2.918032E-01
1.703772E-02
7.101906E-01
1.009209E-01
7.732582E-01
1.197247E-01
5.787213E-02
6.705930E-04
1.408492E-01
3.972179E-03
4.223488E-01
3.605283E-02
6.818063E-03
9.396342E-06
1.659380E-02
5.565812E-05
5.958551E-01
7.215230E-02
9.932163E-03
2.004773E-05
2.417290E-02
1.187504E-04
1.685106E+00
5.752729E-01
9.834826E-03
1.960236E-05
2.393629E-02
1.161151E-04
4.400457E+00
3.886563E+00
3.318560E-03
2.211175E-06
8.211544E-03
1.353859E-05
2.814971E+00
1.585202E+00
1.876190E-02
7.040316E-05
5.218528E-02
5.446713E-04
1.970406E+00
7.765020E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.643086E-01
1.494764E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.421562E-01
3.961223E-02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.421077E-01
8.385255E-02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.717683E+00
5.974687E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.024715E+00
3.251629E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.507885E+00
1.258277E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.237600E+00
3.066006E-01
4.632318E-01
4.294997E-02
1.127413E+00
2.544090E-01
7.085300E-01
1.005838E-01
1.066845E-01
2.280854E-03
2.596486E-01
1.351037E-02
4.080999E-01
3.359223E-02
1.334606E-02
3.590848E-05
3.248163E-02
2.126996E-04
5.587686E-01
6.339264E-02
1.868162E-02
7.084921E-05
4.546734E-02
4.196670E-04
1.647852E+00
5.503007E-01
1.940765E-02
7.635704E-05
4.723492E-02
4.523030E-04
4.672431E+00
4.381400E+00
7.228567E-03
1.048470E-05
1.788658E-02
6.419579E-05
3.054158E+00
1.866081E+00
4.202691E-02
3.534383E-04
1.168957E-01
2.734362E-03
1.313678E+00
3.454374E-01
4.964081E-01
4.934325E-02
1.208158E+00
2.922789E-01
7.298176E-01
1.067017E-01
1.106319E-01
2.452869E-03
2.692559E-01
1.452928E-02
4.129038E-01
3.441066E-02
1.357107E-02
3.713729E-05
3.302926E-02
2.199783E-04
5.677464E-01
6.543807E-02
1.908992E-02
7.390190E-05
4.646105E-02
4.377492E-04
1.651067E+00
5.522453E-01
1.956802E-02
7.754346E-05
4.762523E-02
4.593308E-04
4.583305E+00
4.217589E+00
7.135836E-03
1.022408E-05
1.765713E-02
6.260005E-05
2.988394E+00
1.786305E+00
4.141550E-02
3.431964E-04
1.151951E-01
2.655125E-03

View file

@ -0,0 +1,36 @@
import os
import openmc
from openmc.examples import random_ray_lattice, random_ray_three_region_cube
from openmc.utility_funcs import change_directory
import pytest
from tests.testing_harness import TolerantPyAPITestHarness
class MGXSTestHarness(TolerantPyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = 'mgxs.h5'
if os.path.exists(f):
os.remove(f)
def test_random_ray_basic():
openmc.reset_auto_ids()
model = random_ray_lattice(second_temp=True)
# Set the temperature of the lower-left pin to 395 K -> cell instances [0, 8).
# All other pins are set to 295.
for id, cell in model.geometry.get_all_cells().items():
if cell.fill.name == "UO2 fuel":
cell.temperature = [(100.0 * (i < 8) + 295.0) for i in range(24)]
model.settings.temperature = {
'method' : 'nearest',
'tolerance' : 10.0,
'range' : (200.0, 400.0)
}
# Gold file was generated with manually scaled fuel cross sections.
harness = MGXSTestHarness('statepoint.10.h5', model)
harness.main()