Convert synchronize_bank to C++

This commit is contained in:
Paul Romano 2018-10-11 12:55:33 -05:00
parent 217f87c085
commit 02fbc1582d
10 changed files with 278 additions and 342 deletions

View file

@ -24,6 +24,9 @@ extern "C" int64_t n_bank;
// Non-member functions
//==============================================================================
//! Sample/redistribute source sites from accumulated fission sites
extern "C" void synchronize_bank();
//! Calculates the Shannon entropy of the fission source distribution to assess
//! source convergence
extern "C" void shannon_entropy();

View file

@ -4,7 +4,7 @@
#ifndef OPENMC_SEARCH_H
#define OPENMC_SEARCH_H
#include <algorithm> // for lower_bound
#include <algorithm> // for lower_bound, upper_bound
namespace openmc {
@ -19,6 +19,14 @@ lower_bound_index(It first, It last, const T& value)
return (index == last) ? -1 : index - first;
}
template<class It, class T>
typename std::iterator_traits<It>::difference_type
upper_bound_index(It first, It last, const T& value)
{
It index = std::upper_bound(first, last, value) - 1;
return (index == last) ? -1 : index - first;
}
} // namespace openmc
#endif // OPENMC_SEARCH_H

View file

@ -181,7 +181,6 @@ contains
err = 0
#ifdef OPENMC_MPI
! Free all MPI types
call MPI_TYPE_FREE(MPI_BANK, err)
call openmc_free_bank()
#endif

View file

