Updates to the LibMesh methods for adding variables and writing.

CMakeLists.txt

@@ -395,7 +395,7 @@ if (MPI_ENABLED)
   target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI)
 endif()
 if (libmesh)
-  target_compile_definitions(libopenmc PUBLIC -DLIBMESH)
+  target_compile_definitions(libopenmc PRIVATE LIBMESH)
 endif()
 
 # Set git SHA1 hash as a compile definition

include/openmc/mesh.h

@@ -282,6 +282,16 @@ public:
 
   std::string bin_label(int bin) const override;
 
+  //! Add a variable to the libmesh mesh instance
+  virtual void add_score(const std::string& var_name) = 0;
+
+  //! Set the value of a bin for a variable on the libmesh mesh instance
+  virtual void set_score_data(const std::string& var_name,
+                              std::vector<double> values,
+                              std::vector<double> std_dev) const = 0;
+
+  virtual void write(std::string filename) const = 0;
+
   std::string filename_;
 };
 
@@ -495,15 +505,17 @@ public:
 
   bool intersects(Position& r0, Position r1, int* ijk) const;
 
-  void write(const std::string& filename) const;
+  void write(std::string filename) const override;
 
   void to_hdf5(hid_t group) const;
 
   //! Add a variable to the libmesh mesh instance
-  void add_variable(const std::string& var_name);
+  void add_score(const std::string& var_name) override;
 
   //! Set the value of a bin for a variable on the libmesh mesh instance
-  void set_variable(const std::string& var_name, int bin, double value);
+  void set_score_data(const std::string& var_name,
+                      std::vector<double> values,
+                      std::vector<double> std_dev) const override;
 
 private:
 

include/openmc/state_point.h

@@ -19,7 +19,7 @@ void read_source_bank(hid_t group_id, std::vector<Particle::Bank>& sites, bool d
 void write_tally_results_nr(hid_t file_id);
 void restart_set_keff();
 
-#ifdef DAGMC
+#ifdef LIBMESH
 void write_unstructured_mesh_results();
 #endif
 

src/mesh.cpp

@@ -2277,30 +2277,60 @@ LibMesh::get_element_from_bin(int bin) const {
 }
 
 void
-LibMesh::add_variable(const std::string& var_name) {
+LibMesh::add_score(const std::string& var_name) {
   // check if this is a new varaible
-  if (!variable_map_.count(var_name)) {
+  std::string value_name = var_name + "_value";
+  if (!variable_map_.count(value_name)) {
     auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
-    auto var_num = eqn_sys.add_variable(var_name, libMesh::CONSTANT, libMesh::MONOMIAL);
-    variable_map_[var_name] = var_num;
+    auto var_num = eqn_sys.add_variable(value_name, libMesh::CONSTANT, libMesh::MONOMIAL);
+    variable_map_[value_name] = var_num;
+    equation_systems_->init();
+  }
+
+  std::string std_dev_name = var_name + "_std_dev";
+  // check if this is a new varaible
+  if (!variable_map_.count(std_dev_name)) {
+    auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
+    auto var_num = eqn_sys.add_variable(std_dev_name, libMesh::CONSTANT, libMesh::MONOMIAL);
+    variable_map_[std_dev_name] = var_num;
     equation_systems_->init();
   }
 }
 
 void
-LibMesh::set_variable(const std::string& var_name, int bin, double value) {
-  // look up the variable
-  unsigned int var_num = variable_map_.at(var_name);
-  auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
-  const libMesh::DofMap& dof_map = eqn_sys.get_dof_map();
-  auto e = m_->elem_ptr(bin);
-  std::vector<libMesh::dof_id_type> dof_indices;
-  dof_map.dof_indices(e, dof_indices, var_num);
-  Ensures(dof_indices.size() == 1);
-  eqn_sys.solution->set(dof_indices[0], value);
+LibMesh::set_score_data(const std::string& var_name,
+                      std::vector<double> values,
+                      std::vector<double> std_dev) const
+{
+    auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
+    const libMesh::DofMap& dof_map = eqn_sys.get_dof_map();
+
+    // look up the value variable
+    std::string value_name = var_name + "_value";
+    unsigned int value_num = variable_map_.at(value_name);
+    // look up the std dev variable
+    std::string std_dev_name = var_name + "_std_dev";
+    unsigned int std_dev_num = variable_map_.at(std_dev_name);
+
+  for (int bin = 0; bin < this->n_bins(); bin++) {
+    auto e = m_->elem_ptr(bin);
+
+    // set value
+    std::vector<libMesh::dof_id_type> value_dof_indices;
+    dof_map.dof_indices(e, value_dof_indices, value_num);
+    Ensures(value_dof_indices.size() == 1);
+    eqn_sys.solution->set(value_dof_indices[0], values[bin]);
+
+    // set std dev
+    std::vector<libMesh::dof_id_type> std_dev_dof_indices;
+    dof_map.dof_indices(e, std_dev_dof_indices, std_dev_num);
+    Ensures(std_dev_dof_indices.size() == 1);
+    eqn_sys.solution->set(std_dev_dof_indices[0], std_dev[bin]);
+
+  }
 }
 
-void LibMesh::write(const std::string& filename) const {
+void LibMesh::write(std::string filename) const {
   libMesh::ExodusII_IO exo(*m_);
   std::set<std::string> systems_out = {eq_system_name_};
   exo.write_discontinuous_exodusII(filename, *equation_systems_, &systems_out);

src/state_point.cpp

@@ -770,31 +770,27 @@ void read_source_bank(hid_t group_id, std::vector<Particle::Bank>& sites, bool d
   H5Tclose(banktype);
 }
 
-#ifdef DAGMC
+#ifdef LIBMESH
 void write_unstructured_mesh_results() {
-
+  std::cout << "Checking for unstructured meshes to write..." << std::endl;
   for (auto& tally : model::tallies) {
 
     std::vector<std::string> tally_scores;
     for (auto filter_idx : tally->filters()) {
       auto& filter = model::tally_filters[filter_idx];
-      if (filter->type() != "mesh") continue;
-
-      // check if the filter uses an unstructured mesh
-      auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
-      auto mesh_idx = mesh_filter->mesh();
-      auto umesh = dynamic_cast<UnstructuredMesh*>(model::meshes[mesh_idx].get());
-
-      if (!umesh) continue;
-
-      // if this tally has more than one filter, print
-      // warning and skip writing the mesh
-      if (tally->filters().size() > 1) {
-        warning(fmt::format("Skipping unstructured mesh writing for tally "
-                            "{}. More than one filter is present on the tally.",
-                            tally->id_));
-        break;
-      }
+      if (filter->type() == "mesh") {
+        auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
+        auto& mesh = model::meshes[mesh_filter->mesh()];
+        auto umesh = dynamic_cast<UnstructuredMeshBase*>(mesh.get());
+        if (umesh) {
+          for (int i = 0; i < tally->scores_.size(); i++) {
+            // get values for this score and create vectors for marking
+            // up the mesh
+            std::vector<double> vals_vec, sum_sq_vec;
+            for (int j = 0; j < tally->results_.shape()[0]; j++) {
+              vals_vec.push_back(tally->results_(j, i, TallyResult::SUM) / tally->n_realizations_);
+              sum_sq_vec.push_back(tally->results_(j , i, TallyResult::SUM_SQ) - std::pow(vals_vec[i], 2)/ tally->n_realizations_);
+            }
 
       int n_realizations = tally->n_realizations_;