crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1164 from paulromano/cpp-volume

Browse source 
Translate volume calculations to C++
This commit is contained in:
Sterling Harper 2019-02-15 10:23:09 -05:00 • committed by GitHub
commit 5af7510c40
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
17 changed files with 527 additions and 827 deletions
Download patch file Download diff file Expand all files Collapse all files

5
CMakeLists.txt
View file

@ -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
@ -338,12 +337,9 @@ 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
src/volume_calc.F90
src/volume_header.F90
src/xml_interface.F90
src/tallies/tally.F90
src/tallies/tally_derivative_header.F90
@ -433,6 +429,7 @@ add_library(libopenmc SHARED
src/tallies/trigger.cpp
src/timer.cpp
src/thermal.cpp
src/volume_calc.cpp
src/wmp.cpp
src/xml_interface.cpp
src/xsdata.cpp)

8
include/openmc/hdf5_interface.h
View file

@ -391,6 +391,14 @@ write_attribute(hid_t obj_id, const char* name, const std::vector<T>& buffer)
write_attr(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data());
}
inline void
write_attribute(hid_t obj_id, const char* name, Position r)
{
std::array<double, 3> buffer {r.x, r.y, r.z};
write_attribute(obj_id, name, buffer);
}
//==============================================================================
// Templates/overloads for write_dataset

74
include/openmc/volume_calc.h Normal file
View file

@ -0,0 +1,74 @@
#ifndef VOLUME_CALC_H
#define VOLUME_CALC_H
#include "openmc/position.h"
#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"
#include <string>
#include <vector>
namespace openmc {
//==============================================================================
// Volume calculation class
//==============================================================================
class VolumeCalculation {
public:
// Aliases, types
struct Result {
std::array<double, 2> volume; //!< Mean/standard deviation of volume
std::vector<int> nuclides; //!< Index of nuclides
std::vector<double> atoms; //!< Number of atoms for each nuclide
std::vector<double> uncertainty; //!< Uncertainty on number of atoms
}; // Results for a single domain
// Constructors
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<Result> execute() 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<Result>& 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<int> domain_ids_; //!< IDs of domains to find volumes of
private:
//! \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<int>& indices,
std::vector<int>& hits) const;
};
//==============================================================================
// Global variables
//==============================================================================
namespace model {
extern std::vector<VolumeCalculation> volume_calcs;
}
} // namespace openmc
#endif // VOLUME_CALC_H

4
src/api.F90
View file

@ -67,7 +67,6 @@ contains
use surface_header
use tally_filter_header
use tally_header
use volume_header
interface
subroutine free_memory_source() bind(C)
@ -88,6 +87,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

5
src/cell.cpp
View file

@ -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);}
}

2
src/geometry.F90
View file

@ -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

6
src/geometry_header.F90
View file

@ -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

32
src/input_xml.F90
View file

@ -29,7 +29,6 @@ module input_xml
use tally_derivative_header
use tally_filter_header
use tally_filter
use volume_header
use xml_interface
implicit none
@ -56,9 +55,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
@ -90,34 +86,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
!===============================================================================

96
src/output.F90
View file

@ -1,96 +0,0 @@
module output
use, intrinsic :: ISO_C_BINDING
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
use math, only: t_percentile
use message_passing, only: master, n_procs
use mgxs_interface
use nuclide_header
use particle_header, only: LocalCoord, 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
!===============================================================================
! 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.
!===============================================================================
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

6
src/particle_restart.F90
View file

@ -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

3
src/reaction_header.F90
View file

@ -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

9
src/settings.cpp
View file

@ -26,6 +26,7 @@
#include "openmc/source.h"
#include "openmc/string_utils.h"
#include "openmc/tallies/trigger.h"
#include "openmc/volume_calc.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -721,7 +722,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")) {
@ -776,9 +780,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());
}
//==============================================================================

3
src/string.F90
View file

@ -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

88
src/timer_header.F90
View file

@ -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

514
src/volume_calc.F90
View file

@ -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

430
src/volume_calc.cpp Normal file
View file

@ -0,0 +1,430 @@
#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 <omp.h>
#endif
#include "xtensor/xadapt.hpp"
#include "xtensor/xview.hpp"
#include <algorithm> // for copy
#include <cmath> // for pow, sqrt
#include <sstream>
#include <unordered_set>
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
namespace model {
std::vector<VolumeCalculation> 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<int>(node, "domain_ids");
lower_left_ = get_node_array<double>(node, "lower_left");
upper_right_ = get_node_array<double>(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<int> 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::Result> VolumeCalculation::execute() const
{
// Shared data that is collected from all threads
int n = domain_ids_.size();
std::vector<std::vector<int>> master_indices(n); // List of material indices for each domain
std::vector<std::vector<int>> 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
std::vector<std::vector<int>> indices(n);
std::vector<std::vector<int>> hits(n);
Particle p;
p.initialize();
prn_set_stream(STREAM_VOLUME);
// Sample 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_material, indices[i_domain], hits[i_domain]);
break;
}
}
}
} 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_material, indices[i_domain], hits[i_domain]);
break;
}
}
}
} 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_material, indices[i_domain], hits[i_domain]);
break;
}
}
}
}
}
// 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
master_indices = indices;
master_hits = hits;
#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<Result> 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<double, 2> atoms({n_nuc, 2}, 0.0);
#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);
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];
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];
double f = static_cast<double>(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<double>(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);
}
}
}
}
return results;
}
void VolumeCalculation::to_hdf5(const std::string& filename,
const std::vector<Result>& 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<std::string> 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<double, 2> 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);
}
void VolumeCalculation::check_hit(int i_material, std::vector<int>& indices,
std::vector<int>& 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
//==============================================================================
// 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.to_hdf5(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;
}
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" void free_memory_volume() { openmc::model::volume_calcs.clear(); }

69
src/volume_header.F90
View file

@ -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