diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 4b76c3dc5..77127db0a 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -59,9 +59,6 @@ extern std::string path_source; extern std::string path_sourcepoint; //!< path to a source file extern std::string path_statepoint; //!< path to a statepoint file -extern int32_t index_entropy_mesh; //!< Index of entropy mesh in global mesh array -extern int32_t index_ufs_mesh; //!< Index of UFS mesh in global mesh array - extern "C" int32_t n_batches; //!< number of (inactive+active) batches extern "C" int32_t n_inactive; //!< number of inactive batches extern "C" int32_t gen_per_batch; //!< number of generations per batch diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index fef6ebbcd..8efac4c7f 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -4,6 +4,7 @@ #ifndef OPENMC_SIMULATION_H #define OPENMC_SIMULATION_H +#include "openmc/mesh.h" #include "openmc/particle.h" #include @@ -39,6 +40,9 @@ extern "C" int total_gen; //!< total number of generations simulated extern double total_weight; //!< Total source weight in a batch extern int64_t work_per_rank; //!< number of particles per MPI rank +extern const RegularMesh* entropy_mesh; +extern const RegularMesh* ufs_mesh; + extern std::vector k_generation; extern std::vector work_index; diff --git a/openmc/filter.py b/openmc/filter.py index 4d14f054e..fa1085b97 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -611,7 +611,9 @@ class MeshFilter(Filter): self.id = filter_id def __hash__(self): - return hash(repr(self)) + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id) + return hash(string) def __repr__(self): string = type(self).__name__ + '\n' diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 0f4b49f46..310279d36 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -532,14 +532,10 @@ int openmc_get_keff(double* k_combined) void shannon_entropy() { - // Get reference to entropy mesh - auto& m = *dynamic_cast( - model::meshes[settings::index_entropy_mesh].get()); - // Get source weight in each mesh bin bool sites_outside; - xt::xtensor p = m.count_sites(simulation::fission_bank, - &sites_outside); + xt::xtensor p = simulation::entropy_mesh->count_sites( + simulation::fission_bank, &sites_outside); // display warning message if there were sites outside entropy box if (sites_outside) { @@ -565,22 +561,19 @@ void shannon_entropy() void ufs_count_sites() { - auto& m = *dynamic_cast( - model::meshes[settings::index_entropy_mesh].get()); - if (simulation::current_batch == 1 && simulation::current_gen == 1) { // On the first generation, just assume that the source is already evenly // distributed so that effectively the production of fission sites is not // biased auto s = xt::view(simulation::source_frac, xt::all()); - s = m.volume_frac_; + s = simulation::ufs_mesh->volume_frac_; } else { // count number of source sites in each ufs mesh cell bool sites_outside; - simulation::source_frac = m.count_sites(simulation::source_bank, - &sites_outside); + simulation::source_frac = simulation::ufs_mesh->count_sites( + simulation::source_bank, &sites_outside); // Check for sites outside of the mesh if (mpi::master && sites_outside) { @@ -589,7 +582,7 @@ void ufs_count_sites() #ifdef OPENMC_MPI // Send source fraction to all processors - int n_bins = xt::prod(m.shape_)(); + int n_bins = xt::prod(simulation::ufs_mesh->shape_)(); MPI_Bcast(simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm); #endif @@ -607,18 +600,16 @@ void ufs_count_sites() double ufs_get_weight(const Particle* p) { - auto& m = *dynamic_cast( - model::meshes[settings::index_entropy_mesh].get()); - // Determine indices on ufs mesh for current location - int mesh_bin = m.get_bin(p->r()); + int mesh_bin = simulation::ufs_mesh->get_bin(p->r()); if (mesh_bin < 0) { p->write_restart(); fatal_error("Source site outside UFS mesh!"); } if (simulation::source_frac(mesh_bin) != 0.0) { - return m.volume_frac_ / simulation::source_frac(mesh_bin); + return simulation::ufs_mesh->volume_frac_ + / simulation::source_frac(mesh_bin); } else { return 1.0; } diff --git a/src/finalize.cpp b/src/finalize.cpp index 34c4b9781..fa6b7628c 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -71,8 +71,6 @@ int openmc_finalize() settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0}; settings::entropy_on = false; settings::gen_per_batch = 1; - settings::index_entropy_mesh = -1; - settings::index_ufs_mesh = -1; settings::legendre_to_tabular = true; settings::legendre_to_tabular_points = -1; settings::n_particles = -1; @@ -114,6 +112,9 @@ int openmc_finalize() simulation::satisfy_triggers = false; simulation::total_gen = 0; + simulation::entropy_mesh = nullptr; + simulation::ufs_mesh = nullptr; + data::energy_max = {INFTY, INFTY}; data::energy_min = {0.0, 0.0}; model::root_universe = -1; diff --git a/src/plot.cpp b/src/plot.cpp index 070726bcb..43e7ed69c 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -19,6 +19,7 @@ #include "openmc/progress_bar.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/string_utils.h" namespace openmc { @@ -497,20 +498,38 @@ Plot::set_meshlines(pugi::xml_node plot_node) // Set mesh based on type if ("ufs" == meshtype) { - if (settings::index_ufs_mesh < 0) { + if (!simulation::ufs_mesh) { std::stringstream err_msg; err_msg << "No UFS mesh for meshlines on plot " << id_; fatal_error(err_msg); } else { - index_meshlines_mesh_ = settings::index_ufs_mesh; + for (int i = 0; i < model::meshes.size(); ++i) { + if (const auto* m + = dynamic_cast(model::meshes[i].get())) { + if (m == simulation::ufs_mesh) { + index_meshlines_mesh_ = i; + } + } + } + if (index_meshlines_mesh_ == -1) + fatal_error("Could not find the UFS mesh for meshlines plot"); } } else if ("entropy" == meshtype) { - if (settings::index_entropy_mesh < 0) { + if (!simulation::entropy_mesh) { std::stringstream err_msg; - err_msg <<"No entropy mesh for meshlines on plot " << id_; + err_msg << "No entropy mesh for meshlines on plot " << id_; fatal_error(err_msg); } else { - index_meshlines_mesh_ = settings::index_entropy_mesh; + for (int i = 0; i < model::meshes.size(); ++i) { + if (const auto* m + = dynamic_cast(model::meshes[i].get())) { + if (m == simulation::entropy_mesh) { + index_meshlines_mesh_ = i; + } + } + } + if (index_meshlines_mesh_ == -1) + fatal_error("Could not find the entropy mesh for meshlines plot"); } } else if ("tally" == meshtype) { // Ensure that there is a mesh id if the type is tally diff --git a/src/settings.cpp b/src/settings.cpp index 27f8c31d2..63827c3f4 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -73,9 +73,6 @@ std::string path_source; std::string path_sourcepoint; std::string path_statepoint; -int32_t index_entropy_mesh {-1}; -int32_t index_ufs_mesh {-1}; - int32_t n_batches; int32_t n_inactive {0}; int32_t gen_per_batch {1}; @@ -481,6 +478,7 @@ void read_settings_xml() read_meshes(root); // Shannon Entropy mesh + int32_t i_entropy_mesh = -1; if (check_for_node(root, "entropy_mesh")) { int temp = std::stoi(get_node_value(root, "entropy_mesh")); if (model::mesh_map.find(temp) == model::mesh_map.end()) { @@ -488,7 +486,7 @@ void read_settings_xml() msg << "Mesh " << temp << " specified for Shannon entropy does not exist."; fatal_error(msg); } - index_entropy_mesh = model::mesh_map.at(temp); + i_entropy_mesh = model::mesh_map.at(temp); } else if (check_for_node(root, "entropy")) { warning("Specifying a Shannon entropy mesh via the element " @@ -500,25 +498,28 @@ void read_settings_xml() model::meshes.push_back(std::make_unique(node_entropy)); // Set entropy mesh index - index_entropy_mesh = model::meshes.size() - 1; + i_entropy_mesh = model::meshes.size() - 1; // Assign ID and set mapping model::meshes.back()->id_ = 10000; - model::mesh_map[10000] = index_entropy_mesh; + model::mesh_map[10000] = i_entropy_mesh; } - if (index_entropy_mesh >= 0) { - auto& m = *dynamic_cast(model::meshes[index_entropy_mesh].get()); - if (m.shape_.dimension() == 0) { + if (i_entropy_mesh >= 0) { + auto* m = dynamic_cast(model::meshes[i_entropy_mesh].get()); + if (!m) fatal_error("Only regular meshes can be used as an entropy mesh"); + simulation::entropy_mesh = m; + + if (m->shape_.dimension() == 0) { // If the user did not specify how many mesh cells are to be used in // each direction, we automatically determine an appropriate number of // cells int n = std::ceil(std::pow(n_particles / 20.0, 1.0/3.0)); - m.shape_ = {n, n, n}; - m.n_dimension_ = 3; + m->shape_ = {n, n, n}; + m->n_dimension_ = 3; // Calculate width - m.width_ = (m.upper_right_ - m.lower_left_) / m.shape_; + m->width_ = (m->upper_right_ - m->lower_left_) / m->shape_; } // Turn on Shannon entropy calculation @@ -526,6 +527,7 @@ void read_settings_xml() } // Uniform fission source weighting mesh + int32_t i_ufs_mesh = -1; if (check_for_node(root, "ufs_mesh")) { auto temp = std::stoi(get_node_value(root, "ufs_mesh")); if (model::mesh_map.find(temp) == model::mesh_map.end()) { @@ -534,7 +536,7 @@ void read_settings_xml() "does not exist."; fatal_error(msg); } - index_ufs_mesh = model::mesh_map.at(temp); + i_ufs_mesh = model::mesh_map.at(temp); } else if (check_for_node(root, "uniform_fs")) { warning("Specifying a UFS mesh via the element " @@ -546,14 +548,18 @@ void read_settings_xml() model::meshes.push_back(std::make_unique(node_ufs)); // Set entropy mesh index - index_ufs_mesh = model::meshes.size() - 1; + i_ufs_mesh = model::meshes.size() - 1; // Assign ID and set mapping model::meshes.back()->id_ = 10001; - model::mesh_map[10001] = index_entropy_mesh; + model::mesh_map[10001] = i_entropy_mesh; } - if (index_ufs_mesh >= 0) { + if (i_ufs_mesh >= 0) { + auto* m = dynamic_cast(model::meshes[i_ufs_mesh].get()); + if (!m) fatal_error("Only regular meshes can be used as a UFS mesh"); + simulation::ufs_mesh = m; + // Turn on uniform fission source weighting ufs_on = true; } diff --git a/src/simulation.cpp b/src/simulation.cpp index fc7a78851..f94ba232a 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -247,6 +247,9 @@ int total_gen {0}; double total_weight; int64_t work_per_rank; +const RegularMesh* entropy_mesh {nullptr}; +const RegularMesh* ufs_mesh {nullptr}; + std::vector k_generation; std::vector work_index;