crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

New multiply_density option on tallies (#2539)

Browse source
This commit is contained in:
Paul Romano 2023-06-09 07:24:00 -05:00 • committed by GitHub
parent a288e42a06
commit c123ce0fa0
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
18 changed files with 286 additions and 48 deletions
Download patch file Download diff file Expand all files Collapse all files

2
docs/source/io_formats/statepoint.rst
View file

@ -139,6 +139,8 @@ The current version of the statepoint file format is 17.0.
- **internal** (*int*) -- Flag indicating the presence of tally
data (0) or absence of tally data (1). All user defined
tallies will have a value of 0 unless otherwise instructed.
- **multiply_density** (*int*) -- Flag indicating whether reaction
rates should be multiplied by atom density (1) or not (0).
:Datasets: - **n_realizations** (*int*) -- Number of realizations.
- **n_filters** (*int*) -- Number of filters used.

6
docs/source/io_formats/tallies.rst
View file

@ -69,6 +69,12 @@ The ``<tally>`` element accepts the following sub-elements:
list of valid scores can be found in the :ref:`user's guide
<usersguide_scores>`.
:multiply_density:
A boolean that indicates whether reaction rate scores should be computed by
multiplying by the atom density of a nuclide present in a material.
*Default*: true
:trigger:
Precision trigger applied to all filter bins and nuclides for this tally.
It must specify the trigger's type, threshold and scores to which it will

18
include/openmc/nuclide.h
View file

@ -29,6 +29,7 @@ namespace openmc {
class Nuclide {
public:
//============================================================================
// Types, aliases
using EmissionMode = ReactionProduct::EmissionMode;
struct EnergyGrid {
@ -36,18 +37,32 @@ public:
vector<double> energy;
};
//============================================================================
// Constructors/destructors
Nuclide(hid_t group, const vector<double>& temperature);
~Nuclide();
//============================================================================
// Methods
//! Initialize logarithmic grid for energy searches
void init_grid();
//! Calculate microscopic cross sections
//
//! \param[in] i_sab Index in data::thermal_scatt
//! \param[in] i_log_union Log-grid search index
//! \param[in] sab_frac S(a,b) table fraction
//! \param[in,out] p Particle object
void calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle& p);
//! Calculate thermal scattering cross section
//
//! \param[in] i_sab Index in data::thermal_scatt
//! \param[in] sab_frac S(a,b) table fraction
//! \param[in,out] p Particle object
void calculate_sab_xs(int i_sab, double sab_frac, Particle& p);
// Methods
double nu(double E, EmissionMode mode, int group = 0) const;
void calculate_elastic_xs(Particle& p) const;
@ -69,6 +84,7 @@ public:
double collapse_rate(int MT, double temperature,
gsl::span<const double> energy, gsl::span<const double> flux) const;
//============================================================================
// Data members
std::string name_; //!< Name of nuclide, e.g. "U235"
int Z_; //!< Atomic number

13
include/openmc/particle.h
View file

@ -35,6 +35,9 @@ public:
Particle() = default;
//==========================================================================
// Methods
double speed() const;
//! create a secondary particle
@ -106,6 +109,16 @@ public:
//! create a particle restart HDF5 file
void write_restart() const;
//! Update microscopic cross section cache
//
//! \param[in] i_nuclide Index in data::nuclides
//! \param[in] i_grid Index on log union grid
//! \param[in] i_sab Index in data::thermal_scatt
//! \param[in] sab_frac S(a,b) table fraction
//! \param[in] ncrystal_xs Thermal scattering xs from NCrystal
void update_neutron_xs(int i_nuclide, int i_grid = C_NONE, int i_sab = C_NONE,
double sab_frac = 0.0, double ncrystal_xs = -1.0);
};
//============================================================================

3
include/openmc/particle_data.h
View file

@ -218,8 +218,9 @@ struct BoundaryInfo {
* https://doi.org/10.1016/j.anucene.2017.11.032.
*/
class ParticleData {
public:
//----------------------------------------------------------------------------
// Constructors
ParticleData();
private:

7
include/openmc/tallies/tally.h
View file

@ -37,6 +37,8 @@ public:
void set_active(bool active) { active_ = active; }
void set_multiply_density(bool value) { multiply_density_ = value; }
void set_writable(bool writable) { writable_ = writable; }
void set_scores(pugi::xml_node node);
@ -62,6 +64,8 @@ public:
int32_t n_filter_bins() const { return n_filter_bins_; }
bool multiply_density() const { return multiply_density_; }
bool writable() const { return writable_; }
//----------------------------------------------------------------------------
@ -139,6 +143,9 @@ private:
int32_t n_filter_bins_ {0};
//! Whether to multiply by atom density for reaction rates
bool multiply_density_ {true};
gsl::index index_;
};

20
openmc/lib/tally.py
View file

@ -39,6 +39,9 @@ _dll.openmc_tally_get_filters.argtypes = [
c_int32, POINTER(POINTER(c_int32)), POINTER(c_size_t)]
_dll.openmc_tally_get_filters.restype = c_int
_dll.openmc_tally_get_filters.errcheck = _error_handler
_dll.openmc_tally_get_multiply_density.argtypes = [c_int32, POINTER(c_bool)]
_dll.openmc_tally_get_multiply_density.restype = c_int
_dll.openmc_tally_get_multiply_density.errcheck = _error_handler
_dll.openmc_tally_get_n_realizations.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_n_realizations.restype = c_int
_dll.openmc_tally_get_n_realizations.errcheck = _error_handler
@ -75,6 +78,9 @@ _dll.openmc_tally_set_estimator.errcheck = _error_handler
_dll.openmc_tally_set_id.argtypes = [c_int32, c_int32]
_dll.openmc_tally_set_id.restype = c_int
_dll.openmc_tally_set_id.errcheck = _error_handler
_dll.openmc_tally_set_multiply_density.argtypes = [c_int32, c_bool]
_dll.openmc_tally_set_multiply_density.restype = c_int
_dll.openmc_tally_set_multiply_density.errcheck = _error_handler
_dll.openmc_tally_set_nuclides.argtypes = [c_int32, c_int, POINTER(c_char_p)]
_dll.openmc_tally_set_nuclides.restype = c_int
_dll.openmc_tally_set_nuclides.errcheck = _error_handler
@ -174,6 +180,10 @@ class Tally(_FortranObjectWithID):
List of tally filters
mean : numpy.ndarray
An array containing the sample mean for each bin
multiply_density : bool
Whether reaction rates should be multiplied by atom density
.. versionadded:: 0.13.4
nuclides : list of str
List of nuclides to score results for
num_realizations : int
@ -363,6 +373,16 @@ class Tally(_FortranObjectWithID):
def writable(self, writable):
_dll.openmc_tally_set_writable(self._index, writable)
@property
def multiply_density(self):
multiply_density = c_bool()
_dll.openmc_tally_get_multiply_density(self._index, multiply_density)
return multiply_density.value
@multiply_density.setter
def multiply_density(self, multiply_density):
_dll.openmc_tally_set_multiply_density(self._index, multiply_density)
def reset(self):
"""Reset results and num_realizations of tally"""
_dll.openmc_tally_reset(self._index)

4
openmc/statepoint.py
View file

@ -398,6 +398,10 @@ class StatePoint:
tally._sp_filename = self._f.filename
tally.name = group['name'][()].decode() if 'name' in group else ''
# Check if tally has multiply_density attribute
if "multiply_density" in group.attrs:
tally.multiply_density = group.attrs["multiply_density"].item() > 0
# Read the number of realizations
n_realizations = group['n_realizations'][()]

25
openmc/tallies.py
View file

@ -14,7 +14,7 @@ import scipy.sparse as sps
import openmc
import openmc.checkvalue as cv
from ._xml import clean_indentation, reorder_attributes
from ._xml import clean_indentation, reorder_attributes, get_text
from .mixin import IDManagerMixin
from .mesh import MeshBase
@ -54,6 +54,10 @@ class Tally(IDManagerMixin):
Unique identifier for the tally
name : str
Name of the tally
multiply_density : bool
Whether reaction rates should be multiplied by atom density
.. versionadded:: 0.13.4
filters : list of openmc.Filter
List of specified filters for the tally
nuclides : list of str
@ -111,6 +115,7 @@ class Tally(IDManagerMixin):
self._estimator = None
self._triggers = cv.CheckedList(openmc.Trigger, 'tally triggers')
self._derivative = None
self._multiply_density = True
self._num_realizations = 0
self._with_summary = False
@ -138,12 +143,17 @@ class Tally(IDManagerMixin):
parts.append('{: <15}=\t{}'.format('Nuclides', nuclides))
parts.append('{: <15}=\t{}'.format('Scores', self.scores))
parts.append('{: <15}=\t{}'.format('Estimator', self.estimator))
parts.append('{: <15}=\t{}'.format('Multiply dens.', self.multiply_density))
return '\n\t'.join(parts)
@property
def name(self):
return self._name
@property
def multiply_density(self):
return self._multiply_density
@property
def filters(self):
return self._filters
@ -323,6 +333,11 @@ class Tally(IDManagerMixin):
cv.check_type('tally name', name, str, none_ok=True)
self._name = name
@multiply_density.setter
def multiply_density(self, value):
cv.check_type('multiply density', value, bool)
self._multiply_density = value
@derivative.setter
def derivative(self, deriv):
cv.check_type('tally derivative', deriv, openmc.TallyDerivative,
@ -829,6 +844,10 @@ class Tally(IDManagerMixin):
if self.name != '':
element.set("name", self.name)
# Multiply by density
if not self.multiply_density:
element.set("multiply_density", str(self.multiply_density).lower())
# Optional Tally filters
if len(self.filters) > 0:
subelement = ET.SubElement(element, "filters")
@ -885,6 +904,10 @@ class Tally(IDManagerMixin):
name = elem.get('name', '')
tally = cls(tally_id=tally_id, name=name)
text = get_text(elem, 'multiply_density')
if text is not None:
tally.multiply_density = text in ('true', '1')
# Read filters
filters_elem = elem.find('filters')
if filters_elem is not None:

19
src/material.cpp
View file

@ -353,7 +353,7 @@ Material::~Material()
model::material_map.erase(id_);
}
Material & Material::clone()
Material& Material::clone()
{
std::unique_ptr<Material> mat = std::make_unique<Material>();
@ -868,21 +868,12 @@ void Material::calculate_neutron_xs(Particle& p) const
// ======================================================================
// CALCULATE MICROSCOPIC CROSS SECTION
// Determine microscopic cross sections for this nuclide
// Get nuclide index
int i_nuclide = nuclide_[i];
// Calculate microscopic cross section for this nuclide
auto& micro {p.neutron_xs(i_nuclide)};
if (p.E() != micro.last_E || p.sqrtkT() != micro.last_sqrtkT ||
i_sab != micro.index_sab || sab_frac != micro.sab_frac) {
data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p);
// If NCrystal is being used, update micro cross section cache
if (ncrystal_xs >= 0.0) {
data::nuclides[i_nuclide]->calculate_elastic_xs(p);
ncrystal_update_micro(ncrystal_xs, micro);
}
}
// Update microscopic cross section for this nuclide
p.update_neutron_xs(i_nuclide, i_grid, i_sab, sab_frac, ncrystal_xs);
auto& micro = p.neutron_xs(i_nuclide);
// ======================================================================
// ADD TO MACROSCOPIC CROSS SECTION

2
src/nuclide.cpp
View file

@ -471,7 +471,7 @@ void Nuclide::create_derived(
xs_cdf_sum +=
(std::sqrt(E[i]) * xs[i] + std::sqrt(E[i + 1]) * xs[i + 1]) / 2.0 *
(E[i + 1] - E[i]);
xs_cdf_[i+1] = xs_cdf_sum;
xs_cdf_[i + 1] = xs_cdf_sum;
}
}
}

