From c758c66ebb53dae8360c7f8c4d877ba135650699 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Thu, 17 Aug 2017 12:32:14 -0500 Subject: [PATCH 1/8] [bug fix] DAGMC Advanced tallies were not producing vector tagged mesh data as expected by PyNE \ \Made it so that\a) the first bin 0-lowest bin is ignored \b) the total bin is no longer tagged as this was throwing off the correct number of bins calculation \c) changed the tally tag to a vector tag, have confirmed this all works in PyNE --- src/mcnp/meshtal_funcs.cpp | 5 ++- src/tally/KDEMeshTally.cpp | 30 +++++++------ src/tally/MeshTally.cpp | 69 ++++++++++++++++-------------- src/tally/MeshTally.hpp | 4 +- src/tally/TrackLengthMeshTally.cpp | 42 +++++++++++++----- 5 files changed, 93 insertions(+), 57 deletions(-) diff --git a/src/mcnp/meshtal_funcs.cpp b/src/mcnp/meshtal_funcs.cpp index 32c6658279..f068527e8e 100644 --- a/src/mcnp/meshtal_funcs.cpp +++ b/src/mcnp/meshtal_funcs.cpp @@ -139,7 +139,10 @@ 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 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]); } diff --git a/src/tally/KDEMeshTally.cpp b/src/tally/KDEMeshTally.cpp index f837eb69f1..f62074223f 100644 --- a/src/tally/KDEMeshTally.cpp +++ b/src/tally/KDEMeshTally.cpp @@ -183,7 +183,14 @@ 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 tally_data = data->get_data(point_index, j); double tally = tally_data.first; double error = tally_data.second; @@ -198,15 +205,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,8 +231,9 @@ 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 output_tags = tally_tags; - output_tags.insert(output_tags.end(), error_tags.begin(), error_tags.end()); + std::vector 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(), @@ -234,7 +241,6 @@ void KDEMeshTally::write_data(double num_histories) { &tally_mesh_set, 1, &(output_tags[0]), output_tags.size()); - assert(moab::MB_SUCCESS == rval); } //---------------------------------------------------------------------------// diff --git a/src/tally/MeshTally.cpp b/src/tally/MeshTally.cpp index 1e11589dc9..2bbb139d7f 100644 --- a/src/tally/MeshTally.cpp +++ b/src/tally/MeshTally.cpp @@ -113,40 +113,43 @@ 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; } diff --git a/src/tally/MeshTally.hpp b/src/tally/MeshTally.hpp index e6a3fb22e8..11ea1e703e 100644 --- a/src/tally/MeshTally.hpp +++ b/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 tally_tags, error_tags; + // std::vector tally_tags, error_tags; + moab::Tag tally_tag, error_tag; + moab::Tag total_tally_tag, total_error_tag; // >>> PROTECTED METHODS diff --git a/src/tally/TrackLengthMeshTally.cpp b/src/tally/TrackLengthMeshTally.cpp index ed2e2146a2..d517b9c0e7 100644 --- a/src/tally/TrackLengthMeshTally.cpp +++ b/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,19 @@ 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 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); + 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 +276,31 @@ 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); + tally_vect[j] = score; + error_vect[j] = rel_err; + + // TODO: It would be convenient for users to be able to get the total + // values tagged onto the dataset automatically, since our typical + // visualiation path is h5m -> vtk. The vtk doesnt understand + // vector tags so it would be useful to be able to the total_tag + // for easy visualisation -> maybe just fix pyne to be able to expand the + // tags of an arbitrary mesh + } + + 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)); } - std::vector output_tags = tally_tags; - output_tags.insert(output_tags.end(), error_tags.begin(), error_tags.end()); + /* + */ + + std::vector output_tags; + output_tags.push_back(tally_tag); + output_tags.push_back(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); From 1c12a7783956bec1a294d0d740cfcccc38b4e68a Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Fri, 15 Sep 2017 13:36:58 -0500 Subject: [PATCH 2/8] Now support scalar total tag, and vector flux tag --- src/tally/TallyData.cpp | 3 +- src/tally/TrackLengthMeshTally.cpp | 46 ++++++++++++++++++++---------- 2 files changed, 33 insertions(+), 16 deletions(-) diff --git a/src/tally/TallyData.cpp b/src/tally/TallyData.cpp index 807e641dd6..6c64abb56a 100644 --- a/src/tally/TallyData.cpp +++ b/src/tally/TallyData.cpp @@ -37,7 +37,8 @@ TallyData::TallyData(unsigned int num_energy_bins, bool total_energy_bin) { //---------------------------------------------------------------------------// std::pair 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; diff --git a/src/tally/TrackLengthMeshTally.cpp b/src/tally/TrackLengthMeshTally.cpp index d517b9c0e7..54059af9e6 100644 --- a/src/tally/TrackLengthMeshTally.cpp +++ b/src/tally/TrackLengthMeshTally.cpp @@ -256,9 +256,10 @@ void TrackLengthMeshTally::write_data(double num_histories) { unsigned int tet_index = get_entity_index(t); 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]; @@ -278,29 +279,44 @@ void TrackLengthMeshTally::write_data(double num_histories) { tally_vect[j] = score; error_vect[j] = rel_err; - - // TODO: It would be convenient for users to be able to get the total - // values tagged onto the dataset automatically, since our typical - // visualiation path is h5m -> vtk. The vtk doesnt understand - // vector tags so it would be useful to be able to the total_tag - // for easy visualisation -> maybe just fix pyne to be able to expand the - // tags of an arbitrary mesh - } - + 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 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. + double rel_err = 0; + if (error != 0) { + rel_err = sqrt((error / (tally * tally)) - (1. / num_histories)); + } + + 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 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); From 8c2d955eb6771e7a7748714abab7885d2b0582c6 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 20 Sep 2017 11:22:00 -0500 Subject: [PATCH 3/8] Unit test modified to add a vector tag test, assert that when egroups > 1 we get a vector tag and a total scalar tag --- src/tally/tests/test_TrackLengthMeshTally.cpp | 57 +++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/src/tally/tests/test_TrackLengthMeshTally.cpp b/src/tally/tests/test_TrackLengthMeshTally.cpp index ef9ba62ed1..2458e4a228 100644 --- a/src/tally/tests/test_TrackLengthMeshTally.cpp +++ b/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 :) +} From 0e79890fc8d23cd9980d5107631df2677059938a Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 20 Sep 2017 12:40:45 -0500 Subject: [PATCH 4/8] Applied SG --- src/mcnp/meshtal_funcs.cpp | 3 +- src/tally/KDEMeshTally.cpp | 5 +- src/tally/MeshTally.cpp | 47 ++++++++++--------- src/tally/TallyData.cpp | 2 +- src/tally/TrackLengthMeshTally.cpp | 39 +++++++-------- src/tally/tests/test_TrackLengthMeshTally.cpp | 32 ++++++------- 6 files changed, 66 insertions(+), 62 deletions(-) diff --git a/src/mcnp/meshtal_funcs.cpp b/src/mcnp/meshtal_funcs.cpp index f068527e8e..afc48b5bce 100644 --- a/src/mcnp/meshtal_funcs.cpp +++ b/src/mcnp/meshtal_funcs.cpp @@ -141,7 +141,8 @@ void dagmc_fmesh_setup_mesh_(int* fm_ipt, int* id, int* fmesh_idx, std::vector energy_boundaries; // if there is more than 1 bin, dont want 0->bottom bin int bottom = 0; - if(*n_energy_mesh > 2) bottom = 1; + if (*n_energy_mesh > 2) + bottom = 1; for (int i = bottom ; i < *n_energy_mesh; ++i) { energy_boundaries.push_back(energy_mesh[i]); } diff --git a/src/tally/KDEMeshTally.cpp b/src/tally/KDEMeshTally.cpp index f62074223f..31c332bedf 100644 --- a/src/tally/KDEMeshTally.cpp +++ b/src/tally/KDEMeshTally.cpp @@ -185,11 +185,12 @@ void KDEMeshTally::write_data(double num_histories) { 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--; + 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 tally_data = data->get_data(point_index, j); double tally = tally_data.first; diff --git a/src/tally/MeshTally.cpp b/src/tally/MeshTally.cpp index 2bbb139d7f..02de7fa3a9 100644 --- a/src/tally/MeshTally.cpp +++ b/src/tally/MeshTally.cpp @@ -113,42 +113,43 @@ moab::ErrorCode MeshTally::setup_tags(moab::Interface* mbi, const char* prefix) int tag_size = 1; unsigned int num_bins = data->get_num_energy_bins(); - if(data->has_total_energy_bin()) 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"); - + 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"); - + 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"); + 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"); + 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"); return moab::MB_SUCCESS; diff --git a/src/tally/TallyData.cpp b/src/tally/TallyData.cpp index 6c64abb56a..21ae1c1dfe 100644 --- a/src/tally/TallyData.cpp +++ b/src/tally/TallyData.cpp @@ -37,7 +37,7 @@ TallyData::TallyData(unsigned int num_energy_bins, bool total_energy_bin) { //---------------------------------------------------------------------------// std::pair TallyData::get_data(unsigned int tally_point_index, unsigned int energy_bin) const { - + assert(energy_bin < num_energy_bins + this->total_energy_bin); assert(tally_point_index < num_tally_points); diff --git a/src/tally/TrackLengthMeshTally.cpp b/src/tally/TrackLengthMeshTally.cpp index 54059af9e6..3e43a65aab 100644 --- a/src/tally/TrackLengthMeshTally.cpp +++ b/src/tally/TrackLengthMeshTally.cpp @@ -236,7 +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; @@ -258,11 +258,12 @@ void TrackLengthMeshTally::write_data(double num_histories) { 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--; - + 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 tally_data = data->get_data(tet_index, j); double tally = tally_data.first; @@ -280,12 +281,12 @@ void TrackLengthMeshTally::write_data(double 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)); - + + 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()) { @@ -293,23 +294,23 @@ void TrackLengthMeshTally::write_data(double num_histories) { std::pair 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. double rel_err = 0; if (error != 0) { - rel_err = sqrt((error / (tally * tally)) - (1. / num_histories)); + rel_err = sqrt((error / (tally * tally)) - (1. / num_histories)); } - - 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)); + + 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 output_tags; output_tags.push_back(tally_tag); output_tags.push_back(error_tag); diff --git a/src/tally/tests/test_TrackLengthMeshTally.cpp b/src/tally/tests/test_TrackLengthMeshTally.cpp index 2458e4a228..b023d3a9a5 100644 --- a/src/tally/tests/test_TrackLengthMeshTally.cpp +++ b/src/tally/tests/test_TrackLengthMeshTally.cpp @@ -1216,7 +1216,7 @@ TEST_F(TrackLengthMeshTallyTest, ComputeScoreTallyManager1RayProton) { } //---------------------------------------------------------------------------// -TEST_F(TrackLengthMeshTallyTest, EnsureVectorTags ) { +TEST_F(TrackLengthMeshTallyTest, EnsureVectorTags) { input.tally_type = "unstr_track"; input.energy_bin_bounds.push_back(20.0); @@ -1237,36 +1237,36 @@ TEST_F(TrackLengthMeshTallyTest, EnsureVectorTags ) { // write all 0's tallyManager->writeData(1.); - moab::Core *MBI = new moab::Core(); + 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); + 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); + 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); - + 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 :) } From 02a61eb8381c67da985fe56abcd1db9369c53b4d Mon Sep 17 00:00:00 2001 From: Lucas John Jacobson Date: Wed, 7 Nov 2018 09:14:36 -0600 Subject: [PATCH 5/8] Enable c++11 for dagtally --- CMakeLists.txt | 2 +- src/tally/CMakeLists.txt | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0b1226c1d6..c80dc63f91 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/src/tally/CMakeLists.txt b/src/tally/CMakeLists.txt index cf2fd9b264..185b94d97a 100644 --- a/src/tally/CMakeLists.txt +++ b/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) From ee89b9a31de1f569dfd0486a95b1852d52d990e4 Mon Sep 17 00:00:00 2001 From: Chelsea D'Angelo Date: Fri, 30 Nov 2018 22:41:39 +0000 Subject: [PATCH 6/8] lots of print msgs to find e group mismatch --- src/mcnp/meshtal_funcs.cpp | 3 +++ src/tally/MeshTally.cpp | 2 ++ src/tally/Tally.cpp | 2 ++ src/tally/TallyData.cpp | 5 +++++ 4 files changed, 12 insertions(+) diff --git a/src/mcnp/meshtal_funcs.cpp b/src/mcnp/meshtal_funcs.cpp index afc48b5bce..9e95ee4f49 100644 --- a/src/mcnp/meshtal_funcs.cpp +++ b/src/mcnp/meshtal_funcs.cpp @@ -124,6 +124,7 @@ void dagmc_fmesh_setup_mesh_(int* fm_ipt, int* id, int* fmesh_idx, std::cout << "DAGMC tally " << *id << " has these " << *n_energy_mesh << " energy bin boundaries: " << std::endl; + std::cout << "meshtal funcs: n_energy_mesh " << *n_energy_mesh << std::endl; for (int i = 0; i < *n_energy_mesh; ++i) { std::cout << " " << energy_mesh[i] << std::endl; } @@ -144,8 +145,10 @@ void dagmc_fmesh_setup_mesh_(int* fm_ipt, int* id, int* fmesh_idx, if (*n_energy_mesh > 2) bottom = 1; for (int i = bottom ; i < *n_energy_mesh; ++i) { + std::cout << "energy mesh[i] " << energy_mesh[i] << std::endl; energy_boundaries.push_back(energy_mesh[i]); } + std::cout << "meshtal funcs: size energy bounds " << energy_boundaries.size() << std::endl; // Parse FC card and create input data for MeshTally std::multimap fc_settings; diff --git a/src/tally/MeshTally.cpp b/src/tally/MeshTally.cpp index 02de7fa3a9..c0cf8523e6 100644 --- a/src/tally/MeshTally.cpp +++ b/src/tally/MeshTally.cpp @@ -116,6 +116,8 @@ moab::ErrorCode MeshTally::setup_tags(moab::Interface* mbi, const char* prefix) if (data->has_total_energy_bin()) num_bins--; + std::cout << "meshtally: num_bins " << std::endl; + std::string tag_name = pfx + "TALLY_TAG"; std::string error_tag_name = pfx + "ERROR_TAG"; diff --git a/src/tally/Tally.cpp b/src/tally/Tally.cpp index 89ddd59f3e..b8e6faadf1 100644 --- a/src/tally/Tally.cpp +++ b/src/tally/Tally.cpp @@ -19,8 +19,10 @@ Tally::Tally(const TallyInput& input) // This is a placeholder for a future option to set t.e.b. false via the // TallyInput bool total_energy_bin = true; + std::cout << "Tally: tot bin " << total_energy_bin << std::endl; unsigned int num_energy_bins = input_data.energy_bin_bounds.size() - 1; + std::cout << "Tally: num_energy_bins " << num_energy_bins << std::endl; data = new TallyData(num_energy_bins, total_energy_bin); } diff --git a/src/tally/TallyData.cpp b/src/tally/TallyData.cpp index 21ae1c1dfe..0f73016889 100644 --- a/src/tally/TallyData.cpp +++ b/src/tally/TallyData.cpp @@ -21,15 +21,20 @@ TallyData::TallyData(unsigned int num_energy_bins, bool total_energy_bin) { // Otherwise, pass through the value this->total_energy_bin = total_energy_bin; } + std::cout << "TD: total bin? " << total_energy_bin << std::endl; + std::cout << "TD: num_energy_bins " << num_energy_bins << std::endl; ////////////////////////////////////////////////// if (this->total_energy_bin) { // Add an extra bin for the total, if required this->num_energy_bins = num_energy_bins + 1; + std::cout << "tally data: total bin? " << total_energy_bin << std::endl; + std::cout << "tally data: num_energy_bins " << num_energy_bins << std::endl; } else { this->num_energy_bins = num_energy_bins; } + std::cout << "tally data: num_energy_bins " << num_energy_bins << std::endl; this->num_tally_points = 0; } //---------------------------------------------------------------------------// From a4a57880ed4ea96cecef502ec81dfb212b253b55 Mon Sep 17 00:00:00 2001 From: chelsea Date: Tue, 4 Dec 2018 15:18:17 -0600 Subject: [PATCH 7/8] print outs to show num ebins and bounds --- src/tally/MeshTally.cpp | 3 ++- src/tally/Tally.cpp | 6 ++++-- src/tally/TallyData.cpp | 6 +++--- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/tally/MeshTally.cpp b/src/tally/MeshTally.cpp index c0cf8523e6..fa46be34fc 100644 --- a/src/tally/MeshTally.cpp +++ b/src/tally/MeshTally.cpp @@ -113,10 +113,11 @@ moab::ErrorCode MeshTally::setup_tags(moab::Interface* mbi, const char* prefix) int tag_size = 1; unsigned int num_bins = data->get_num_energy_bins(); + std::cout << "MT: num bins " << num_bins << std::endl; if (data->has_total_energy_bin()) num_bins--; - std::cout << "meshtally: num_bins " << std::endl; + std::cout << "meshtally: num_bins " << num_bins << std::endl; std::string tag_name = pfx + "TALLY_TAG"; std::string error_tag_name = pfx + "ERROR_TAG"; diff --git a/src/tally/Tally.cpp b/src/tally/Tally.cpp index b8e6faadf1..9656846604 100644 --- a/src/tally/Tally.cpp +++ b/src/tally/Tally.cpp @@ -21,7 +21,8 @@ Tally::Tally(const TallyInput& input) bool total_energy_bin = true; std::cout << "Tally: tot bin " << total_energy_bin << std::endl; - unsigned int num_energy_bins = input_data.energy_bin_bounds.size() - 1; +// unsigned int num_energy_bins = input_data.energy_bin_bounds.size() - 1; + unsigned int num_energy_bins = input_data.energy_bin_bounds.size(); std::cout << "Tally: num_energy_bins " << num_energy_bins << std::endl; data = new TallyData(num_energy_bins, total_energy_bin); @@ -91,7 +92,8 @@ std::string Tally::get_tally_type() { //---------------------------------------------------------------------------// bool Tally::get_energy_bin(double energy, unsigned int& ebin) { bool bin_exists = false; - + + //std::cout << "TD: input data ebin bound 174 " << input_data.energy_bin_bounds.at(174) << std::endl; if (energy_in_bounds(energy)) { // in bounds, energy bin index must exist bin_exists = true; diff --git a/src/tally/TallyData.cpp b/src/tally/TallyData.cpp index 0f73016889..2cbe09fcaa 100644 --- a/src/tally/TallyData.cpp +++ b/src/tally/TallyData.cpp @@ -21,8 +21,8 @@ TallyData::TallyData(unsigned int num_energy_bins, bool total_energy_bin) { // Otherwise, pass through the value this->total_energy_bin = total_energy_bin; } - std::cout << "TD: total bin? " << total_energy_bin << std::endl; - std::cout << "TD: num_energy_bins " << num_energy_bins << std::endl; + std::cout << "TD: total bin? " << this->total_energy_bin << std::endl; + std::cout << "TD: num_energy_bins " << this->num_energy_bins << std::endl; ////////////////////////////////////////////////// if (this->total_energy_bin) { @@ -34,7 +34,7 @@ TallyData::TallyData(unsigned int num_energy_bins, bool total_energy_bin) { this->num_energy_bins = num_energy_bins; } - std::cout << "tally data: num_energy_bins " << num_energy_bins << std::endl; + std::cout << "tally data: num_energy_bins " << this->num_energy_bins << std::endl; this->num_tally_points = 0; } //---------------------------------------------------------------------------// From 697a1d55be81153de18f8608319655e9d96ae49b Mon Sep 17 00:00:00 2001 From: chelsea Date: Tue, 4 Dec 2018 15:23:24 -0600 Subject: [PATCH 8/8] removed subtraction from energy bin bounds in Tally.cpp --- src/mcnp/meshtal_funcs.cpp | 3 --- src/tally/MeshTally.cpp | 3 --- src/tally/Tally.cpp | 4 ---- src/tally/TallyData.cpp | 5 ----- 4 files changed, 15 deletions(-) diff --git a/src/mcnp/meshtal_funcs.cpp b/src/mcnp/meshtal_funcs.cpp index 9e95ee4f49..afc48b5bce 100644 --- a/src/mcnp/meshtal_funcs.cpp +++ b/src/mcnp/meshtal_funcs.cpp @@ -124,7 +124,6 @@ void dagmc_fmesh_setup_mesh_(int* fm_ipt, int* id, int* fmesh_idx, std::cout << "DAGMC tally " << *id << " has these " << *n_energy_mesh << " energy bin boundaries: " << std::endl; - std::cout << "meshtal funcs: n_energy_mesh " << *n_energy_mesh << std::endl; for (int i = 0; i < *n_energy_mesh; ++i) { std::cout << " " << energy_mesh[i] << std::endl; } @@ -145,10 +144,8 @@ void dagmc_fmesh_setup_mesh_(int* fm_ipt, int* id, int* fmesh_idx, if (*n_energy_mesh > 2) bottom = 1; for (int i = bottom ; i < *n_energy_mesh; ++i) { - std::cout << "energy mesh[i] " << energy_mesh[i] << std::endl; energy_boundaries.push_back(energy_mesh[i]); } - std::cout << "meshtal funcs: size energy bounds " << energy_boundaries.size() << std::endl; // Parse FC card and create input data for MeshTally std::multimap fc_settings; diff --git a/src/tally/MeshTally.cpp b/src/tally/MeshTally.cpp index fa46be34fc..02de7fa3a9 100644 --- a/src/tally/MeshTally.cpp +++ b/src/tally/MeshTally.cpp @@ -113,12 +113,9 @@ moab::ErrorCode MeshTally::setup_tags(moab::Interface* mbi, const char* prefix) int tag_size = 1; unsigned int num_bins = data->get_num_energy_bins(); - std::cout << "MT: num bins " << num_bins << std::endl; if (data->has_total_energy_bin()) num_bins--; - std::cout << "meshtally: num_bins " << num_bins << std::endl; - std::string tag_name = pfx + "TALLY_TAG"; std::string error_tag_name = pfx + "ERROR_TAG"; diff --git a/src/tally/Tally.cpp b/src/tally/Tally.cpp index 9656846604..b2d10a419f 100644 --- a/src/tally/Tally.cpp +++ b/src/tally/Tally.cpp @@ -19,11 +19,8 @@ Tally::Tally(const TallyInput& input) // This is a placeholder for a future option to set t.e.b. false via the // TallyInput bool total_energy_bin = true; - std::cout << "Tally: tot bin " << total_energy_bin << std::endl; -// unsigned int num_energy_bins = input_data.energy_bin_bounds.size() - 1; unsigned int num_energy_bins = input_data.energy_bin_bounds.size(); - std::cout << "Tally: num_energy_bins " << num_energy_bins << std::endl; data = new TallyData(num_energy_bins, total_energy_bin); } @@ -93,7 +90,6 @@ std::string Tally::get_tally_type() { bool Tally::get_energy_bin(double energy, unsigned int& ebin) { bool bin_exists = false; - //std::cout << "TD: input data ebin bound 174 " << input_data.energy_bin_bounds.at(174) << std::endl; if (energy_in_bounds(energy)) { // in bounds, energy bin index must exist bin_exists = true; diff --git a/src/tally/TallyData.cpp b/src/tally/TallyData.cpp index 2cbe09fcaa..21ae1c1dfe 100644 --- a/src/tally/TallyData.cpp +++ b/src/tally/TallyData.cpp @@ -21,20 +21,15 @@ TallyData::TallyData(unsigned int num_energy_bins, bool total_energy_bin) { // Otherwise, pass through the value this->total_energy_bin = total_energy_bin; } - std::cout << "TD: total bin? " << this->total_energy_bin << std::endl; - std::cout << "TD: num_energy_bins " << this->num_energy_bins << std::endl; ////////////////////////////////////////////////// if (this->total_energy_bin) { // Add an extra bin for the total, if required this->num_energy_bins = num_energy_bins + 1; - std::cout << "tally data: total bin? " << total_energy_bin << std::endl; - std::cout << "tally data: num_energy_bins " << num_energy_bins << std::endl; } else { this->num_energy_bins = num_energy_bins; } - std::cout << "tally data: num_energy_bins " << this->num_energy_bins << std::endl; this->num_tally_points = 0; } //---------------------------------------------------------------------------//