From 382bcb2e8ea63b5ffd9ead818c651e5aa9d4c7ce Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 30 Jun 2023 22:29:47 -0500 Subject: [PATCH] Updates to object lifecycle for WeightWindows and WeightWindowGenerators (#2582) --- include/openmc/weight_windows.h | 318 +++++++++--------- openmc/weight_windows.py | 13 +- src/initialize.cpp | 5 + src/weight_windows.cpp | 121 ++++--- .../weightwindows/generators/test.py | 13 +- tests/unit_tests/weightwindows/test.py | 16 +- 6 files changed, 275 insertions(+), 211 deletions(-) diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 75073e917..992c8fdc1 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -15,211 +15,215 @@ #include "openmc/tallies/tally.h" #include "openmc/vector.h" - namespace openmc -{ +namespace openmc { - enum class WeightWindowUpdateMethod { - MAGIC, - }; +enum class WeightWindowUpdateMethod { + MAGIC, +}; - //============================================================================== - // Constants - //============================================================================== +//============================================================================== +// Constants +//============================================================================== - constexpr double DEFAULT_WEIGHT_CUTOFF {1.0e-38}; // default low weight cutoff +constexpr double DEFAULT_WEIGHT_CUTOFF {1.0e-38}; // default low weight cutoff - //============================================================================== - // Non-member functions - //============================================================================== +//============================================================================== +// Non-member functions +//============================================================================== - //! Apply weight windows to a particle - //! \param[in] p Particle to apply weight windows to - void apply_weight_windows(Particle & p); +//! Apply weight windows to a particle +//! \param[in] p Particle to apply weight windows to +void apply_weight_windows(Particle& p); - //! Free memory associated with weight windows - void free_memory_weight_windows(); +//! Free memory associated with weight windows +void free_memory_weight_windows(); - //============================================================================== - // Global variables - //============================================================================== +//============================================================================== +// Global variables +//============================================================================== - class WeightWindows; - class WeightWindowsGenerator; +class WeightWindows; +class WeightWindowsGenerator; - namespace variance_reduction { +namespace variance_reduction { - extern std::unordered_map ww_map; - extern vector> weight_windows; - extern vector> weight_windows_generators; +extern std::unordered_map ww_map; +extern vector> weight_windows; +extern vector> weight_windows_generators; - } // namespace variance_reduction +} // namespace variance_reduction - //============================================================================== - //! Individual weight window information - //============================================================================== +//============================================================================== +//! Individual weight window information +//============================================================================== - struct WeightWindow { - double lower_weight {-1}; // -1 indicates invalid state - double upper_weight {1}; - double max_lb_ratio {1}; - double survival_weight {0.5}; - double weight_cutoff {DEFAULT_WEIGHT_CUTOFF}; - int max_split {1}; +struct WeightWindow { + double lower_weight {-1}; // -1 indicates invalid state + double upper_weight {1}; + double max_lb_ratio {1}; + double survival_weight {0.5}; + double weight_cutoff {DEFAULT_WEIGHT_CUTOFF}; + int max_split {10}; - //! Whether the weight window is in a valid state - bool is_valid() const { return lower_weight >= 0.0; } + //! Whether the weight window is in a valid state + bool is_valid() const { return lower_weight >= 0.0; } - //! Adjust the weight window by a constant factor - void scale(double factor) - { - lower_weight *= factor; - upper_weight *= factor; - } - }; + //! Adjust the weight window by a constant factor + void scale(double factor) + { + lower_weight *= factor; + upper_weight *= factor; + } +}; - //============================================================================== - //! Weight window settings - //============================================================================== +//============================================================================== +//! Weight window settings +//============================================================================== - class WeightWindows { - public: - //---------------------------------------------------------------------------- - // Constructors - WeightWindows(int32_t id = -1); - WeightWindows(pugi::xml_node node); - ~WeightWindows(); - static WeightWindows* create(int32_t id = -1); - static WeightWindows* from_hdf5( - hid_t wws_group, const std::string& group_name); +class WeightWindows { +public: + //---------------------------------------------------------------------------- + // Constructors + WeightWindows(int32_t id = -1); + WeightWindows(pugi::xml_node node); + ~WeightWindows(); + static WeightWindows* create(int32_t id = -1); + static WeightWindows* from_hdf5( + hid_t wws_group, const std::string& group_name); - //---------------------------------------------------------------------------- - // Methods - private: - template - void check_bounds(const T& lower, const T& upper) const; + //---------------------------------------------------------------------------- + // Methods +private: + template + void check_bounds(const T& lower, const T& upper) const; - template - void check_bounds(const T& lower) const; + template + void check_bounds(const T& lower) const; - void check_tally_update_compatibility(const Tally* tally); + void check_tally_update_compatibility(const Tally* tally); - public: - //! Set the weight window ID - void set_id(int32_t id = -1); +public: + //! Set the weight window ID + void set_id(int32_t id = -1); - void set_energy_bounds(gsl::span bounds); + void set_energy_bounds(gsl::span bounds); - void set_mesh(const std::unique_ptr& mesh); + void set_mesh(const std::unique_ptr& mesh); - void set_mesh(const Mesh* mesh); + void set_mesh(const Mesh* mesh); - void set_mesh(int32_t mesh_idx); + void set_mesh(int32_t mesh_idx); - //! Ready the weight window class for use - void set_defaults(); + //! Ready the weight window class for use + void set_defaults(); - //! Update weight window boundaries using tally results - //! \param[in] tally Pointer to the tally whose results will be used to - //! update weight windows \param[in] value String representing the type of - //! value to use for weight window generation (one of "mean" or "rel_err") - //! \param[in] threshold Relative error threshold. Results over this - //! threshold will be ignored \param[in] ratio Ratio of upper to lower - //! weight window bounds - void update_magic(const Tally* tally, const std::string& value = "mean", - double threshold = 1.0, double ratio = 5.0); + //! Ensure the weight window lower bounds are properly allocated + void allocate_ww_bounds(); - // NOTE: This is unused for now but may be used in the future - //! Write weight window settings to an HDF5 file - //! \param[in] group HDF5 group to write to - void to_hdf5(hid_t group) const; + //! Update weight window boundaries using tally results + //! \param[in] tally Pointer to the tally whose results will be used to + //! update weight windows \param[in] value String representing the type of + //! value to use for weight window generation (one of "mean" or "rel_err") + //! \param[in] threshold Relative error threshold. Results over this + //! threshold will be ignored \param[in] ratio Ratio of upper to lower + //! weight window bounds + void update_magic(const Tally* tally, const std::string& value = "mean", + double threshold = 1.0, double ratio = 5.0); - //! Retrieve the weight window for a particle - //! \param[in] p Particle to get weight window for - WeightWindow get_weight_window(const Particle& p) const; + // NOTE: This is unused for now but may be used in the future + //! Write weight window settings to an HDF5 file + //! \param[in] group HDF5 group to write to + void to_hdf5(hid_t group) const; - std::array bounds_size() const; + //! Retrieve the weight window for a particle + //! \param[in] p Particle to get weight window for + WeightWindow get_weight_window(const Particle& p) const; - const vector& energy_bounds() const { return energy_bounds_; } + std::array bounds_size() const; - void set_bounds(const xt::xtensor& lower_ww_bounds, - const xt::xtensor& upper_bounds); + const vector& energy_bounds() const { return energy_bounds_; } - void set_bounds(const xt::xtensor& lower_bounds, double ratio); + void set_bounds(const xt::xtensor& lower_ww_bounds, + const xt::xtensor& upper_bounds); - void set_bounds(gsl::span lower_bounds, - gsl::span upper_bounds); + void set_bounds(const xt::xtensor& lower_bounds, double ratio); - void set_bounds(gsl::span lower_bounds, double ratio); + void set_bounds( + gsl::span lower_bounds, gsl::span upper_bounds); - void set_particle_type(ParticleType p_type); + void set_bounds(gsl::span lower_bounds, double ratio); - //---------------------------------------------------------------------------- - // Accessors - int32_t id() const { return id_; } - int32_t& id() { return id_; } + void set_particle_type(ParticleType p_type); - int32_t index() const { return index_; } + //---------------------------------------------------------------------------- + // Accessors + int32_t id() const { return id_; } + int32_t& id() { return id_; } - vector& energy_bounds() { return energy_bounds_; } + int32_t index() const { return index_; } - const std::unique_ptr& mesh() const - { - return model::meshes[mesh_idx_]; - } + vector& energy_bounds() { return energy_bounds_; } - const xt::xtensor& lower_ww_bounds() const { return lower_ww_; } - xt::xtensor& lower_ww_bounds() { return lower_ww_; } + const std::unique_ptr& mesh() const { return model::meshes[mesh_idx_]; } - const xt::xtensor& upper_ww_bounds() const { return upper_ww_; } - xt::xtensor& upper_ww_bounds() { return upper_ww_; } + const xt::xtensor& lower_ww_bounds() const { return lower_ww_; } + xt::xtensor& lower_ww_bounds() { return lower_ww_; } - ParticleType particle_type() const { return particle_type_; } + const xt::xtensor& upper_ww_bounds() const { return upper_ww_; } + xt::xtensor& upper_ww_bounds() { return upper_ww_; } - private: - //---------------------------------------------------------------------------- - // Data members - int32_t id_; //!< Unique ID - gsl::index index_; //!< Index into weight windows vector - ParticleType particle_type_ { - ParticleType::neutron}; //!< Particle type to apply weight windows to - vector energy_bounds_; //!< Energy boundaries [eV] - xt::xtensor lower_ww_; //!< Lower weight window bounds (shape: - //!< energy_bins, mesh_bins (k, j, i)) - xt::xtensor - upper_ww_; //!< Upper weight window bounds (shape: energy_bins, mesh_bins) - double survival_ratio_ {3.0}; //!< Survival weight ratio - double max_lb_ratio_ { - 1.0}; //!< Maximum lower bound to particle weight ratio - double weight_cutoff_ {DEFAULT_WEIGHT_CUTOFF}; //!< Weight cutoff - int max_split_ {10}; //!< Maximum value for particle splitting - int32_t mesh_idx_; //!< Index in meshes vector - }; + ParticleType particle_type() const { return particle_type_; } - class WeightWindowsGenerator { - public: - // Constructors - WeightWindowsGenerator(pugi::xml_node node); +private: + //---------------------------------------------------------------------------- + // Data members + int32_t id_; //!< Unique ID + gsl::index index_; //!< Index into weight windows vector + ParticleType particle_type_ { + ParticleType::neutron}; //!< Particle type to apply weight windows to + vector energy_bounds_; //!< Energy boundaries [eV] + xt::xtensor lower_ww_; //!< Lower weight window bounds (shape: + //!< energy_bins, mesh_bins (k, j, i)) + xt::xtensor + upper_ww_; //!< Upper weight window bounds (shape: energy_bins, mesh_bins) + double survival_ratio_ {3.0}; //!< Survival weight ratio + double max_lb_ratio_ {1.0}; //!< Maximum lower bound to particle weight ratio + double weight_cutoff_ {DEFAULT_WEIGHT_CUTOFF}; //!< Weight cutoff + int max_split_ {10}; //!< Maximum value for particle splitting + int32_t mesh_idx_ {-1}; //!< Index in meshes vector +}; - // Methods - void update() const; +class WeightWindowsGenerator { +public: + // Constructors + WeightWindowsGenerator(pugi::xml_node node); - // Data members - int32_t - tally_idx_; //!< Index of the tally used to update the weight windows - int32_t ww_idx_; //!< Index of the weight windows object being generated - std::string method_; //!< Method used to update weight window. Only "magic" - //!< is valid for now. - int32_t max_realizations_; //!< Maximum number of tally realizations - int32_t update_interval_; //!< Determines how often updates occur - bool on_the_fly_; //!< Whether or not weight windows + // Methods + void update() const; - // MAGIC update parameters - std::string tally_value_ { - "mean"}; // int: if self.energy_bounds is None: - raise ValueError('Energy bounds are not set') + return 1 return self.energy_bounds.size - 1 @property @@ -650,7 +653,8 @@ class WeightWindowGenerator: Mesh used to represent the weight windows spatially energy_bounds : Iterable of Real A list of values for which each successive pair constitutes a range of - energies in [eV] for a single bin + energies in [eV] for a single bin. If no energy bins are provided, the + maximum and minimum energy for the data available at runtime. particle_type : {'neutron', 'photon'} Particle type the weight windows apply to @@ -681,6 +685,7 @@ class WeightWindowGenerator: def __init__(self, mesh, energy_bounds=None, particle_type='neutron'): self.mesh = mesh + self._energy_bounds = None if energy_bounds is not None: self.energy_bounds = energy_bounds self.particle_type = particle_type diff --git a/src/initialize.cpp b/src/initialize.cpp index 323da2cc0..d667f4b03 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -33,6 +33,7 @@ #include "openmc/thermal.h" #include "openmc/timer.h" #include "openmc/vector.h" +#include "openmc/weight_windows.h" #ifdef LIBMESH #include "libmesh/libmesh.h" @@ -404,6 +405,8 @@ bool read_model_xml() } } + finalize_variance_reduction(); + return true; } @@ -427,6 +430,8 @@ void read_separate_xml_files() // Read the plots.xml regardless of plot mode in case plots are requested // via the API read_plots_xml(); + + finalize_variance_reduction(); } void initial_output() diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index dcdf13abe..a1f8731f8 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -257,21 +257,26 @@ WeightWindows* WeightWindows::from_hdf5( void WeightWindows::set_defaults() { - // ensure default values are set + // set energy bounds to the min/max energy supported by the data if (energy_bounds_.size() == 0) { int p_type = static_cast(particle_type_); energy_bounds_.push_back(data::energy_min[p_type]); energy_bounds_.push_back(data::energy_max[p_type]); } +} - // some constructors won't allocate space for the bounds - // do that here so the object is valid - if (lower_ww_.size() == 0 || upper_ww_.size() == 0) { - lower_ww_ = xt::empty(bounds_size()); - lower_ww_.fill(-1); - upper_ww_ = xt::empty(bounds_size()); - upper_ww_.fill(-1); +void WeightWindows::allocate_ww_bounds() +{ + auto shape = bounds_size(); + if (shape[0] * shape[1] == 0) { + auto msg = fmt::format( + "Size of weight window bounds is zero for WeightWindows {}", id()); + warning(msg); } + lower_ww_ = xt::empty(shape); + lower_ww_.fill(-1); + upper_ww_ = xt::empty(shape); + upper_ww_.fill(-1); } void WeightWindows::set_id(int32_t id) @@ -308,6 +313,9 @@ void WeightWindows::set_energy_bounds(gsl::span bounds) { energy_bounds_.clear(); energy_bounds_.insert(energy_bounds_.begin(), bounds.begin(), bounds.end()); + // if the mesh is set, allocate space for weight window bounds + if (mesh_idx_ != C_NONE) + allocate_ww_bounds(); } void WeightWindows::set_particle_type(ParticleType p_type) @@ -325,6 +333,7 @@ void WeightWindows::set_mesh(int32_t mesh_idx) fatal_error(fmt::format("Could not find a mesh for index {}", mesh_idx)); mesh_idx_ = mesh_idx; + allocate_ww_bounds(); } void WeightWindows::set_mesh(const std::unique_ptr& mesh) @@ -743,43 +752,6 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) e_bounds.push_back(data::energy_max[p_type]); } - // create a tally based on the WWG information - Tally* ww_tally = Tally::create(); - tally_idx_ = model::tally_map[ww_tally->id()]; - ww_tally->set_scores({"flux"}); - - // see if there's already a mesh filter using this mesh - bool found_mesh_filter = false; - for (const auto& f : model::tally_filters) { - if (f->type() == FilterType::MESH) { - const auto* mesh_filter = dynamic_cast(f.get()); - if (mesh_filter->mesh() == mesh_idx && !mesh_filter->translated()) { - ww_tally->add_filter(f.get()); - found_mesh_filter = true; - break; - } - } - } - - if (!found_mesh_filter) { - auto mesh_filter = Filter::create("mesh"); - openmc_mesh_filter_set_mesh(mesh_filter->index(), model::mesh_map[mesh_id]); - ww_tally->add_filter(mesh_filter); - } - - if (e_bounds.size() > 0) { - auto energy_filter = Filter::create("energy"); - openmc_energy_filter_set_bins( - energy_filter->index(), e_bounds.size(), e_bounds.data()); - ww_tally->add_filter(energy_filter); - } - - // add a particle filter - auto particle_filter = Filter::create("particle"); - auto pf = dynamic_cast(particle_filter); - pf->set_particles({&particle_type, 1}); - ww_tally->add_filter(particle_filter); - // set method and parameters for updates method_ = get_node_value(node, "method"); if (method_ == "magic") { @@ -815,13 +787,59 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) // create a matching weight windows object auto wws = WeightWindows::create(); ww_idx_ = wws->index(); + wws->set_mesh(mesh_idx); if (e_bounds.size() > 0) wws->set_energy_bounds(e_bounds); - wws->set_mesh(model::mesh_map[mesh_id]); wws->set_particle_type(particle_type); wws->set_defaults(); } +void WeightWindowsGenerator::create_tally() +{ + const auto& wws = variance_reduction::weight_windows[ww_idx_]; + + // create a tally based on the WWG information + Tally* ww_tally = Tally::create(); + tally_idx_ = model::tally_map[ww_tally->id()]; + ww_tally->set_scores({"flux"}); + + int32_t mesh_id = wws->mesh()->id(); + int32_t mesh_idx = model::mesh_map.at(mesh_id); + // see if there's already a mesh filter using this mesh + bool found_mesh_filter = false; + for (const auto& f : model::tally_filters) { + if (f->type() == FilterType::MESH) { + const auto* mesh_filter = dynamic_cast(f.get()); + if (mesh_filter->mesh() == mesh_idx && !mesh_filter->translated()) { + ww_tally->add_filter(f.get()); + found_mesh_filter = true; + break; + } + } + } + + if (!found_mesh_filter) { + auto mesh_filter = Filter::create("mesh"); + openmc_mesh_filter_set_mesh(mesh_filter->index(), model::mesh_map[mesh_id]); + ww_tally->add_filter(mesh_filter); + } + + const auto& e_bounds = wws->energy_bounds(); + if (e_bounds.size() > 0) { + auto energy_filter = Filter::create("energy"); + openmc_energy_filter_set_bins( + energy_filter->index(), e_bounds.size(), e_bounds.data()); + ww_tally->add_filter(energy_filter); + } + + // add a particle filter + auto particle_type = wws->particle_type(); + auto particle_filter = Filter::create("particle"); + auto pf = dynamic_cast(particle_filter); + pf->set_particles({&particle_type, 1}); + ww_tally->add_filter(particle_filter); +} + void WeightWindowsGenerator::update() const { const auto& wws = variance_reduction::weight_windows[ww_idx_]; @@ -845,6 +863,17 @@ void WeightWindowsGenerator::update() const // complete } +//============================================================================== +// Non-member functions +//============================================================================== + +void finalize_variance_reduction() +{ + for (const auto& wwg : variance_reduction::weight_windows_generators) { + wwg->create_tally(); + } +} + //============================================================================== // C API //============================================================================== diff --git a/tests/regression_tests/weightwindows/generators/test.py b/tests/regression_tests/weightwindows/generators/test.py index f5bfc4dcd..1d516db8d 100644 --- a/tests/regression_tests/weightwindows/generators/test.py +++ b/tests/regression_tests/weightwindows/generators/test.py @@ -28,13 +28,20 @@ def test_ww_generator(run_in_tmpdir): energy_bounds = np.linspace(0.0, 1e6, 70) particle = 'neutron' + # include another tally to make sure user-specified tallies and those automaticaly + # created by weight window generators can coexist + tally = openmc.Tally() + ef = openmc.EnergyFilter(energy_bounds) + tally.filters = [ef] + tally.scores = ['flux'] + model.tallies = [tally] + wwg = openmc.WeightWindowGenerator(mesh, energy_bounds, particle) - wwg.update_parameters = {'ratio' : 5.0, 'threshold': 0.8, 'value' : 'mean'} + wwg.update_parameters = {'ratio': 5.0, 'threshold': 0.8, 'value': 'mean'} model.settings.weight_window_generators = wwg - model.export_to_xml() - model.run() + # we test the effectiveness of the update method elsewhere, so # just test that the generation happens successfully here assert os.path.exists('weight_windows.h5') diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index f8db66820..8af843cbd 100644 --- a/tests/unit_tests/weightwindows/test.py +++ b/tests/unit_tests/weightwindows/test.py @@ -241,7 +241,21 @@ def test_roundtrip(run_in_tmpdir, model, wws): assert(ww_out == ww_in) -def test_ww_attrs(run_in_tmpdir, model): +def test_ww_attrs_python(model): + mesh = openmc.RegularMesh.from_domain(model.geometry) + lower_bounds = np.ones(mesh.dimension) + + # ensure that creation of weight window objects with default arg values + # is successful + wws = openmc.WeightWindows(mesh, lower_bounds, upper_bound_ratio=10.0) + + assert wws.energy_bounds == None + + wwg = openmc.WeightWindowGenerator(mesh) + + assert wwg.energy_bounds == None + +def test_ww_attrs_capi(run_in_tmpdir, model): model.export_to_xml() openmc.lib.init()