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