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

Convert Particle type to C++, add templated HDF5 functions

Browse source
This commit is contained in:
Paul Romano 2018-07-06 06:36:05 -05:00
parent 6ea533a61f
commit 764f638d66
27 changed files with 734 additions and 465 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

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

2
include/openmc.h
View file

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

4
src/api.F90
View file

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

3
src/cell.cpp
View file

@ -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()) {

12
src/geometry.F90
View file

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

6
src/geometry.h Normal file
View file

@ -0,0 +1,6 @@
#ifndef GEOMETRY_H
#define GEOMETRY_H
extern "C" int openmc_root_universe;
#endif // GEOMETRY_H

2
src/geometry_header.F90
View file

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

24
src/hdf5_interface.cpp
View file

@ -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<int>::type_id = H5T_NATIVE_INT;
template<>
const hid_t H5TypeMap<int64_t>::type_id = H5T_NATIVE_INT64;
template<>
const hid_t H5TypeMap<double>::type_id = H5T_NATIVE_DOUBLE;
template <>
const hid_t H5TypeMap<char>::type_id = H5T_NATIVE_CHAR;
} // namespace openmc

183
src/hdf5_interface.h
View file

@ -5,6 +5,7 @@
#include "hdf5_hl.h"
#include <array>
#include <cstddef>
#include <string>
#include <sstream>
#include <vector>
@ -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<double>& result,
@ -89,51 +70,115 @@ read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > >& 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<std::size_t array_len> void
write_int(hid_t group_id, char const *name,
const std::array<int, array_len> &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<typename T>
struct H5TypeMap { static const hid_t type_id; };
//==============================================================================
// Template functions used to provide simple interface to lower-level functions
//==============================================================================
template<typename T>
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<T>::type_id, &buffer);
}
template<std::size_t array_len> void
write_double(hid_t group_id, char const *name,
const std::array<double, array_len> &buffer, bool indep)
template<> inline void
write_attribute<const char*>(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<typename T, std::size_t N> inline void
write_attribute(hid_t obj_id, const char* name, const std::array<T, N>& buffer)
{
hsize_t dims[] {N};
write_attr(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data());
}
template<typename T> inline void
write_dataset(hid_t obj_id, const char* name, T buffer)
{
write_dataset(obj_id, 0, nullptr, name, H5TypeMap<T>::type_id, &buffer, false);
}
template<> inline void
write_dataset<const char*>(hid_t obj_id, const char* name, const char* buffer)
{
write_string(obj_id, name, buffer, false);
}
template<typename T, std::size_t N> inline void
write_dataset(hid_t obj_id, const char* name, const std::array<T, N>& buffer)
{
hsize_t dims[] {N};
write_dataset(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data(), false);
}
} // namespace openmc

14
src/input_xml.F90
View file

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

28
src/lattice.cpp
View file

@ -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<double, 2> pitch_short {{pitch[0], pitch[1]}};
write_double(lat_group, "pitch", pitch_short, false);
write_dataset(lat_group, "pitch", pitch_short);
std::array<double, 2> 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<int, 2> 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<double, 1> pitch_short {{pitch[0]}};
write_double(lat_group, "pitch", pitch_short, false);
write_dataset(lat_group, "pitch", pitch_short);
std::array<double, 2> 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.

225
src/particle.cpp Normal file
View file

@ -0,0 +1,225 @@
#include "particle.h"
#include <algorithm>
#include <sstream>
#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; i<MAX_COORD; ++i) coord[i].reset();
}
void
Particle::create_secondary(const double* uvw, double E, int type, bool run_CE)
{
if (n_secondary == MAX_SECONDARY) {
fatal_error("Too many secondary particles created.");
}
int64_t n = n_secondary;
secondary_bank[n].particle = type;
secondary_bank[n].wgt = wgt;
std::copy(coord[0].xyz, coord[0].xyz + 3, secondary_bank[n].xyz);
std::copy(uvw, uvw + 3, secondary_bank[n].uvw);
secondary_bank[n].E = E;
if (!run_CE) secondary_bank[n].E = g;
n_secondary += 1;
}
void
Particle::initialize()
{
// Clear coordinate lists
clear();
// Set particle to neutron that's alive
type = NEUTRON;
alive = true;
// clear attributes
surface = 0;
cell_born = 0;
material = 0;
last_material = 0;
last_sqrtkT = 0;
wgt = 1.0;
last_wgt = 1.0;
absorb_wgt = 0.0;
n_bank = 0;
wgt_bank = 0.0;
sqrtkT = -1.0;
n_collision = 0;
fission = false;
delayed_group = 0;
for (int i=0; i<MAX_DELAYED_GROUPS; ++i) {
n_delayed_bank[i] = 0;
}
g = 0;
// Set up base level coordinates
coord[0].universe = C_NONE;
n_coord = 1;
last_n_coord = 1;
}
void
Particle::from_source(const Bank* src, bool run_CE, const double* energy_bin_avg)
{
// set defaults
initialize();
// copy attributes from source bank site
type = src->particle;
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<int>(src->E);
last_g = static_cast<int>(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

160
src/particle.h Normal file
View file

@ -0,0 +1,160 @@
#ifndef PARTICLE_H
#define PARTICLE_H
#include <cstdint>
#include <array>
#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

398
src/particle_header.F90
View file

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

4
src/particle_restart.F90
View file

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

11
src/photon_physics.F90
View file

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

33
src/physics.F90
View file

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

4
src/physics_mg.F90
View file

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

6
src/plot.F90
View file

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

2
src/settings.F90
View file

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

4
src/simulation.F90
View file

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

13
src/simulation.h Normal file
View file

@ -0,0 +1,13 @@
#ifndef SIMULATION_H
#define SIMULATION_H
#include <cstdint>
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

8
src/simulation_header.F90
View file

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

6
src/source_header.F90
View file

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

24
src/surface.cpp
View file

@ -318,7 +318,7 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-plane", false);
std::array<double, 1> 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<double, 1> 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<double, 1> 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<double, 4> 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<double, 3> 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<double, 3> 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<double, 3> 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<double, 4> 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<double, 4> 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<double, 4> 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<double, 4> 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<double, 10> coeffs {{A, B, C, D, E, F, G, H, J, K}};
write_double(group_id, "coefficients", coeffs, false);
write_dataset(group_id, "coefficients", coeffs);
}
//==============================================================================

18
src/tracking.F90
View file

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

4
src/volume_calc.F90
View file

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