From 9f6a7466e24f6bd309bfb37265965eb7180402cd Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 26 Jan 2019 13:53:01 -0500 Subject: [PATCH] Use interoperable TallyDerivative structs --- include/openmc/tallies/derivative.h | 2 +- include/openmc/tallies/tally.h | 4 +- src/api.F90 | 2 - src/input_xml.F90 | 49 ++++------- src/output.F90 | 8 +- src/state_point.F90 | 17 ++-- src/tallies/derivative.cpp | 35 ++++++-- src/tallies/tally.F90 | 37 +++++---- src/tallies/tally.cpp | 6 ++ src/tallies/tally_derivative_header.F90 | 106 +++++------------------- src/tallies/tally_header.F90 | 31 ++++++- 11 files changed, 139 insertions(+), 158 deletions(-) diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index 585e9702c..d9202433c 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -12,7 +12,7 @@ namespace openmc { -struct TallyDerivative { +extern "C" struct TallyDerivative { int id; //!< User-defined identifier int variable; //!< Independent variable (like temperature) int diff_material; //!< Material this derivative is applied to diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 5f4deb34c..77e5958ef 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -39,8 +39,10 @@ public: int surface_filter_ {C_NONE}; int mesh_filter_ {C_NONE}; + int deriv_ {F90_NONE}; //!< Index of a TallyDerivative object for diff tallies. + private: - std::vector filters_; //< Filter indices in global filters array + std::vector filters_; //!< Filter indices in global filters array //! Index strides assigned to each filter to support 1D indexing. std::vector strides_; diff --git a/src/api.F90 b/src/api.F90 index 640ad910a..cf5670bc9 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -66,7 +66,6 @@ contains use settings use simulation_header use surface_header - use tally_derivative_header use tally_filter_header use tally_header use trigger_header @@ -102,7 +101,6 @@ contains call free_memory_mesh() call free_memory_tally() call free_memory_tally_filter() - call free_memory_tally_derivative() call free_memory_bank() call free_memory_cmfd() #ifdef DAGMC diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 181e52843..87eeac6c3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -913,6 +913,7 @@ contains 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 :: i_elem ! index of entry in dictionary integer :: n ! size of arrays in mesh specification @@ -948,6 +949,7 @@ contains type(XMLNode), allocatable :: node_trigger_list(:) type(XMLNode), allocatable :: node_deriv_list(:) type(DictEntryCI) :: elem + type(TallyDerivative), pointer :: deriv interface subroutine read_tally_derivatives(node_ptr) bind(C) @@ -960,12 +962,6 @@ contains filename = trim(path_input) // "tallies.xml" inquire(FILE=filename, EXIST=file_exists) if (.not. file_exists) then - ! We need to allocate tally_derivs to avoid segfaults. Also needs to be - ! done in parallel because tally derivs are threadprivate. -!$omp parallel - allocate(tally_derivs(0)) -!$omp end parallel - ! Since a tallies.xml file is optional, no error is issued here return end if @@ -1008,28 +1004,10 @@ contains call read_tally_derivatives(root % ptr) - ! Get pointer list to XML nodes and allocate global array. - ! The array is threadprivate so it must be allocated in parallel. call get_node_list(root, "derivative", node_deriv_list) -!$omp parallel - allocate(tally_derivs(size(node_deriv_list))) -!$omp end parallel - - ! Make sure this is not an MG run. - if (.not. run_CE .and. size(node_deriv_list) > 0) then - call fatal_error("Differential tallies not supported in multi-group mode") - end if - - ! Read derivative attributes. do i = 1, size(node_deriv_list) -!$omp parallel -!$omp critical (ReadTallyDeriv) - call tally_derivs(i) % from_xml(node_deriv_list(i)) -!$omp end critical (ReadTallyDeriv) -!$omp end parallel - - ! Update tally derivative dictionary - call tally_deriv_dict % set(tally_derivs(i) % id, i) + deriv => tally_deriv_c(i - 1) + call tally_deriv_dict % set(deriv % id, i) end do ! ========================================================================== @@ -1618,13 +1596,13 @@ contains ! Check for a tally derivative. if (check_for_node(node_tal, "derivative")) then - ! Temporarily store the derivative id. - call get_node_value(node_tal, "derivative", t % deriv) + call get_node_value(node_tal, "derivative", deriv_id) ! Find the derivative with the given id, and store it's index. - do j = 1, size(tally_derivs) - if (tally_derivs(j) % id == t % deriv) then - t % deriv = j + do j = 1, n_tally_derivs() + deriv => tally_deriv_c(j - 1) + 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 @@ -1633,15 +1611,16 @@ contains ! We found the derivative we were looking for; exit the do loop. exit end if - if (j == size(tally_derivs)) then + if (j == n_tally_derivs()) then call fatal_error("Could not find derivative " & - // trim(to_str(t % deriv)) // " specified on tally " & + // trim(to_str(deriv_id)) // " specified on tally " & // trim(to_str(t % id))) end if end do - if (tally_derivs(t % deriv) % variable == DIFF_NUCLIDE_DENSITY & - .or. tally_derivs(t % deriv) % variable == DIFF_TEMPERATURE) then + deriv => tally_deriv_c(t % deriv() - 1) + if (deriv % variable == DIFF_NUCLIDE_DENSITY & + .or. deriv % variable == DIFF_TEMPERATURE) then if (any(t % nuclide_bins == -1)) then if (has_energyout) then call fatal_error("Error on tally " // trim(to_str(t % id)) & diff --git a/src/output.F90 b/src/output.F90 index 6efd2c6e6..04e35c17d 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -471,6 +471,7 @@ contains ! to be applied at write-time integer, allocatable :: filter_bins(:) character(MAX_WORD_LEN) :: temp_name + type(TallyDerivative), pointer :: deriv ! Skip if there are no tallies if (n_tallies == 0) return @@ -530,8 +531,9 @@ contains endif ! Write derivative information. - if (t % deriv /= NONE) then - associate(deriv => tally_derivs(t % deriv)) + if (t % deriv() /= NONE) then + !associate(deriv => tally_derivs(t % deriv())) + deriv => tally_deriv_c(t % deriv() - 1) select case (deriv % variable) case (DIFF_DENSITY) write(unit=unit_tally, fmt="(' Density derivative Material ',A)") & @@ -548,7 +550,7 @@ contains call fatal_error("Differential tally dependent variable for tally "& // trim(to_str(t % id)) // " not defined in output.F90.") end select - end associate + !end associate end if ! WARNING: Admittedly, the logic for moving for printing results is diff --git a/src/state_point.F90 b/src/state_point.F90 index fba79dd84..f9fde0191 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -28,7 +28,7 @@ module state_point use string, only: to_str, count_digits, zero_padded, to_f_string use tally_header use tally_filter_header - use tally_derivative_header, only: tally_derivs + use tally_derivative_header use timer_header implicit none @@ -69,6 +69,7 @@ contains character(len=:, kind=C_CHAR), allocatable :: filename_ character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name logical :: parallel + type(TallyDerivative), pointer :: deriv interface subroutine meshes_to_hdf5(group) bind(C) @@ -174,10 +175,11 @@ contains call meshes_to_hdf5(tallies_group) ! Write information for derivatives. - if (size(tally_derivs) > 0) then + if (n_tally_derivs() > 0) then derivs_group = create_group(tallies_group, "derivatives") - do i = 1, size(tally_derivs) - associate(deriv => tally_derivs(i)) + do i = 1, n_tally_derivs() + !associate(deriv => tally_derivs(i)) + deriv => tally_deriv_c(i - 1) deriv_group = create_group(derivs_group, "derivative " & // trim(to_str(deriv % id))) select case (deriv % variable) @@ -200,7 +202,7 @@ contains &state_point.F90.") end select call close_group(deriv_group) - end associate + !end associate end do call close_group(derivs_group) end if @@ -303,9 +305,10 @@ contains deallocate(str_array) ! Write derivative information. - if (tally % deriv /= NONE) then + if (tally % deriv() /= NONE) then + deriv => tally_deriv_c(tally % deriv() - 1) call write_dataset(tally_group, "derivative", & - tally_derivs(tally % deriv) % id) + deriv % id) end if ! Write scores. diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index c78be5973..66889708a 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -33,12 +33,12 @@ TallyDerivative::TallyDerivative(pugi::xml_node node) if (id <= 0) fatal_error(" IDs must be an integer greater than zero"); - std::string variable = get_node_value(node, "variable"); + std::string variable_str = get_node_value(node, "variable"); - if (variable == "density") { + if (variable_str == "density") { variable = DIFF_DENSITY; - } else if (variable == "nuclide_density") { + } else if (variable_str == "nuclide_density") { variable = DIFF_NUCLIDE_DENSITY; std::string nuclide_name = get_node_value(node, "nuclide"); @@ -46,7 +46,8 @@ TallyDerivative::TallyDerivative(pugi::xml_node node) for (auto i = 0; i < data::nuclides.size(); ++i) { if (data::nuclides[i]->name_ == nuclide_name) { found = true; - diff_nuclide = i; + //TODO: off-by-one + diff_nuclide = i + 1; } } if (!found) { @@ -56,12 +57,13 @@ TallyDerivative::TallyDerivative(pugi::xml_node node) fatal_error(out); } - } else if (variable == "temperature") { + } else if (variable_str == "temperature") { variable = DIFF_TEMPERATURE; } else { std::stringstream out; - out << "Unrecognized variable \"" << variable << "\" on derivative " << id; + out << "Unrecognized variable \"" << variable_str + << "\" on derivative " << id; fatal_error(out); } @@ -100,4 +102,25 @@ read_tally_derivatives(pugi::xml_node* node) fatal_error("Differential tallies not supported in multi-group mode"); } +/* +//! Set the flux derivatives on differential tallies to zero. +extern "C" void +zero_flux_derivs_c() +{ + for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.; +} +*/ + +//============================================================================== +// Fortran interop +//============================================================================== + +extern "C" int n_tally_derivs() {return model::tally_derivs.size();} + +extern "C" TallyDerivative* +tally_deriv_c(int i) +{ + return &(model::tally_derivs[i]); +} + }// namespace openmc diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 2b274364e..3122bf866 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -17,7 +17,7 @@ module tally use settings use simulation_header use string, only: to_str - use tally_derivative_header, only: tally_derivs + use tally_derivative_header use tally_filter use tally_header @@ -1184,7 +1184,7 @@ contains !######################################################################### ! Add derivative information on score for differential tallies. - if (t % deriv /= NONE) then + if (t % deriv() /= NONE) then call apply_derivative_to_score(p, t, i_nuclide, atom_density, & score_bin, score) end if @@ -2375,7 +2375,7 @@ contains ! Add derivative information for differential tallies. Note that the ! i_nuclide and atom_density arguments do not matter since this is an ! analog estimator. - if (t % deriv /= NONE) then + if (t % deriv() /= NONE) then call apply_derivative_to_score(p, t, 0, ZERO, SCORE_NU_FISSION, score) end if @@ -2992,6 +2992,7 @@ contains integer, intent(in) :: score_bin real(8), intent(inout) :: score + type(TallyDerivative), pointer :: deriv integer :: l logical :: scoring_diff_nuclide real(8) :: flux_deriv @@ -3004,10 +3005,12 @@ contains ! where (1/f * d_f/d_p) is the (logarithmic) flux derivative and p is the ! perturbated variable. - associate(deriv => tally_derivs(t % deriv)) + !associate(deriv => tally_derivs(t % deriv())) + deriv => tally_deriv_c(t % deriv() - 1) flux_deriv = deriv % flux_deriv - select case (tally_derivs(t % deriv) % variable) + !select case (tally_derivs(t % deriv()) % variable) + select case (deriv % variable) !========================================================================= ! Density derivative: @@ -3546,7 +3549,7 @@ contains &analog and collision estimators.") end select end select - end associate + !end associate end subroutine apply_derivative_to_score !=============================================================================== @@ -3558,14 +3561,16 @@ contains type(Particle), intent(in) :: p real(8), intent(in) :: distance ! Neutron flight distance + type(TallyDerivative), pointer :: deriv integer :: i, l real(8) :: dsig_s, dsig_a, dsig_f ! A void material cannot be perturbed so it will not affect flux derivatives if (p % material == MATERIAL_VOID) return - do i = 1, size(tally_derivs) - associate(deriv => tally_derivs(i)) + do i = 1, n_tally_derivs() + !associate(deriv => tally_derivs(i)) + deriv => tally_deriv_c(i - 1) select case (deriv % variable) case (DIFF_DENSITY) @@ -3609,7 +3614,7 @@ contains end if end associate end select - end associate + !end associate end do end subroutine score_track_derivative @@ -3632,14 +3637,16 @@ contains subroutine score_collision_derivative(p) type(Particle), intent(in) :: p + type(TallyDerivative), pointer :: deriv integer :: i, j, l real(8) :: dsig_s, dsig_a, dsig_f ! A void material cannot be perturbed so it will not affect flux derivatives if (p % material == MATERIAL_VOID) return - do i = 1, size(tally_derivs) - associate(deriv => tally_derivs(i)) + do i = 1, n_tally_derivs() + !associate(deriv => tally_derivs(i)) + deriv => tally_deriv_c(i - 1) select case (deriv % variable) case (DIFF_DENSITY) @@ -3702,7 +3709,7 @@ contains end if end associate end select - end associate + !end associate end do end subroutine score_collision_derivative @@ -3712,8 +3719,10 @@ contains subroutine zero_flux_derivs() integer :: i - do i = 1, size(tally_derivs) - tally_derivs(i) % flux_deriv = ZERO + type(TallyDerivative), pointer :: deriv + do i = 1, n_tally_derivs() + deriv => tally_deriv_c(i - 1) + deriv % flux_deriv = ZERO end do end subroutine zero_flux_derivs diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 38ead0f08..169e886d6 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -3,6 +3,7 @@ #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/error.h" +#include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/filter_energy.h" #include "openmc/tallies/filter_delayedgroup.h" @@ -175,6 +176,7 @@ free_memory_tally_c() #pragma omp parallel { simulation::filter_matches.clear(); + model::tally_derivs.clear(); } model::tally_filters.clear(); @@ -260,6 +262,10 @@ extern "C" { int tally_get_mesh_filter_c(Tally* tally) {return tally->mesh_filter_;} + + int tally_get_deriv_c(Tally* tally) {return tally->deriv_;} + + int tally_set_deriv_c(Tally* tally, int deriv) {tally->deriv_ = deriv;} } } // namespace openmc diff --git a/src/tallies/tally_derivative_header.F90 b/src/tallies/tally_derivative_header.F90 index 304ae5e45..efeddbb36 100644 --- a/src/tallies/tally_derivative_header.F90 +++ b/src/tallies/tally_derivative_header.F90 @@ -1,104 +1,38 @@ module tally_derivative_header - use constants - use dict_header, only: DictIntInt, EMPTY - use error, only: fatal_error - use nuclide_header, only: nuclide_dict - use string, only: to_str, to_lower - use xml_interface + use, intrinsic :: ISO_C_BINDING + + use dict_header, only: DictIntInt implicit none - private - public :: free_memory_tally_derivative !=============================================================================== ! TALLYDERIVATIVE describes a first-order derivative that can be applied to ! tallies. !=============================================================================== - type, public :: TallyDerivative - integer :: id - integer :: variable - integer :: diff_material - integer :: diff_nuclide - real(8) :: flux_deriv - contains - procedure :: from_xml + type, bind(C), public :: TallyDerivative + integer(C_INT) :: id + integer(C_INT) :: variable + integer(C_INT) :: diff_material + integer(C_INT) :: diff_nuclide + real(C_DOUBLE) :: flux_deriv end type TallyDerivative - type(TallyDerivative), public, allocatable :: tally_derivs(:) -!$omp threadprivate(tally_derivs) - ! Dictionary that maps user IDs to indices in 'tally_derivs' type(DictIntInt), public :: tally_deriv_dict -contains + interface + function n_tally_derivs() result(n) bind(C) + import C_INT + integer(C_INT) :: n + end function - subroutine from_xml(this, node) - class(TallyDerivative), intent(inout) :: this - type(XMLNode), intent(in) :: node - - character(MAX_WORD_LEN) :: temp_str - character(MAX_WORD_LEN) :: word - integer :: val - - ! Copy the derivative id. - if (check_for_node(node, "id")) then - call get_node_value(node, "id", this % id) - else - call fatal_error("Must specify an ID for elements in the& - & tally XML file") - end if - - ! Make sure the id is > 0. - if (this % id <= 0) then - call fatal_error(" IDs must be an integer greater than & - &zero") - end if - - ! Make sure this id has not already been used. - if (tally_deriv_dict % has(this % id)) then - call fatal_error("Two or more 's use the same unique & - &ID: " // trim(to_str(this % id))) - end if - - ! Read the independent variable name. - call get_node_value(node, "variable", temp_str) - temp_str = to_lower(temp_str) - - select case(temp_str) - case("density") - this % variable = DIFF_DENSITY - - case("nuclide_density") - this % variable = DIFF_NUCLIDE_DENSITY - - call get_node_value(node, "nuclide", word) - word = trim(to_lower(word)) - val = nuclide_dict % get(word) - if (val == EMPTY) then - call fatal_error("Could not find the nuclide " & - // trim(word) // " specified in derivative " & - // trim(to_str(this % id)) // " in any material.") - end if - this % diff_nuclide = val - - case("temperature") - this % variable = DIFF_TEMPERATURE - end select - - call get_node_value(node, "material", this % diff_material) - - end subroutine from_xml - -!=============================================================================== -! FREE_MEMORY_TALLY_DERIVATIVE deallocates global arrays defined in this module -!=============================================================================== - - subroutine free_memory_tally_derivative() -!$omp parallel - if (allocated(tally_derivs)) deallocate(tally_derivs) -!$omp end parallel - end subroutine free_memory_tally_derivative + function tally_deriv_c(i) result(deriv) bind(C) + import C_INT, TallyDerivative + integer(C_INT), value :: i + type(TallyDerivative), pointer :: deriv + end function + end interface end module tally_derivative_header diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 0b216e3ad..f4300789e 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -93,9 +93,6 @@ module tally_header integer :: n_triggers = 0 ! # of triggers type(TriggerObject), allocatable :: triggers(:) ! Array of triggers - ! Index for the TallyDerivative for differential tallies. - integer :: deriv = NONE - contains procedure :: accumulate => tally_accumulate procedure :: allocate_results => tally_allocate_results @@ -110,6 +107,8 @@ module tally_header procedure :: delayedgroup_filter => tally_get_delayedgroup_filter procedure :: surface_filter => tally_get_surface_filter procedure :: mesh_filter => tally_get_mesh_filter + procedure :: deriv => tally_get_deriv + procedure :: set_deriv => tally_set_deriv end type TallyObject type, public :: TallyContainer @@ -404,6 +403,32 @@ contains filt = tally_get_mesh_filter_c(this % ptr) end function + function tally_get_deriv(this) result(deriv) + class(TallyObject) :: this + integer(C_INT) :: deriv + interface + function tally_get_deriv_c(tally) result(deriv) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: tally + integer(C_INT) :: deriv + end function + end interface + deriv = tally_get_deriv_c(this % ptr) + end function + + subroutine tally_set_deriv(this, deriv) + class(TallyObject) :: this + integer(C_INT) :: deriv + interface + subroutine tally_set_deriv_c(tally, deriv) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: tally + integer(C_INT), value :: deriv + end subroutine + end interface + call tally_set_deriv_c(this % ptr, deriv) + end subroutine + !=============================================================================== ! CONFIGURE_TALLIES initializes several data structures related to tallies. This ! is called after the basic tally data has already been read from the