@ -21,277 +21,6 @@ module eigenvalue
contains
!===============================================================================
! SYNCHRONIZE_BANK samples source sites from the fission sites that were
! accumulated during the generation. This routine is what allows this Monte
! Carlo to scale to large numbers of processors where other codes cannot.
!===============================================================================
subroutine synchronize_bank()
integer :: i ! loop indices
integer :: j ! loop indices
integer(8) :: start ! starting index in global bank
integer(8) :: finish ! ending index in global bank
integer(8) :: total ! total sites in global fission bank
integer(8) :: index_temp ! index in temporary source bank
integer(8) :: sites_needed ! # of sites to be sampled
real(8) :: p_sample ! probability of sampling a site
type(Bank), save, allocatable :: &
& temp_sites(:) ! local array of extra sites on each node
#ifdef OPENMC_MPI
integer :: mpi_err ! MPI error code
integer(8) :: n ! number of sites to send/recv
integer :: neighbor ! processor to send/recv data from
#ifdef OPENMC_MPIF08
type(MPI_Request) :: request(20)
#else
integer :: request(20) ! communication request for send/recving sites
#endif
integer :: n_request ! number of communication requests
integer(8) :: index_local ! index in local source bank
integer(8), save, allocatable :: &
& bank_position(:) ! starting positions in global source bank
#endif
! 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
! of fission bank sites -- each processor needs to know the total number of
! sites in order to figure out the probability for selecting
! sites. Furthermore, each proc also needs to know where in the 'global'
! fission bank its own sites starts in order to ensure reproducibility by
! skipping ahead to the proper seed.
#ifdef OPENMC_MPI
start = 0_8
call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, &
mpi_intracomm, mpi_err)
! While we would expect the value of start on rank 0 to be 0, the MPI
! standard says that the receive buffer on rank 0 is undefined and not
! significant
if (rank == 0) start = 0_8
finish = start + n_bank
total = finish
call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, &
mpi_intracomm, mpi_err)
#else
start = 0_8
finish = n_bank
total = n_bank
#endif
! If there are not that many particles per generation, it's possible that no
! fission sites were created at all on a single processor. Rather than add
! extra logic to treat this circumstance, we really want to ensure the user
! runs enough particles to avoid this in the first place.
if (n_bank == 0) then
call fatal_error("No fission sites banked on processor " // to_str(rank))
end if
! Make sure all processors start at the same point for random sampling. Then
! skip ahead in the sequence using the starting index in the 'global'
! fission bank for each processor.
call set_particle_seed(int(total_gen + overall_generation(), 8))
call advance_prn_seed(start)
! Determine how many fission sites we need to sample from the source bank
! and the probability for selecting a site.
if (total < n_particles) then
sites_needed = mod(n_particles,total)
else
sites_needed = n_particles
end if
p_sample = real(sites_needed,8)/real(total,8)
call time_bank_sample % start()
! ==========================================================================
! SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
! Allocate temporary source bank
index_temp = 0_8
if (.not. allocated(temp_sites)) allocate(temp_sites(3*work))
do i = 1, int(n_bank,4)
! If there are less than n_particles particles banked, automatically add
! int(n_particles/total) sites to temp_sites. For example, if you need
! 1000 and 300 were banked, this would add 3 source sites per banked site
! and the remaining 100 would be randomly sampled.
if (total < n_particles) then
do j = 1, int(n_particles/total)
index_temp = index_temp + 1
temp_sites(index_temp) = fission_bank(i)
end do
end if
! Randomly sample sites needed
if (prn() < p_sample) then
index_temp = index_temp + 1
temp_sites(index_temp) = fission_bank(i)
end if
end do
! At this point, the sampling of source sites is done and now we need to
! figure out where to send source sites. Since it is possible that one
! processor's share of the source bank spans more than just the immediate
! neighboring processors, we have to perform an ALLGATHER to determine the
! indices for all processors
#ifdef OPENMC_MPI
! First do an exclusive scan to get the starting indices for
start = 0_8
call MPI_EXSCAN(index_temp, start, 1, MPI_INTEGER8, MPI_SUM, &
mpi_intracomm, mpi_err)
finish = start + index_temp
! Allocate space for bank_position if this hasn't been done yet
if (.not. allocated(bank_position)) allocate(bank_position(n_procs))
call MPI_ALLGATHER(start, 1, MPI_INTEGER8, bank_position, 1, &
MPI_INTEGER8, mpi_intracomm, mpi_err)
#else
start = 0_8
finish = index_temp
#endif
! Now that the sampling is complete, we need to ensure that we have exactly
! n_particles source sites. The way this is done in a reproducible manner is
! to adjust only the source sites on the last processor.
if (rank == n_procs - 1) then
if (finish > n_particles) then
! If we have extra sites sampled, we will simply discard the extra
! ones on the last processor
index_temp = n_particles - start
elseif (finish < n_particles) then
! If we have too few sites, repeat sites from the very end of the
! fission bank
sites_needed = n_particles - finish
do i = 1, int(sites_needed,4)
index_temp = index_temp + 1
temp_sites(index_temp) = fission_bank(n_bank - sites_needed + i)
end do
end if
! the last processor should not be sending sites to right
finish = work_index(rank + 1)
end if
call time_bank_sample % stop()
call time_bank_sendrecv % start()
#ifdef OPENMC_MPI
! ==========================================================================
! SEND BANK SITES TO NEIGHBORS
index_local = 1
n_request = 0
if (start < n_particles) then
! Determine the index of the processor which has the first part of the
! source_bank for the local processor
neighbor = binary_search(work_index, n_procs + 1, start) - 1
SEND_SITES: do while (start < finish)
! Determine the number of sites to send
n = min(work_index(neighbor + 1), finish) - start
! Initiate an asynchronous send of source sites to the neighboring
! process
if (neighbor /= rank) then
n_request = n_request + 1
call MPI_ISEND(temp_sites(index_local), int(n), MPI_BANK, neighbor, &
rank, mpi_intracomm, request(n_request), mpi_err)
end if
! Increment all indices
start = start + n
index_local = index_local + n
neighbor = neighbor + 1
! Check for sites out of bounds -- this only happens in the rare
! circumstance that a processor close to the end has so many sites that
! it would exceed the bank on the last processor
if (neighbor > n_procs - 1) exit
end do SEND_SITES
end if
! ==========================================================================
! RECEIVE BANK SITES FROM NEIGHBORS OR TEMPORARY BANK
start = work_index(rank)
index_local = 1
! Determine what process has the source sites that will need to be stored at
! the beginning of this processor's source bank.
if (start >= bank_position(n_procs)) then
neighbor = n_procs - 1
else
neighbor = binary_search(bank_position, n_procs, start) - 1
end if
RECV_SITES: do while (start < work_index(rank + 1))
! Determine how many sites need to be received
if (neighbor == n_procs - 1) then
n = work_index(rank + 1) - start
else
n = min(bank_position(neighbor + 2), work_index(rank + 1)) - start
end if
if (neighbor /= rank) then
! If the source sites are not on this processor, initiate an
! asynchronous receive for the source sites
n_request = n_request + 1
call MPI_IRECV(source_bank(index_local), int(n), MPI_BANK, &
neighbor, neighbor, mpi_intracomm, request(n_request), mpi_err)
else
! If the source sites are on this procesor, we can simply copy them
! from the temp_sites bank
index_temp = start - bank_position(rank+1) + 1
source_bank(index_local:index_local+n-1) = &
temp_sites(index_temp:index_temp+n-1)
end if
! Increment all indices
start = start + n
index_local = index_local + n
neighbor = neighbor + 1
end do RECV_SITES
! Since we initiated a series of asynchronous ISENDs and IRECVs, now we have
! to ensure that the data has actually been communicated before moving on to
! the next generation
call MPI_WAITALL(n_request, request, MPI_STATUSES_IGNORE, mpi_err)
! Deallocate space for bank_position on the very last generation
if (current_batch == n_max_batches .and. current_gen == gen_per_batch) &
deallocate(bank_position)
#else
source_bank = temp_sites(1:n_particles)
#endif
call time_bank_sendrecv % stop()
! Deallocate space for the temporary source bank on the last generation
if (current_batch == n_max_batches .and. current_gen == gen_per_batch) &
deallocate(temp_sites)
end subroutine synchronize_bank
!===============================================================================
! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto
! the master processor and normalizes them. This should work whether or not the

View file

@ -9,9 +9,14 @@
#include "openmc/hdf5_interface.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include <algorithm> // for min
#include <string>
namespace openmc {
//==============================================================================
@ -25,6 +30,245 @@ xt::xtensor<double, 1> source_frac;
// Non-member functions
//==============================================================================
void synchronize_bank()
{
// 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
// of fission bank sites -- each processor needs to know the total number of
// sites in order to figure out the probability for selecting
// sites. Furthermore, each proc also needs to know where in the 'global'
// fission bank its own sites starts in order to ensure reproducibility by
// skipping ahead to the proper seed.
#ifdef OPENMC_MPI
int64_t start = 0;
MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
// While we would expect the value of start on rank 0 to be 0, the MPI
// standard says that the receive buffer on rank 0 is undefined and not
// significant
if (mpi::rank == 0) start = 0;
int64_t finish = start + n_bank;
int64_t total = finish;
MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm);
#else
int64_t start = 0;
int64_t finish = n_bank;
int64_t total = n_bank;
#endif
// If there are not that many particles per generation, it's possible that no
// fission sites were created at all on a single processor. Rather than add
// extra logic to treat this circumstance, we really want to ensure the user
// runs enough particles to avoid this in the first place.
if (n_bank == 0) {
fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
}
// Make sure all processors start at the same point for random sampling. Then
// skip ahead in the sequence using the starting index in the 'global'
// fission bank for each processor.
set_particle_seed(simulation::total_gen + overall_generation());
advance_prn_seed(start);
// Determine how many fission sites we need to sample from the source bank
// and the probability for selecting a site.
int64_t sites_needed;
if (total < settings::n_particles) {
sites_needed = settings::n_particles % total;
} else {
sites_needed = settings::n_particles;
}
double p_sample = static_cast<double>(sites_needed) / total;
//time_bank_sample % start()
// ==========================================================================
// SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
// Allocate temporary source bank
int64_t index_temp = 0;
Bank temp_sites[3*simulation::work];
for (int64_t i = 0; i < n_bank; ++i) {
// If there are less than n_particles particles banked, automatically add
// int(n_particles/total) sites to temp_sites. For example, if you need
// 1000 and 300 were banked, this would add 3 source sites per banked site
// 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];
++index_temp;
}
}
// Randomly sample sites needed
if (prn() < p_sample) {
temp_sites[index_temp] = fission_bank[i];
++index_temp;
}
}
// At this point, the sampling of source sites is done and now we need to
// figure out where to send source sites. Since it is possible that one
// processor's share of the source bank spans more than just the immediate
// neighboring processors, we have to perform an ALLGATHER to determine the
// indices for all processors
#ifdef OPENMC_MPI
// First do an exclusive scan to get the starting indices for
start = 0;
MPI_Exscan(&index_temp, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
finish = start + index_temp;
// Allocate space for bank_position if this hasn't been done yet
int64_t bank_position[mpi::n_procs];
MPI_Allgather(&start, 1, MPI_INT64_T, bank_position, 1,
MPI_INT64_T, mpi::intracomm);
#else
start = 0;
finish = index_temp;
#endif
// Now that the sampling is complete, we need to ensure that we have exactly
// n_particles source sites. The way this is done in a reproducible manner is
// to adjust only the source sites on the last processor.
if (mpi::rank == mpi::n_procs - 1) {
if (finish > settings::n_particles) {
// If we have extra sites sampled, we will simply discard the extra
// ones on the last processor
index_temp = settings::n_particles - start;
} else if (finish < settings::n_particles) {
// If we have too few sites, repeat sites from the very end of the
// fission bank
sites_needed = settings::n_particles - finish;
for (int i = 0; i < sites_needed; ++i) {
temp_sites[index_temp] = fission_bank[n_bank - sites_needed + i];
++index_temp;
}
}
// the last processor should not be sending sites to right
finish = simulation::work_index[mpi::rank + 1];
}
//time_bank_sample % stop()
//time_bank_sendrecv % start()
#ifdef OPENMC_MPI
// ==========================================================================
// SEND BANK SITES TO NEIGHBORS
int64_t index_local = 0;
std::vector<MPI_Request> requests;
if (start < settings::n_particles) {
// Determine the index of the processor which has the first part of the
// source_bank for the local processor
int neighbor = upper_bound_index(simulation::work_index.begin(),
simulation::work_index.end(), start);
while (start < finish) {
// Determine the number of sites to send
int64_t n = std::min(simulation::work_index[neighbor + 1], finish) - start;
// Initiate an asynchronous send of source sites to the neighboring
// process
if (neighbor != mpi::rank) {
requests.emplace_back();
MPI_Isend(&temp_sites[index_local], static_cast<int>(n), mpi::bank,
neighbor, mpi::rank, mpi::intracomm, &requests.back());
}
// Increment all indices
start += n;
index_local += n;
++neighbor;
// Check for sites out of bounds -- this only happens in the rare
// circumstance that a processor close to the end has so many sites that
// it would exceed the bank on the last processor
if (neighbor > mpi::n_procs - 1) break;
}
}
// ==========================================================================
// RECEIVE BANK SITES FROM NEIGHBORS OR TEMPORARY BANK
start = simulation::work_index[mpi::rank];
index_local = 0;
// Determine what process has the source sites that will need to be stored at
// the beginning of this processor's source bank.
int neighbor;
if (start >= bank_position[mpi::n_procs - 1]) {
neighbor = mpi::n_procs - 1;
} else {
neighbor = upper_bound_index(bank_position, bank_position + mpi::n_procs, start);
}
while (start < simulation::work_index[mpi::rank + 1]) {
// Determine how many sites need to be received
int64_t n;
if (neighbor == mpi::n_procs - 1) {
n = simulation::work_index[mpi::rank + 1] - start;
} else {
n = std::min(bank_position[neighbor + 1], simulation::work_index[mpi::rank + 1]) - start;
}
if (neighbor != mpi::rank) {
// If the source sites are not on this processor, initiate an
// asynchronous receive for the source sites
requests.emplace_back();
MPI_Irecv(&source_bank[index_local], static_cast<int>(n), mpi::bank,
neighbor, neighbor, mpi::intracomm, &requests.back());
} else {
// If the source sites are on this procesor, we can simply copy them
// from the temp_sites bank
index_temp = start - bank_position[mpi::rank];
std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n],
&source_bank[index_local]);
}
// Increment all indices
start += n;
index_local += n;
++neighbor;
}
// Since we initiated a series of asynchronous ISENDs and IRECVs, now we have
// to ensure that the data has actually been communicated before moving on to
// the next generation
int n_request = requests.size();
MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
#else
std::copy(temp_sites, temp_sites + settings::n_particles, source_bank);
#endif
//time_bank_sendrecv % stop()
}
void shannon_entropy()
{
// Get pointer to entropy mesh

View file

@ -83,8 +83,9 @@ contains
call time_initialize%start()
#ifdef OPENMC_MPI
! Setup MPI
call initialize_mpi(comm)
! Indicate that MPI is turned on
mpi_enabled = .true.
mpi_intracomm = intracomm
#endif
#ifdef _OPENMP
@ -114,57 +115,6 @@ contains
err = 0
end function openmc_init_f
#ifdef OPENMC_MPI
!===============================================================================
! INITIALIZE_MPI starts up the Message Passing Interface (MPI) and determines
! the number of processors the problem is being run with as well as the rank of
! each processor.
!===============================================================================
subroutine initialize_mpi(intracomm)
#ifdef OPENMC_MPIF08
type(MPI_Comm), intent(in) :: intracomm ! MPI intracommunicator
#else
integer, intent(in) :: intracomm ! MPI intracommunicator
#endif
integer :: mpi_err ! MPI error code
integer :: bank_blocks(5) ! Count for each datatype
#ifdef OPENMC_MPIF08
type(MPI_Datatype) :: bank_types(5)
#else
integer :: bank_types(5) ! Datatypes
#endif
integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements
type(Bank) :: b
! Indicate that MPI is turned on
mpi_enabled = .true.
mpi_intracomm = intracomm
! ==========================================================================
! CREATE MPI_BANK TYPE
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(b % wgt, bank_disp(1), mpi_err)
call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err)
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err)
call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err)
call MPI_GET_ADDRESS(b % delayed_group, bank_disp(5), mpi_err)
! Adjust displacements
bank_disp = bank_disp - bank_disp(1)
! Define MPI_BANK for fission sites
bank_blocks = (/ 1, 3, 3, 1, 1 /)
bank_types = (/ MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_INTEGER /)
call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, &
bank_types, MPI_BANK, mpi_err)
call MPI_TYPE_COMMIT(MPI_BANK, mpi_err)
end subroutine initialize_mpi
#endif
!===============================================================================
! READ_COMMAND_LINE reads all parameters from the command line
!===============================================================================

