Move source and fission banks to C++

This commit is contained in:
Paul Romano 2018-11-15 08:13:14 -06:00
parent 887636d398
commit b8187571b3
23 changed files with 385 additions and 299 deletions

View file

@ -381,6 +381,7 @@ add_library(libopenmc SHARED
src/tallies/tally_header.F90
src/tallies/trigger.F90
src/tallies/trigger_header.F90
src/bank.cpp
src/dagmc.cpp
src/cell.cpp
src/cmfd_execute.cpp

47
include/openmc/bank.h Normal file
View file

@ -0,0 +1,47 @@
#ifndef OPENMC_BANK_H
#define OPENMC_BANK_H
#include <cstdint>
#include <vector>
namespace openmc {
struct Bank {
double wgt;
double xyz[3];
double uvw[3];
double E;
int delayed_group;
int particle;
};
} // namespace openmc
// Without an explicit instantiation of vector<Bank>, the Intel compiler
// will complain about the threadprivate directive on filter_matches. Note that
// this has to happen *outside* of the openmc namespace
extern template class std::vector<openmc::Bank>;
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
namespace simulation {
extern "C" int64_t n_bank;
extern std::vector<Bank> source_bank;
extern std::vector<Bank> fission_bank;
#ifdef _OPENMP
extern std::vector<Bank> master_fission_bank;
#endif
#pragma omp threadprivate(fission_bank, n_bank)
} // namespace simulation
} // namespace openmc
#endif // OPENMC_BANK_H

View file

