From d0b6d9fbd640e086d68af15fd0b4444f0b3e08b3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Feb 2019 23:34:36 -0600 Subject: [PATCH] Move read_tallies_xml to C++ --- CMakeLists.txt | 1 - include/openmc/mesh.h | 2 +- include/openmc/settings.h | 2 +- include/openmc/tallies/derivative.h | 3 + include/openmc/tallies/filter.h | 14 +- include/openmc/tallies/tally.h | 8 +- src/initialize.cpp | 2 +- src/input_xml.F90 | 471 ---------------------------- src/mesh.cpp | 4 +- src/settings.cpp | 2 +- src/tallies/derivative.cpp | 6 +- src/tallies/filter.cpp | 116 +++---- src/tallies/tally.cpp | 334 +++++++++++++++++++- src/tallies/tally_filter.F90 | 10 +- 14 files changed, 413 insertions(+), 562 deletions(-) delete mode 100644 src/input_xml.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index b89b8a9825..342dc46361 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -310,7 +310,6 @@ add_library(libopenmc SHARED src/geometry.F90 src/hdf5_interface.F90 src/initialize.F90 - src/input_xml.F90 src/material_header.F90 src/message_passing.F90 src/particle_header.F90 diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 1ca4d57b4a..389590e0d0 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -124,7 +124,7 @@ private: //! Read meshes from either settings/tallies //! \param[in] root XML node -extern "C" void read_meshes(pugi::xml_node* root); +void read_meshes(pugi::xml_node root); //! Write mesh data to an HDF5 group //! \param[in] group HDF5 group diff --git a/include/openmc/settings.h b/include/openmc/settings.h index e8ff62972e..49d06f7f68 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -23,7 +23,7 @@ namespace openmc { namespace settings { // Boolean flags -extern "C" bool assume_separate; //!< assume tallies are spatially separate? +extern bool assume_separate; //!< assume tallies are spatially separate? extern bool check_overlaps; //!< check overlaps in geometry? extern bool confidence_intervals; //!< use confidence intervals for results? extern bool create_fission_neutrons; //!< create fission neutrons (fixed source)? diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index 02532afc08..673afb0c94 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -29,6 +29,9 @@ struct TallyDerivative { // Non-method functions //============================================================================== +//! Read tally derivatives from a tallies.xml file +void read_tally_derivatives(pugi::xml_node node); + //! Scale the given score by its logarithmic derivative void diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index b601685943..1ea3ad8fc8 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -4,6 +4,7 @@ #include #include #include +#include #include #include "openmc/hdf5_interface.h" @@ -89,13 +90,16 @@ extern std::vector filter_matches; } // namespace simulation namespace model { - -extern "C" int32_t n_filters; -extern std::vector> tally_filters; - -} // namespace model + extern "C" int32_t n_filters; + extern std::vector> tally_filters; + extern std::unordered_map filter_map; +} //============================================================================== +// Non-member functions +//============================================================================== + +Filter* allocate_filter(const std::string& type); // Filter-related Fortran functions that will be called from C++ extern "C" int verify_filter(int32_t index); diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 75fa3d10a1..7719668191 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -9,6 +9,7 @@ #include "xtensor/xtensor.hpp" #include // for unique_ptr +#include #include #include @@ -46,7 +47,7 @@ public: //---------------------------------------------------------------------------- // Other methods. - void init_triggers(pugi::xml_node node, int i_tally); + void init_triggers(pugi::xml_node node); void init_results(); @@ -123,6 +124,8 @@ namespace model { extern std::vector active_collision_tallies; extern std::vector active_meshsurf_tallies; extern std::vector active_surface_tallies; + + extern std::unordered_map tally_map; } namespace simulation { @@ -149,6 +152,9 @@ extern double global_tally_leakage; // Non-member functions //============================================================================== +//! Read tally specification from tallies.xml +void read_tallies_xml(); + //! \brief Accumulate the sum of the contributions from each history within the //! batch to a new random variable void accumulate_tallies(); diff --git a/src/initialize.cpp b/src/initialize.cpp index b1caccecbd..553ad5d3db 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -28,12 +28,12 @@ #include "openmc/simulation.h" #include "openmc/string_utils.h" #include "openmc/summary.h" +#include "openmc/tallies/tally.h" #include "openmc/thermal.h" #include "openmc/timer.h" // data/functions from Fortran side extern "C" void read_command_line(); -extern "C" void read_tallies_xml(); // Paths to various files extern "C" { diff --git a/src/input_xml.F90 b/src/input_xml.F90 deleted file mode 100644 index b28380ba35..0000000000 --- a/src/input_xml.F90 +++ /dev/null @@ -1,471 +0,0 @@ -module input_xml - - use, intrinsic :: ISO_C_BINDING - - use constants - use error, only: fatal_error, warning, write_message, openmc_err_msg - use material_header - use message_passing - use settings - use stl_vector, only: VectorInt, VectorReal, VectorChar - use string, only: to_lower, to_str, str_to_int, & - starts_with, ends_with, to_c_string - use tally - use tally_header, only: openmc_extend_tallies - use tally_derivative_header - use tally_filter_header - use tally_filter - use xml_interface - - implicit none - save - -contains - -!=============================================================================== -! READ_TALLIES_XML reads data from a tallies.xml file and parses it, checking -! for errors and placing properly-formatted data in the right data structures -!=============================================================================== - - subroutine read_tallies_xml() bind(C) - - integer :: i ! loop over user-specified tallies - integer :: j ! loop over words - integer :: l ! loop over bins - integer :: filter_id ! user-specified identifier for filter - integer :: tally_id ! user-specified identifier for filter - integer :: deriv_id - integer :: i_filt ! index in filters array - integer :: n ! size of arrays in mesh specification - integer :: n_words ! number of words read - integer :: n_filter ! number of filters - integer :: i_start, i_end - integer(C_INT) :: err - logical :: file_exists ! does tallies.xml file exist? - integer, allocatable :: temp_filter(:) ! temporary filter indices - logical :: has_energyout - integer :: particle_filter_index - character(MAX_LINE_LEN) :: filename - character(MAX_WORD_LEN) :: temp_str - type(TallyFilterContainer), pointer :: f - type(XMLDocument) :: doc - type(XMLNode) :: root - type(XMLNode) :: node_tal - type(XMLNode) :: node_filt - type(XMLNode), allocatable :: node_tal_list(:) - type(XMLNode), allocatable :: node_filt_list(:) - type(TallyDerivative), pointer :: deriv - - interface - subroutine tally_init_from_xml(tally_ptr, xml_node) bind(C) - import C_PTR - type(C_PTR), value :: tally_ptr - type(C_PTR) :: xml_node - end subroutine - - subroutine tally_set_scores(tally_ptr, xml_node) bind(C) - import C_PTR - type(C_PTR), value :: tally_ptr - type(C_PTR) :: xml_node - end subroutine - - subroutine tally_set_nuclides(tally_ptr, xml_node) bind(C) - import C_PTR - type(C_PTR), value :: tally_ptr - type(C_PTR) :: xml_node - end subroutine - - subroutine tally_init_triggers(tally_ptr, i_tally, xml_node) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: tally_ptr - integer(C_INT), value :: i_tally - type(C_PTR) :: xml_node - end subroutine - - subroutine read_tally_derivatives(node_ptr) bind(C) - import C_PTR - type(C_PTR) :: node_ptr - end subroutine - - subroutine read_meshes(node_ptr) bind(C) - import C_PTR - type(C_PTR) :: node_ptr - end subroutine - end interface - - ! Check if tallies.xml exists - filename = trim(path_input) // "tallies.xml" - inquire(FILE=filename, EXIST=file_exists) - if (.not. file_exists) then - ! Since a tallies.xml file is optional, no error is issued here - return - end if - - ! Display output message - call write_message("Reading tallies XML file...", 5) - - ! Parse tallies.xml file - call doc % load_file(filename) - root = doc % document_element() - - ! ========================================================================== - ! DETERMINE SIZE OF ARRAYS AND ALLOCATE - - ! Get pointer list to XML - call get_node_list(root, "filter", node_filt_list) - - ! Get pointer list to XML - call get_node_list(root, "tally", node_tal_list) - - ! Check for setting - if (check_for_node(root, "assume_separate")) then - call get_node_value(root, "assume_separate", assume_separate) - end if - - ! ========================================================================== - ! READ MESH DATA - - ! Check for user meshes and allocate - call read_meshes(root % ptr) - - ! We only need the mesh info for plotting - if (run_mode == MODE_PLOTTING) then - call doc % clear() - return - end if - - ! ========================================================================== - ! READ DATA FOR DERIVATIVES - - call read_tally_derivatives(root % ptr) - - ! ========================================================================== - ! READ FILTER DATA - - ! Check for user filters and allocate - n = size(node_filt_list) - if (n > 0) then - err = openmc_extend_filters(n, i_start, i_end) - end if - - READ_FILTERS: do i = 1, n - f => filters(i_start + i - 1) - - ! Get pointer to filter xml node - node_filt = node_filt_list(i) - - ! Copy filter id - if (check_for_node(node_filt, "id")) then - call get_node_value(node_filt, "id", filter_id) - else - call fatal_error("Must specify id for filter in tally XML file.") - end if - - ! Check to make sure 'id' hasn't been used - if (filter_dict % has(filter_id)) then - call fatal_error("Two or more filters use the same unique ID: " & - // to_str(filter_id)) - end if - - ! Convert filter type to lower case - temp_str = '' - if (check_for_node(node_filt, "type")) & - call get_node_value(node_filt, "type", temp_str) - temp_str = to_lower(temp_str) - - ! Make sure bins have been set - select case(temp_str) - case ("energy", "energyout", "mu", "polar", "azimuthal") - if (.not. check_for_node(node_filt, "bins")) then - call fatal_error("Bins not set in filter " // trim(to_str(filter_id))) - end if - case ("mesh", "meshsurface", "universe", "material", "cell", "distribcell", & - "cellborn", "cellfrom", "surface", "delayedgroup") - if (.not. check_for_node(node_filt, "bins")) then - call fatal_error("Bins not set in filter " // trim(to_str(filter_id))) - end if - end select - - ! Allocate according to the filter type - err = openmc_filter_set_type(i_start + i - 1, to_c_string(temp_str)) - if (err /= 0) call fatal_error(to_f_string(openmc_err_msg)) - - ! Read filter data from XML - call f % obj % from_xml(node_filt) - - ! Set filter id - err = openmc_filter_set_id(i_start + i - 1, filter_id) - - ! Initialize filter - call f % obj % initialize() - end do READ_FILTERS - - ! ========================================================================== - ! READ TALLY DATA - - ! Check for user tallies - n = size(node_tal_list) - if (n == 0) then - if (master) call warning("No tallies present in tallies.xml file!") - end if - - ! Allocate user tallies - if (n > 0 .and. run_mode /= MODE_PLOTTING) then - err = openmc_extend_tallies(n, i_start, i_end) - end if - - READ_TALLIES: do i = 1, n - ! Allocate tally - err = openmc_tally_allocate(i_start + i - 1, & - C_CHAR_'generic' // C_NULL_CHAR) - - ! Get pointer to tally - associate (t => tallies(i_start + i - 1) % obj) - - ! Get pointer to tally xml node - node_tal = node_tal_list(i) - - call tally_init_from_xml(t % ptr, node_tal % ptr) - - ! Copy and set tally id - if (check_for_node(node_tal, "id")) then - call get_node_value(node_tal, "id", tally_id) - err = openmc_tally_set_id(i_start + i - 1, tally_id) - if (err /= 0) call fatal_error(to_f_string(openmc_err_msg)) - else - call fatal_error("Must specify id for tally in tally XML file.") - end if - - ! Copy tally name - if (check_for_node(node_tal, "name")) & - call get_node_value(node_tal, "name", t % name) - - ! ======================================================================= - ! READ DATA FOR FILTERS - - ! Check if user is using old XML format and throw an error if so - if (check_for_node(node_tal, "filter")) then - call fatal_error("Tally filters must be specified independently of & - &tallies in a element. The element itself should & - &have a list of filters that apply, e.g., 1 2 & - &where 1 and 2 are the IDs of filters specified outside of & - &.") - end if - - ! Determine number of filters - if (check_for_node(node_tal, "filters")) then - n_filter = node_word_count(node_tal, "filters") - else - n_filter = 0 - end if - - ! Allocate and store filter user ids - allocate(temp_filter(n_filter)) - if (n_filter > 0) then - call get_node_array(node_tal, "filters", temp_filter) - - do j = 1, n_filter - ! Get pointer to filter - if (filter_dict % has(temp_filter(j))) then - i_filt = filter_dict % get(temp_filter(j)) - f => filters(i_filt) - else - call fatal_error("Could not find filter " & - // trim(to_str(temp_filter(j))) // " specified on tally " & - // trim(to_str(t % id()))) - end if - - ! Store the index of the filter - temp_filter(j) = i_filt - 1 - end do - - ! Set the filters - err = openmc_tally_set_filters(i_start + i - 1, n_filter, temp_filter) - end if - deallocate(temp_filter) - - ! Check for the presence of certain filter types - has_energyout = (t % energyout_filter() > 0) - particle_filter_index = 0 - do j = 1, t % n_filters() - select type (filt => filters(t % filter(j) + 1) % obj) - type is (ParticleFilter) - particle_filter_index = j - end select - end do - - ! Change the tally estimator if a filter demands it - do j = 1, t % n_filters() - select type (filt => filters(t % filter(j) + 1) % obj) - type is (EnergyoutFilter) - call t % set_estimator(ESTIMATOR_ANALOG) - type is (LegendreFilter) - call t % set_estimator(ESTIMATOR_ANALOG) - type is (SphericalHarmonicsFilter) - if (filt % cosine() == COSINE_SCATTER) then - call t % set_estimator(ESTIMATOR_ANALOG) - end if - type is (SpatialLegendreFilter) - call t % set_estimator(ESTIMATOR_COLLISION) - type is (ZernikeFilter) - call t % set_estimator(ESTIMATOR_COLLISION) - type is (ZernikeRadialFilter) - call t % set_estimator(ESTIMATOR_COLLISION) - end select - end do - - ! ======================================================================= - ! READ DATA FOR NUCLIDES - - call tally_set_nuclides(t % ptr, node_tal % ptr) - - ! ======================================================================= - ! READ DATA FOR SCORES - - call tally_set_scores(t % ptr, node_tal % ptr) - - if (check_for_node(node_tal, "scores")) then - n_words = node_word_count(node_tal, "scores") - - ! Check if tally is compatible with particle type - if (photon_transport) then - if (particle_filter_index == 0) then - do j = 1, t % n_score_bins() - select case (t % score_bins(j)) - case (SCORE_INVERSE_VELOCITY) - call fatal_error("Particle filter must be used with photon & - &transport on and inverse velocity score") - case (SCORE_FLUX, SCORE_TOTAL, SCORE_SCATTER, SCORE_NU_SCATTER, & - SCORE_ABSORPTION, SCORE_FISSION, SCORE_NU_FISSION, & - SCORE_CURRENT, SCORE_EVENTS, SCORE_DELAYED_NU_FISSION, & - SCORE_PROMPT_NU_FISSION, SCORE_DECAY_RATE) - call warning("Particle filter is not used with photon transport& - & on and " // trim(to_str(t % score_bins(j))) // " score") - end select - end do - else - select type(filt => filters(particle_filter_index) % obj) - type is (ParticleFilter) - do l = 1, filt % n_bins - if (filt % particles(l) == ELECTRON .or. filt % particles(l) == POSITRON) then - call t % set_estimator(ESTIMATOR_ANALOG) - end if - end do - end select - end if - else - if (particle_filter_index > 0) then - select type(filt => filters(particle_filter_index) % obj) - type is (ParticleFilter) - do l = 1, filt % n_bins - if (filt % particles(l) /= NEUTRON) then - call warning("Particle filter other than NEUTRON used with & - &photon transport turned off. All tallies for particle & - &type " // trim(to_str(filt % particles(l))) // " will have no scores") - end if - end do - end select - end if - end if - else - call fatal_error("No specified on tally " & - // trim(to_str(t % id())) // ".") - end if - - ! Check for a tally derivative. - if (check_for_node(node_tal, "derivative")) then - call get_node_value(node_tal, "derivative", deriv_id) - - ! Find the derivative with the given id, and store it's index. - do j = 0, n_tally_derivs() - 1 - deriv => tally_deriv_c(j) - if (deriv % id == deriv_id) then - call t % set_deriv(j) - ! Only analog or collision estimators are supported for differential - ! tallies. - if (t % estimator() == ESTIMATOR_TRACKLENGTH) then - call t % set_estimator(ESTIMATOR_COLLISION) - end if - ! We found the derivative we were looking for; exit the do loop. - exit - end if - if (j == n_tally_derivs()) then - call fatal_error("Could not find derivative " & - // trim(to_str(deriv_id)) // " specified on tally " & - // trim(to_str(t % id()))) - end if - end do - - deriv => tally_deriv_c(t % deriv()) - if (deriv % variable == DIFF_NUCLIDE_DENSITY & - .or. deriv % variable == DIFF_TEMPERATURE) then - do j = 1, t % n_nuclide_bins() - if (has_energyout .and. t % nuclide_bins(j) == -1) then - call fatal_error("Error on tally " // trim(to_str(t % id())) & - // ": Cannot use a 'nuclide_density' or 'temperature' & - &derivative on a tally with an outgoing energy filter and & - &'total' nuclide rate. Instead, tally each nuclide in the & - &material individually.") - ! Note that diff tallies with these characteristics would work - ! correctly if no tally events occur in the perturbed material - ! (e.g. pertrubing moderator but only tallying fuel), but this - ! case would be hard to check for by only reading inputs. - end if - end do - end if - end if - - ! If settings.xml trigger is turned on, create tally triggers - if (trigger_on) then - !TODO: off-by-one - call tally_init_triggers(t % ptr, i_start + i - 1 - 1, node_tal % ptr) - end if - - ! ======================================================================= - ! SET TALLY ESTIMATOR - - ! Check if user specified estimator - if (check_for_node(node_tal, "estimator")) then - temp_str = '' - call get_node_value(node_tal, "estimator", temp_str) - select case(trim(temp_str)) - case ('analog') - call t % set_estimator(ESTIMATOR_ANALOG) - - case ('tracklength', 'track-length', 'pathlength', 'path-length') - ! If the estimator was set to an analog estimator, this means the - ! tally needs post-collision information - if (t % estimator() == ESTIMATOR_ANALOG) then - call fatal_error("Cannot use track-length estimator for tally " & - // to_str(t % id())) - end if - - ! Set estimator to track-length estimator - call t % set_estimator(ESTIMATOR_TRACKLENGTH) - - case ('collision') - ! If the estimator was set to an analog estimator, this means the - ! tally needs post-collision information - if (t % estimator() == ESTIMATOR_ANALOG) then - call fatal_error("Cannot use collision estimator for tally " & - // to_str(t % id())) - end if - - ! Set estimator to collision estimator - call t % set_estimator(ESTIMATOR_COLLISION) - - case default - call fatal_error("Invalid estimator '" // trim(temp_str) & - // "' on tally " // to_str(t % id())) - end select - end if - - end associate - end do READ_TALLIES - - ! Close XML document - call doc % clear() - - end subroutine read_tallies_xml - -end module input_xml diff --git a/src/mesh.cpp b/src/mesh.cpp index 6ec298a78a..1681c23619 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -863,9 +863,9 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, // Non-member functions //============================================================================== -void read_meshes(pugi::xml_node* root) +void read_meshes(pugi::xml_node root) { - for (auto node : root->children("mesh")) { + for (auto node : root.children("mesh")) { // Read mesh and add to vector model::meshes.emplace_back(new RegularMesh{node}); diff --git a/src/settings.cpp b/src/settings.cpp index ec62c4aa62..0d4878f3d8 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -490,7 +490,7 @@ void read_settings_xml() } // Read meshes - read_meshes(&root); + read_meshes(root); // Shannon Entropy mesh if (check_for_node(root, "entropy_mesh")) { diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index b28b78d0ae..e510335347 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -79,14 +79,14 @@ TallyDerivative::TallyDerivative(pugi::xml_node node) // Non-method functions //============================================================================== -extern "C" void -read_tally_derivatives(pugi::xml_node* node) +void +read_tally_derivatives(pugi::xml_node node) { // Populate the derivatives array. This must be done in parallel because // the derivatives are threadprivate. #pragma omp parallel { - for (auto deriv_node : node->children("derivative")) + for (auto deriv_node : node.children("derivative")) model::tally_derivs.emplace_back(deriv_node); } diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 7670fddf42..66a26a5bc8 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -36,12 +36,67 @@ namespace openmc { //============================================================================== namespace simulation { -std::vector filter_matches; -} // namespace simulation + std::vector filter_matches; +} namespace model { -std::vector> tally_filters; -} // namespace model + std::vector> tally_filters; + std::unordered_map filter_map; +} + +Filter* +allocate_filter(const std::string& type) +{ + if (type == "azimuthal") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "cell") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "cellborn") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "cellfrom") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "distribcell") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "delayedgroup") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "energyfunction") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "energy") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "energyout") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "legendre") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "material") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "mesh") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "meshsurface") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "mu") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "particle") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "polar") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "surface") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "spatiallegendre") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "sphericalharmonics") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "universe") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "zernike") { + model::tally_filters.push_back(std::make_unique()); + } else if (type == "zernikeradial") { + model::tally_filters.push_back(std::make_unique()); + } else { + throw std::runtime_error{"Unknown filter type: " + type}; + } + return model::tally_filters.back().get(); +} + //============================================================================== // Fortran compatibility functions @@ -50,59 +105,6 @@ std::vector> tally_filters; extern "C" { // filter_match_point moved to simulation.cpp - Filter* - allocate_filter(const char* type) - { - std::string type_ {type}; - if (type_ == "azimuthal") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "cell") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "cellborn") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "cellfrom") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "distribcell") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "delayedgroup") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "energyfunction") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "energy") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "energyout") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "legendre") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "material") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "mesh") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "meshsurface") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "mu") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "particle") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "polar") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "surface") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "spatiallegendre") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "sphericalharmonics") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "universe") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "zernike") { - model::tally_filters.push_back(std::make_unique()); - } else if (type_ == "zernikeradial") { - model::tally_filters.push_back(std::make_unique()); - } else { - return nullptr; - } - return model::tally_filters.back().get(); - } int32_t filter_get_id(Filter* filt) {return filt->id_;} diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 31806b5b4e..19aa567bb3 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -3,9 +3,12 @@ #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/error.h" +#include "openmc/file_utils.h" #include "openmc/message_passing.h" +#include "openmc/mesh.h" #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" +#include "openmc/particle.h" #include "openmc/reaction_product.h" #include "openmc/settings.h" #include "openmc/simulation.h" @@ -19,6 +22,8 @@ #include "openmc/tallies/filter_legendre.h" #include "openmc/tallies/filter_mesh.h" #include "openmc/tallies/filter_meshsurface.h" +#include "openmc/tallies/filter_particle.h" +#include "openmc/tallies/filter_sph_harm.h" #include "openmc/tallies/filter_surface.h" #include "openmc/xml_interface.h" @@ -45,6 +50,7 @@ namespace model { std::vector active_collision_tallies; std::vector active_meshsurf_tallies; std::vector active_surface_tallies; + std::unordered_map tally_map; } namespace simulation { @@ -439,7 +445,7 @@ Tally::set_nuclides(pugi::xml_node node) } void -Tally::init_triggers(pugi::xml_node node, int i_tally) +Tally::init_triggers(pugi::xml_node node) { for (auto trigger_node: node.children("trigger")) { // Read the trigger type. @@ -484,19 +490,18 @@ Tally::init_triggers(pugi::xml_node node, int i_tally) } // Parse the trigger scores and populate the triggers_ vector. - const auto& tally {*model::tallies[i_tally]}; for (auto score_str : trigger_scores) { if (score_str == "all") { - triggers_.reserve(triggers_.size() + tally.scores_.size()); - for (auto i_score = 0; i_score < tally.scores_.size(); ++i_score) { + triggers_.reserve(triggers_.size() + this->scores_.size()); + for (auto i_score = 0; i_score < this->scores_.size(); ++i_score) { triggers_.push_back({metric, threshold, i_score}); } } else { int i_score = 0; - for (; i_score < tally.scores_.size(); ++i_score) { - if (reaction_name(tally.scores_[i_score]) == score_str) break; + for (; i_score < this->scores_.size(); ++i_score) { + if (reaction_name(this->scores_[i_score]) == score_str) break; } - if (i_score == tally.scores_.size()) { + if (i_score == this->scores_.size()) { std::stringstream msg; msg << "Could not find the score \"" << score_str << "\" in tally " << id_ << " but it was listed in a trigger on that tally"; @@ -547,6 +552,317 @@ void Tally::accumulate() // Non-member functions //============================================================================== +void read_tallies_xml() +{ + // Check if tallies.xml exists. If not, just return since it is optional + std::string filename = settings::path_input + "tallies.xml"; + if (!file_exists(filename)) return; + + write_message("Reading tallies XML file...", 5); + + // Parse tallies.xml file + pugi::xml_document doc; + doc.load_file(filename.c_str()); + pugi::xml_node root = doc.document_element(); + + // Check for setting + if (check_for_node(root, "assume_separate")) { + settings::assume_separate = get_node_value_bool(root, "assume_separate"); + } + + // Check for user meshes and allocate + read_meshes(root); + + // We only need the mesh info for plotting + if (settings::run_mode == RUN_MODE_PLOTTING) return; + + // Read data for tally derivatives + read_tally_derivatives(root); + + // ========================================================================== + // READ FILTER DATA + + // Check for user filters and allocate + for (auto node_filt : root.children("filter")) { + // Copy filter id + if (!check_for_node(node_filt, "id")) { + fatal_error("Must specify id for filter in tally XML file."); + } + int filter_id = std::stoi(get_node_value(node_filt, "id")); + + // Check to make sure 'id' hasn't been used + if (model::filter_map.find(filter_id) != model::filter_map.end()) { + fatal_error("Two or more filters use the same unique ID: " + + std::to_string(filter_id)); + } + + // Convert filter type to lower case + std::string s; + if (check_for_node(node_filt, "type")) { + s = get_node_value(node_filt, "type", true); + } + + // Make sure bins have been set + if (s == "energy" || s == "energyout" || s == "mu" || s == "polar" + || s == "azimuthal" || s == "mesh" || s == "meshsurface" || s == "universe" + || s == "material" || s == "cell" || s == "distribcell" || s == "cellborn" + || s == "cellfrom" || s == "surface" || s == "delayedgroup") { + if (!check_for_node(node_filt, "bins")) { + fatal_error("Bins not set in filter " + std::to_string(filter_id)); + } + } + + // Allocate according to the filter type + Filter* f = allocate_filter(s); + + // Read filter data from XML + f->from_xml(node_filt); + + // Set filter id + f->id_ = filter_id; + model::filter_map[filter_id] = model::tally_filters.size() - 1; + + // Initialize filter + f->initialize(); + } + + // ========================================================================== + // READ TALLY DATA + + // Check for user tallies + int n = 0; + for (auto node : root.children("tally")) ++n; + if (n == 0 && mpi::master) { + warning("No tallies present in tallies.xml file."); + } + + for (auto node_tal : root.children("tally")) { + model::tallies.push_back(std::make_unique()); + + auto& t {model::tallies.back()}; + t->init_from_xml(node_tal); + + // Copy and set tally id + if (!check_for_node(node_tal, "id")) { + fatal_error("Must specify id for tally in tally XML file."); + } + t->id_ = std::stoi(get_node_value(node_tal, "id")); + model::tally_map[t->id_] = model::tallies.size() - 1; + + // Copy tally name + if (check_for_node(node_tal, "name")) { + t->name_ = get_node_value(node_tal, "name"); + } + + // ======================================================================= + // READ DATA FOR FILTERS + + // Check if user is using old XML format and throw an error if so + if (check_for_node(node_tal, "filter")) { + fatal_error("Tally filters must be specified independently of " + "tallies in a element. The element itself should " + "have a list of filters that apply, e.g., 1 2 " + "where 1 and 2 are the IDs of filters specified outside of " + "."); + } + + // Determine number of filters + std::vector filters; + if (check_for_node(node_tal, "filters")) { + filters = get_node_array(node_tal, "filters"); + } + + // Allocate and store filter user ids + if (!filters.empty()) { + std::vector filter_indices; + for (int filter_id : filters) { + // Determine if filter ID is valid + auto it = model::filter_map.find(filter_id); + if (it == model::filter_map.end()) { + fatal_error("Could not find filter " + std::to_string(filter_id) + + " specified on tally " + std::to_string(t->id_)); + } + + // Store the index of the filter + filter_indices.push_back(it->second); + } + + // Set the filters + t->set_filters(filter_indices.data(), filter_indices.size()); + } + + // Check for the presence of certain filter types + bool has_energyout = t->energyout_filter_ >= 0; + int particle_filter_index = C_NONE; + for (int j = 0; j < t->filters().size(); ++j) { + int i_filter = t->filters(j); + const auto& f = model::tally_filters[i_filter].get(); + + auto pf = dynamic_cast(f); + if (pf) particle_filter_index = j; + + // Change the tally estimator if a filter demands it + std::string filt_type = f->type(); + if (filt_type == "energyout" || filt_type == "legendre") { + t->estimator_ = ESTIMATOR_ANALOG; + } else if (filt_type == "sphericalharmonics") { + auto sf = dynamic_cast(f); + if (sf->cosine_ == SphericalHarmonicsCosine::scatter) { + t->estimator_ = ESTIMATOR_ANALOG; + } + } else if (filt_type == "spatiallegendre" || filt_type == "zernike" + || filt_type == "zernikeradial") { + t->estimator_ = ESTIMATOR_COLLISION; + } + } + + // ======================================================================= + // READ DATA FOR NUCLIDES + + t->set_nuclides(node_tal); + + // ======================================================================= + // READ DATA FOR SCORES + + t->set_scores(node_tal); + + if (!check_for_node(node_tal, "scores")) { + fatal_error("No scores specified on tally " + std::to_string(t->id_) + + "."); + } + + // Check if tally is compatible with particle type + if (settings::photon_transport) { + if (particle_filter_index == C_NONE) { + for (int score : t->scores_) { + switch (score) { + case SCORE_INVERSE_VELOCITY: + fatal_error("Particle filter must be used with photon " + "transport on and inverse velocity score"); + break; + case SCORE_FLUX: + case SCORE_TOTAL: + case SCORE_SCATTER: + case SCORE_NU_SCATTER: + case SCORE_ABSORPTION: + case SCORE_FISSION: + case SCORE_NU_FISSION: + case SCORE_CURRENT: + case SCORE_EVENTS: + case SCORE_DELAYED_NU_FISSION: + case SCORE_PROMPT_NU_FISSION: + case SCORE_DECAY_RATE: + warning("Particle filter is not used with photon transport" + " on and " + std::to_string(score) + " score."); + break; + } + } + } else { + const auto& f = model::tally_filters[particle_filter_index].get(); + auto pf = dynamic_cast(f); + for (int p : pf->particles_) { + if (p == static_cast(ParticleType::electron) || + p == static_cast(ParticleType::positron)) { + t->estimator_ = ESTIMATOR_ANALOG; + } + } + } + } else { + if (particle_filter_index >= 0) { + const auto& f = model::tally_filters[particle_filter_index].get(); + auto pf = dynamic_cast(f); + for (int p : pf->particles_) { + if (p != static_cast(ParticleType::neutron)) { + warning("Particle filter other than NEUTRON used with photon " + "transport turned off. All tallies for particle type " + + std::to_string(p) + " will have no scores"); + } + } + } + } + + // Check for a tally derivative. + if (check_for_node(node_tal, "derivative")) { + int deriv_id = std::stoi(get_node_value(node_tal, "derivative")); + + // Find the derivative with the given id, and store it's index. + auto it = model::tally_deriv_map.find(deriv_id); + if (it == model::tally_deriv_map.end()) { + fatal_error("Could not find derivative " + std::to_string(deriv_id) + + " specified on tally " + std::to_string(t->id_)); + } + + t->deriv_ = it->second; + + // Only analog or collision estimators are supported for differential + // tallies. + if (t->estimator_ == ESTIMATOR_TRACKLENGTH) { + t->estimator_ = ESTIMATOR_COLLISION; + } + + const auto& deriv = model::tally_derivs[t->deriv_]; + if (deriv.variable == DIFF_NUCLIDE_DENSITY + || deriv.variable == DIFF_TEMPERATURE) { + for (int i_nuc : t->nuclides_) { + if (has_energyout && i_nuc == -1) { + fatal_error("Error on tally " + std::to_string(t->id_) + + ": Cannot use a 'nuclide_density' or 'temperature' " + "derivative on a tally with an outgoing energy filter and " + "'total' nuclide rate. Instead, tally each nuclide in the " + "material individually."); + // Note that diff tallies with these characteristics would work + // correctly if no tally events occur in the perturbed material + // (e.g. pertrubing moderator but only tallying fuel), but this + // case would be hard to check for by only reading inputs. + } + } + } + } + + // If settings.xml trigger is turned on, create tally triggers + if (settings::trigger_on) { + t->init_triggers(node_tal); + } + + // ======================================================================= + // SET TALLY ESTIMATOR + + // Check if user specified estimator + if (check_for_node(node_tal, "estimator")) { + std::string est = get_node_value(node_tal, "estimator"); + if (est == "analog") { + t->estimator_ = ESTIMATOR_ANALOG; + } else if (est == "tracklength" || est == "track-length" + || est == "pathlength" || est == "path-length") { + // If the estimator was set to an analog estimator, this means the + // tally needs post-collision information + if (t->estimator_ == ESTIMATOR_ANALOG) { + fatal_error("Cannot use track-length estimator for tally " + + std::to_string(t->id_)); + } + + // Set estimator to track-length estimator + t->estimator_ = ESTIMATOR_TRACKLENGTH; + + } else if (est == "collision") { + // If the estimator was set to an analog estimator, this means the + // tally needs post-collision information + if (t->estimator_ == ESTIMATOR_ANALOG) { + fatal_error("Cannot use collision estimator for tally " + + std::to_string(t->id_)); + } + + // Set estimator to collision estimator + t->estimator_ = ESTIMATOR_COLLISION; + + } else { + fatal_error("Invalid estimator '" + est + "' on tally " + + std::to_string(t->id_)); + } + } + } +} + #ifdef OPENMC_MPI void reduce_tally_results() { @@ -1041,8 +1357,8 @@ extern "C" { void tally_set_nuclides(Tally* tally, pugi::xml_node* node) {tally->set_nuclides(*node);} - void tally_init_triggers(Tally* tally, int i_tally, pugi::xml_node* node) - {tally->init_triggers(*node, i_tally);} + // void tally_init_triggers(Tally* tally, int i_tally, pugi::xml_node* node) + // {tally->init_triggers(*node, i_tally);} int tally_get_deriv_c(Tally* tally) {return tally->deriv_;} diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index 58b6a196e0..9a8eff9d52 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -97,14 +97,6 @@ contains character(:), allocatable :: type_ - interface - function allocate_filter(type) result(ptr) bind(C) - import C_CHAR, C_PTR - character(kind=C_CHAR), intent(in) :: type(*) - type(C_PTR) :: ptr - end function allocate_filter - end interface - ! Convert C string to Fortran string type_ = to_f_string(type) @@ -164,7 +156,7 @@ contains call set_errmsg("Unknown filter type: " // trim(type_)) end select - filters(index) % obj % ptr = allocate_filter(type) + !filters(index) % obj % ptr = allocate_filter(type) if (.not. c_associated(filters(index) % obj % ptr)) then err = E_UNASSIGNED call set_errmsg("Could not allocate C++ tally filter")