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()