crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1374 from pshriwise/volume_trigger

Browse source 
Volume Calculation Trigger
This commit is contained in:
Paul Romano 2019-10-15 15:58:20 +07:00 • committed by GitHub
commit 6799da13f1
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
8 changed files with 452 additions and 150 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -23,6 +23,8 @@ The current version of the volume file format is 1.0.
bounding box
- **upper_right** (*double[3]*) -- Upper-right coordinates of
bounding box
- **threshold** (*double*) -- Threshold used for volume uncertainty
- **trigger_type** (*char[]*) -- Trigger type used for volume uncertainty
**/domain_<id>/**

10
docs/source/usersguide/volume.rst
View file

@ -37,6 +37,16 @@ arguments are not necessary. For example,
Of course, the volumes that you *need* this capability for are often the ones
with complex definitions.
A threshold can be applied for the calculation's variance, standard deviation,
or relative error of volume estimates using :meth:`openmc.VolumeCalculation.set_trigger`::
vol_calc.set_trigger(1e-05, 'std_dev')
If a threshold is provided, calculations will be performed iteratively using the
number of samples specified on the calculation until all volume uncertainties are below
the threshold value. If no threshold is provided, the calculation will run the number of
samples specified once and return the result.
Once you have one or more :class:`openmc.VolumeCalculation` objects created, you
can then assign then to :attr:`Settings.volume_calculations`::

8
include/openmc/volume_calc.h
View file

@ -2,12 +2,14 @@
#define OPENMC_VOLUME_CALC_H
#include "openmc/position.h"
#include "openmc/tallies/trigger.h"
#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"
#include <string>
#include <vector>
#include <gsl/gsl>
namespace openmc {
@ -16,6 +18,7 @@ namespace openmc {
//==============================================================================
class VolumeCalculation {
public:
// Aliases, types
struct Result {
@ -23,6 +26,7 @@ public:
std::vector<int> nuclides; //!< Index of nuclides
std::vector<double> atoms; //!< Number of atoms for each nuclide
std::vector<double> uncertainty; //!< Uncertainty on number of atoms
int iterations; //!< Number of iterations needed to obtain the results
}; // Results for a single domain
// Constructors
@ -44,7 +48,9 @@ public:
// Data members
int domain_type_; //!< Type of domain (cell, material, etc.)
int n_samples_; //!< Number of samples to use
size_t n_samples_; //!< Number of samples to use
double threshold_ {-1.0}; //!< Error threshold for domain volumes
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
std::vector<int> domain_ids_; //!< IDs of domains to find volumes of

72
openmc/volume.py
View file

@ -45,6 +45,8 @@ class VolumeCalculation(object):
Lower-left coordinates of bounding box used to sample points
upper_right : Iterable of float
Upper-right coordinates of bounding box used to sample points
threshold : float
Threshold for the maximum standard deviation of volume in the calculation
atoms : dict
Dictionary mapping unique IDs of domains to a mapping of nuclides to
total number of atoms for each nuclide present in the domain. For
@ -54,12 +56,20 @@ class VolumeCalculation(object):
in each domain specified.
volumes : dict
Dictionary mapping unique IDs of domains to estimated volumes in cm^3.
threshold : float
Threshold for the maxmimum standard deviation of volumes.
trigger_type : {'variance', 'std_dev', 'rel_err'}
Value type used to halt volume calculation
iterations : int
Number of iterations over samples (for calculations with a trigger).
"""
def __init__(self, domains, samples, lower_left=None,
upper_right=None):
def __init__(self, domains, samples, lower_left=None, upper_right=None):
self._atoms = {}
self._volumes = {}
self._threshold = None
self._trigger_type = None
self._iterations = None
cv.check_type('domains', domains, Iterable,
(openmc.Cell, openmc.Material, openmc.Universe))
@ -123,6 +133,18 @@ class VolumeCalculation(object):
def upper_right(self):
return self._upper_right
@property
def threshold(self):
return self._threshold
@property
def trigger_type(self):
return self._trigger_type
@property
def iterations(self):
return self._iterations
@property
def domain_type(self):
return self._domain_type
@ -170,6 +192,26 @@ class VolumeCalculation(object):
cv.check_length(name, upper_right, 3)
self._upper_right = upper_right
@threshold.setter
def threshold(self, threshold):
name = 'volume std. dev. threshold'
cv.check_type(name, threshold, Real)
cv.check_greater_than(name, threshold, 0.0)
self._threshold = threshold
@trigger_type.setter
def trigger_type(self, trigger_type):
cv.check_value('tally trigger type', trigger_type,
('variance', 'std_dev', 'rel_err'))
self._trigger_type = trigger_type
@iterations.setter
def iterations(self, iterations):
name = 'volume calculation iterations'
cv.check_type(name, iterations, Integral)
cv.check_greater_than(name, iterations, 0)
self._iterations = iterations
@volumes.setter
def volumes(self, volumes):
cv.check_type('volumes', volumes, Mapping)
@ -180,6 +222,19 @@ class VolumeCalculation(object):
cv.check_type('atoms', atoms, Mapping)
self._atoms = atoms
def set_trigger(self, threshold, trigger_type):
"""Set a trigger on the voulme calculation
Parameters
----------
threshold : float
Threshold for the maxmimum standard deviation of volumes
trigger_type : {'variance', 'std_dev', 'rel_err'}
Value type used to halt volume calculation
"""
self.trigger_type = trigger_type
self.threshold = threshold
@classmethod
def from_hdf5(cls, filename):
"""Load stochastic volume calculation results from HDF5 file.
@ -203,6 +258,10 @@ class VolumeCalculation(object):
lower_left = f.attrs['lower_left']
upper_right = f.attrs['upper_right']
threshold = f.attrs.get('threshold')
trigger_type = f.attrs.get('trigger_type')
iterations = f.attrs.get('iterations', 1)
volumes = {}
atoms = {}
ids = []
@ -233,6 +292,11 @@ class VolumeCalculation(object):
# Instantiate the class and assign results
vol = cls(domains, samples, lower_left, upper_right)
if trigger_type is not None:
vol.set_trigger(threshold, trigger_type.decode())
vol.iterations = iterations
vol.volumes = volumes
vol.atoms = atoms
return vol
@ -277,4 +341,8 @@ class VolumeCalculation(object):
ll_elem.text = ' '.join(str(x) for x in self.lower_left)
ur_elem = ET.SubElement(element, "upper_right")
ur_elem.text = ' '.join(str(x) for x in self.upper_right)
if self.threshold:
trigger_elem = ET.SubElement(element, "threshold")
trigger_elem.set("type", self.trigger_type)
trigger_elem.set("threshold", str(self.threshold))
return element

387
src/volume_calc.cpp
View file

@ -59,7 +59,32 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node)
domain_ids_ = get_node_array<int>(node, "domain_ids");
lower_left_ = get_node_array<double>(node, "lower_left");
upper_right_ = get_node_array<double>(node, "upper_right");
n_samples_ = std::stoi(get_node_value(node, "samples"));
n_samples_ = std::stoull(get_node_value(node, "samples"));
if (check_for_node(node, "threshold")) {
pugi::xml_node threshold_node = node.child("threshold");
threshold_ = std::stod(get_node_value(threshold_node, "threshold"));
if (threshold_ <= 0.0) {
std::stringstream msg;
msg << "Invalid error threshold " << threshold_ << " provided for a volume calculation.";
fatal_error(msg);
}
std::string tmp = get_node_value(threshold_node, "type");
if (tmp == "variance") {
trigger_type_ = TriggerMetric::variance;
} else if (tmp == "std_dev") {
trigger_type_ = TriggerMetric::standard_deviation;
} else if ( tmp == "rel_err") {
trigger_type_ = TriggerMetric::relative_error;
} else {
std::stringstream msg;
msg << "Invalid volume calculation trigger type '" << tmp << "' provided.";
fatal_error(msg);
}
}
// Ensure there are no duplicates by copying elements to a set and then
// comparing the length with the original vector
@ -77,11 +102,12 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
int n = domain_ids_.size();
std::vector<std::vector<int>> master_indices(n); // List of material indices for each domain
std::vector<std::vector<int>> master_hits(n); // Number of hits for each material in each domain
int iterations = 0;
// Divide work over MPI processes
int min_samples = n_samples_ / mpi::n_procs;
int remainder = n_samples_ % mpi::n_procs;
int i_start, i_end;
size_t min_samples = n_samples_ / mpi::n_procs;
size_t remainder = n_samples_ % mpi::n_procs;
size_t i_start, i_end;
if (mpi::rank < remainder) {
i_start = (min_samples + 1)*mpi::rank;
i_end = i_start + min_samples + 1;
@ -90,182 +116,230 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
i_end = i_start + min_samples;
}
#pragma omp parallel
{
// Variables that are private to each thread
std::vector<std::vector<int>> indices(n);
std::vector<std::vector<int>> hits(n);
Particle p;
while (true) {
prn_set_stream(STREAM_VOLUME);
#pragma omp parallel
{
// Variables that are private to each thread
std::vector<std::vector<int>> indices(n);
std::vector<std::vector<int>> hits(n);
Particle p;
// Sample locations and count hits
#pragma omp for
for (int i = i_start; i < i_end; i++) {
set_particle_seed(i);
prn_set_stream(STREAM_VOLUME);
p.n_coord_ = 1;
Position xi {prn(), prn(), prn()};
p.r() = lower_left_ + xi*(upper_right_ - lower_left_);
p.u() = {0.5, 0.5, 0.5};
// Sample locations and count hits
#pragma omp for
for (size_t i = i_start; i < i_end; i++) {
set_particle_seed(iterations * n_samples_ + i);
// If this location is not in the geometry at all, move on to next block
if (!find_cell(&p, false)) continue;
p.n_coord_ = 1;
Position xi {prn(), prn(), prn()};
p.r() = lower_left_ + xi*(upper_right_ - lower_left_);
p.u() = {0.5, 0.5, 0.5};
if (domain_type_ == FILTER_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;
// If this location is not in the geometry at all, move on to next block
if (!find_cell(&p, false)) continue;
if (domain_type_ == FILTER_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_ == FILTER_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_ == FILTER_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_ == FILTER_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;
} else if (domain_type_ == FILTER_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;
}
}
}
}
}
}
// 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
// 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
#ifdef _OPENMP
#pragma omp for ordered schedule(static)
for (int i = 0; i < omp_get_num_threads(); ++i) {
#pragma omp ordered
for (int i_domain = 0; i_domain < n; ++i_domain) {
for (int j = 0; j < indices[i_domain].size(); ++j) {
// Check if this material has been added to the master list and if so,
// accumulate the number of hits
bool already_added = false;
for (int k = 0; k < master_indices[i_domain].size(); k++) {
if (indices[i_domain][j] == master_indices[i_domain][k]) {
master_hits[i_domain][k] += hits[i_domain][j];
already_added = true;
#ifdef _OPENMP
int n_threads = omp_get_num_threads();
#else
int n_threads = 1;
#endif
#pragma omp for ordered schedule(static)
for (int i = 0; i < n_threads; ++i) {
#pragma omp ordered
for (int i_domain = 0; i_domain < n; ++i_domain) {
for (int j = 0; j < indices[i_domain].size(); ++j) {
// Check if this material has been added to the master list and if so,
// accumulate the number of hits
bool already_added = false;
for (int k = 0; k < master_indices[i_domain].size(); k++) {
if (indices[i_domain][j] == master_indices[i_domain][k]) {
master_hits[i_domain][k] += hits[i_domain][j];
already_added = true;
}
}
if (!already_added) {
// If we made it here, the material hasn't yet been added to the master
// list, so add entries to the master indices and master hits lists
master_indices[i_domain].push_back(indices[i_domain][j]);
master_hits[i_domain].push_back(hits[i_domain][j]);
}
}
if (!already_added) {
// If we made it here, the material hasn't yet been added to the master
// list, so add entries to the master indices and master hits lists
master_indices[i_domain].push_back(indices[i_domain][j]);
master_hits[i_domain].push_back(hits[i_domain][j]);
}
}
}
}
#else
master_indices = indices;
master_hits = hits;
#endif
prn_set_stream(STREAM_TRACKING);
} // omp parallel
prn_set_stream(STREAM_TRACKING);
} // omp parallel
// Reduce hits onto master process
// Reduce hits onto master process
// Determine volume of bounding box
Position d {upper_right_ - lower_left_};
double volume_sample = d.x*d.y*d.z;
// Determine volume of bounding box
Position d {upper_right_ - lower_left_};
double volume_sample = d.x*d.y*d.z;
// bump iteration counter and get total number
// of samples at this point
iterations++;
size_t total_samples = iterations * n_samples_;
// Set size for members of the Result struct
std::vector<Result> results(n);
// reset
double trigger_val = -INFTY;
for (int i_domain = 0; i_domain < n; ++i_domain) {
// Get reference to result for this domain
auto& result {results[i_domain]};
// Set size for members of the Result struct
std::vector<Result> results(n);
// 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 = data::nuclides.size();
xt::xtensor<double, 2> atoms({n_nuc, 2}, 0.0);
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 = data::nuclides.size();
xt::xtensor<double, 2> atoms({n_nuc, 2}, 0.0);
#ifdef OPENMC_MPI
if (mpi::master) {
for (int j = 1; j < mpi::n_procs; j++) {
int q;
MPI_Recv(&q, 1, MPI_INTEGER, j, 0, mpi::intracomm, MPI_STATUS_IGNORE);
int buffer[2*q];
MPI_Recv(&buffer[0], 2*q, MPI_INTEGER, j, 1, mpi::intracomm, MPI_STATUS_IGNORE);
for (int k = 0; k < q; ++k) {
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];
break;
if (mpi::master) {
for (int j = 1; j < mpi::n_procs; j++) {
int q;
MPI_Recv(&q, 1, MPI_INTEGER, j, 0, mpi::intracomm, MPI_STATUS_IGNORE);
int buffer[2*q];
MPI_Recv(&buffer[0], 2*q, MPI_INTEGER, j, 1, mpi::intracomm, MPI_STATUS_IGNORE);
for (int k = 0; k < q; ++k) {
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];
break;
}
}
}
}
}
} else {
int q = master_indices[i_domain].size();
int 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];
}
} else {
int q = master_indices[i_domain].size();
int 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_INTEGER, 0, 0, mpi::intracomm);
MPI_Send(&buffer[0], 2*q, MPI_INTEGER, 0, 1, mpi::intracomm);
MPI_Send(&q, 1, MPI_INTEGER, 0, 0, mpi::intracomm);
MPI_Send(&buffer[0], 2*q, MPI_INTEGER, 0, 1, mpi::intracomm);
}
#endif
if (mpi::master) {
int 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;
int i_material = master_indices[i_domain][j];
if (i_material == MATERIAL_VOID) 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;
}
// update max if entry is valid
if (val > 0.0) { trigger_val = std::max(trigger_val, val); }
}
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) {
result.nuclides.push_back(j);
result.atoms.push_back(mean);
result.uncertainty.push_back(stdev);
}
}
}
} // end domain loop
// if no trigger is applied, we're done
if (trigger_type_ == TriggerMetric::not_active) { return results; }
#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
if (mpi::master) {
int 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]) / n_samples_;
double var_f = f*(1.0 - f)/n_samples_;
int i_material = master_indices[i_domain][j];
if (i_material == MATERIAL_VOID) 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) / n_samples_ * volume_sample;
result.volume[1] = std::sqrt(result.volume[0]
* (volume_sample - result.volume[0]) / n_samples_);
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) {
result.nuclides.push_back(j);
result.atoms.push_back(mean);
result.uncertainty.push_back(stdev);
}
}
}
}
return results;
} // end while
}
void VolumeCalculation::to_hdf5(const std::string& filename,
@ -289,6 +363,27 @@ void VolumeCalculation::to_hdf5(const std::string& filename,
write_attribute(file_id, "samples", n_samples_);
write_attribute(file_id, "lower_left", lower_left_);
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, "threshold", threshold_);
std::string trigger_str;
switch(trigger_type_) {
case TriggerMetric::variance:
trigger_str = "variance";
break;
case TriggerMetric::standard_deviation:
trigger_str = "std_dev";
break;
case TriggerMetric::relative_error:
trigger_str = "rel_err";
break;
}
write_attribute(file_id, "trigger_type", trigger_str);
} else {
write_attribute(file_id, "iterations", 1);
}
if (domain_type_ == FILTER_CELL) {
write_attribute(file_id, "domain_type", "cell");
}

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

@ -48,4 +48,28 @@
<lower_left>-1.0 -1.0 -6.0</lower_left>
<upper_right>1.0 1.0 6.0</upper_right>
</volume_calc>
<volume_calc>
<domain_type>cell</domain_type>
<domain_ids>1 2 3</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 threshold="0.1" type="std_dev" />
</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 threshold="0.1" type="rel_err" />
</volume_calc>
<volume_calc>
<domain_type>cell</domain_type>
<domain_ids>1 2 3</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 threshold="0.05" type="variance" />
</volume_calc>
</settings>

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

@ -1,4 +1,7 @@
Volume calculation 0
Trigger Type: None
Trigger threshold: None
Iterations: 1
Domain 1: 31.47+/-0.07 cm^3
Domain 2: 2.093+/-0.031 cm^3
Domain 3: 2.049+/-0.031 cm^3
@ -12,6 +15,9 @@ Domain 3: 2.049+/-0.031 cm^3
6 3 O16 (6.85+/-0.10)e+22
7 3 B10 (6.85+/-0.10)e+18
Volume calculation 1
Trigger Type: None
Trigger threshold: None
Iterations: 1
Domain 1: 4.14+/-0.04 cm^3
Domain 2: 31.47+/-0.07 cm^3
Material Nuclide Atoms
@ -21,6 +27,9 @@ Domain 2: 31.47+/-0.07 cm^3
3 2 U235 (3.482+/-0.008)e+23
4 2 Mo99 (3.482+/-0.008)e+22
Volume calculation 2
Trigger Type: None
Trigger threshold: None
Iterations: 1
Domain 0: 35.61+/-0.07 cm^3
Universe Nuclide Atoms
0 0 H1 (2.770+/-0.029)e+23
@ -28,3 +37,47 @@ Domain 0: 35.61+/-0.07 cm^3
2 0 B10 (1.385+/-0.014)e+19
3 0 U235 (3.482+/-0.008)e+23
4 0 Mo99 (3.482+/-0.008)e+22
Volume calculation 3
Trigger Type: std_dev
Trigger threshold: 0.1
Iterations: 521
Domain 1: 31.47+/-0.10 cm^3
Domain 2: 2.10+/-0.04 cm^3
Domain 3: 2.11+/-0.04 cm^3
Cell Nuclide Atoms
0 1 U235 (3.481+/-0.011)e+23
1 1 Mo99 (3.481+/-0.011)e+22
2 2 H1 (1.403+/-0.029)e+23
3 2 O16 (7.01+/-0.14)e+22
4 2 B10 (7.01+/-0.14)e+18
5 3 H1 (1.412+/-0.029)e+23
6 3 O16 (7.06+/-0.14)e+22
7 3 B10 (7.06+/-0.14)e+18
Volume calculation 4
Trigger Type: rel_err
Trigger threshold: 0.1
Iterations: 10
Domain 1: 4.5+/-0.4 cm^3
Domain 2: 30.5+/-0.7 cm^3
Material Nuclide Atoms
0 1 H1 (3.02+/-0.30)e+23
1 1 O16 (1.51+/-0.15)e+23
2 1 B10 (1.51+/-0.15)e+19
3 2 U235 (3.38+/-0.08)e+23
4 2 Mo99 (3.38+/-0.08)e+22
Volume calculation 5
Trigger Type: variance
Trigger threshold: 0.05
Iterations: 105
Domain 1: 31.51+/-0.22 cm^3
Domain 2: 2.13+/-0.10 cm^3
Domain 3: 2.11+/-0.10 cm^3
Cell Nuclide Atoms
0 1 U235 (3.486+/-0.025)e+23
1 1 Mo99 (3.486+/-0.025)e+22
2 2 H1 (1.42+/-0.06)e+23
3 2 O16 (7.11+/-0.32)e+22
4 2 B10 (7.11+/-0.32)e+18
5 3 H1 (1.41+/-0.06)e+23
6 3 O16 (7.05+/-0.32)e+22
7 3 B10 (7.05+/-0.32)e+18

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

@ -8,6 +8,14 @@ from tests.testing_harness import PyAPITestHarness
class VolumeTest(PyAPITestHarness):
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self.exp_std_dev = 1e-01
self.exp_rel_err = 1e-01
self.exp_variance = 5e-02
def _build_inputs(self):
# Define materials
water = openmc.Material(1)
@ -45,9 +53,18 @@ class VolumeTest(PyAPITestHarness):
vol_calcs = [
openmc.VolumeCalculation(list(root.cells.values()), 100000),
openmc.VolumeCalculation([water, fuel], 100000, ll, ur),
openmc.VolumeCalculation([root], 100000, ll, ur)
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)
]
vol_calcs[3].set_trigger(self.exp_std_dev, 'std_dev')
vol_calcs[4].set_trigger(self.exp_rel_err, 'rel_err')
vol_calcs[5].set_trigger(self.exp_variance, 'variance')
# Define settings
settings = openmc.Settings()
settings.run_mode = 'volume'
@ -62,6 +79,33 @@ class VolumeTest(PyAPITestHarness):
# Read volume calculation results
volume_calc = openmc.VolumeCalculation.from_hdf5(filename)
outstr += 'Trigger Type: {}\n'.format(volume_calc.trigger_type)
outstr += 'Trigger threshold: {}\n'.format(volume_calc.threshold)
outstr += 'Iterations: {}\n'.format(volume_calc.iterations)
if i == 3:
assert(volume_calc.trigger_type == 'std_dev')
assert(volume_calc.threshold == self.exp_std_dev)
elif i == 4:
assert(volume_calc.trigger_type == 'rel_err')
assert(volume_calc.threshold == self.exp_rel_err)
elif i == 5:
assert(volume_calc.trigger_type == 'variance')
assert(volume_calc.threshold == self.exp_variance)
else:
assert(volume_calc.trigger_type == None)
assert(volume_calc.threshold == None)
assert(volume_calc.iterations == 1)
# if a trigger is applied, make sure the calculation satisfies the trigger
for vol in volume_calc.volumes.values():
if volume_calc.trigger_type == 'std_dev':
assert(vol.std_dev <= self.exp_std_dev)
if volume_calc.trigger_type == 'rel_err':
assert(vol.std_dev/vol.nominal_value <= self.exp_rel_err)
if volume_calc.trigger_type == 'variance':
assert(vol.std_dev * vol.std_dev <= self.exp_variance)
# Write cell volumes and total # of atoms for each nuclide
for uid, volume in sorted(volume_calc.volumes.items()):
outstr += 'Domain {}: {} cm^3\n'.format(uid, volume)