@ -6,9 +6,11 @@
#include <stddef.h>
#ifdef __cplusplus
#include "openmc/bank.h"
extern "C" {
#endif
int openmc_fission_bank(openmc::Bank** ptr, int64_t* n);
int openmc_source_bank(openmc::Bank** ptr, int64_t* n);
#else
struct Bank {
double wgt;
double xyz[3];
@ -18,6 +20,10 @@ extern "C" {
int particle;
};
int openmc_fission_bank(struct Bank** ptr, int64_t* n);
int openmc_source_bank(struct Bank** ptr, int64_t* n);
#endif
int openmc_calculate_volumes();
int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n);
int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n);
@ -38,7 +44,6 @@ extern "C" {
int openmc_filter_set_type(int32_t index, const char* type);
int openmc_finalize();
int openmc_find_cell(double* xyz, int32_t* index, int32_t* instance);
int openmc_fission_bank(struct Bank** ptr, int64_t* n);
int openmc_get_cell_index(int32_t id, int32_t* index);
int openmc_get_filter_index(int32_t id, int32_t* index);
void openmc_get_filter_next_id(int32_t* id);
@ -85,7 +90,6 @@ extern "C" {
void openmc_set_seed(int64_t new_seed);
int openmc_simulation_finalize();
int openmc_simulation_init();
int openmc_source_bank(struct Bank** ptr, int64_t* n);
int openmc_spatial_legendre_filter_get_order(int32_t index, int* order);
int openmc_spatial_legendre_filter_get_params(int32_t index, int* axis, double* min, double* max);
int openmc_spatial_legendre_filter_set_order(int32_t index, int order);

View file

@ -25,9 +25,6 @@ extern std::array<double, 2> k_sum; //!< Used to reduce sum and sum_sq
extern std::vector<double> entropy; //!< Shannon entropy at each generation
extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
extern "C" int64_t n_bank;
#pragma omp threadprivate(n_bank)
} // namespace simulation
//==============================================================================
@ -35,15 +32,22 @@ extern "C" int64_t n_bank;
//==============================================================================
//! Collect/normalize the tracklength keff from each process
extern "C" void calculate_generation_keff();
void calculate_generation_keff();
//! Calculate mean/standard deviation of keff during active generations
//!
//! This function sets the global variables keff and keff_std which represent
//! the mean and standard deviation of the mean of k-effective over active
//! generations. It also broadcasts the value from the master process.
extern "C" void calculate_average_keff();
void calculate_average_keff();
#ifdef _OPENMP
//! Join threadprivate fission banks into a single fission bank
//!
//! Note that this operation is necessarily sequential to preserve the order of
//! the bank when using varying numbers of threads.
void join_bank_from_threads();
#endif
//! Calculates a minimum variance estimate of k-effective
//!
@ -60,16 +64,16 @@ extern "C" void calculate_average_keff();
extern "C" int openmc_get_keff(double* k_combined);
//! Sample/redistribute source sites from accumulated fission sites
extern "C" void synchronize_bank();
void synchronize_bank();
//! Calculates the Shannon entropy of the fission source distribution to assess
//! source convergence
extern "C" void shannon_entropy();
void shannon_entropy();
//! Determines the source fraction in each UFS mesh cell and reweights the
//! source bank so that the sum of the weights is equal to n_particles. The
//! 'source_frac' variable is used later to bias the production of fission sites
extern "C" void ufs_count_sites();
void ufs_count_sites();
//! Get UFS weight corresponding to particle's location
extern "C" double ufs_get_weight(const Particle* p);

View file

@ -56,6 +56,9 @@ extern "C" int thread_id; //!< ID of a given thread
// Functions
//==============================================================================
//! Allocate space for source and fission banks
void allocate_banks();
//! Determine number of particles to transport per process
void calculate_work();

View file

@ -9,9 +9,9 @@
#include "pugixml.hpp"
#include "openmc/bank.h"
#include "openmc/distribution_multi.h"
#include "openmc/distribution_spatial.h"
#include "openmc/capi.h"
#include "openmc/particle.h"
namespace openmc {
@ -62,7 +62,7 @@ extern "C" void initialize_source();
//! Sample a site from all external source distributions in proportion to their
//! source strength
//! \return Sampled source site
extern "C" Bank sample_external_source();
Bank sample_external_source();
//! Fill source bank at end of generation for fixed source simulations
void fill_source_bank_fixedsource();

View file

@ -10,8 +10,8 @@
namespace openmc {
void write_source_point(const char* filename);
extern "C" void write_source_bank(hid_t group_id, Bank* source_bank);
extern "C" void read_source_bank(hid_t group_id, Bank* source_bank);
extern "C" void write_source_bank(hid_t group_id);
extern "C" void read_source_bank(hid_t group_id);
extern "C" void write_tally_results_nr(hid_t file_id);
extern "C" void restart_set_keff();

View file

@ -82,6 +82,9 @@ contains
subroutine free_memory_settings() bind(C)
end subroutine free_memory_settings
subroutine free_memory_bank() bind(C)
end subroutine free_memory_bank
end interface
call free_memory_geometry()

107
src/bank.cpp Normal file
View file

@ -0,0 +1,107 @@
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/error.h"
#include <cstdint>
// Explicit template instantiation definition
template class std::vector<openmc::Bank>;
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
namespace simulation {
int64_t n_bank;
std::vector<Bank> source_bank;
std::vector<Bank> fission_bank;
#ifdef _OPENMP
std::vector<Bank> master_fission_bank;
#endif
} // namespace simulation
//==============================================================================
// C API
//==============================================================================
extern "C" int openmc_source_bank(Bank** ptr, int64_t* n)
{
if (simulation::source_bank.size() == 0) {
set_errmsg("Source bank has not been allocated.");
return OPENMC_E_ALLOCATE;
} else {
*ptr = simulation::source_bank.data();
*n = simulation::source_bank.size();
return 0;
}
}
extern "C" int openmc_fission_bank(Bank** ptr, int64_t* n)
{
if (simulation::fission_bank.size() == 0) {
set_errmsg("Fission bank has not been allocated.");
return OPENMC_E_ALLOCATE;
} else {
*ptr = simulation::fission_bank.data();
*n = simulation::fission_bank.size();
return 0;
}
}
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" void free_memory_bank()
{
simulation::source_bank.clear();
#pragma omp parallel
{
simulation::fission_bank.clear();
}
#ifdef _OPENMP
simulation::master_fission_bank.clear();
#endif
}
extern "C" int fission_bank_delayed_group(int64_t i) {
return simulation::fission_bank[i-1].delayed_group;
}
extern "C" double fission_bank_E(int64_t i) {
return simulation::fission_bank[i-1].E;
}
extern "C" double fission_bank_wgt(int64_t i) {
return simulation::fission_bank[i-1].wgt;
}
extern "C" void source_bank_xyz(int64_t i, double* xyz)
{
xyz[0] = simulation::source_bank[i-1].xyz[0];
xyz[1] = simulation::source_bank[i-1].xyz[1];
xyz[2] = simulation::source_bank[i-1].xyz[2];
}
extern "C" double source_bank_E(int64_t i)
{
return simulation::source_bank[i-1].E;
}
extern "C" double source_bank_wgt(int64_t i)
{
return simulation::source_bank[i-1].wgt;
}
extern "C" void source_bank_set_wgt(int64_t i, double wgt)
{
simulation::source_bank[i-1].wgt = wgt;
}
} // namespace openmc

View file

@ -2,8 +2,6 @@ module bank_header
use, intrinsic :: ISO_C_BINDING
use error, only: E_ALLOCATE, set_errmsg
implicit none
!===============================================================================
@ -21,70 +19,58 @@ module bank_header
integer(C_INT) :: particle ! particle type (neutron, photon, etc.)
end type Bank
! Source and fission bank
type(Bank), allocatable, target :: source_bank(:)
type(Bank), allocatable, target :: fission_bank(:)
#ifdef _OPENMP
type(Bank), allocatable, target :: master_fission_bank(:)
#endif
integer(C_INT64_T), bind(C) :: n_bank ! # of sites in fission bank
!$omp threadprivate(n_bank)
!$omp threadprivate(fission_bank, n_bank)
interface
function openmc_fission_bank(ptr, n) result(err) bind(C)
import C_PTR, C_INT64_T, C_INT
type(C_PTR), intent(out) :: ptr
integer(C_INT64_T), intent(out) :: n
integer(C_INT) :: err
end function
contains
function fission_bank_delayed_group(i) result(g) bind(C)
import C_INT64_T, C_INT
integer(C_INT64_T), value :: i
integer(C_INT) :: g
end function
!===============================================================================
! FREE_MEMORY_BANK deallocates global arrays defined in this module
!===============================================================================
function fission_bank_E(i) result(E) bind(C, name='fission_bank_E')
import C_INT64_T, C_DOUBLE
integer(C_INT64_T), value :: i
real(C_DOUBLE) :: E
end function
subroutine free_memory_bank()
function fission_bank_wgt(i) result(wgt) bind(C)
import C_INT64_T, C_DOUBLE
integer(C_INT64_T), value :: i
real(C_DOUBLE) :: wgt
end function
! Deallocate fission and source bank and entropy
!$omp parallel
if (allocated(fission_bank)) deallocate(fission_bank)
!$omp end parallel
#ifdef _OPENMP
if (allocated(master_fission_bank)) deallocate(master_fission_bank)
#endif
if (allocated(source_bank)) deallocate(source_bank)
subroutine source_bank_xyz(i, xyz) bind(C)
import C_INT64_T, C_DOUBLE
integer(C_INT64_T), value :: i
real(C_DOUBLE), intent(in) :: xyz(*)
end subroutine
end subroutine free_memory_bank
function source_bank_E(i) result(E) bind(C, name='source_bank_E')
import C_INT64_T, C_DOUBLE
integer(C_INT64_T), value :: i
real(C_DOUBLE) :: E
end function
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function source_bank_wgt(i) result(wgt) bind(C)
import C_INT64_T, C_DOUBLE
integer(C_INT64_T), value :: i
real(C_DOUBLE) :: wgt
end function
function openmc_source_bank(ptr, n) result(err) bind(C)
! Return a pointer to the source bank
type(C_PTR), intent(out) :: ptr
integer(C_INT64_T), intent(out) :: n
integer(C_INT) :: err
if (.not. allocated(source_bank)) then
err = E_ALLOCATE
call set_errmsg("Source bank has not been allocated.")
else
err = 0
ptr = C_LOC(source_bank)
n = size(source_bank)
end if
end function openmc_source_bank
function openmc_fission_bank(ptr, n) result(err) bind(C)
! Return a pointer to the source bank
type(C_PTR), intent(out) :: ptr
integer(C_INT64_T), intent(out) :: n
integer(C_INT) :: err
if (.not. allocated(fission_bank)) then
err = E_ALLOCATE
call set_errmsg("Fission bank has not been allocated.")
else
err = 0
ptr = C_LOC(fission_bank)
n = size(fission_bank)
end if
end function openmc_fission_bank
subroutine source_bank_set_wgt(i, wgt) bind(C)
import C_INT64_T, C_DOUBLE
integer(C_INT64_T), value :: i
real(C_DOUBLE), value :: wgt
end subroutine
end interface
end module bank_header

View file

@ -218,7 +218,7 @@ contains
subroutine cmfd_reweight(new_weights)
use algorithm, only: binary_search
use bank_header, only: source_bank
use bank_header
use constants, only: ZERO, ONE
use error, only: warning, fatal_error
use message_passing
@ -230,13 +230,14 @@ contains
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer(C_INT) :: ng ! maximum number of energy groups
integer :: i ! iteration counter
integer(8) :: i ! iteration counter
integer :: g ! index for group
integer :: ijk(3) ! spatial bin location
integer :: e_bin ! energy bin of source particle
integer :: mesh_bin ! mesh bin of soruce particle
integer :: n_groups ! number of energy groups
real(8) :: norm ! normalization factor
real(C_DOUBLE) :: xyz(3)
logical(C_BOOL) :: outside ! any source sites outside mesh
logical :: in_mesh ! source site is inside mesh
@ -309,21 +310,22 @@ contains
end if
! begin loop over source bank
do i = 1, int(work,4)
do i = 1, work
! Determine spatial bin
call cmfd_mesh % get_indices(source_bank(i) % xyz, ijk, in_mesh)
call source_bank_xyz(i, xyz)
call cmfd_mesh % get_indices(xyz, ijk, in_mesh)
! Determine energy bin
n_groups = size(cmfd % egrid) - 1
if (source_bank(i) % E < cmfd % egrid(1)) then
if (source_bank_E(i) < cmfd % egrid(1)) then
e_bin = 1
if (master) call warning('Source pt below energy grid')
elseif (source_bank(i) % E > cmfd % egrid(n_groups + 1)) then
elseif (source_bank_E(i) > cmfd % egrid(n_groups + 1)) then
e_bin = n_groups
if (master) call warning('Source pt above energy grid')
else
e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank(i) % E)
e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank_E(i))
end if
! Reverese energy bin (lowest grp is highest energy bin)
@ -335,8 +337,8 @@ contains
end if
! Reweight particle
source_bank(i) % wgt = source_bank(i) % wgt * &
cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3))
call source_bank_set_wgt(i, source_bank_wgt(i) * &
cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3)))
end do
end subroutine cmfd_reweight

