R2S tally fix-- correcting number of energy groups

CMakeLists.txt

@@ -1,5 +1,5 @@
 project(DAGMC)
-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 3.1)
 enable_language(Fortran)
 
 # Set DAGMC version

src/mcnp/meshtal_funcs.cpp

@@ -139,7 +139,11 @@ void dagmc_fmesh_setup_mesh_(int* fm_ipt, int* id, int* fmesh_idx,
 
   // Copy emesh bin boundaries from MCNP (includes 0.0 MeV)
   std::vector<double> energy_boundaries;
-  for (int i = 0; i < *n_energy_mesh; ++i) {
+  // if there is more than 1 bin, dont want 0->bottom bin
+  int bottom = 0;
+  if (*n_energy_mesh > 2)
+    bottom = 1;
+  for (int i = bottom ; i < *n_energy_mesh; ++i) {
     energy_boundaries.push_back(energy_mesh[i]);
   }
 

src/tally/CMakeLists.txt

@@ -6,6 +6,8 @@ file(GLOB PUB_HEADERS "*.hpp")
 set(LINK_LIBS lapack dagmc)
 set(LINK_LIBS_EXTERN_NAMES)
 
+set(CMAKE_CXX_STANDARD 11)
+
 dagmc_install_library(dagtally)
 
 if (BUILD_TESTS)

src/tally/KDEMeshTally.cpp

@@ -183,7 +183,15 @@ void KDEMeshTally::write_data(double num_histories) {
     moab::EntityHandle point = *i;
     unsigned int point_index = get_entity_index(point);
 
-    for (unsigned int j = 0; j < data->get_num_energy_bins(); ++ j) {
+    unsigned int num_ebins = data->get_num_energy_bins();
+    // if there is a total, dont do anything with it
+    if (data->has_total_energy_bin())
+      num_ebins--;
+
+    double tally_vect[num_ebins];
+    double error_vect[num_ebins];
+
+    for (unsigned int j = 0; j < num_ebins; ++ j) {
       std::pair <double, double> tally_data = data->get_data(point_index, j);
       double tally = tally_data.first;
       double error = tally_data.second;
@@ -198,15 +206,14 @@ void KDEMeshTally::write_data(double num_histories) {
       // normalize mesh tally result by the number of source particles
       tally /= num_histories;
 
-      // set tally and error tag values for this entity
-      rval = mbi->tag_set_data(tally_tags[j], &point, 1, &tally);
-
-      assert(moab::MB_SUCCESS == rval);
-
-      rval = mbi->tag_set_data(error_tags[j], &point, 1, &rel_error);
-
-      assert(moab::MB_SUCCESS == rval);
+      tally_vect[j] = tally;
+      error_vect[j] = rel_error;
     }
+    // set tally and error tag values for this entity
+    rval = mbi->tag_set_data(tally_tag, &point, 1, tally_vect);
+    MB_CHK_SET_ERR_RET(rval, "Failed to set the tally_tag data");
+    rval = mbi->tag_set_data(error_tag, &point, 1, error_vect);
+    MB_CHK_SET_ERR_RET(rval, "Failed to set the error_tag data");
   }
 
   // create a global tag to store the bandwidth value
@@ -225,16 +232,16 @@ void KDEMeshTally::write_data(double num_histories) {
   assert(moab::MB_SUCCESS == rval);
 
   // define list of tags to include and write mesh to output file
-  std::vector<moab::Tag> output_tags = tally_tags;
-  output_tags.insert(output_tags.end(), error_tags.begin(), error_tags.end());
+  std::vector<moab::Tag> output_tags;
+  output_tags.push_back(tally_tag);
+  output_tags.push_back(error_tag);
   output_tags.push_back(bandwidth_tag);
 
   rval = mbi->write_file(output_filename.c_str(),
                          NULL, NULL,
                          &tally_mesh_set, 1,
                          &(output_tags[0]),
                          output_tags.size());
-
   assert(moab::MB_SUCCESS == rval);
 }
 //---------------------------------------------------------------------------//

src/tally/MeshTally.cpp

@@ -113,40 +113,44 @@ moab::ErrorCode MeshTally::setup_tags(moab::Interface* mbi, const char* prefix)
   int tag_size = 1;
 
   unsigned int num_bins = data->get_num_energy_bins();
-  tally_tags.resize(num_bins);
-  error_tags.resize(num_bins);
+  if (data->has_total_energy_bin())
+    num_bins--;
+
+  std::string tag_name = pfx + "TALLY_TAG";
+  std::string error_tag_name = pfx + "ERROR_TAG";
+
+  // create vector tag to score particle spectra
+  rval = mbi->tag_get_handle(tag_name.c_str(),
+                             num_bins,
+                             moab::MB_TYPE_DOUBLE,
+                             tally_tag,
+                             moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
+  MB_CHK_SET_ERR(rval, "Failed to get the tag handle");
+
+  rval = mbi->tag_get_handle(error_tag_name.c_str(),
+                             num_bins,
+                             moab::MB_TYPE_DOUBLE,
+                             error_tag,
+                             moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
+  MB_CHK_SET_ERR(rval, "Failed to get the tag handle");
+
+  // create single tag to store the total
+  std::string total_tally_tag_name = tag_name + "_TOTAL";
+  rval = mbi->tag_get_handle(total_tally_tag_name.c_str(),
+                             tag_size,
+                             moab::MB_TYPE_DOUBLE,
+                             total_tally_tag,
+                             moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
+  MB_CHK_SET_ERR(rval, "Failed to get the tag handle");
+
+  std::string total_error_tag_name = error_tag_name + "_TOTAL";
+  rval = mbi->tag_get_handle(total_error_tag_name.c_str(),
+                             tag_size,
+                             moab::MB_TYPE_DOUBLE,
+                             total_error_tag,
+                             moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
+  MB_CHK_SET_ERR(rval, "Failed to get the tag handle");
 
-  // Create separate MOAB tag handles for every energy bin
-  for (unsigned i = 0; i < num_bins; ++i) {
-    std::string t_name = pfx + "TALLY_TAG", e_name = pfx + "ERROR_TAG";
-    std::stringstream str;
-
-    if (i + 1 != num_bins) {
-      str << "_" << input_data.energy_bin_bounds[i]
-          << '-' << input_data.energy_bin_bounds[i + 1];
-    }
-
-    t_name += str.str();
-    e_name += str.str();
-
-    rval = mbi->tag_get_handle(t_name.c_str(),
-                               tag_size,
-                               moab::MB_TYPE_DOUBLE,
-                               tally_tags[i],
-                               moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
-
-    if (rval != moab::MB_SUCCESS)
-      return rval;
-
-    rval = mbi->tag_get_handle(e_name.c_str(),
-                               tag_size,
-                               moab::MB_TYPE_DOUBLE,
-                               error_tags[i],
-                               moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
-
-    if (rval != moab::MB_SUCCESS)
-      return rval;
-  }
 
   return moab::MB_SUCCESS;
 }

src/tally/MeshTally.hpp

@@ -81,7 +81,9 @@ class MeshTally : public Tally {
   moab::Range tally_points;
 
   /// Tag arrays for storing energy bin labels
-  std::vector<moab::Tag> tally_tags, error_tags;
+  //  std::vector<moab::Tag> tally_tags, error_tags;
+  moab::Tag tally_tag, error_tag;
+  moab::Tag total_tally_tag, total_error_tag;
 
   // >>> PROTECTED METHODS
 

src/tally/Tally.cpp

@@ -20,7 +20,7 @@ Tally::Tally(const TallyInput& input)
   // TallyInput
   bool total_energy_bin = true;
 
-  unsigned int num_energy_bins = input_data.energy_bin_bounds.size() - 1;
+  unsigned int num_energy_bins = input_data.energy_bin_bounds.size();
 
   data = new TallyData(num_energy_bins, total_energy_bin);
 }
@@ -89,7 +89,7 @@ std::string Tally::get_tally_type() {
 //---------------------------------------------------------------------------//
 bool Tally::get_energy_bin(double energy, unsigned int& ebin) {
   bool bin_exists = false;
-
+ 
   if (energy_in_bounds(energy)) {
     // in bounds, energy bin index must exist
     bin_exists = true;

src/tally/TallyData.cpp

@@ -37,7 +37,8 @@ TallyData::TallyData(unsigned int num_energy_bins, bool total_energy_bin) {
 //---------------------------------------------------------------------------//
 std::pair <double, double> TallyData::get_data(unsigned int tally_point_index,
                                                unsigned int energy_bin) const {
-  assert(energy_bin < num_energy_bins);
+
+  assert(energy_bin < num_energy_bins + this->total_energy_bin);
   assert(tally_point_index < num_tally_points);
 
   int index = tally_point_index * num_energy_bins + energy_bin;

src/tally/TrackLengthMeshTally.cpp

@@ -236,6 +236,7 @@ void TrackLengthMeshTally::write_data(double num_histories) {
   }
   assert(rval == MB_SUCCESS);
 
+
   for (Range::const_iterator i = all_tets.begin(); i != all_tets.end(); ++i) {
     EntityHandle t = *i;
 
@@ -254,16 +255,47 @@ void TrackLengthMeshTally::write_data(double num_histories) {
     double volume = tet_volume(v[0], v[1], v[2], v[3]);
     unsigned int tet_index = get_entity_index(t);
 
-    for (unsigned j = 0; j < data->get_num_energy_bins(); ++j) {
+    unsigned int num_ebins = data->get_num_energy_bins();
+
+    // if there is a total, dont do anything with it
+    if (data->has_total_energy_bin())
+      num_ebins--;
+
+    double tally_vect[num_ebins];
+    double error_vect[num_ebins];
+
+    for (unsigned j = 0; j < num_ebins ; ++j) {
       std::pair <double, double> tally_data = data->get_data(tet_index, j);
       double tally = tally_data.first;
       double error = tally_data.second;
-      // double tally = get_data( tally_data, t, j );
-      // double error = get_data( error_data, t, j );
+
       double score = (tally / (volume * num_histories));
 
-      rval = mb->tag_set_data(tally_tags[j], &t, 1, &score);
-      assert(rval == MB_SUCCESS);
+      // Use 0 as the error output value if nothing has been computed for this mesh cell;
+      // this reflects MCNP's approach to avoiding a divide-by-zero situation.
+      double rel_err = 0;
+      if (error != 0) {
+        rel_err = sqrt((error / (tally * tally)) - (1. / num_histories));
+      }
+
+      tally_vect[j] = score;
+      error_vect[j] = rel_err;
+    }
+
+    rval = mb->tag_set_data(tally_tag, &t, 1, tally_vect);
+    MB_CHK_SET_ERR_RET(rval, "Failed to set tally_tag " + std::to_string(rval) + " " + std::to_string(t));
+    rval = mb->tag_set_data(error_tag, &t, 1, error_vect);
+    MB_CHK_SET_ERR_RET(rval, "Failed to set error_tag " + std::to_string(rval) + " " + std::to_string(t));
+
+
+    // if we have a total bin, write it out
+    if (data->has_total_energy_bin()) {
+      int j = num_ebins++;
+      std::pair <double, double> tally_data = data->get_data(tet_index, j);
+      double tally = tally_data.first;
+      double error = tally_data.second;
+
+      double score = (tally / (volume * num_histories));
 
       // Use 0 as the error output value if nothing has been computed for this mesh cell;
       // this reflects MCNP's approach to avoiding a divide-by-zero situation.
@@ -272,13 +304,20 @@ void TrackLengthMeshTally::write_data(double num_histories) {
         rel_err = sqrt((error / (tally * tally)) - (1. / num_histories));
       }
 
-      rval = mb->tag_set_data(error_tags[j], &t, 1, &rel_err);
-      assert(rval == MB_SUCCESS);
+      rval = mb->tag_set_data(total_tally_tag, &t, 1, &score);
+      MB_CHK_SET_ERR_RET(rval, "Failed to set tally_tag " + std::to_string(rval) + " " + std::to_string(t));
+      rval = mb->tag_set_data(total_error_tag, &t, 1, &rel_err);
+      MB_CHK_SET_ERR_RET(rval, "Failed to set error_tag " + std::to_string(rval) + " " + std::to_string(t));
     }
   }
 
-  std::vector<Tag> output_tags = tally_tags;
-  output_tags.insert(output_tags.end(), error_tags.begin(), error_tags.end());
+  std::vector<Tag> output_tags;
+  output_tags.push_back(tally_tag);
+  output_tags.push_back(error_tag);
+  if (data->has_total_energy_bin()) {
+    output_tags.push_back(total_tally_tag);
+    output_tags.push_back(total_error_tag);
+  }
 
   rval = mb->write_file(output_filename.c_str(), NULL, NULL, &tally_mesh_set, 1, &(output_tags[0]), output_tags.size());
   assert(rval == MB_SUCCESS);

src/tally/tests/test_TrackLengthMeshTally.cpp

@@ -1,6 +1,7 @@
 // test_TrackLengthMeshTally.cpp
 #include "gtest/gtest.h"
 
+#include "moab/Core.hpp"
 #include "moab/CartVect.hpp"
 
 #include "../Tally.hpp"
@@ -1213,3 +1214,59 @@ TEST_F(TrackLengthMeshTallyTest, ComputeScoreTallyManager1RayProton) {
 
   EXPECT_EQ(total, event.track_length);
 }
+
+//---------------------------------------------------------------------------//
+TEST_F(TrackLengthMeshTallyTest, EnsureVectorTags) {
+  input.tally_type = "unstr_track";
+  input.energy_bin_bounds.push_back(20.0);
+
+  input.options.insert(std::make_pair("inp", "unstructured_mesh.h5m"));
+  input.options.insert(std::make_pair("out", "mesh.out.h5m"));
+
+  // dummy variable to prevent segfault during teardown
+  mesh_tally = Tally::create_tally(input);
+
+  TallyEvent event;
+  make_event(event);
+  event.particle = 9; // proton
+
+  TallyManager* tallyManager = new TallyManager();
+  tallyManager->addNewTally(input.tally_id, input.tally_type, event.particle,
+                            input.energy_bin_bounds, input.options);
+
+  // write all 0's
+  tallyManager->writeData(1.);
+
+  moab::Core* MBI = new moab::Core();
+  // make a meshset
+  moab::EntityHandle mesh_set;
+  moab::ErrorCode rval = MBI->create_meshset(moab::MESHSET_SET, mesh_set);
+  EXPECT_EQ(rval, moab::MB_SUCCESS);
+  // load the mesh
+  rval = MBI->load_file("mesh.out.h5m");
+  EXPECT_EQ(rval, moab::MB_SUCCESS);
+  // look for the tag
+  moab::Tag tally_tag;
+  std::string tag_name = "TALLY_TAG";
+  // get the tag handle for the vector tag
+  rval = MBI->tag_get_handle(tag_name.c_str(),
+                             2,
+                             moab::MB_TYPE_DOUBLE,
+                             tally_tag,
+                             moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
+  // expect that this succeeded
+  EXPECT_EQ(rval, moab::MB_SUCCESS);
+
+  // now look for the scalar tag
+  tag_name = "TALLY_TAG_TOTAL";
+  // get the tag handle for the vector tag
+  rval = MBI->tag_get_handle(tag_name.c_str(),
+                             1,
+                             moab::MB_TYPE_DOUBLE,
+                             tally_tag,
+                             moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
+  // expect that this succeeded
+  EXPECT_EQ(rval, moab::MB_SUCCESS);
+
+  // all done :)
+}