From 68e65ce2d1ecc42af3307e06dc7303d8761e90ec Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 6 Feb 2019 21:47:08 -0600 Subject: [PATCH 1/7] Translate volume calculations to C++ --- CMakeLists.txt | 3 +- include/openmc/hdf5_interface.h | 8 + include/openmc/volume_calc.h | 61 ++++ src/api.F90 | 4 +- src/input_xml.F90 | 32 -- src/settings.cpp | 9 +- src/volume_calc.F90 | 514 -------------------------------- src/volume_calc.cpp | 425 ++++++++++++++++++++++++++ src/volume_header.F90 | 69 ----- 9 files changed, 503 insertions(+), 622 deletions(-) create mode 100644 include/openmc/volume_calc.h delete mode 100644 src/volume_calc.F90 create mode 100644 src/volume_calc.cpp delete mode 100644 src/volume_header.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 3285041995..48e9f622b5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -342,8 +342,6 @@ add_library(libopenmc SHARED src/tracking.F90 src/track_output.F90 src/vector_header.F90 - src/volume_calc.F90 - src/volume_header.F90 src/xml_interface.F90 src/tallies/tally.F90 src/tallies/tally_derivative_header.F90 @@ -436,6 +434,7 @@ add_library(libopenmc SHARED src/tallies/tally.cpp src/timer.cpp src/thermal.cpp + src/volume_calc.cpp src/wmp.cpp src/xml_interface.cpp src/xsdata.cpp) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 398f03b09d..1fc378bdeb 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -391,6 +391,14 @@ write_attribute(hid_t obj_id, const char* name, const std::vector& buffer) write_attr(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data()); } +inline void +write_attribute(hid_t obj_id, const char* name, Position r) +{ + std::array buffer {r.x, r.y, r.z}; + write_attribute(obj_id, name, buffer); +} + + //============================================================================== // Templates/overloads for write_dataset diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h new file mode 100644 index 0000000000..5f905b8fa4 --- /dev/null +++ b/include/openmc/volume_calc.h @@ -0,0 +1,61 @@ +#ifndef VOLUME_CALC_H +#define VOLUME_CALC_H + +#include "openmc/position.h" + +#include "pugixml.hpp" +#include "xtensor/xtensor.hpp" + +#include +#include + +namespace openmc { + +//============================================================================== +// Volume calculation class +//============================================================================== + +class VolumeCalculation { +public: + // Aliases, types + using int_2dvec = std::vector>; + + struct Result { + std::array volume; //!< Mean/standard deviation of volume + std::vector nuclides; //!< Index of nuclides + std::vector atoms; //!< Number of atoms for each nuclide + std::vector uncertainty; //!< Uncertainty on number of atoms + }; // Results for a single domain + + // Constructors + VolumeCalculation() = default; + VolumeCalculation(pugi::xml_node node); + + // Methods + std::vector execute() const; + void write_volume(const std::string& filename, const std::vector& results) const; + + // Data members + int domain_type_; //!< Type of domain (cell, material, etc.) + int n_samples_; //!< Number of samples to use + Position lower_left_; //!< Lower-left position of bounding box + Position upper_right_; //!< Upper-right position of bounding box + std::vector domain_ids_; //!< IDs of domains to find volumes of + +private: + void check_hit(int i_domain, int i_material, int_2dvec& indices, + int_2dvec& hits) const; + +}; + +//============================================================================== +// Global variables +//============================================================================== + +namespace model { +extern std::vector volume_calcs; +} + +} // namespace openmc + +#endif // VOLUME_CALC_H diff --git a/src/api.F90 b/src/api.F90 index 08112d752f..1412e38e22 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -69,7 +69,6 @@ contains use tally_filter_header use tally_header use trigger_header - use volume_header interface subroutine free_memory_source() bind(C) @@ -90,6 +89,9 @@ contains subroutine free_memory_cmfd() bind(C) end subroutine free_memory_cmfd + subroutine free_memory_volume() bind(C) + end subroutine + subroutine sab_clear() bind(C) end subroutine end interface diff --git a/src/input_xml.F90 b/src/input_xml.F90 index bf661fcdbc..faefbdcddc 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -30,7 +30,6 @@ module input_xml use tally_filter_header use tally_filter use trigger_header - use volume_header use xml_interface implicit none @@ -57,9 +56,6 @@ module input_xml type(C_PTR) :: node_ptr end subroutine read_lattices - subroutine read_settings_xml() bind(C) - end subroutine read_settings_xml - subroutine read_materials(node_ptr) bind(C) import C_PTR type(C_PTR) :: node_ptr @@ -91,34 +87,6 @@ module input_xml contains -!=============================================================================== -! READ_SETTINGS_XML reads data from a settings.xml file and parses it, checking -! for errors and placing properly-formatted data in the right data structures -!=============================================================================== - - subroutine read_settings_xml_f(root_ptr) bind(C) - type(C_PTR), value :: root_ptr - - integer :: i - integer :: n - type(XMLNode) :: root - type(XMLNode) :: node_vol - type(XMLNode), allocatable :: node_vol_list(:) - - ! Get proper XMLNode type given pointer - root % ptr = root_ptr - - call get_node_list(root, "volume_calc", node_vol_list) - n = size(node_vol_list) - allocate(volume_calcs(n)) - do i = 1, n - node_vol = node_vol_list(i) - call volume_calcs(i) % from_xml(node_vol) - end do - - end subroutine read_settings_xml_f - - #ifdef DAGMC !=============================================================================== diff --git a/src/settings.cpp b/src/settings.cpp index 3d40eaea30..479a85624d 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -25,6 +25,7 @@ #include "openmc/simulation.h" #include "openmc/source.h" #include "openmc/string_utils.h" +#include "openmc/volume_calc.h" #include "openmc/xml_interface.h" namespace openmc { @@ -727,7 +728,10 @@ void read_settings_xml() } } - // TODO: Get volume calculations + // Get volume calculations + for (pugi::xml_node node_vol : root.children("volume_calc")) { + model::volume_calcs.emplace_back(node_vol); + } // Get temperature settings if (check_for_node(root, "temperature_default")) { @@ -782,9 +786,6 @@ void read_settings_xml() create_fission_neutrons = get_node_value_bool(root, "create_fission_neutrons"); } } - - // Read remaining settings from Fortran side - read_settings_xml_f(root.internal_object()); } //============================================================================== diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 deleted file mode 100644 index 55c872cb33..0000000000 --- a/src/volume_calc.F90 +++ /dev/null @@ -1,514 +0,0 @@ -module volume_calc - - use, intrinsic :: ISO_C_BINDING - -#ifdef _OPENMP - use omp_lib -#endif - - use constants - use error, only: write_message - use geometry, only: find_cell - use geometry_header, only: cells, universe_id - use hdf5_interface, only: file_open, file_close, write_attribute, & - create_group, close_group, write_dataset, HID_T - use output, only: header, time_stamp - use material_header - use message_passing - use nuclide_header, only: nuclides - use particle_header - use random_lcg, only: prn, prn_set_stream, set_particle_seed - use settings, only: path_output - use stl_vector, only: VectorInt, VectorReal - use string, only: to_str - use timer_header, only: Timer - use volume_header - - implicit none - private - - public :: openmc_calculate_volumes - -contains - -!=============================================================================== -! OPENMC_CALCULATE_VOLUMES runs each of the stochastic volume calculations that -! the user has specified and writes results to HDF5 files -!=============================================================================== - - function openmc_calculate_volumes() result(err) bind(C) - integer :: i, j - integer :: n - integer(C_INT) :: err - real(8), allocatable :: volume(:,:) ! volume mean/stdev in each domain - character(10) :: domain_type - character(MAX_FILE_LEN) :: filename ! filename for HDF5 file - type(Timer) :: time_volume ! timer for volume calculation - type(VectorInt), allocatable :: nuclide_vec(:) ! indices in nuclides array - type(VectorReal), allocatable :: atoms_vec(:) ! total # of atoms of each nuclide - type(VectorReal), allocatable :: uncertainty_vec(:) ! uncertainty of total # of atoms - - if (master) then - call header("STOCHASTIC VOLUME CALCULATION", 3) - call time_volume % start() - end if - - do i = 1, size(volume_calcs) - n = size(volume_calcs(i) % domain_id) - allocate(nuclide_vec(n)) - allocate(atoms_vec(n), uncertainty_vec(n)) - allocate(volume(2,n)) - - if (master) then - call write_message("Running volume calculation " // trim(to_str(i)) & - // "...", 4) - end if - - call get_volume(volume_calcs(i), volume, nuclide_vec, atoms_vec, & - uncertainty_vec) - - if (master) then - select case (volume_calcs(i) % domain_type) - case (FILTER_CELL) - domain_type = ' Cell' - case (FILTER_MATERIAL) - domain_type = ' Material' - case (FILTER_UNIVERSE) - domain_type = ' Universe' - end select - - ! Display domain volumes - do j = 1, size(volume_calcs(i) % domain_id) - call write_message(trim(domain_type) // " " // trim(to_str(& - volume_calcs(i) % domain_id(j))) // ": " // trim(to_str(& - volume(1,j))) // " +/- " // trim(to_str(volume(2,j))) // & - " cm^3", 4) - end do - call write_message("", 4) - - filename = trim(path_output) // 'volume_' // trim(to_str(i)) // '.h5' - call write_volume(volume_calcs(i), filename, volume, nuclide_vec, & - atoms_vec, uncertainty_vec) - end if - - deallocate(nuclide_vec, atoms_vec, uncertainty_vec, volume) - end do - - ! Show elapsed time - if (master) then - call time_volume % stop() - call write_message("Elapsed time: " // trim(to_str(time_volume % & - get_value())) // " s", 6) - end if - err = 0 - end function openmc_calculate_volumes - -!=============================================================================== -! GET_VOLUME stochastically determines the volume of a set of domains along with -! the average number densities of nuclides within the domain -!=============================================================================== - - subroutine get_volume(this, volume, nuclide_vec, atoms_vec, uncertainty_vec) - type(VolumeCalculation), intent(in) :: this - real(8), intent(out) :: volume(:,:) ! volume mean/stdev in each domain - type(VectorInt), intent(out) :: nuclide_vec(:) ! indices in nuclides array - type(VectorReal), intent(out) :: atoms_vec(:) ! total # of atoms of each nuclide - type(VectorReal), intent(out) :: uncertainty_vec(:) ! uncertainty of total # of atoms - - ! Variables that are private to each thread - integer(8) :: i - integer :: j, k - integer :: i_domain ! index in domain_id array - integer :: i_material ! index in materials array - integer :: level ! local coordinate level - integer :: n_mat(size(this % domain_id)) ! Number of materials for each domain - integer, allocatable :: indices(:,:) ! List of material indices for each domain - integer, allocatable :: hits(:,:) ! Number of hits for each material in each domain - logical :: found_cell - type(Particle) :: p - - ! Shared variables - integer :: i_start, i_end ! Starting/ending sample for each process - type(VectorInt) :: master_indices(size(this % domain_id)) - type(VectorInt) :: master_hits(size(this % domain_id)) - - ! Variables used outside of parallel region - integer :: i_nuclide ! index in nuclides array - integer :: total_hits ! total hits for a single domain (summed over materials) - integer :: min_samples ! minimum number of samples per process - integer :: remainder ! leftover samples from uneven divide -#ifdef OPENMC_MPI - integer :: m ! index over materials - integer(C_INT) :: n ! number of materials - integer(C_INT), allocatable :: data(:) ! array used to send number of hits -#endif - real(8) :: f ! fraction of hits - real(8) :: var_f ! variance of fraction of hits - real(8) :: volume_sample ! total volume of sampled region - real(8) :: atoms(2, size(nuclides)) - -#ifdef OPENMC_MPI - interface - subroutine send_int(buffer, count, dest, tag) bind(C) - import C_INT - integer(C_INT), intent(in) :: buffer - integer(C_INT), value :: count - integer(C_INT), value :: dest - integer(C_INT), value :: tag - end subroutine - - subroutine recv_int(buffer, count, source, tag) bind(C) - import C_INT - integer(C_INT), intent(out) :: buffer - integer(C_INT), value :: count - integer(C_INT), value :: source - integer(C_INT), value :: tag - end subroutine - end interface -#endif - - ! Divide work over MPI processes - min_samples = this % samples / n_procs - remainder = mod(this % samples, n_procs) - if (rank < remainder) then - i_start = (min_samples + 1)*rank - i_end = i_start + min_samples - else - i_start = (min_samples + 1)*remainder + (rank - remainder)*min_samples - i_end = i_start + min_samples - 1 - end if - - call particle_initialize(p) - -!$omp parallel private(i, j, k, i_domain, i_material, level, found_cell, & -!$omp& indices, hits, n_mat) firstprivate(p) - - ! Create space for material indices and number of hits for each - allocate(indices(size(this % domain_id), 8)) - allocate(hits(size(this % domain_id), 8)) - n_mat(:) = 0 - - call prn_set_stream(STREAM_VOLUME) - - ! ========================================================================== - ! SAMPLES LOCATIONS AND COUNT HITS - -!$omp do - SAMPLE_LOOP: do i = i_start, i_end - call set_particle_seed(i) - - p % n_coord = 1 - p % coord(1) % xyz(1) = this % lower_left(1) + prn()*(& - this % upper_right(1) - this % lower_left(1)) - p % coord(1) % xyz(2) = this % lower_left(2) + prn()*(& - this % upper_right(2) - this % lower_left(2)) - p % coord(1) % xyz(3) = this % lower_left(3) + prn()*(& - this % upper_right(3) - this % lower_left(3)) - p % coord(1) % uvw(:) = [HALF, HALF, HALF] - - ! If this location is not in the geometry at all, move on to the next - ! block - call find_cell(p, found_cell) - if (.not. found_cell) cycle - - if (this % domain_type == FILTER_MATERIAL) then - i_material = p % material - if (i_material /= MATERIAL_VOID) then - do i_domain = 1, size(this % domain_id) - if (material_id(i_material) == this % domain_id(i_domain)) then - call check_hit(i_domain, i_material, indices, hits, n_mat) - end if - end do - end if - - elseif (this % domain_type == FILTER_CELL) THEN - do level = 1, p % n_coord - do i_domain = 1, size(this % domain_id) - if (cells(p % coord(level) % cell + 1) % id() & - == this % domain_id(i_domain)) then - i_material = p % material - call check_hit(i_domain, i_material, indices, hits, n_mat) - end if - end do - end do - - elseif (this % domain_type == FILTER_UNIVERSE) then - do level = 1, p % n_coord - do i_domain = 1, size(this % domain_id) - if (universe_id(p % coord(level) % universe) == & - this % domain_id(i_domain)) then - i_material = p % material - call check_hit(i_domain, i_material, indices, hits, n_mat) - end if - end do - end do - - end if - end do SAMPLE_LOOP -!$omp end do - - ! ========================================================================== - ! REDUCE HITS ONTO MASTER THREAD - - ! At this point, each thread has its own pair of index/hits lists and we now - ! need to reduce them. OpenMP is not nearly smart enough to do this on its - ! own, so we have to manually reduce them. - -#ifdef _OPENMP -!$omp do ordered schedule(static) - THREAD_LOOP: do i = 1, omp_get_num_threads() -!$omp ordered - do i_domain = 1, size(this % domain_id) - INDEX_LOOP: do j = 1, n_mat(i_domain) - ! Check if this material has been added to the master list and if so, - ! accumulate the number of hits - do k = 1, master_indices(i_domain) % size() - if (indices(i_domain, j) == master_indices(i_domain) % data(k)) then - master_hits(i_domain) % data(k) = & - master_hits(i_domain) % data(k) + hits(i_domain, j) - cycle INDEX_LOOP - end if - end do - - ! If we made it here, this means the material hasn't yet been added to - ! the master list, so add an entry to both the master indices and master - ! hits lists - call master_indices(i_domain) % push_back(indices(i_domain, j)) - call master_hits(i_domain) % push_back(hits(i_domain, j)) - end do INDEX_LOOP - end do -!$omp end ordered - end do THREAD_LOOP -!$omp end do -#else - do i_domain = 1, size(this % domain_id) - do j = 1, n_mat(i_domain) - call master_indices(i_domain) % push_back(indices(i_domain, j)) - call master_hits(i_domain) % push_back(hits(i_domain, j)) - end do - end do -#endif - - call prn_set_stream(STREAM_TRACKING) -!$omp end parallel - - ! ========================================================================== - ! REDUCE HITS ONTO MASTER PROCESS - - volume_sample = product(this % upper_right - this % lower_left) - - do i_domain = 1, size(this % domain_id) - atoms(:, :) = ZERO - total_hits = 0 - - if (master) then -#ifdef OPENMC_MPI - do j = 1, n_procs - 1 - call recv_int(n, 1, j, 0) - - allocate(data(2*n)) - call recv_int(data(1), 2*n, j, 1) - do k = 0, n - 1 - do m = 1, master_indices(i_domain) % size() - if (data(2*k + 1) == master_indices(i_domain) % data(m)) then - master_hits(i_domain) % data(m) = master_hits(i_domain) % data(m) + & - data(2*k + 2) - end if - end do - end do - deallocate(data) - end do -#endif - - do j = 1, master_indices(i_domain) % size() - total_hits = total_hits + master_hits(i_domain) % data(j) - f = real(master_hits(i_domain) % data(j), 8) / this % samples - var_f = f*(ONE - f) / this % samples - - i_material = master_indices(i_domain) % data(j) - if (i_material == MATERIAL_VOID) cycle - - do k = 1, material_nuclide_size(i_material) - ! Accumulate nuclide density - i_nuclide = material_nuclide(i_material, k) - atoms(1, i_nuclide) = atoms(1, i_nuclide) + & - material_atom_density(i_material, k) * f - atoms(2, i_nuclide) = atoms(2, i_nuclide) + & - material_atom_density(i_material, k)**2 * var_f - end do - end do - - ! Determine volume - volume(1, i_domain) = real(total_hits, 8) / this % samples * volume_sample - volume(2, i_domain) = sqrt(volume(1, i_domain) * (volume_sample - & - volume(1, i_domain)) / this % samples) - - ! Determine total number of atoms. At this point, we have values in - ! atoms/b-cm. To get to atoms we multiple by 10^24 V. - do j = 1, size(atoms, 2) - atoms(1, j) = 1.0e24_8 * volume_sample * atoms(1, j) - atoms(2, j) = 1.0e24_8 * volume_sample * sqrt(atoms(2, j)) - end do - - ! Convert full arrays to vectors - do j = 1, size(nuclides) - if (atoms(1, j) > ZERO) then - call nuclide_vec(i_domain) % push_back(j) - call atoms_vec(i_domain) % push_back(atoms(1, j)) - call uncertainty_vec(i_domain) % push_back(atoms(2, j)) - end if - end do - - else -#ifdef OPENMC_MPI - n = master_indices(i_domain) % size() - allocate(data(2*n)) - do k = 0, n - 1 - data(2*k + 1) = master_indices(i_domain) % data(k + 1) - data(2*k + 2) = master_hits(i_domain) % data(k + 1) - end do - - call send_int(n, 1, 0, 0) - call send_int(data(1), 2*n, 0, 1) - deallocate(data) -#endif - end if - end do - - contains - - !=========================================================================== - ! CHECK_HIT is an internal subroutine that checks for whether a material has - ! already been hit for a given domain. If not, it increases the list size by - ! one (taking care of re-allocation if needed). - !=========================================================================== - - subroutine check_hit(i_domain, i_material, indices, hits, n_mat) - integer :: i_domain - integer :: i_material - integer, allocatable :: indices(:,:) - integer, allocatable :: hits(:,:) - integer :: n_mat(:) - - integer, allocatable :: temp(:,:) - logical :: already_hit - integer :: j, k, nm - - ! Check if we've already had a hit in this material and if so, - ! simply add one - already_hit = .false. - nm = n_mat(i_domain) - do j = 1, nm - if (indices(i_domain, j) == i_material) then - hits(i_domain, j) = hits(i_domain, j) + 1 - already_hit = .true. - end if - end do - - if (.not. already_hit) then - ! If we make it here, that means we haven't yet had a hit in this - ! material. First check if the indices and hits arrays are large enough - ! and if not, double them. - if (nm == size(indices, 2)) then - k = 2*size(indices, 2) - allocate(temp(size(this % domain_id), k)) - temp(:, 1:nm) = indices(:, 1:nm) - call move_alloc(FROM=temp, TO=indices) - - allocate(temp(size(this % domain_id), k)) - temp(:, 1:nm) = hits(:, 1:nm) - call move_alloc(FROM=temp, TO=hits) - end if - - ! Add an entry to both the indices list and the hits list - n_mat(i_domain) = n_mat(i_domain) + 1 - indices(i_domain, n_mat(i_domain)) = i_material - hits(i_domain, n_mat(i_domain)) = 1 - end if - end subroutine check_hit - - end subroutine get_volume - -!=============================================================================== -! WRITE_VOLUME writes the results of a single stochastic volume calculation to -! an HDF5 file -!=============================================================================== - - subroutine write_volume(this, filename, volume, nuclide_vec, atoms_vec, & - uncertainty_vec) - type(VolumeCalculation), intent(in) :: this - character(*), intent(in) :: filename ! filename for HDF5 file - real(8), intent(in) :: volume(:,:) ! volume mean/stdev in each domain - type(VectorInt), intent(in) :: nuclide_vec(:) ! indices in nuclides array - type(VectorReal), intent(in) :: atoms_vec(:) ! total # of atoms of each nuclide - type(VectorReal), intent(in) :: uncertainty_vec(:) ! uncertainty of total # of atoms - - integer :: i, j - integer :: n - integer(HID_T) :: file_id - integer(HID_T) :: group_id - real(8), allocatable :: atom_data(:,:) ! mean/stdev of total # of atoms for - ! each nuclide - character(MAX_WORD_LEN), allocatable :: nucnames(:) ! names of nuclides - - ! Create HDF5 file - file_id = file_open(filename, 'w') - - ! Write header info - call write_attribute(file_id, "filetype", "volume") - call write_attribute(file_id, "version", VERSION_VOLUME) - call write_attribute(file_id, "openmc_version", VERSION) -#ifdef GIT_SHA1 - call write_attribute(file_id, "git_sha1", GIT_SHA1) -#endif - - ! Write current date and time - call write_attribute(file_id, "date_and_time", time_stamp()) - - ! Write basic metadata - select case (this % domain_type) - case (FILTER_CELL) - call write_attribute(file_id, "domain_type", "cell") - case (FILTER_MATERIAL) - call write_attribute(file_id, "domain_type", "material") - case (FILTER_UNIVERSE) - call write_attribute(file_id, "domain_type", "universe") - end select - call write_attribute(file_id, "samples", this % samples) - call write_attribute(file_id, "lower_left", this % lower_left) - call write_attribute(file_id, "upper_right", this % upper_right) - - do i = 1, size(this % domain_id) - group_id = create_group(file_id, "domain_" // trim(to_str(& - this % domain_id(i)))) - - ! Write volume for domain - call write_dataset(group_id, "volume", volume(:, i)) - - ! Create array of nuclide names from the vector - n = nuclide_vec(i) % size() - if (n > 0) then - allocate(nucnames(n)) - do j = 1, n - nucnames(j) = nuclides(nuclide_vec(i) % data(j)) % name - end do - - ! Create array of total # of atoms with uncertainty for each nuclide - allocate(atom_data(2, n)) - atom_data(1, :) = atoms_vec(i) % data(1:n) - atom_data(2, :) = uncertainty_vec(i) % data(1:n) - - ! Write results - call write_dataset(group_id, "nuclides", nucnames) - call write_dataset(group_id, "atoms", atom_data) - - deallocate(nucnames) - deallocate(atom_data) - end if - - call close_group(group_id) - end do - call file_close(file_id) - end subroutine write_volume - -end module volume_calc diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp new file mode 100644 index 0000000000..2599a80966 --- /dev/null +++ b/src/volume_calc.cpp @@ -0,0 +1,425 @@ +#include "openmc/volume_calc.h" + +#include "openmc/capi.h" +#include "openmc/cell.h" +#include "openmc/constants.h" +#include "openmc/error.h" +#include "openmc/geometry.h" +#include "openmc/hdf5_interface.h" +#include "openmc/material.h" +#include "openmc/message_passing.h" +#include "openmc/nuclide.h" +#include "openmc/output.h" +#include "openmc/random_lcg.h" +#include "openmc/settings.h" +#include "openmc/timer.h" +#include "openmc/xml_interface.h" + +#ifdef _OPENMP +#include +#endif +#include "xtensor/xadapt.hpp" +#include "xtensor/xview.hpp" + +#include // for copy +#include // for pow, sqrt +#include + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +namespace model { +std::vector volume_calcs; +} + +//============================================================================== +// VolumeCalculation implementation +//============================================================================== + +VolumeCalculation::VolumeCalculation(pugi::xml_node node) { + // Read domain type (cell, material or universe) + std::string domain_type = get_node_value(node, "domain_type"); + if (domain_type == "cell") { + domain_type_ = FILTER_CELL; + } + else if (domain_type == "material") { + domain_type_ = FILTER_MATERIAL; + } + else if (domain_type == "universe") { + domain_type_ = FILTER_UNIVERSE; + } + else { + fatal_error(std::string("Unrecognized domain type for stochastic " + "volume calculation:" + domain_type)); + } + + // Read domain IDs, bounding corodinates and number of samples + domain_ids_ = get_node_array(node, "domain_ids"); + lower_left_ = get_node_array(node, "lower_left"); + upper_right_ = get_node_array(node, "upper_right"); + n_samples_ = std::stoi(get_node_value(node, "samples")); +} + +void VolumeCalculation::check_hit(int i_domain, int i_material, + int_2dvec& indices, int_2dvec& hits) const { + + // Check if this material was previously hit and if so, increment count + bool already_hit = false; + for (int j = 0; j < indices[i_domain].size(); j++) { + if (indices[i_domain][j] == i_material) { + hits[i_domain][j]++; + already_hit = true; + } + } + + // If the material was not previously hit, append an entry to the material + // indices and hits lists + if (!already_hit) { + indices[i_domain].push_back(i_material); + hits[i_domain].push_back(1); + } +} + +// Stochastically determin the volume of a set of domains along with the +// average number densities of nuclides within the domain +// @param volumes volume mean/stdev in each domain +// @param i_nuclides indices in nuclides array +// @param n_atoms total # of atoms of each nuclide +// @param n_atoms_uncertainty uncertainty of total # of atoms +std::vector VolumeCalculation::execute() const +{ + // Shared data that is collected from all threads + int n = domain_ids_.size(); + int_2dvec master_indices(n); // List of material indices for each domain + int_2dvec master_hits(n); // Number of hits for each material in each domain + + // Divide work over MPI processes + int min_samples = n_samples_ / mpi::n_procs; + int remainder = n_samples_ % mpi::n_procs; + int i_start, i_end; + if (mpi::rank < remainder) { + i_start = (min_samples + 1)*mpi::rank; + i_end = i_start + min_samples + 1; + } + else { + i_start = (min_samples + 1)*remainder + (mpi::rank - remainder)*min_samples; + i_end = i_start + min_samples; + } + +#pragma omp parallel + { + // Variables that are private to each thread + int_2dvec indices(n); + int_2dvec hits(n); + Particle p; + p.initialize(); + + prn_set_stream(STREAM_VOLUME); + + // Samples locations and count hits + #pragma omp for + for (int i = i_start; i < i_end; i++) { + set_particle_seed(i); + + p.n_coord = 1; + Position xi {prn(), prn(), prn()}; + Position r {lower_left_ + xi*(upper_right_ - lower_left_)}; + // TODO: assign directly when xyz is Position + std::copy(&r.x, &r.x + 3, p.coord[0].xyz); + p.coord[0].uvw[0] = 0.5; + p.coord[1].uvw[1] = 0.5; + p.coord[2].uvw[2] = 0.5; + + // If this location is not in the geometry at all, move on to next block + if (!find_cell(&p, false)) continue; + + // TODO: off-by-one + int i_material = p.material == MATERIAL_VOID ? p.material : p.material - 1; + + if (domain_type_ == FILTER_MATERIAL) { + if (i_material != MATERIAL_VOID) { + for (int i_domain = 0; i_domain < n; i_domain++) { + if (model::materials[i_material]->id_ == domain_ids_[i_domain]) { + this->check_hit(i_domain, i_material, indices, hits); + } + } + } + } else if (domain_type_ == FILTER_CELL) { + for (int level = 0; level < p.n_coord; ++level) { + for (int i_domain=0; i_domain < n; i_domain++) { + if (model::cells[p.coord[level].cell]->id_ == domain_ids_[i_domain]) { + this->check_hit(i_domain, i_material, indices, hits); + } + } + } + } else if (domain_type_ == FILTER_UNIVERSE) { + for (int level = 0; level < p.n_coord; ++level) { + for (int i_domain = 0; i_domain < n; ++i_domain) { + if (model::universes[p.coord[level].universe]->id_ == domain_ids_[i_domain]) { + check_hit(i_domain, i_material, indices, hits); + } + } + } + } + } + + // Reduce hits onto master thread + // At this point, each thread has its own pair of index/hits lists and we now + // need to reduce them. OpenMP is not nearly smart enough to do this on its own, + // so we have to manually reduce them + +#ifdef _OPENMP +#pragma omp for ordered schedule(static) + for (int i = 0; i < omp_get_num_threads(); ++i) { +#pragma omp ordered + for (int i_domain = 0; i_domain < n; ++i_domain) { + for (int j = 0; j < indices[i_domain].size(); ++j) { + // Check if this material has been added to the master list and if so, + // accumulate the number of hits + bool already_added = false; + for (int k = 0; k < master_indices[i_domain].size(); k++) { + if (indices[i_domain][j] == master_indices[i_domain][k]) { + master_hits[i_domain][k] += hits[i_domain][j]; + already_added = true; + } + } + if (!already_added) { + // If we made it here, the material hasn't yet been added to the master + // list, so add entries to the master indices and master hits lists + master_indices[i_domain].push_back(indices[i_domain][j]); + master_hits[i_domain].push_back(hits[i_domain][j]); + } + } + } + } +#else + for (int i_domain = 0; i_domain < n; ++i_domain) { + for (int j = 0; j < indices[i_domain].size(); ++j) { + master_indices[i_domain].push_back(indices[i_domain][j]); + master_hits[i_domain].push_back(hits[i_domain][j]); + } + } +#endif + + prn_set_stream(STREAM_TRACKING); + } // omp parallel + + // Reduce hits onto master process + + // Determine volume of bounding box + Position d {upper_right_ - lower_left_}; + double volume_sample = d.x*d.y*d.z; + + // Set size for members of the Result struct + std::vector results(n); + + for (int i_domain = 0; i_domain < n; ++i_domain) { + // Get reference to result for this domain + auto& result {results[i_domain]}; + + // Create 2D array to store atoms/uncertainty for each nuclide. Later this + // is compressed into vectors storing only those nuclides that are non-zero + auto n_nuc = data::nuclides.size(); + xt::xtensor atoms({n_nuc, 2}, 0.0); + + if (mpi::master) { +#ifdef OPENMC_MPI + for (int j = 1; j < mpi::n_procs; j++) { + int q; + MPI_Recv(&q, 1, MPI_INTEGER, j, 0, mpi::intracomm, MPI_STATUS_IGNORE); + int buffer[2*q]; + MPI_Recv(&buffer[0], 2*q, MPI_INTEGER, j, 1, mpi::intracomm, MPI_STATUS_IGNORE); + for (int k = 0; k < q; ++k) { + for (int m = 0; m < master_indices[i_domain].size(); ++m) { + if (buffer[2*k] == master_indices[i_domain][m]) { + master_hits[i_domain][m] += buffer[2*k + 1]; + } + } + } + } +#endif + + int total_hits = 0; + for (int j = 0; j < master_indices[i_domain].size(); ++j) { + total_hits += master_hits[i_domain][j]; + double f = static_cast(master_hits[i_domain][j]) / n_samples_; + double var_f = f*(1.0 - f)/n_samples_; + + int i_material = master_indices[i_domain][j]; + if (i_material == MATERIAL_VOID) continue; + + const auto& mat = model::materials[i_material]; + for (int k = 0; k < mat->nuclide_.size(); ++k) { + // Accumulate nuclide density + int i_nuclide = mat->nuclide_[k]; + atoms(i_nuclide, 0) += mat->atom_density_[k] * f; + atoms(i_nuclide, 1) += std::pow(mat->atom_density_[k], 2) * var_f; + } + } + + // Determine volume + result.volume[0] = static_cast(total_hits) / n_samples_ * volume_sample; + result.volume[1] = std::sqrt(result.volume[0] + * (volume_sample - result.volume[0]) / n_samples_); + + for (int j = 0; j < n_nuc; ++j) { + // Determine total number of atoms. At this point, we have values in + // atoms/b-cm. To get to atoms we multiply by 10^24 V. + double mean = 1.0e24 * volume_sample * atoms(j, 0); + double stdev = 1.0e24 * volume_sample * std::sqrt(atoms(j, 1)); + + // Convert full arrays to vectors + if (mean > 0.0) { + result.nuclides.push_back(j); + result.atoms.push_back(mean); + result.uncertainty.push_back(stdev); + } + } + } else { +#ifdef OPENMC_MPI + int q = master_indices[i_domain].size(); + int buffer[2*q]; + for (int k = 0; k < q; ++k) { + buffer[2*k] = master_indices[i_domain][k]; + buffer[2*k + 1] = master_hits[i_domain][k]; + } + + MPI_Send(&q, 1, MPI_INTEGER, 0, 0, mpi::intracomm); + MPI_Send(&buffer[0], 2*q, MPI_INTEGER, 0, 1, mpi::intracomm); +#endif + } + } + + return results; +} + +void VolumeCalculation::write_volume(const std::string& filename, + const std::vector& results) const +{ + // Create HDF5 file + hid_t file_id = file_open(filename, 'w'); + + // Write header info + write_attribute(file_id, "filetype", "volume"); + write_attribute(file_id, "version", VERSION_VOLUME); + write_attribute(file_id, "openmc_version", VERSION); +#ifdef GIT_SHA1 + write_attribute(file_id, "git_sha1", GIT_SHA1); +#endif + + // Write current date and time + write_attribute(file_id, "date_and_time", time_stamp()); + + // Write basic metadata + write_attribute(file_id, "samples", n_samples_); + write_attribute(file_id, "lower_left", lower_left_); + write_attribute(file_id, "upper_right", upper_right_); + if (domain_type_ == FILTER_CELL) { + write_attribute(file_id, "domain_type", "cell"); + } + else if (domain_type_ == FILTER_MATERIAL) { + write_attribute(file_id, "domain_type", "material"); + } + else if (domain_type_ == FILTER_UNIVERSE) { + write_attribute(file_id, "domain_type", "universe"); + } + + for (int i = 0; i < domain_ids_.size(); ++i) + { + hid_t group_id = create_group(file_id, "domain_" + + std::to_string(domain_ids_[i])); + + // Write volume for domain + const auto& result {results[i]}; + write_dataset(group_id, "volume", result.volume); + + // Create array of nuclide names from the vector + auto n_nuc = result.nuclides.size(); + + if (!result.nuclides.empty()) { + std::vector nucnames; + for (int i_nuc : result.nuclides) { + nucnames.push_back(data::nuclides[i_nuc]->name_); + } + + // Create array of total # of atoms with uncertainty for each nuclide + xt::xtensor atom_data({n_nuc, 2}); + xt::view(atom_data, xt::all(), 0) = xt::adapt(result.atoms); + xt::view(atom_data, xt::all(), 1) = xt::adapt(result.uncertainty); + + // Write results + write_dataset(group_id, "nuclides", nucnames); + write_dataset(group_id, "atoms", atom_data); + } + + close_group(group_id); + } + + file_close(file_id); +} + +} // namespace openmc + +//============================================================================== +// OPENMC_CALCULATE_VOLUMES runs each of the stochastic volume calculations +// that the user has specified and writes results to HDF5 files +//============================================================================== + +int openmc_calculate_volumes() { + using namespace openmc; + + if (mpi::master) { + header("STOCHASTIC VOLUME CALCULATION", 3); + } + Timer time_volume; + time_volume.start(); + + for (int i = 0; i < model::volume_calcs.size(); ++i) { + if (mpi::master) { + write_message("Running volume calculation " + std::to_string(i+1) + "...", 4); + } + + // Run volume calculation + const auto& vol_calc {model::volume_calcs[i]}; + auto results = vol_calc.execute(); + + if (mpi::master) { + std::string domain_type; + if (vol_calc.domain_type_ == FILTER_CELL) { + domain_type = " Cell "; + } else if (vol_calc.domain_type_ == FILTER_MATERIAL) { + domain_type = " Material "; + } else { + domain_type = " Universe "; + } + + // Display domain volumes + for (int j = 0; j < vol_calc.domain_ids_.size(); j++) { + std::stringstream msg; + msg << domain_type << vol_calc.domain_ids_[j] << ": " << + results[j].volume[0] << " +/- " << results[j].volume[1] << " cm^3"; + write_message(msg, 4); + } + + // Write volumes to HDF5 file + std::string filename = settings::path_output + "volume_" + + std::to_string(i+1) + ".h5"; + vol_calc.write_volume(filename, results); + } + + } + + // Show elapsed time + time_volume.stop(); + if (mpi::master) { + write_message("Elapsed time: " + std::to_string(time_volume.elapsed()) + + " s", 6); + } + + return 0; +} + +extern "C" void free_memory_volume() { openmc::model::volume_calcs.clear(); } diff --git a/src/volume_header.F90 b/src/volume_header.F90 deleted file mode 100644 index 59a579a5fb..0000000000 --- a/src/volume_header.F90 +++ /dev/null @@ -1,69 +0,0 @@ -module volume_header - - use constants, only: FILTER_CELL, FILTER_MATERIAL, FILTER_UNIVERSE - use error, only: fatal_error - use xml_interface - - implicit none - - type VolumeCalculation - integer :: domain_type - integer, allocatable :: domain_id(:) - real(8) :: lower_left(3) - real(8) :: upper_right(3) - integer :: samples - contains - procedure :: from_xml => volume_from_xml - end type VolumeCalculation - - type(VolumeCalculation), allocatable :: volume_calcs(:) - -contains - - subroutine volume_from_xml(this, node_vol) - class(VolumeCalculation), intent(out) :: this - type(XMLNode), intent(in) :: node_vol - - integer :: num_domains - character(10) :: temp_str - - ! Check domain type - call get_node_value(node_vol, "domain_type", temp_str) - select case (temp_str) - case ('cell') - this % domain_type = FILTER_CELL - case ('material') - this % domain_type = FILTER_MATERIAL - case ('universe') - this % domain_type = FILTER_UNIVERSE - case default - call fatal_error("Unrecognized domain type for stochastic volume & - &calculation: " // trim(temp_str)) - end select - - ! Read cell IDs - if (check_for_node(node_vol, "domain_ids")) then - num_domains = node_word_count(node_vol, "domain_ids") - else - call fatal_error("Must specify at least one cell for a volume calculation") - end if - allocate(this % domain_id(num_domains)) - call get_node_array(node_vol, "domain_ids", this % domain_id) - - ! Read lower-left and upper-right bounding coordinates - call get_node_array(node_vol, "lower_left", this % lower_left) - call get_node_array(node_vol, "upper_right", this % upper_right) - - ! Read number of samples - call get_node_value(node_vol, "samples", this % samples) - end subroutine volume_from_xml - -!=============================================================================== -! FREE_MEMORY_VOLUME deallocates global arrays defined in this module -!=============================================================================== - - subroutine free_memory_volume() - if (allocated(volume_calcs)) deallocate(volume_calcs) - end subroutine free_memory_volume - -end module volume_header From 9bce2c357b4e090799b10ac26ccfdaf1caaa82f7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 6 Feb 2019 22:18:00 -0600 Subject: [PATCH 2/7] Remove unused Fortran code --- CMakeLists.txt | 1 - src/cell.cpp | 5 --- src/geometry.F90 | 2 - src/geometry_header.F90 | 6 --- src/output.F90 | 19 +-------- src/reaction_header.F90 | 3 -- src/string.F90 | 3 +- src/timer_header.F90 | 88 ----------------------------------------- 8 files changed, 2 insertions(+), 125 deletions(-) delete mode 100644 src/timer_header.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 48e9f622b5..d62c07308b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -338,7 +338,6 @@ add_library(libopenmc SHARED src/stl_vector.F90 src/string.F90 src/surface_header.F90 - src/timer_header.F90 src/tracking.F90 src/track_output.F90 src/vector_header.F90 diff --git a/src/cell.cpp b/src/cell.cpp index 1dec11f5dd..6420595067 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -845,11 +845,6 @@ extern "C" { } model::n_cells = model::cells.size(); } - - int32_t universe_id(int i_univ) {return model::universes[i_univ]->id_;} - - void universes_to_hdf5(hid_t universes_group) - {for (Universe* u : model::universes) u->to_hdf5(universes_group);} } diff --git a/src/geometry.F90 b/src/geometry.F90 index 83bf2325f8..7456644e67 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -1,13 +1,11 @@ module geometry use constants - use error, only: fatal_error, warning, write_message use geometry_header use particle_header use simulation_header use settings use surface_header - use stl_vector, only: VectorInt use string, only: to_str use, intrinsic :: ISO_C_BINDING diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 8c18a40308..b6c454b9dc 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -8,12 +8,6 @@ module geometry_header implicit none interface - function universe_id(universe_ind) bind(C) result(id) - import C_INT, C_INT32_T - integer(C_INT), intent(in), value :: universe_ind - integer(C_INT32_T) :: id - end function universe_id - function cell_pointer(cell_ind) bind(C) result(ptr) import C_PTR, C_INT32_T integer(C_INT32_T), intent(in), value :: cell_ind diff --git a/src/output.F90 b/src/output.F90 index bfab4ddbff..bb04bdae3f 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -4,7 +4,6 @@ module output use, intrinsic :: ISO_FORTRAN_ENV use constants - use eigenvalue, only: openmc_get_keff use endf, only: reaction_name use error, only: fatal_error, warning use geometry_header @@ -12,7 +11,7 @@ module output use message_passing, only: master, n_procs use mgxs_interface use nuclide_header - use particle_header, only: LocalCoord, Particle + use particle_header, only: Particle use settings use simulation_header use surface_header, only: surfaces @@ -38,22 +37,6 @@ module output contains -!=============================================================================== -! TIME_STAMP returns the current date and time in a formatted string -!=============================================================================== - - function time_stamp() result(current_time) - - character(19) :: current_time ! ccyy-mm-dd hh:mm:ss - character(8) :: date_ ! ccyymmdd - character(10) :: time_ ! hhmmss.sss - - call date_and_time(DATE=date_, TIME=time_) - current_time = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) // & - " " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) - - end function time_stamp - !=============================================================================== ! HEADER displays a header block according to a specified level. If no level is ! specified, it is assumed to be a minor header block. diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index 9d57db5ccf..bc0fc8f326 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -2,10 +2,7 @@ module reaction_header use, intrinsic :: ISO_C_BINDING - use constants, only: MAX_WORD_LEN use hdf5_interface - use stl_vector, only: VectorInt - use string, only: to_str, starts_with implicit none private diff --git a/src/string.F90 b/src/string.F90 index f0419e2c14..20fd36c469 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -3,8 +3,7 @@ module string use, intrinsic :: ISO_C_BINDING use constants, only: ERROR_INT, ERROR_REAL - use error, only: fatal_error, warning - use stl_vector, only: VectorInt + use error, only: fatal_error implicit none diff --git a/src/timer_header.F90 b/src/timer_header.F90 deleted file mode 100644 index 3bbc207620..0000000000 --- a/src/timer_header.F90 +++ /dev/null @@ -1,88 +0,0 @@ -module timer_header - - use, intrinsic :: ISO_C_BINDING - - use constants, only: ZERO - -!=============================================================================== -! TIMER represents a timer that can be started and stopped to measure how long -! different routines run. The intrinsic routine system_clock is used to measure -! time rather than cpu_time. -!=============================================================================== - - type Timer - private - logical :: running = .false. ! is timer running? - integer(8) :: start_counts = 0 ! counts when started - real(8), public :: elapsed = ZERO ! total time elapsed in seconds - contains - procedure :: start => timer_start - procedure :: get_value => timer_get_value - procedure :: stop => timer_stop - procedure :: reset => timer_reset - end type Timer - -contains - -!=============================================================================== -! TIMER_START starts running a timer and measures the current time -!=============================================================================== - - subroutine timer_start(self) - class(Timer), intent(inout) :: self - - ! Turn timer on and measure starting time - self % running = .true. - call system_clock(self % start_counts) - end subroutine timer_start - -!=============================================================================== -! TIMER_GET_VALUE returns the current value of the timer -!=============================================================================== - - function timer_get_value(self) result(elapsed) - class(Timer), intent(in) :: self ! the timer - real(8) :: elapsed ! total elapsed time - - integer(8) :: end_counts ! current number of counts - integer(8) :: count_rate ! system-dependent counting rate - real(8) :: elapsed_time ! elapsed time since last start - - if (self % running) then - call system_clock(end_counts, count_rate) - elapsed_time = real(end_counts - self % start_counts, 8) / & - real(count_rate, 8) - elapsed = self % elapsed + elapsed_time - else - elapsed = self % elapsed - end if - end function timer_get_value - -!=============================================================================== -! TIMER_STOP stops the timer and sets the elapsed time -!=============================================================================== - - subroutine timer_stop(self) - class(Timer), intent(inout) :: self - - ! Check to make sure timer was running - if (.not. self % running) return - - ! Stop timer and add time - self % elapsed = self % get_value() - self % running = .false. - end subroutine timer_stop - -!=============================================================================== -! TIMER_RESET resets a timer to have a zero value -!=============================================================================== - - pure subroutine timer_reset(self) - class(Timer), intent(inout) :: self - - self % running = .false. - self % start_counts = 0 - self % elapsed = ZERO - end subroutine timer_reset - -end module timer_header From bd2d96568bc30e5cfb9cf74d5a816f41fcc19b91 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 6 Feb 2019 23:46:45 -0600 Subject: [PATCH 3/7] A little more cleanup of volume calculations --- include/openmc/volume_calc.h | 25 +++++++++++++++------ src/volume_calc.cpp | 42 +++++++++++++++++------------------- 2 files changed, 39 insertions(+), 28 deletions(-) diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 5f905b8fa4..6bc52be5ee 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -18,8 +18,6 @@ namespace openmc { class VolumeCalculation { public: // Aliases, types - using int_2dvec = std::vector>; - struct Result { std::array volume; //!< Mean/standard deviation of volume std::vector nuclides; //!< Index of nuclides @@ -28,12 +26,21 @@ public: }; // Results for a single domain // Constructors - VolumeCalculation() = default; VolumeCalculation(pugi::xml_node node); // Methods + + //! \brief Stochastically determine the volume of a set of domains along with the + //! average number densities of nuclides within the domain + // + //! \return Vector of results for each user-specified domain std::vector execute() const; - void write_volume(const std::string& filename, const std::vector& results) const; + + //! \brief Write volume calculation results to HDF5 file + // + //! \param[in] filename Path to HDF5 file to write + //! \param[in] results Vector of results for each domain + void to_hdf5(const std::string& filename, const std::vector& results) const; // Data members int domain_type_; //!< Type of domain (cell, material, etc.) @@ -43,8 +50,14 @@ public: std::vector domain_ids_; //!< IDs of domains to find volumes of private: - void check_hit(int i_domain, int i_material, int_2dvec& indices, - int_2dvec& hits) const; + //! \brief Check whether a material has already been hit for a given domain. + //! If not, add new entries to the vectors + // + //! \param[in] i_material Index in global materials vector + //! \param[in,out] indices Vector of material indices + //! \param[in,out] hits Number of hits corresponding to each material + void check_hit(int i_material, std::vector& indices, + std::vector& hits) const; }; diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 2599a80966..2d6c10eec2 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -63,14 +63,14 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) { n_samples_ = std::stoi(get_node_value(node, "samples")); } -void VolumeCalculation::check_hit(int i_domain, int i_material, - int_2dvec& indices, int_2dvec& hits) const { +void VolumeCalculation::check_hit(int i_material, std::vector& indices, + std::vector& hits) const { // Check if this material was previously hit and if so, increment count bool already_hit = false; - for (int j = 0; j < indices[i_domain].size(); j++) { - if (indices[i_domain][j] == i_material) { - hits[i_domain][j]++; + for (int j = 0; j < indices.size(); j++) { + if (indices[j] == i_material) { + hits[j]++; already_hit = true; } } @@ -78,23 +78,17 @@ void VolumeCalculation::check_hit(int i_domain, int i_material, // If the material was not previously hit, append an entry to the material // indices and hits lists if (!already_hit) { - indices[i_domain].push_back(i_material); - hits[i_domain].push_back(1); + indices.push_back(i_material); + hits.push_back(1); } } -// Stochastically determin the volume of a set of domains along with the -// average number densities of nuclides within the domain -// @param volumes volume mean/stdev in each domain -// @param i_nuclides indices in nuclides array -// @param n_atoms total # of atoms of each nuclide -// @param n_atoms_uncertainty uncertainty of total # of atoms std::vector VolumeCalculation::execute() const { // Shared data that is collected from all threads int n = domain_ids_.size(); - int_2dvec master_indices(n); // List of material indices for each domain - int_2dvec master_hits(n); // Number of hits for each material in each domain + std::vector> master_indices(n); // List of material indices for each domain + std::vector> master_hits(n); // Number of hits for each material in each domain // Divide work over MPI processes int min_samples = n_samples_ / mpi::n_procs; @@ -112,8 +106,8 @@ std::vector VolumeCalculation::execute() const #pragma omp parallel { // Variables that are private to each thread - int_2dvec indices(n); - int_2dvec hits(n); + std::vector> indices(n); + std::vector> hits(n); Particle p; p.initialize(); @@ -143,7 +137,7 @@ std::vector VolumeCalculation::execute() const if (i_material != MATERIAL_VOID) { for (int i_domain = 0; i_domain < n; i_domain++) { if (model::materials[i_material]->id_ == domain_ids_[i_domain]) { - this->check_hit(i_domain, i_material, indices, hits); + this->check_hit(i_material, indices[i_domain], hits[i_domain]); } } } @@ -151,7 +145,7 @@ std::vector VolumeCalculation::execute() const for (int level = 0; level < p.n_coord; ++level) { for (int i_domain=0; i_domain < n; i_domain++) { if (model::cells[p.coord[level].cell]->id_ == domain_ids_[i_domain]) { - this->check_hit(i_domain, i_material, indices, hits); + this->check_hit(i_material, indices[i_domain], hits[i_domain]); } } } @@ -159,7 +153,7 @@ std::vector VolumeCalculation::execute() const for (int level = 0; level < p.n_coord; ++level) { for (int i_domain = 0; i_domain < n; ++i_domain) { if (model::universes[p.coord[level].universe]->id_ == domain_ids_[i_domain]) { - check_hit(i_domain, i_material, indices, hits); + check_hit(i_material, indices[i_domain], hits[i_domain]); } } } @@ -296,7 +290,7 @@ std::vector VolumeCalculation::execute() const return results; } -void VolumeCalculation::write_volume(const std::string& filename, +void VolumeCalculation::to_hdf5(const std::string& filename, const std::vector& results) const { // Create HDF5 file @@ -407,7 +401,7 @@ int openmc_calculate_volumes() { // Write volumes to HDF5 file std::string filename = settings::path_output + "volume_" + std::to_string(i+1) + ".h5"; - vol_calc.write_volume(filename, results); + vol_calc.to_hdf5(filename, results); } } @@ -422,4 +416,8 @@ int openmc_calculate_volumes() { return 0; } +//============================================================================== +// Fortran compatibility +//============================================================================== + extern "C" void free_memory_volume() { openmc::model::volume_calcs.clear(); } From 57bd735fe812f552e259357ba78084a2b6d97fac Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Feb 2019 13:39:13 -0600 Subject: [PATCH 4/7] Address @smharper comments on #1164 --- src/volume_calc.cpp | 94 +++++++++++++++++++++------------------------ 1 file changed, 43 insertions(+), 51 deletions(-) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 2d6c10eec2..d13057799a 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -44,16 +44,13 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) { std::string domain_type = get_node_value(node, "domain_type"); if (domain_type == "cell") { domain_type_ = FILTER_CELL; - } - else if (domain_type == "material") { + } else if (domain_type == "material") { domain_type_ = FILTER_MATERIAL; - } - else if (domain_type == "universe") { + } else if (domain_type == "universe") { domain_type_ = FILTER_UNIVERSE; - } - else { + } else { fatal_error(std::string("Unrecognized domain type for stochastic " - "volume calculation:" + domain_type)); + "volume calculation: " + domain_type)); } // Read domain IDs, bounding corodinates and number of samples @@ -63,26 +60,6 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) { n_samples_ = std::stoi(get_node_value(node, "samples")); } -void VolumeCalculation::check_hit(int i_material, std::vector& indices, - std::vector& hits) const { - - // Check if this material was previously hit and if so, increment count - bool already_hit = false; - for (int j = 0; j < indices.size(); j++) { - if (indices[j] == i_material) { - hits[j]++; - already_hit = true; - } - } - - // If the material was not previously hit, append an entry to the material - // indices and hits lists - if (!already_hit) { - indices.push_back(i_material); - hits.push_back(1); - } -} - std::vector VolumeCalculation::execute() const { // Shared data that is collected from all threads @@ -97,8 +74,7 @@ std::vector VolumeCalculation::execute() const if (mpi::rank < remainder) { i_start = (min_samples + 1)*mpi::rank; i_end = i_start + min_samples + 1; - } - else { + } else { i_start = (min_samples + 1)*remainder + (mpi::rank - remainder)*min_samples; i_end = i_start + min_samples; } @@ -113,8 +89,8 @@ std::vector VolumeCalculation::execute() const prn_set_stream(STREAM_VOLUME); - // Samples locations and count hits - #pragma omp for + // Sample locations and count hits +#pragma omp for for (int i = i_start; i < i_end; i++) { set_particle_seed(i); @@ -160,7 +136,6 @@ std::vector VolumeCalculation::execute() const } } - // Reduce hits onto master thread // At this point, each thread has its own pair of index/hits lists and we now // need to reduce them. OpenMP is not nearly smart enough to do this on its own, // so we have to manually reduce them @@ -190,12 +165,7 @@ std::vector VolumeCalculation::execute() const } } #else - for (int i_domain = 0; i_domain < n; ++i_domain) { - for (int j = 0; j < indices[i_domain].size(); ++j) { - master_indices[i_domain].push_back(indices[i_domain][j]); - master_hits[i_domain].push_back(hits[i_domain][j]); - } - } + master_indices = indices; #endif prn_set_stream(STREAM_TRACKING); @@ -219,8 +189,8 @@ std::vector VolumeCalculation::execute() const auto n_nuc = data::nuclides.size(); xt::xtensor atoms({n_nuc, 2}, 0.0); - if (mpi::master) { #ifdef OPENMC_MPI + if (mpi::master) { for (int j = 1; j < mpi::n_procs; j++) { int q; MPI_Recv(&q, 1, MPI_INTEGER, j, 0, mpi::intracomm, MPI_STATUS_IGNORE); @@ -230,12 +200,25 @@ std::vector VolumeCalculation::execute() const for (int m = 0; m < master_indices[i_domain].size(); ++m) { if (buffer[2*k] == master_indices[i_domain][m]) { master_hits[i_domain][m] += buffer[2*k + 1]; + break; } } } } + } else { + int q = master_indices[i_domain].size(); + int buffer[2*q]; + for (int k = 0; k < q; ++k) { + buffer[2*k] = master_indices[i_domain][k]; + buffer[2*k + 1] = master_hits[i_domain][k]; + } + + MPI_Send(&q, 1, MPI_INTEGER, 0, 0, mpi::intracomm); + MPI_Send(&buffer[0], 2*q, MPI_INTEGER, 0, 1, mpi::intracomm); + } #endif + if (mpi::master) { int total_hits = 0; for (int j = 0; j < master_indices[i_domain].size(); ++j) { total_hits += master_hits[i_domain][j]; @@ -272,18 +255,6 @@ std::vector VolumeCalculation::execute() const result.uncertainty.push_back(stdev); } } - } else { -#ifdef OPENMC_MPI - int q = master_indices[i_domain].size(); - int buffer[2*q]; - for (int k = 0; k < q; ++k) { - buffer[2*k] = master_indices[i_domain][k]; - buffer[2*k + 1] = master_hits[i_domain][k]; - } - - MPI_Send(&q, 1, MPI_INTEGER, 0, 0, mpi::intracomm); - MPI_Send(&buffer[0], 2*q, MPI_INTEGER, 0, 1, mpi::intracomm); -#endif } } @@ -355,6 +326,27 @@ void VolumeCalculation::to_hdf5(const std::string& filename, file_close(file_id); } +void VolumeCalculation::check_hit(int i_material, std::vector& indices, + std::vector& hits) const +{ + + // Check if this material was previously hit and if so, increment count + bool already_hit = false; + for (int j = 0; j < indices.size(); j++) { + if (indices[j] == i_material) { + hits[j]++; + already_hit = true; + } + } + + // If the material was not previously hit, append an entry to the material + // indices and hits lists + if (!already_hit) { + indices.push_back(i_material); + hits.push_back(1); + } +} + } // namespace openmc //============================================================================== From ca7a6adcd132640c6acb28364e1bba1def16f887 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Feb 2019 13:59:30 -0600 Subject: [PATCH 5/7] Get rid of output.F90 --- CMakeLists.txt | 1 - src/output.F90 | 79 ---------------------------------------- src/particle_restart.F90 | 6 ++- 3 files changed, 5 insertions(+), 81 deletions(-) delete mode 100644 src/output.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index d956b8d7e1..5a204fb66a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -323,7 +323,6 @@ add_library(libopenmc SHARED src/message_passing.F90 src/mgxs_interface.F90 src/nuclide_header.F90 - src/output.F90 src/particle_header.F90 src/particle_restart.F90 src/photon_header.F90 diff --git a/src/output.F90 b/src/output.F90 deleted file mode 100644 index a9430deb34..0000000000 --- a/src/output.F90 +++ /dev/null @@ -1,79 +0,0 @@ -module output - - use, intrinsic :: ISO_C_BINDING - use, intrinsic :: ISO_FORTRAN_ENV - - use constants - use endf, only: reaction_name - use error, only: fatal_error, warning - use geometry_header - use math, only: t_percentile - use message_passing, only: master, n_procs - use mgxs_interface - use nuclide_header - use particle_header, only: Particle - use settings - use simulation_header - use surface_header, only: surfaces - use string, only: to_upper, to_str - use tally_header - use tally_derivative_header - use tally_filter - - implicit none - - ! Short names for output and error units - integer :: ou = OUTPUT_UNIT - integer :: eu = ERROR_UNIT - - interface - subroutine print_particle(p) bind(C) - import Particle - type(Particle), intent(in) :: p - end subroutine - - subroutine write_tallies() bind(C) - end subroutine - end interface - -contains - -!=============================================================================== -! HEADER displays a header block according to a specified level. If no level is -! specified, it is assumed to be a minor header block. -!=============================================================================== - - subroutine header(msg, level, unit) - character(*), intent(in) :: msg ! header message - integer, intent(in) :: level - integer, intent(in), optional :: unit ! unit to write to - - integer :: n ! number of = signs on left - integer :: m ! number of = signs on right - integer :: unit_ ! unit to write to - character(MAX_LINE_LEN) :: line - - ! set default unit - if (present(unit)) then - unit_ = unit - else - unit_ = OUTPUT_UNIT - end if - - ! determine how many times to repeat '=' character - n = (63 - len_trim(msg))/2 - m = n - if (mod(len_trim(msg),2) == 0) m = m + 1 - - ! convert line to upper case - line = to_upper(msg) - - ! print header based on level - if (verbosity >= level) then - write(UNIT=unit_, FMT='(/1X,A/)') repeat('=', n) // '> ' // & - trim(line) // ' <' // repeat('=', m) - end if - - end subroutine header - -end module output diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 9474f58896..dd1c89669d 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -8,7 +8,6 @@ module particle_restart use hdf5_interface, only: file_open, file_close, read_dataset, HID_T use mgxs_interface, only: energy_bin_avg use nuclide_header, only: micro_xs, n_nuclides - use output, only: print_particle use particle_header use photon_header, only: micro_photon_xs, n_elements use random_lcg, only: set_particle_seed @@ -37,6 +36,11 @@ contains interface subroutine set_micro_xs() bind(C) end subroutine + + subroutine print_particle(p) bind(C) + import Particle + type(Particle), intent(in) :: p + end subroutine end interface err = 0 From 2d7c649ed904445a42a10d6f56d2a319834c30bf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Feb 2019 13:59:40 -0600 Subject: [PATCH 6/7] Make sure volume calc domain IDs are unique --- src/volume_calc.cpp | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index d13057799a..a35a7bd19b 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -24,6 +24,7 @@ #include // for copy #include // for pow, sqrt #include +#include namespace openmc { @@ -58,6 +59,15 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node) { lower_left_ = get_node_array(node, "lower_left"); upper_right_ = get_node_array(node, "upper_right"); n_samples_ = std::stoi(get_node_value(node, "samples")); + + // Ensure there are no duplicates by copying elements to a set and then + // comparing the length with the original vector + std::unordered_set unique_ids(domain_ids_.cbegin(), domain_ids_.cend()); + if (unique_ids.size() != domain_ids_.size()) { + throw std::runtime_error{"Domain IDs for a volume calculation " + "must be unique."}; + } + } std::vector VolumeCalculation::execute() const @@ -79,7 +89,7 @@ std::vector VolumeCalculation::execute() const i_end = i_start + min_samples; } -#pragma omp parallel + #pragma omp parallel { // Variables that are private to each thread std::vector> indices(n); @@ -90,7 +100,7 @@ std::vector VolumeCalculation::execute() const prn_set_stream(STREAM_VOLUME); // Sample locations and count hits -#pragma omp for + #pragma omp for for (int i = i_start; i < i_end; i++) { set_particle_seed(i); @@ -114,6 +124,7 @@ std::vector VolumeCalculation::execute() const for (int i_domain = 0; i_domain < n; i_domain++) { if (model::materials[i_material]->id_ == domain_ids_[i_domain]) { this->check_hit(i_material, indices[i_domain], hits[i_domain]); + break; } } } @@ -122,6 +133,7 @@ std::vector VolumeCalculation::execute() const for (int i_domain=0; i_domain < n; i_domain++) { if (model::cells[p.coord[level].cell]->id_ == domain_ids_[i_domain]) { this->check_hit(i_material, indices[i_domain], hits[i_domain]); + break; } } } @@ -130,6 +142,7 @@ std::vector VolumeCalculation::execute() const for (int i_domain = 0; i_domain < n; ++i_domain) { if (model::universes[p.coord[level].universe]->id_ == domain_ids_[i_domain]) { check_hit(i_material, indices[i_domain], hits[i_domain]); + break; } } } @@ -141,9 +154,9 @@ std::vector VolumeCalculation::execute() const // so we have to manually reduce them #ifdef _OPENMP -#pragma omp for ordered schedule(static) + #pragma omp for ordered schedule(static) for (int i = 0; i < omp_get_num_threads(); ++i) { -#pragma omp ordered + #pragma omp ordered for (int i_domain = 0; i_domain < n; ++i_domain) { for (int j = 0; j < indices[i_domain].size(); ++j) { // Check if this material has been added to the master list and if so, From 248a0fbc8014b6d6a28502ac473ec051b32f834c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Feb 2019 14:51:41 -0600 Subject: [PATCH 7/7] Fix non-OpenMP volume calcs, move a brace --- src/volume_calc.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index a35a7bd19b..67c22e2447 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -40,7 +40,8 @@ std::vector volume_calcs; // VolumeCalculation implementation //============================================================================== -VolumeCalculation::VolumeCalculation(pugi::xml_node node) { +VolumeCalculation::VolumeCalculation(pugi::xml_node node) +{ // Read domain type (cell, material or universe) std::string domain_type = get_node_value(node, "domain_type"); if (domain_type == "cell") { @@ -179,6 +180,7 @@ std::vector VolumeCalculation::execute() const } #else master_indices = indices; + master_hits = hits; #endif prn_set_stream(STREAM_TRACKING);