Use RegularMesh* for entropy and ufs

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

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 <cstdint>
@@ -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<double> k_generation;
 extern std::vector<int64_t> work_index;
 

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'

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<const RegularMesh*>(
-    model::meshes[settings::index_entropy_mesh].get());
-
   // Get source weight in each mesh bin
   bool sites_outside;
-  xt::xtensor<double, 1> p = m.count_sites(simulation::fission_bank,
-    &sites_outside);
+  xt::xtensor<double, 1> 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<const RegularMesh*>(
-    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<const RegularMesh*>(
-    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;
   }

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;

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<const RegularMesh*>(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<const RegularMesh*>(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

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 <entropy> element "
@@ -500,25 +498,28 @@ void read_settings_xml()
     model::meshes.push_back(std::make_unique<RegularMesh>(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<RegularMesh*>(model::meshes[index_entropy_mesh].get());
-    if (m.shape_.dimension() == 0) {
+  if (i_entropy_mesh >= 0) {
+    auto* m = dynamic_cast<RegularMesh*>(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 <uniform_fs> element "
@@ -546,14 +548,18 @@ void read_settings_xml()
     model::meshes.push_back(std::make_unique<RegularMesh>(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<RegularMesh*>(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;
   }

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<double> k_generation;
 std::vector<int64_t> work_index;