From 764f638d66915a606a64f04170b3a8bee501f622 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 6 Jul 2018 06:36:05 -0500 Subject: [PATCH] Convert Particle type to C++, add templated HDF5 functions --- CMakeLists.txt | 1 + include/openmc.h | 2 + src/api.F90 | 4 +- src/cell.cpp | 3 +- src/geometry.F90 | 12 +- src/geometry.h | 6 + src/geometry_header.F90 | 2 +- src/hdf5_interface.cpp | 24 ++- src/hdf5_interface.h | 183 +++++++++++------- src/input_xml.F90 | 14 +- src/lattice.cpp | 28 +-- src/particle.cpp | 225 +++++++++++++++++++++ src/particle.h | 160 +++++++++++++++ src/particle_header.F90 | 398 ++++++++++---------------------------- src/particle_restart.F90 | 4 +- src/photon_physics.F90 | 11 +- src/physics.F90 | 33 ++-- src/physics_mg.F90 | 4 +- src/plot.F90 | 6 +- src/settings.F90 | 2 +- src/simulation.F90 | 4 +- src/simulation.h | 13 ++ src/simulation_header.F90 | 8 +- src/source_header.F90 | 6 +- src/surface.cpp | 24 +-- src/tracking.F90 | 18 +- src/volume_calc.F90 | 4 +- 27 files changed, 734 insertions(+), 465 deletions(-) create mode 100644 src/geometry.h create mode 100644 src/particle.cpp create mode 100644 src/particle.h create mode 100644 src/simulation.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 73b3a3851..ac3b3c814 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -393,6 +393,7 @@ add_library(libopenmc SHARED src/message_passing.cpp src/mgxs.cpp src/mgxs_interface.cpp + src/particle.cpp src/plot.cpp src/pugixml/pugixml_c.cpp src/random_lcg.cpp diff --git a/include/openmc.h b/include/openmc.h index f211b4930..1f919d7c6 100644 --- a/include/openmc.h +++ b/include/openmc.h @@ -14,6 +14,7 @@ extern "C" { double uvw[3]; double E; int delayed_group; + int particle; }; int openmc_calculate_volumes(); @@ -128,6 +129,7 @@ extern "C" { extern char openmc_err_msg[256]; extern double openmc_keff; extern double openmc_keff_std; + extern int32_t gen_per_batch; extern int32_t n_batches; extern int32_t n_cells; extern int32_t n_filters; diff --git a/src/api.F90 b/src/api.F90 index 01bd565b6..cc9120c7b 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -15,7 +15,7 @@ module openmc_api use message_passing use nuclide_header use initialize, only: openmc_init_f - use particle_header, only: Particle + use particle_header use plot, only: openmc_plot_geometry use random_lcg, only: openmc_get_seed, openmc_set_seed use settings @@ -198,7 +198,7 @@ contains logical :: found type(Particle) :: p - call p % initialize() + call particle_initialize(p) p % coord(1) % xyz(:) = xyz p % coord(1) % uvw(:) = [ZERO, ZERO, ONE] call find_cell(p, found) diff --git a/src/cell.cpp b/src/cell.cpp index e4388e701..263191b31 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -329,8 +329,7 @@ Cell::to_hdf5(hid_t cell_group) const } //TODO: Fix the off-by-one indexing. - write_int(cell_group, 0, nullptr, "universe", - &global_universes[universe-1]->id, false); + write_dataset(cell_group, "universe", global_universes[universe-1]->id); // Write the region specification. if (!region.empty()) { diff --git a/src/geometry.F90 b/src/geometry.F90 index a5365dbf6..5bf0012a3 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -3,7 +3,7 @@ module geometry use constants use error, only: fatal_error, warning, write_message use geometry_header - use particle_header, only: LocalCoord, Particle + use particle_header use simulation_header use settings use surface_header @@ -113,13 +113,13 @@ contains logical :: use_search_cells ! use cells provided as argument do j = p % n_coord + 1, MAX_COORD - call p % coord(j) % reset() + call reset_coord(p % coord(j)) end do j = p % n_coord ! Determine universe (if not yet set, use root universe) i_universe = p % coord(j) % universe - if (i_universe == NONE) then + if (i_universe == C_NONE) then p % coord(j) % universe = root_universe i_universe = root_universe end if @@ -336,7 +336,7 @@ contains call find_cell(p, found) if (.not. found) then if (p % alive) then ! Particle may have been killed in find_cell - call p % mark_as_lost("Could not locate particle " & + call particle_mark_as_lost(p, "Could not locate particle " & // trim(to_str(p % id)) // " after crossing a lattice boundary.") return end if @@ -360,7 +360,7 @@ contains ! Search for particle call find_cell(p, found) if (.not. found) then - call p % mark_as_lost("Could not locate particle " // & + call particle_mark_as_lost(p, "Could not locate particle " // & trim(to_str(p % id)) // " after crossing a lattice boundary.") return end if @@ -440,7 +440,7 @@ contains end select LAT_TYPE if (d_lat < ZERO) then - call p % mark_as_lost("Particle " // trim(to_str(p % id)) & + call particle_mark_as_lost(p, "Particle " // trim(to_str(p % id)) & //" had a negative distance to a lattice boundary. d = " & //trim(to_str(d_lat))) end if diff --git a/src/geometry.h b/src/geometry.h new file mode 100644 index 000000000..8837c961a --- /dev/null +++ b/src/geometry.h @@ -0,0 +1,6 @@ +#ifndef GEOMETRY_H +#define GEOMETRY_H + +extern "C" int openmc_root_universe; + +#endif // GEOMETRY_H diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 2c943706d..71c458105 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -288,7 +288,7 @@ module geometry_header end type Cell ! array index of the root universe - integer :: root_universe = -1 + integer(C_INT), bind(C, name='openmc_root_universe') :: root_universe = -1 integer(C_INT32_T), bind(C) :: n_cells ! # of cells integer(C_INT32_T), bind(C) :: n_universes ! # of universes diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index 7133ad51a..9226cb9ad 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -169,9 +169,9 @@ file_open(const char* filename, char mode, bool parallel) } hid_t -file_open(const std::string& filename, char mode, bool parallel=false) +file_open(const std::string& filename, char mode, bool parallel) { - file_open(filename.c_str(), mode, parallel); + return file_open(filename.c_str(), mode, parallel); } void file_close(hid_t file_id) @@ -302,26 +302,24 @@ object_exists(hid_t object_id, const char* name) hid_t open_dataset(hid_t group_id, const char* name) { - if (object_exists(group_id, name)) { - return H5Dopen(group_id, name, H5P_DEFAULT); - } else { + if (!object_exists(group_id, name)) { std::stringstream err_msg; err_msg << "Group \"" << name << "\" does not exist"; fatal_error(err_msg); } + return H5Dopen(group_id, name, H5P_DEFAULT); } hid_t open_group(hid_t group_id, const char* name) { - if (object_exists(group_id, name)) { - return H5Gopen(group_id, name, H5P_DEFAULT); - } else { + if (!object_exists(group_id, name)) { std::stringstream err_msg; err_msg << "Group \"" << name << "\" does not exist"; fatal_error(err_msg); } + return H5Gopen(group_id, name, H5P_DEFAULT); } void @@ -821,4 +819,14 @@ using_mpio_device(hid_t obj_id) return driver == H5FD_MPIO; } +// Specializations of the H5TypeMap template struct +template<> +const hid_t H5TypeMap::type_id = H5T_NATIVE_INT; +template<> +const hid_t H5TypeMap::type_id = H5T_NATIVE_INT64; +template<> +const hid_t H5TypeMap::type_id = H5T_NATIVE_DOUBLE; +template <> +const hid_t H5TypeMap::type_id = H5T_NATIVE_CHAR; + } // namespace openmc diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index b00fe7c87..2face49ea 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -5,6 +5,7 @@ #include "hdf5_hl.h" #include +#include #include #include #include @@ -13,48 +14,28 @@ namespace openmc { -extern "C" bool attribute_exists(hid_t obj_id, const char* name); -extern "C" size_t attribute_typesize(hid_t obj_id, const char* name); -extern "C" hid_t create_group(hid_t parent_id, const char* name); -hid_t create_group(hid_t parent_id, const std::string& name); -extern "C" void close_dataset(hid_t dataset_id); -extern "C" void close_group(hid_t group_id); -extern "C" int dataset_ndims(hid_t dset); -extern "C" size_t dataset_typesize(hid_t dset); -extern "C" hid_t file_open(const char* filename, char mode, bool parallel); -hid_t file_open(const std::string& filename, char mode, bool parallel); -extern "C" void file_close(hid_t file_id); -extern "C" void get_name(hid_t obj_id, char* name); -extern "C" int get_num_datasets(hid_t group_id); -extern "C" int get_num_groups(hid_t group_id); -extern "C" void get_datasets(hid_t group_id, char* name[]); -extern "C" void get_groups(hid_t group_id, char* name[]); -extern "C" void get_shape(hid_t obj_id, hsize_t* dims); -extern "C" void get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims); -extern "C" bool object_exists(hid_t object_id, const char* name); -extern "C" hid_t open_dataset(hid_t group_id, const char* name); -extern "C" hid_t open_group(hid_t group_id, const char* name); -bool using_mpio_device(hid_t obj_id); +//============================================================================== +// Low-level internal functions +//============================================================================== void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id, const void* buffer); -extern "C" void read_attr_double(hid_t obj_id, const char* name, double* buffer); -extern "C" void read_attr_int(hid_t obj_id, const char* name, int* buffer); -extern "C" void read_attr_string(hid_t obj_id, const char* name, size_t slen, - char* buffer); - +void write_attr(hid_t obj_id, int ndim, const hsize_t* dims, const char* name, + hid_t mem_type_id, const void* buffer); void read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id, void* buffer, bool indep); -extern "C" void read_double(hid_t obj_id, const char* name, double* buffer, - bool indep); -extern "C" void read_int(hid_t obj_id, const char* name, int* buffer, - bool indep); -extern "C" void read_llong(hid_t obj_id, const char* name, long long* buffer, - bool indep); -extern "C" void read_string(hid_t obj_id, const char* name, size_t slen, - char* buffer, bool indep); -extern "C" void read_complex(hid_t obj_id, const char* name, - double _Complex* buffer, bool indep); +void write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name, + hid_t mem_type_id, const void* buffer, bool indep); +bool using_mpio_device(hid_t obj_id); + +//============================================================================== +// Normal functions that are used to read/write files +//============================================================================== + +hid_t create_group(hid_t parent_id, const std::string& name); +hid_t file_open(const std::string& filename, char mode, bool parallel=false); +void write_string(hid_t group_id, const char* name, const std::string& buffer, + bool indep); void read_nd_vector(hid_t obj_id, const char* name, std::vector& result, @@ -89,51 +70,115 @@ read_nd_vector(hid_t obj_id, const char* name, std::vector > > > >& result, bool must_have = false); -extern "C" void read_tally_results(hid_t group_id, hsize_t n_filter, +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" { + bool attribute_exists(hid_t obj_id, const char* name); + size_t attribute_typesize(hid_t obj_id, const char* name); + hid_t create_group(hid_t parent_id, const char* name); + void close_dataset(hid_t dataset_id); + void close_group(hid_t group_id); + int dataset_ndims(hid_t dset); + size_t dataset_typesize(hid_t dset); + hid_t file_open(const char* filename, char mode, bool parallel); + void file_close(hid_t file_id); + void get_name(hid_t obj_id, char* name); + int get_num_datasets(hid_t group_id); + int get_num_groups(hid_t group_id); + void get_datasets(hid_t group_id, char* name[]); + void get_groups(hid_t group_id, char* name[]); + void get_shape(hid_t obj_id, hsize_t* dims); + void get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims); + bool object_exists(hid_t object_id, const char* name); + hid_t open_dataset(hid_t group_id, const char* name); + hid_t open_group(hid_t group_id, const char* name); + void read_attr_double(hid_t obj_id, const char* name, double* buffer); + void read_attr_int(hid_t obj_id, const char* name, int* buffer); + void read_attr_string(hid_t obj_id, const char* name, size_t slen, + char* buffer); + void read_complex(hid_t obj_id, const char* name, + double _Complex* buffer, bool indep); + void read_double(hid_t obj_id, const char* name, double* buffer, + bool indep); + void read_int(hid_t obj_id, const char* name, int* buffer, + bool indep); + void read_llong(hid_t obj_id, const char* name, long long* buffer, + bool indep); + void read_string(hid_t obj_id, const char* name, size_t slen, + char* buffer, bool indep); + + + void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results); - -void write_attr(hid_t obj_id, int ndim, const hsize_t* dims, const char* name, - hid_t mem_type_id, const void* buffer); -extern "C" void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims, + void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims, const char* name, const double* buffer); -extern "C" void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims, + void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims, const char* name, const int* buffer); -extern "C" void write_attr_string(hid_t obj_id, const char* name, const char* buffer); - - -void write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name, - hid_t mem_type_id, const void* buffer, bool indep); -extern "C" void write_double(hid_t group_id, int ndim, const hsize_t* dims, + void write_attr_string(hid_t obj_id, const char* name, const char* buffer); + void write_double(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const double* buffer, bool indep); -extern "C" void write_int(hid_t group_id, int ndim, const hsize_t* dims, + void write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const int* buffer, bool indep); -extern "C" void write_llong(hid_t group_id, int ndim, const hsize_t* dims, + void write_llong(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const long long* buffer, bool indep); - -extern "C" void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen, + void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen, const char* name, char const* buffer, bool indep); -void write_string(hid_t group_id, const char* name, const std::string& buffer, bool indep); - - -extern "C" void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, + void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results); +} // extern "C" -template void -write_int(hid_t group_id, char const *name, - const std::array &buffer, bool indep) +//============================================================================== +// Template struct used to map types to HDF5 datatype IDs, which are stored +// using the type hid_t. By having a single static data member, the template can +// be specialized for each type we know of. The specializations appear in the +// .cpp file since they are definitions. +//============================================================================== + +template +struct H5TypeMap { static const hid_t type_id; }; + +//============================================================================== +// Template functions used to provide simple interface to lower-level functions +//============================================================================== + +template +void write_attribute(hid_t obj_id, const char* name, T buffer) { - hsize_t dims[1] {array_len}; - write_dataset(group_id, 1, dims, name, H5T_NATIVE_INT, buffer.data(), indep); + write_attr(obj_id, name, 0, nullptr, H5TypeMap::type_id, &buffer); } - -template void -write_double(hid_t group_id, char const *name, - const std::array &buffer, bool indep) +template<> inline void +write_attribute(hid_t obj_id, const char* name, const char* buffer) { - hsize_t dims[1] {array_len}; - write_dataset(group_id, 1, dims, name, H5T_NATIVE_DOUBLE, - buffer.data(), indep); + write_attr_string(obj_id, name, buffer); +} + +template inline void +write_attribute(hid_t obj_id, const char* name, const std::array& buffer) +{ + hsize_t dims[] {N}; + write_attr(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data()); +} + +template inline void +write_dataset(hid_t obj_id, const char* name, T buffer) +{ + write_dataset(obj_id, 0, nullptr, name, H5TypeMap::type_id, &buffer, false); +} + +template<> inline void +write_dataset(hid_t obj_id, const char* name, const char* buffer) +{ + write_string(obj_id, name, buffer, false); +} + +template inline void +write_dataset(hid_t obj_id, const char* name, const std::array& buffer) +{ + hsize_t dims[] {N}; + write_dataset(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data(), false); } } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index b9884bcde..e2feb6848 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1232,8 +1232,8 @@ contains end if ! Read cell temperatures. If the temperature is not specified, set it to - ! ERROR_REAL for now. During initialization we'll replace ERROR_REAL with - ! the temperature from the material data. + ! a negative number for now. During initialization we'll replace + ! negatives with the temperature from the material data. if (check_for_node(node_cell, "temperature")) then n = node_word_count(node_cell, "temperature") if (n > 0) then @@ -1258,11 +1258,11 @@ contains c % sqrtkT(:) = sqrt(K_BOLTZMANN * c % sqrtkT(:)) else allocate(c % sqrtkT(1)) - c % sqrtkT(1) = ERROR_REAL + c % sqrtkT(1) = -1.0 end if else allocate(c % sqrtkT(1)) - c % sqrtkT = ERROR_REAL + c % sqrtkT = -1.0 end if ! Add cell to dictionary @@ -1580,7 +1580,7 @@ contains if (check_for_node(node_mat, "temperature")) then call get_node_value(node_mat, "temperature", material_temps(i)) else - material_temps(i) = ERROR_REAL + material_temps(i) = -1.0 end if ! Get pointer to density element @@ -3847,7 +3847,7 @@ contains do i = 1, n_cells ! Ignore non-normal cells and cells with defined temperature. if (cells(i) % material(1) == NONE) cycle - if (cells(i) % sqrtkT(1) /= ERROR_REAL) cycle + if (cells(i) % sqrtkT(1) >= ZERO) cycle ! Set the number of temperatures equal to the number of materials. deallocate(cells(i) % sqrtkT) @@ -3863,7 +3863,7 @@ contains ! Use material default or global default temperature i_material = cells(i) % material(j) - if (material_temps(i_material) /= ERROR_REAL) then + if (material_temps(i_material) >= ZERO) then cells(i) % sqrtkT(j) = sqrt(K_BOLTZMANN * & material_temps(i_material)) else diff --git a/src/lattice.cpp b/src/lattice.cpp index 959ab932e..c95762a63 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -119,9 +119,9 @@ Lattice::to_hdf5(hid_t lattices_group) const if (outer != NO_OUTER_UNIVERSE) { int32_t outer_id = global_universes[outer]->id; - write_int(lat_group, 0, nullptr, "outer", &outer_id, false); + write_dataset(lat_group, "outer", outer_id); } else { - write_int(lat_group, 0, nullptr, "outer", &outer, false); + write_dataset(lat_group, "outer", outer); } // Call subclass-overriden function to fill in other details. @@ -348,16 +348,16 @@ RectLattice::to_hdf5_inner(hid_t lat_group) const // Write basic lattice information. write_string(lat_group, "type", "rectangular", false); if (is_3d) { - write_double(lat_group, "pitch", pitch, false); - write_double(lat_group, "lower_left", lower_left, false); - write_int(lat_group, "dimension", n_cells, false); + write_dataset(lat_group, "pitch", pitch); + write_dataset(lat_group, "lower_left", lower_left); + write_dataset(lat_group, "dimension", n_cells); } else { std::array pitch_short {{pitch[0], pitch[1]}}; - write_double(lat_group, "pitch", pitch_short, false); + write_dataset(lat_group, "pitch", pitch_short); std::array ll_short {{lower_left[0], lower_left[1]}}; - write_double(lat_group, "lower_left", ll_short, false); + write_dataset(lat_group, "lower_left", ll_short); std::array nc_short {{n_cells[0], n_cells[1]}}; - write_int(lat_group, "dimension", nc_short, false); + write_dataset(lat_group, "dimension", nc_short); } // Write the universe ids. The convention here is to switch the ordering on @@ -826,16 +826,16 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const { // Write basic lattice information. write_string(lat_group, "type", "hexagonal", false); - write_int(lat_group, 0, nullptr, "n_rings", &n_rings, false); - write_int(lat_group, 0, nullptr, "n_axial", &n_axial, false); + write_dataset(lat_group, "n_rings", n_rings); + write_dataset(lat_group, "n_axial", n_axial); if (is_3d) { - write_double(lat_group, "pitch", pitch, false); - write_double(lat_group, "center", center, false); + write_dataset(lat_group, "pitch", pitch); + write_dataset(lat_group, "center", center); } else { std::array pitch_short {{pitch[0]}}; - write_double(lat_group, "pitch", pitch_short, false); + write_dataset(lat_group, "pitch", pitch_short); std::array center_short {{center[0], center[1]}}; - write_double(lat_group, "center", center_short, false); + write_dataset(lat_group, "center", center_short); } // Write the universe ids. diff --git a/src/particle.cpp b/src/particle.cpp new file mode 100644 index 000000000..a7f91e0f3 --- /dev/null +++ b/src/particle.cpp @@ -0,0 +1,225 @@ +#include "particle.h" + +#include +#include + +#include "constants.h" +#include "error.h" +#include "hdf5_interface.h" +#include "openmc.h" +#include "simulation.h" + +namespace openmc { + +//============================================================================== +// LocalCoord implementation +//============================================================================== + +void +LocalCoord::reset() +{ + cell = 0; + universe = 0; + lattice = 0; + lattice_x = 0; + lattice_y = 0; + rotated = false; +} + +//============================================================================== +// Particle implementation +//============================================================================== + +void +Particle::clear() +{ + // reset any coordinate levels + for (int i=0; iparticle; + wgt = src->wgt; + last_wgt = src->wgt; + std::copy(src->xyz, src->xyz + 3, coord[0].xyz); + std::copy(src->uvw, src->uvw + 3, coord[0].uvw); + std::copy(src->xyz, src->xyz + 3, last_xyz_current); + std::copy(src->xyz, src->xyz + 3, last_xyz); + std::copy(src->uvw, src->uvw + 3, last_uvw); + if (run_CE) { + E = src->E; + g = 0; + } else { + g = static_cast(src->E); + last_g = static_cast(src->E); + E = energy_bin_avg[g - 1]; + } + last_E = E; +} + +void +Particle::mark_as_lost(const char* message) +{ + // Print warning and write lost particle file + warning(message); + write_restart(); + + // Increment number of lost particles + alive = false; +#pragma omp atomic + openmc_n_lost_particles += 1; + + // Count the total number of simulated particles (on this processor) + auto n = openmc_current_batch * gen_per_batch * openmc_work; + + // Abort the simulation if the maximum number of lost particles has been + // reached + if (openmc_n_lost_particles >= MAX_LOST_PARTICLES && + openmc_n_lost_particles >= REL_MAX_LOST_PARTICLES*n) { + fatal_error("Maximum number of lost particles has been reached."); + } +} + +void +Particle::write_restart() +{ + // Dont write another restart file if in particle restart mode + if (openmc_run_mode == RUN_MODE_PARTICLE) return; + + // Set up file name + // TODO: Add path_output when available on C++ side + std::stringstream filename; + filename << "particle_" << openmc_current_batch << '_' << id << ".h5"; + +#pragma omp critical (WriteParticleRestart) + { + // Create file + hid_t file_id = file_open(filename.str(), 'w'); + + // Write filetype and version info + write_attribute(file_id, "filetype", "particle restart"); + write_attribute(file_id, "version", VERSION_PARTICLE_RESTART); + write_attribute(file_id, "openmc_version", VERSION); +#ifdef GIT_SHA1 + write_attr_string(file_id, "git_sha1", GIT_SHA1); +#endif + + // Write data to file + write_dataset(file_id, "current_batch", openmc_current_batch); + write_dataset(file_id, "generations_per_batch", gen_per_batch); + write_dataset(file_id, "current_generation", openmc_current_gen); + write_dataset(file_id, "n_particles", n_particles); + switch (openmc_run_mode) { + case RUN_MODE_FIXEDSOURCE: + write_dataset(file_id, "run_mode", "fixed source"); + break; + case RUN_MODE_EIGENVALUE: + write_dataset(file_id, "run_mode", "eigenvalue"); + break; + case RUN_MODE_PARTICLE: + write_dataset(file_id, "run_mode", "particle restart"); + break; + } + write_dataset(file_id, "id", id); + write_dataset(file_id, "type", type); + + // Get pointer to source bank + Bank* src; + int64_t n; + openmc_source_bank(&src, &n); + + int64_t i = openmc_current_work; + write_dataset(file_id, "weight", src[i-1].wgt); + write_dataset(file_id, "energy", src[i-1].E); + hsize_t dims[] {3}; + write_double(file_id, 1, dims, "xyz", src[i-1].xyz, false); + write_double(file_id, 1, dims, "uvw", src[i-1].uvw, false); + + // Close file + file_close(file_id); + } // #pragma omp critical +} + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +void reset_coord(LocalCoord* c) { c->reset(); } +void particle_clear(Particle* p) { p->clear(); } +void particle_create_secondary(Particle* p, const double* uvw, double E, + int type, bool run_CE) +{ + p->create_secondary(uvw, E, type, run_CE); +} +void particle_initialize(Particle* p) { p->initialize(); } +void particle_from_source(Particle* p, const Bank* src, bool run_CE, + const double* energy_bin_avg) +{ + p->from_source(src, run_CE, energy_bin_avg); +} +void particle_mark_as_lost(Particle* p, const char* message) +{ + p->mark_as_lost(message); +} +void particle_write_restart(Particle* p) { p->write_restart(); } + +} // namespace openmc diff --git a/src/particle.h b/src/particle.h new file mode 100644 index 000000000..81089f606 --- /dev/null +++ b/src/particle.h @@ -0,0 +1,160 @@ +#ifndef PARTICLE_H +#define PARTICLE_H + +#include +#include + +#include "openmc.h" + +namespace openmc { + +//============================================================================== +// Constants +//============================================================================== + +constexpr int MAX_DELAYED_GROUPS {8}; +constexpr int MAX_SECONDARY {1000}; +constexpr int NEUTRON {1}; +constexpr int MAX_LOST_PARTICLES {10}; +constexpr double REL_MAX_LOST_PARTICLES {1.0e-6}; + +extern "C" { + + struct LocalCoord { + int cell; + int universe; + int lattice; + int lattice_x; + int lattice_y; + int lattice_z; + double xyz[3]; //!< particle position + double uvw[3]; //!< particle direction + bool rotated {false}; //!< Is the level rotated? + + //! clear data from a single coordinate level + void reset(); + }; + + //============================================================================ + //! State of a particle being transported through geometry + //============================================================================ + + struct Particle { + int64_t id; //!< Unique ID + int type; //!< Particle type (n, p, e, etc.) + + int n_coord; //!< number of current coordinate levels + int cell_instance; //!< offset for distributed properties + LocalCoord coord[MAX_COORD]; //!< coordinates for all levels + + // Particle coordinates before crossing a surface + int last_n_coord; //!< number of current coordinates + int last_cell[MAX_COORD]; //!< coordinates for all levels + + // Energy data + double E; //!< post-collision energy in eV + double last_E; //!< pre-collision energy in eV + int g; //!< post-collision energy group (MG only) + int last_g; //!< pre-collision energy group (MG only) + + // Other physical data + double wgt; //!< particle weight + double mu; //!< angle of scatter + bool alive; //!< is particle alive? + + // Other physical data + double last_xyz_current[3]; //!< coordinates of the last collision or + //!< reflective/periodic surface crossing for + //!< current tallies + double last_xyz[3]; //!< previous coordinates + double last_uvw[3]; //!< previous direction coordinates + double last_wgt; //!< pre-collision particle weight + double absorb_wgt; //!< weight absorbed for survival biasing + + // What event took place + bool fission; //!< did particle cause implicit fission + int event; //!< scatter, absorption + int event_nuclide; //!< index in nuclides array + int event_MT; //!< reaction MT + int delayed_group; //!< delayed group + + // Post-collision physical data + int n_bank; //!< number of fission sites banked + double wgt_bank; //!< weight of fission sites banked + int n_delayed_bank[MAX_DELAYED_GROUPS]; //!< number of delayed fission + //!< sites banked + + // Indices for various arrays + int surface; //!< index for surface particle is on + int cell_born; //!< index for cell particle was born in + int material; //!< index for current material + int last_material; //!< index for last material + + // Temperature of current cell + double sqrtkT; //!< sqrt(k_Boltzmann * temperature) in eV + double last_sqrtkT; //!< last temperature + + // Statistical data + int n_collision; //!< # of collisions + + // Track output + bool write_track {false}; + + // Secondary particles created + int64_t n_secondary {}; + Bank secondary_bank[MAX_SECONDARY]; + + //! resets all coordinate levels for the particle + void clear(); + + //! create a secondary particle + // + //! stores the current phase space attributes of the particle in the + //! secondary bank and increments the number of sites in the secondary bank. + //! @param uvw Direction of the secondary particle + //! @param E Energy of the secondary particle in [eV] + //! @param type Particle type + //! @param run_CE Whether continuous-energy data is being used + void create_secondary(const double* uvw, double E, int type, bool run_CE); + + //! sets default attributes for a particle + void initialize(); + + //! initialize from a source site + // + //! initializes a particle from data stored in a source site. The source + //! site may have been produced from an external source, from fission, or + //! simply as a secondary particle. + //! @param src Source site data + //! @param run_CE Whether continuous-energy data is being used + //! @param energy_bin_avg An array of energy group bin averages + void from_source(const Bank* src, bool run_CE, const double* energy_bin_avg); + + //! mark a particle as lost and create a particle restart file + //! @param message A warning message to display + void mark_as_lost(const char* message); + + //! create a particle restart HDF5 file + void write_restart(); + }; + + + //============================================================================ + // Fortran compatibility functions + //============================================================================ + + void reset_coord(LocalCoord* c); + void particle_clear(Particle* p); + void particle_create_secondary(Particle* p, const double* uvw, double E, + int type, bool run_CE); + void particle_initialize(Particle* p); + void particle_from_source(Particle* p, const Bank* src, bool run_CE, + const double* energy_bin_avg); + void particle_mark_as_lost(Particle* p, const char* message); + void particle_write_restart(Particle* p); + +} // extern "C" + +} // namespace openmc + +#endif // PARTICLE_H diff --git a/src/particle_header.F90 b/src/particle_header.F90 index a04379c7e..df5deff8c 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -1,41 +1,34 @@ module particle_header - use bank_header, only: Bank, source_bank + use, intrinsic :: ISO_C_BINDING + + use bank_header, only: Bank use constants - use error, only: fatal_error, warning - use hdf5_interface - use settings - use simulation_header - use string, only: to_str + use string, only: to_c_string implicit none - private - !=============================================================================== ! LOCALCOORD describes the location of a particle local to a single ! universe. When the geometry consists of nested universes, a particle will have ! a list of coordinates in each level !=============================================================================== - type, public :: LocalCoord - + type, bind(C) :: LocalCoord ! Indices in various arrays for this level - integer :: cell = NONE - integer :: universe = NONE - integer :: lattice = NONE - integer :: lattice_x = NONE - integer :: lattice_y = NONE - integer :: lattice_z = NONE + integer(C_INT) :: cell = NONE + integer(C_INT) :: universe = NONE + integer(C_INT) :: lattice = NONE + integer(C_INT) :: lattice_x = NONE + integer(C_INT) :: lattice_y = NONE + integer(C_INT) :: lattice_z = NONE ! Particle position and direction for this level - real(8) :: xyz(3) - real(8) :: uvw(3) + real(C_DOUBLE) :: xyz(3) + real(C_DOUBLE) :: uvw(3) ! Is this level rotated? - logical :: rotated = .false. - contains - procedure :: reset => reset_coord + logical(C_BOOL) :: rotated = .false. end type LocalCoord !=============================================================================== @@ -43,311 +36,126 @@ module particle_header ! geometry !=============================================================================== - type, public :: Particle + type, bind(C) :: Particle ! Basic data - integer(8) :: id ! Unique ID - integer :: type ! Particle type (n, p, e, etc) + integer(C_INT64_T) :: id ! Unique ID + integer(C_INT) :: type ! Particle type (n, p, e, etc) ! Particle coordinates - integer :: n_coord ! number of current coordinates - integer :: cell_instance ! offset for distributed properties + integer(C_INT) :: n_coord ! number of current coordinates + integer(C_INT) :: cell_instance ! offset for distributed properties type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels ! Particle coordinates before crossing a surface - integer :: last_n_coord ! number of current coordinates - integer :: last_cell(MAX_COORD) ! coordinates for all levels + integer(C_INT) :: last_n_coord ! number of current coordinates + integer(C_INT) :: last_cell(MAX_COORD) ! coordinates for all levels ! Energy Data - real(8) :: E ! post-collision energy - real(8) :: last_E ! pre-collision energy - integer :: g ! post-collision energy group (MG only) - integer :: last_g ! pre-collision energy group (MG only) + real(C_DOUBLE) :: E ! post-collision energy + real(C_DOUBLE) :: last_E ! pre-collision energy + integer(C_INT) :: g ! post-collision energy group (MG only) + integer(C_INT) :: last_g ! pre-collision energy group (MG only) ! Other physical data - real(8) :: wgt ! particle weight - real(8) :: mu ! angle of scatter - logical :: alive ! is particle alive? + real(C_DOUBLE) :: wgt ! particle weight + real(C_DOUBLE) :: mu ! angle of scatter + logical(C_BOOL) :: alive ! is particle alive? ! Pre-collision physical data - real(8) :: last_xyz_current(3) ! coordinates of the last collision or - ! reflective/periodic surface crossing - ! for current tallies - real(8) :: last_xyz(3) ! previous coordinates - real(8) :: last_uvw(3) ! previous direction coordinates - real(8) :: last_wgt ! pre-collision particle weight - real(8) :: absorb_wgt ! weight absorbed for survival biasing + real(C_DOUBLE) :: last_xyz_current(3) ! coordinates of the last collision or + ! reflective/periodic surface crossing + ! for current tallies + real(C_DOUBLE) :: last_xyz(3) ! previous coordinates + real(C_DOUBLE) :: last_uvw(3) ! previous direction coordinates + real(C_DOUBLE) :: last_wgt ! pre-collision particle weight + real(C_DOUBLE) :: absorb_wgt ! weight absorbed for survival biasing ! What event last took place - logical :: fission ! did the particle cause implicit fission - integer :: event ! scatter, absorption - integer :: event_nuclide ! index in nuclides array - integer :: event_MT ! reaction MT - integer :: delayed_group ! delayed group + logical(C_BOOL) :: fission ! did the particle cause implicit fission + integer(C_INT) :: event ! scatter, absorption + integer(C_INT) :: event_nuclide ! index in nuclides array + integer(C_INT) :: event_MT ! reaction MT + integer(C_INT) :: delayed_group ! delayed group ! Post-collision physical data - integer :: n_bank ! number of fission sites banked - real(8) :: wgt_bank ! weight of fission sites banked - integer :: n_delayed_bank(MAX_DELAYED_GROUPS) ! number of delayed fission + integer(C_INT) :: n_bank ! number of fission sites banked + real(C_DOUBLE) :: wgt_bank ! weight of fission sites banked + integer(C_INT) :: n_delayed_bank(MAX_DELAYED_GROUPS) ! number of delayed fission ! sites banked ! Indices for various arrays - integer :: surface ! index for surface particle is on - integer :: cell_born ! index for cell particle was born in - integer :: material ! index for current material - integer :: last_material ! index for last material + integer(C_INT) :: surface ! index for surface particle is on + integer(C_INT) :: cell_born ! index for cell particle was born in + integer(C_INT) :: material ! index for current material + integer(C_INT) :: last_material ! index for last material ! Temperature of the current cell - real(8) :: sqrtkT ! sqrt(k_Boltzmann * temperature) in eV - real(8) :: last_sqrtKT ! last temperature + real(C_DOUBLE) :: sqrtkT ! sqrt(k_Boltzmann * temperature) in eV + real(C_DOUBLE) :: last_sqrtKT ! last temperature ! Statistical data - integer :: n_collision ! # of collisions + integer(C_INT) :: n_collision ! # of collisions ! Track output - logical :: write_track = .false. + logical(C_BOOL) :: write_track = .false. ! Secondary particles created - integer(8) :: n_secondary = 0 - type(Bank) :: secondary_bank(MAX_SECONDARY) - - contains - procedure :: clear - procedure :: create_secondary - procedure :: initialize - procedure :: initialize_from_source - procedure :: mark_as_lost - procedure :: write_restart + integer(C_INT64_T) :: n_secondary = 0 + type(Bank) :: secondary_bank(MAX_SECONDARY) end type Particle + interface + subroutine reset_coord(c) bind(C) + import LocalCoord + type(LocalCoord), intent(inout) :: c + end subroutine reset_coord + + subroutine particle_clear(p) bind(C) + import Particle + type(Particle), intent(inout) :: p + end subroutine particle_clear + + subroutine particle_create_secondary(p, uvw, E, type, run_CE) bind(C) + import Particle, C_DOUBLE, C_INT, C_BOOL + type(Particle), intent(inout) :: p + real(C_DOUBLE), intent(in) :: uvw(3) + real(C_DOUBLE), value :: E + integer(C_INT), value :: type + logical(C_BOOL), value :: run_CE + end subroutine particle_create_secondary + + subroutine particle_initialize(p) bind(C) + import Particle + type(Particle), intent(inout) :: p + end subroutine particle_initialize + + subroutine particle_from_source(p, src, run_CE, energy_bin_avg) bind(C) + import Particle, Bank, C_BOOL, C_DOUBLE + type(Particle), intent(inout) :: p + type(Bank), intent(in) :: src + logical(C_BOOL), value :: run_CE + real(C_DOUBLE), intent(in) :: energy_bin_avg(*) + end subroutine particle_from_source + + subroutine particle_mark_as_lost_c(p, message) bind(C, name='particle_mark_as_lost') + import Particle, C_CHAR + type(Particle), intent(in) :: p + character(kind=C_CHAR), intent(in) :: message(*) + end subroutine particle_mark_as_lost_c + + subroutine particle_write_restart(p) bind(C) + import Particle + type(Particle), intent(in) :: p + end subroutine particle_write_restart + end interface + contains -!=============================================================================== -! RESET_COORD clears data from a single coordinate level -!=============================================================================== + subroutine particle_mark_as_lost(this, message) + type(Particle), intent(inout) :: this + character(*) :: message - elemental subroutine reset_coord(this) - class(LocalCoord), intent(inout) :: this - - this % cell = NONE - this % universe = NONE - this % lattice = NONE - this % lattice_x = NONE - this % lattice_y = NONE - this % lattice_z = NONE - this % rotated = .false. - - end subroutine reset_coord - -!=============================================================================== -! CLEAR_PARTICLE resets all coordinate levels for the particle -!=============================================================================== - - subroutine clear(this) - class(Particle) :: this - - integer :: i - - ! remove any coordinate levels - do i = 1, MAX_COORD - call this % coord(i) % reset() - end do - end subroutine clear - -!=============================================================================== -! CREATE_SECONDARY stores the current phase space attributes of the particle in -! the secondary bank and increments the number of sites in the secondary bank. -!=============================================================================== - - subroutine create_secondary(this, uvw, E, type, run_CE) - class(Particle), intent(inout) :: this - real(8), intent(in) :: uvw(3) - real(8), intent(in) :: E - integer, intent(in) :: type - logical, intent(in) :: run_CE - - integer(8) :: n - - ! Check to make sure that the hard-limit on secondary particles is not - ! exceeded. - if (this % n_secondary == MAX_SECONDARY) then - call fatal_error("Too many secondary particles created.") - end if - - n = this % n_secondary + 1 - this % secondary_bank(n) % particle = type - this % secondary_bank(n) % wgt = this % wgt - this % secondary_bank(n) % xyz(:) = this % coord(1) % xyz - this % secondary_bank(n) % uvw(:) = uvw - this % secondary_bank(n) % E = E - if (.not. run_CE) then - this % secondary_bank(n) % E = real(this % g, 8) - end if - this % n_secondary = n - - end subroutine create_secondary - -!=============================================================================== -! INITIALIZE sets default attributes for a particle from the source bank -!=============================================================================== - - subroutine initialize(this) - - class(Particle) :: this - - ! Clear coordinate lists - call this % clear() - - ! Set particle to neutron that's alive - this % type = NEUTRON - this % alive = .true. - - ! clear attributes - this % surface = NONE - this % cell_born = NONE - this % material = NONE - this % last_material = NONE - this % last_sqrtkT = NONE - this % wgt = ONE - this % last_wgt = ONE - this % absorb_wgt = ZERO - this % n_bank = 0 - this % wgt_bank = ZERO - this % sqrtkT = ERROR_REAL - this % n_collision = 0 - this % fission = .false. - this % delayed_group = 0 - this % n_delayed_bank(:) = 0 - this % g = NONE - - ! Set up base level coordinates - this % coord(1) % universe = NONE - this % n_coord = 1 - this % last_n_coord = 1 - - end subroutine initialize - -!=============================================================================== -! INITIALIZE_FROM_SOURCE initializes a particle from data stored in a source -! site. The source site may have been produced from an external source, from -! fission, or simply as a secondary particle. -!=============================================================================== - - subroutine initialize_from_source(this, src, run_CE, energy_bin_avg) - class(Particle), intent(inout) :: this - type(Bank), intent(in) :: src - logical, intent(in) :: run_CE - real(8), allocatable, intent(in) :: energy_bin_avg(:) - - ! set defaults - call this % initialize() - - ! copy attributes from source bank site - this % type = src % particle - this % wgt = src % wgt - this % last_wgt = src % wgt - this % coord(1) % xyz = src % xyz - this % coord(1) % uvw = src % uvw - this % last_xyz_current = src % xyz - this % last_xyz = src % xyz - this % last_uvw = src % uvw - if (run_CE) then - this % E = src % E - this % g = NONE - else - this % g = int(src % E) - this % last_g = int(src % E) - this % E = energy_bin_avg(this % g) - end if - this % last_E = this % E - - end subroutine initialize_from_source - -!=============================================================================== -! MARK_AS_LOST -!=============================================================================== - - subroutine mark_as_lost(this, message) - class(Particle), intent(inout) :: this - character(*) :: message - - integer(8) :: tot_n_particles - - ! Print warning and write lost particle file - call warning(message) - call this % write_restart() - - ! Increment number of lost particles - this % alive = .false. -!$omp atomic - n_lost_particles = n_lost_particles + 1 - - ! Count the total number of simulated particles (on this processor) - tot_n_particles = current_batch * gen_per_batch * work - - ! Abort the simulation if the maximum number of lost particles has been - ! reached - if (n_lost_particles >= MAX_LOST_PARTICLES .and. & - n_lost_particles >= REL_MAX_LOST_PARTICLES * tot_n_particles) then - call fatal_error("Maximum number of lost particles has been reached.") - end if - - end subroutine mark_as_lost - -!=============================================================================== -! WRITE_RESTART creates a particle restart file -!=============================================================================== - - subroutine write_restart(this) - class(Particle), intent(in) :: this - - integer(HID_T) :: file_id - character(MAX_FILE_LEN) :: filename - - ! Dont write another restart file if in particle restart mode - if (run_mode == MODE_PARTICLE) return - - ! Set up file name - filename = trim(path_output) // 'particle_' // trim(to_str(current_batch)) & - // '_' // trim(to_str(this % id)) // '.h5' - -!$omp critical (WriteParticleRestart) - ! Create file - file_id = file_open(filename, 'w') - - associate (src => source_bank(current_work)) - ! Write filetype and version info - call write_attribute(file_id, 'filetype', 'particle restart') - call write_attribute(file_id, 'version', VERSION_PARTICLE_RESTART) - call write_attribute(file_id, "openmc_version", VERSION) -#ifdef GIT_SHA1 - call write_attribute(file_id, "git_sha1", GIT_SHA1) -#endif - - ! Write data to file - call write_dataset(file_id, 'current_batch', current_batch) - call write_dataset(file_id, 'generations_per_batch', gen_per_batch) - call write_dataset(file_id, 'current_generation', current_gen) - call write_dataset(file_id, 'n_particles', n_particles) - select case(run_mode) - case (MODE_FIXEDSOURCE) - call write_dataset(file_id, 'run_mode', 'fixed source') - case (MODE_EIGENVALUE) - call write_dataset(file_id, 'run_mode', 'eigenvalue') - case (MODE_PARTICLE) - call write_dataset(file_id, 'run_mode', 'particle restart') - end select - call write_dataset(file_id, 'id', this % id) - call write_dataset(file_id, 'type', this % type) - call write_dataset(file_id, 'weight', src % wgt) - call write_dataset(file_id, 'energy', src % E) - call write_dataset(file_id, 'xyz', src % xyz) - call write_dataset(file_id, 'uvw', src % uvw) - end associate - - ! Close file - call file_close(file_id) -!$omp end critical (WriteParticleRestart) - - end subroutine write_restart + call particle_mark_as_lost_c(this, to_c_string(message)) + end subroutine particle_mark_as_lost end module particle_header diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 2ded5e75e..0190c6e97 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -9,7 +9,7 @@ module particle_restart use mgxs_interface, only: energy_bin_avg use nuclide_header, only: micro_xs, n_nuclides use output, only: print_particle - use particle_header, only: Particle + use particle_header use random_lcg, only: set_particle_seed use settings use simulation_header @@ -41,7 +41,7 @@ contains allocate(micro_xs(n_nuclides)) ! Initialize the particle to be tracked - call p % initialize() + call particle_initialize(p) ! Read in the restart information call read_particle_restart(p, previous_run_mode) diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 3a84a9a01..e17e0f47d 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -2,7 +2,7 @@ module photon_physics use algorithm, only: binary_search use constants - use particle_header, only: Particle + use particle_header use photon_header, only: PhotonInteraction, BremsstrahlungData, & compton_profile_pz, ttb_e_grid, ttb use random_lcg, only: prn @@ -332,7 +332,7 @@ contains uvw(2) = sqrt(ONE - mu*mu)*cos(phi) uvw(3) = sqrt(ONE - mu*mu)*sin(phi) E = elm % shells(i_shell) % binding_energy - call p % create_secondary(uvw, E, PHOTON, run_ce=.true.) + call particle_create_secondary(p, uvw, E, PHOTON, run_ce=.true._C_BOOL) return end if @@ -363,7 +363,7 @@ contains ! Non-radiative transition -- Auger/Coster-Kronig effect ! Create auger electron - call p % create_secondary(uvw, E, ELECTRON, run_ce=.true.) + call particle_create_secondary(p, uvw, E, ELECTRON, run_ce=.true._C_BOOL) ! Fill hole left by emitted auger electron i_hole = elm % shell_dict % get(secondary) @@ -372,7 +372,7 @@ contains ! Radiative transition -- get X-ray energy ! Create fluorescent photon - call p % create_secondary(uvw, E, PHOTON, run_ce=.true.) + call particle_create_secondary(p, uvw, E, PHOTON, run_ce=.true._C_BOOL) end if @@ -630,7 +630,8 @@ contains if (w > energy_cutoff(PHOTON)) then ! Create secondary photon - call p % create_secondary(p % coord(1) % uvw, w, PHOTON, run_ce=.true.) + call particle_create_secondary(p, p % coord(1) % uvw, w, PHOTON, & + run_ce=.true._C_BOOL) E_lost = E_lost + w end if end do diff --git a/src/physics.F90 b/src/physics.F90 index c0a813ba2..0eed83418 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -9,7 +9,7 @@ module physics use mesh_header, only: meshes use message_passing use nuclide_header - use particle_header, only: Particle + use particle_header use photon_header use photon_physics, only: rayleigh_scatter, compton_scatter, & atomic_relaxation, pair_production, & @@ -231,7 +231,7 @@ contains / sqrt(alpha**2 + alpha_out**2 - TWO*alpha*alpha_out*mu) phi = TWO*PI*prn() uvw = rotate_angle(p % coord(1) % uvw, mu_electron, phi) - call p % create_secondary(uvw, E_electron, ELECTRON, .true.) + call particle_create_secondary(p, uvw, E_electron, ELECTRON, .true._C_BOOL) ! TODO: Compton subshell data does not match atomic relaxation data ! Allow electrons to fill orbital and produce auger electrons @@ -288,7 +288,8 @@ contains uvw(3) = sqrt(ONE - mu*mu)*sin(phi) ! Create secondary electron - call p % create_secondary(uvw, E_electron, ELECTRON, run_CE=.true.) + call particle_create_secondary(p, uvw, E_electron, ELECTRON, & + run_CE=.true._C_BOOL) ! Allow electrons to fill orbital and produce auger electrons ! and fluorescent photons @@ -311,11 +312,11 @@ contains ! Create secondary electron uvw = rotate_angle(p % coord(1) % uvw, mu_electron) - call p % create_secondary(uvw, E_electron, ELECTRON, .true.) + call particle_create_secondary(p, uvw, E_electron, ELECTRON, .true._C_BOOL) ! Create secondary positron uvw = rotate_angle(p % coord(1) % uvw, mu_positron) - call p % create_secondary(uvw, E_positron, POSITRON, .true.) + call particle_create_secondary(p, uvw, E_positron, POSITRON, .true._C_BOOL) p % event_MT = PAIR_PROD p % alive = .false. @@ -380,8 +381,8 @@ contains uvw(3) = sqrt(ONE - mu*mu)*sin(phi) ! Create annihilation photon pair traveling in opposite directions - call p % create_secondary( uvw, MASS_ELECTRON_EV, PHOTON, .true.) - call p % create_secondary(-uvw, MASS_ELECTRON_EV, PHOTON, .true.) + call particle_create_secondary(p, uvw, MASS_ELECTRON_EV, PHOTON, .true._C_BOOL) + call particle_create_secondary(p, -uvw, MASS_ELECTRON_EV, PHOTON, .true._C_BOOL) p % E = ZERO p % alive = .false. @@ -425,7 +426,7 @@ contains ! Check to make sure that a nuclide was sampled if (i_nuc_mat > mat % n_nuclides) then - call p % write_restart() + call particle_write_restart(p) call fatal_error("Did not sample any nuclide during collision.") end if @@ -475,7 +476,7 @@ contains ! Check to make sure that a nuclide was sampled if (i > mat % n_nuclides) then - call p % write_restart() + call particle_write_restart(p) call fatal_error("Did not sample any element during collision.") end if @@ -744,7 +745,7 @@ contains ! Check to make sure inelastic scattering reaction sampled if (i > size(nuc % reactions)) then - call p % write_restart() + call particle_write_restart(p) call fatal_error("Did not sample any reaction for nuclide " & &// trim(nuc % name)) end if @@ -1451,7 +1452,7 @@ contains ! Determine indices on ufs mesh for current location call m % get_bin(p % coord(1) % xyz, mesh_bin) if (mesh_bin == NO_BIN_FOUND) then - call p % write_restart() + call particle_write_restart(p) call fatal_error("Source site outside UFS mesh!") end if @@ -1607,7 +1608,7 @@ contains ! check for large number of resamples n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then - ! call p % write_restart() + ! call particle_write_restart(p) call fatal_error("Resampled energy distribution maximum number of " & // "times for nuclide " // nuc % name) end if @@ -1631,7 +1632,7 @@ contains ! check for large number of resamples n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then - ! call p % write_restart() + ! call particle_write_restart(p) call fatal_error("Resampled energy distribution maximum number of " & // "times for nuclide " // nuc % name) end if @@ -1695,8 +1696,8 @@ contains if (mod(yield, ONE) == ZERO) then ! If yield is integral, create exactly that many secondary particles do i = 1, nint(yield) - 1 - call p % create_secondary(p % coord(1) % uvw, p % E, & - NEUTRON, run_CE=.true.) + call particle_create_secondary(p, p % coord(1) % uvw, p % E, & + NEUTRON, run_CE=.true._C_BOOL) end do else ! Otherwise, change weight of particle based on yield @@ -1746,7 +1747,7 @@ contains uvw = rotate_angle(p % coord(1) % uvw, mu) ! Create the secondary photon - call p % create_secondary(uvw, E, PHOTON, run_CE=.true.) + call particle_create_secondary(p, uvw, E, PHOTON, run_CE=.true._C_BOOL) end do end subroutine sample_secondary_photons diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index fd7f7d9a9..40df6e1cb 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -11,7 +11,7 @@ module physics_mg use mgxs_interface use message_passing use nuclide_header, only: material_xs - use particle_header, only: Particle + use particle_header use physics_common use random_lcg, only: prn use settings @@ -185,7 +185,7 @@ contains call m % get_bin(p % coord(1) % xyz, mesh_bin) if (mesh_bin == NO_BIN_FOUND) then - call p % write_restart() + call particle_write_restart(p) call fatal_error("Source site outside UFS mesh!") end if diff --git a/src/plot.F90 b/src/plot.F90 index 789ad0ee9..d8db97ec6 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -9,7 +9,7 @@ module plot use hdf5_interface use output, only: time_stamp use material_header, only: materials - use particle_header, only: LocalCoord, Particle + use particle_header use plot_header use progress_header, only: ProgressBar use settings, only: check_overlaps @@ -161,7 +161,7 @@ contains end if ! allocate and initialize particle - call p % initialize() + call particle_initialize(p) p % coord(1) % xyz = xyz p % coord(1) % uvw = [ HALF, HALF, HALF ] p % coord(1) % universe = root_universe @@ -388,7 +388,7 @@ contains ll = pl % origin - pl % width / TWO ! allocate and initialize particle - call p % initialize() + call particle_initialize(p) p % coord(1) % xyz = ll p % coord(1) % uvw = [ HALF, HALF, HALF ] p % coord(1) % universe = root_universe diff --git a/src/settings.F90 b/src/settings.F90 index 1fbcb9209..972eba1d5 100644 --- a/src/settings.F90 +++ b/src/settings.F90 @@ -9,7 +9,7 @@ module settings ! ============================================================================ ! ENERGY TREATMENT RELATED VARIABLES - logical :: run_CE = .true. ! Run in CE mode? + logical(C_BOOL) :: run_CE = .true. ! Run in CE mode? ! ============================================================================ ! CONTINUOUS-ENERGY CROSS SECTION RELATED VARIABLES diff --git a/src/simulation.F90 b/src/simulation.F90 index e721036f2..c8c9dfcf8 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -25,7 +25,7 @@ module simulation use output, only: header, print_columns, & print_batch_keff, print_generation, print_runtime, & print_results, print_overlap_check, write_tallies - use particle_header, only: Particle + use particle_header use photon_header, only: micro_photon_xs, n_elements use random_lcg, only: set_particle_seed use settings @@ -142,7 +142,7 @@ contains integer :: i ! set defaults - call p % initialize_from_source(source_bank(index_source), run_CE, & + call particle_from_source(p, source_bank(index_source), run_CE, & energy_bin_avg) ! set identifier for particle diff --git a/src/simulation.h b/src/simulation.h new file mode 100644 index 000000000..b7544a6da --- /dev/null +++ b/src/simulation.h @@ -0,0 +1,13 @@ +#ifndef SIMULATION_H +#define SIMULATION_H + +#include + +extern "C" int openmc_current_batch; +extern "C" int openmc_current_gen; +extern "C" int64_t openmc_current_work; +extern "C" int openmc_n_lost_particles; + +#pragma omp threadprivate(openmc_current_work) + +#endif // SIMULATION_H diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 806bcc65e..60bed6425 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -13,15 +13,15 @@ module simulation_header ! GEOMETRY-RELATED VARIABLES ! Number of lost particles - integer :: n_lost_particles = 0 + integer(C_INT), bind(C, name='openmc_n_lost_particles') :: n_lost_particles = 0 real(8) :: log_spacing ! spacing on logarithmic grid ! ============================================================================ ! SIMULATION VARIABLES - integer :: current_batch ! current batch - integer :: current_gen ! current generation within a batch + integer(C_INT), bind(C, name='openmc_current_batch') :: current_batch ! current batch + integer(C_INT), bind(C, name='openmc_current_gen') :: current_gen ! current generation within a batch integer :: total_gen = 0 ! total number of generations simulated logical(C_BOOL), bind(C, name='openmc_simulation_initialized') :: & simulation_initialized = .false. @@ -34,7 +34,7 @@ module simulation_header integer(C_INT64_T), bind(C, name='openmc_work') :: work ! number of particles per processor integer(C_INT64_T), allocatable :: work_index(:) ! starting index in source bank for each process - integer(8) :: current_work ! index in source bank of current history simulated + integer(C_INT64_T), bind(C, name='openmc_current_work') :: current_work ! index in source bank of current history simulated ! ============================================================================ ! K-EIGENVALUE SIMULATION VARIABLES diff --git a/src/source_header.F90 b/src/source_header.F90 index 3cac26aa7..163ae1c19 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -10,7 +10,7 @@ module source_header use geometry, only: find_cell use material_header, only: materials use nuclide_header, only: energy_min, energy_max - use particle_header, only: Particle + use particle_header use settings, only: photon_transport use string, only: to_lower use xml_interface @@ -247,7 +247,7 @@ contains found = .false. do while (.not. found) ! Set particle defaults - call p % initialize() + call particle_initialize(p) ! Set particle type site % particle = this % particle @@ -290,7 +290,7 @@ contains ! Increment number of accepted samples n_accept = n_accept + 1 - call p % clear() + call particle_clear(p) ! Sample angle site % uvw(:) = this % angle % sample() diff --git a/src/surface.cpp b/src/surface.cpp index 56e57a854..6400da2f7 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -318,7 +318,7 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-plane", false); std::array coeffs {{x0}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -383,7 +383,7 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-plane", false); std::array coeffs {{y0}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -449,7 +449,7 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-plane", false); std::array coeffs {{z0}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -510,7 +510,7 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "plane", false); std::array coeffs {{A, B, C, D}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -642,7 +642,7 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cylinder", false); std::array coeffs {{y0, z0, r}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } //============================================================================== @@ -676,7 +676,7 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-cylinder", false); std::array coeffs {{x0, z0, r}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } //============================================================================== @@ -710,7 +710,7 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-cylinder", false); std::array coeffs {{x0, y0, r}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } //============================================================================== @@ -781,7 +781,7 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "sphere", false); std::array coeffs {{x0, y0, z0, r}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } //============================================================================== @@ -898,7 +898,7 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cone", false); std::array coeffs {{x0, y0, z0, r_sq}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } //============================================================================== @@ -932,7 +932,7 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-cone", false); std::array coeffs {{x0, y0, z0, r_sq}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } //============================================================================== @@ -966,7 +966,7 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-cone", false); std::array coeffs {{x0, y0, z0, r_sq}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } //============================================================================== @@ -1060,7 +1060,7 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "quadric", false); std::array coeffs {{A, B, C, D, E, F, G, H, J, K}}; - write_double(group_id, "coefficients", coeffs, false); + write_dataset(group_id, "coefficients", coeffs); } //============================================================================== diff --git a/src/tracking.F90 b/src/tracking.F90 index f443d71f8..dd0c26b43 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -11,7 +11,7 @@ module tracking use message_passing use mgxs_interface use nuclide_header - use particle_header, only: LocalCoord, Particle + use particle_header use physics, only: collision use physics_mg, only: collision_mg use random_lcg, only: prn, prn_set_stream @@ -94,7 +94,7 @@ contains if (p % coord(p % n_coord) % cell == NONE) then call find_cell(p, found_cell) if (.not. found_cell) then - call p % mark_as_lost("Could not find the cell containing" & + call particle_mark_as_lost(p, "Could not find the cell containing" & // " particle " // trim(to_str(p %id))) return end if @@ -274,8 +274,8 @@ contains ! Check for secondary particles if this particle is dead if (.not. p % alive) then if (p % n_secondary > 0) then - call p % initialize_from_source(p % secondary_bank(p % n_secondary), & - run_CE, energy_bin_avg) + call particle_from_source(p, p % secondary_bank(p % n_secondary), & + run_CE, energy_bin_avg) p % n_secondary = p % n_secondary - 1 n_event = 0 @@ -355,7 +355,7 @@ contains ! Do not handle reflective boundary conditions on lower universes if (p % n_coord /= 1) then - call p % mark_as_lost("Cannot reflect particle " & + call particle_mark_as_lost(p, "Cannot reflect particle " & // trim(to_str(p % id)) // " off surface in a lower universe.") return end if @@ -392,7 +392,7 @@ contains p % n_coord = 1 call find_cell(p, found) if (.not. found) then - call p % mark_as_lost("Couldn't find particle after reflecting& + call particle_mark_as_lost(p, "Couldn't find particle after reflecting& & from surface " // trim(to_str(surf % id())) // ".") return end if @@ -412,7 +412,7 @@ contains ! Do not handle periodic boundary conditions on lower universes if (p % n_coord /= 1) then - call p % mark_as_lost("Cannot transfer particle " & + call particle_mark_as_lost(p, "Cannot transfer particle " & // trim(to_str(p % id)) // " across surface in a lower universe.& & Boundary conditions must be applied to universe 0.") return @@ -447,7 +447,7 @@ contains p % n_coord = 1 call find_cell(p, found) if (.not. found) then - call p % mark_as_lost("Couldn't find particle after hitting & + call particle_mark_as_lost(p, "Couldn't find particle after hitting & &periodic boundary on surface " // trim(to_str(surf % id())) & // ".") return @@ -505,7 +505,7 @@ contains ! undefined region in the geometry. if (.not. found) then - call p % mark_as_lost("After particle " // trim(to_str(p % id)) & + call particle_mark_as_lost(p, "After particle " // trim(to_str(p % id)) & // " crossed surface " // trim(to_str(surf % id())) & // " it could not be located in any cell and it did not leak.") return diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index f793d4689..93087df89 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -16,7 +16,7 @@ module volume_calc use material_header, only: materials use message_passing use nuclide_header, only: nuclides - use particle_header, only: Particle + 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 @@ -159,7 +159,7 @@ contains i_end = i_start + min_samples - 1 end if - call p % initialize() + call particle_initialize(p) !$omp parallel private(i, j, k, i_domain, i_material, level, found_cell, & !$omp& indices, hits, n_mat) firstprivate(p)