View file

@ -5,6 +5,7 @@
#include "xtensor/xarray.hpp"
#include "xtensor/xio.hpp"
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
@ -18,14 +19,10 @@ extern "C" void
cmfd_populate_sourcecounts(int n_energy, const double* energies,
double* source_counts, bool* outside)
{
// Get pointer to source bank
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
// Get source counts in each mesh bin / energy bin
auto& m = model::meshes.at(settings::index_cmfd_mesh);
xt::xarray<double> counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside);
xt::xarray<double> counts = m->count_sites(simulation::work,
simulation::source_bank.data(), n_energy, energies, outside);
// Copy data from the xarray into the source counts array
std::copy(counts.begin(), counts.end(), source_counts);

View file

@ -2,9 +2,6 @@ module eigenvalue
use, intrinsic :: ISO_C_BINDING
use bank_header
use simulation_header
implicit none
interface
@ -15,49 +12,4 @@ module eigenvalue
end function
end interface
contains
#ifdef _OPENMP
!===============================================================================
! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission
! bank that can be sampled. Note that this operation is necessarily sequential
! to preserve the order of the bank when using varying numbers of threads.
!===============================================================================
subroutine join_bank_from_threads() bind(C)
integer(8) :: total ! total number of fission bank sites
integer :: i ! loop index for threads
! Initialize the total number of fission bank sites
total = 0
!$omp parallel
! Copy thread fission bank sites to one shared copy
!$omp do ordered schedule(static)
do i = 1, n_threads
!$omp ordered
master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank)
total = total + n_bank
!$omp end ordered
end do
!$omp end do
! Make sure all threads have made it to this point
!$omp barrier
! Now copy the shared fission bank sites back to the master thread's copy.
if (thread_id == 0) then
n_bank = total
fission_bank(1:n_bank) = master_fission_bank(1:n_bank)
else
n_bank = 0
end if
!$omp end parallel
end subroutine join_bank_from_threads
#endif
end module eigenvalue

