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

Merge c217ba86b3 into d3bc1669d3

Browse source
This commit is contained in:
Vitaly Mogulian 2026-07-17 20:51:40 +08:00 • committed by GitHub
commit c5ddd408e8
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
11 changed files with 765 additions and 327 deletions
Download patch file Download diff file Expand all files Collapse all files

2
include/openmc/message_passing.h
View file

@ -19,6 +19,8 @@ extern bool master;
#ifdef OPENMC_MPI
extern MPI_Datatype source_site;
extern MPI_Datatype collision_track_site;
extern MPI_Datatype volume_results;
extern MPI_Datatype volume_tally;
extern MPI_Comm intracomm;
#endif

193
include/openmc/volume_calc.h
View file

@ -7,16 +7,15 @@
#include <vector>
#include "openmc/array.h"
#include "openmc/cell.h"
#include "openmc/openmp_interface.h"
#include "openmc/position.h"
#include "openmc/tallies/trigger.h"
#include "openmc/timer.h"
#include "openmc/vector.h"
#include "openmc/tensor.h"
#include "pugixml.hpp"
#ifdef _OPENMP
#include <omp.h>
#endif
namespace openmc {
@ -25,60 +24,171 @@ namespace openmc {
//==============================================================================
class VolumeCalculation {
public:
private:
// Aliases, types
struct Result {
array<double, 2> volume; //!< Mean/standard deviation of volume
struct NuclResult {
vector<int> nuclides; //!< Index of nuclides
vector<double> atoms; //!< Number of atoms for each nuclide
vector<double> uncertainty; //!< Uncertainty on number of atoms
int iterations; //!< Number of iterations needed to obtain the results
}; // Results for a single domain
}; // Results of nuclides calculation for a single domain
public:
//! \brief Tally corresponding to a single material (AoS)
struct VolTally {
double score {0.0}; //!< Current batch scores accumulator
array<double, 2> score_acc {}; //!< Scores and squared scores accumulator
int32_t index {0}; //!< Material ID
VolTally() = default;
inline VolTally(int i, double s = 0.0)
{
index = i;
score = s;
}
//! \brief Add batch scores means to a tally
//! \param[in] batch_size_1 Inversed batch size
inline void finalize_batch(const double batch_size_1);
//! \brief Pass counters data from a given tally to this
//! \param[in] vol_tally Data source
inline void assign_tally(const VolTally& vol_tally);
//! \brief Add counters data from a given tally to this
//! \param[in] vol_tally Data source
inline void append_tally(const VolTally& vol_tally);
//! \brief Determines given trigger condition satisfaction for this tally
//
//! \param[in] trigger_type Type of trigger condition
//! \param[in] threshold Value for trigger condition (either volume
//! fraction variance or squared rel. err. dependent on the trigger type)
//! \param[in] n_samples Statistics size
//! \return True if the trigger condition is satisfied
inline bool trigger_state(const TriggerMetric trigger_type,
const double threshold, const size_t& n_samples) const;
}; // public just because it is used externally for the MPI struct definition
//! \brief Online results of calculation specific for each thread
struct CalcResults {
uint64_t n_samples; //!< Number of samples
int iterations; //!< Number of iterations needed to obtain the results
double cost; //!< Product of spent time and number of threads/processes
Timer sampling_time; // Timer for measurment of the simulation
vector<vector<VolumeCalculation::VolTally>>
vol_tallies; //!< Volume tallies for each domain
vector<VolumeCalculation::NuclResult>
nuc_results; //!< Nuclides of each domain
CalcResults(const VolumeCalculation& vol_calc);
CalcResults() = default;
//! \brief Reset all counters
void reset();
//! \brief Append another counters to this
void append(const CalcResults& other);
#ifdef OPENMC_MPI
//! \brief Collects results from all MPI processes to this
void collect_MPI();
#endif
};
// Constructors
VolumeCalculation(pugi::xml_node node);
VolumeCalculation() = default;
// Methods
//! \brief Stochastically determine the volume of a set of domains along with
//! the
//! average number densities of nuclides within the domain
//! the average number densities of nuclides within the domain
//
//! \param[in,out] results Results of calculation entity for filling
//! \return Vector of results for each user-specified domain
vector<Result> execute() const;
void execute(CalcResults& results) const;
//! \brief Print volume calculation results
//
//! \param[in] results Full volume calculation results
void show_results(const CalcResults& results) const;
//! \brief Write volume calculation results to HDF5 file
//
//! \param[in] filename Path to HDF5 file to write
//! \param[in] results Vector of results for each domain
void to_hdf5(
const std::string& filename, const vector<Result>& results) const;
//! \param[in] results Results entity
void to_hdf5(const std::string& filename, const CalcResults& results) const;
private:
//! \brief Rejection estimator
inline void score_hit(const Particle& p, CalcResults& results) const;
//! \brief Check whether a material has already been hit for a given domain.
//! If not, add new entries to the vectors
//
//! \param[in] i_material Index in global materials vector
//! \param[in] contrib Scoring value
//! \param[in,out] vol_tallies Vector of tallies corresponding to each
//! material
void check_hit(const int32_t i_material, const double contrib,
vector<VolTally>& vol_tallies) const;
//! \brief Reduce vector of volumetric tallies from each thread to a single
//! copy
//
//! \param[in] local_results Results specific to each thread
//! \param[out] results Reduced results
void reduce_results(
const CalcResults& local_results, CalcResults& results) const;
//! \brief Prints a statistics parameter
//
//! \param[in] label Name of parameter
//! \param[in] units Units of measure
//! \param[in] value Value of parameter
void show_vol_stat(
const std::string label, const std::string units, const double value) const;
//! \brief Prints volume result for a domain
//
//! \param[in] domain_type Either material or cell, ect.
//! \param[in] domain_id Number of this domain
//! \param[in] region_name Domain description
//! \param[in] mean Center of confidence interval
//! \param[in] stddev Half-width of confidence interval
void show_volume(const std::string domain_type, const int domain_id,
const std::string region_name, const double mean,
const double stddev) const;
// Tally filter and map types
enum class TallyDomain { UNIVERSE, MATERIAL, CELL };
// Data members
TallyDomain domain_type_; //!< Type of domain (cell, material, etc.)
size_t n_samples_; //!< Number of samples to use
double threshold_ {-1.0}; //!< Error threshold for domain volumes
TallyDomain domain_type_; //!< Type of domain (cell, material, etc.)
size_t n_samples_; //!< Number of samples to use
double volume_sample_; //!< Volume of bounding primitive
double threshold_ {-1.0}; //!< Error threshold for domain volumes
double threshold_cnd_; //!< Pre-computed value for trigger condition
int max_iterations_ {INT_MAX}; //!< Limit of iterations number (necessary
//!< maximum value of data type by default)
TriggerMetric trigger_type_ {
TriggerMetric::not_active}; //!< Trigger metric for the volume calculation
Position lower_left_; //!< Lower-left position of bounding box
Position upper_right_; //!< Upper-right position of bounding box
vector<int> domain_ids_; //!< IDs of domains to find volumes of
private:
//! \brief Check whether a material has already been hit for a given domain.
//! If not, add new entries to the vectors
constexpr static int _INDEX_TOTAL =
-999; //!< Index of zero-element tally for entire domain totals should be
//!< out of material ID range
//! \brief Computes estimated mean and std.dev. for a tally
//
//! \param[in] i_material Index in global materials vector
//! \param[in,out] indices Vector of material indices
//! \param[in,out] hits Number of hits corresponding to each material
void check_hit(
int i_material, vector<uint64_t>& indices, vector<uint64_t>& hits) const;
//! \param[in] n_samples Statistic's size
//! \param[in] coeff_norm Normalization coefficient to multiply
//! \param[in] vol_tally Tally
//! \return Array of mean and stddev
array<double, 2> get_tally_results(const size_t& n_samples,
const double coeff_norm, const VolTally& vol_tally) const;
};
//==============================================================================
@ -93,35 +203,6 @@ extern vector<VolumeCalculation> volume_calcs;
// Non-member functions
//==============================================================================
//! Reduce vector of indices and hits from each thread to a single copy
//
//! \param[in] local_indices Indices specific to each thread
//! \param[in] local_hits Hit count specific to each thread
//! \param[out] indices Reduced vector of indices
//! \param[out] hits Reduced vector of hits
template<typename T, typename T2>
void reduce_indices_hits(const vector<T>& local_indices,
const vector<T2>& local_hits, vector<T>& indices, vector<T2>& hits)
{
const int n_threads = num_threads();
#pragma omp for ordered schedule(static, 1)
for (int i = 0; i < n_threads; ++i) {
#pragma omp ordered
for (int j = 0; j < local_indices.size(); ++j) {
// Check if this material has been added to the master list and if
// so, accumulate the number of hits
auto it = std::find(indices.begin(), indices.end(), local_indices[j]);
if (it == indices.end()) {
indices.push_back(local_indices[j]);
hits.push_back(local_hits[j]);
} else {
hits[it - indices.begin()] += local_hits[j];
}
}
}
}
void free_memory_volume();
} // namespace openmc

31
openmc/volume.py
View file

@ -66,6 +66,11 @@ class VolumeCalculation:
Number of iterations over samples (for calculations with a trigger).
.. versionadded:: 0.12
max_iterations : int
Limit of the maximal allowed iterations number (optional, for
calculations with a trigger).
.. versionadded:: 0.15.x
"""
def __init__(self, domains, samples, lower_left=None, upper_right=None):
@ -73,6 +78,7 @@ class VolumeCalculation:
self._volumes = {}
self._threshold = None
self._trigger_type = None
self._max_iterations = None
self._iterations = None
cv.check_type('domains', domains, Iterable,
@ -187,6 +193,18 @@ class VolumeCalculation:
('variance', 'std_dev', 'rel_err'))
self._trigger_type = trigger_type
@property
def max_iterations(self):
return self._max_iterations
@max_iterations.setter
def max_iterations(self, max_iterations):
name = 'volume calculation iterations limit'
cv.check_type(name, max_iterations, Integral, none_ok=True)
if max_iterations is not None:
cv.check_greater_than(name, max_iterations, 0)
self._max_iterations = max_iterations
@property
def iterations(self):
return self._iterations
@ -230,7 +248,7 @@ class VolumeCalculation:
return pd.DataFrame.from_records(items, columns=columns)
def set_trigger(self, threshold, trigger_type):
def set_trigger(self, threshold, trigger_type, max_iterations=None):
"""Set a trigger on the volume calculation
.. versionadded:: 0.12
@ -241,9 +259,12 @@ class VolumeCalculation:
Threshold for the maximum standard deviation of volumes
trigger_type : {'variance', 'std_dev', 'rel_err'}
Value type used to halt volume calculation
max_iterations : int
Maximal allowed number of iterations (optional)
"""
self.trigger_type = trigger_type
self.threshold = threshold
self.max_iterations = max_iterations
@classmethod
def from_hdf5(cls, filename):
@ -270,6 +291,7 @@ class VolumeCalculation:
threshold = f.attrs.get('threshold')
trigger_type = f.attrs.get('trigger_type')
max_iterations = f.attrs.get('max_iterations')
iterations = f.attrs.get('iterations', 1)
volumes = {}
@ -304,7 +326,7 @@ class VolumeCalculation:
vol = cls(domains, samples, lower_left, upper_right)
if trigger_type is not None:
vol.set_trigger(threshold, trigger_type.decode())
vol.set_trigger(threshold, trigger_type.decode(), max_iterations)
vol.iterations = iterations
vol.volumes = volumes
@ -355,6 +377,8 @@ class VolumeCalculation:
trigger_elem = ET.SubElement(element, "threshold")
trigger_elem.set("type", self.trigger_type)
trigger_elem.set("threshold", str(self.threshold))
if self.max_iterations is not None:
trigger_elem.set("max_iterations", str(self.max_iterations))
return element
@classmethod
@ -398,6 +422,7 @@ class VolumeCalculation:
if trigger_elem is not None:
trigger_type = get_text(trigger_elem, "type")
threshold = float(get_text(trigger_elem, "threshold"))
vol.set_trigger(threshold, trigger_type)
max_iterations = Integral(get_text(trigger_elem, "max_iterations"))
vol.set_trigger(threshold, trigger_type, max_iterations)
return vol

6
src/finalize.cpp
View file

@ -192,6 +192,12 @@ int openmc_finalize()
if (mpi::collision_track_site != MPI_DATATYPE_NULL) {
MPI_Type_free(&mpi::collision_track_site);
}
if (mpi::volume_results != MPI_DATATYPE_NULL) {
MPI_Type_free(&mpi::volume_results);
}
if (mpi::volume_tally != MPI_DATATYPE_NULL) {
MPI_Type_free(&mpi::volume_tally);
}
#endif
openmc_finalize_random_ray();

31
src/initialize.cpp
View file

@ -36,6 +36,7 @@
#include "openmc/thermal.h"
#include "openmc/timer.h"
#include "openmc/vector.h"
#include "openmc/volume_calc.h"
#include "openmc/weight_windows.h"
#ifdef OPENMC_LIBMESH_ENABLED
@ -234,6 +235,36 @@ void initialize_mpi(MPI_Comm intracomm)
MPI_Type_create_struct(
16, blocksc, dispc, typesc, &mpi::collision_track_site);
MPI_Type_commit(&mpi::collision_track_site);
// Volume Calculation data types
VolumeCalculation::CalcResults cr;
MPI_Aint cr_disp[2], cr_d;
MPI_Get_address(&cr, &cr_d);
MPI_Get_address(&cr.n_samples, &cr_disp[0]);
MPI_Get_address(&cr.cost, &cr_disp[1]);
for (int i = 0; i < 2; i++) {
cr_disp[i] -= cr_d;
}
int cr_blocks[] {1, 1};
MPI_Datatype cr_types[] {MPI_UINT64_T, MPI_DOUBLE};
MPI_Type_create_struct(2, cr_blocks, cr_disp, cr_types, &mpi::volume_results);
MPI_Type_commit(&mpi::volume_results);
VolumeCalculation::VolTally vt = VolumeCalculation::VolTally();
MPI_Aint vt_disp[3], vt_d;
MPI_Get_address(&vt, &vt_d);
MPI_Get_address(&vt.score, &vt_disp[0]);
MPI_Get_address(&vt.score_acc, &vt_disp[1]);
MPI_Get_address(&vt.index, &vt_disp[2]);
for (int i = 0; i < 3; i++) {
vt_disp[i] -= vt_d;
}
int vt_blocks[] {1, 2, 1};
MPI_Datatype vt_types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_INT32_T};
MPI_Type_create_struct(3, vt_blocks, vt_disp, vt_types, &mpi::volume_tally);
MPI_Type_commit(&mpi::volume_tally);
}
#endif // OPENMC_MPI

2
src/message_passing.cpp
View file

@ -11,6 +11,8 @@ bool master {true};
MPI_Comm intracomm {MPI_COMM_NULL};
MPI_Datatype source_site {MPI_DATATYPE_NULL};
MPI_Datatype collision_track_site {MPI_DATATYPE_NULL};
MPI_Datatype volume_results {MPI_DATATYPE_NULL};
MPI_Datatype volume_tally {MPI_DATATYPE_NULL};
#endif
extern "C" bool openmc_master()

741
src/volume_calc.cpp
View file

@ -3,15 +3,19 @@
#include "openmc/capi.h"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/distribution_multi.h"
#include "openmc/error.h"
#include "openmc/finalize.h"
#include "openmc/geometry.h"
#include "openmc/hdf5_interface.h"
#include "openmc/initialize.h"
#include "openmc/material.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/openmp_interface.h"
#include "openmc/output.h"
#include "openmc/random_dist.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/timer.h"
@ -32,7 +36,7 @@ namespace openmc {
namespace model {
vector<VolumeCalculation> volume_calcs;
}
} // namespace model
//==============================================================================
// VolumeCalculation implementation
@ -60,6 +64,10 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node)
upper_right_ = get_node_array<double>(node, "upper_right");
n_samples_ = std::stoull(get_node_value(node, "samples"));
// Determine volume of bounding box
const Position d {upper_right_ - lower_left_};
volume_sample_ = d.x * d.y * d.z;
if (check_for_node(node, "threshold")) {
pugi::xml_node threshold_node = node.child("threshold");
@ -70,17 +78,41 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node)
threshold_));
}
// Prior calculation of the trigger volume fractions-related values t'
// deriviated from the given volumes-related values t to prevent excessive
// repeated computations during Monte Carlo execution. The values of t' are
// computed via the sample mean \bar{x} and the adjusted sample variance
// s^2.
std::string tmp = get_node_value(threshold_node, "type");
if (tmp == "variance") {
trigger_type_ = TriggerMetric::variance;
// Condition: s^2 < t' / N
// t' = s^2 = t / (V_b)^2, in sq. volume fraction units
threshold_cnd_ = threshold_ / std::pow(volume_sample_, 2);
;
} else if (tmp == "std_dev") {
trigger_type_ = TriggerMetric::standard_deviation;
// Condition: s^2 < t' / N
// t' = s^2 = (t / V_b)^2, in sq. volume fraction units
threshold_cnd_ = std::pow(threshold_ / volume_sample_, 2);
} else if (tmp == "rel_err") {
trigger_type_ = TriggerMetric::relative_error;
// Condition: s^2 / \bar{x}^2 < t' / N
// t' = s^2 / \bar{x}^2 = t^2, in relative units
threshold_cnd_ = threshold_ * threshold_;
} else {
fatal_error(fmt::format(
"Invalid volume calculation trigger type '{}' provided.", tmp));
}
if (check_for_node(threshold_node, "max_iterations")) {
max_iterations_ =
std::stoi(get_node_value(threshold_node, "max_iterations"));
if (max_iterations_ <= 0) {
fatal_error(fmt::format(
"Invalid error max_iterations {} provided.", max_iterations_));
}
}
}
// Ensure there are no duplicates by copying elements to a set and then
@ -92,7 +124,7 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node)
}
}
vector<VolumeCalculation::Result> VolumeCalculation::execute() const
void VolumeCalculation::execute(CalcResults& master_results) const
{
// Check to make sure domain IDs are valid
for (auto uid : domain_ids_) {
@ -120,12 +152,7 @@ vector<VolumeCalculation::Result> VolumeCalculation::execute() const
}
// Shared data that is collected from all threads
int n = domain_ids_.size();
vector<vector<uint64_t>> master_indices(
n); // List of material indices for each domain
vector<vector<uint64_t>> master_hits(
n); // Number of hits for each material in each domain
int iterations = 0;
const int n_domains = domain_ids_.size();
// Divide work over MPI processes
uint64_t min_samples = n_samples_ / mpi::n_procs;
@ -144,246 +171,202 @@ vector<VolumeCalculation::Result> VolumeCalculation::execute() const
#pragma omp parallel
{
// Variables that are private to each thread
vector<vector<uint64_t>> indices(n);
vector<vector<uint64_t>> hits(n);
// Temporary variables that are private to each thread
CalcResults results(*this);
results.sampling_time.start();
Particle p;
// Sample locations and count hits
// Sample locations and count scores
#pragma omp for
for (size_t i = i_start; i < i_end; i++) {
uint64_t id = iterations * n_samples_ + i;
uint64_t id = master_results.iterations * n_samples_ + i;
uint64_t seed = init_seed(id, STREAM_VOLUME);
p.n_coord() = 1;
Position xi {prn(&seed), prn(&seed), prn(&seed)};
p.r() = lower_left_ + xi * (upper_right_ - lower_left_);
p.u() = {1. / std::sqrt(3.), 1. / std::sqrt(3.), 1. / std::sqrt(3.)};
p.r() = {uniform_distribution(lower_left_.x, upper_right_.x, &seed),
uniform_distribution(lower_left_.y, upper_right_.y, &seed),
uniform_distribution(lower_left_.z, upper_right_.z, &seed)};
const double sqrt3_1 = 1. / std::sqrt(3.);
p.u() = {sqrt3_1, sqrt3_1, sqrt3_1};
// If this location is not in the geometry at all, move on to next block
if (!exhaustive_find_cell(p))
continue;
// TO REVIEWER: THE SWITCH IS TRANSFERED TO score_hit()
if (exhaustive_find_cell(p))
this->score_hit(p, results);
if (domain_type_ == TallyDomain::MATERIAL) {
if (p.material() != MATERIAL_VOID) {
for (int i_domain = 0; i_domain < n; i_domain++) {
if (model::materials[p.material()]->id_ ==
domain_ids_[i_domain]) {
this->check_hit(
p.material(), indices[i_domain], hits[i_domain]);
break;
}
}
}
} else if (domain_type_ == TallyDomain::CELL) {
for (int level = 0; level < p.n_coord(); ++level) {
for (int i_domain = 0; i_domain < n; i_domain++) {
if (model::cells[p.coord(level).cell()]->id_ ==
domain_ids_[i_domain]) {
this->check_hit(
p.material(), indices[i_domain], hits[i_domain]);
break;
}
}
}
} else if (domain_type_ == TallyDomain::UNIVERSE) {
for (int level = 0; level < p.n_coord(); ++level) {
for (int i_domain = 0; i_domain < n; ++i_domain) {
if (model::universes[p.coord(level).universe()]->id_ ==
domain_ids_[i_domain]) {
check_hit(p.material(), indices[i_domain], hits[i_domain]);
break;
}
}
results.n_samples++;
// This passing across all tallies after each sample can be
// inefficient for the case of large number of domains and small
// size of batch, but the batch size is currently assigned to 1
// for keep the input format being unchanged
const double batch_size_1 = 1. / static_cast<double>(1);
for (auto& vt : results.vol_tallies) {
for (auto& vol_tally : vt) {
vol_tally.finalize_batch(batch_size_1);
}
}
}
} // sample/batch loop
// At this point, each thread has its own pair of index/hits lists and we
// now need to reduce them. OpenMP is not nearly smart enough to do this
// on its own, so we have to manually reduce them
for (int i_domain = 0; i_domain < n; ++i_domain) {
reduce_indices_hits(indices[i_domain], hits[i_domain],
master_indices[i_domain], master_hits[i_domain]);
}
results.sampling_time.stop();
results.cost = results.sampling_time.elapsed();
// At this point, each thread has its own volume tallies lists and we
// now need to reduce them. OpenMP is not nearly smart enough to do
// this on its own, so we have to manually reduce them
reduce_results(results, master_results);
} // omp parallel
// Reduce hits onto master process
// bump iteration counter
master_results.iterations++;
// Determine volume of bounding box
Position d {upper_right_ - lower_left_};
double volume_sample = d.x * d.y * d.z;
// TO REVIEWER: THE MPI-RELATED PART IS MOVED TO collect_MPI()
#ifdef OPENMC_MPI
master_results.collect_MPI(); // collect results to master process
#endif
// Process volume estimation results in master process for the trigger state
// determination
bool stop_calc =
mpi::master && (trigger_type_ == TriggerMetric::not_active ||
master_results.iterations == max_iterations_);
// bump iteration counter and get total number
// of samples at this point
iterations++;
uint64_t total_samples = iterations * n_samples_;
// warn user if total sample size is greater than what the uin64_t type can
// represent
if (total_samples == std::numeric_limits<uint64_t>::max()) {
warning("The number of samples has exceeded the type used to track hits. "
"Volume results may be inaccurate.");
if (!stop_calc) {
// Compute current trigger state across totals (0th elements) only
for (const auto& vt : master_results.vol_tallies) {
stop_calc = vt[0].trigger_state(
trigger_type_, threshold_cnd_, master_results.n_samples);
if (!stop_calc)
break;
}
}
// reset
double trigger_val = -INFTY;
// Set size for members of the Result struct
vector<Result> results(n);
for (int i_domain = 0; i_domain < n; ++i_domain) {
// Get reference to result for this domain
auto& result {results[i_domain]};
// Create 2D array to store atoms/uncertainty for each nuclide. Later this
// is compressed into vectors storing only those nuclides that are
// non-zero
auto n_nuc =
settings::run_CE ? data::nuclides.size() : data::mg.nuclides_.size();
auto atoms =
tensor::zeros<double>({static_cast<size_t>(n_nuc), size_t {2}});
#ifdef OPENMC_MPI
if (mpi::master) {
for (int j = 1; j < mpi::n_procs; j++) {
int q;
// retrieve results
MPI_Recv(
&q, 1, MPI_UINT64_T, j, 2 * j, mpi::intracomm, MPI_STATUS_IGNORE);
vector<uint64_t> buffer(2 * q);
MPI_Recv(buffer.data(), 2 * q, MPI_UINT64_T, j, 2 * j + 1,
mpi::intracomm, MPI_STATUS_IGNORE);
for (int k = 0; k < q; ++k) {
bool already_added = false;
for (int m = 0; m < master_indices[i_domain].size(); ++m) {
if (buffer[2 * k] == master_indices[i_domain][m]) {
master_hits[i_domain][m] += buffer[2 * k + 1];
already_added = true;
break;
}
}
if (!already_added) {
master_indices[i_domain].push_back(buffer[2 * k]);
master_hits[i_domain].push_back(buffer[2 * k + 1]);
}
}
}
} else {
int q = master_indices[i_domain].size();
vector<uint64_t> buffer(2 * q);
for (int k = 0; k < q; ++k) {
buffer[2 * k] = master_indices[i_domain][k];
buffer[2 * k + 1] = master_hits[i_domain][k];
}
MPI_Send(&q, 1, MPI_UINT64_T, 0, 2 * mpi::rank, mpi::intracomm);
MPI_Send(buffer.data(), 2 * q, MPI_UINT64_T, 0, 2 * mpi::rank + 1,
mpi::intracomm);
}
// Send the state of calculation continuation just obtained in master
// process to all processes
MPI_Bcast(&stop_calc, 1, MPI_CXX_BOOL, 0, mpi::intracomm);
#endif
if (mpi::master) {
size_t total_hits = 0;
for (int j = 0; j < master_indices[i_domain].size(); ++j) {
total_hits += master_hits[i_domain][j];
double f =
static_cast<double>(master_hits[i_domain][j]) / total_samples;
double var_f = f * (1.0 - f) / total_samples;
if (!stop_calc)
continue; // while loop
// No trigger is applied or the trigger condition is satisfied, we're
// done
int i_material = master_indices[i_domain][j];
if (i_material == MATERIAL_VOID)
if (mpi::master) {
// Normalize all results on the bounding primitive volume and compute
// stddev
for (auto& vt : master_results.vol_tallies) {
for (auto& vol_tally : vt) {
vol_tally.score_acc = get_tally_results(
master_results.n_samples, volume_sample_, vol_tally);
}
}
// Compute nuclides
for (int i_domain = 0; i_domain < n_domains; ++i_domain) {
// Create 2D array to store atoms/uncertainty for each nuclide. Later
// this is compressed into vectors storing only those nuclides that are
// non-zero
auto n_nuc =
settings::run_CE ? data::nuclides.size() : data::mg.nuclides_.size();
auto atoms =
tensor::zeros<double>({static_cast<size_t>(n_nuc), size_t {2}});
for (int j = 0; j < master_results.vol_tallies[i_domain].size(); ++j) {
const int i_material = master_results.vol_tallies[i_domain][j].index;
if (i_material == MATERIAL_VOID || i_material == _INDEX_TOTAL)
continue;
const auto& mat = model::materials[i_material];
for (int k = 0; k < mat->nuclide_.size(); ++k) {
// Accumulate nuclide density
int i_nuclide = mat->nuclide_[k];
atoms(i_nuclide, 0) += mat->atom_density_[k] * f;
atoms(i_nuclide, 1) += std::pow(mat->atom_density_[k], 2) * var_f;
}
}
// Determine volume
result.volume[0] =
static_cast<double>(total_hits) / total_samples * volume_sample;
result.volume[1] =
std::sqrt(result.volume[0] * (volume_sample - result.volume[0]) /
total_samples);
result.iterations = iterations;
// update threshold value if needed
if (trigger_type_ != TriggerMetric::not_active) {
double val = 0.0;
switch (trigger_type_) {
case TriggerMetric::standard_deviation:
val = result.volume[1];
break;
case TriggerMetric::relative_error:
val = result.volume[0] == 0.0 ? INFTY
: result.volume[1] / result.volume[0];
break;
case TriggerMetric::variance:
val = result.volume[1] * result.volume[1];
break;
default:
break;
}
// update max if entry is valid
if (val > 0.0) {
trigger_val = std::max(trigger_val, val);
auto& volume = master_results.vol_tallies[i_domain][j].score_acc;
// Collect calculated nuclide amounts N [atoms] and stddev as
// N = V [cm^3] * \rho [atoms/b-cm] * 1.e24 [b-cm/cm^3]
atoms(mat->nuclide_[k], 0) +=
volume[0] * mat->atom_density(k) * 1.0e24;
atoms(mat->nuclide_[k], 1) +=
volume[1] * mat->atom_density(k) * 1.0e24;
}
}
// Get reference to result for this domain
auto& result {master_results.nuc_results[i_domain]};
// Convert full arrays to vectors
for (int j = 0; j < n_nuc; ++j) {
// Determine total number of atoms. At this point, we have values in
// atoms/b-cm. To get to atoms we multiply by 10^24 V.
double mean = 1.0e24 * volume_sample * atoms(j, 0);
double stdev = 1.0e24 * volume_sample * std::sqrt(atoms(j, 1));
// Convert full arrays to vectors
if (mean > 0.0) {
if (atoms(j, 0) > 0.0) {
result.nuclides.push_back(j);
result.atoms.push_back(mean);
result.uncertainty.push_back(stdev);
result.atoms.push_back(atoms(j, 0));
result.uncertainty.push_back(atoms(j, 1));
}
}
}
} // end domain loop
// if no trigger is applied, we're done
if (trigger_type_ == TriggerMetric::not_active) {
return results;
} // end domains loop
}
#ifdef OPENMC_MPI
// update maximum error value on all processes
MPI_Bcast(&trigger_val, 1, MPI_DOUBLE, 0, mpi::intracomm);
#endif
// return results of the calculation
if (trigger_val < threshold_) {
return results;
}
#ifdef OPENMC_MPI
// if iterating in an MPI run, need to zero indices and hits so they aren't
// counted twice
if (!mpi::master) {
for (auto& v : master_indices) {
std::fill(v.begin(), v.end(), 0);
}
for (auto& v : master_hits) {
std::fill(v.begin(), v.end(), 0);
}
}
#endif
return;
} // end while
}
void VolumeCalculation::show_vol_stat(
const std::string label, const std::string units, const double value) const
{
fmt::print("{0:<20} = {2:10.4e} {1:<}\n", label, units, value);
}
void VolumeCalculation::show_volume(const std::string domain_type,
const int domain_id, const std::string region_name, const double mean,
const double stddev) const
{
fmt::print(" {0:>9}{1:>6}: {2:10.4e} +/- {3:10.4e} cm^3", domain_type,
domain_id, mean, stddev);
if (!region_name.empty()) {
fmt::print(" //{:<}", region_name);
}
fmt::print("\n");
}
// TO REVIEWER: I/O BLOCK BEGINS, I WOULD PREFFER TO PUT IT AT THE END OF FILE
// BUT KEEPT THE FILE STRUCTURE UNCHANGED
void VolumeCalculation::show_results(const CalcResults& results) const
{
// Show tracing statistics
write_message(5, " ");
show_vol_stat(
"Total sample size", "", static_cast<double>(results.n_samples));
show_vol_stat("Running cost", "thread-sec", results.cost);
show_vol_stat("Cost of hitting", "thread-sec/hit",
static_cast<double>(results.cost) / static_cast<double>(results.n_samples));
write_message(5, " ");
std::string domain_type;
if (domain_type_ == TallyDomain::CELL) {
domain_type = "Cell";
} else if (domain_type_ == TallyDomain::MATERIAL) {
domain_type = "Material";
} else {
domain_type = "Universe";
}
// Display domain volumes
for (int j = 0; j < domain_ids_.size(); j++) {
std::string region_name {""};
if (domain_type_ == TallyDomain::CELL) {
int cell_idx = model::cell_map[domain_ids_[j]];
region_name = model::cells[cell_idx]->name();
} else if (domain_type_ == TallyDomain::MATERIAL) {
int mat_idx = model::material_map[domain_ids_[j]];
region_name = model::materials[mat_idx]->name();
}
if (region_name.size())
region_name.insert(0, " "); // prepend space for formatting
show_volume(domain_type, domain_ids_[j], region_name,
results.vol_tallies[j][0].score_acc[0],
results.vol_tallies[j][0].score_acc[1]);
}
write_message(4, " "); // Blank line afer results printed
}
void VolumeCalculation::to_hdf5(
const std::string& filename, const vector<Result>& results) const
const std::string& filename, const CalcResults& results) const
{
// Create HDF5 file
hid_t file_id = file_open(filename, 'w');
@ -405,7 +388,7 @@ void VolumeCalculation::to_hdf5(
write_attribute(file_id, "upper_right", upper_right_);
// Write trigger info
if (trigger_type_ != TriggerMetric::not_active) {
write_attribute(file_id, "iterations", results[0].iterations);
write_attribute(file_id, "iterations", results.iterations);
write_attribute(file_id, "threshold", threshold_);
std::string trigger_str;
switch (trigger_type_) {
@ -422,6 +405,10 @@ void VolumeCalculation::to_hdf5(
break;
}
write_attribute(file_id, "trigger_type", trigger_str);
// check max_iterations on default value
if (max_iterations_ !=
std::numeric_limits<decltype(max_iterations_)>::max())
write_attribute(file_id, "max_iterations", max_iterations_);
} else {
write_attribute(file_id, "iterations", 1);
}
@ -439,8 +426,8 @@ void VolumeCalculation::to_hdf5(
create_group(file_id, fmt::format("domain_{}", domain_ids_[i]));
// Write volume for domain
const auto& result {results[i]};
write_dataset(group_id, "volume", result.volume);
const auto& result {results.nuc_results[i]};
write_dataset(group_id, "volume", results.vol_tallies[i][0].score_acc);
// Create array of nuclide names from the vector
auto n_nuc = result.nuclides.size();
@ -468,24 +455,296 @@ void VolumeCalculation::to_hdf5(
file_close(file_id);
}
void VolumeCalculation::check_hit(
int i_material, vector<uint64_t>& indices, vector<uint64_t>& hits) const
// TO REVIEWER: THIS IS THE SAME check_hit()
void VolumeCalculation::check_hit(const int32_t i_material,
const double contrib, vector<VolTally>& vol_tallies) const
{
// Contribute to entire domain result tally
vol_tallies[0].score += contrib;
// Check if this material was previously hit and if so, increment count
bool already_hit = false;
for (int j = 0; j < indices.size(); j++) {
if (indices[j] == i_material) {
hits[j]++;
already_hit = true;
// Check if this material was previously hit and if so, contribute score
for (int j = 1; j < vol_tallies.size(); j++) {
if (vol_tallies[j].index == i_material) {
vol_tallies[j].score += contrib;
return;
}
}
// If the material was not previously hit, append an entry to the material
// The material was not previously hit, append an entry to the material
// indices and hits lists
if (!already_hit) {
indices.push_back(i_material);
hits.push_back(1);
vol_tallies.push_back(VolTally(i_material, contrib));
}
// TO REVIEWER: THIS IS AN EQUIVALENT OF THE reduce_indicies_hits() TEMPLATE
// FROM volume_calc.h
void VolumeCalculation::reduce_results(
const CalcResults& local_results, CalcResults& results) const
{
auto n_threads = num_threads();
// Collect scalar variables
#pragma omp for ordered schedule(static, 1)
for (int i = 0; i < n_threads; ++i) {
#pragma omp ordered
{
results.n_samples += local_results.n_samples;
results.cost += local_results.cost;
}
}
// Collect vectored domain-wise results
for (int i_domain = 0; i_domain < domain_ids_.size(); ++i_domain) {
const vector<VolTally>& local_vol_tall =
local_results.vol_tallies[i_domain];
vector<VolTally>& vol_tall = results.vol_tallies[i_domain];
#pragma omp for ordered schedule(static, 1)
for (int i = 0; i < n_threads; ++i) {
#pragma omp ordered
for (int j = 0; j < local_vol_tall.size(); ++j) {
// Check if this material has been added to the master list and if
// so, accumulate scores
const auto ind {local_vol_tall[j].index};
const auto it = std::find_if(vol_tall.begin(), vol_tall.end(),
[ind](const VolTally& vt) { return vt.index == ind; });
if (it == vol_tall.end()) {
vol_tall.push_back(local_vol_tall[j]);
} else {
vol_tall[it - vol_tall.begin()].append_tally(local_vol_tall[j]);
}
}
}
}
}
#ifdef OPENMC_MPI
// TO REVIEWER: THIS IS AN EQUIVALENT OF THE RELATED MPI-INTERCHANGE BLOCK OF
// execute()
void VolumeCalculation::CalcResults::collect_MPI()
{
vector<int> domain_sizes(vol_tallies.size());
if (!mpi::master) {
// n_domain + 2 MPI messages will be send in total to node mpi::master
// To determination of an unique tag for each MPI message (as below)
int mpi_offset = mpi::rank * (vol_tallies.size() + 2) + 2;
// Pass root data of the struct
MPI_Send(
(void*)this, 1, mpi::volume_results, 0, mpi_offset - 2, mpi::intracomm);
// Pass sizes of domain-wise data
for (int i_domain = 0; i_domain < vol_tallies.size(); i_domain++) {
domain_sizes[i_domain] = vol_tallies[i_domain].size();
}
MPI_Send(domain_sizes.data(), domain_sizes.size(), MPI_INT, 0,
mpi_offset - 1, mpi::intracomm);
// Pass domain-wise data of struct
for (int i_domain = 0; i_domain < vol_tallies.size(); i_domain++) {
MPI_Send(vol_tallies[i_domain].data(), domain_sizes[i_domain],
mpi::volume_tally, 0, mpi_offset + i_domain, mpi::intracomm);
}
this->reset(); // Delete passed to main process data
} else {
// n_domain + 2 MPI messages will be recieved in total from node mpi::master
for (int i_proc = 1; i_proc < mpi::n_procs; i_proc++) {
CalcResults res_buff(*this); // temporary storage for recived data
// To determination of an unique tag for each MPI message (as above)
int mpi_offset = i_proc * (vol_tallies.size() + 2) + 2;
// Pass root data of struct
MPI_Recv(&res_buff, 1, mpi::volume_results, i_proc, mpi_offset - 2,
mpi::intracomm, MPI_STATUS_IGNORE);
// Pass sizes of domain-wise data
MPI_Recv(domain_sizes.data(), domain_sizes.size(), MPI_INT, i_proc,
mpi_offset - 1, mpi::intracomm, MPI_STATUS_IGNORE);
// Pass domain-wise data of struct
for (int i_domain = 0; i_domain < vol_tallies.size(); i_domain++) {
res_buff.vol_tallies[i_domain].resize(domain_sizes[i_domain]);
MPI_Recv(res_buff.vol_tallies[i_domain].data(), domain_sizes[i_domain],
mpi::volume_tally, i_proc, mpi_offset + i_domain, mpi::intracomm,
MPI_STATUS_IGNORE);
}
this->append(res_buff);
}
}
}
#endif
// TO REVIEWER: STATISTICS AND TALLY MANIPULATION BLOCK BEGIN
VolumeCalculation::CalcResults::CalcResults(const VolumeCalculation& vol_calc)
{
n_samples = 0;
iterations = 0;
cost = 0.;
for (int i = 0; i < vol_calc.domain_ids_.size(); i++) {
vector<VolTally> vt_vect; // Tally group for a domain
vt_vect.push_back(
VolTally(_INDEX_TOTAL)); // Zero-element tally for entire domain totals
vol_tallies.push_back(vt_vect);
nuc_results.push_back(NuclResult());
}
}
void VolumeCalculation::CalcResults::reset()
{
n_samples = 0;
cost = 0.;
for (auto& vt : vol_tallies) {
std::fill(vt.begin(), vt.end(), VolTally());
}
}
// TO REVIEWER: THIS IS A PART OF THE RELATED TO MPI-INTERCHANGE
// BLOCK OF execute()
void VolumeCalculation::CalcResults::append(const CalcResults& other)
{
n_samples += other.n_samples;
cost += other.cost;
// The domain-wise vectors this.vol_tallies and other.vol_tallies are
// conformed each to other by definition
for (auto id = 0; id < vol_tallies.size(); id++) {
// Merging current domain vector from other.vol_tallies into this via
// pair-wise comparisons
for (const auto& vt_other : other.vol_tallies[id]) {
bool already_appended = false;
for (auto& vt : vol_tallies[id]) {
if (vt.index == vt_other.index) {
vt.append_tally(vt_other);
already_appended = true;
break;
}
}
if (!already_appended)
vol_tallies[id].push_back(vt_other);
}
}
}
inline void VolumeCalculation::VolTally::finalize_batch(
const double batch_size_1)
{
if (score != 0.) {
score *= batch_size_1;
score_acc[0] += score;
score_acc[1] += score * score;
score = 0.;
}
}
inline void VolumeCalculation::VolTally::assign_tally(const VolTally& vol_tally)
{
score = vol_tally.score;
score_acc = vol_tally.score_acc;
index = vol_tally.index;
}
inline void VolumeCalculation::VolTally::append_tally(const VolTally& vol_tally)
{
score += vol_tally.score;
score_acc[0] += vol_tally.score_acc[0];
score_acc[1] += vol_tally.score_acc[1];
}
// TO REVIEWER: IF THIS APPROACH WILL BE FOUND OBFUSCATED, IT CAN BE ADJUSTED TO
// LITERAL FORMULAE TRANSLATION WITH REPEATING MULTIPLY CALCULATIONS OF
// THRESHOLD SQUARE ROOTS AND DIVISIONS FOR EACH VOLUME DOMAIN IN THE END OF
// EACH BATCH
inline bool VolumeCalculation::VolTally::trigger_state(
const TriggerMetric trigger_type, const double threshold,
const size_t& n_samples) const
{
// For sample contribution to volume fraction limited by 1, the maximal
// allowed n_samples value is around 1.e102, but this is still much larger
// than the size_t limit equal to ~1.8e19
const double ns1 = static_cast<double>(n_samples - 1);
const double ns = static_cast<double>(n_samples);
const double mean_xi_sq = score_acc[0] * score_acc[0];
// Adjusted sample variance: s^2 = (\bar{x^2} - \bar{x}^2) / (N-1)
// \bar{x}^2 = mean_xi_sq / N^2, \bar{\x^2} = score_acc[1] / N, N = ns
switch (trigger_type) {
case TriggerMetric::variance:
case TriggerMetric::standard_deviation:
// Condition: s^2 / N < t'
// Equivalent implementation:
// N^2 * (\bar{x^2} - \bar{x}^2) < t' * (N-1) * N^2
return score_acc[1] * ns - mean_xi_sq < threshold * ns1 * ns * ns;
case TriggerMetric::relative_error:
// Condition: (s^2 / \mu^2) / N < t'
// Equivalent implementation:
// N^2 * (\bar{x^2} - \bar{x}^2) < t' * (N-1) * (N * \bar{x})^2
return score_acc[1] * ns - mean_xi_sq < threshold * ns1 * mean_xi_sq;
default:
return true;
}
}
array<double, 2> VolumeCalculation::get_tally_results(const size_t& n_samples,
const double coeff_norm, const VolTally& vol_tally) const
{
array<double, 2> volume;
const double ns_1 = 1. / static_cast<double>(n_samples);
volume[0] = vol_tally.score_acc[0] * ns_1;
volume[1] = vol_tally.score_acc[1] * ns_1;
volume[1] = std::sqrt(
(volume[1] - volume[0] * volume[0]) / static_cast<double>(n_samples - 1));
volume[0] *= coeff_norm;
volume[1] *= coeff_norm;
return volume;
}
// TO REVIEWER: THIS IS A FULL EQUIVALENT OF THE SWITCH FROM THE HISTORY LOOP
void VolumeCalculation::score_hit(const Particle& p, CalcResults& results) const
{
const auto n_domains = domain_ids_.size();
const auto id_mat = p.material();
const auto score = 1.; // Floating-point score value
switch (domain_type_) {
case TallyDomain::MATERIAL:
if (id_mat != MATERIAL_VOID) {
for (auto i_domain = 0; i_domain < n_domains; i_domain++) {
if (model::materials[id_mat]->id_ == domain_ids_[i_domain]) {
this->check_hit(id_mat, score, results.vol_tallies[i_domain]);
break;
}
}
}
break;
case TallyDomain::CELL:
for (auto level = 0; level < p.n_coord_last(); ++level) {
for (auto i_domain = 0; i_domain < n_domains; i_domain++) {
if (model::cells[p.coord(level).cell()]->id_ == domain_ids_[i_domain]) {
this->check_hit(id_mat, score, results.vol_tallies[i_domain]);
break;
}
}
}
break;
case TallyDomain::UNIVERSE:
for (auto level = 0; level < p.n_coord_last(); ++level) {
for (auto i_domain = 0; i_domain < n_domains; ++i_domain) {
if (model::universes[model::cells[p.coord(level).cell()]->universe_]
->id_ == domain_ids_[i_domain]) {
this->check_hit(id_mat, score, results.vol_tallies[i_domain]);
break;
}
}
}
}
}
@ -516,43 +775,19 @@ int openmc_calculate_volumes()
// Run volume calculation
const auto& vol_calc {model::volume_calcs[i]};
std::vector<VolumeCalculation::Result> results;
VolumeCalculation::CalcResults results(vol_calc);
try {
results = vol_calc.execute();
vol_calc.execute(results);
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
if (mpi::master) {
std::string domain_type;
if (vol_calc.domain_type_ == VolumeCalculation::TallyDomain::CELL) {
domain_type = " Cell ";
} else if (vol_calc.domain_type_ ==
VolumeCalculation::TallyDomain::MATERIAL) {
domain_type = " Material ";
} else {
domain_type = " Universe ";
}
// Display domain volumes
for (int j = 0; j < vol_calc.domain_ids_.size(); j++) {
std::string region_name {""};
if (vol_calc.domain_type_ == VolumeCalculation::TallyDomain::CELL) {
int cell_idx = model::cell_map[vol_calc.domain_ids_[j]];
region_name = model::cells[cell_idx]->name();
} else if (vol_calc.domain_type_ ==
VolumeCalculation::TallyDomain::MATERIAL) {
int mat_idx = model::material_map[vol_calc.domain_ids_[j]];
region_name = model::materials[mat_idx]->name();
}
if (region_name.size())
region_name.insert(0, " "); // prepend space for formatting
write_message(4, "{}{}{}: {} +/- {} cm^3", domain_type,
vol_calc.domain_ids_[j], region_name, results[j].volume[0],
results[j].volume[1]);
}
// Output volume calculation results and statistics
vol_calc.show_results(results);
// Write volumes to HDF5 file
std::string filename =
@ -563,7 +798,7 @@ int openmc_calculate_volumes()
// Show elapsed time
time_volume.stop();
write_message(6, "Elapsed time: {} s", time_volume.elapsed());
write_message(6, "Elapsed time: {} sec", time_volume.elapsed());
return 0;
}

10
tests/regression_tests/volume_calc/inputs_true.dat
View file

@ -17,7 +17,7 @@
<geometry>
<cell id="1" material="2" region="-1 -3 5" universe="0"/>
<cell id="2" material="1" region="-2 3" universe="0"/>
<cell id="3" material="1" region="-4 -3" universe="0"/>
<cell id="3" material="1" region="-4 -5" universe="0"/>
<surface id="1" type="z-cylinder" boundary="vacuum" coeffs="0.0 0.0 1.0"/>
<surface id="2" type="sphere" boundary="vacuum" coeffs="0.0 0.0 5.0 1.0"/>
<surface id="3" type="z-plane" coeffs="5.0"/>
@ -71,5 +71,13 @@
<upper_right>1.0 1.0 6.0</upper_right>
<threshold type="variance" threshold="0.05"/>
</volume_calc>
<volume_calc>
<domain_type>material</domain_type>
<domain_ids>1 2</domain_ids>
<samples>100</samples>
<lower_left>-1.0 -1.0 -6.0</lower_left>
<upper_right>1.0 1.0 6.0</upper_right>
<threshold type="rel_err" threshold="1e-10" max_iterations="9"/>
</volume_calc>
</settings>
</model>

10
tests/regression_tests/volume_calc/inputs_true_mg.dat
View file

@ -17,7 +17,7 @@
<geometry>
<cell id="1" material="2" region="-1 -3 5" universe="0"/>
<cell id="2" material="1" region="-2 3" universe="0"/>
<cell id="3" material="1" region="-4 -3" universe="0"/>
<cell id="3" material="1" region="-4 -5" universe="0"/>
<surface id="1" type="z-cylinder" boundary="vacuum" coeffs="0.0 0.0 1.0"/>
<surface id="2" type="sphere" boundary="vacuum" coeffs="0.0 0.0 5.0 1.0"/>
<surface id="3" type="z-plane" coeffs="5.0"/>
@ -72,5 +72,13 @@
<upper_right>1.0 1.0 6.0</upper_right>
<threshold type="variance" threshold="0.05"/>
</volume_calc>
<volume_calc>
<domain_type>material</domain_type>
<domain_ids>1 2</domain_ids>
<samples>100</samples>
<lower_left>-1.0 -1.0 -6.0</lower_left>
<upper_right>1.0 1.0 6.0</upper_right>
<threshold type="rel_err" threshold="1e-10" max_iterations="9"/>
</volume_calc>
</settings>
</model>

31
tests/regression_tests/volume_calc/results_true.dat
View file

@ -1,7 +1,8 @@
Volume calculation 0
Trigger Type: None
Trigger threshold: None
Iterations: 1
Iterations limit: None
Iterations completed: 1
Domain 1: 31.36+/-0.07 cm^3
Domain 2: 2.107+/-0.031 cm^3
Domain 3: 2.164+/-0.031 cm^3
@ -17,7 +18,8 @@ Domain 3: 2.164+/-0.031 cm^3
Volume calculation 1
Trigger Type: None
Trigger threshold: None
Iterations: 1
Iterations limit: None
Iterations completed: 1
Domain 1: 4.27+/-0.04 cm^3
Domain 2: 31.36+/-0.07 cm^3
Material Nuclide Atoms
@ -29,7 +31,8 @@ Domain 2: 31.36+/-0.07 cm^3
Volume calculation 2
Trigger Type: None
Trigger threshold: None
Iterations: 1
Iterations limit: None
Iterations completed: 1
Domain 0: 35.63+/-0.07 cm^3
Universe Nuclide Atoms
0 0 H1 (2.856+/-0.029)e+23
@ -40,7 +43,8 @@ Domain 0: 35.63+/-0.07 cm^3
Volume calculation 3
Trigger Type: std_dev
Trigger threshold: 0.1
Iterations: 523
Iterations limit: None
Iterations completed: 523
Domain 1: 31.31+/-0.10 cm^3
Domain 2: 2.10+/-0.04 cm^3
Domain 3: 2.19+/-0.04 cm^3
@ -56,7 +60,8 @@ Domain 3: 2.19+/-0.04 cm^3
Volume calculation 4
Trigger Type: rel_err
Trigger threshold: 0.1
Iterations: 9
Iterations limit: None
Iterations completed: 9
Domain 1: 5.3+/-0.5 cm^3
Domain 2: 29.9+/-0.8 cm^3
Material Nuclide Atoms
@ -68,7 +73,8 @@ Domain 2: 29.9+/-0.8 cm^3
Volume calculation 5
Trigger Type: variance
Trigger threshold: 0.05
Iterations: 105
Iterations limit: None
Iterations completed: 105
Domain 1: 31.16+/-0.22 cm^3
Domain 2: 1.91+/-0.09 cm^3
Domain 3: 2.32+/-0.10 cm^3
@ -81,3 +87,16 @@ Domain 3: 2.32+/-0.10 cm^3
5 3 H1 (1.55+/-0.07)e+23
6 3 O16 (7.76+/-0.34)e+22
7 3 B10 (7.76+/-0.34)e+18
Volume calculation 6
Trigger Type: rel_err
Trigger threshold: 1e-10
Iterations limit: 9
Iterations completed: 9
Domain 1: 5.3+/-0.5 cm^3
Domain 2: 29.9+/-0.8 cm^3
Material Nuclide Atoms
0 1 H1 (3.53+/-0.33)e+23
1 1 O16 (1.77+/-0.17)e+23
2 1 B10 (1.77+/-0.17)e+19
3 2 U235 (3.31+/-0.09)e+23
4 2 Mo99 (3.31+/-0.09)e+22

35
tests/regression_tests/volume_calc/test.py
View file

@ -16,7 +16,9 @@ class VolumeTest(PyAPITestHarness):
self.exp_std_dev = 1e-01
self.exp_rel_err = 1e-01
self.tiny_rel_err = 1e-10
self.exp_variance = 5e-02
self.max_iterations = 9
self.is_ce = is_ce
if not is_ce:
self.inputs_true = 'inputs_true_mg.dat'
@ -49,7 +51,7 @@ class VolumeTest(PyAPITestHarness):
# Define geometry
inside_cyl = openmc.Cell(1, fill=fuel, region=-cyl & -top_plane & +bottom_plane)
top_hemisphere = openmc.Cell(2, fill=water, region=-top_sphere & +top_plane)
bottom_hemisphere = openmc.Cell(3, fill=water, region=-bottom_sphere & -top_plane)
bottom_hemisphere = openmc.Cell(3, fill=water, region=-bottom_sphere & -bottom_plane)
root = openmc.Universe(0, cells=(inside_cyl, top_hemisphere, bottom_hemisphere))
self._model.geometry = openmc.Geometry(root)
@ -62,7 +64,8 @@ class VolumeTest(PyAPITestHarness):
openmc.VolumeCalculation([root], 100000, ll, ur),
openmc.VolumeCalculation(list(root.cells.values()), 100),
openmc.VolumeCalculation([water, fuel], 100, ll, ur),
openmc.VolumeCalculation(list(root.cells.values()), 100)
openmc.VolumeCalculation(list(root.cells.values()), 100),
openmc.VolumeCalculation([water, fuel], 100, ll, ur)
]
vol_calcs[3].set_trigger(self.exp_std_dev, 'std_dev')
@ -71,6 +74,9 @@ class VolumeTest(PyAPITestHarness):
vol_calcs[5].set_trigger(self.exp_variance, 'variance')
vol_calcs[6].set_trigger(self.tiny_rel_err, 'rel_err',
self.max_iterations)
# Define settings
settings = openmc.Settings()
settings.run_mode = 'volume'
@ -130,30 +136,45 @@ class VolumeTest(PyAPITestHarness):
outstr += 'Trigger Type: {}\n'.format(volume_calc.trigger_type)
outstr += 'Trigger threshold: {}\n'.format(volume_calc.threshold)
outstr += 'Iterations: {}\n'.format(volume_calc.iterations)
outstr += 'Iterations limit: {}\n'.format(volume_calc.max_iterations)
outstr += 'Iterations completed: {}\n'.format(volume_calc.iterations)
# make sure the volume calculation results contains true
# values of trigger and iteration limit
if i == 3:
assert volume_calc.trigger_type == 'std_dev'
assert volume_calc.threshold == self.exp_std_dev
assert volume_calc.max_iterations is None
elif i == 4:
assert volume_calc.trigger_type == 'rel_err'
assert volume_calc.threshold == self.exp_rel_err
assert volume_calc.max_iterations is None
elif i == 5:
assert volume_calc.trigger_type == 'variance'
assert volume_calc.threshold == self.exp_variance
assert volume_calc.max_iterations is None
elif i == 6:
assert volume_calc.trigger_type == 'rel_err'
assert volume_calc.threshold == self.tiny_rel_err
assert volume_calc.max_iterations == self.max_iterations
else:
assert volume_calc.trigger_type is None
assert volume_calc.threshold is None
assert volume_calc.max_iterations is None
assert volume_calc.iterations == 1
# if a trigger is applied, make sure the calculation satisfies the trigger
# if a trigger is applied, make sure the calculation satisfies
# the trigger and iteration limit
for vol in volume_calc.volumes.values():
if volume_calc.trigger_type == 'std_dev':
assert vol.std_dev <= self.exp_std_dev
assert (vol.std_dev <= volume_calc.threshold or
volume_calc.iterations == volume_calc.max_iterations)
if volume_calc.trigger_type == 'rel_err':
assert vol.std_dev/vol.nominal_value <= self.exp_rel_err
assert (vol.std_dev/vol.nominal_value <= volume_calc.threshold or
volume_calc.iterations == volume_calc.max_iterations)
if volume_calc.trigger_type == 'variance':
assert vol.std_dev * vol.std_dev <= self.exp_variance
assert (vol.std_dev * vol.std_dev <= volume_calc.threshold or
volume_calc.iterations == volume_calc.max_iterations)
# Write cell volumes and total # of atoms for each nuclide
for uid, volume in sorted(volume_calc.volumes.items()):