View file

@ -80,16 +80,18 @@ void initialize_mpi(MPI_Comm intracomm)
// Create bank datatype
Bank b;
MPI_Aint disp[] {
offsetof(Bank, wgt),
offsetof(Bank, xyz),
offsetof(Bank, uvw),
offsetof(Bank, E),
offsetof(Bank, delayed_group)
};
int blocks[] {1, 3, 3, 1, 1};
MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT};
MPI_Type_create_struct(5, blocks, disp, types, &mpi::bank);
MPI_Aint disp[6];
MPI_Get_address(&b.wgt, &disp[0]);
MPI_Get_address(&b.xyz, &disp[1]);
MPI_Get_address(&b.uvw, &disp[2]);
MPI_Get_address(&b.E, &disp[3]);
MPI_Get_address(&b.delayed_group, &disp[4]);
MPI_Get_address(&b.particle, &disp[5]);
for (int i = 5; i >= 0; --i) disp[i] -= disp[0];
int blocks[] {1, 3, 3, 1, 1, 1};
MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT};
MPI_Type_create_struct(6, blocks, disp, types, &mpi::bank);
MPI_Type_commit(&mpi::bank);
}
#endif // OPENMC_MPI

View file

@ -19,10 +19,8 @@ module message_passing
logical(C_BOOL), bind(C, name='openmc_master') :: master = .true. ! master process?
logical :: mpi_enabled = .false. ! is MPI in use and initialized?
#ifdef OPENMC_MPIF08
type(MPI_Datatype) :: MPI_BANK ! MPI datatype for fission bank
type(MPI_Comm) :: mpi_intracomm ! MPI intra-communicator
#else
integer :: MPI_BANK ! MPI datatype for fission bank
integer :: mpi_intracomm ! MPI intra-communicator
#endif