View file

@ -5,6 +5,7 @@
#include "xtensor/xtensor.hpp"
#include "xtensor/xview.hpp"
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
@ -69,13 +70,6 @@ void synchronize_bank()
{
simulation::time_bank.start();
// Get pointers to source/fission bank
Bank* source_bank;
Bank* fission_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
openmc_fission_bank(&fission_bank, &n);
// In order to properly understand the fission bank algorithm, you need to
// think of the fission and source bank as being one global array divided
// over multiple processors. At the start, each processor has a random amount
@ -147,14 +141,14 @@ void synchronize_bank()
// and the remaining 100 would be randomly sampled.
if (total < settings::n_particles) {
for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
temp_sites[index_temp] = fission_bank[i];
temp_sites[index_temp] = simulation::fission_bank[i];
++index_temp;
}
}
// Randomly sample sites needed
if (prn() < p_sample) {
temp_sites[index_temp] = fission_bank[i];
temp_sites[index_temp] = simulation::fission_bank[i];
++index_temp;
}
}
@ -195,7 +189,8 @@ void synchronize_bank()
// fission bank
sites_needed = settings::n_particles - finish;
for (int i = 0; i < sites_needed; ++i) {
temp_sites[index_temp] = fission_bank[simulation::n_bank - sites_needed + i];
int i_bank = simulation::n_bank - sites_needed + i;
temp_sites[index_temp] = simulation::fission_bank[i_bank];
++index_temp;
}
}
@ -300,7 +295,8 @@ void synchronize_bank()
MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
#else
std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, source_bank);
std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles,
simulation::source_bank.begin());
#endif
simulation::time_bank_sendrecv.stop();
@ -347,6 +343,46 @@ void calculate_average_keff()
}
}
#ifdef _OPENMP
void join_bank_from_threads()
{
// Initialize the total number of fission bank sites
int64_t total = 0;
#pragma omp parallel
{
// Copy thread fission bank sites to one shared copy
#pragma omp for ordered schedule(static)
for (int i = 0; i < simulation::n_threads; ++i) {
#pragma omp ordered
{
std::copy(
&simulation::fission_bank[0],
&simulation::fission_bank[0] + simulation::n_bank,
&simulation::master_fission_bank[total]
);
total += simulation::n_bank;
}
}
// Make sure all threads have made it to this point
#pragma omp barrier
// Now copy the shared fission bank sites back to the master thread's copy.
if (simulation::thread_id == 0) {
simulation::n_bank = total;
std::copy(
&simulation::master_fission_bank[0],
&simulation::master_fission_bank[0] + simulation::n_bank,
&simulation::fission_bank[0]
);
} else {
simulation::n_bank = 0;
}
}
}
#endif
int openmc_get_keff(double* k_combined)
{
k_combined[0] = 0.0;
@ -499,15 +535,10 @@ void shannon_entropy()
// Get pointer to entropy mesh
auto& m = model::meshes[settings::index_entropy_mesh];
// Get pointer to fission bank
Bank* fission_bank;
int64_t n;
openmc_fission_bank(&fission_bank, &n);
// Get source weight in each mesh bin
bool sites_outside;
xt::xtensor<double, 1> p = m->count_sites(
simulation::n_bank, fission_bank, 0, nullptr, &sites_outside);
xt::xtensor<double, 1> p = m->count_sites(simulation::n_bank,
simulation::fission_bank.data(), 0, nullptr, &sites_outside);
// display warning message if there were sites outside entropy box
if (sites_outside) {
@ -544,15 +575,10 @@ void ufs_count_sites()
s = m->volume_frac_;
} else {
// Get pointer to source bank
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
// count number of source sites in each ufs mesh cell
bool sites_outside;
simulation::source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr,
&sites_outside);
simulation::source_frac = m->count_sites(simulation::work,
simulation::source_bank.data(), 0, nullptr, &sites_outside);
// Check for sites outside of the mesh
if (mpi::master && sites_outside) {
@ -572,7 +598,7 @@ void ufs_count_sites()
// Since the total starting weight is not equal to n_particles, we need to
// renormalize the weight of the source sites
for (int i = 0; i < simulation::work; ++i) {
source_bank[i].wgt *= settings::n_particles / total;
simulation::source_bank[i].wgt *= settings::n_particles / total;
}
}
}

