diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index fa6460551..cf7fb732f 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -368,6 +368,15 @@ write_dataset(hid_t obj_id, const char* name, const xt::xarray& arr) arr.data(), false); } +template inline void +write_dataset(hid_t obj_id, const char* name, const xt::xtensor& arr) +{ + auto s = arr.shape(); + std::vector dims {s.cbegin(), s.cend()}; + write_dataset(obj_id, dims.size(), dims.data(), name, H5TypeMap::type_id, + arr.data(), false); +} + inline void write_dataset(hid_t obj_id, const char* name, Position r) { diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 64713ae96..6541a67fc 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -14,5 +14,7 @@ extern "C" void write_source_bank(hid_t group_id, int64_t* work_index, extern "C" void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank); +extern "C" void write_tally_results_nr(hid_t file_id); + } // namespace openmc #endif // OPENMC_STATE_POINT_H diff --git a/src/state_point.F90 b/src/state_point.F90 index 91d0f1090..11e4ad621 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -87,6 +87,10 @@ contains import HID_T integer(HID_T), value :: group end subroutine + subroutine write_tally_results_nr(file_id) bind(C) + import HID_T + integer(HID_T), value :: file_id + end subroutine end interface err = 0 @@ -490,129 +494,6 @@ contains end subroutine write_source_point -!=============================================================================== -! WRITE_TALLY_RESULTS_NR -!=============================================================================== - - subroutine write_tally_results_nr(file_id) - integer(HID_T), intent(in) :: file_id - - integer :: i ! loop index - integer :: n ! number of filter bins - integer :: m ! number of score bins - integer :: n_bins ! total number of bins - integer(HID_T) :: tallies_group, tally_group - real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results - real(8), target :: global_temp(3,N_GLOBAL_TALLIES) -#ifdef OPENMC_MPI - integer :: mpi_err ! MPI error code - real(8) :: dummy ! temporary receive buffer for non-root reduces -#endif - type(TallyObject) :: dummy_tally - - ! ========================================================================== - ! COLLECT AND WRITE GLOBAL TALLIES - - if (master) then - ! Write number of realizations - call write_dataset(file_id, "n_realizations", n_realizations) - - ! Write number of global tallies - call write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES) - - tallies_group = open_group(file_id, "tallies") - end if - - -#ifdef OPENMC_MPI - ! Reduce global tallies - n_bins = size(global_tallies) - call MPI_REDUCE(global_tallies, global_temp, n_bins, MPI_REAL8, MPI_SUM, & - 0, mpi_intracomm, mpi_err) -#endif - - if (master) then - ! Transfer values to value on master - if (current_batch == n_max_batches .or. satisfy_triggers) then - global_tallies(:,:) = global_temp(:,:) - end if - - ! Write out global tallies sum and sum_sq - call write_dataset(file_id, "global_tallies", global_temp) - end if - - if (active_tallies % size() > 0) then - ! Indicate that tallies are on - if (master) then - call write_attribute(file_id, "tallies_present", 1) - end if - - ! Write all tally results - TALLY_RESULTS: do i = 1, n_tallies - associate (t => tallies(i) % obj) - ! Determine size of tally results array - m = size(t % results, 2) - n = size(t % results, 3) - n_bins = m*n*2 - - ! Allocate array for storing sums and sums of squares, but - ! contiguously in memory for each - allocate(tally_temp(2,m,n)) - tally_temp(1,:,:) = t % results(RESULT_SUM,:,:) - tally_temp(2,:,:) = t % results(RESULT_SUM_SQ,:,:) - - if (master) then - tally_group = open_group(tallies_group, "tally " // & - trim(to_str(t % id))) - - ! The MPI_IN_PLACE specifier allows the master to copy values into - ! a receive buffer without having a temporary variable -#ifdef OPENMC_MPI - call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & - MPI_SUM, 0, mpi_intracomm, mpi_err) -#endif - - ! At the end of the simulation, store the results back in the - ! regular TallyResults array - if (current_batch == n_max_batches .or. satisfy_triggers) then - t % results(RESULT_SUM,:,:) = tally_temp(1,:,:) - t % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:) - end if - - ! Put in temporary tally result - allocate(dummy_tally % results(3,m,n)) - dummy_tally % results(RESULT_SUM,:,:) = tally_temp(1,:,:) - dummy_tally % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:) - - ! Write reduced tally results to file - call dummy_tally % write_results_hdf5(tally_group) - - ! Deallocate temporary tally result - deallocate(dummy_tally % results) - else - ! Receive buffer not significant at other processors -#ifdef OPENMC_MPI - call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, & - 0, mpi_intracomm, mpi_err) -#endif - end if - - ! Deallocate temporary copy of tally results - deallocate(tally_temp) - - if (master) call close_group(tally_group) - end associate - end do TALLY_RESULTS - - if (master) call close_group(tallies_group) - else - ! Indicate that tallies are off - if (master) call write_dataset(file_id, "tallies_present", 0) - end if - - - end subroutine write_tally_results_nr - !=============================================================================== ! LOAD_STATE_POINT !=============================================================================== diff --git a/src/state_point.cpp b/src/state_point.cpp index e24993b6c..4bd91e121 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -1,17 +1,20 @@ #include "openmc/state_point.h" -#include -#include - -#ifdef OPENMC_MPI -#include "mpi.h" -#endif - #include "openmc/capi.h" +#include "openmc/constants.h" #include "openmc/error.h" +#include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/tallies/tally.h" + +#include "xtensor/xbuilder.hpp" // for empty_like +#include "xtensor/xview.hpp" + +#include +#include +#include namespace openmc { @@ -165,4 +168,110 @@ void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank) H5Tclose(banktype); } +void write_tally_results_nr(hid_t file_id) +{ + // ========================================================================== + // COLLECT AND WRITE GLOBAL TALLIES + + hid_t tallies_group; + if (mpi::master) { + // Write number of realizations + write_dataset(file_id, "n_realizations", n_realizations); + + // Write number of global tallies + write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES); + + tallies_group = open_group(file_id, "tallies"); + } + + // Get pointer to global tallies + auto gt = global_tallies(); + +#ifdef OPENMC_MPI + // Reduce global tallies + xt::xtensor gt_reduced = xt::empty_like(gt); + MPI_Reduce(gt.data(), gt_reduced.data(), gt.size(), MPI_DOUBLE, + MPI_SUM, 0, mpi::intracomm); + + // Transfer values to value on master + if (mpi::master) { + if (simulation::current_batch == settings::n_max_batches || + simulation::satisfy_triggers) { + std::copy(gt_reduced.begin(), gt_reduced.end(), gt.begin()); + } + } +#endif + + // Write out global tallies sum and sum_sq + if (mpi::master) { + write_dataset(file_id, "global_tallies", gt); + } + + for (int i = 1; i <= n_tallies; ++i) { + // Skip any tallies that are not active + bool active; + openmc_tally_get_active(i, &active); + if (!active) continue; + + if (mpi::master && !object_exists(file_id, "tallies_present")) { + write_attribute(file_id, "tallies_present", 1); + } + + // Get view of accumulated tally values + auto results = tally_results(i); + auto values_view = xt::view(results, xt::all(), xt::all(), + xt::range(RESULT_SUM, RESULT_SUM_SQ + 1)); + + // Make copy of tally values in contiguous array + xt::xtensor values = values_view; + + if (mpi::master) { + // Open group for tally + int id; + openmc_tally_get_id(i, &id); + std::string groupname {"tally " + std::to_string(id)}; + hid_t tally_group = open_group(tallies_group, groupname.c_str()); + + // The MPI_IN_PLACE specifier allows the master to copy values into + // a receive buffer without having a temporary variable +#ifdef OPENMC_MPI + MPI_Reduce(MPI_IN_PLACE, values.data(), values.size(), MPI_DOUBLE, + MPI_SUM, 0, mpi::intracomm); +#endif + + // At the end of the simulation, store the results back in the + // regular TallyResults array + if (simulation::current_batch == settings::n_max_batches || + simulation::satisfy_triggers) { + values_view = values; + } + + // Put in temporary tally result + xt::xtensor results_copy = xt::zeros_like(results); + auto copy_view = xt::view(results_copy, xt::all(), xt::all(), + xt::range(RESULT_SUM, RESULT_SUM_SQ + 1)); + copy_view = values; + + // Write reduced tally results to file + auto shape = results_copy.shape(); + write_tally_results(tally_group, shape[0], shape[1], results_copy.data()); + + close_group(tally_group); + } else { + // Receive buffer not significant at other processors +#ifdef OPENMC_MPI + MPI_Reduce(values.data(), nullptr, values.size(), MPI_REAL8, MPI_SUM, + 0, mpi::intracomm); +#endif + } + } + + if (mpi::master) { + if (!object_exists(file_id, "tallies_present")) { + // Indicate that tallies are off + write_dataset(file_id, "tallies_present", 0); + } + } +} + } // namespace openmc