From 02fbc1582d837e40e01c7aa4c82283440cfa1599 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 11 Oct 2018 12:55:33 -0500 Subject: [PATCH] Convert synchronize_bank to C++ --- include/openmc/eigenvalue.h | 3 + include/openmc/search.h | 10 +- src/api.F90 | 1 - src/eigenvalue.F90 | 271 ------------------------------------ src/eigenvalue.cpp | 244 ++++++++++++++++++++++++++++++++ src/initialize.F90 | 56 +------- src/initialize.cpp | 22 +-- src/message_passing.F90 | 2 - src/simulation.F90 | 5 +- src/simulation.cpp | 6 +- 10 files changed, 278 insertions(+), 342 deletions(-) diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index cc879461e..59b448e29 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -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(); diff --git a/include/openmc/search.h b/include/openmc/search.h index c91ad18fe..446e2916a 100644 --- a/include/openmc/search.h +++ b/include/openmc/search.h @@ -4,7 +4,7 @@ #ifndef OPENMC_SEARCH_H #define OPENMC_SEARCH_H -#include // for lower_bound +#include // 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 +typename std::iterator_traits::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 diff --git a/src/api.F90 b/src/api.F90 index 38329ad80..d79aa2450 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -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 diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index bd2a20b41..df108f6d7 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -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 diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index f03820492..459204ea3 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -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 // for min +#include + namespace openmc { //============================================================================== @@ -25,6 +30,245 @@ xt::xtensor 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(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 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(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(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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 0dbaa83bf..697fc1235 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 !=============================================================================== diff --git a/src/initialize.cpp b/src/initialize.cpp index ac7a5ef15..31fa8fb25 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -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 diff --git a/src/message_passing.F90 b/src/message_passing.F90 index 6ade0895a..05219f0f8 100644 --- a/src/message_passing.F90 +++ b/src/message_passing.F90 @@ -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 diff --git a/src/simulation.F90 b/src/simulation.F90 index 767ec5244..2359306d2 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -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 diff --git a/src/simulation.cpp b/src/simulation.cpp index 386e00510..2935d71ac 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -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; } }