diff --git a/docs/source/io_formats/volume.rst b/docs/source/io_formats/volume.rst index 456a06eef..c79255a1c 100644 --- a/docs/source/io_formats/volume.rst +++ b/docs/source/io_formats/volume.rst @@ -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_/** diff --git a/docs/source/usersguide/volume.rst b/docs/source/usersguide/volume.rst index e4bdbf715..c040eb2d7 100644 --- a/docs/source/usersguide/volume.rst +++ b/docs/source/usersguide/volume.rst @@ -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`:: diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 7ee660d8b..1b49da814 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -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 #include +#include namespace openmc { @@ -16,6 +18,7 @@ namespace openmc { //============================================================================== class VolumeCalculation { + public: // Aliases, types struct Result { @@ -23,6 +26,7 @@ public: std::vector nuclides; //!< Index of nuclides std::vector atoms; //!< Number of atoms for each nuclide std::vector 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 domain_ids_; //!< IDs of domains to find volumes of diff --git a/openmc/volume.py b/openmc/volume.py index 91ff829cb..a59fb2acf 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -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 diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 6c1c48932..419bf0ea5 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -59,7 +59,32 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) domain_ids_ = get_node_array(node, "domain_ids"); lower_left_ = get_node_array(node, "lower_left"); upper_right_ = get_node_array(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::execute() const int n = domain_ids_.size(); std::vector> master_indices(n); // List of material indices for each domain std::vector> 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::execute() const i_end = i_start + min_samples; } - #pragma omp parallel - { - // Variables that are private to each thread - std::vector> indices(n); - std::vector> hits(n); - Particle p; + while (true) { - prn_set_stream(STREAM_VOLUME); + #pragma omp parallel + { + // Variables that are private to each thread + std::vector> indices(n); + std::vector> 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 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 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 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 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(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(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(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(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"); } diff --git a/tests/regression_tests/volume_calc/inputs_true.dat b/tests/regression_tests/volume_calc/inputs_true.dat index 607921af2..aaf6d8b00 100644 --- a/tests/regression_tests/volume_calc/inputs_true.dat +++ b/tests/regression_tests/volume_calc/inputs_true.dat @@ -48,4 +48,28 @@ -1.0 -1.0 -6.0 1.0 1.0 6.0 + + cell + 1 2 3 + 100 + -1.0 -1.0 -6.0 + 1.0 1.0 6.0 + + + + material + 1 2 + 100 + -1.0 -1.0 -6.0 + 1.0 1.0 6.0 + + + + cell + 1 2 3 + 100 + -1.0 -1.0 -6.0 + 1.0 1.0 6.0 + + diff --git a/tests/regression_tests/volume_calc/results_true.dat b/tests/regression_tests/volume_calc/results_true.dat index 466139cd6..800356e35 100644 --- a/tests/regression_tests/volume_calc/results_true.dat +++ b/tests/regression_tests/volume_calc/results_true.dat @@ -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 diff --git a/tests/regression_tests/volume_calc/test.py b/tests/regression_tests/volume_calc/test.py index 8ea4ab04e..f7b9a27a0 100644 --- a/tests/regression_tests/volume_calc/test.py +++ b/tests/regression_tests/volume_calc/test.py @@ -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)