From 610a5d8c3485034626ce54ce87a88eaa8071cc23 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 9 Jun 2023 10:47:27 -0500 Subject: [PATCH] Global Weight Window Generation via CAPI (#2359) Co-authored-by: Jonathan Shimwell Co-authored-by: Paul Romano --- docs/source/io_formats/index.rst | 1 + docs/source/io_formats/settings.rst | 73 +- docs/source/io_formats/statepoint.rst | 2 + docs/source/io_formats/weight_windows.rst | 39 + docs/source/pythonapi/base.rst | 14 +- docs/source/pythonapi/capi.rst | 2 + include/openmc/capi.h | 33 + include/openmc/constants.h | 3 +- include/openmc/mesh.h | 52 +- include/openmc/output.h | 22 + include/openmc/settings.h | 6 +- include/openmc/tallies/filter_mesh.h | 6 +- include/openmc/tallies/tally.h | 44 + include/openmc/weight_windows.h | 259 +++-- openmc/lib/__init__.py | 2 + openmc/lib/core.py | 50 +- openmc/lib/filter.py | 18 +- openmc/lib/mesh.py | 9 +- openmc/lib/settings.py | 1 + openmc/lib/tally.py | 24 + openmc/lib/weight_windows.py | 323 ++++++ openmc/settings.py | 68 +- openmc/source.py | 47 +- openmc/tallies.py | 1 - openmc/tracks.py | 3 +- openmc/weight_windows.py | 302 +++++- src/mesh.cpp | 64 +- src/settings.cpp | 33 +- src/simulation.cpp | 17 + src/tallies/filter.cpp | 9 + src/tallies/filter_particle.cpp | 19 + src/tallies/tally.cpp | 61 ++ src/weight_windows.cpp | 930 +++++++++++++++++- .../track_output/results_true.dat | 6 +- .../weightwindows/generators/__init__.py | 0 .../weightwindows/generators/test.py | 54 + tests/unit_tests/weightwindows/test.py | 47 +- tests/unit_tests/weightwindows/test_ww_gen.py | 316 ++++++ 38 files changed, 2789 insertions(+), 171 deletions(-) create mode 100644 docs/source/io_formats/weight_windows.rst create mode 100644 openmc/lib/weight_windows.py create mode 100644 tests/regression_tests/weightwindows/generators/__init__.py create mode 100644 tests/regression_tests/weightwindows/generators/test.py create mode 100644 tests/unit_tests/weightwindows/test_ww_gen.py diff --git a/docs/source/io_formats/index.rst b/docs/source/io_formats/index.rst index 8c89bd2a1..4bbaa961a 100644 --- a/docs/source/io_formats/index.rst +++ b/docs/source/io_formats/index.rst @@ -51,3 +51,4 @@ Output Files track voxel volume + weight_windows diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 0431c95ab..2a85ebc8a 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -1066,7 +1066,7 @@ sub-elements/attributes: :particle_type: The particle that the weight windows will apply to (e.g., 'neutron') - *Default*: None + *Default*: 'neutron' :energy_bins: Monotonically increasing list of bounding energies in [eV] to be used for @@ -1104,3 +1104,74 @@ sub-elements/attributes: Threshold below which particles will be terminated *Default*: :math:`10^{-38}` + +-------------------------------------- +```` Element +-------------------------------------- + +The ```` element provides information for creating a set of +mesh-based weight windows. + + :mesh: + ID of a mesh that is to be used for the weight windows spatial bins + + *Default*: None + + :energy_bounds: + The weight window energy bounds. If not present, the max/min energy of the + cross section data is applied as a single energy bin. + + *Default*: None + + :particle_type: + The particle that the weight windows will apply to (e.g., 'neutron') + + *Default*: neutron + + :max_realizations: + The number of tally realizations after which the weight windows will stop updating. + + *Default*: 1 + + :update_interval: + The number of tally realizations between weight window updates. + + *Default*: 1 + + :on_the_fly: + Controls whether or not the tally results are reset after a weight window update. + + *Default*: true + + :method: + Method used to update weight window values (currently only 'magic' is supported) + + *Default*: magic + + :update_parameters: + Method-specific update parameters used when generating/updating weight windows. + + For MAGIC: + + :value: + The type of tally value to use when creating weight windows (one of 'mean' or 'rel_err') + + *Default*: 'mean' + + :threshold: + The relative error threshold above which tally results will be ignored. + + *Default*: 1.0 + + :ratio: + The ratio of the lower to upper weight window bounds. + + *Default*: 5.0 + +-------------------------------------- +```` Element +-------------------------------------- + + The ``weight_windows_file`` element has no attributes and contains the path to + a weight windows HDF5 file to load during simulation initialization. + diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 89d973a2b..37a379b52 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -68,6 +68,8 @@ The current version of the statepoint file format is 17.0. :Attributes: - **n_meshes** (*int*) -- Number of meshes in the problem. - **ids** (*int[]*) -- User-defined unique ID of each mesh. +.. _mesh-spec-hdf5: + **/tallies/meshes/mesh /** :Datasets: - **type** (*char[]*) -- Type of mesh. diff --git a/docs/source/io_formats/weight_windows.rst b/docs/source/io_formats/weight_windows.rst new file mode 100644 index 000000000..302e89318 --- /dev/null +++ b/docs/source/io_formats/weight_windows.rst @@ -0,0 +1,39 @@ +.. _io_weight_windows: + +==================== +Weight Window Format +==================== + +The current revision of the weight window file format is 1.0. + +**/** + +:Attributes: - **filetype** (*char[]*) -- String indicating the type of file. + - **version** (*int[2]*) -- Major and minor version of the weight + window file format. + +**/weight_windows/** + +:Attributes: - **n_weight_windows** (*int*) -- Number of weight window objects in the file. + - **ids** (*int[]*) -- Unique IDs of weight window objects in the file. + +**/weight_windows/weight_windows_/** + +:Datasets: - **mesh** (*int*) -- ID of the mesh associated with the weight window object. + - **particle_type** (*char[]*) -- Particle type to which the weight windows apply. + - **energy_bounds** (*double[]*) -- Energy bounds of the weight windows in [eV] + - **lower_ww_bounds** (*double[]*) -- Weight window lower bounds. + - **upper_ww_bounds** (*double[]*) -- Weight window upper bounds. + - **survival_ratio** (*double*) -- Weight window survival ratio. + - **max_lower_bound_ratio** (*double*) -- Maximum particle weight to lower weight window bound ratio. + - **max_split** (*int*) -- Maximum number of splits per weight window check. + - **weight_cutoff** (*double*) -- Particle weight cutoff. + +**/meshes/** + +:Attributes: - **n_meshes** (*int*) -- Number of meshes in the file. + - **ids** (*int[]*) -- User-defined unique ID of each mesh. + +**/meshes/mesh /** + +Please see the section on **/tallies/meshes/** in the :doc:`statepoint`. diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index c7304f462..ffd521a75 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -24,7 +24,6 @@ Simulation Settings openmc.Source openmc.SourceParticle openmc.VolumeCalculation - openmc.WeightWindows openmc.Settings .. autosummary:: @@ -230,6 +229,19 @@ Various classes may be created when performing tally slicing and/or arithmetic: openmc.arithmetic.AggregateNuclide openmc.arithmetic.AggregateFilter +Variance Reduction +------------------ + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass + + openmc.WeightWindows + openmc.WeightWindowGenerator + openmc.hdf5_to_wws + + Coarse Mesh Finite Difference Acceleration ------------------------------------------ diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 89fa167ff..bce647ccb 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -15,6 +15,8 @@ Functions calculate_volumes current_batch export_properties + export_weight_windows + import_weight_windows finalize find_cell find_material diff --git a/include/openmc/capi.h b/include/openmc/capi.h index ce58bed12..ac69c652f 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -50,6 +50,7 @@ int openmc_extend_meshes( int openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_filter_get_id(int32_t index, int32_t* id); int openmc_filter_get_type(int32_t index, char* type); +int openmc_filter_get_num_bins(int32_t index, int* n_bins); int openmc_filter_set_id(int32_t index, int32_t id); int openmc_finalize(); int openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance); @@ -163,6 +164,36 @@ int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides); int openmc_tally_set_scores(int32_t index, int n, const char** scores); int openmc_tally_set_type(int32_t index, const char* type); int openmc_tally_set_writable(int32_t index, bool writable); +int openmc_get_weight_windows_index(int32_t id, int32_t* idx); +int openmc_weight_windows_get_id(int32_t index, int32_t* id); +int openmc_weight_windows_set_id(int32_t index, int32_t id); + +//! Updates weight window values using the specified tally +//! \param[in] ww_idx Index of the weight window object +//! \param[in] tally_idx Index of the tally to use for the update +//! \param[in] value Tally value to use for the update (one of 'mean', +//! 'rel_err') \param[in] threshold Relative error threshold for applied results +//! \param[in] ratio Upper to lower weight window bound ratio +int openmc_weight_windows_update_magic(int32_t ww_idx, int32_t tally_idx, + const char* value, double threshold, double ratio); + +int openmc_extend_weight_windows( + int32_t n, int32_t* index_start, int32_t* index_end); +int openmc_weight_windows_get_mesh(int32_t index, int32_t* mesh_idx); +int openmc_weight_windows_set_mesh(int32_t index, int32_t mesh_idx); +int openmc_weight_windows_set_energy_bounds( + int32_t index, double* e_bounds, size_t e_bounds_size); +int openmc_weight_windows_get_energy_bounds( + int32_t index, const double** e_bounds, size_t* e_bounds_size); +int openmc_weight_windows_set_particle(int32_t index, int particle); +int openmc_weight_windows_get_particle(int32_t index, int* particle); +int openmc_weight_windows_get_bounds(int32_t index, const double** lower_bounds, + const double** upper_bounds, size_t* size); +int openmc_weight_windows_set_bounds(int32_t index, const double* lower_bounds, + const double* upper_bounds, size_t size); +size_t openmc_weight_windows_size(); +int openmc_weight_windows_export(const char* filename = nullptr); +int openmc_weight_windows_import(const char* filename = nullptr); int openmc_zernike_filter_get_order(int32_t index, int* order); int openmc_zernike_filter_get_params( int32_t index, double* x, double* y, double* r); @@ -170,6 +201,8 @@ int openmc_zernike_filter_set_order(int32_t index, int order); int openmc_zernike_filter_set_params( int32_t index, const double* x, const double* y, const double* r); +int openmc_particle_filter_get_bins(int32_t idx, int bins[]); + //! Sets the mesh and energy grid for CMFD reweight //! \param[in] meshtyally_id id of CMFD Mesh Tally //! \param[in] cmfd_indices indices storing spatial and energy dimensions of diff --git a/include/openmc/constants.h b/include/openmc/constants.h index b0031fae3..6a77396b7 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -33,6 +33,7 @@ constexpr array VERSION_VOLUME {1, 0}; constexpr array VERSION_VOXEL {2, 0}; constexpr array VERSION_MGXS_LIBRARY {1, 0}; constexpr array VERSION_PROPERTIES {1, 0}; +constexpr array VERSION_WEIGHT_WINDOWS {1, 0}; // ============================================================================ // ADJUSTABLE PARAMETERS @@ -276,7 +277,7 @@ enum class MgxsType { // ============================================================================ // TALLY-RELATED CONSTANTS -enum class TallyResult { VALUE, SUM, SUM_SQ }; +enum class TallyResult { VALUE, SUM, SUM_SQ, SIZE }; enum class TallyType { VOLUME, MESH_SURFACE, SURFACE }; diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index cef4e0518..47f54c023 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -88,12 +88,6 @@ public: //! \return sampled position within tet virtual Position sample(uint64_t* seed, int32_t bin) const = 0; - //! Get the volume of a mesh bin - // - //! \param[in] bin Bin to return the volume for - //! \return Volume of the bin - virtual double volume(int bin) const = 0; - //! Determine which bins were crossed by a particle // //! \param[in] r0 Previous position of the particle @@ -125,6 +119,8 @@ public: //! Get the number of mesh cell surfaces. virtual int n_surface_bins() const = 0; + int32_t id() const { return id_; } + //! Set the mesh ID void set_id(int32_t id = -1); @@ -149,7 +145,15 @@ public: //! \param[in] bin Mesh bin to generate a label for virtual std::string bin_label(int bin) const = 0; - //! Return the mesh type + //! Get the volume of a mesh bin + // + //! \param[in] bin Bin to return the volume for + //! \return Volume of the bin + virtual double volume(int bin) const = 0; + + //! Volumes of all elements in the mesh in bin ordering + vector volumes() const; + virtual std::string get_mesh_type() const = 0; // Data members @@ -181,8 +185,6 @@ public: Position sample(uint64_t* seed, int32_t bin) const override; - double volume(int bin) const override; - int get_bin(Position r) const override; int n_bins() const override; @@ -260,6 +262,16 @@ public: //! Get shape as xt::xtensor xt::xtensor get_x_shape() const; + double volume(int bin) const override + { + return this->volume(get_indices_from_bin(bin)); + } + + //! Get the volume of a specified element + //! \param[in] ijk Mesh index to return the volume for + //! \return Volume of the bin + virtual double volume(const MeshIndex& ijk) const = 0; + // Data members xt::xtensor lower_left_; //!< Lower-left coordinates of mesh xt::xtensor upper_right_; //!< Upper-right coordinates of mesh @@ -331,8 +343,12 @@ public: xt::xtensor count_sites( const SourceSite* bank, int64_t length, bool* outside) const; + //! Return the volume for a given mesh index + double volume(const MeshIndex& ijk) const override; + // Data members double volume_frac_; //!< Volume fraction of each mesh element + double element_volume_; //!< Volume of each mesh element xt::xtensor width_; //!< Width of each mesh element }; @@ -370,9 +386,13 @@ public: //! \param[in] i Direction index double negative_grid_boundary(const MeshIndex& ijk, int i) const; - array, 3> grid_; + //! Return the volume for a given mesh index + double volume(const MeshIndex& ijk) const override; int set_grid(); + + // Data members + array, 3> grid_; }; class CylindricalMesh : public PeriodicStructuredMesh { @@ -398,6 +418,8 @@ public: void to_hdf5(hid_t group) const override; + double volume(const MeshIndex& ijk) const override; + array, 3> grid_; int set_grid(); @@ -478,6 +500,8 @@ private: } } + double volume(const MeshIndex& ijk) const override; + inline int sanitize_theta(int idx) const { return sanitize_angular_index(idx, full_theta_, shape_[1]); @@ -634,6 +658,10 @@ public: std::vector connectivity(int id) const override; + //! Get the volume of a mesh bin + // + //! \param[in] bin Bin to return the volume for + //! \return Volume of the bin double volume(int bin) const override; private: @@ -791,6 +819,10 @@ public: std::vector connectivity(int id) const override; + //! Get the volume of a mesh bin + // + //! \param[in] bin Bin to return the volume for + //! \return Volume of the bin double volume(int bin) const override; libMesh::MeshBase* mesh_ptr() const { return m_; }; diff --git a/include/openmc/output.h b/include/openmc/output.h index cf5b49607..1ef95a887 100644 --- a/include/openmc/output.h +++ b/include/openmc/output.h @@ -59,3 +59,25 @@ void write_tallies(); } // namespace openmc #endif // OPENMC_OUTPUT_H + +////////////////////////////////////// +// Custom formatters +////////////////////////////////////// +namespace fmt { + +template +struct formatter> { + template + constexpr auto parse(ParseContext& ctx) + { + return ctx.begin(); + } + + template + auto format(const std::array& arr, FormatContext& ctx) + { + return format_to(ctx.out(), "({}, {})", arr[0], arr[1]); + } +}; + +} // namespace fmt \ No newline at end of file diff --git a/include/openmc/settings.h b/include/openmc/settings.h index f6cbe3138..3ce13c89f 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -29,7 +29,7 @@ extern bool confidence_intervals; //!< use confidence intervals for results? extern bool create_fission_neutrons; //!< create fission neutrons (fixed source)? extern bool create_delayed_neutrons; //!< create delayed fission neutrons? -extern "C" bool cmfd_run; //!< is a CMFD run? +extern "C" bool cmfd_run; //!< is a CMFD run? extern bool delayed_photon_scaling; //!< Scale fission photon yield to include delayed extern "C" bool entropy_on; //!< calculate Shannon entropy? @@ -58,7 +58,7 @@ extern "C" bool trigger_on; //!< tally triggers enabled? extern bool trigger_predict; //!< predict batches for triggers? extern bool ufs_on; //!< uniform fission site method on? extern bool urr_ptables_on; //!< use unresolved resonance prob. tables? -extern bool weight_windows_on; //!< are weight windows are enabled? +extern "C" bool weight_windows_on; //!< are weight windows are enabled? extern bool write_all_tracks; //!< write track files for every particle? extern bool write_initial_source; //!< write out initial source file? @@ -69,6 +69,8 @@ extern std::string path_output; //!< directory where output files are written extern std::string path_particle_restart; //!< path to a particle restart file extern std::string path_sourcepoint; //!< path to a source file extern std::string path_statepoint; //!< path to a statepoint file +extern std::string weight_windows_file; //!< Location of weight window file to + //!< load on simulation initialization // This is required because the c_str() may not be the first thing in // std::string. Sometimes it is, but it seems libc++ may not be like that diff --git a/include/openmc/tallies/filter_mesh.h b/include/openmc/tallies/filter_mesh.h index e3bcd7c20..2f48b6d4d 100644 --- a/include/openmc/tallies/filter_mesh.h +++ b/include/openmc/tallies/filter_mesh.h @@ -55,9 +55,9 @@ protected: //---------------------------------------------------------------------------- // Data members - int32_t mesh_; - bool translated_ {false}; - Position translation_ {0.0, 0.0, 0.0}; + int32_t mesh_; //!< Index of the mesh + bool translated_ {false}; //!< Whether or not the filter is translated + Position translation_ {0.0, 0.0, 0.0}; //!< Filter translation }; } // namespace openmc diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 6ad3bef6f..303b1ebb1 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -35,6 +35,8 @@ public: void set_id(int32_t id); + int id() const { return id_; } + void set_active(bool active) { active_ = active; } void set_multiply_density(bool value) { multiply_density_ = value; } @@ -45,16 +47,52 @@ public: void set_scores(const vector& scores); + std::vector scores() const; + + int32_t n_scores() const { return scores_.size(); } + void set_nuclides(pugi::xml_node node); void set_nuclides(const vector& nuclides); + const xt::xtensor& results() const { return results_; } + //! returns vector of indices corresponding to the tally this is called on const vector& filters() const { return filters_; } + //! returns a vector of filter types for the tally + std::vector filter_types() const; + + //! returns a mapping of filter types to index into the tally's filters + std::unordered_map filter_indices() const; + //! \brief Returns the tally filter at index i int32_t filters(int i) const { return filters_[i]; } + //! \brief Return a const pointer to a filter instance based on type. Always + //! returns the first matching filter type + template + const T* get_filter() const + { + const T* out; + for (auto filter_idx : filters_) { + if ((out = dynamic_cast(model::tally_filters[filter_idx].get()))) + return out; + } + return nullptr; + } + + template + const T* get_filter(int idx) const + { + if (const T* out = dynamic_cast(model::tally_filters[filters_.at(idx)])) + return out; + return nullptr; + } + + //! \brief Check if this tally has a specified type of filter + bool has_filter(FilterType filter_type) const; + void set_filters(gsl::span filters); //! Given already-set filters, set the stride lengths @@ -81,6 +119,12 @@ public: void accumulate(); + //! return the index of a score specified by name + int score_index(const std::string& score) const; + + //! Tally results reshaped according to filter sizes + xt::xarray get_reshaped_data() const; + //! A string representing the i-th score on this tally std::string score_name(int score_idx) const; diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 45c1de919..75073e917 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -4,6 +4,7 @@ #include #include +#include #include #include @@ -11,104 +12,214 @@ #include "openmc/memory.h" #include "openmc/mesh.h" #include "openmc/particle.h" +#include "openmc/tallies/tally.h" #include "openmc/vector.h" -namespace openmc { + namespace openmc +{ -//============================================================================== -// Constants -//============================================================================== + enum class WeightWindowUpdateMethod { + MAGIC, + }; -constexpr double DEFAULT_WEIGHT_CUTOFF {1.0e-38}; // default low weight cutoff + //============================================================================== + // Constants + //============================================================================== -//============================================================================== -// Non-member functions -//============================================================================== + constexpr double DEFAULT_WEIGHT_CUTOFF {1.0e-38}; // default low weight cutoff -//! Apply weight windows to a particle -//! \param[in] p Particle to apply weight windows to -void apply_weight_windows(Particle& p); + //============================================================================== + // Non-member functions + //============================================================================== -//! Free memory associated with weight windows -void free_memory_weight_windows(); + //! Apply weight windows to a particle + //! \param[in] p Particle to apply weight windows to + void apply_weight_windows(Particle & p); -//============================================================================== -// Global variables -//============================================================================== + //! Free memory associated with weight windows + void free_memory_weight_windows(); -class WeightWindows; + //============================================================================== + // Global variables + //============================================================================== -namespace variance_reduction { + class WeightWindows; + class WeightWindowsGenerator; -extern std::unordered_map ww_map; -extern vector> weight_windows; + namespace variance_reduction { -} // namespace variance_reduction + extern std::unordered_map ww_map; + extern vector> weight_windows; + extern vector> weight_windows_generators; -//============================================================================== -//! Individual weight window information -//============================================================================== + } // namespace variance_reduction -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}; + //============================================================================== + //! Individual weight window information + //============================================================================== - //! Whether the weight window is in a valid state - bool is_valid() const { return lower_weight >= 0.0; } + 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}; - //! Adjust the weight window by a constant factor - void scale(double factor) - { - lower_weight *= factor; - upper_weight *= factor; - } -}; + //! Whether the weight window is in a valid state + bool is_valid() const { return lower_weight >= 0.0; } -//============================================================================== -//! Weight window settings -//============================================================================== + //! Adjust the weight window by a constant factor + void scale(double factor) + { + lower_weight *= factor; + upper_weight *= factor; + } + }; -class WeightWindows { -public: - // Constructors - WeightWindows(); - WeightWindows(pugi::xml_node node); + //============================================================================== + //! Weight window settings + //============================================================================== - // Methods + 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); - //! Set the weight window ID - void set_id(int32_t id = -1); + //---------------------------------------------------------------------------- + // Methods + private: + template + void check_bounds(const T& lower, const T& upper) 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; + template + void check_bounds(const T& lower) 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; + void check_tally_update_compatibility(const Tally* tally); - // Accessors - int32_t id() const { return id_; } - const Mesh& mesh() const { return *model::meshes[mesh_idx_]; } + public: + //! Set the weight window ID + void set_id(int32_t id = -1); -private: - // Data members - int32_t id_; //!< Unique ID - ParticleType particle_type_; //!< Particle type to apply weight windows to - vector energy_bounds_; //!< Energy boundaries [eV] - vector lower_ww_; //!< Lower weight window bounds - vector upper_ww_; //!< Upper weight window bounds - 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 -}; + void set_energy_bounds(gsl::span bounds); + + void set_mesh(const std::unique_ptr& mesh); + + void set_mesh(const Mesh* mesh); + + void set_mesh(int32_t mesh_idx); + + //! 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); + + // 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; + + //! Retrieve the weight window for a particle + //! \param[in] p Particle to get weight window for + WeightWindow get_weight_window(const Particle& p) const; + + std::array bounds_size() const; + + const vector& energy_bounds() const { return energy_bounds_; } + + void set_bounds(const xt::xtensor& lower_ww_bounds, + const xt::xtensor& upper_bounds); + + void set_bounds(const xt::xtensor& lower_bounds, double ratio); + + void set_bounds(gsl::span lower_bounds, + gsl::span upper_bounds); + + void set_bounds(gsl::span lower_bounds, double ratio); + + void set_particle_type(ParticleType p_type); + + //---------------------------------------------------------------------------- + // Accessors + int32_t id() const { return id_; } + int32_t& id() { return id_; } + + int32_t index() const { return index_; } + + vector& energy_bounds() { return energy_bounds_; } + + const std::unique_ptr& mesh() const + { + return model::meshes[mesh_idx_]; + } + + const xt::xtensor& lower_ww_bounds() const { return lower_ww_; } + xt::xtensor& lower_ww_bounds() { return lower_ww_; } + + const xt::xtensor& upper_ww_bounds() const { return upper_ww_; } + xt::xtensor& upper_ww_bounds() { return upper_ww_; } + + ParticleType particle_type() const { return particle_type_; } + + 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 + }; + + class WeightWindowsGenerator { + public: + // Constructors + WeightWindowsGenerator(pugi::xml_node node); + + // Methods + void update() const; + + // 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 + + // MAGIC update parameters + std::string tally_value_ { + "mean"}; // bool: return self._weight_windows_on + @property + def weight_windows_file(self) -> Optional[PathLike]: + return self._weight_windows_file + + @weight_windows_file.setter + def weight_windows_file(self, value: PathLike): + cv.check_type('weight windows file', value, (str, Path)) + self._weight_windows_file = value + + @property + def weight_window_generators(self) -> typing.List[WeightWindowGenerator]: + return self._weight_window_generators + + @weight_window_generators.setter + def weight_window_generators(self, wwgs): + if not isinstance(wwgs, MutableSequence): + wwgs = [wwgs] + self._weight_window_generators = cv.CheckedList(WeightWindowGenerator, 'weight window generators', wwgs) + @property def max_splits(self) -> int: return self._max_splits @@ -922,7 +954,7 @@ class Settings: self._weight_windows = cv.CheckedList(WeightWindows, 'weight windows', value) @weight_windows_on.setter - def weight_windows_on(self, value): + def weight_windows_on(self, value: bool): cv.check_type('weight windows on', value, bool) self._weight_windows_on = value @@ -1283,6 +1315,28 @@ class Settings: elem = ET.SubElement(root, "weight_windows_on") elem.text = str(self._weight_windows_on).lower() + def _create_weight_window_generators_subelement(self, root, mesh_memo=None): + if not self.weight_window_generators: + return + elem = ET.SubElement(root, 'weight_window_generators') + for wwg in self.weight_window_generators: + elem.append(wwg.to_xml_element()) + + # ensure that mesh elements are created if needed + for wwg in self.weight_window_generators: + if mesh_memo is not None and wwg.mesh.id in mesh_memo: + continue + + root.append(wwg.mesh.to_xml_element()) + if mesh_memo is not None: + mesh_memo.add(wwg.mesh) + + def _create_weight_windows_file_element(self, root): + if self.weight_windows_file is not None: + element = ET.Element("weight_windows_file") + element.text = self.weight_windows_file + root.append(element) + def _create_max_splits_subelement(self, root): if self._max_splits is not None: elem = ET.SubElement(root, "max_splits") @@ -1595,6 +1649,11 @@ class Settings: if text is not None: self.write_initial_source = text in ('true', '1') + def _weight_window_generators_from_xml_element(self, root, meshes=None): + for elem in root.iter('weight_windows_generator'): + wwg = WeightWindowGenerator.from_xml_element(elem, meshes) + self.weight_window_generators.append(wwg) + def _weight_windows_from_xml_element(self, root, meshes=None): for elem in root.findall('weight_windows'): ww = WeightWindows.from_xml_element(elem, root) @@ -1671,6 +1730,8 @@ class Settings: self._create_log_grid_bins_subelement(element) self._create_write_initial_source_subelement(element) self._create_weight_windows_subelement(element, mesh_memo) + self._create_weight_window_generators_subelement(element, mesh_memo) + self._create_weight_windows_file_element(element) self._create_max_splits_subelement(element) self._create_max_tracks_subelement(element) @@ -1764,6 +1825,7 @@ class Settings: settings._log_grid_bins_from_xml_element(elem) settings._write_initial_source_from_xml_element(elem) settings._weight_windows_from_xml_element(elem, meshes) + settings._weight_window_generators_from_xml_element(elem, meshes) settings._max_splits_from_xml_element(elem) settings._max_tracks_from_xml_element(elem) diff --git a/openmc/source.py b/openmc/source.py index 9851ea59d..cc68c9535 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -1,11 +1,11 @@ from collections.abc import Iterable -from enum import Enum +from enum import IntEnum from numbers import Real import warnings import typing # imported separately as py3.8 requires typing.Iterable # also required to prevent typing.Union namespace overwriting Union from typing import Optional, Sequence -import lxml.etree as ET +import lxml.etree as ET import numpy as np import h5py @@ -333,12 +333,53 @@ class Source: return source -class ParticleType(Enum): +class ParticleType(IntEnum): + """ + IntEnum class representing a particle type. Type + values mirror those found in the C++ class. + """ NEUTRON = 0 PHOTON = 1 ELECTRON = 2 POSITRON = 3 + @classmethod + def from_string(cls, value: str): + """ + Constructs a ParticleType instance from a string. + + Parameters + ---------- + value : str + The string representation of the particle type. + + Returns + ------- + The corresponding ParticleType instance. + """ + try: + return cls[value.upper()] + except KeyError: + raise ValueError(f"Invalid string for creation of {cls.__name__}: {value}") + + def __repr__(self) -> str: + """ + Returns a string representation of the ParticleType instance. + + Returns: + str: The lowercase name of the ParticleType instance. + """ + return self.name.lower() + + # needed for < Python 3.11 + def __str__(self) -> str: + return self.__repr__() + + # needed for <= 3.7, IntEnum will use the mixed-in type's `__format__` method otherwise + # this forces it to default to the standard object format, relying on __str__ under the hood + def __format__(self, spec): + return object.__format__(self, spec) + class SourceParticle: """Source particle diff --git a/openmc/tallies.py b/openmc/tallies.py index 85daaa17a..6f56fbb64 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3222,7 +3222,6 @@ class Tallies(cv.CheckedList): return element - def export_to_xml(self, path='tallies.xml'): """Create a tallies.xml file that can be used for a simulation. diff --git a/openmc/tracks.py b/openmc/tracks.py index 9b4b55247..61e5a7244 100644 --- a/openmc/tracks.py +++ b/openmc/tracks.py @@ -25,8 +25,7 @@ states : numpy.ndarray """ def _particle_track_repr(self): - name = self.particle.name.lower() - return f"" + return f"" ParticleTrack.__repr__ = _particle_track_repr diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index c09e1e862..2411b83f8 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -2,19 +2,20 @@ from __future__ import annotations from collections.abc import Iterable from numbers import Real, Integral import pathlib -import typing from typing import Iterable, List, Optional, Union, Dict +import warnings import lxml.etree as ET import numpy as np import h5py +import openmc from openmc.filter import _PARTICLES from openmc.mesh import MeshBase, RectilinearMesh, UnstructuredMesh import openmc.checkvalue as cv from openmc.checkvalue import PathLike -from ._xml import get_text +from ._xml import get_text, clean_indentation from .mixin import IDManagerMixin @@ -108,8 +109,8 @@ class WeightWindows(IDManagerMixin): used_ids = set() def __init__( - self, - mesh: MeshBase, + self, + mesh: MeshBase, lower_ww_bounds: Iterable[float], upper_ww_bounds: Optional[Iterable[float]] = None, upper_bound_ratio: Optional[float] = None, @@ -159,7 +160,7 @@ class WeightWindows(IDManagerMixin): def __repr__(self) -> str: string = type(self).__name__ + '\n' string += '{: <16}=\t{}\n'.format('\tID', self._id) - string += '{: <16}=\t{}\n'.format('\tMesh:', self.mesh) + string += '{: <16}=\t{}\n'.format('\tMesh', self.mesh) string += '{: <16}=\t{}\n'.format('\tParticle Type', self._particle_type) string += '{: <16}=\t{}\n'.format('\tEnergy Bounds', self._energy_bounds) string += '{: <16}=\t{}\n'.format('\tLower WW Bounds', self._lower_ww_bounds) @@ -286,7 +287,7 @@ class WeightWindows(IDManagerMixin): @max_lower_bound_ratio.setter def max_lower_bound_ratio(self, val: float): cv.check_type('Maximum lower bound ratio', val, Real) - cv.check_greater_than('Maximum lower bound ratio', val, 1.0) + cv.check_greater_than('Maximum lower bound ratio', val, 1.0, equality=True) self._max_lower_bound_ratio = val @property @@ -424,10 +425,15 @@ class WeightWindows(IDManagerMixin): id = int(group.name.split('/')[-1].lstrip('weight_windows')) mesh_id = group['mesh'][()] + mesh = meshes[mesh_id] + ptype = group['particle_type'][()].decode() e_bounds = group['energy_bounds'][()] - lower_ww_bounds = group['lower_ww_bounds'][()] - upper_ww_bounds = group['upper_ww_bounds'][()] + # weight window bounds are stored with the shape (e, k, j, i) + # in C++ and HDF5 -- the opposite of how they are stored here + shape = (e_bounds.size - 1, *mesh.dimension[::-1]) + lower_ww_bounds = group['lower_ww_bounds'][()].reshape(shape).T + upper_ww_bounds = group['upper_ww_bounds'][()].reshape(shape).T survival_ratio = group['survival_ratio'][()] max_lower_bound_ratio = None @@ -438,7 +444,7 @@ class WeightWindows(IDManagerMixin): weight_cutoff = group['weight_cutoff'][()] return cls( - mesh=meshes[mesh_id], + mesh=mesh, lower_ww_bounds=lower_ww_bounds, upper_ww_bounds=upper_ww_bounds, energy_bounds=e_bounds, @@ -622,3 +628,281 @@ def wwinp_to_wws(path: PathLike) -> List[WeightWindows]: wws.append(ww) return wws + + +class WeightWindowGenerator: + """Class passed to setting to govern weight window generation + using the OpenMC executable + + Parameters + ---------- + mesh : :class:`openmc.MeshBase` + 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 + particle_type : {'neutron', 'photon'} + Particle type the weight windows apply to + + Attributes + ---------- + mesh : openmc.MeshBase + 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 + particle_type : {'neutron', 'photon'} + Particle type the weight windows apply to + method : {'magic'} + The weight window generation methodology applied during an update. Only + 'magic' is currently supported. + max_realizations : int + The upper limit for number of tally realizations when generating weight + windows. + update_interval : int + The number of tally realizations between updates. (default: 1) + update_parameters : dict + A set of parameters related to the update. + on_the_fly : bool + Whether or not to apply weight windows on the fly. (default: True) + """ + + _MAGIC_PARAMS = {'value': str, 'threshold': float, 'ratio': float} + + def __init__(self, mesh, energy_bounds=None, particle_type='neutron'): + self.mesh = mesh + if energy_bounds is not None: + self.energy_bounds = energy_bounds + self.particle_type = particle_type + self.max_realizations = 1 + + self._update_parameters = None + + self.method = 'magic' + self.particle_type = particle_type + self.update_interval = 1 + self.on_the_fly = True + + + def __repr__(self): + string = type(self).__name__ + '\n' + string += f'\t{"Mesh":<20}=\t{self.mesh.id}\n' + string += f'\t{"Particle:":<20}=\t{self.particle_type}\n' + string += f'\t{"Energy Bounds:":<20}=\t{self.energy_bounds}\n' + string += f'\t{"Method":<20}=\t{self.method}\n' + string += f'\t{"Max Realizations:":<20}=\t{self.max_realizations}\n' + string += f'\t{"Update Interval:":<20}=\t{self.update_interval}\n' + string += f'\t{"On The Fly:":<20}=\t{self.on_the_fly}\n' + if self.update_parameters is not None: + string += f'\t{"Update Parameters:":<20}\n\t\t\t{self.update_parameters}\n' + string + + return string + + @property + def mesh(self) -> openmc.MeshBase: + return self._mesh + + @mesh.setter + def mesh(self, m: openmc.MeshBase): + cv.check_type('mesh', m, openmc.MeshBase) + self._mesh = m + + @property + def energy_bounds(self) -> Iterable[Real]: + return self._energy_bounds + + @energy_bounds.setter + def energy_bounds(self, eb: Iterable[float]): + cv.check_type('energy bounds', eb, Iterable, Real) + self._energy_bounds = eb + + @property + def particle_type(self) -> str: + return self._particle_type + + @particle_type.setter + def particle_type(self, pt: str): + cv.check_value('particle type', pt, ('neutron', 'photon')) + self._particle_type = pt + + @property + def method(self) -> str: + return self._method + + @method.setter + def method(self, m: str): + cv.check_type('generation method', m, str) + cv.check_value('generation method', m, {'magic'}) + self._method = m + if self._update_parameters is not None: + try: + self._check_update_parameters() + except (TypeError, KeyError): + warnings.warn(f'Update parameters are invalid for the "{m}" method.') + + @property + def max_realizations(self) -> int: + return self._max_realizations + + @max_realizations.setter + def max_realizations(self, m: int): + cv.check_type('max tally realizations', m, Integral) + cv.check_greater_than('max tally realizations', m, 0) + self._max_realizations = m + + @property + def update_interval(self) -> int: + return self._update_interval + + @update_interval.setter + def update_interval(self, ui: int): + cv.check_type('update interval', ui, Integral) + cv.check_greater_than('update interval', ui , 0) + self._update_interval = ui + + @property + def update_parameters(self) -> dict: + return self._update_parameters + + def _check_update_parameters(self, params: dict): + if self.method == 'magic': + check_params = self._MAGIC_PARAMS + + for key, val in params.items(): + if key not in check_params: + raise ValueError(f'Invalid param "{key}" for {self.method} ' + 'weight window generation') + cv.check_type(f'weight window generation param: "{key}"', val, self._MAGIC_PARAMS[key]) + + @update_parameters.setter + def update_parameters(self, params: dict): + self._check_update_parameters(params) + self._update_parameters = params + + @property + def on_the_fly(self) -> bool: + return self._on_the_fly + + @on_the_fly.setter + def on_the_fly(self, otf: bool): + cv.check_type('on the fly generation', otf, bool) + self._on_the_fly = otf + + def _update_parameters_subelement(self, element: ET.Element): + if not self.update_parameters: + return + params_element = ET.SubElement(element, 'update_parameters') + for pname, value in self.update_parameters.items(): + param_element = ET.SubElement(params_element, pname) + param_element.text = str(value) + + @classmethod + def _sanitize_update_parameters(cls, method: str, update_parameters: dict): + """ + Attempt to convert update parameters to their appropriate types + + Parameters + ---------- + method : str + The update method for which these update parameters should comply + update_parameters : dict + The update parameters as-read from the XML node (keys: str, values: str) + """ + if method == 'magic': + check_params = cls._MAGIC_PARAMS + + for param, param_type in check_params.items(): + if param in update_parameters: + update_parameters[param] = param_type(update_parameters[param]) + + def to_xml_element(self): + """Creates a 'weight_window_generator' element to be written to an XML file. + """ + element = ET.Element('weight_windows_generator') + + mesh_elem = ET.SubElement(element, 'mesh') + mesh_elem.text = str(self.mesh.id) + if self.energy_bounds is not None: + subelement = ET.SubElement(element, 'energy_bounds') + subelement.text = ' '.join(str(e) for e in self.energy_bounds) + particle_elem = ET.SubElement(element, 'particle_type') + particle_elem.text = self.particle_type + realizations_elem = ET.SubElement(element, 'max_realizations') + realizations_elem.text = str(self.max_realizations) + update_interval_elem = ET.SubElement(element, 'update_interval') + update_interval_elem.text = str(self.update_interval) + otf_elem = ET.SubElement(element, 'on_the_fly') + otf_elem.text = str(self.on_the_fly).lower() + method_elem = ET.SubElement(element, 'method') + method_elem.text = self.method + if self.update_parameters is not None: + self._update_parameters_subelement(element) + + clean_indentation(element) + + return element + + @classmethod + def from_xml_element(cls, elem: ET.Element, meshes: dict) -> WeightWindowGenerator: + """ + Create a weight window generation object from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + meshes : dict + A dictionary with IDs as keys and openmc.MeshBase instances as values + + Returns + ------- + openmc.WeightWindowGenerator + """ + + mesh_id = int(get_text(elem, 'mesh')) + mesh = meshes[mesh_id] + + energy_bounds = [float(x) for x in get_text(elem, 'energy_bounds').split()] + particle_type = get_text(elem, 'particle_type') + + wwg = cls(mesh, energy_bounds, particle_type) + + wwg.max_realizations = int(get_text(elem, 'max_realizations')) + wwg.update_interval = int(get_text(elem, 'update_interval')) + wwg.on_the_fly = bool(get_text(elem, 'on_the_fly')) + wwg.method = get_text(elem, 'method') + + if elem.find('update_parameters'): + update_parameters = {} + params_elem = elem.find('update_parameters') + for entry in params_elem: + update_parameters[entry.tag] = entry.text + + cls._sanitize_update_parameters(wwg.method, update_parameters) + wwg.update_parameters = update_parameters + + return wwg + +def hdf5_to_wws(path='weight_windows.h5'): + """Create WeightWindows instances from a weight windows HDF5 file + + .. versionadded:: 0.13.4 + + Parameters + ---------- + path : cv.PathLike + Path to the weight windows hdf5 file + + Returns + ------- + list of openmc.WeightWindows + """ + + with h5py.File(path) as h5_file: + # read in all of the meshes in the mesh node + meshes = {} + for mesh_group in h5_file['meshes']: + mesh = MeshBase.from_hdf5(h5_file['meshes'][mesh_group]) + meshes[mesh.id] = mesh + return [WeightWindows.from_hdf5(ww, meshes) for ww in h5_file['weight_windows'].values()] \ No newline at end of file diff --git a/src/mesh.cpp b/src/mesh.cpp index e72a5cf7d..12052bed6 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -132,6 +132,15 @@ void Mesh::set_id(int32_t id) model::mesh_map[id] = model::meshes.size() - 1; } +vector Mesh::volumes() const +{ + vector volumes(n_bins()); + for (int i = 0; i < n_bins(); i++) { + volumes[i] = this->volume(i); + } + return volumes; +} + //============================================================================== // Structured Mesh implementation //============================================================================== @@ -377,11 +386,6 @@ Position StructuredMesh::sample(uint64_t* seed, int32_t bin) const fatal_error("Position sampling on structured meshes is not yet implemented"); } -double StructuredMesh::volume(int bin) const -{ - fatal_error("Unable to get volume of structured mesh, not yet implemented"); -} - int StructuredMesh::get_bin(Position r) const { // Determine indices @@ -722,6 +726,11 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node} // Set volume fraction volume_frac_ = 1.0 / xt::prod(shape)(); + + element_volume_ = 1.0; + for (int i = 0; i < n_dimension_; i++) { + element_volume_ *= width_[i]; + } } int RegularMesh::get_index_in_direction(double r, int i) const @@ -875,6 +884,11 @@ xt::xtensor RegularMesh::count_sites( return counts; } +double RegularMesh::volume(const MeshIndex& ijk) const +{ + return element_volume_; +} + //============================================================================== // RectilinearMesh implementation //============================================================================== @@ -1004,6 +1018,16 @@ void RectilinearMesh::to_hdf5(hid_t group) const close_group(mesh_group); } +double RectilinearMesh::volume(const MeshIndex& ijk) const +{ + double vol {1.0}; + + for (int i = 0; i < n_dimension_; i++) { + vol *= grid_[i][ijk[i] + 1] - grid_[i][ijk[i]]; + } + return vol; +} + //============================================================================== // CylindricalMesh implementation //============================================================================== @@ -1012,7 +1036,6 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : PeriodicStructuredMesh {node} { n_dimension_ = 3; - grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "phi_grid"); grid_[2] = get_node_array(node, "z_grid"); @@ -1249,6 +1272,20 @@ void CylindricalMesh::to_hdf5(hid_t group) const close_group(mesh_group); } +double CylindricalMesh::volume(const MeshIndex& ijk) const +{ + double r_i = grid_[0][ijk[0] - 1]; + double r_o = grid_[0][ijk[0]]; + + double phi_i = grid_[1][ijk[1] - 1]; + double phi_o = grid_[1][ijk[1]]; + + double z_i = grid_[2][ijk[2] - 1]; + double z_o = grid_[2][ijk[2]]; + + return 0.5 * (r_o * r_o - r_i * r_i) * (phi_o - phi_i) * (z_o - z_i); +} + //============================================================================== // SphericalMesh implementation //============================================================================== @@ -1523,6 +1560,21 @@ void SphericalMesh::to_hdf5(hid_t group) const close_group(mesh_group); } +double SphericalMesh::volume(const MeshIndex& ijk) const +{ + double r_i = grid_[0][ijk[0] - 1]; + double r_o = grid_[0][ijk[0]]; + + double theta_i = grid_[1][ijk[1] - 1]; + double theta_o = grid_[1][ijk[1]]; + + double phi_i = grid_[2][ijk[2] - 1]; + double phi_o = grid_[2][ijk[2]]; + + return (1.0 / 3.0) * (r_o * r_o * r_o - r_i * r_i * r_i) * + (std::cos(theta_i) - std::cos(theta_o)) * (phi_o - phi_i); +} + //============================================================================== // Helper functions for the C API //============================================================================== diff --git a/src/settings.cpp b/src/settings.cpp index 6d1444122..e16afbedc 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -83,6 +83,7 @@ std::string path_particle_restart; std::string path_sourcepoint; std::string path_statepoint; const char* path_statepoint_c {path_statepoint.c_str()}; +std::string weight_windows_file; int32_t n_inactive {0}; int32_t max_lost_particles {10}; @@ -904,11 +905,19 @@ void read_settings_xml(pugi::xml_node root) for (pugi::xml_node node_ww : root.children("weight_windows")) { variance_reduction::weight_windows.emplace_back( std::make_unique(node_ww)); - - // Enable weight windows by default if one or more are present - settings::weight_windows_on = true; } + // Enable weight windows by default if one or more are present + if (variance_reduction::weight_windows.size() > 0) + settings::weight_windows_on = true; + + // read weight windows from file + if (check_for_node(root, "weight_windows_file")) { + weight_windows_file = get_node_value(root, "weight_windows_file"); + } + + // read settings for weight windows value, this will override + // the automatic setting even if weight windows are present if (check_for_node(root, "weight_windows_on")) { weight_windows_on = get_node_value_bool(root, "weight_windows_on"); } @@ -920,6 +929,24 @@ void read_settings_xml(pugi::xml_node root) if (check_for_node(root, "max_tracks")) { settings::max_tracks = std::stoi(get_node_value(root, "max_tracks")); } + + // Create weight window generator objects + if (check_for_node(root, "weight_window_generators")) { + auto wwgs_node = root.child("weight_window_generators"); + for (pugi::xml_node node_wwg : + wwgs_node.children("weight_windows_generator")) { + variance_reduction::weight_windows_generators.emplace_back( + std::make_unique(node_wwg)); + } + // if any of the weight windows are intended to be generated otf, make sure + // they're applied + for (const auto& wwg : variance_reduction::weight_windows_generators) { + if (wwg->on_the_fly_) { + settings::weight_windows_on = true; + break; + } + } + } } void free_memory_settings() diff --git a/src/simulation.cpp b/src/simulation.cpp index d46787c28..310629936 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -24,6 +24,7 @@ #include "openmc/tallies/trigger.h" #include "openmc/timer.h" #include "openmc/track_output.h" +#include "openmc/weight_windows.h" #ifdef _OPENMP #include @@ -137,6 +138,11 @@ int openmc_simulation_init() } } + // load weight windows from file + if (!settings::weight_windows_file.empty()) { + openmc_weight_windows_import(settings::weight_windows_file.c_str()); + } + // Set flag indicating initialization is done simulation::initialized = true; return 0; @@ -175,6 +181,12 @@ int openmc_simulation_finalize() if (settings::output_tallies && mpi::master) write_tallies(); + // If weight window generators are present in this simulation, + // write a weight windows file + if (variance_reduction::weight_windows_generators.size() > 0) { + openmc_weight_windows_export(); + } + // Deactivate all tallies for (auto& t : model::tallies) { t->active_ = false; @@ -356,6 +368,11 @@ void finalize_batch() accumulate_tallies(); simulation::time_tallies.stop(); + // update weight windows if needed + for (const auto& wwg : variance_reduction::weight_windows_generators) { + wwg->update(); + } + // Reset global tally results if (simulation::current_batch <= settings::n_inactive) { xt::view(simulation::global_tallies, xt::all()) = 0.0; diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 5e8204b2a..d8efd590f 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -221,6 +221,15 @@ extern "C" int openmc_filter_get_type(int32_t index, char* type) return 0; } +extern "C" int openmc_filter_get_num_bins(int32_t index, int* n_bins) +{ + if (int err = verify_filter(index)) + return err; + + *n_bins = model::tally_filters[index]->n_bins(); + return 0; +} + extern "C" int openmc_get_filter_index(int32_t id, int32_t* index) { auto it = model::filter_map.find(id); diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index 3865824f6..56d93a47f 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -58,4 +58,23 @@ std::string ParticleFilter::text_label(int bin) const return fmt::format("Particle: {}", particle_type_to_str(p)); } +extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[]) +{ + if (int err = verify_filter(idx)) + return err; + + const auto& f = model::tally_filters[idx]; + auto pf = dynamic_cast(f.get()); + if (pf) { + const auto& particles = pf->particles(); + for (int i = 0; i < particles.size(); i++) { + bins[i] = static_cast(particles[i]); + } + } else { + set_errmsg("The filter at the specified index is not a ParticleFilter"); + return OPENMC_E_INVALID_ARGUMENT; + } + return 0; +} + } // namespace openmc diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 85d999fd4..9a7e52b95 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -358,6 +358,34 @@ void Tally::set_id(int32_t id) model::tally_map[id] = index_; } +std::vector Tally::filter_types() const +{ + std::vector filter_types; + for (auto idx : this->filters()) + filter_types.push_back(model::tally_filters[idx]->type()); + return filter_types; +} + +std::unordered_map Tally::filter_indices() const +{ + std::unordered_map filter_indices; + for (int i = 0; i < this->filters().size(); i++) { + const auto& f = model::tally_filters[this->filters(i)]; + + filter_indices[f->type()] = i; + } + return filter_indices; +} + +bool Tally::has_filter(FilterType filter_type) const +{ + for (auto idx : this->filters()) { + if (model::tally_filters[idx]->type() == filter_type) + return true; + } + return false; +} + void Tally::set_filters(gsl::span filters) { // Clear old data. @@ -692,6 +720,31 @@ void Tally::accumulate() } } +int Tally::score_index(const std::string& score) const +{ + for (int i = 0; i < scores_.size(); i++) { + if (this->score_name(i) == score) + return i; + } + return -1; +} + +xt::xarray Tally::get_reshaped_data() const +{ + std::vector shape; + for (auto f : filters()) { + shape.push_back(model::tally_filters[f]->n_bins()); + } + + // add number of scores and nuclides to tally + shape.push_back(results_.shape()[1]); + shape.push_back(results_.shape()[2]); + + xt::xarray reshaped_results = results_; + reshaped_results.reshape(shape); + return reshaped_results; +} + std::string Tally::score_name(int score_idx) const { if (score_idx < 0 || score_idx >= scores_.size()) { @@ -700,6 +753,14 @@ std::string Tally::score_name(int score_idx) const return reaction_name(scores_[score_idx]); } +std::vector Tally::scores() const +{ + std::vector score_names; + for (int score : scores_) + score_names.push_back(reaction_name(score)); + return score_names; +} + std::string Tally::nuclide_name(int nuclide_idx) const { if (nuclide_idx < 0 || nuclide_idx >= nuclides_.size()) { diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index e8e8f4b05..104c38815 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -1,12 +1,33 @@ #include "openmc/weight_windows.h" +#include +#include +#include +#include + +#include "xtensor/xindex_view.hpp" +#include "xtensor/xio.hpp" +#include "xtensor/xmasked_view.hpp" +#include "xtensor/xnoalias.hpp" +#include "xtensor/xstrided_view.hpp" +#include "xtensor/xview.hpp" + #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/hdf5_interface.h" +#include "openmc/mesh.h" +#include "openmc/message_passing.h" +#include "openmc/nuclide.h" +#include "openmc/output.h" #include "openmc/particle.h" #include "openmc/particle_data.h" #include "openmc/physics_common.h" #include "openmc/search.h" +#include "openmc/settings.h" +#include "openmc/tallies/filter_energy.h" +#include "openmc/tallies/filter_mesh.h" +#include "openmc/tallies/filter_particle.h" +#include "openmc/tallies/tally.h" #include "openmc/xml_interface.h" #include @@ -22,6 +43,7 @@ namespace variance_reduction { std::unordered_map ww_map; openmc::vector> weight_windows; +openmc::vector> weight_windows_generators; } // namespace variance_reduction @@ -108,6 +130,13 @@ void free_memory_weight_windows() // WeightWindowSettings implementation //============================================================================== +WeightWindows::WeightWindows(int32_t id) +{ + index_ = variance_reduction::weight_windows.size(); + set_id(id); + set_defaults(); +} + WeightWindows::WeightWindows(pugi::xml_node node) { // Make sure required elements are present @@ -132,11 +161,8 @@ WeightWindows::WeightWindows(pugi::xml_node node) mesh_idx_ = model::mesh_map.at(mesh_id); // energy bounds - energy_bounds_ = get_node_array(node, "energy_bounds"); - - // read the lower/upper weight bounds - lower_ww_ = get_node_array(node, "lower_ww_bounds"); - upper_ww_ = get_node_array(node, "upper_ww_bounds"); + if (check_for_node(node, "energy_bounds")) + energy_bounds_ = get_node_array(node, "energy_bounds"); // get the survival value - optional if (check_for_node(node, "survival_ratio")) { @@ -170,22 +196,81 @@ WeightWindows::WeightWindows(pugi::xml_node node) fatal_error("weight_cutoff must be less than 1"); } - // make sure that the upper and lower bounds have the same size - if (upper_ww_.size() != lower_ww_.size()) { - fatal_error("The upper and lower weight window lengths do not match."); + // read the lower/upper weight bounds + this->set_bounds(get_node_array(node, "lower_ww_bounds"), + get_node_array(node, "upper_ww_bounds")); + + set_defaults(); +} + +WeightWindows::~WeightWindows() +{ + variance_reduction::ww_map.erase(id()); +} + +WeightWindows* WeightWindows::create(int32_t id) +{ + variance_reduction::weight_windows.push_back(make_unique()); + auto wws = variance_reduction::weight_windows.back().get(); + variance_reduction::ww_map[wws->id()] = + variance_reduction::weight_windows.size() - 1; + return wws; +} + +WeightWindows* WeightWindows::from_hdf5( + hid_t wws_group, const std::string& group_name) +{ + // collect ID from the name of this group + hid_t ww_group = open_group(wws_group, group_name); + + auto wws = WeightWindows::create(); + + std::string particle_type; + read_dataset(ww_group, "particle_type", particle_type); + wws->particle_type_ = openmc::str_to_particle_type(particle_type); + + read_dataset(ww_group, "energy_bounds", wws->energy_bounds_); + + int32_t mesh_id; + read_dataset(ww_group, "mesh", mesh_id); + + if (model::mesh_map.count(mesh_id) == 0) { + fatal_error( + fmt::format("Mesh {} used in weight windows does not exist.", mesh_id)); + } + wws->set_mesh(model::mesh_map[mesh_id]); + + wws->lower_ww_ = xt::empty(wws->bounds_size()); + wws->upper_ww_ = xt::empty(wws->bounds_size()); + + read_dataset(ww_group, "lower_ww_bounds", wws->lower_ww_); + read_dataset(ww_group, "upper_ww_bounds", wws->upper_ww_); + read_dataset(ww_group, "survival_ratio", wws->survival_ratio_); + read_dataset(ww_group, "max_lower_bound_ratio", wws->max_lb_ratio_); + read_dataset(ww_group, "max_split", wws->max_split_); + read_dataset(ww_group, "weight_cutoff", wws->weight_cutoff_); + + close_group(ww_group); + + return wws; +} + +void WeightWindows::set_defaults() +{ + // ensure default values are set + 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]); } - // num spatial*energy bins must match num weight bins - int num_spatial_bins = this->mesh().n_bins(); - int num_energy_bins = energy_bounds_.size() - 1; - int num_weight_bins = lower_ww_.size(); - if (num_weight_bins != num_spatial_bins * num_energy_bins) { - auto err_msg = - fmt::format("In weight window domain {} the number of " - "energy/spatial bins ({}) does not match the number " - "of weight bins provided ({})", - id_, num_energy_bins * num_spatial_bins, num_weight_bins); - fatal_error(err_msg); + // 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); } } @@ -216,8 +301,40 @@ void WeightWindows::set_id(int32_t id) // Update ID and entry in the mesh map id_ = id; - variance_reduction::ww_map[id] = - variance_reduction::weight_windows.size() - 1; + variance_reduction::ww_map[id] = index_; +} + +void WeightWindows::set_energy_bounds(gsl::span bounds) +{ + energy_bounds_.clear(); + energy_bounds_.insert(energy_bounds_.begin(), bounds.begin(), bounds.end()); +} + +void WeightWindows::set_particle_type(ParticleType p_type) +{ + if (p_type != ParticleType::neutron && p_type != ParticleType::photon) + fatal_error( + fmt::format("Particle type '{}' cannot be applied to weight windows.", + particle_type_to_str(p_type))); + particle_type_ = p_type; +} + +void WeightWindows::set_mesh(int32_t mesh_idx) +{ + if (mesh_idx < 0 || mesh_idx >= model::meshes.size()) + fatal_error(fmt::format("Could not find a mesh for index {}", mesh_idx)); + + mesh_idx_ = mesh_idx; +} + +void WeightWindows::set_mesh(const std::unique_ptr& mesh) +{ + set_mesh(mesh.get()); +} + +void WeightWindows::set_mesh(const Mesh* mesh) +{ + set_mesh(model::mesh_map[mesh->id_]); } WeightWindow WeightWindows::get_weight_window(const Particle& p) const @@ -229,10 +346,10 @@ WeightWindow WeightWindows::get_weight_window(const Particle& p) const // Get mesh index for particle's position const auto& mesh = this->mesh(); - int ww_index = mesh.get_bin(p.r()); + int mesh_bin = mesh->get_bin(p.r()); // particle is outside the weight window mesh - if (ww_index < 0) + if (mesh_bin < 0) return {}; // particle energy @@ -246,13 +363,11 @@ WeightWindow WeightWindows::get_weight_window(const Particle& p) const int energy_bin = lower_bound_index(energy_bounds_.begin(), energy_bounds_.end(), E); - // indices now points to the correct weight for the given energy - ww_index += energy_bin * mesh.n_bins(); - + // mesh_bin += energy_bin * mesh->n_bins(); // Create individual weight window WeightWindow ww; - ww.lower_weight = lower_ww_[ww_index]; - ww.upper_weight = upper_ww_[ww_index]; + ww.lower_weight = lower_ww_(energy_bin, mesh_bin); + ww.upper_weight = upper_ww_(energy_bin, mesh_bin); ww.survival_weight = ww.lower_weight * survival_ratio_; ww.max_lb_ratio = max_lb_ratio_; ww.max_split = max_split_; @@ -260,22 +375,751 @@ WeightWindow WeightWindows::get_weight_window(const Particle& p) const return ww; } -void WeightWindows::to_hdf5(hid_t group) const +std::array WeightWindows::bounds_size() const { - hid_t ww_group = create_group(group, fmt::format("weight_windows {}", id_)); - - write_dataset( - ww_group, "particle_type", openmc::particle_type_to_str(particle_type_)); - write_dataset(ww_group, "energy_bounds", energy_bounds_); - write_dataset(ww_group, "lower_ww_bounds", lower_ww_); - write_dataset(ww_group, "upper_ww_bounds", upper_ww_); - write_dataset(ww_group, "survival_ratio", survival_ratio_); - write_dataset(ww_group, "max_lower_bound_ratio", max_lb_ratio_); - write_dataset(ww_group, "max_split", max_split_); - write_dataset(ww_group, "weight_cutoff", weight_cutoff_); - write_dataset(ww_group, "mesh", this->mesh().id_); - - close_group(ww_group); + int num_spatial_bins = this->mesh()->n_bins(); + int num_energy_bins = + energy_bounds_.size() > 0 ? energy_bounds_.size() - 1 : 1; + return {num_energy_bins, num_spatial_bins}; } +template +void WeightWindows::check_bounds(const T& lower, const T& upper) const +{ + // make sure that the upper and lower bounds have the same size + if (lower.size() != upper.size()) { + auto msg = fmt::format("The upper and lower weight window lengths do not " + "match.\n Lower size: {}\n Upper size: {}", + lower.size(), upper.size()); + fatal_error(msg); + } + this->check_bounds(lower); +} + +template +void WeightWindows::check_bounds(const T& bounds) const +{ + // check that the number of weight window entries is correct + auto dims = this->bounds_size(); + if (bounds.size() != dims[0] * dims[1]) { + auto err_msg = + fmt::format("In weight window domain {} the number of spatial " + "energy/spatial bins ({}) does not match the number " + "of weight bins ({})", + id_, dims, bounds.size()); + fatal_error(err_msg); + } +} + +void WeightWindows::set_bounds(const xt::xtensor& lower_bounds, + const xt::xtensor& upper_bounds) +{ + + this->check_bounds(lower_bounds, upper_bounds); + + // set new weight window values + lower_ww_ = lower_bounds; + upper_ww_ = upper_bounds; +} + +void WeightWindows::set_bounds( + const xt::xtensor& lower_bounds, double ratio) +{ + this->check_bounds(lower_bounds); + + // set new weight window values + lower_ww_ = lower_bounds; + upper_ww_ = lower_bounds; + upper_ww_ *= ratio; +} + +void WeightWindows::set_bounds( + gsl::span lower_bounds, gsl::span upper_bounds) +{ + check_bounds(lower_bounds, upper_bounds); + auto shape = this->bounds_size(); + lower_ww_ = xt::empty(shape); + upper_ww_ = xt::empty(shape); + + // set new weight window values + xt::view(lower_ww_, xt::all()) = + xt::adapt(lower_bounds.data(), lower_ww_.shape()); + xt::view(upper_ww_, xt::all()) = + xt::adapt(upper_bounds.data(), upper_ww_.shape()); +} + +void WeightWindows::set_bounds( + gsl::span lower_bounds, double ratio) +{ + this->check_bounds(lower_bounds); + + auto shape = this->bounds_size(); + lower_ww_ = xt::empty(shape); + upper_ww_ = xt::empty(shape); + + // set new weight window values + xt::view(lower_ww_, xt::all()) = + xt::adapt(lower_bounds.data(), lower_ww_.shape()); + xt::view(upper_ww_, xt::all()) = + xt::adapt(lower_bounds.data(), upper_ww_.shape()); + upper_ww_ *= ratio; +} + +void WeightWindows::update_magic( + const Tally* tally, const std::string& value, double threshold, double ratio) +{ + /////////////////////////// + // Setup and checks + /////////////////////////// + this->check_tally_update_compatibility(tally); + + lower_ww_.fill(-1); + upper_ww_.fill(-1); + + // determine which value to use + const std::set allowed_values = {"mean", "rel_err"}; + if (allowed_values.count(value) == 0) { + fatal_error(fmt::format("Invalid value '{}' specified for weight window " + "generation. Must be one of: 'mean' or 'rel_err'", + value)); + } + + // determine the index of the specified score + int score_index = tally->score_index("flux"); + if (score_index == C_NONE) { + fatal_error( + fmt::format("A 'flux' score required for weight window generation " + "is not present on tally {}.", + tally->id())); + } + + /////////////////////////// + // Extract tally data + // + // At the end of this section, the mean and rel_err array + // is a 2D view of tally data (n_e_groups, n_mesh_bins) + // + /////////////////////////// + + // build a shape for a view of the tally results, this will always be + // dimension 5 (3 filter dimensions, 1 score dimension, 1 results dimension) + std::array shape = { + 1, 1, 1, tally->n_scores(), static_cast(TallyResult::SIZE)}; + + // set the shape for the filters applied on the tally + for (int i = 0; i < tally->filters().size(); i++) { + const auto& filter = model::tally_filters[tally->filters(i)]; + shape[i] = filter->n_bins(); + } + + // build the transpose information to re-order data according to filter type + std::array transpose = {0, 1, 2, 3, 4}; + + // track our filter types and where we've added new ones + std::vector filter_types = tally->filter_types(); + + // assign other filter types to dummy positions if needed + if (!tally->has_filter(FilterType::PARTICLE)) + filter_types.push_back(FilterType::PARTICLE); + + if (!tally->has_filter(FilterType::ENERGY)) + filter_types.push_back(FilterType::ENERGY); + + // particle axis mapping + transpose[0] = + std::find(filter_types.begin(), filter_types.end(), FilterType::PARTICLE) - + filter_types.begin(); + + // energy axis mapping + transpose[1] = + std::find(filter_types.begin(), filter_types.end(), FilterType::ENERGY) - + filter_types.begin(); + + // mesh axis mapping + transpose[2] = + std::find(filter_types.begin(), filter_types.end(), FilterType::MESH) - + filter_types.begin(); + + // get a fully reshaped view of the tally according to tally ordering of + // filters + auto tally_values = xt::reshape_view(tally->results(), shape); + + // get a that is (particle, energy, mesh, scores, values) + auto transposed_view = xt::transpose(tally_values, transpose); + + // determine the dimension and index of the particle data + int particle_idx = 0; + if (tally->has_filter(FilterType::PARTICLE)) { + // get the particle filter + auto pf = tally->get_filter(); + const auto& particles = pf->particles(); + + // find the index of the particle that matches these weight windows + auto p_it = + std::find(particles.begin(), particles.end(), this->particle_type_); + // if the particle filter doesn't have particle data for the particle + // used on this weight windows instance, report an error + if (p_it == particles.end()) { + auto msg = fmt::format("Particle type '{}' not present on Filter {} for " + "Tally {} used to update WeightWindows {}", + particle_type_to_str(this->particle_type_), pf->id(), tally->id(), + this->id()); + fatal_error(msg); + } + + // use the index of the particle in the filter to down-select data later + particle_idx = p_it - particles.begin(); + } + + // down-select data based on particle and score + auto sum = xt::view(transposed_view, particle_idx, xt::all(), xt::all(), + score_index, static_cast(TallyResult::SUM)); + auto sum_sq = xt::view(transposed_view, particle_idx, xt::all(), xt::all(), + score_index, static_cast(TallyResult::SUM_SQ)); + int n = tally->n_realizations_; + + ////////////////////////////////////////////// + // + // Assign new weight windows + // + // Use references to the existing weight window data + // to store and update the values + // + ////////////////////////////////////////////// + + // up to this point the data arrays are views into the tally results (no + // computation has been performed) now we'll switch references to the tally's + // bounds to avoid allocating additional memory + auto& new_bounds = this->lower_ww_; + auto& rel_err = this->upper_ww_; + + // noalias avoids memory allocation here + xt::noalias(new_bounds) = sum / n; + + xt::noalias(rel_err) = + xt::sqrt(((sum_sq / n) - xt::square(new_bounds)) / (n - 1)) / new_bounds; + xt::filter(rel_err, sum <= 0.0).fill(INFTY); + + if (value == "rel_err") + xt::noalias(new_bounds) = 1 / rel_err; + + // get mesh volumes + auto mesh_vols = this->mesh()->volumes(); + + int e_bins = new_bounds.shape()[0]; + for (int e = 0; e < e_bins; e++) { + // select all + auto group_view = xt::view(new_bounds, e); + + double group_max = *std::max_element(group_view.begin(), group_view.end()); + // normalize values in this energy group by the maximum value for this + // group + if (group_max > 0.0) + group_view /= group_max; + + // divide by volume of mesh elements + for (int i = 0; i < group_view.size(); i++) { + group_view[i] /= mesh_vols[i]; + } + } + + // make sure that values where the mean is zero are set s.t. the weight window + // value will be ignored + xt::filter(new_bounds, sum <= 0.0).fill(-1.0); + + // make sure the weight windows are ignored for any locations where the + // relative error is higher than the specified relative error threshold + xt::filter(new_bounds, rel_err > threshold).fill(-1.0); + + // update the bounds of this weight window class + // noalias avoids additional memory allocation + xt::noalias(upper_ww_) = ratio * lower_ww_; +} + +void WeightWindows::check_tally_update_compatibility(const Tally* tally) +{ + // define the set of allowed filters for the tally + const std::set allowed_filters = { + FilterType::MESH, FilterType::ENERGY, FilterType::PARTICLE}; + + // retrieve a mapping of filter type to filter index for the tally + auto filter_indices = tally->filter_indices(); + + // a mesh filter is required for a tally used to update weight windows + if (!filter_indices.count(FilterType::MESH)) { + fatal_error( + "A mesh filter is required for a tally to update weight window bounds"); + } + + // ensure the mesh filter is using the same mesh as this weight window object + auto mesh_filter = tally->get_filter(); + + // make sure that all of the filters present on the tally are allowed + for (auto filter_pair : filter_indices) { + if (allowed_filters.find(filter_pair.first) == allowed_filters.end()) { + fatal_error(fmt::format("Invalid filter type '{}' found on tally " + "used for weight window generation.", + model::tally_filters[filter_pair.second]->type_str())); + } + } + + if (mesh_filter->mesh() != mesh_idx_) { + int32_t mesh_filter_id = model::meshes[mesh_filter->mesh()]->id(); + int32_t ww_mesh_id = model::meshes[this->mesh_idx_]->id(); + fatal_error(fmt::format("Mesh filter {} uses a different mesh ({}) than " + "weight window {} mesh ({})", + mesh_filter->id(), mesh_filter_id, id_, ww_mesh_id)); + } + + // if an energy filter exists, make sure the energy grid matches that of this + // weight window object + if (auto energy_filter = tally->get_filter()) { + std::vector filter_bins = energy_filter->bins(); + std::set filter_e_bounds( + energy_filter->bins().begin(), energy_filter->bins().end()); + if (filter_e_bounds.size() != energy_bounds().size()) { + fatal_error( + fmt::format("Energy filter {} does not have the same number of energy " + "bounds ({}) as weight window object {} ({})", + energy_filter->id(), filter_e_bounds.size(), id_, + energy_bounds().size())); + } + + for (auto e : energy_bounds()) { + if (filter_e_bounds.count(e) == 0) { + fatal_error(fmt::format( + "Energy bounds of filter {} and weight windows {} do not match", + energy_filter->id(), id_)); + } + } + } +} + + void WeightWindows::to_hdf5(hid_t group) const + { + hid_t ww_group = + create_group(group, fmt::format("weight_windows_{}", id())); + + write_dataset(ww_group, "mesh", this->mesh()->id()); + write_dataset( + ww_group, "particle_type", openmc::particle_type_to_str(particle_type_)); + write_dataset(ww_group, "energy_bounds", energy_bounds_); + write_dataset(ww_group, "lower_ww_bounds", lower_ww_); + write_dataset(ww_group, "upper_ww_bounds", upper_ww_); + write_dataset(ww_group, "survival_ratio", survival_ratio_); + write_dataset(ww_group, "max_lower_bound_ratio", max_lb_ratio_); + write_dataset(ww_group, "max_split", max_split_); + write_dataset(ww_group, "weight_cutoff", weight_cutoff_); + + close_group(ww_group); + } + + WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) + { + // read information from the XML node + int32_t mesh_id = std::stoi(get_node_value(node, "mesh")); + int32_t mesh_idx = model::mesh_map[mesh_id]; + max_realizations_ = std::stoi(get_node_value(node, "max_realizations")); + + int active_batches = settings::n_batches - settings::n_inactive; + if (max_realizations_ > active_batches) { + auto msg = fmt::format( + "The maximum number of specified tally realizations ({}) is " + "greater than the number of active batches ({}).", + max_realizations_, active_batches); + warning(msg); + } + auto tmp_str = get_node_value(node, "particle_type", true, true); + auto particle_type = str_to_particle_type(tmp_str); + + update_interval_ = std::stoi(get_node_value(node, "update_interval")); + on_the_fly_ = get_node_value_bool(node, "on_the_fly"); + + std::vector e_bounds; + if (check_for_node(node, "energy_bounds")) { + e_bounds = get_node_array(node, "energy_bounds"); + } else { + int p_type = static_cast(particle_type); + e_bounds.push_back(data::energy_min[p_type]); + 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") { + // parse non-default update parameters if specified + if (check_for_node(node, "update_parameters")) { + pugi::xml_node params_node = node.child("update_parameters"); + if (check_for_node(params_node, "value")) + tally_value_ = get_node_value(params_node, "value"); + if (check_for_node(params_node, "threshold")) + threshold_ = std::stod(get_node_value(params_node, "threshold")); + if (check_for_node(params_node, "ratio")) { + ratio_ = std::stod(get_node_value(params_node, "ratio")); + } + } + // check update parameter values + if (tally_value_ != "mean" && tally_value_ != "rel_err") { + fatal_error(fmt::format("Unsupported tally value '{}' specified for " + "weight window generation.", + tally_value_)); + } + if (threshold_ <= 0.0) + fatal_error( + fmt::format("Invalid relative error threshold '{}' (<= 0.0) " + "specified for weight window generation", + ratio_)); + if (ratio_ <= 1.0) + fatal_error(fmt::format("Invalid weight window ratio '{}' (<= 1.0) " + "specified for weight window generation")); + } else { + fatal_error(fmt::format( + "Unknown weight window update method '{}' specified", method_)); + } + + // create a matching weight windows object + auto wws = WeightWindows::create(); + ww_idx_ = wws->index(); + 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::update() const + { + const auto& wws = variance_reduction::weight_windows[ww_idx_]; + + Tally* tally = model::tallies[tally_idx_].get(); + + // if we're beyond the number of max realizations or not at the corrrect + // update interval, skip the update + if (max_realizations_ < tally->n_realizations_ || + tally->n_realizations_ % update_interval_ != 0) + return; + + wws->update_magic(tally, tally_value_, threshold_, ratio_); + + // if we're not doing on the fly generation, reset the tally results once + // we're done with the update + if (!on_the_fly_) + tally->reset(); + + // TODO: deactivate or remove tally once weight window generation is + // complete + } + + //============================================================================== + // C API + //============================================================================== + + int verify_ww_index(int32_t index) + { + if (index < 0 || index >= variance_reduction::weight_windows.size()) { + set_errmsg( + fmt::format("Index '{}' for weight windows is invalid", index)); + return OPENMC_E_OUT_OF_BOUNDS; + } + return 0; + } + + extern "C" int openmc_get_weight_windows_index(int32_t id, int32_t * idx) + { + auto it = variance_reduction::ww_map.find(id); + if (it == variance_reduction::ww_map.end()) { + set_errmsg(fmt::format("No weight windows exist with ID={}", id)); + return OPENMC_E_INVALID_ID; + } + + *idx = it->second; + return 0; + } + + extern "C" int openmc_weight_windows_get_id(int32_t index, int32_t * id) + { + if (int err = verify_ww_index(index)) + return err; + + const auto& wws = variance_reduction::weight_windows.at(index); + *id = wws->id(); + return 0; + } + + extern "C" int openmc_weight_windows_set_id(int32_t index, int32_t id) + { + if (int err = verify_ww_index(index)) + return err; + + const auto& wws = variance_reduction::weight_windows.at(index); + wws->set_id(id); + return 0; + } + + extern "C" int openmc_weight_windows_update_magic(int32_t ww_idx, + int32_t tally_idx, const char* value, double threshold, double ratio) + { + if (int err = verify_ww_index(ww_idx)) + return err; + + if (tally_idx < 0 || tally_idx >= model::tallies.size()) { + set_errmsg(fmt::format("Index '{}' for tally is invalid", tally_idx)); + return OPENMC_E_OUT_OF_BOUNDS; + } + + // get the requested tally + const Tally* tally = model::tallies.at(tally_idx).get(); + + // get the WeightWindows object + const auto& wws = variance_reduction::weight_windows.at(ww_idx); + + wws->update_magic(tally, value, threshold, ratio); + + return 0; + } + + extern "C" int openmc_weight_windows_set_mesh( + int32_t ww_idx, int32_t mesh_idx) + { + if (int err = verify_ww_index(ww_idx)) + return err; + const auto& wws = variance_reduction::weight_windows.at(ww_idx); + wws->set_mesh(mesh_idx); + return 0; + } + + extern "C" int openmc_weight_windows_get_mesh( + int32_t ww_idx, int32_t * mesh_idx) + { + if (int err = verify_ww_index(ww_idx)) + return err; + const auto& wws = variance_reduction::weight_windows.at(ww_idx); + *mesh_idx = model::mesh_map.at(wws->mesh()->id()); + return 0; + } + + extern "C" int openmc_weight_windows_set_energy_bounds( + int32_t ww_idx, double* e_bounds, size_t e_bounds_size) + { + if (int err = verify_ww_index(ww_idx)) + return err; + const auto& wws = variance_reduction::weight_windows.at(ww_idx); + wws->set_energy_bounds({e_bounds, e_bounds_size}); + return 0; + } + + extern "C" int openmc_weight_windows_get_energy_bounds( + int32_t ww_idx, const double** e_bounds, size_t* e_bounds_size) + { + if (int err = verify_ww_index(ww_idx)) + return err; + const auto& wws = variance_reduction::weight_windows[ww_idx].get(); + *e_bounds = wws->energy_bounds().data(); + *e_bounds_size = wws->energy_bounds().size(); + return 0; + } + + extern "C" int openmc_weight_windows_set_particle(int32_t index, int particle) + { + if (int err = verify_ww_index(index)) + return err; + + const auto& wws = variance_reduction::weight_windows.at(index); + wws->set_particle_type(static_cast(particle)); + return 0; + } + + extern "C" int openmc_weight_windows_get_particle( + int32_t index, int* particle) + { + if (int err = verify_ww_index(index)) + return err; + + const auto& wws = variance_reduction::weight_windows.at(index); + *particle = static_cast(wws->particle_type()); + return 0; + } + + extern "C" int openmc_weight_windows_get_bounds(int32_t index, + const double** lower_bounds, const double** upper_bounds, size_t* size) + { + if (int err = verify_ww_index(index)) + return err; + + const auto& wws = variance_reduction::weight_windows[index]; + *size = wws->lower_ww_bounds().size(); + *lower_bounds = wws->lower_ww_bounds().data(); + *upper_bounds = wws->upper_ww_bounds().data(); + return 0; + } + + extern "C" int openmc_weight_windows_set_bounds(int32_t index, + const double* lower_bounds, const double* upper_bounds, size_t size) + { + if (int err = verify_ww_index(index)) + return err; + + const auto& wws = variance_reduction::weight_windows[index]; + wws->set_bounds({lower_bounds, size}, {upper_bounds, size}); + return 0; + } + + extern "C" int openmc_extend_weight_windows( + int32_t n, int32_t * index_start, int32_t * index_end) + { + if (index_start) + *index_start = variance_reduction::weight_windows.size(); + if (index_end) + *index_end = variance_reduction::weight_windows.size() + n - 1; + for (int i = 0; i < n; ++i) + variance_reduction::weight_windows.push_back( + make_unique()); + return 0; + } + + extern "C" size_t openmc_weight_windows_size() + { + return variance_reduction::weight_windows.size(); + } + + extern "C" int openmc_weight_windows_export(const char* filename) + { + + if (!mpi::master) + return 0; + + std::string name = filename ? filename : "weight_windows.h5"; + + write_message(fmt::format("Exporting weight windows to {}...", name), 5); + + hid_t ww_file = file_open(name, 'w'); + + // Write file type + write_attribute(ww_file, "filetype", "weight_windows"); + + // Write revisiion number for state point file + write_attribute(ww_file, "version", VERSION_WEIGHT_WINDOWS); + + hid_t weight_windows_group = create_group(ww_file, "weight_windows"); + + hid_t mesh_group = create_group(ww_file, "meshes"); + + std::vector mesh_ids; + std::vector ww_ids; + for (const auto& ww : variance_reduction::weight_windows) { + + ww->to_hdf5(weight_windows_group); + ww_ids.push_back(ww->id()); + + // if the mesh has already been written, move on + int32_t mesh_id = ww->mesh()->id(); + if (std::find(mesh_ids.begin(), mesh_ids.end(), mesh_id) != + mesh_ids.end()) + continue; + + mesh_ids.push_back(mesh_id); + ww->mesh()->to_hdf5(mesh_group); + } + + write_attribute(mesh_group, "n_meshes", mesh_ids.size()); + write_attribute(mesh_group, "ids", mesh_ids); + close_group(mesh_group); + + write_attribute(weight_windows_group, "n_weight_windows", ww_ids.size()); + write_attribute(weight_windows_group, "ids", ww_ids); + close_group(weight_windows_group); + + file_close(ww_file); + + return 0; + } + + extern "C" int openmc_weight_windows_import(const char* filename) + { + std::string name = filename ? filename : "weight_windows.h5"; + + if (mpi::master) + write_message( + fmt::format("Importing weight windows from {}...", name), 5); + + if (!file_exists(name)) { + set_errmsg(fmt::format("File '{}' does not exist", name)); + } + + hid_t ww_file = file_open(name, 'r'); + + // Check that filetype is correct + std::string filetype; + read_attribute(ww_file, "filetype", filetype); + if (filetype != "weight_windows") { + file_close(ww_file); + set_errmsg(fmt::format("File '{}' is not a weight windows file.", name)); + return OPENMC_E_INVALID_ARGUMENT; + } + + // Check that the file version is compatible + std::array file_version; + read_attribute(ww_file, "version", file_version); + if (file_version[0] != VERSION_WEIGHT_WINDOWS[0]) { + std::string err_msg = + fmt::format("File '{}' has version {} which is incompatible with the " + "expected version ({}).", + name, file_version, VERSION_WEIGHT_WINDOWS); + set_errmsg(err_msg); + return OPENMC_E_INVALID_ARGUMENT; + } + + hid_t weight_windows_group = open_group(ww_file, "weight_windows"); + + std::vector names = group_names(weight_windows_group); + + for (const auto& name : names) { + WeightWindows::from_hdf5(weight_windows_group, name); + } + + close_group(weight_windows_group); + + file_close(ww_file); + + return 0; + } + } // namespace openmc diff --git a/tests/regression_tests/track_output/results_true.dat b/tests/regression_tests/track_output/results_true.dat index 37f4a0f98..746f1448f 100644 --- a/tests/regression_tests/track_output/results_true.dat +++ b/tests/regression_tests/track_output/results_true.dat @@ -1,4 +1,4 @@ -ParticleType.NEUTRON [((-9.663085e-01, -6.616522e-01, -9.863690e-01), (-7.513921e-01, 4.140970e-01, 5.137447e-01), 5.293873e+06, 0.000000e+00, 1.000000e+00, 23, 2403, 1) +neutron [((-9.663085e-01, -6.616522e-01, -9.863690e-01), (-7.513921e-01, 4.140970e-01, 5.137447e-01), 5.293873e+06, 0.000000e+00, 1.000000e+00, 23, 2403, 1) ((-1.281491e+00, -4.879529e-01, -7.708709e-01), (-7.513921e-01, 4.140970e-01, 5.137447e-01), 5.293873e+06, 1.323629e-10, 1.000000e+00, 22, 2403, 3) ((-1.368452e+00, -4.400285e-01, -7.114141e-01), (-7.513921e-01, 4.140970e-01, 5.137447e-01), 5.293873e+06, 1.688824e-10, 1.000000e+00, 21, 2403, 2) ((-2.150539e+00, -9.015151e-03, -1.766823e-01), (-7.513921e-01, 4.140970e-01, 5.137447e-01), 5.293873e+06, 4.973246e-10, 1.000000e+00, 22, 2403, 3) @@ -117,7 +117,7 @@ ParticleType.NEUTRON [((-9.663085e-01, -6.616522e-01, -9.863690e-01), (-7.513921 ((-1.682182e+01, 1.413107e+00, -2.343954e+01), (6.002458e-01, 4.967948e-01, 6.268172e-01), 8.602342e-03, 2.164421e-05, 1.000000e+00, 21, 2129, 2) ((-1.636364e+01, 1.792326e+00, -2.296107e+01), (9.641593e-01, -1.162836e-01, 2.384846e-01), 7.837617e-03, 2.759442e-05, 1.000000e+00, 21, 2129, 2) ((-1.598564e+01, 1.746738e+00, -2.286758e+01), (9.641593e-01, -1.162836e-01, 2.384846e-01), 7.837617e-03, 3.079609e-05, 0.000000e+00, 21, 2129, 2)] -ParticleType.NEUTRON [((-9.469716e-01, -2.580266e-01, 1.414357e-01), (1.438873e-01, 2.365819e-01, 9.608983e-01), 9.724461e+05, 0.000000e+00, 1.000000e+00, 23, 2403, 1) +neutron [((-9.469716e-01, -2.580266e-01, 1.414357e-01), (1.438873e-01, 2.365819e-01, 9.608983e-01), 9.724461e+05, 0.000000e+00, 1.000000e+00, 23, 2403, 1) ((-9.064818e-01, -1.914524e-01, 4.118326e-01), (9.231638e-01, 3.100833e-01, 2.271937e-01), 1.746594e+05, 2.064696e-10, 1.000000e+00, 23, 2403, 1) ((-8.178800e-01, -1.616918e-01, 4.336377e-01), (9.231638e-01, 3.100833e-01, 2.271937e-01), 1.746594e+05, 3.725262e-10, 1.000000e+00, 11, 1756, 1) ((-3.834875e-01, -1.578285e-02, 5.405432e-01), (-5.547498e-01, -6.358598e-02, 8.295839e-01), 1.474011e+05, 1.186660e-09, 1.000000e+00, 11, 1756, 1) @@ -143,7 +143,7 @@ ParticleType.NEUTRON [((-9.469716e-01, -2.580266e-01, 1.414357e-01), (1.438873e- ((-3.920090e+00, -9.211794e-01, 2.908406e+00), (9.094673e-01, -3.659267e-01, -1.974002e-01), 9.547085e+00, 1.260840e-07, 1.000000e+00, 23, 2387, 1) ((-3.729948e+00, -9.976834e-01, 2.867135e+00), (-2.916959e-01, -6.494657e-01, -7.022164e-01), 8.590591e+00, 1.750036e-07, 1.000000e+00, 23, 2387, 1) ((-3.744933e+00, -1.031047e+00, 2.831062e+00), (-2.916959e-01, -6.494657e-01, -7.022164e-01), 8.590591e+00, 1.876753e-07, 0.000000e+00, 23, 2387, 1)] -ParticleType.NEUTRON [((6.474155e+00, -5.192870e+00, 5.413003e+00), (-9.112559e-01, 3.950037e-01, 1.165536e-01), 2.052226e+06, 4.450859e-05, 1.000000e+00, 21, 2367, 2) +neutron [((6.474155e+00, -5.192870e+00, 5.413003e+00), (-9.112559e-01, 3.950037e-01, 1.165536e-01), 2.052226e+06, 4.450859e-05, 1.000000e+00, 21, 2367, 2) ((6.037250e+00, -5.003484e+00, 5.468885e+00), (-9.112559e-01, 3.950037e-01, 1.165536e-01), 2.052226e+06, 4.450883e-05, 1.000000e+00, 22, 2367, 3) ((5.942573e+00, -4.962444e+00, 5.480995e+00), (-9.112559e-01, 3.950037e-01, 1.165536e-01), 2.052226e+06, 4.450888e-05, 1.000000e+00, 23, 2367, 1) ((5.861800e+00, -4.927431e+00, 5.491326e+00), (-9.518772e-01, -3.064714e-01, -2.259335e-03), 1.141417e+06, 4.450892e-05, 1.000000e+00, 23, 2367, 1) diff --git a/tests/regression_tests/weightwindows/generators/__init__.py b/tests/regression_tests/weightwindows/generators/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/weightwindows/generators/test.py b/tests/regression_tests/weightwindows/generators/test.py new file mode 100644 index 000000000..f5bfc4dcd --- /dev/null +++ b/tests/regression_tests/weightwindows/generators/test.py @@ -0,0 +1,54 @@ +import os + +import numpy as np +import openmc +import pytest + + +def test_ww_generator(run_in_tmpdir): + # create a simple spherical shield model + model = openmc.Model() + + water = openmc.Material() + water.set_density('g/cc', 1.0) + water.add_nuclide('H1', 0.66) + water.add_nuclide('O16', 0.34) + + s = openmc.Sphere(r=50, boundary_type='vacuum') + c = openmc.Cell(fill=water, region=-s) + + model.geometry = openmc.Geometry([c]) + + model.settings.particles = 500 + model.settings.batches = 5 + model.settings.run_mode = 'fixed source' + model.settings.max_splits = 100 + + mesh = openmc.RegularMesh.from_domain(model.geometry.root_universe) + energy_bounds = np.linspace(0.0, 1e6, 70) + particle = 'neutron' + + wwg = openmc.WeightWindowGenerator(mesh, energy_bounds, particle) + 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') + + wws_mean = openmc.hdf5_to_wws() + assert len(wws_mean) == 1 + + # check that generation using the relative error works too + wwg.update_parameters['value'] = 'rel_err' + model.run() + + wws_rel_err = openmc.hdf5_to_wws() + assert len(wws_rel_err) == 1 + + # we should not get the same set of weight windows when switching to use of + # rel. err. + assert (wws_mean[0].lower_ww_bounds != wws_rel_err[0].lower_ww_bounds).any() diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index 1dfad56aa..dce75c740 100644 --- a/tests/unit_tests/weightwindows/test.py +++ b/tests/unit_tests/weightwindows/test.py @@ -6,6 +6,7 @@ import pytest from uncertainties import ufloat import openmc +import openmc.lib from openmc.stats import Discrete, Point from tests import cdtemp @@ -15,8 +16,6 @@ from tests import cdtemp def wws(): # weight windows - - ww_files = ('ww_n.txt', 'ww_p.txt') cwd = Path(__file__).parent.absolute() ww_n_file, ww_p_file = [cwd / Path(f) for f in ww_files] @@ -240,3 +239,47 @@ def test_roundtrip(run_in_tmpdir, model, wws): # ensure the lower bounds read in from the XML match those of the for ww_out, ww_in in zipped_wws: assert(ww_out == ww_in) + + +def test_ww_attrs(run_in_tmpdir, model): + model.export_to_xml() + + openmc.lib.init() + + tally = openmc.lib.tallies[model.tallies[0].id] + + wws = openmc.lib.WeightWindows.from_tally(tally) + + # this is the first weight window object created + assert wws.id == 1 + + with pytest.raises(ValueError): + tally.find_filter(openmc.lib.AzimuthalFilter) + + mesh_filter = tally.find_filter(openmc.lib.MeshFilter) + mesh = mesh_filter.mesh + + assert wws.mesh.id == mesh.id + + assert wws.particle == openmc.ParticleType.NEUTRON + + wws.particle = 1 + assert wws.particle == openmc.ParticleType.PHOTON + wws.particle = 'photon' + assert wws.particle == openmc.ParticleType.PHOTON + + with pytest.raises(ValueError): + wws.particle = '🌠' + + energy_filter = tally.find_filter(openmc.lib.EnergyFilter) + np.testing.assert_allclose(np.unique(energy_filter.bins), wws.energy_bounds) + + # at this point the weight window bounds are uninitialized + assert all(wws.bounds[0] == -1) + assert all(wws.bounds[1] == -1) + + wws = openmc.lib.WeightWindows.from_tally(tally, particle='photon') + assert wws.id == 2 + assert wws.particle == openmc.ParticleType.PHOTON + + openmc.lib.finalize() \ No newline at end of file diff --git a/tests/unit_tests/weightwindows/test_ww_gen.py b/tests/unit_tests/weightwindows/test_ww_gen.py new file mode 100644 index 000000000..74bbcfada --- /dev/null +++ b/tests/unit_tests/weightwindows/test_ww_gen.py @@ -0,0 +1,316 @@ +from itertools import permutations +from pathlib import Path + +import numpy as np +import openmc +import openmc.lib +import pytest + + +@pytest.fixture +def model(): + openmc.reset_auto_ids() + + # create a simple spherical shell shielding model + + ### Materials ### + water = openmc.Material() + water.set_density('g/cc', 1.0) + water.add_nuclide("H1", 2) + water.add_nuclide("O16", 1) + + steel = openmc.Material() + steel.set_density('g/cc', 8.0) + steel.add_nuclide("Fe56", 1.0) + + air = openmc.Material() + air.set_density('g/cc', 0.001205) + air.add_nuclide("N14", 0.781557629247) + air.add_nuclide("O16", 0.210668126508) + + boron = openmc.Material() + boron.set_density('g/cc', 2.52) + boron.add_nuclide("B10", 0.15856) + boron.add_nuclide("B11", 0.64144) + boron.add_nuclide("C0", 0.2) + + ### Geometry ### + radii = [5.0, 10.0, 30.0, 31.0, 50.0] + + surfs = [openmc.Sphere(r=r) for r in radii] + + surfs[-1].boundary_type = 'vacuum' + + regions = openmc.model.subdivide(surfs) + + mats = [air, water, steel, boron, air] + + cells = [openmc.Cell(fill=m, region=r) for r, m in zip(regions, mats)] + + geometry = openmc.Geometry(cells) + + ### Settings ### + + settings = openmc.Settings( + run_mode='fixed source', + particles=100, + batches=10, + max_splits=10, + survival_biasing=False + ) + + # 10 keV neutron point source at the origin + space = openmc.stats.Point() + energy = openmc.stats.Discrete(x=[1e4], p=[1.0]) + settings.source = openmc.Source(space=space, energy=energy) + + return openmc.Model(geometry=geometry, settings=settings) + + +# create a tally used for weight window generation +mesh = openmc.RegularMesh() +mesh.lower_left = [-50.0] * 3 +mesh.upper_right = [50.0] * 3 +# use un-equal mesh widths in each dimension to more robustly check +# use of tally data +mesh.dimension = (19, 20, 21) + +mf = openmc.MeshFilter(mesh) + +ef = openmc.EnergyFilter([0.0, 1e7]) + +pf = openmc.ParticleFilter(['neutron', 'photon']) + +filters = [mf, ef, pf] + +test_cases = list(permutations(filters)) +test_cases += list(permutations(filters[:-1])) +test_cases += list(permutations(filters[::2])) + + +def labels(params): + out = [] + for p in params: + if isinstance(p, openmc.ParticleFilter): + out.append('particle') + elif isinstance(p, openmc.MeshFilter): + out.append('mesh') + elif isinstance(p, openmc.EnergyFilter): + out.append('energy') + return "filters:" + '-'.join(out) + + +@pytest.mark.parametrize("filters", test_cases, ids=labels) +def test_ww_gen(filters, model): + + tally = openmc.Tally() + tally.filters = list(filters) + tally.scores = ['flux'] + model.tallies = openmc.Tallies([tally]) + + model.export_to_model_xml() + + ref_lower = None + ref_upper = None + # test weight window generation capability + with openmc.lib.run_in_memory(): + + # retrieve the tally we created above in memory + lib_tally = openmc.lib.tallies[tally.id] + + # create a new weight windows object + ww = openmc.lib.WeightWindows.from_tally(lib_tally) + + # run particle transport + openmc.lib.run() + + # capture analog data + analog_mean = np.copy(lib_tally.mean) + + # update the weight window values using tally results + ww.update_magic(lib_tally) + + assert any(ww.bounds[0] != -1) + assert any(ww.bounds[1] != -1) + + # make sure that the weight window update doesn't change tally values + np.testing.assert_equal(lib_tally.mean, analog_mean) + + # check against weight windows from the previous iteration + # the order of filters should not change the weight window values + if ref_lower is None: + ref_lower = ww.bounds[0].copy() + else: + np.testing.assert_equal(ref_lower, ww.bounds[0]) + + if ref_upper is None: + ref_upper = ww.bounds[1].copy() + else: + np.testing.assert_equal(ref_upper, ww.bounds[1]) + + # turn on weight windows for the subsequent run + openmc.lib.settings.weight_windows_on = True + + openmc.lib.hard_reset() + + openmc.lib.run() + + ww_mean = np.copy(lib_tally.mean) + + # we expect that the application of weight windows will populate more tally + # bins than the analog run for the meshes in this test model + assert any(ww_mean != analog_mean) + assert np.count_nonzero(ww_mean) > np.count_nonzero(analog_mean) + + +def test_ww_import_export(run_in_tmpdir, model): + # create a tally for weight windows + mesh = openmc.RegularMesh() + mesh.lower_left = [-50.0] * 3 + mesh.upper_right = [50.0] * 3 + # use un-equal mesh widths in each dimension to more robustly check + # use of tally data + mesh.dimension = (3, 4, 5) + + mf = openmc.MeshFilter(mesh) + + e_groups = np.logspace(0, 7, 8) + ef = openmc.EnergyFilter(e_groups) + + pf = openmc.ParticleFilter(['neutron', 'photon']) + + tally = openmc.Tally() + tally.filters = [mf, ef, pf] + tally.scores = ['flux'] + + model.tallies = openmc.Tallies([tally]) + + # first, generate some weight windows + model.export_to_model_xml() + + openmc.lib.init() + + tally = openmc.lib.tallies[tally.id] + + ww = openmc.lib.WeightWindows.from_tally(tally) + + openmc.lib.run() + + mean_before = np.array(tally.mean) + + ww.update_magic(tally) + + mean_after = np.array(tally.mean) + + assert (mean_before == mean_after).all() + + lb_before, up_before = ww.bounds + + openmc.lib.export_weight_windows() + + assert Path('weight_windows.h5').exists() + + openmc.lib.import_weight_windows('weight_windows.h5') + + ww = openmc.lib.weight_windows[2] + + lb_after, up_after = ww.bounds + + assert np.allclose(lb_before, lb_after) + assert np.allclose(up_before, up_after) + + openmc.lib.finalize() + + +def test_ww_gen_roundtrip(run_in_tmpdir, model): + + mesh = openmc.RegularMesh.from_domain(model.geometry.root_universe) + energy_bounds = np.linspace(0.0, 1e8, 11) + particle_type = 'neutron' + + wwg = openmc.WeightWindowGenerator(mesh, energy_bounds, particle_type) + wwg.update_params = {'ratio' : 5.0, + 'threshold': 0.8, + 'value' : 'mean'} + + model.settings.weight_window_generators = wwg + model.export_to_xml() + + model_in = openmc.Model.from_xml() + + assert len(model_in.settings.weight_window_generators) == 1 + + wwg_in = model.settings.weight_window_generators[0] + + # rountrip tests + model_in = openmc.Model.from_xml() + assert len(model_in.settings.weight_window_generators) == 1 + wwg_in = model_in.settings.weight_window_generators[0] + assert wwg_in.max_realizations == 1 + assert wwg_in.on_the_fly == True + assert wwg_in.update_interval == 1 + assert wwg_in.update_parameters == wwg.update_parameters + + with pytest.raises(ValueError): + wwg.method = '🦍🐒' + + with pytest.raises(TypeError): + wwg.update_parameters = {'ratio' : 'one-to-one'} + + with pytest.raises(ValueError): + wwg.max_realizations = -1 + + +def test_python_hdf5_roundtrip(run_in_tmpdir, model): + + # add a tally to the model + mesh = openmc.RegularMesh.from_domain(model.geometry) + + # some arbitrary energy groups + e_groups = np.logspace(0, 6, 4) + energy_filter = openmc.EnergyFilter(e_groups) + + bounds = np.arange(energy_filter.num_bins * np.prod(mesh.dimension)) + + wws = openmc.WeightWindows(mesh, bounds, bounds, energy_bounds=e_groups) + + model.settings.weight_windows = [wws] + + model.export_to_xml() + + # initialize and export wws to HDF5 + openmc.lib.init() + + openmc.lib.export_weight_windows() + + openmc.lib.finalize() + + wws_hdf5 = openmc.hdf5_to_wws()[0] + + # ensure + assert all(wws.energy_bounds == wws_hdf5.energy_bounds) + assert wws.id == wws_hdf5.id + assert wws.mesh.id == wws_hdf5.mesh.id + np.testing.assert_array_equal(wws.lower_ww_bounds, wws_hdf5.lower_ww_bounds) + np.testing.assert_array_equal(wws.upper_ww_bounds, wws_hdf5.upper_ww_bounds) + + +def test_ww_bounds_set_in_memory(run_in_tmpdir, model): + tally = openmc.Tally() + tally.filters = filters + tally.scores = ['flux'] + model.tallies = [tally] + + bounds = np.arange(ef.num_bins * np.prod(mf.mesh.dimension)) + + model.export_to_xml() + + openmc.lib.init() + + lib_tally = openmc.lib.tallies[tally.id] + + wws = openmc.lib.WeightWindows.from_tally(lib_tally) + + wws.bounds = (bounds, bounds) + + openmc.lib.finalize()