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/mcnp/meshtal_funcs.cpp b/src/mcnp/meshtal_funcs.cpp index 32c6658279..afc48b5bce 100644 --- a/src/mcnp/meshtal_funcs.cpp +++ b/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 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/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) diff --git a/src/tally/KDEMeshTally.cpp b/src/tally/KDEMeshTally.cpp index f837eb69f1..31c332bedf 100644 --- a/src/tally/KDEMeshTally.cpp +++ b/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 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,8 +232,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 +242,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..02de7fa3a9 100644 --- a/src/tally/MeshTally.cpp +++ b/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; } 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/Tally.cpp b/src/tally/Tally.cpp index 89ddd59f3e..b2d10a419f 100644 --- a/src/tally/Tally.cpp +++ b/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; diff --git a/src/tally/TallyData.cpp b/src/tally/TallyData.cpp index 807e641dd6..21ae1c1dfe 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 ed2e2146a2..3e43a65aab 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,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 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 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 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); + 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); diff --git a/src/tally/tests/test_TrackLengthMeshTally.cpp b/src/tally/tests/test_TrackLengthMeshTally.cpp index ef9ba62ed1..b023d3a9a5 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 :) +}