View file

@ -3,6 +3,7 @@
#include <algorithm>
#include <sstream>
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
@ -184,17 +185,12 @@ Particle::write_restart() const
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 = simulation::current_work;
write_dataset(file_id, "weight", src[i-1].wgt);
write_dataset(file_id, "energy", src[i-1].E);
write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt);
write_dataset(file_id, "energy", simulation::source_bank[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);
write_double(file_id, 1, dims, "xyz", simulation::source_bank[i-1].xyz, false);
write_double(file_id, 1, dims, "uvw", simulation::source_bank[i-1].uvw, false);
// Close file
file_close(file_id);

View file

@ -1,6 +1,7 @@
module physics
use algorithm, only: binary_search
use bank_header
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message
@ -86,7 +87,11 @@ contains
integer :: i_nuclide ! index in nuclides array
integer :: i_nuc_mat ! index in material's nuclides array
integer :: i_reaction ! index in nuc % reactions array
integer(C_INT) :: err
integer(C_INT64_T) :: n
type(Nuclide), pointer :: nuc
type(C_PTR) :: ptr
type(Bank), pointer :: fission_bank(:)
call sample_nuclide(p, 'total ', i_nuclide, i_nuc_mat)
@ -103,6 +108,10 @@ contains
if (nuc % fissionable) then
if (run_mode == MODE_EIGENVALUE) then
! Get fission bank pointer
err = openmc_fission_bank(ptr, n)
call c_f_pointer(ptr, fission_bank, [n])
call sample_fission(i_nuclide, p % E, i_reaction)
call create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank)
elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then

View file