View file

@ -11,7 +11,7 @@ module simulation
use cmfd_header, only: cmfd_on
use constants, only: ZERO
use eigenvalue, only: calculate_average_keff, calculate_generation_keff, &
synchronize_bank, keff_generation, k_sum
keff_generation, k_sum
#ifdef _OPENMP
use eigenvalue, only: join_bank_from_threads
#endif
@ -258,6 +258,9 @@ contains
subroutine shannon_entropy() bind(C)
end subroutine
subroutine synchronize_bank() bind(C)
end subroutine
end interface
! Update global tallies with the omp private accumulation variables

View file

@ -86,8 +86,8 @@ void calculate_work()
int64_t remainder = settings::n_particles % mpi::n_procs;
int64_t i_bank = 0;
simulation::work_index.reserve(mpi::n_procs);
simulation::work_index.push_back(0);
simulation::work_index.resize(mpi::n_procs + 1);
simulation::work_index[0] = 0;
for (int i = 0; i < mpi::n_procs; ++i) {
// Number of particles for rank i
int64_t work_i = i < remainder ? min_work + 1 : min_work;
@ -97,7 +97,7 @@ void calculate_work()
// Set index into source bank for rank i
i_bank += work_i;
simulation::work_index.push_back(i_bank);
simulation::work_index[i + 1] = i_bank;
}
}