31
src/particle.cpp
View file

@ -36,6 +36,10 @@
namespace openmc {
//==============================================================================
// Particle implementation
//==============================================================================
double Particle::speed() const
{
// Determine mass in eV/c^2
@ -433,7 +437,9 @@ void Particle::cross_surface()
#ifdef DAGMC
// in DAGMC, we know what the next cell should be
if (surf->geom_type_ == GeometryType::DAG) {
int32_t i_cell = next_cell(i_surface, cell_last(n_coord() - 1), lowest_coord().universe) - 1;
int32_t i_cell =
next_cell(i_surface, cell_last(n_coord() - 1), lowest_coord().universe) -
1;
// save material and temp
material_last() = material();
sqrtkT_last() = sqrtkT();
@ -703,6 +709,29 @@ void Particle::write_restart() const
} // #pragma omp critical
}
void Particle::update_neutron_xs(
int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)
{
// Get microscopic cross section cache
auto& micro = this->neutron_xs(i_nuclide);
// If the cache doesn't match, recalculate micro xs
if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||
i_sab != micro.index_sab || sab_frac != micro.sab_frac) {
data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);
// If NCrystal is being used, update micro cross section cache
if (ncrystal_xs >= 0.0) {
data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
ncrystal_update_micro(ncrystal_xs, micro);
}
}
}
//==============================================================================
// Non-method functions
//==============================================================================
std::string particle_type_to_str(ParticleType type)
{
switch (type) {

6
src/state_point.cpp
View file

@ -193,6 +193,12 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
continue;
}
if (tally->multiply_density()) {
write_attribute(tally_group, "multiply_density", 1);
} else {
write_attribute(tally_group, "multiply_density", 0);
}
if (tally->estimator_ == TallyEstimator::ANALOG) {
write_dataset(tally_group, "estimator", "analog");
} else if (tally->estimator_ == TallyEstimator::TRACKLENGTH) {

37
src/tallies/tally.cpp
View file

@ -92,6 +92,10 @@ Tally::Tally(pugi::xml_node node)
if (check_for_node(node, "name"))
name_ = get_node_value(node, "name");
if (check_for_node(node, "multiply_density")) {
multiply_density_ = get_node_value_bool(node, "multiply_density");
}
// =======================================================================
// READ DATA FOR FILTERS
@ -564,11 +568,12 @@ void Tally::set_nuclides(const vector<std::string>& nuclides)
nuclides_.push_back(-1);
} else {
auto search = data::nuclide_map.find(nuc);
if (search == data::nuclide_map.end())
fatal_error(fmt::format("Could not find the nuclide {} specified in "
"tally {} in any material",
nuc, id_));
nuclides_.push_back(search->second);
if (search == data::nuclide_map.end()) {
int err = openmc_load_nuclide(nuc.c_str(), nullptr, 0);
if (err < 0)
throw std::runtime_error {openmc_err_msg};
}
nuclides_.push_back(data::nuclide_map.at(nuc));
}
}
}
@ -1108,6 +1113,28 @@ extern "C" int openmc_tally_set_writable(int32_t index, bool writable)
return 0;
}
extern "C" int openmc_tally_get_multiply_density(int32_t index, bool* value)
{
if (index < 0 || index >= model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*value = model::tallies[index]->multiply_density();
return 0;
}
extern "C" int openmc_tally_set_multiply_density(int32_t index, bool value)
{
if (index < 0 || index >= model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
model::tallies[index]->set_multiply_density(value);
return 0;
}
extern "C" int openmc_tally_get_scores(int32_t index, int** scores, int* n)
{
if (index < 0 || index >= model::tallies.size()) {

76
src/tallies/tally_scoring.cpp
View file

@ -939,21 +939,17 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
if (i_nuclide >= 0) {
const auto& micro = p.photon_xs(i_nuclide);
double xs = (score_bin == COHERENT)
? micro.coherent
: (score_bin == INCOHERENT) ? micro.incoherent
: (score_bin == PHOTOELECTRIC)
? micro.photoelectric
: micro.pair_production;
double xs = (score_bin == COHERENT) ? micro.coherent
: (score_bin == INCOHERENT) ? micro.incoherent
: (score_bin == PHOTOELECTRIC) ? micro.photoelectric
: micro.pair_production;
score = xs * atom_density * flux;
} else {
double xs = (score_bin == COHERENT)
? p.macro_xs().coherent
: (score_bin == INCOHERENT)
? p.macro_xs().incoherent
: (score_bin == PHOTOELECTRIC)
? p.macro_xs().photoelectric
: p.macro_xs().pair_production;
double xs = (score_bin == COHERENT) ? p.macro_xs().coherent
: (score_bin == INCOHERENT) ? p.macro_xs().incoherent
: (score_bin == PHOTOELECTRIC)
? p.macro_xs().photoelectric
: p.macro_xs().pair_production;
score = xs * flux;
}
break;
@ -2308,6 +2304,9 @@ void score_tracklength_tally(Particle& p, double distance)
// Determine the tracklength estimate of the flux
double flux = p.wgt() * distance;
// Set 'none' value for log union grid index
int i_log_union = C_NONE;
for (auto i_tally : model::active_tracklength_tallies) {
const Tally& tally {*model::tallies[i_tally]};
@ -2331,11 +2330,25 @@ void score_tracklength_tally(Particle& p, double distance)
double atom_density = 0.;
if (i_nuclide >= 0) {
if (p.material() != MATERIAL_VOID) {
auto j =
model::materials[p.material()]->mat_nuclide_index_[i_nuclide];
if (j == C_NONE)
continue;
atom_density = model::materials[p.material()]->atom_density_(j);
const auto& mat = model::materials[p.material()];
auto j = mat->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) {
// Determine log union grid index
if (i_log_union == C_NONE) {
int neutron = static_cast<int>(ParticleType::neutron);
i_log_union = std::log(p.E() / data::energy_min[neutron]) /
simulation::log_spacing;
}
// Update micro xs cache
if (!tally.multiply_density()) {
p.update_neutron_xs(i_nuclide, i_log_union);
atom_density = 1.0;
}
} else {
atom_density =
tally.multiply_density() ? mat->atom_density_(j) : 1.0;
}
}
}
@ -2371,6 +2384,9 @@ void score_collision_tally(Particle& p)
flux = p.wgt_last() / p.macro_xs().total;
}
// Set 'none value for log union grid index
int i_log_union = C_NONE;
for (auto i_tally : model::active_collision_tallies) {
const Tally& tally {*model::tallies[i_tally]};
@ -2393,11 +2409,25 @@ void score_collision_tally(Particle& p)
double atom_density = 0.;
if (i_nuclide >= 0) {
auto j =
model::materials[p.material()]->mat_nuclide_index_[i_nuclide];
if (j == C_NONE)
continue;
atom_density = model::materials[p.material()]->atom_density_(j);
const auto& mat = model::materials[p.material()];
auto j = mat->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) {
// Determine log union grid index
if (i_log_union == C_NONE) {
int neutron = static_cast<int>(ParticleType::neutron);
i_log_union = std::log(p.E() / data::energy_min[neutron]) /
simulation::log_spacing;
}
// Update micro xs cache
if (!tally.multiply_density()) {
p.update_neutron_xs(i_nuclide, i_log_union);
atom_density = 1.0;
}
} else {
atom_density =
tally.multiply_density() ? mat->atom_density_(j) : 1.0;
}
}
// TODO: consider replacing this "if" with pointers or templates

2
tests/regression_tests/tallies/results_true.dat
View file

@ -1 +1 @@
6b1d8d6f4d7a70af6c39cc76b9267a61ba9a9d0c14b75a4df6f83e72a2bfaf1533caaf936c2185f0a158443eb9da266bd5a77b7996020fd638551ea841d821e0
d1decdbec6cb59df91ba5c42cb37a04f413a34fa7faf7ad1eecfd7d53a14af8cb65b58335c4ede4e88f6d9ab35a1746251983cc991d74eab3e678887c84183bd

13
tests/unit_tests/test_lib.py
View file

@ -366,6 +366,19 @@ def test_tally_activate(lib_simulation_init):
assert t.active
def test_tally_multiply_density(lib_simulation_init):
# multiply_density is True by default
t = openmc.lib.tallies[1]
assert t.multiply_density
# Make sure setting multiply_density works
t.multiply_density = False
assert not t.multiply_density
# Reset to True
t.multiply_density = True
def test_tally_writable(lib_simulation_init):
t = openmc.lib.tallies[1]
assert t.writable

50
tests/unit_tests/test_tally_multiply_density.py Normal file
View file

@ -0,0 +1,50 @@
import numpy as np
import openmc
import pytest
def test_micro_macro_compare():
# Create simple sphere model with H1 and H2
mat = openmc.Material()
mat.add_components({'H1': 1.0, 'H2': 1.0})
mat.set_density('g/cm3', 1.0)
sph = openmc.Sphere(r=10.0, boundary_type='vacuum')
cell = openmc.Cell(fill=mat, region=-sph)
model = openmc.Model()
model.geometry = openmc.Geometry([cell])
model.settings.run_mode = 'fixed source'
model.settings.particles = 1000
model.settings.batches = 10
# Set up two reaction rate tallies, one that multplies by density and the
# other that doesn't
tally_macro = openmc.Tally()
tally_macro.nuclides = ['H1', 'H2', 'H3']
tally_macro.scores = ['total', 'elastic']
tally_micro = openmc.Tally()
tally_micro.nuclides = ['H1', 'H2', 'H3']
tally_micro.scores = ['total', 'elastic']
tally_micro.multiply_density = False
model.tallies = [tally_macro, tally_micro]
sp_filename = model.run()
with openmc.StatePoint(sp_filename) as sp:
tally_macro = sp.tallies[tally_macro.id]
tally_micro = sp.tallies[tally_micro.id]
# Make sure multply_density attribute from statepoint is set correctly
assert tally_macro.multiply_density
assert not tally_micro.multiply_density
# Dividing macro by density should give micro
density = mat.get_nuclide_atom_densities()
for nuc in ('H1', 'H2'):
micro_derived = tally_macro.get_values(nuclides=[nuc]) / density[nuc]
micro = tally_micro.get_values(nuclides=[nuc])
assert micro_derived == pytest.approx(micro)
# For macro tally, H3 scores should be zero
assert np.all(tally_macro.get_values(nuclides=['H3']) == 0.0)
# For micro tally, H3 scores should be positive
assert np.all(tally_micro.get_values(nuclides=['H3']) > 0.0)