@ -5,6 +5,7 @@
#include "xtensor/xarray.hpp"
#include "openmc/bank.h"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
@ -47,12 +48,8 @@ sample_reaction(Particle* p, const MaterialMacroXS* material_xs)
if (model::materials[p->material - 1]->fissionable) {
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
Bank* result_bank;
int64_t result_bank_size;
// Get pointer to fission bank from Fortran side
openmc_fission_bank(&result_bank, &result_bank_size);
create_fission_sites(p, result_bank, &simulation::n_bank, result_bank_size,
material_xs);
create_fission_sites(p, simulation::fission_bank.data(), &simulation::n_bank,
simulation::fission_bank.size(), material_xs);
} else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
(settings::create_fission_neutrons)) {
create_fission_sites(p, p->secondary_bank, &(p->n_secondary),

View file

@ -2,20 +2,9 @@ module simulation
use, intrinsic :: ISO_C_BINDING
#ifdef _OPENMP
use omp_lib
#endif
use bank_header, only: source_bank
use constants, only: ZERO
use error, only: fatal_error
use material_header, only: n_materials, materials
use message_passing
use nuclide_header, only: micro_xs, n_nuclides
use photon_header, only: micro_photon_xs, n_elements
use settings
use simulation_header
use tally_header
use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer
implicit none
@ -69,56 +58,4 @@ contains
end subroutine
!===============================================================================
! ALLOCATE_BANKS allocates memory for the fission and source banks
!===============================================================================
subroutine allocate_banks() bind(C)
integer :: alloc_err ! allocation error code
! Allocate source bank
if (allocated(source_bank)) deallocate(source_bank)
allocate(source_bank(work), STAT=alloc_err)
! Check for allocation errors
if (alloc_err /= 0) then
call fatal_error("Failed to allocate source bank.")
end if
if (run_mode == MODE_EIGENVALUE) then
#ifdef _OPENMP
! If OpenMP is being used, each thread needs its own private fission
! bank. Since the private fission banks need to be combined at the end of
! a generation, there is also a 'master_fission_bank' that is used to
! collect the sites from each thread.
n_threads = omp_get_max_threads()
!$omp parallel
thread_id = omp_get_thread_num()
if (allocated(fission_bank)) deallocate(fission_bank)
if (thread_id == 0) then
allocate(fission_bank(3*work))
else
allocate(fission_bank(3*work/n_threads))
end if
!$omp end parallel
if (allocated(master_fission_bank)) deallocate(master_fission_bank)
allocate(master_fission_bank(3*work), STAT=alloc_err)
#else
if (allocated(fission_bank)) deallocate(fission_bank)
allocate(fission_bank(3*work), STAT=alloc_err)
#endif
! Check for allocation errors
if (alloc_err /= 0) then
call fatal_error("Failed to allocate fission bank.")
end if
end if
end subroutine allocate_banks
end module simulation

View file

@ -1,5 +1,6 @@
#include "openmc/simulation.h"
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/container_util.h"
#include "openmc/eigenvalue.h"
@ -15,6 +16,8 @@
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include <omp.h>
#include <algorithm>
#include <string>
@ -24,14 +27,12 @@ namespace openmc {
extern "C" bool cmfd_on;
extern "C" void accumulate_tallies();
extern "C" void allocate_banks();
extern "C" void allocate_tally_results();
extern "C" void check_triggers();
extern "C" void cmfd_init_batch();
extern "C" void cmfd_tally_init();
extern "C" void execute_cmfd();
extern "C" void init_tally_routines();
extern "C" void join_bank_from_threads();
extern "C" void load_state_point();
extern "C" void print_batch_keff();
extern "C" void print_columns();
@ -280,6 +281,36 @@ int thread_id; //!< ID of a given thread
// Non-member functions
//==============================================================================
void allocate_banks()
{
// Allocate source bank
simulation::source_bank.resize(simulation::work);
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
#ifdef _OPENMP
// If OpenMP is being used, each thread needs its own private fission
// bank. Since the private fission banks need to be combined at the end of
// a generation, there is also a 'master_fission_bank' that is used to
// collect the sites from each thread.
simulation::n_threads = omp_get_max_threads();
#pragma omp parallel
{
simulation::thread_id = omp_get_thread_num();
if (simulation::thread_id == 0) {
simulation::fission_bank.resize(3*simulation::work);
} else {
simulation::fission_bank.resize(3*simulation::work / simulation::n_threads);
}
}
simulation::master_fission_bank.resize(3*simulation::work);
#else
simulation::fission_bank.resize(3*simulation::work);
#endif
}
}
void initialize_batch()
{
// Increment current batch
@ -453,13 +484,8 @@ void finalize_generation()
void initialize_history(Particle* p, int64_t index_source)
{
// Get pointer to source bank
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
// set defaults
p->from_source(&source_bank[index_source - 1]);
p->from_source(&simulation::source_bank[index_source - 1]);
// set identifier for particle
p->id = simulation::work_index[mpi::rank] + index_source;

View file

@ -5,6 +5,7 @@
#include "xtensor/xadapt.hpp"
#include "openmc/bank.h"
#include "openmc/cell.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
@ -242,11 +243,6 @@ void initialize_source()
{
write_message("Initializing source particles...", 5);
// Get pointer to source bank
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
if (settings::path_source != "") {
// Read the source from a binary file instead of sampling from some
// assumed source distribution
@ -268,7 +264,7 @@ void initialize_source()
}
// Read in the source bank
read_source_bank(file_id, source_bank);
read_source_bank(file_id);
// Close file
file_close(file_id);
@ -282,7 +278,7 @@ void initialize_source()
set_particle_seed(id);
// sample external source distribution
source_bank[i] = sample_external_source();
simulation::source_bank[i] = sample_external_source();
}
}
@ -291,7 +287,7 @@ void initialize_source()
write_message("Writing out initial source...", 5);
std::string filename = settings::path_output + "initial_source.h5";
hid_t file_id = file_open(filename, 'w', true);
write_source_bank(file_id, source_bank);
write_source_bank(file_id);
file_close(file_id);
}
}
@ -354,11 +350,6 @@ extern "C" double total_source_strength()
void fill_source_bank_fixedsource()
{
if (settings::path_source.empty()) {
// Get pointer to source bank
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
for (int64_t i = 0; i < simulation::work; ++i) {
// initialize random number seed
int64_t id = (simulation::total_gen + overall_generation()) *
@ -366,7 +357,7 @@ void fill_source_bank_fixedsource()
set_particle_seed(id);
// sample external source distribution
source_bank[i] = sample_external_source();
simulation::source_bank[i] = sample_external_source();
}
}
}

View file

@ -35,16 +35,14 @@ module state_point
implicit none
interface
subroutine write_source_bank(group_id, bank_) bind(C)
import HID_T, C_INT64_T, Bank
subroutine write_source_bank(group_id) bind(C)
import HID_T
integer(HID_T), value :: group_id
type(Bank), intent(in) :: bank_(*)
end subroutine write_source_bank
subroutine read_source_bank(group_id, bank_) bind(C)
import HID_T, C_INT64_T, Bank
subroutine read_source_bank(group_id) bind(C)
import HID_T
integer(HID_T), value :: group_id
type(Bank), intent(out) :: bank_(*)
end subroutine read_source_bank
end interface
@ -437,7 +435,7 @@ contains
if (master .or. parallel) then
file_id = file_open(filename_, 'a', parallel=.true.)
end if
call write_source_bank(file_id, source_bank)
call write_source_bank(file_id)
if (master .or. parallel) call file_close(file_id)
end if
end function openmc_statepoint_write
@ -634,8 +632,8 @@ contains
file_id = file_open(path_source_point, 'r', parallel=.true.)
end if
! Write out source
call read_source_bank(file_id, source_bank)
! Read source
call read_source_bank(file_id)
end if

View file

@ -8,6 +8,7 @@
#include "xtensor/xbuilder.hpp" // for empty_like
#include "xtensor/xview.hpp"
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
@ -70,16 +71,13 @@ write_source_point(const char* filename)
}
// Get pointer to source bank and write to file
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
write_source_bank(file_id, source_bank);
write_source_bank(file_id);
if (mpi::master || parallel) file_close(file_id);
}
void
write_source_bank(hid_t group_id, Bank* source_bank)
write_source_bank(hid_t group_id)
{
hid_t banktype = h5banktype();
@ -103,7 +101,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
// Write data to file in parallel
H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank);
H5Dwrite(dset, banktype, memspace, dspace, plist, simulation::source_bank.data());
// Free resources
H5Sclose(dspace);
@ -122,7 +120,8 @@ write_source_bank(hid_t group_id, Bank* source_bank)
// Save source bank sites since the souce_bank array is overwritten below
#ifdef OPENMC_MPI
std::vector<Bank> temp_source {source_bank, source_bank + simulation::work};
std::vector<Bank> temp_source {simulation::source_bank.begin(),
simulation::source_bank.begin() + simulation::work};
#endif
for (int i = 0; i < mpi::n_procs; ++i) {
@ -134,7 +133,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
#ifdef OPENMC_MPI
// Receive source sites from other processes
if (i > 0)
MPI_Recv(source_bank, count[0], mpi::bank, i, i,
MPI_Recv(source_bank.data(), count[0], mpi::bank, i, i,
mpi::intracomm, MPI_STATUS_IGNORE);
#endif
@ -144,7 +143,8 @@ write_source_bank(hid_t group_id, Bank* source_bank)
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Write data to hyperslab
H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT,
simulation::source_bank.data());
H5Sclose(memspace);
H5Sclose(dspace);
@ -155,12 +155,12 @@ write_source_bank(hid_t group_id, Bank* source_bank)
#ifdef OPENMC_MPI
// Restore state of source bank
std::copy(temp_source.begin(), temp_source.end(), source_bank);
std::copy(temp_source.begin(), temp_source.end(), source_bank.begin());
#endif
} else {
#ifdef OPENMC_MPI
MPI_Send(source_bank, simulation::work, mpi::bank, 0, mpi::rank,
mpi::intracomm);
MPI_Send(simulation::source_bank.data(), simulation::work, mpi::bank,
0, mpi::rank, mpi::intracomm);
#endif
}
#endif
@ -169,7 +169,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
}
void read_source_bank(hid_t group_id, Bank* source_bank)
void read_source_bank(hid_t group_id)
{
hid_t banktype = h5banktype();
@ -197,10 +197,10 @@ void read_source_bank(hid_t group_id, Bank* source_bank)
// Read data in parallel
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
H5Dread(dset, banktype, memspace, dspace, plist, source_bank);
H5Dread(dset, banktype, memspace, dspace, plist, simulation::source_bank.data());
H5Pclose(plist);
#else
H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, simulation::source_bank.data());
#endif
// Close all ids

