Convert particle restart to C++

This commit is contained in:
Paul Romano 2019-02-14 12:43:05 -06:00
parent 24c84e394c
commit 968e9a95d4
10 changed files with 137 additions and 162 deletions

View file

@ -319,10 +319,8 @@ add_library(libopenmc SHARED
src/mgxs_interface.F90
src/nuclide_header.F90
src/particle_header.F90
src/particle_restart.F90
src/photon_header.F90
src/pugixml/pugixml_f.F90
src/random_lcg.F90
src/relaxng
src/settings.F90
src/simulation_header.F90
@ -370,6 +368,7 @@ add_library(libopenmc SHARED
src/nuclide.cpp
src/output.cpp
src/particle.cpp
src/particle_restart.cpp
src/photon.cpp
src/physics.cpp
src/physics_common.cpp

View file

@ -84,7 +84,6 @@ extern "C" {
int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_next_batch(int* status);
int openmc_nuclide_name(int index, const char** name);
int openmc_particle_restart();
int openmc_plot_geometry();
int openmc_reset();
int openmc_run();

View file

@ -223,8 +223,9 @@ read_attribute(hid_t obj_id, const char* name, std::vector<std::string>& vec)
// Templates/overloads for read_dataset and related methods
//==============================================================================
template<typename T>
void read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)
template<typename T> inline
std::enable_if_t<std::is_scalar<std::decay_t<T>>::value>
read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)
{
read_dataset(obj_id, name, H5TypeMap<T>::type_id, &buffer, indep);
}
@ -242,6 +243,14 @@ read_dataset(hid_t obj_id, const char* name, std::string& str, bool indep=false)
str = std::string{buffer, n};
}
// array version
template<typename T, std::size_t N> inline void
read_dataset(hid_t dset, const char* name, std::array<T, N>& buffer, bool indep=false)
{
read_dataset(dset, name, H5TypeMap<T>::type_id, buffer.data(), indep);
}
// vector version
template <typename T>
void read_dataset(hid_t dset, std::vector<T>& vec, bool indep=false)
{

View file

@ -0,0 +1,10 @@
#ifndef OPENMC_PARTICLE_RESTART_H
#define OPENMC_PARTICLE_RESTART_H
namespace openmc {
void run_particle_restart();
} // namespace openmc
#endif // OPENMC_PARTICLE_RESTART_H

View file

@ -83,6 +83,9 @@ void finalize_generation();
//! Determine overall generation number
extern "C" int overall_generation();
extern "C" void simulation_init_f();
extern "C" void simulation_finalize_f();
#ifdef OPENMC_MPI
void broadcast_results();
#endif

View file

@ -147,18 +147,6 @@ module constants
DIFF_NUCLIDE_DENSITY = 2, &
DIFF_TEMPERATURE = 3
! ============================================================================
! RANDOM NUMBER STREAM CONSTANTS
integer(C_INT), bind(C, name='N_STREAMS') :: N_STREAMS
integer(C_INT), bind(C, name='STREAM_TRACKING') :: STREAM_TRACKING
integer(C_INT), bind(C, name='STREAM_TALLIES') :: STREAM_TALLIES
integer(C_INT), bind(C, name='STREAM_SOURCE') :: STREAM_SOURCE
integer(C_INT), bind(C, name='STREAM_URR_PTABLE') :: STREAM_URR_PTABLE
integer(C_INT), bind(C, name='STREAM_VOLUME') :: STREAM_VOLUME
integer(C_INT), bind(C, name='STREAM_PHOTON') :: STREAM_PHOTON
integer(C_INT64_T), parameter :: DEFAULT_SEED = 1_8
! ============================================================================
! MISCELLANEOUS CONSTANTS

View file

@ -5,6 +5,7 @@
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/message_passing.h"
#include "openmc/particle_restart.h"
#include "openmc/settings.h"
@ -36,7 +37,8 @@ int main(int argc, char* argv[]) {
err = openmc_plot_geometry();
break;
case RUN_MODE_PARTICLE:
if (mpi::master) err = openmc_particle_restart();
if (mpi::master) run_particle_restart();
err = 0;
break;
case RUN_MODE_VOLUME:
err = openmc_calculate_volumes();

View file

@ -1,142 +0,0 @@
module particle_restart
use, intrinsic :: ISO_FORTRAN_ENV
use bank_header, only: Bank
use constants
use error, only: write_message
use hdf5_interface, only: file_open, file_close, read_dataset, HID_T
use mgxs_interface, only: energy_bin_avg
use nuclide_header, only: micro_xs, nuclides_size
use particle_header
use photon_header, only: micro_photon_xs, n_elements
use random_lcg, only: set_particle_seed
use settings
use simulation_header
use tally_header, only: n_tallies
implicit none
private
public :: openmc_particle_restart
contains
!===============================================================================
! OPENMC_PARTICLE_RESTART is the main routine that runs the particle restart
!===============================================================================
function openmc_particle_restart() result(err) bind(C)
integer(C_INT) :: err
integer(8) :: particle_seed
integer :: previous_run_mode
type(Particle) :: p
interface
subroutine set_micro_xs() bind(C)
end subroutine
subroutine print_particle(p) bind(C)
import Particle
type(Particle), intent(in) :: p
end subroutine
end interface
err = 0
! Set verbosity high
verbosity = 10
!$omp parallel
allocate(micro_xs(nuclides_size()))
allocate(micro_photon_xs(n_elements))
!$omp end parallel
call set_micro_xs()
! Initialize the particle to be tracked
call particle_initialize(p)
! Read in the restart information
call read_particle_restart(p, previous_run_mode)
! Set all tallies to 0 for now (just tracking errors)
n_tallies = 0
! Compute random number seed
select case (previous_run_mode)
case (MODE_EIGENVALUE)
particle_seed = (total_gen + overall_generation() - 1)*n_particles + p % id
case (MODE_FIXEDSOURCE)
particle_seed = p % id
end select
call set_particle_seed(particle_seed)
! Transport neutron
! FIXME
!call transport(p)
! Write output if particle made it
call print_particle(p)
deallocate(micro_xs)
end function openmc_particle_restart
!===============================================================================
! READ_PARTICLE_RESTART reads the particle restart file
!===============================================================================
subroutine read_particle_restart(p, previous_run_mode)
type(Particle), intent(inout) :: p
integer, intent(inout) :: previous_run_mode
integer(HID_T) :: file_id
character(MAX_WORD_LEN) :: tempstr
! Write meessage
call write_message("Loading particle restart file " &
// trim(path_particle_restart) // "...", 5)
! Open file
file_id = file_open(path_particle_restart, 'r')
! Read data from file
call read_dataset(current_batch, file_id, 'current_batch')
call read_dataset(gen_per_batch, file_id, 'generations_per_batch')
call read_dataset(current_gen, file_id, 'current_generation')
call read_dataset(n_particles, file_id, 'n_particles')
call read_dataset(tempstr, file_id, 'run_mode')
select case (tempstr)
case ('eigenvalue')
previous_run_mode = MODE_EIGENVALUE
case ('fixed source')
previous_run_mode = MODE_FIXEDSOURCE
end select
call read_dataset(p % id, file_id, 'id')
call read_dataset(p % type, file_id, 'type')
call read_dataset(p % wgt, file_id, 'weight')
call read_dataset(p % E, file_id, 'energy')
call read_dataset(p % coord(1) % xyz, file_id, 'xyz')
call read_dataset(p % coord(1) % uvw, file_id, 'uvw')
! Set energy group and average energy in multi-group mode
if (.not. run_CE) then
p % g = int(p % E)
p % E = energy_bin_avg(p % g)
end if
! Set particle last attributes
p % last_wgt = p % wgt
p % last_xyz_current = p % coord(1)%xyz
p % last_xyz = p % coord(1)%xyz
p % last_uvw = p % coord(1)%uvw
p % last_E = p % E
p % last_g = p % g
! Close hdf5 file
call file_close(file_id)
end subroutine read_particle_restart
end module particle_restart

109
src/particle_restart.cpp Normal file
View file

@ -0,0 +1,109 @@
#include "openmc/particle_restart.h"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/output.h"
#include "openmc/particle.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/tally.h"
#include <algorithm> // for copy
#include <array>
#include <string>
namespace openmc {
void read_particle_restart(Particle& p, int& previous_run_mode)
{
// Write meessage
write_message("Loading particle restart file " +
settings::path_particle_restart + "...", 5);
// Open file
hid_t file_id = file_open(settings::path_particle_restart, 'r');
// Read data from file
read_dataset(file_id, "current_batch", simulation::current_batch);
read_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
read_dataset(file_id, "current_generation", simulation::current_gen);
read_dataset(file_id, "n_particles", settings::n_particles);
std::string mode;
read_dataset(file_id, "run_mode", mode);
if (mode == "eigenvalue") {
previous_run_mode = RUN_MODE_EIGENVALUE;
} else if (mode == "fixed source") {
previous_run_mode = RUN_MODE_FIXEDSOURCE;
}
read_dataset(file_id, "id", p.id);
read_dataset(file_id, "type", p.type);
read_dataset(file_id, "weight", p.wgt);
read_dataset(file_id, "energy", p.E);
std::array<double, 3> x;
read_dataset(file_id, "xyz", x);
std::copy(x.data(), x.data() + 3, p.coord[0].xyz);
read_dataset(file_id, "uvw", x);
std::copy(x.data(), x.data() + 3, p.coord[0].uvw);
// Set energy group and average energy in multi-group mode
if (!settings::run_CE) {
p.g = p.E;
p.E = data::energy_bin_avg[p.g - 1];
}
// Set particle last attributes
p.last_wgt = p.wgt;
std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz_current);
std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz);
std::copy(p.coord[0].uvw, p.coord[0].uvw + 3, p.last_uvw);
p.last_E = p.E;
p.last_g = p.g;
// Close hdf5 file
file_close(file_id);
}
void run_particle_restart()
{
// Set verbosity high
settings::verbosity = 10;
simulation_init_f();
set_micro_xs();
// Initialize the particle to be tracked
Particle p;
p.initialize();
// Read in the restart information
int previous_run_mode;
read_particle_restart(p, previous_run_mode);
// Set all tallies to 0 for now (just tracking errors)
model::tallies.clear();
// Compute random number seed
int64_t particle_seed;
switch (previous_run_mode) {
case RUN_MODE_EIGENVALUE:
particle_seed = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + p.id;
break;
case RUN_MODE_FIXEDSOURCE:
particle_seed = p.id;
break;
}
set_particle_seed(particle_seed);
// Transport neutron
p.transport();
// Write output if particle made it
print_particle(&p);
simulation_finalize_f();
}
} // namespace openmc

View file

@ -30,8 +30,6 @@ namespace openmc {
extern "C" void accumulate_tallies();
extern "C" void allocate_tally_results();
extern "C" void setup_active_tallies();
extern "C" void simulation_init_f();
extern "C" void simulation_finalize_f();
extern "C" void write_tallies();
} // namespace openmc