From 130b7a5eebea279206ac84138608fa990fe0f6dd Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 01:14:27 -0500 Subject: [PATCH 01/32] Setting up new function implement a stopping criterion for volume calculations. --- include/openmc/volume_calc.h | 29 +++++++++++++++++++++++++++ src/volume_calc.cpp | 38 ++++++++++++++++++++++++++++++++++-- 2 files changed, 65 insertions(+), 2 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 7ee660d8b..c7f6413b0 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -8,6 +8,7 @@ #include #include +#include namespace openmc { @@ -23,6 +24,31 @@ public: std::vector nuclides; //!< Index of nuclides std::vector atoms; //!< Number of atoms for each nuclide std::vector uncertainty; //!< Uncertainty on number of atoms + size_t num_samples; + + Result& operator +=( const Result& other) { + Expects(volume.size() == other.volume.size()); + Expects(atoms.size() == atoms.size()); + + size_t total_samples = num_samples + other.num_samples; + + for (int i = 0; i < volume.size(); i++) { + // average volume results + volume[0] = (num_samples * volume[0] + other.num_samples * other.volume[0]) / total_samples; + // propagate error + volume[1] = std::sqrt(num_samples * volume[1] *volume[1] + other.num_samples * other.volume[1] * other.volume[1]) / total_samples; + } + + for (int i = 0; i < atoms.size(); i++) { + atoms[i] = (num_samples * atoms[i] + other.num_samples * other.atoms[i]) / total_samples; + uncertainty[i] = std::sqrt(num_samples * uncertainty[i] * uncertainty[i] + other.num_samples * other.uncertainty[i] * other.uncertainty[i]) / total_samples; + } + + num_samples = total_samples; + + return *this; + + } }; // Results for a single domain // Constructors @@ -36,6 +62,8 @@ public: //! \return Vector of results for each user-specified domain std::vector execute() const; + std::vector _execute(size_t seed_offset = 0) const; + //! \brief Write volume calculation results to HDF5 file // //! \param[in] filename Path to HDF5 file to write @@ -45,6 +73,7 @@ public: // Data members int domain_type_; //!< Type of domain (cell, material, etc.) int n_samples_; //!< Number of samples to use + int seed_offset_; 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/src/volume_calc.cpp b/src/volume_calc.cpp index 6c1c48932..bd6ada96a 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -71,7 +71,40 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) } -std::vector VolumeCalculation::execute() const +std::vector VolumeCalculation::execute() const { + + std::vector results; + size_t offset = 0; + + results = _execute(offset); + offset += n_samples_; + + double max_err = -INFTY; + for (int i = 0; i < results.size(); i++) { + max_err = std::max(max_err, results[i].volume[1]); + } + + double error_limit = 1E-05; + int iters = 1; + while (max_err > error_limit) { + std::cout << "Iter " << iters++ << std::endl; + std::vector tmp = _execute(offset); + max_err = -INFTY; + for (int i = 0; i < results.size(); i++) { + auto& result = results[i]; + result += tmp[i]; + max_err = std::max(max_err, result.volume[1]); + } + + offset += n_samples_; + + std::cout << "Max error: " << max_err << std::endl; + } + + return results; +} + +std::vector VolumeCalculation::_execute(size_t seed_offset) const { // Shared data that is collected from all threads int n = domain_ids_.size(); @@ -102,7 +135,7 @@ std::vector VolumeCalculation::execute() const // Sample locations and count hits #pragma omp for for (int i = i_start; i < i_end; i++) { - set_particle_seed(i); + set_particle_seed(seed_offset + i); p.n_coord_ = 1; Position xi {prn(), prn(), prn()}; @@ -248,6 +281,7 @@ std::vector VolumeCalculation::execute() const 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_); + result.num_samples = n_samples_; for (int j = 0; j < n_nuc; ++j) { // Determine total number of atoms. At this point, we have values in From 1c5218aad6cc14772c454360ced1da08ad69c645 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 01:26:21 -0500 Subject: [PATCH 02/32] Cleaning up operator a bit. --- include/openmc/volume_calc.h | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index c7f6413b0..b5e847742 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -30,25 +30,41 @@ public: Expects(volume.size() == other.volume.size()); Expects(atoms.size() == atoms.size()); + auto& a_samples = num_samples; + auto& b_samples = other.num_samples; + size_t total_samples = num_samples + other.num_samples; for (int i = 0; i < volume.size(); i++) { - // average volume results - volume[0] = (num_samples * volume[0] + other.num_samples * other.volume[0]) / total_samples; + // calculate weighted average of volume results + auto& a_vol = volume[0]; + auto& b_vol = other.volume[1]; + volume[0] = (a_samples * a_vol + b_samples * b_vol) / total_samples; + // propagate error - volume[1] = std::sqrt(num_samples * volume[1] *volume[1] + other.num_samples * other.volume[1] * other.volume[1]) / total_samples; + auto& a_err = volume[1]; + auto& b_err = other.volume[1]; + volume[1] = std::sqrt(a_samples * a_err * a_err + b_samples * b_err * b_err) / total_samples; } for (int i = 0; i < atoms.size(); i++) { - atoms[i] = (num_samples * atoms[i] + other.num_samples * other.atoms[i]) / total_samples; - uncertainty[i] = std::sqrt(num_samples * uncertainty[i] * uncertainty[i] + other.num_samples * other.uncertainty[i] * other.uncertainty[i]) / total_samples; + // calculate weighted average of atom results + auto& a_atoms = atoms[i]; + auto& b_atoms = other.atoms[i]; + atoms[i] = (a_samples * a_atoms + b_samples * b_atoms) / total_samples; + + // propagate error + auto& a_err = uncertainty[i]; + auto& b_err = other.uncertainty[i]; + uncertainty[i] = std::sqrt(a_samples * a_err * a_err + b_samples * b_err * b_err) / total_samples; } + // update number of samples on the returned set of results; num_samples = total_samples; - - return *this; + return *this; } + }; // Results for a single domain // Constructors From d3bcfeda0f667650c2e6d5fdb948f5b0db8f26af Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 01:41:07 -0500 Subject: [PATCH 03/32] Cleaning up while loop. --- src/volume_calc.cpp | 36 ++++++++++++++---------------------- 1 file changed, 14 insertions(+), 22 deletions(-) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index bd6ada96a..40b5f08dd 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -73,32 +73,24 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) std::vector VolumeCalculation::execute() const { - std::vector results; - size_t offset = 0; - - results = _execute(offset); - offset += n_samples_; - - double max_err = -INFTY; - for (int i = 0; i < results.size(); i++) { - max_err = std::max(max_err, results[i].volume[1]); - } - + std::vector results = _execute(); + size_t offset = n_samples_; double error_limit = 1E-05; - int iters = 1; - while (max_err > error_limit) { - std::cout << "Iter " << iters++ << std::endl; - std::vector tmp = _execute(offset); + double max_err; + + while (true) { + // check maximum error value for all domains max_err = -INFTY; - for (int i = 0; i < results.size(); i++) { - auto& result = results[i]; - result += tmp[i]; - max_err = std::max(max_err, result.volume[1]); - } - + for (const auto& result : results) { max_err = std::max(max_err, result.volume[1]); } + + // exit once we're below our error limit + if (max_err <= error_limit) { break; } + + std::vector tmp = _execute(offset); offset += n_samples_; - std::cout << "Max error: " << max_err << std::endl; + // update current results + for (int i = 0; i < results.size(); i++) { results[i] += tmp[i]; } } return results; From 5109bd0de1222f3096f9f511b194c35e25a69696 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 02:13:12 -0500 Subject: [PATCH 04/32] Adding error trigger attribute to volume calc. --- include/openmc/volume_calc.h | 4 ++-- src/volume_calc.cpp | 32 +++++++++++++++++++++++++------- 2 files changed, 27 insertions(+), 9 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index b5e847742..2a8273f9b 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -88,8 +88,8 @@ public: // Data members int domain_type_; //!< Type of domain (cell, material, etc.) - int n_samples_; //!< Number of samples to use - int seed_offset_; + size_t n_samples_; //!< Number of samples to use + double error_trigger_ {-1.0}; //!< Error below which the calculation will stop 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/src/volume_calc.cpp b/src/volume_calc.cpp index 40b5f08dd..a57305db9 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -59,7 +59,23 @@ 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")); + std::stringstream size_t_stream(get_node_value(node, "samples")); + size_t_stream >> n_samples_; + if (size_t_stream.fail()) { + std::stringstream msg; + msg << "Could not read number of samples (" + << size_t_stream.str() << ")\n"; + fatal_error(msg); + } + + if (check_for_node(node, "error_trigger")) { + error_trigger_ = std::stod(get_node_value(node, "error_trigger")); + if (error_trigger_ <= 0.0) { + std::stringstream msg; + msg << "Invalid error trigger " << error_trigger_ << " provided for a volume calculation."; + fatal_error(msg); + } + } // Ensure there are no duplicates by copying elements to a set and then // comparing the length with the original vector @@ -74,8 +90,10 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) std::vector VolumeCalculation::execute() const { std::vector results = _execute(); + + if (error_trigger_ <= 0.0) { return results; } + size_t offset = n_samples_; - double error_limit = 1E-05; double max_err; while (true) { @@ -84,7 +102,7 @@ std::vector VolumeCalculation::execute() const { for (const auto& result : results) { max_err = std::max(max_err, result.volume[1]); } // exit once we're below our error limit - if (max_err <= error_limit) { break; } + if (max_err <= error_trigger_) { break; } std::vector tmp = _execute(offset); offset += n_samples_; @@ -104,9 +122,9 @@ std::vector VolumeCalculation::_execute(size_t seed_o std::vector> master_hits(n); // Number of hits for each material in each domain // 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; @@ -126,7 +144,7 @@ std::vector VolumeCalculation::_execute(size_t seed_o // Sample locations and count hits #pragma omp for - for (int i = i_start; i < i_end; i++) { + for (size_t i = i_start; i < i_end; i++) { set_particle_seed(seed_offset + i); p.n_coord_ = 1; From 4b690ab7518e6de2376ab3ebd01a12c7e463cba1 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 02:18:16 -0500 Subject: [PATCH 05/32] Fixing a bug. --- include/openmc/volume_calc.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 2a8273f9b..18aa2e2d1 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -38,7 +38,7 @@ public: for (int i = 0; i < volume.size(); i++) { // calculate weighted average of volume results auto& a_vol = volume[0]; - auto& b_vol = other.volume[1]; + auto& b_vol = other.volume[0]; volume[0] = (a_samples * a_vol + b_samples * b_vol) / total_samples; // propagate error From 507bc0d9339a896dd00ed80ed00f0e5f6724eca7 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 10:23:11 -0500 Subject: [PATCH 06/32] Exposing volume trigger through Python API. --- include/openmc/volume_calc.h | 12 ++++++--- openmc/volume.py | 31 +++++++++++++++++++--- src/volume_calc.cpp | 16 ++++++----- tests/regression_tests/volume_calc/test.py | 2 +- 4 files changed, 48 insertions(+), 13 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 18aa2e2d1..4e2b2a2a6 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -78,8 +78,6 @@ public: //! \return Vector of results for each user-specified domain std::vector execute() const; - std::vector _execute(size_t seed_offset = 0) const; - //! \brief Write volume calculation results to HDF5 file // //! \param[in] filename Path to HDF5 file to write @@ -89,7 +87,7 @@ public: // Data members int domain_type_; //!< Type of domain (cell, material, etc.) size_t n_samples_; //!< Number of samples to use - double error_trigger_ {-1.0}; //!< Error below which the calculation will stop + double trigger_ {-1.0}; //!< Error threshold for domain volumes 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 @@ -104,6 +102,14 @@ private: void check_hit(int i_material, std::vector& indices, std::vector& hits) const; + //! \brief Perform calculation for domain volumes and average nuclide density + //! using n_samples_ + // + //! \param[in] seed_offset Seed offset used for independent calculations + //! \return Vector of results for each user-specified domain + std::vector _execute(size_t seed_offset = 0) const; + + }; //============================================================================== diff --git a/openmc/volume.py b/openmc/volume.py index 91ff829cb..189287bef 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -12,7 +12,7 @@ from uncertainties import ufloat import openmc import openmc.checkvalue as cv -_VERSION_VOLUME = 1 +_VERSION_VOLUME = 2 class VolumeCalculation(object): @@ -32,6 +32,8 @@ class VolumeCalculation(object): Upper-right coordinates of bounding box used to sample points. If this argument is not supplied, an attempt is made to automatically determine a bounding box. + trigger : float + Threshold for the maxmimum standard deviation of volumes Attributes ---------- @@ -45,6 +47,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 + trigger : 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 @@ -57,9 +61,10 @@ class VolumeCalculation(object): """ def __init__(self, domains, samples, lower_left=None, - upper_right=None): + upper_right=None, trigger=None): self._atoms = {} self._volumes = {} + self._trigger = None cv.check_type('domains', domains, Iterable, (openmc.Cell, openmc.Material, openmc.Universe)) @@ -72,6 +77,9 @@ class VolumeCalculation(object): self.ids = [d.id for d in domains] self.samples = samples + + if trigger is not None: + self.trigger = trigger if lower_left is not None: if upper_right is None: @@ -123,6 +131,10 @@ class VolumeCalculation(object): def upper_right(self): return self._upper_right + @property + def trigger(self): + return self._trigger + @property def domain_type(self): return self._domain_type @@ -170,6 +182,12 @@ class VolumeCalculation(object): cv.check_length(name, upper_right, 3) self._upper_right = upper_right + @trigger.setter + def trigger(self, trigger): + name = 'Volume std. dev. trigger' + cv.check_type(name, trigger, Real) + self._trigger = trigger + @volumes.setter def volumes(self, volumes): cv.check_type('volumes', volumes, Mapping) @@ -202,6 +220,10 @@ class VolumeCalculation(object): samples = f.attrs['samples'] lower_left = f.attrs['lower_left'] upper_right = f.attrs['upper_right'] + trigger = f.attrs['trigger'] + + if trigger == -1.0: + trigger = None volumes = {} atoms = {} @@ -232,7 +254,7 @@ class VolumeCalculation(object): domains = [openmc.Universe(uid) for uid in ids] # Instantiate the class and assign results - vol = cls(domains, samples, lower_left, upper_right) + vol = cls(domains, samples, lower_left, upper_right, trigger) vol.volumes = volumes vol.atoms = atoms return vol @@ -277,4 +299,7 @@ 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.trigger: + trigger_elem = ET.SubElement(element, "trigger") + trigger_elem.text = str(self.trigger) return element diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index a57305db9..358d8cfd4 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -68,11 +68,11 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) fatal_error(msg); } - if (check_for_node(node, "error_trigger")) { - error_trigger_ = std::stod(get_node_value(node, "error_trigger")); - if (error_trigger_ <= 0.0) { + if (check_for_node(node, "trigger")) { + trigger_ = std::stod(get_node_value(node, "trigger")); + if (trigger_ <= 0.0) { std::stringstream msg; - msg << "Invalid error trigger " << error_trigger_ << " provided for a volume calculation."; + msg << "Invalid error threshold " << trigger_ << " provided for a volume calculation."; fatal_error(msg); } } @@ -89,9 +89,11 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) std::vector VolumeCalculation::execute() const { + // execute the calculation once std::vector results = _execute(); - if (error_trigger_ <= 0.0) { return results; } + // if no std. dev. threshold is set, return these resuls + if (trigger_ == -1.0) { return results; } size_t offset = n_samples_; double max_err; @@ -102,8 +104,9 @@ std::vector VolumeCalculation::execute() const { for (const auto& result : results) { max_err = std::max(max_err, result.volume[1]); } // exit once we're below our error limit - if (max_err <= error_trigger_) { break; } + if (max_err <= trigger_) { break; } + // perform the calculation std::vector tmp = _execute(offset); offset += n_samples_; @@ -333,6 +336,7 @@ 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_attribute(file_id, "trigger", trigger_); if (domain_type_ == FILTER_CELL) { write_attribute(file_id, "domain_type", "cell"); } diff --git a/tests/regression_tests/volume_calc/test.py b/tests/regression_tests/volume_calc/test.py index 8ea4ab04e..2affc4085 100644 --- a/tests/regression_tests/volume_calc/test.py +++ b/tests/regression_tests/volume_calc/test.py @@ -74,4 +74,4 @@ class VolumeTest(PyAPITestHarness): def test_volume_calc(): harness = VolumeTest('') - harness.main() + harness._build_inputs() From 48fab6dcac29a9935a575d26d158503de420379b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 10:50:32 -0500 Subject: [PATCH 07/32] Only write threshold to file if used. --- include/openmc/volume_calc.h | 2 +- openmc/volume.py | 46 ++++++++++++++++++------------------ src/volume_calc.cpp | 17 +++++++------ 3 files changed, 34 insertions(+), 31 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 4e2b2a2a6..33c8cbfea 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -87,7 +87,7 @@ public: // Data members int domain_type_; //!< Type of domain (cell, material, etc.) size_t n_samples_; //!< Number of samples to use - double trigger_ {-1.0}; //!< Error threshold for domain volumes + double threshold_ {-1.0}; //!< Error threshold for domain volumes 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 189287bef..bcc46d19e 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -12,7 +12,7 @@ from uncertainties import ufloat import openmc import openmc.checkvalue as cv -_VERSION_VOLUME = 2 +_VERSION_VOLUME = 1 class VolumeCalculation(object): @@ -32,7 +32,7 @@ class VolumeCalculation(object): Upper-right coordinates of bounding box used to sample points. If this argument is not supplied, an attempt is made to automatically determine a bounding box. - trigger : float + threshold : float Threshold for the maxmimum standard deviation of volumes Attributes @@ -47,7 +47,7 @@ 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 - trigger : float + 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 @@ -61,10 +61,10 @@ class VolumeCalculation(object): """ def __init__(self, domains, samples, lower_left=None, - upper_right=None, trigger=None): + upper_right=None, threshold=None): self._atoms = {} self._volumes = {} - self._trigger = None + self._threshold = None cv.check_type('domains', domains, Iterable, (openmc.Cell, openmc.Material, openmc.Universe)) @@ -78,8 +78,8 @@ class VolumeCalculation(object): self.samples = samples - if trigger is not None: - self.trigger = trigger + if threshold is not None: + self.threshold = threshold if lower_left is not None: if upper_right is None: @@ -132,8 +132,8 @@ class VolumeCalculation(object): return self._upper_right @property - def trigger(self): - return self._trigger + def threshold(self): + return self._threshold @property def domain_type(self): @@ -182,11 +182,11 @@ class VolumeCalculation(object): cv.check_length(name, upper_right, 3) self._upper_right = upper_right - @trigger.setter - def trigger(self, trigger): - name = 'Volume std. dev. trigger' - cv.check_type(name, trigger, Real) - self._trigger = trigger + @threshold.setter + def threshold(self, threshold): + name = 'volume std. dev. threshold' + cv.check_type(name, threshold, Real) + self._threshold = threshold @volumes.setter def volumes(self, volumes): @@ -220,11 +220,11 @@ class VolumeCalculation(object): samples = f.attrs['samples'] lower_left = f.attrs['lower_left'] upper_right = f.attrs['upper_right'] - trigger = f.attrs['trigger'] - - if trigger == -1.0: - trigger = None - + try: + threshold = f.attrs['threshold'] + except KeyError: + threshold = None + volumes = {} atoms = {} ids = [] @@ -254,7 +254,7 @@ class VolumeCalculation(object): domains = [openmc.Universe(uid) for uid in ids] # Instantiate the class and assign results - vol = cls(domains, samples, lower_left, upper_right, trigger) + vol = cls(domains, samples, lower_left, upper_right, threshold) vol.volumes = volumes vol.atoms = atoms return vol @@ -299,7 +299,7 @@ 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.trigger: - trigger_elem = ET.SubElement(element, "trigger") - trigger_elem.text = str(self.trigger) + if self.threshold: + threshold_elem = ET.SubElement(element, "threshold") + threshold_elem.text = str(self.threshold) return element diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 358d8cfd4..9b2ef7c51 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -68,11 +68,11 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) fatal_error(msg); } - if (check_for_node(node, "trigger")) { - trigger_ = std::stod(get_node_value(node, "trigger")); - if (trigger_ <= 0.0) { + if (check_for_node(node, "threshold")) { + threshold_ = std::stod(get_node_value(node, "threshold")); + if (threshold_ <= 0.0) { std::stringstream msg; - msg << "Invalid error threshold " << trigger_ << " provided for a volume calculation."; + msg << "Invalid error threshold " << threshold_ << " provided for a volume calculation."; fatal_error(msg); } } @@ -93,7 +93,7 @@ std::vector VolumeCalculation::execute() const { std::vector results = _execute(); // if no std. dev. threshold is set, return these resuls - if (trigger_ == -1.0) { return results; } + if (threshold_ == -1.0) { return results; } size_t offset = n_samples_; double max_err; @@ -104,7 +104,7 @@ std::vector VolumeCalculation::execute() const { for (const auto& result : results) { max_err = std::max(max_err, result.volume[1]); } // exit once we're below our error limit - if (max_err <= trigger_) { break; } + if (max_err <= threshold_) { break; } // perform the calculation std::vector tmp = _execute(offset); @@ -336,7 +336,10 @@ 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_attribute(file_id, "trigger", trigger_); + if (threshold_ != -1.0) { + write_attribute(file_id, "threshold", threshold_); + } + if (domain_type_ == FILTER_CELL) { write_attribute(file_id, "domain_type", "cell"); } From 0914330f8d773b62f321b6ceb21a218b95247bcb Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 10:58:11 -0500 Subject: [PATCH 08/32] Updating volume documentation. --- docs/source/io_formats/volume.rst | 1 + docs/source/usersguide/volume.rst | 12 +++++++++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/docs/source/io_formats/volume.rst b/docs/source/io_formats/volume.rst index 456a06eef..e421de566 100644 --- a/docs/source/io_formats/volume.rst +++ b/docs/source/io_formats/volume.rst @@ -23,6 +23,7 @@ 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 **/domain_/** diff --git a/docs/source/usersguide/volume.rst b/docs/source/usersguide/volume.rst index e4bdbf715..f34710232 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 for the uncertainty in volume estimates can be specified using +::attr::`openmc.VolumeCalculation.threshold` :: + + vol_calc.threshold = 1E-05 + +If a threshold is provided, calculations will be performed iteratively using the +number of samples specified on the calculation until all volume estimates have a +standard deviation lower than this 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`:: @@ -66,4 +76,4 @@ After the results are loaded, volume estimates will be stored in :attr:`VolumeCalculation.volumes`. There is also a :attr:`VolumeCalculation.atoms_dataframe` attribute that shows stochastic estimates of the number of atoms of each type of nuclide within the specified -domains along with their uncertainties. +domains along with their uncertainties. \ No newline at end of file From 9415bbb71550772c7152836ab1452b635ad5f720 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 11:02:23 -0500 Subject: [PATCH 09/32] Adding check for valid threshold value in Python API. --- openmc/volume.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/openmc/volume.py b/openmc/volume.py index bcc46d19e..f185620cd 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -186,6 +186,9 @@ class VolumeCalculation(object): def threshold(self, threshold): name = 'volume std. dev. threshold' cv.check_type(name, threshold, Real) + if (threshold <= 0.0): + raise ValueError("Invalid value '{}' (<= 0.0) provided for volume " + "calculation threshold.".format(threshold)) self._threshold = threshold @volumes.setter @@ -224,7 +227,7 @@ class VolumeCalculation(object): threshold = f.attrs['threshold'] except KeyError: threshold = None - + volumes = {} atoms = {} ids = [] From 8d1e8711f564f4f7e97fb79265bf4fbbaa0174ed Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 12:33:00 -0500 Subject: [PATCH 10/32] Adding volume trigger via different methods. --- include/openmc/volume_calc.h | 8 +++++++ src/volume_calc.cpp | 44 ++++++++++++++++++++++++++++++++---- 2 files changed, 47 insertions(+), 5 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 33c8cbfea..e2cf1b3ef 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -12,11 +12,18 @@ namespace openmc { +enum class ThresholdType { + VARIANCE = 0, + STD_DEV = 1, + REL_ERR = 2 +}; + //============================================================================== // Volume calculation class //============================================================================== class VolumeCalculation { + public: // Aliases, types struct Result { @@ -88,6 +95,7 @@ public: int 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 + ThresholdType threshold_type_; 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/src/volume_calc.cpp b/src/volume_calc.cpp index 9b2ef7c51..870721f3e 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -75,6 +75,21 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) msg << "Invalid error threshold " << threshold_ << " provided for a volume calculation."; fatal_error(msg); } + + pugi::xml_node threshold_node = node.child("threshold"); + std::string tmp = get_node_value(threshold_node, "type"); + if (tmp == "variance") { + threshold_type_ = ThresholdType::VARIANCE; + } else if (tmp == "std_dev") { + threshold_type_ = ThresholdType::STD_DEV; + } else if ( tmp == "rel_err") { + threshold_type_ = ThresholdType::REL_ERR; + } 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 @@ -96,15 +111,30 @@ std::vector VolumeCalculation::execute() const { if (threshold_ == -1.0) { return results; } size_t offset = n_samples_; - double max_err; + double max_val; while (true) { // check maximum error value for all domains - max_err = -INFTY; - for (const auto& result : results) { max_err = std::max(max_err, result.volume[1]); } + max_val = -INFTY; + for (const auto& result : results) { + double val; + switch (threshold_type_) { + case ThresholdType::STD_DEV: + val = result.volume[1]; + break; + case ThresholdType::REL_ERR: + val = result.volume[1] / result.volume[0]; + break; + case ThresholdType::VARIANCE: + val = result.volume[1] * result.volume[1]; + break; + } + // update max + max_val = std::max(max_val, val); + } // exit once we're below our error limit - if (max_err <= threshold_) { break; } + if (max_val <= threshold_) { break; } // perform the calculation std::vector tmp = _execute(offset); @@ -307,6 +337,10 @@ std::vector VolumeCalculation::_execute(size_t seed_o result.nuclides.push_back(j); result.atoms.push_back(mean); result.uncertainty.push_back(stdev); + } else { + result.nuclides.push_back(j); + result.atoms.push_back(0.0); + result.uncertainty.push_back(INFTY); } } } @@ -339,7 +373,7 @@ void VolumeCalculation::to_hdf5(const std::string& filename, if (threshold_ != -1.0) { write_attribute(file_id, "threshold", threshold_); } - + if (domain_type_ == FILTER_CELL) { write_attribute(file_id, "domain_type", "cell"); } From 85e25d79eafd65a064c3afeea3a77e30e50b0047 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 13:12:31 -0500 Subject: [PATCH 11/32] Finishing addition of a trigger type. --- include/openmc/volume_calc.h | 2 +- src/volume_calc.cpp | 21 +++++++++++++++++---- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index e2cf1b3ef..a7dd75009 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -95,7 +95,7 @@ public: int 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 - ThresholdType threshold_type_; + ThresholdType trigger_type_; 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/src/volume_calc.cpp b/src/volume_calc.cpp index 870721f3e..f352e71f8 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -79,11 +79,11 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) pugi::xml_node threshold_node = node.child("threshold"); std::string tmp = get_node_value(threshold_node, "type"); if (tmp == "variance") { - threshold_type_ = ThresholdType::VARIANCE; + trigger_type_ = ThresholdType::VARIANCE; } else if (tmp == "std_dev") { - threshold_type_ = ThresholdType::STD_DEV; + trigger_type_ = ThresholdType::STD_DEV; } else if ( tmp == "rel_err") { - threshold_type_ = ThresholdType::REL_ERR; + trigger_type_ = ThresholdType::REL_ERR; } else { std::stringstream msg; msg << "Invalid volume calculation trigger type '" << tmp << "' provided."; @@ -118,7 +118,7 @@ std::vector VolumeCalculation::execute() const { max_val = -INFTY; for (const auto& result : results) { double val; - switch (threshold_type_) { + switch (trigger_type_) { case ThresholdType::STD_DEV: val = result.volume[1]; break; @@ -372,6 +372,19 @@ void VolumeCalculation::to_hdf5(const std::string& filename, write_attribute(file_id, "upper_right", upper_right_); if (threshold_ != -1.0) { write_attribute(file_id, "threshold", threshold_); + + switch(trigger_type_) { + case ThresholdType::VARIANCE: + write_attribute(file_id, "trigger_type", "variance"); + break; + case ThresholdType::STD_DEV: + write_attribute(file_id, "trigger_type", "std_dev"); + break; + case ThresholdType::REL_ERR: + write_attribute(file_id, "trigger_type", "rel_err"); + break; + } + } if (domain_type_ == FILTER_CELL) { From 8aed7b019abc7666d553b4f508deee2a5fb9cf94 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 13:12:56 -0500 Subject: [PATCH 12/32] Adding support for trigger type in the Python API. --- openmc/volume.py | 50 ++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 42 insertions(+), 8 deletions(-) diff --git a/openmc/volume.py b/openmc/volume.py index f185620cd..a72f32ab3 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -32,8 +32,7 @@ class VolumeCalculation(object): Upper-right coordinates of bounding box used to sample points. If this argument is not supplied, an attempt is made to automatically determine a bounding box. - threshold : float - Threshold for the maxmimum standard deviation of volumes + Attributes ---------- @@ -58,13 +57,17 @@ 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 """ - def __init__(self, domains, samples, lower_left=None, - upper_right=None, threshold=None): + def __init__(self, domains, samples, lower_left=None, upper_right=None): self._atoms = {} self._volumes = {} self._threshold = None + self._trigger_type = None cv.check_type('domains', domains, Iterable, (openmc.Cell, openmc.Material, openmc.Universe)) @@ -78,9 +81,6 @@ class VolumeCalculation(object): self.samples = samples - if threshold is not None: - self.threshold = threshold - if lower_left is not None: if upper_right is None: raise ValueError('Both lower-left and upper-right coordinates ' @@ -135,6 +135,10 @@ class VolumeCalculation(object): def threshold(self): return self._threshold + @property + def trigger_type(self): + return self._trigger_type + @property def domain_type(self): return self._domain_type @@ -191,6 +195,12 @@ class VolumeCalculation(object): "calculation threshold.".format(threshold)) 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 + @volumes.setter def volumes(self, volumes): cv.check_type('volumes', volumes, Mapping) @@ -201,6 +211,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. @@ -223,11 +246,17 @@ class VolumeCalculation(object): samples = f.attrs['samples'] lower_left = f.attrs['lower_left'] upper_right = f.attrs['upper_right'] + try: threshold = f.attrs['threshold'] except KeyError: threshold = None + try: + trigger_type = f.attrs['trigger_type'].decode() + except KeyError: + trigger_type = None + volumes = {} atoms = {} ids = [] @@ -257,7 +286,11 @@ class VolumeCalculation(object): domains = [openmc.Universe(uid) for uid in ids] # Instantiate the class and assign results - vol = cls(domains, samples, lower_left, upper_right, threshold) + vol = cls(domains, samples, lower_left, upper_right) + + if threshold is not None: + vol.set_trigger(threshold, trigger_type) + vol.volumes = volumes vol.atoms = atoms return vol @@ -305,4 +338,5 @@ class VolumeCalculation(object): if self.threshold: threshold_elem = ET.SubElement(element, "threshold") threshold_elem.text = str(self.threshold) + threshold_elem.set("type", self.trigger_type) return element From 1d27972d5040ad869818492f1e3efbaafe83800e Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 5 Oct 2019 13:25:25 -0500 Subject: [PATCH 13/32] Updating volume documentation to include info about volume trigger types. --- docs/source/io_formats/volume.rst | 1 + docs/source/usersguide/volume.rst | 12 ++++++------ 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/docs/source/io_formats/volume.rst b/docs/source/io_formats/volume.rst index e421de566..c79255a1c 100644 --- a/docs/source/io_formats/volume.rst +++ b/docs/source/io_formats/volume.rst @@ -24,6 +24,7 @@ The current version of the volume file format is 1.0. - **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 f34710232..d35386b38 100644 --- a/docs/source/usersguide/volume.rst +++ b/docs/source/usersguide/volume.rst @@ -37,15 +37,15 @@ 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 for the uncertainty in volume estimates can be specified using -::attr::`openmc.VolumeCalculation.threshold` :: +A threshold can be applied for the calculation's variance, standard deviation, +or relative error of volume estimates using ::attr::`openmc.VolumeCalculation.set_trigger`:: - vol_calc.threshold = 1E-05 + 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 estimates have a -standard deviation lower than this value. If no threshold is provided, the -calculation will run the number of samples specified once and return the result. +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`:: From f3fa1f25782fa1ebc8c8427f7f950be4abf053b7 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 7 Oct 2019 00:08:47 -0500 Subject: [PATCH 14/32] Fixing issue with INF error values. Adding test case. --- openmc/volume.py | 6 +-- src/volume_calc.cpp | 22 ++++---- .../volume_calc/inputs_true.dat | 8 +++ .../volume_calc/results_true.dat | 54 +++++++++++++++---- tests/regression_tests/volume_calc/test.py | 12 ++++- 5 files changed, 76 insertions(+), 26 deletions(-) diff --git a/openmc/volume.py b/openmc/volume.py index a72f32ab3..3996ce9dd 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -336,7 +336,7 @@ class VolumeCalculation(object): ur_elem = ET.SubElement(element, "upper_right") ur_elem.text = ' '.join(str(x) for x in self.upper_right) if self.threshold: - threshold_elem = ET.SubElement(element, "threshold") - threshold_elem.text = str(self.threshold) - threshold_elem.set("type", self.trigger_type) + 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 f352e71f8..f65cb89ed 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -69,19 +69,20 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) } if (check_for_node(node, "threshold")) { - threshold_ = std::stod(get_node_value(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); } - pugi::xml_node threshold_node = node.child("threshold"); std::string tmp = get_node_value(threshold_node, "type"); if (tmp == "variance") { trigger_type_ = ThresholdType::VARIANCE; } else if (tmp == "std_dev") { - trigger_type_ = ThresholdType::STD_DEV; + trigger_type_ = ThresholdType::STD_DEV; } else if ( tmp == "rel_err") { trigger_type_ = ThresholdType::REL_ERR; } else { @@ -116,11 +117,11 @@ std::vector VolumeCalculation::execute() const { while (true) { // check maximum error value for all domains max_val = -INFTY; - for (const auto& result : results) { + for (const auto& result : results) { double val; switch (trigger_type_) { case ThresholdType::STD_DEV: - val = result.volume[1]; + val = result.volume[1]; break; case ThresholdType::REL_ERR: val = result.volume[1] / result.volume[0]; @@ -129,8 +130,9 @@ std::vector VolumeCalculation::execute() const { val = result.volume[1] * result.volume[1]; break; } - // update max - max_val = std::max(max_val, val); + // update max if entry is valid + if (val > 0.0) { max_val = std::max(max_val, val); } + } // exit once we're below our error limit @@ -306,7 +308,7 @@ std::vector VolumeCalculation::_execute(size_t seed_o 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_; + double var_f = f*(1.0 - f) / n_samples_; int i_material = master_indices[i_domain][j]; if (i_material == MATERIAL_VOID) continue; @@ -340,7 +342,7 @@ std::vector VolumeCalculation::_execute(size_t seed_o } else { result.nuclides.push_back(j); result.atoms.push_back(0.0); - result.uncertainty.push_back(INFTY); + result.uncertainty.push_back(0.0); } } } @@ -384,7 +386,7 @@ void VolumeCalculation::to_hdf5(const std::string& filename, write_attribute(file_id, "trigger_type", "rel_err"); break; } - + } if (domain_type_ == FILTER_CELL) { diff --git a/tests/regression_tests/volume_calc/inputs_true.dat b/tests/regression_tests/volume_calc/inputs_true.dat index 607921af2..1af345405 100644 --- a/tests/regression_tests/volume_calc/inputs_true.dat +++ b/tests/regression_tests/volume_calc/inputs_true.dat @@ -48,4 +48,12 @@ -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..4fe3b4b6e 100644 --- a/tests/regression_tests/volume_calc/results_true.dat +++ b/tests/regression_tests/volume_calc/results_true.dat @@ -2,15 +2,22 @@ Volume calculation 0 Domain 1: 31.47+/-0.07 cm^3 Domain 2: 2.093+/-0.031 cm^3 Domain 3: 2.049+/-0.031 cm^3 - Cell Nuclide Atoms -0 1 U235 (3.482+/-0.008)e+23 -1 1 Mo99 (3.482+/-0.008)e+22 -2 2 H1 (1.400+/-0.021)e+23 -3 2 O16 (7.00+/-0.10)e+22 -4 2 B10 (7.00+/-0.10)e+18 -5 3 H1 (1.370+/-0.021)e+23 -6 3 O16 (6.85+/-0.10)e+22 -7 3 B10 (6.85+/-0.10)e+18 + Cell Nuclide Atoms +0 1 H1 0.0+/-0 +1 1 O16 0.0+/-0 +2 1 B10 0.0+/-0 +3 1 U235 (3.482+/-0.008)e+23 +4 1 Mo99 (3.482+/-0.008)e+22 +5 2 H1 (1.400+/-0.021)e+23 +6 2 O16 (7.00+/-0.10)e+22 +7 2 B10 (7.00+/-0.10)e+18 +8 2 U235 0.0+/-0 +9 2 Mo99 0.0+/-0 +10 3 H1 (1.370+/-0.021)e+23 +11 3 O16 (6.85+/-0.10)e+22 +12 3 B10 (6.85+/-0.10)e+18 +13 3 U235 0.0+/-0 +14 3 Mo99 0.0+/-0 Volume calculation 1 Domain 1: 4.14+/-0.04 cm^3 Domain 2: 31.47+/-0.07 cm^3 @@ -18,8 +25,13 @@ Domain 2: 31.47+/-0.07 cm^3 0 1 H1 (2.770+/-0.029)e+23 1 1 O16 (1.385+/-0.014)e+23 2 1 B10 (1.385+/-0.014)e+19 -3 2 U235 (3.482+/-0.008)e+23 -4 2 Mo99 (3.482+/-0.008)e+22 +3 1 U235 0.0+/-0 +4 1 Mo99 0.0+/-0 +5 2 H1 0.0+/-0 +6 2 O16 0.0+/-0 +7 2 B10 0.0+/-0 +8 2 U235 (3.482+/-0.008)e+23 +9 2 Mo99 (3.482+/-0.008)e+22 Volume calculation 2 Domain 0: 35.61+/-0.07 cm^3 Universe Nuclide Atoms @@ -28,3 +40,23 @@ 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 +Domain 1: 31.35476+/-0.00010 cm^3 +Domain 2: 2.09147+/-0.00005 cm^3 +Domain 3: 2.110402+/-0.000033 cm^3 + Cell Nuclide Atoms +0 1 H1 0.0+/-0 +1 1 O16 0.0+/-0 +2 1 B10 0.0+/-0 +3 1 U235 (3.474110+/-0.000011)e+23 +4 1 Mo99 (3.474110+/-0.000011)e+22 +5 2 H1 (1.398833+/-0.000031)e+23 +6 2 O16 (6.99416+/-0.00015)e+22 +7 2 B10 (6.99416+/-0.00015)e+18 +8 2 U235 0.0+/-0 +9 2 Mo99 0.0+/-0 +10 3 H1 (1.401947+/-0.000022)e+23 +11 3 O16 (7.00974+/-0.00011)e+22 +12 3 B10 (7.00974+/-0.00011)e+18 +13 3 U235 0.0+/-0 +14 3 Mo99 0.0+/-0 diff --git a/tests/regression_tests/volume_calc/test.py b/tests/regression_tests/volume_calc/test.py index 2affc4085..6bcedd397 100644 --- a/tests/regression_tests/volume_calc/test.py +++ b/tests/regression_tests/volume_calc/test.py @@ -45,9 +45,12 @@ 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) ] + vol_calcs[-1].set_trigger(1e-04, 'std_dev') + # Define settings settings = openmc.Settings() settings.run_mode = 'volume' @@ -62,6 +65,11 @@ class VolumeTest(PyAPITestHarness): # Read volume calculation results volume_calc = openmc.VolumeCalculation.from_hdf5(filename) + if volume_calc.samples == 100: + assert(volume_calc.trigger_type == 'std_dev') + assert(volume_calc.threshold == 1e-04) + + # 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) @@ -74,4 +82,4 @@ class VolumeTest(PyAPITestHarness): def test_volume_calc(): harness = VolumeTest('') - harness._build_inputs() + harness.main() From 51ea427b25214a31033bad675069ee00f3e9d20c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 8 Oct 2019 09:44:24 -0500 Subject: [PATCH 15/32] Moved trigger logic into main func. --- include/openmc/volume_calc.h | 9 +- src/volume_calc.cpp | 329 +++++++++++++++++++---------------- 2 files changed, 186 insertions(+), 152 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index a7dd75009..982e6a9e7 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -13,9 +13,10 @@ namespace openmc { enum class ThresholdType { - VARIANCE = 0, - STD_DEV = 1, - REL_ERR = 2 + NONE = 0, + VARIANCE = 1, + STD_DEV = 2, + REL_ERR = 3 }; //============================================================================== @@ -115,7 +116,7 @@ private: // //! \param[in] seed_offset Seed offset used for independent calculations //! \return Vector of results for each user-specified domain - std::vector _execute(size_t seed_offset = 0) const; + std::vector _execute() const; }; diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index f65cb89ed..f2cb7a7d1 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -107,6 +107,7 @@ std::vector VolumeCalculation::execute() const { // execute the calculation once std::vector results = _execute(); + return results; // if no std. dev. threshold is set, return these resuls if (threshold_ == -1.0) { return results; } @@ -139,7 +140,7 @@ std::vector VolumeCalculation::execute() const { if (max_val <= threshold_) { break; } // perform the calculation - std::vector tmp = _execute(offset); + std::vector tmp = _execute(); offset += n_samples_; // update current results @@ -149,12 +150,13 @@ std::vector VolumeCalculation::execute() const { return results; } -std::vector VolumeCalculation::_execute(size_t seed_offset) const +std::vector VolumeCalculation::_execute() const { // Shared data that is collected from all threads 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 size_t min_samples = n_samples_ / mpi::n_procs; @@ -168,187 +170,218 @@ std::vector VolumeCalculation::_execute(size_t seed_o 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) { + #pragma omp parallel + { + // Variables that are private to each thread + std::vector> indices(n); + std::vector> hits(n); + Particle p; - prn_set_stream(STREAM_VOLUME); + prn_set_stream(STREAM_VOLUME); - // Sample locations and count hits - #pragma omp for - for (size_t i = i_start; i < i_end; i++) { - set_particle_seed(seed_offset + i); + // 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); - 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}; + 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 this location is not in the geometry at all, move on to next block - if (!find_cell(&p, false)) continue; + // 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; + 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); + double max_vol_err = -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]}; -#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); + // 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; + } + } + } + } + } else { + int q = master_indices[i_domain].size(); 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]; + 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); + } + #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.num_samples = total_samples; + + // update threshold value if needed + if (trigger_type_ != ThresholdType::NONE) { + double val = 0.0; + switch (trigger_type_) { + case ThresholdType::STD_DEV: + val = result.volume[1]; break; - } + case ThresholdType::REL_ERR: + val = result.volume[1] / result.volume[0]; + break; + case ThresholdType::VARIANCE: + val = result.volume[1] * result.volume[1]; + break; + } + // update max if entry is valid + if (val > 0.0) { max_vol_err = std::max(max_vol_err, 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); + } else { + result.nuclides.push_back(j); + result.atoms.push_back(0.0); + result.uncertainty.push_back(0.0); } } } - } 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); } -#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_; + // return results of the calculation + if (trigger_type_ == ThresholdType::NONE || max_vol_err < threshold_) { + return results; + } - 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_); - result.num_samples = 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); - } else { - result.nuclides.push_back(j); - result.atoms.push_back(0.0); - result.uncertainty.push_back(0.0); - } - } - } - } - - return results; + } // end while } void VolumeCalculation::to_hdf5(const std::string& filename, From 4c344e962617b193130ddc39260c6e2757ca67fd Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 8 Oct 2019 09:46:52 -0500 Subject: [PATCH 16/32] Removing unused code. --- include/openmc/volume_calc.h | 52 ++---------------------------------- src/volume_calc.cpp | 51 ++--------------------------------- 2 files changed, 4 insertions(+), 99 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 982e6a9e7..74331de87 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -31,48 +31,8 @@ public: std::array volume; //!< Mean/standard deviation of volume std::vector nuclides; //!< Index of nuclides std::vector atoms; //!< Number of atoms for each nuclide - std::vector uncertainty; //!< Uncertainty on number of atoms - size_t num_samples; - - Result& operator +=( const Result& other) { - Expects(volume.size() == other.volume.size()); - Expects(atoms.size() == atoms.size()); - - auto& a_samples = num_samples; - auto& b_samples = other.num_samples; - - size_t total_samples = num_samples + other.num_samples; - - for (int i = 0; i < volume.size(); i++) { - // calculate weighted average of volume results - auto& a_vol = volume[0]; - auto& b_vol = other.volume[0]; - volume[0] = (a_samples * a_vol + b_samples * b_vol) / total_samples; - - // propagate error - auto& a_err = volume[1]; - auto& b_err = other.volume[1]; - volume[1] = std::sqrt(a_samples * a_err * a_err + b_samples * b_err * b_err) / total_samples; - } - - for (int i = 0; i < atoms.size(); i++) { - // calculate weighted average of atom results - auto& a_atoms = atoms[i]; - auto& b_atoms = other.atoms[i]; - atoms[i] = (a_samples * a_atoms + b_samples * b_atoms) / total_samples; - - // propagate error - auto& a_err = uncertainty[i]; - auto& b_err = other.uncertainty[i]; - uncertainty[i] = std::sqrt(a_samples * a_err * a_err + b_samples * b_err * b_err) / total_samples; - } - - // update number of samples on the returned set of results; - num_samples = total_samples; - - return *this; - } - + std::vector uncertainty; //!< Uncertainty on number of atoms + int iterations; }; // Results for a single domain // Constructors @@ -111,14 +71,6 @@ private: void check_hit(int i_material, std::vector& indices, std::vector& hits) const; - //! \brief Perform calculation for domain volumes and average nuclide density - //! using n_samples_ - // - //! \param[in] seed_offset Seed offset used for independent calculations - //! \return Vector of results for each user-specified domain - std::vector _execute() const; - - }; //============================================================================== diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index f2cb7a7d1..21d9c3cb0 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -103,54 +103,7 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) } -std::vector VolumeCalculation::execute() const { - - // execute the calculation once - std::vector results = _execute(); - return results; - - // if no std. dev. threshold is set, return these resuls - if (threshold_ == -1.0) { return results; } - - size_t offset = n_samples_; - double max_val; - - while (true) { - // check maximum error value for all domains - max_val = -INFTY; - for (const auto& result : results) { - double val; - switch (trigger_type_) { - case ThresholdType::STD_DEV: - val = result.volume[1]; - break; - case ThresholdType::REL_ERR: - val = result.volume[1] / result.volume[0]; - break; - case ThresholdType::VARIANCE: - val = result.volume[1] * result.volume[1]; - break; - } - // update max if entry is valid - if (val > 0.0) { max_val = std::max(max_val, val); } - - } - - // exit once we're below our error limit - if (max_val <= threshold_) { break; } - - // perform the calculation - std::vector tmp = _execute(); - offset += n_samples_; - - // update current results - for (int i = 0; i < results.size(); i++) { results[i] += tmp[i]; } - } - - return results; -} - -std::vector VolumeCalculation::_execute() const +std::vector VolumeCalculation::execute() const { // Shared data that is collected from all threads int n = domain_ids_.size(); @@ -336,7 +289,7 @@ std::vector VolumeCalculation::_execute() const 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.num_samples = total_samples; + result.iterations = iterations; // update threshold value if needed if (trigger_type_ != ThresholdType::NONE) { From 24bfa0051a09ef06286b36b61537c3be921100ac Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 8 Oct 2019 10:55:04 -0500 Subject: [PATCH 17/32] Updating default trigger value and how trigger value is written to file --- include/openmc/volume_calc.h | 2 +- src/volume_calc.cpp | 18 +++---- .../volume_calc/results_true.dat | 54 ++++--------------- 3 files changed, 20 insertions(+), 54 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 74331de87..33d10209e 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -56,7 +56,7 @@ public: int 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 - ThresholdType trigger_type_; + ThresholdType trigger_type_ {ThresholdType::NONE}; 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/src/volume_calc.cpp b/src/volume_calc.cpp index 21d9c3cb0..23bb08f82 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -320,10 +320,6 @@ std::vector VolumeCalculation::execute() const result.nuclides.push_back(j); result.atoms.push_back(mean); result.uncertainty.push_back(stdev); - } else { - result.nuclides.push_back(j); - result.atoms.push_back(0.0); - result.uncertainty.push_back(0.0); } } } @@ -358,21 +354,23 @@ 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_); - if (threshold_ != -1.0) { + // Write trigger info + if (trigger_type_ != ThresholdType::NONE) { + write_attribute(file_id, "iterations", results[0].iterations); write_attribute(file_id, "threshold", threshold_); - + std::string trigger_str; switch(trigger_type_) { case ThresholdType::VARIANCE: - write_attribute(file_id, "trigger_type", "variance"); + trigger_str = "variance"; break; case ThresholdType::STD_DEV: - write_attribute(file_id, "trigger_type", "std_dev"); + trigger_str = "std_dev"; break; case ThresholdType::REL_ERR: - write_attribute(file_id, "trigger_type", "rel_err"); + trigger_str = "rel_err"; break; } - + write_attribute(file_id, "trigger_type", trigger_str); } if (domain_type_ == FILTER_CELL) { diff --git a/tests/regression_tests/volume_calc/results_true.dat b/tests/regression_tests/volume_calc/results_true.dat index 4fe3b4b6e..466139cd6 100644 --- a/tests/regression_tests/volume_calc/results_true.dat +++ b/tests/regression_tests/volume_calc/results_true.dat @@ -2,22 +2,15 @@ Volume calculation 0 Domain 1: 31.47+/-0.07 cm^3 Domain 2: 2.093+/-0.031 cm^3 Domain 3: 2.049+/-0.031 cm^3 - Cell Nuclide Atoms -0 1 H1 0.0+/-0 -1 1 O16 0.0+/-0 -2 1 B10 0.0+/-0 -3 1 U235 (3.482+/-0.008)e+23 -4 1 Mo99 (3.482+/-0.008)e+22 -5 2 H1 (1.400+/-0.021)e+23 -6 2 O16 (7.00+/-0.10)e+22 -7 2 B10 (7.00+/-0.10)e+18 -8 2 U235 0.0+/-0 -9 2 Mo99 0.0+/-0 -10 3 H1 (1.370+/-0.021)e+23 -11 3 O16 (6.85+/-0.10)e+22 -12 3 B10 (6.85+/-0.10)e+18 -13 3 U235 0.0+/-0 -14 3 Mo99 0.0+/-0 + Cell Nuclide Atoms +0 1 U235 (3.482+/-0.008)e+23 +1 1 Mo99 (3.482+/-0.008)e+22 +2 2 H1 (1.400+/-0.021)e+23 +3 2 O16 (7.00+/-0.10)e+22 +4 2 B10 (7.00+/-0.10)e+18 +5 3 H1 (1.370+/-0.021)e+23 +6 3 O16 (6.85+/-0.10)e+22 +7 3 B10 (6.85+/-0.10)e+18 Volume calculation 1 Domain 1: 4.14+/-0.04 cm^3 Domain 2: 31.47+/-0.07 cm^3 @@ -25,13 +18,8 @@ Domain 2: 31.47+/-0.07 cm^3 0 1 H1 (2.770+/-0.029)e+23 1 1 O16 (1.385+/-0.014)e+23 2 1 B10 (1.385+/-0.014)e+19 -3 1 U235 0.0+/-0 -4 1 Mo99 0.0+/-0 -5 2 H1 0.0+/-0 -6 2 O16 0.0+/-0 -7 2 B10 0.0+/-0 -8 2 U235 (3.482+/-0.008)e+23 -9 2 Mo99 (3.482+/-0.008)e+22 +3 2 U235 (3.482+/-0.008)e+23 +4 2 Mo99 (3.482+/-0.008)e+22 Volume calculation 2 Domain 0: 35.61+/-0.07 cm^3 Universe Nuclide Atoms @@ -40,23 +28,3 @@ 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 -Domain 1: 31.35476+/-0.00010 cm^3 -Domain 2: 2.09147+/-0.00005 cm^3 -Domain 3: 2.110402+/-0.000033 cm^3 - Cell Nuclide Atoms -0 1 H1 0.0+/-0 -1 1 O16 0.0+/-0 -2 1 B10 0.0+/-0 -3 1 U235 (3.474110+/-0.000011)e+23 -4 1 Mo99 (3.474110+/-0.000011)e+22 -5 2 H1 (1.398833+/-0.000031)e+23 -6 2 O16 (6.99416+/-0.00015)e+22 -7 2 B10 (6.99416+/-0.00015)e+18 -8 2 U235 0.0+/-0 -9 2 Mo99 0.0+/-0 -10 3 H1 (1.401947+/-0.000022)e+23 -11 3 O16 (7.00974+/-0.00011)e+22 -12 3 B10 (7.00974+/-0.00011)e+18 -13 3 U235 0.0+/-0 -14 3 Mo99 0.0+/-0 From 884a083a50eaf28ac86a3ba35ad96572534724c0 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 8 Oct 2019 10:55:45 -0500 Subject: [PATCH 18/32] Adding iteration property to volume calculation in Python API. --- openmc/volume.py | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/openmc/volume.py b/openmc/volume.py index 3996ce9dd..99eadb57d 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -58,7 +58,9 @@ class VolumeCalculation(object): volumes : dict Dictionary mapping unique IDs of domains to estimated volumes in cm^3. threshold : float - Threshold for the maxmimum standard deviation of volumes + Threshold for the maxmimum standard deviation of volumes. + iterations : int + Number of iterations over samples (for calculations with a trigger). trigger_type : {'variance', 'std_dev', 'rel_err'} Value type used to halt volume calculation @@ -68,6 +70,7 @@ class VolumeCalculation(object): self._volumes = {} self._threshold = None self._trigger_type = None + self._iterations = None cv.check_type('domains', domains, Iterable, (openmc.Cell, openmc.Material, openmc.Universe)) @@ -139,6 +142,10 @@ class VolumeCalculation(object): def trigger_type(self): return self._trigger_type + @property + def iterations(self): + return self._iterations + @property def domain_type(self): return self._domain_type @@ -189,10 +196,8 @@ class VolumeCalculation(object): @threshold.setter def threshold(self, threshold): name = 'volume std. dev. threshold' - cv.check_type(name, threshold, Real) - if (threshold <= 0.0): - raise ValueError("Invalid value '{}' (<= 0.0) provided for volume " - "calculation threshold.".format(threshold)) + cv.check_type(name, threshold, Real) + cv.check_greater_than(name, threshold, 0.0) self._threshold = threshold @trigger_type.setter @@ -201,6 +206,13 @@ class VolumeCalculation(object): ['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) @@ -247,15 +259,9 @@ class VolumeCalculation(object): lower_left = f.attrs['lower_left'] upper_right = f.attrs['upper_right'] - try: - threshold = f.attrs['threshold'] - except KeyError: - threshold = None - - try: - trigger_type = f.attrs['trigger_type'].decode() - except KeyError: - trigger_type = None + threshold = f.attrs.get('threshold') + trigger_type = f.attrs.get('trigger_type') + iterations = f.attrs.get('iterations', 1) volumes = {} atoms = {} @@ -288,9 +294,10 @@ class VolumeCalculation(object): # Instantiate the class and assign results vol = cls(domains, samples, lower_left, upper_right) - if threshold is not None: - vol.set_trigger(threshold, trigger_type) + if trigger_type is not None: + vol.set_trigger(threshold, trigger_type.decode()) + vol.iterations = iterations vol.volumes = volumes vol.atoms = atoms return vol From 17ebd34b452b94c8d9431bdae3a3d99af835e878 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 8 Oct 2019 10:55:54 -0500 Subject: [PATCH 19/32] Updating tests. --- .../volume_calc/inputs_true.dat | 18 ++++++- .../volume_calc/results_true.dat | 35 ++++++++++++++ tests/regression_tests/volume_calc/test.py | 48 +++++++++++++++++-- 3 files changed, 95 insertions(+), 6 deletions(-) diff --git a/tests/regression_tests/volume_calc/inputs_true.dat b/tests/regression_tests/volume_calc/inputs_true.dat index 1af345405..aaf6d8b00 100644 --- a/tests/regression_tests/volume_calc/inputs_true.dat +++ b/tests/regression_tests/volume_calc/inputs_true.dat @@ -54,6 +54,22 @@ 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..a31360abc 100644 --- a/tests/regression_tests/volume_calc/results_true.dat +++ b/tests/regression_tests/volume_calc/results_true.dat @@ -28,3 +28,38 @@ 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 +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 +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 +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 6bcedd397..8f003aa68 100644 --- a/tests/regression_tests/volume_calc/test.py +++ b/tests/regression_tests/volume_calc/test.py @@ -8,6 +8,19 @@ 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.std_dev_iters = 521 + + self.exp_rel_err = 1e-01 + self.rel_err_iters = 10 + + self.exp_variance = 5e-02 + self.variance_iters = 105 + def _build_inputs(self): # Define materials water = openmc.Material(1) @@ -46,10 +59,17 @@ class VolumeTest(PyAPITestHarness): openmc.VolumeCalculation(list(root.cells.values()), 100000), openmc.VolumeCalculation([water, fuel], 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[-1].set_trigger(1e-04, 'std_dev') + + 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() @@ -65,10 +85,28 @@ class VolumeTest(PyAPITestHarness): # Read volume calculation results volume_calc = openmc.VolumeCalculation.from_hdf5(filename) - if volume_calc.samples == 100: + if i == 3: assert(volume_calc.trigger_type == 'std_dev') - assert(volume_calc.threshold == 1e-04) - + assert(volume_calc.threshold == self.exp_std_dev) + assert(volume_calc.iterations == self.std_dev_iters) + for vol in volume_calc.volumes.values(): + assert(vol.std_dev <= self.exp_std_dev) + elif i == 4: + assert(volume_calc.trigger_type == 'rel_err') + assert(volume_calc.threshold == self.exp_rel_err) + assert(volume_calc.iterations == self.rel_err_iters) + for vol in volume_calc.volumes.values(): + assert(vol.std_dev/vol.nominal_value <= self.exp_rel_err) + elif i == 5: + assert(volume_calc.trigger_type == 'variance') + assert(volume_calc.threshold == self.exp_variance) + assert(volume_calc.iterations == self.variance_iters) + for vol in volume_calc.volumes.values(): + assert(vol.std_dev * vol.std_dev <= self.exp_variance) + else: + assert(volume_calc.trigger_type == None) + assert(volume_calc.threshold == None) + assert(volume_calc.iterations == 1) # Write cell volumes and total # of atoms for each nuclide for uid, volume in sorted(volume_calc.volumes.items()): @@ -82,4 +120,4 @@ class VolumeTest(PyAPITestHarness): def test_volume_calc(): harness = VolumeTest('') - harness.main() + harness.main() \ No newline at end of file From b2c064cc357415b4699302dabf7e0ea6d1f40d9d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 8 Oct 2019 11:34:26 -0500 Subject: [PATCH 20/32] Updates for MPI runs. --- src/volume_calc.cpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 23bb08f82..13cfcf569 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -325,11 +325,30 @@ std::vector VolumeCalculation::execute() const } } +#ifdef OPENMC_MPI + // update maximum error value on all processes + if (mpi::master) { + for (int i = 1; i < mpi::n_procs; i++) { + MPI_Send(&max_vol_err, 1, MPI_DOUBLE, i, 0, mpi::intracomm); + } + } else { + MPI_Recv(&max_vol_err, 1, MPI_DOUBLE, 0, 0, mpi::intracomm, MPI_STATUS_IGNORE); + } +#endif + // return results of the calculation if (trigger_type_ == ThresholdType::NONE || max_vol_err < threshold_) { return results; } +#ifdef OPENMC_MPI + // if iterating in MPI, need to zero indices and hits to they aren't counted twice + if (!mpi::master) { + for (auto& v : master_indices) { std::fill(v.begin(), v.end(), 0.0); } + for (auto& v : master_hits) { std::fill(v.begin(), v.end(), 0.0); } + } +#endif + } // end while } From eeb8b3cbb337635387de225c0d4c31819ac4e0df Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 8 Oct 2019 13:15:30 -0500 Subject: [PATCH 21/32] Updating name VolumeCalculation attribute. --- include/openmc/volume_calc.h | 4 ++-- src/volume_calc.cpp | 24 ++++++++++++------------ 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 33d10209e..dade67f5a 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -12,7 +12,7 @@ namespace openmc { -enum class ThresholdType { +enum class VolumeTriggerMetric { NONE = 0, VARIANCE = 1, STD_DEV = 2, @@ -56,7 +56,7 @@ public: int 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 - ThresholdType trigger_type_ {ThresholdType::NONE}; + VolumeTriggerMetric trigger_type_ {VolumeTriggerMetric::NONE}; //!< 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/src/volume_calc.cpp b/src/volume_calc.cpp index 13cfcf569..0b5c0204c 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -80,11 +80,11 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) std::string tmp = get_node_value(threshold_node, "type"); if (tmp == "variance") { - trigger_type_ = ThresholdType::VARIANCE; + trigger_type_ = VolumeTriggerMetric::VARIANCE; } else if (tmp == "std_dev") { - trigger_type_ = ThresholdType::STD_DEV; + trigger_type_ = VolumeTriggerMetric::STD_DEV; } else if ( tmp == "rel_err") { - trigger_type_ = ThresholdType::REL_ERR; + trigger_type_ = VolumeTriggerMetric::REL_ERR; } else { std::stringstream msg; msg << "Invalid volume calculation trigger type '" << tmp << "' provided."; @@ -292,16 +292,16 @@ std::vector VolumeCalculation::execute() const result.iterations = iterations; // update threshold value if needed - if (trigger_type_ != ThresholdType::NONE) { + if (trigger_type_ != VolumeTriggerMetric::NONE) { double val = 0.0; switch (trigger_type_) { - case ThresholdType::STD_DEV: + case VolumeTriggerMetric::STD_DEV: val = result.volume[1]; break; - case ThresholdType::REL_ERR: + case VolumeTriggerMetric::REL_ERR: val = result.volume[1] / result.volume[0]; break; - case ThresholdType::VARIANCE: + case VolumeTriggerMetric::VARIANCE: val = result.volume[1] * result.volume[1]; break; } @@ -337,7 +337,7 @@ std::vector VolumeCalculation::execute() const #endif // return results of the calculation - if (trigger_type_ == ThresholdType::NONE || max_vol_err < threshold_) { + if (trigger_type_ == VolumeTriggerMetric::NONE || max_vol_err < threshold_) { return results; } @@ -374,18 +374,18 @@ void VolumeCalculation::to_hdf5(const std::string& filename, write_attribute(file_id, "lower_left", lower_left_); write_attribute(file_id, "upper_right", upper_right_); // Write trigger info - if (trigger_type_ != ThresholdType::NONE) { + if (trigger_type_ != VolumeTriggerMetric::NONE) { write_attribute(file_id, "iterations", results[0].iterations); write_attribute(file_id, "threshold", threshold_); std::string trigger_str; switch(trigger_type_) { - case ThresholdType::VARIANCE: + case VolumeTriggerMetric::VARIANCE: trigger_str = "variance"; break; - case ThresholdType::STD_DEV: + case VolumeTriggerMetric::STD_DEV: trigger_str = "std_dev"; break; - case ThresholdType::REL_ERR: + case VolumeTriggerMetric::REL_ERR: trigger_str = "rel_err"; break; } From f251822b9f16af12f18f08ec7fa94c1fae1be59e Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 8 Oct 2019 13:31:37 -0500 Subject: [PATCH 22/32] Removing trailing whitespace in doc file. --- docs/source/usersguide/volume.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/usersguide/volume.rst b/docs/source/usersguide/volume.rst index d35386b38..b946858ea 100644 --- a/docs/source/usersguide/volume.rst +++ b/docs/source/usersguide/volume.rst @@ -42,8 +42,8 @@ or relative error of volume estimates using ::attr::`openmc.VolumeCalculation.se 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 +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. From 919c19c9094add032881fe07ad7d12d1091d6424 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 8 Oct 2019 13:53:22 -0500 Subject: [PATCH 23/32] Cleanup from self-review. --- include/openmc/volume_calc.h | 4 +- openmc/volume.py | 13 +++---- src/volume_calc.cpp | 43 ++++++++++++---------- tests/regression_tests/volume_calc/test.py | 6 +-- 4 files changed, 35 insertions(+), 31 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index dade67f5a..bf85a435d 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -31,8 +31,8 @@ public: std::array volume; //!< Mean/standard deviation of volume 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; + std::vector uncertainty; //!< Uncertainty on number of atoms + int iterations; //!< Number of iterations needed to obtain the results }; // Results for a single domain // Constructors diff --git a/openmc/volume.py b/openmc/volume.py index 99eadb57d..a59fb2acf 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -32,7 +32,6 @@ class VolumeCalculation(object): Upper-right coordinates of bounding box used to sample points. If this argument is not supplied, an attempt is made to automatically determine a bounding box. - Attributes ---------- @@ -59,10 +58,10 @@ class VolumeCalculation(object): Dictionary mapping unique IDs of domains to estimated volumes in cm^3. threshold : float Threshold for the maxmimum standard deviation of volumes. - iterations : int - Number of iterations over samples (for calculations with a trigger). 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): @@ -83,7 +82,7 @@ class VolumeCalculation(object): self.ids = [d.id for d in domains] self.samples = samples - + if lower_left is not None: if upper_right is None: raise ValueError('Both lower-left and upper-right coordinates ' @@ -193,7 +192,7 @@ class VolumeCalculation(object): cv.check_length(name, upper_right, 3) self._upper_right = upper_right - @threshold.setter + @threshold.setter def threshold(self, threshold): name = 'volume std. dev. threshold' cv.check_type(name, threshold, Real) @@ -203,7 +202,7 @@ class VolumeCalculation(object): @trigger_type.setter def trigger_type(self, trigger_type): cv.check_value('tally trigger type', trigger_type, - ['variance', 'std_dev', 'rel_err']) + ('variance', 'std_dev', 'rel_err')) self._trigger_type = trigger_type @iterations.setter @@ -293,7 +292,7 @@ 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()) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 0b5c0204c..32b6f03e2 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -64,13 +64,13 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) if (size_t_stream.fail()) { std::stringstream msg; msg << "Could not read number of samples (" - << size_t_stream.str() << ")\n"; + << size_t_stream.str() << ")"; fatal_error(msg); } 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; @@ -124,6 +124,7 @@ std::vector VolumeCalculation::execute() const } while (true) { + #pragma omp parallel { // Variables that are private to each thread @@ -223,7 +224,8 @@ std::vector VolumeCalculation::execute() const iterations++; size_t total_samples = iterations * n_samples_; - double max_vol_err = -INFTY; + // reset + double trigger_val = -INFTY; // Set size for members of the Result struct std::vector results(n); @@ -237,7 +239,7 @@ std::vector VolumeCalculation::execute() const auto n_nuc = data::nuclides.size(); xt::xtensor atoms({n_nuc, 2}, 0.0); - #ifdef OPENMC_MPI +#ifdef OPENMC_MPI if (mpi::master) { for (int j = 1; j < mpi::n_procs; j++) { int q; @@ -264,7 +266,7 @@ std::vector VolumeCalculation::execute() const MPI_Send(&q, 1, MPI_INTEGER, 0, 0, mpi::intracomm); MPI_Send(&buffer[0], 2*q, MPI_INTEGER, 0, 1, mpi::intracomm); } - #endif +#endif if (mpi::master) { int total_hits = 0; @@ -299,14 +301,14 @@ std::vector VolumeCalculation::execute() const val = result.volume[1]; break; case VolumeTriggerMetric::REL_ERR: - val = result.volume[1] / result.volume[0]; + val = result.volume[0] == 0.0 ? 0.0 : result.volume[1] / result.volume[0]; break; case VolumeTriggerMetric::VARIANCE: val = result.volume[1] * result.volume[1]; break; } // update max if entry is valid - if (val > 0.0) { max_vol_err = std::max(max_vol_err, val); } + if (val > 0.0) { trigger_val = std::max(trigger_val, val); } } for (int j = 0; j < n_nuc; ++j) { @@ -323,29 +325,30 @@ std::vector VolumeCalculation::execute() const } } } - } + } // end domain loop + + // if no trigger is applied, we're done + if (trigger_type_ == VolumeTriggerMetric::NONE) { return results; } #ifdef OPENMC_MPI // update maximum error value on all processes if (mpi::master) { for (int i = 1; i < mpi::n_procs; i++) { - MPI_Send(&max_vol_err, 1, MPI_DOUBLE, i, 0, mpi::intracomm); + MPI_Send(&trigger_val, 1, MPI_DOUBLE, i, 0, mpi::intracomm); } } else { - MPI_Recv(&max_vol_err, 1, MPI_DOUBLE, 0, 0, mpi::intracomm, MPI_STATUS_IGNORE); - } + MPI_Recv(&trigger_val, 1, MPI_DOUBLE, 0, 0, mpi::intracomm, MPI_STATUS_IGNORE); + } #endif - + // return results of the calculation - if (trigger_type_ == VolumeTriggerMetric::NONE || max_vol_err < threshold_) { - return results; - } + if (trigger_val < threshold_) { return results; } #ifdef OPENMC_MPI - // if iterating in MPI, need to zero indices and hits to they aren't counted twice + // 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.0); } - for (auto& v : master_hits) { std::fill(v.begin(), v.end(), 0.0); } + 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 @@ -373,7 +376,7 @@ 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 + // Write trigger info if (trigger_type_ != VolumeTriggerMetric::NONE) { write_attribute(file_id, "iterations", results[0].iterations); write_attribute(file_id, "threshold", threshold_); @@ -390,6 +393,8 @@ void VolumeCalculation::to_hdf5(const std::string& filename, break; } write_attribute(file_id, "trigger_type", trigger_str); + } else { + write_attribute(file_id, "iterations", 1); } if (domain_type_ == FILTER_CELL) { diff --git a/tests/regression_tests/volume_calc/test.py b/tests/regression_tests/volume_calc/test.py index 8f003aa68..96acf59b9 100644 --- a/tests/regression_tests/volume_calc/test.py +++ b/tests/regression_tests/volume_calc/test.py @@ -64,7 +64,7 @@ class VolumeTest(PyAPITestHarness): 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') @@ -102,7 +102,7 @@ class VolumeTest(PyAPITestHarness): assert(volume_calc.threshold == self.exp_variance) assert(volume_calc.iterations == self.variance_iters) for vol in volume_calc.volumes.values(): - assert(vol.std_dev * vol.std_dev <= self.exp_variance) + assert(vol.std_dev * vol.std_dev <= self.exp_variance) else: assert(volume_calc.trigger_type == None) assert(volume_calc.threshold == None) @@ -120,4 +120,4 @@ class VolumeTest(PyAPITestHarness): def test_volume_calc(): harness = VolumeTest('') - harness.main() \ No newline at end of file + harness.main() From a4f98661c1bab079323f58ebc0fa6aa54906b703 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 14 Oct 2019 18:45:30 -0500 Subject: [PATCH 24/32] Update docs/source/usersguide/volume.rst Co-Authored-By: Paul Romano --- docs/source/usersguide/volume.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/usersguide/volume.rst b/docs/source/usersguide/volume.rst index b946858ea..c040eb2d7 100644 --- a/docs/source/usersguide/volume.rst +++ b/docs/source/usersguide/volume.rst @@ -38,7 +38,7 @@ 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 ::attr::`openmc.VolumeCalculation.set_trigger`:: +or relative error of volume estimates using :meth:`openmc.VolumeCalculation.set_trigger`:: vol_calc.set_trigger(1e-05, 'std_dev') @@ -76,4 +76,4 @@ After the results are loaded, volume estimates will be stored in :attr:`VolumeCalculation.volumes`. There is also a :attr:`VolumeCalculation.atoms_dataframe` attribute that shows stochastic estimates of the number of atoms of each type of nuclide within the specified -domains along with their uncertainties. \ No newline at end of file +domains along with their uncertainties. From 143fe6405757b15f34b1396ee63e3be2dfb668a6 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 14 Oct 2019 20:26:23 -0500 Subject: [PATCH 25/32] Using trigger enum --- include/openmc/volume_calc.h | 10 ++-------- src/volume_calc.cpp | 24 ++++++++++++------------ 2 files changed, 14 insertions(+), 20 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index bf85a435d..1b49da814 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -2,6 +2,7 @@ #define OPENMC_VOLUME_CALC_H #include "openmc/position.h" +#include "openmc/tallies/trigger.h" #include "pugixml.hpp" #include "xtensor/xtensor.hpp" @@ -12,13 +13,6 @@ namespace openmc { -enum class VolumeTriggerMetric { - NONE = 0, - VARIANCE = 1, - STD_DEV = 2, - REL_ERR = 3 -}; - //============================================================================== // Volume calculation class //============================================================================== @@ -56,7 +50,7 @@ public: int 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 - VolumeTriggerMetric trigger_type_ {VolumeTriggerMetric::NONE}; //!< Trigger metric for the volume calculation + 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/src/volume_calc.cpp b/src/volume_calc.cpp index 32b6f03e2..b354b1db6 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -80,11 +80,11 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) std::string tmp = get_node_value(threshold_node, "type"); if (tmp == "variance") { - trigger_type_ = VolumeTriggerMetric::VARIANCE; + trigger_type_ = TriggerMetric::variance; } else if (tmp == "std_dev") { - trigger_type_ = VolumeTriggerMetric::STD_DEV; + trigger_type_ = TriggerMetric::standard_deviation; } else if ( tmp == "rel_err") { - trigger_type_ = VolumeTriggerMetric::REL_ERR; + trigger_type_ = TriggerMetric::relative_error; } else { std::stringstream msg; msg << "Invalid volume calculation trigger type '" << tmp << "' provided."; @@ -294,16 +294,16 @@ std::vector VolumeCalculation::execute() const result.iterations = iterations; // update threshold value if needed - if (trigger_type_ != VolumeTriggerMetric::NONE) { + if (trigger_type_ != TriggerMetric::not_active) { double val = 0.0; switch (trigger_type_) { - case VolumeTriggerMetric::STD_DEV: + case TriggerMetric::standard_deviation: val = result.volume[1]; break; - case VolumeTriggerMetric::REL_ERR: + case TriggerMetric::relative_error: val = result.volume[0] == 0.0 ? 0.0 : result.volume[1] / result.volume[0]; break; - case VolumeTriggerMetric::VARIANCE: + case TriggerMetric::variance: val = result.volume[1] * result.volume[1]; break; } @@ -328,7 +328,7 @@ std::vector VolumeCalculation::execute() const } // end domain loop // if no trigger is applied, we're done - if (trigger_type_ == VolumeTriggerMetric::NONE) { return results; } + if (trigger_type_ == TriggerMetric::not_active) { return results; } #ifdef OPENMC_MPI // update maximum error value on all processes @@ -377,18 +377,18 @@ void VolumeCalculation::to_hdf5(const std::string& filename, write_attribute(file_id, "lower_left", lower_left_); write_attribute(file_id, "upper_right", upper_right_); // Write trigger info - if (trigger_type_ != VolumeTriggerMetric::NONE) { + 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 VolumeTriggerMetric::VARIANCE: + case TriggerMetric::variance: trigger_str = "variance"; break; - case VolumeTriggerMetric::STD_DEV: + case TriggerMetric::standard_deviation: trigger_str = "std_dev"; break; - case VolumeTriggerMetric::REL_ERR: + case TriggerMetric::relative_error: trigger_str = "rel_err"; break; } From 5c723bd360f8b8a73ac7c39b9a2887bdb353ad4f Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 14 Oct 2019 20:27:18 -0500 Subject: [PATCH 26/32] Keeping Python API checks, but writing number of iterations to results file. --- .../volume_calc/results_true.dat | 18 +++++++++++++ tests/regression_tests/volume_calc/test.py | 27 +++++++++---------- 2 files changed, 31 insertions(+), 14 deletions(-) diff --git a/tests/regression_tests/volume_calc/results_true.dat b/tests/regression_tests/volume_calc/results_true.dat index a31360abc..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 @@ -29,6 +38,9 @@ Domain 0: 35.61+/-0.07 cm^3 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 @@ -42,6 +54,9 @@ Domain 3: 2.11+/-0.04 cm^3 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 @@ -51,6 +66,9 @@ Domain 2: 30.5+/-0.7 cm^3 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 diff --git a/tests/regression_tests/volume_calc/test.py b/tests/regression_tests/volume_calc/test.py index 96acf59b9..ae4987f1c 100644 --- a/tests/regression_tests/volume_calc/test.py +++ b/tests/regression_tests/volume_calc/test.py @@ -13,13 +13,8 @@ class VolumeTest(PyAPITestHarness): super().__init__(*args, **kwargs) self.exp_std_dev = 1e-01 - self.std_dev_iters = 521 - self.exp_rel_err = 1e-01 - self.rel_err_iters = 10 - self.exp_variance = 5e-02 - self.variance_iters = 105 def _build_inputs(self): # Define materials @@ -85,29 +80,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) - assert(volume_calc.iterations == self.std_dev_iters) - for vol in volume_calc.volumes.values(): - assert(vol.std_dev <= self.exp_std_dev) elif i == 4: assert(volume_calc.trigger_type == 'rel_err') assert(volume_calc.threshold == self.exp_rel_err) - assert(volume_calc.iterations == self.rel_err_iters) - for vol in volume_calc.volumes.values(): - assert(vol.std_dev/vol.nominal_value <= self.exp_rel_err) elif i == 5: assert(volume_calc.trigger_type == 'variance') assert(volume_calc.threshold == self.exp_variance) - assert(volume_calc.iterations == self.variance_iters) - for vol in volume_calc.volumes.values(): - assert(vol.std_dev * vol.std_dev <= 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) From e91a21add239108c9f950125661b0d6fc626eb8d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 14 Oct 2019 20:30:25 -0500 Subject: [PATCH 27/32] Simplyfing MPI send to broadcast. --- src/volume_calc.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index b354b1db6..a89c67368 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -333,9 +333,7 @@ std::vector VolumeCalculation::execute() const #ifdef OPENMC_MPI // update maximum error value on all processes if (mpi::master) { - for (int i = 1; i < mpi::n_procs; i++) { - MPI_Send(&trigger_val, 1, MPI_DOUBLE, i, 0, mpi::intracomm); - } + MPI_Bcast(&trigger_val, 1, MPI_DOUBLE, mpi:intracomm); } else { MPI_Recv(&trigger_val, 1, MPI_DOUBLE, 0, 0, mpi::intracomm, MPI_STATUS_IGNORE); } From 6c387b66c56620bf5fa506d7b325e0297368370b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 14 Oct 2019 20:39:52 -0500 Subject: [PATCH 28/32] Setting volume variance to INF if the domain is unsampled to avoid early trigger exit. --- src/volume_calc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index a89c67368..2c4aa8287 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -301,7 +301,7 @@ std::vector VolumeCalculation::execute() const val = result.volume[1]; break; case TriggerMetric::relative_error: - val = result.volume[0] == 0.0 ? 0.0 : result.volume[1] / result.volume[0]; + val = result.volume[0] == 0.0 ? INFTY : result.volume[1] / result.volume[0]; break; case TriggerMetric::variance: val = result.volume[1] * result.volume[1]; From f47e35da75e6c6a8522e1234a4f57669680777ea Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 14 Oct 2019 20:53:04 -0500 Subject: [PATCH 29/32] Removing blank line. --- tests/regression_tests/volume_calc/test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/regression_tests/volume_calc/test.py b/tests/regression_tests/volume_calc/test.py index ae4987f1c..f7b9a27a0 100644 --- a/tests/regression_tests/volume_calc/test.py +++ b/tests/regression_tests/volume_calc/test.py @@ -59,7 +59,6 @@ class VolumeTest(PyAPITestHarness): 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') From d213373a4ba18d4bea3a13b8d576b823d758d9d9 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 14 Oct 2019 20:55:52 -0500 Subject: [PATCH 30/32] Converting to unsigned long long for number of samples in volume calc. --- src/volume_calc.cpp | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 2c4aa8287..aa685e353 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -59,14 +59,7 @@ 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"); - std::stringstream size_t_stream(get_node_value(node, "samples")); - size_t_stream >> n_samples_; - if (size_t_stream.fail()) { - std::stringstream msg; - msg << "Could not read number of samples (" - << size_t_stream.str() << ")"; - fatal_error(msg); - } + n_samples_ = std::stoull(get_node_value(node, "samples")); if (check_for_node(node, "threshold")) { pugi::xml_node threshold_node = node.child("threshold"); From dc5fe42af3deda3f43263203e9e7aca9fc158bd1 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 15 Oct 2019 00:40:58 -0500 Subject: [PATCH 31/32] Update src/volume_calc.cpp Co-Authored-By: Paul Romano --- src/volume_calc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index aa685e353..0c91c463e 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -326,7 +326,7 @@ std::vector VolumeCalculation::execute() const #ifdef OPENMC_MPI // update maximum error value on all processes if (mpi::master) { - MPI_Bcast(&trigger_val, 1, MPI_DOUBLE, mpi:intracomm); + MPI_Bcast(&trigger_val, 1, MPI_DOUBLE, 0, mpi::intracomm); } else { MPI_Recv(&trigger_val, 1, MPI_DOUBLE, 0, 0, mpi::intracomm, MPI_STATUS_IGNORE); } From 67907e53d9279eefb9a138e82e449bcb965550c9 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 15 Oct 2019 00:43:21 -0500 Subject: [PATCH 32/32] Removing if block in MPI ifdef. --- src/volume_calc.cpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 0c91c463e..419bf0ea5 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -325,11 +325,7 @@ std::vector VolumeCalculation::execute() const #ifdef OPENMC_MPI // update maximum error value on all processes - if (mpi::master) { - MPI_Bcast(&trigger_val, 1, MPI_DOUBLE, 0, mpi::intracomm); - } else { - MPI_Recv(&trigger_val, 1, MPI_DOUBLE, 0, 0, mpi::intracomm, MPI_STATUS_IGNORE); - } + MPI_Bcast(&trigger_val, 1, MPI_DOUBLE, 0, mpi::intracomm); #endif // return results of the calculation