View file

@ -3,6 +3,7 @@ module tally
use, intrinsic :: ISO_C_BINDING
use algorithm, only: binary_search
use bank_header
use constants
use dict_header, only: EMPTY
use error, only: fatal_error
@ -687,7 +688,7 @@ contains
do k = 1, p % n_bank
! get the delayed group
g = fission_bank(n_bank - p % n_bank + k) % delayed_group
g = fission_bank_delayed_group(n_bank - p % n_bank + k)
! Case for tallying delayed emissions
if (g /= 0) then
@ -697,8 +698,8 @@ contains
reactions(nuclides(p % event_nuclide) % index_fission(1)))
! determine score based on bank site weight and keff.
score = score + keff * fission_bank(n_bank - p % n_bank + k) &
% wgt * rxn % product_decay_rate(1 + g) * flux
score = score + keff * fission_bank_wgt(n_bank - p % n_bank + k) &
* rxn % product_decay_rate(1 + g) * flux
end associate
! if the delayed group filter is present, tally to corresponding
@ -1843,7 +1844,7 @@ contains
do k = 1, p % n_bank
! get the delayed group
g = fission_bank(n_bank - p % n_bank + k) % delayed_group
g = fission_bank_delayed_group(n_bank - p % n_bank + k)
! Case for tallying delayed emissions
if (g /= 0) then
@ -1851,13 +1852,13 @@ contains
! determine score based on bank site weight and keff.
if (i_nuclide > 0) then
score = score + keff * atom_density * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
fission_bank_wgt(n_bank - p % n_bank + k) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
else
score = score + keff * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
fission_bank_wgt(n_bank - p % n_bank + k) * &
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
end if
@ -2375,10 +2376,10 @@ contains
do k = 1, p % n_bank
! get the delayed group
g = fission_bank(n_bank - p % n_bank + k) % delayed_group
g = fission_bank_delayed_group(n_bank - p % n_bank + k)
! determine score based on bank site weight and keff
score = keff * fission_bank(n_bank - p % n_bank + k) % wgt
score = keff * fission_bank_wgt(n_bank - p % n_bank + k)
! Add derivative information for differential tallies. Note that the
! i_nuclide and atom_density arguments do not matter since this is an
@ -2390,7 +2391,7 @@ contains
if (.not. run_CE .and. eo_filt % matches_transport_groups) then
! determine outgoing energy group from fission bank
g_out = int(fission_bank(n_bank - p % n_bank + k) % E)
g_out = int(fission_bank_E(n_bank - p % n_bank + k))
! modify the value so that g_out = 1 corresponds to the highest
! energy bin
@ -2403,10 +2404,9 @@ contains
! determine outgoing energy from fission bank
if (run_CE) then
E_out = fission_bank(n_bank - p % n_bank + k) % E
E_out = fission_bank_E(n_bank - p % n_bank + k)
else
E_out = energy_bin_avg(int(fission_bank(n_bank - p % n_bank + k) &
% E))
E_out = energy_bin_avg(int(fission_bank_E(n_bank - p % n_bank + k)))
end if
! If this outgoing energy falls within the energyout filter's range,