mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge pull request #1094 from paulromano/cpp-mpi
Move MPI-related code to C++
This commit is contained in:
Sterling Harper
GitHub
commit
04fa0b27f9
48 changed files with 1177 additions and 947 deletions
|
|
@ -18,7 +18,6 @@ option(profile "Compile with profiling flags" OFF)
|
|||
option(debug "Compile with debug flags" OFF)
|
||||
option(optimize "Turn on all compiler optimization flags" OFF)
|
||||
option(coverage "Compile with coverage analysis flags" OFF)
|
||||
option(mpif08 "Use Fortran 2008 MPI interface" OFF)
|
||||
option(dagmc "Enable support for DAGMC (CAD) geometry" OFF)
|
||||
|
||||
# Maximum number of nested coordinates levels
|
||||
|
|
@ -407,6 +406,7 @@ add_library(libopenmc SHARED
|
|||
src/distribution_spatial.cpp
|
||||
src/eigenvalue.cpp
|
||||
src/endf.cpp
|
||||
src/error.cpp
|
||||
src/initialize.cpp
|
||||
src/finalize.cpp
|
||||
src/geometry.cpp
|
||||
|
|
@ -442,9 +442,11 @@ add_library(libopenmc SHARED
|
|||
src/string_functions.cpp
|
||||
src/summary.cpp
|
||||
src/surface.cpp
|
||||
src/timer.cpp
|
||||
src/thermal.cpp
|
||||
src/xml_interface.cpp
|
||||
src/xsdata.cpp)
|
||||
src/xsdata.cpp
|
||||
src/tallies/tally.cpp)
|
||||
|
||||
set_target_properties(libopenmc PROPERTIES
|
||||
OUTPUT_NAME openmc
|
||||
|
|
@ -471,9 +473,6 @@ if (HDF5_IS_PARALLEL)
|
|||
endif()
|
||||
if (MPI_ENABLED)
|
||||
target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI)
|
||||
if (mpif08)
|
||||
target_compile_definitions(libopenmc PRIVATE -DOPENMC_MPIF08)
|
||||
endif()
|
||||
endif()
|
||||
|
||||
# Set git SHA1 hash as a compile definition
|
||||
|
|
|
|||
|
|
@ -3,6 +3,7 @@
|
|||
|
||||
#include <stdint.h>
|
||||
#include <stdbool.h>
|
||||
#include <stddef.h>
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
|
|
@ -48,9 +49,9 @@ extern "C" {
|
|||
int64_t openmc_get_seed();
|
||||
int openmc_get_tally_index(int32_t id, int32_t* index);
|
||||
void openmc_get_tally_next_id(int32_t* id);
|
||||
int openmc_global_tallies(double** ptr);
|
||||
int openmc_hard_reset();
|
||||
int openmc_init(int argc, char* argv[], const void* intracomm);
|
||||
int openmc_init_f(const int* intracomm);
|
||||
int openmc_legendre_filter_get_order(int32_t index, int* order);
|
||||
int openmc_legendre_filter_set_order(int32_t index, int order);
|
||||
int openmc_load_nuclide(const char name[]);
|
||||
|
|
@ -105,7 +106,7 @@ extern "C" {
|
|||
int openmc_tally_get_scores(int32_t index, int** scores, int* n);
|
||||
int openmc_tally_get_type(int32_t index, int32_t* type);
|
||||
int openmc_tally_reset(int32_t index);
|
||||
int openmc_tally_results(int32_t index, double** ptr, int shape_[3]);
|
||||
int openmc_tally_results(int32_t index, double** ptr, size_t shape_[3]);
|
||||
int openmc_tally_set_active(int32_t index, bool active);
|
||||
int openmc_tally_set_estimator(int32_t index, const char* estimator);
|
||||
int openmc_tally_set_filters(int32_t index, int n, const int32_t* indices);
|
||||
|
|
@ -134,8 +135,6 @@ extern "C" {
|
|||
|
||||
// Global variables
|
||||
extern char openmc_err_msg[256];
|
||||
extern double openmc_keff;
|
||||
extern double openmc_keff_std;
|
||||
extern int32_t n_cells;
|
||||
extern int32_t n_filters;
|
||||
extern int32_t n_lattices;
|
||||
|
|
@ -148,17 +147,12 @@ extern "C" {
|
|||
extern int32_t n_surfaces;
|
||||
extern int32_t n_tallies;
|
||||
extern int32_t n_universes;
|
||||
extern bool openmc_simulation_initialized;
|
||||
extern double global_tally_absorption;
|
||||
#pragma omp threadprivate(global_tally_absorption)
|
||||
|
||||
// Variables that are shared by necessity (can be removed from public header
|
||||
// later)
|
||||
extern bool openmc_master;
|
||||
extern int openmc_n_procs;
|
||||
extern int openmc_n_threads;
|
||||
extern int openmc_rank;
|
||||
extern int64_t openmc_work;
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
|
|
|
|||
|
|
@ -314,9 +314,9 @@ constexpr int MG_GET_XS_CHI_DELAYED {14};
|
|||
// TALLY-RELATED CONSTANTS
|
||||
|
||||
// Tally result entries
|
||||
constexpr int RESULT_VALUE {1};
|
||||
constexpr int RESULT_SUM {2};
|
||||
constexpr int RESULT_SUM_SQ {3};
|
||||
constexpr int RESULT_VALUE {0};
|
||||
constexpr int RESULT_SUM {1};
|
||||
constexpr int RESULT_SUM_SQ {2};
|
||||
|
||||
// Tally type
|
||||
// TODO: Convert to enum
|
||||
|
|
@ -407,10 +407,11 @@ constexpr int RELATIVE_ERROR {2};
|
|||
constexpr int STANDARD_DEVIATION {3};
|
||||
|
||||
// Global tally parameters
|
||||
constexpr int K_COLLISION {1};
|
||||
constexpr int K_ABSORPTION {2};
|
||||
constexpr int K_TRACKLENGTH {3};
|
||||
constexpr int LEAKAGE {4};
|
||||
constexpr int N_GLOBAL_TALLIES {4};
|
||||
constexpr int K_COLLISION {0};
|
||||
constexpr int K_ABSORPTION {1};
|
||||
constexpr int K_TRACKLENGTH {2};
|
||||
constexpr int LEAKAGE {3};
|
||||
|
||||
// Differential tally independent variables
|
||||
constexpr int DIFF_DENSITY {1};
|
||||
|
|
|
|||
|
|
@ -14,6 +14,7 @@ namespace openmc {
|
|||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
extern double keff_generation; //!< Single-generation k on each processor
|
||||
extern std::vector<double> entropy; //!< Shannon entropy at each generation
|
||||
extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
|
||||
|
||||
|
|
@ -24,6 +25,12 @@ extern "C" int64_t n_bank;
|
|||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
//! Collect/normalize the tracklength keff from each process
|
||||
extern "C" void calculate_generation_keff();
|
||||
|
||||
//! 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();
|
||||
|
|
|
|||
|
|
@ -9,7 +9,6 @@
|
|||
|
||||
namespace openmc {
|
||||
|
||||
|
||||
extern "C" void fatal_error_from_c(const char* message, int message_len);
|
||||
extern "C" void warning_from_c(const char* message, int message_len);
|
||||
extern "C" void write_message_from_c(const char* message, int message_len,
|
||||
|
|
@ -81,5 +80,9 @@ void write_message(const std::stringstream& message, int level)
|
|||
write_message(message.str(), level);
|
||||
}
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
extern "C" void abort_mpi(int code);
|
||||
#endif
|
||||
|
||||
} // namespace openmc
|
||||
#endif // OPENMC_ERROR_H
|
||||
|
|
|
|||
|
|
@ -1,6 +1,10 @@
|
|||
#ifndef OPENMC_FINALIZE_H
|
||||
#define OPENMC_FINALIZE_H
|
||||
|
||||
namespace openmc {
|
||||
|
||||
extern "C" void openmc_free_bank();
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_FINALIZE_H
|
||||
|
|
|
|||
|
|
@ -8,6 +8,7 @@
|
|||
#include <cstring> // for strlen
|
||||
#include <string>
|
||||
#include <sstream>
|
||||
#include <type_traits>
|
||||
#include <vector>
|
||||
|
||||
#include "hdf5.h"
|
||||
|
|
@ -333,7 +334,9 @@ write_attribute(hid_t obj_id, const char* name, const std::vector<T>& buffer)
|
|||
// Templates/overloads for write_dataset
|
||||
//==============================================================================
|
||||
|
||||
template<typename T> inline void
|
||||
// Template for scalars (ensured by SFINAE)
|
||||
template<typename T> inline
|
||||
std::enable_if_t<std::is_scalar<std::decay_t<T>>::value>
|
||||
write_dataset(hid_t obj_id, const char* name, T buffer)
|
||||
{
|
||||
write_dataset(obj_id, 0, nullptr, name, H5TypeMap<T>::type_id, &buffer, false);
|
||||
|
|
@ -359,9 +362,11 @@ write_dataset(hid_t obj_id, const char* name, const std::vector<T>& buffer)
|
|||
write_dataset(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data(), false);
|
||||
}
|
||||
|
||||
template<typename T> inline void
|
||||
write_dataset(hid_t obj_id, const char* name, const xt::xarray<T>& arr)
|
||||
// Template for xarray, xtensor, etc.
|
||||
template<typename D> inline void
|
||||
write_dataset(hid_t obj_id, const char* name, const xt::xcontainer<D>& arr)
|
||||
{
|
||||
using T = typename D::value_type;
|
||||
auto s = arr.shape();
|
||||
std::vector<hsize_t> dims {s.cbegin(), s.cend()};
|
||||
write_dataset(obj_id, dims.size(), dims.data(), name, H5TypeMap<T>::type_id,
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@
|
|||
#define OPENMC_MESSAGE_PASSING_H
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
#include "mpi.h"
|
||||
#include <mpi.h>
|
||||
#endif
|
||||
|
||||
namespace openmc {
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -10,29 +10,62 @@
|
|||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
// Global variable declarations
|
||||
//==============================================================================
|
||||
|
||||
extern "C" int openmc_current_batch;
|
||||
extern "C" int openmc_current_gen;
|
||||
extern "C" int64_t openmc_current_work;
|
||||
extern "C" int openmc_n_lost_particles;
|
||||
extern "C" int openmc_total_gen;
|
||||
extern "C" bool openmc_trace;
|
||||
namespace simulation {
|
||||
|
||||
extern "C" int current_batch; //!< current batch
|
||||
extern "C" int current_gen; //!< current fission generation
|
||||
extern "C" int64_t current_work; //!< index in source back of current particle
|
||||
extern "C" double keff; //!< average k over batches
|
||||
extern "C" double keff_std; //!< standard deviation of average k
|
||||
extern "C" double k_col_abs; //!< sum over batches of k_collision * k_absorption
|
||||
extern "C" double k_col_tra; //!< sum over batches of k_collision * k_tracklength
|
||||
extern "C" double k_abs_tra; //!< sum over batches of k_absorption * k_tracklength
|
||||
extern "C" double log_spacing; //!< lethargy spacing for energy grid searches
|
||||
extern "C" int n_lost_particles; //!< cumulative number of lost particles
|
||||
extern "C" bool need_depletion_rx; //!< need to calculate depletion rx?
|
||||
extern "C" int restart_batch; //!< batch at which a restart job resumed
|
||||
extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
|
||||
extern "C" bool simulation_initialized; //!< has simulation been initialized?
|
||||
extern "C" int total_gen; //!< total number of generations simulated
|
||||
extern "C" int64_t work; //!< number of particles per process
|
||||
|
||||
extern std::vector<double> k_generation;
|
||||
extern std::vector<int64_t> work_index;
|
||||
|
||||
#pragma omp threadprivate(openmc_current_work, openmc_trace)
|
||||
// Threadprivate variables
|
||||
extern "C" bool trace; //!< flag to show debug information
|
||||
#ifdef _OPENMP
|
||||
extern "C" int n_threads; //!< number of OpenMP threads
|
||||
extern "C" int thread_id; //!< ID of a given thread
|
||||
#endif
|
||||
|
||||
#pragma omp threadprivate(current_work, thread_id, trace)
|
||||
|
||||
} // namespace simulation
|
||||
|
||||
//==============================================================================
|
||||
// Functions
|
||||
//==============================================================================
|
||||
|
||||
//! Determine number of particles to transport per process
|
||||
void calculate_work();
|
||||
|
||||
//! Initialize simulation
|
||||
extern "C" void openmc_simulation_init_c();
|
||||
|
||||
//! Determine number of particles to transport per process
|
||||
void calculate_work();
|
||||
//! Initialize a fission generation
|
||||
extern "C" void initialize_generation();
|
||||
|
||||
//! Determine overall generation number
|
||||
extern "C" int overall_generation();
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
extern "C" void broadcast_results();
|
||||
extern "C" void broadcast_triggers();
|
||||
#endif
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
|
|
|
|||
|
|
@ -9,10 +9,10 @@
|
|||
|
||||
namespace openmc {
|
||||
|
||||
extern "C" void write_source_bank(hid_t group_id, int64_t* work_index,
|
||||
Bank* source_bank);
|
||||
extern "C" void read_source_bank(hid_t group_id, int64_t* work_index,
|
||||
Bank* source_bank);
|
||||
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_tally_results_nr(hid_t file_id);
|
||||
|
||||
} // namespace openmc
|
||||
#endif // OPENMC_STATE_POINT_H
|
||||
|
|
|
|||
46
include/openmc/tallies/tally.h
Normal file
46
include/openmc/tallies/tally.h
Normal file
|
|
@ -0,0 +1,46 @@
|
|||
#ifndef OPENMC_TALLIES_TALLY_H
|
||||
#define OPENMC_TALLIES_TALLY_H
|
||||
|
||||
#include "openmc/constants.h"
|
||||
|
||||
#include "xtensor/xtensor.hpp"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// Global variable declarations
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double total_weight;
|
||||
|
||||
// Threadprivate variables
|
||||
|
||||
extern "C" double global_tally_absorption;
|
||||
#pragma omp threadprivate(global_tally_absorption)
|
||||
|
||||
//==============================================================================
|
||||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
// Alias for the type returned by xt::adapt(...). N is the dimension of the
|
||||
// multidimensional array
|
||||
template <std::size_t N>
|
||||
using adaptor_type = xt::xtensor_adaptor<xt::xbuffer_adaptor<double*&, xt::no_ownership>, N>;
|
||||
|
||||
//! Get the global tallies as a multidimensional array
|
||||
//! \return Global tallies array
|
||||
adaptor_type<2> global_tallies();
|
||||
|
||||
//! Get tally results as a multidimensional array
|
||||
//! \param idx Index in tallies array
|
||||
//! \return Tally results array
|
||||
adaptor_type<3> tally_results(int idx);
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
//! Collect all tally results onto master process
|
||||
extern "C" void reduce_tally_results();
|
||||
#endif
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_TALLIES_TALLY_H
|
||||
47
include/openmc/timer.h
Normal file
47
include/openmc/timer.h
Normal file
|
|
@ -0,0 +1,47 @@
|
|||
#ifndef OPENMC_TIMER_H
|
||||
#define OPENMC_TIMER_H
|
||||
|
||||
#include <chrono>
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
//! Class for measuring time elapsed
|
||||
//==============================================================================
|
||||
|
||||
class Timer {
|
||||
public:
|
||||
using clock = std::chrono::high_resolution_clock;
|
||||
|
||||
Timer() {};
|
||||
|
||||
//! Start running the timer
|
||||
void start ();
|
||||
|
||||
//! Get total elapsed time in seconds
|
||||
//! \return Elapsed time in [s]
|
||||
double elapsed();
|
||||
|
||||
//! Stop running the timer
|
||||
void stop();
|
||||
|
||||
//! Stop the timer and reset its elapsed time
|
||||
void reset();
|
||||
|
||||
private:
|
||||
bool running_ {false}; //!< is timer running?
|
||||
std::chrono::time_point<clock> start_; //!< starting point for clock
|
||||
double elapsed_; //!< elasped time in [s]
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
extern Timer time_bank;
|
||||
extern Timer time_bank_sample;
|
||||
extern Timer time_bank_sendrecv;
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_TIMER_H
|
||||
|
|
@ -71,7 +71,7 @@ def current_batch():
|
|||
Current batch of the simulation
|
||||
|
||||
"""
|
||||
return c_int.in_dll(_dll, 'openmc_current_batch').value
|
||||
return c_int.in_dll(_dll, 'current_batch').value
|
||||
|
||||
|
||||
def finalize():
|
||||
|
|
@ -220,8 +220,8 @@ def keff():
|
|||
return tuple(k)
|
||||
else:
|
||||
# Otherwise, return the tracklength estimator
|
||||
mean = c_double.in_dll(_dll, 'openmc_keff').value
|
||||
std_dev = c_double.in_dll(_dll, 'openmc_keff_std').value \
|
||||
mean = c_double.in_dll(_dll, 'keff').value
|
||||
std_dev = c_double.in_dll(_dll, 'keff_std').value \
|
||||
if n > 1 else np.inf
|
||||
return (mean, std_dev)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
from collections.abc import Mapping
|
||||
from ctypes import c_int, c_int32, c_double, c_char_p, c_bool, POINTER
|
||||
from ctypes import c_int, c_int32, c_size_t, c_double, c_char_p, c_bool, POINTER
|
||||
from weakref import WeakValueDictionary
|
||||
|
||||
import numpy as np
|
||||
|
|
@ -60,7 +60,7 @@ _dll.openmc_tally_reset.argtypes = [c_int32]
|
|||
_dll.openmc_tally_reset.restype = c_int
|
||||
_dll.openmc_tally_reset.errcheck = _error_handler
|
||||
_dll.openmc_tally_results.argtypes = [
|
||||
c_int32, POINTER(POINTER(c_double)), POINTER(c_int*3)]
|
||||
c_int32, POINTER(POINTER(c_double)), POINTER(c_size_t*3)]
|
||||
_dll.openmc_tally_results.restype = c_int
|
||||
_dll.openmc_tally_results.errcheck = _error_handler
|
||||
_dll.openmc_tally_set_active.argtypes = [c_int32, c_bool]
|
||||
|
|
@ -296,9 +296,9 @@ class Tally(_FortranObjectWithID):
|
|||
@property
|
||||
def results(self):
|
||||
data = POINTER(c_double)()
|
||||
shape = (c_int*3)()
|
||||
shape = (c_size_t*3)()
|
||||
_dll.openmc_tally_results(self._index, data, shape)
|
||||
return as_array(data, tuple(shape[::-1]))
|
||||
return as_array(data, tuple(shape))
|
||||
|
||||
@property
|
||||
def scores(self):
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
@ -277,14 +276,12 @@ contains
|
|||
call time_initialize % reset()
|
||||
call time_read_xs % reset()
|
||||
call time_unionize % reset()
|
||||
call time_bank % reset()
|
||||
call time_bank_sample % reset()
|
||||
call time_bank_sendrecv % reset()
|
||||
call time_tallies % reset()
|
||||
call time_inactive % reset()
|
||||
call time_active % reset()
|
||||
call time_transport % reset()
|
||||
call time_finalize % reset()
|
||||
call reset_timers()
|
||||
|
||||
err = 0
|
||||
end function openmc_reset
|
||||
|
|
|
|||
|
|
@ -13,6 +13,16 @@ module cmfd_execute
|
|||
private
|
||||
public :: execute_cmfd, cmfd_init_batch, cmfd_tally_init
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
interface
|
||||
subroutine cmfd_broadcast(n, buffer) bind(C)
|
||||
import C_DOUBLE, C_INT
|
||||
integer(C_INT), value :: n
|
||||
real(C_DOUBLE), intent(out) :: buffer
|
||||
end subroutine
|
||||
end interface
|
||||
#endif
|
||||
|
||||
contains
|
||||
|
||||
!==============================================================================
|
||||
|
|
@ -105,9 +115,6 @@ contains
|
|||
real(8) :: hxyz(3) ! cell dimensions of current ijk cell
|
||||
real(8) :: vol ! volume of cell
|
||||
real(8),allocatable :: source(:,:,:,:) ! tmp source array for entropy
|
||||
#ifdef OPENMC_MPI
|
||||
integer :: mpi_err ! MPI error code
|
||||
#endif
|
||||
|
||||
! Get maximum of spatial and group indices
|
||||
nx = cmfd % indices(1)
|
||||
|
|
@ -199,7 +206,7 @@ contains
|
|||
|
||||
#ifdef OPENMC_MPI
|
||||
! Broadcast full source to all procs
|
||||
call MPI_BCAST(cmfd % cmfd_src, n, MPI_REAL8, 0, mpi_intracomm, mpi_err)
|
||||
call cmfd_broadcast(n, cmfd % cmfd_src(1,1,1,1))
|
||||
#endif
|
||||
|
||||
end subroutine calc_fission_source
|
||||
|
|
@ -232,9 +239,6 @@ contains
|
|||
real(8) :: norm ! normalization factor
|
||||
logical(C_BOOL) :: outside ! any source sites outside mesh
|
||||
logical :: in_mesh ! source site is inside mesh
|
||||
#ifdef OPENMC_MPI
|
||||
integer :: mpi_err
|
||||
#endif
|
||||
|
||||
interface
|
||||
subroutine cmfd_populate_sourcecounts(ng, energies, source_counts, outside) bind(C)
|
||||
|
|
@ -300,8 +304,7 @@ contains
|
|||
|
||||
! Broadcast weight factors to all procs
|
||||
#ifdef OPENMC_MPI
|
||||
call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, &
|
||||
mpi_intracomm, mpi_err)
|
||||
call cmfd_broadcast(ng*nx*ny*nz, cmfd % weightfactors(1,1,1,1))
|
||||
#endif
|
||||
end if
|
||||
|
||||
|
|
|
|||
|
|
@ -7,6 +7,8 @@
|
|||
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/mesh.h"
|
||||
#include "openmc/message_passing.h"
|
||||
#include "openmc/simulation.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -23,10 +25,19 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies,
|
|||
|
||||
// Get source counts in each mesh bin / energy bin
|
||||
auto& m = meshes.at(index_cmfd_mesh);
|
||||
xt::xarray<double> counts = m->count_sites(openmc_work, source_bank, n_energy, energies, outside);
|
||||
xt::xarray<double> counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside);
|
||||
|
||||
// Copy data from the xarray into the source counts array
|
||||
std::copy(counts.begin(), counts.end(), source_counts);
|
||||
}
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
extern "C" void
|
||||
cmfd_broadcast(int n, double* buffer)
|
||||
{
|
||||
MPI_Bcast(buffer, n, MPI_DOUBLE, 0, mpi::intracomm);
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -51,8 +51,8 @@ module cmfd_header
|
|||
real(8), allocatable :: hxyz(:,:,:,:)
|
||||
|
||||
! Source distributions
|
||||
real(8), allocatable :: cmfd_src(:,:,:,:)
|
||||
real(8), allocatable :: openmc_src(:,:,:,:)
|
||||
real(C_DOUBLE), allocatable :: cmfd_src(:,:,:,:)
|
||||
real(C_DOUBLE), allocatable :: openmc_src(:,:,:,:)
|
||||
|
||||
! Source sites in each mesh box
|
||||
real(C_DOUBLE), allocatable :: sourcecounts(:,:)
|
||||
|
|
|
|||
|
|
@ -16,313 +16,15 @@ module eigenvalue
|
|||
|
||||
implicit none
|
||||
|
||||
real(8) :: keff_generation ! Single-generation k on each processor
|
||||
real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq
|
||||
|
||||
interface
|
||||
subroutine calculate_generation_keff() bind(C)
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
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
|
||||
! no-reduce method is being used.
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_generation_keff()
|
||||
|
||||
real(8) :: keff_reduced
|
||||
#ifdef OPENMC_MPI
|
||||
integer :: mpi_err ! MPI error code
|
||||
#endif
|
||||
|
||||
! Get keff for this generation by subtracting off the starting value
|
||||
keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
! Combine values across all processors
|
||||
call MPI_ALLREDUCE(keff_generation, keff_reduced, 1, MPI_REAL8, &
|
||||
MPI_SUM, mpi_intracomm, mpi_err)
|
||||
#else
|
||||
keff_reduced = keff_generation
|
||||
#endif
|
||||
|
||||
! Normalize single batch estimate of k
|
||||
! TODO: This should be normalized by total_weight, not by n_particles
|
||||
keff_reduced = keff_reduced / n_particles
|
||||
call k_generation % push_back(keff_reduced)
|
||||
|
||||
end subroutine calculate_generation_keff
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_AVERAGE_KEFF calculates the mean and standard deviation of the mean
|
||||
! of k-effective during active generations and broadcasts the mean to all
|
||||
|
|
@ -347,12 +49,12 @@ contains
|
|||
if (n <= 0) then
|
||||
! For inactive generations, use current generation k as estimate for next
|
||||
! generation
|
||||
keff = k_generation % data(i)
|
||||
keff = k_generation(i)
|
||||
|
||||
else
|
||||
! Sample mean of keff
|
||||
k_sum(1) = k_sum(1) + k_generation % data(i)
|
||||
k_sum(2) = k_sum(2) + k_generation % data(i)**2
|
||||
k_sum(1) = k_sum(1) + k_generation(i)
|
||||
k_sum(2) = k_sum(2) + k_generation(i)**2
|
||||
|
||||
! Determine mean
|
||||
keff = k_sum(1) / n
|
||||
|
|
|
|||
|
|
@ -5,12 +5,20 @@
|
|||
#include "xtensor/xview.hpp"
|
||||
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/error.h"
|
||||
#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 "openmc/timer.h"
|
||||
#include "openmc/tallies/tally.h"
|
||||
|
||||
#include <algorithm> // for min
|
||||
#include <string>
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -18,6 +26,7 @@ namespace openmc {
|
|||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
double keff_generation;
|
||||
std::vector<double> entropy;
|
||||
xt::xtensor<double, 1> source_frac;
|
||||
|
||||
|
|
@ -25,6 +34,270 @@ xt::xtensor<double, 1> source_frac;
|
|||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
void calculate_generation_keff()
|
||||
{
|
||||
auto gt = global_tallies();
|
||||
|
||||
// Get keff for this generation by subtracting off the starting value
|
||||
keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - keff_generation;
|
||||
|
||||
double keff_reduced;
|
||||
#ifdef OPENMC_MPI
|
||||
// Combine values across all processors
|
||||
MPI_Allreduce(&keff_generation, &keff_reduced, 1, MPI_DOUBLE,
|
||||
MPI_SUM, mpi::intracomm);
|
||||
#else
|
||||
keff_reduced = keff_generation;
|
||||
#endif
|
||||
|
||||
// Normalize single batch estimate of k
|
||||
// TODO: This should be normalized by total_weight, not by n_particles
|
||||
keff_reduced /= settings::n_particles;
|
||||
simulation::k_generation.push_back(keff_reduced);
|
||||
}
|
||||
|
||||
void synchronize_bank()
|
||||
{
|
||||
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
|
||||
// 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();
|
||||
time_bank.stop();
|
||||
}
|
||||
|
||||
void shannon_entropy()
|
||||
{
|
||||
// Get pointer to entropy mesh
|
||||
|
|
@ -66,7 +339,7 @@ void ufs_count_sites()
|
|||
{
|
||||
auto &m = meshes[settings::index_ufs_mesh];
|
||||
|
||||
if (openmc_current_batch == 1 && openmc_current_gen == 1) {
|
||||
if (simulation::current_batch == 1 && simulation::current_gen == 1) {
|
||||
// On the first generation, just assume that the source is already evenly
|
||||
// distributed so that effectively the production of fission sites is not
|
||||
// biased
|
||||
|
|
@ -82,7 +355,7 @@ void ufs_count_sites()
|
|||
|
||||
// count number of source sites in each ufs mesh cell
|
||||
bool sites_outside;
|
||||
source_frac = m->count_sites(openmc_work, source_bank, 0, nullptr,
|
||||
source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr,
|
||||
&sites_outside);
|
||||
|
||||
// Check for sites outside of the mesh
|
||||
|
|
@ -102,7 +375,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 < openmc_work; ++i) {
|
||||
for (int i = 0; i < simulation::work; ++i) {
|
||||
source_bank[i].wgt *= settings::n_particles / total;
|
||||
}
|
||||
}
|
||||
|
|
@ -127,18 +400,34 @@ double ufs_get_weight(const Particle* p)
|
|||
}
|
||||
}
|
||||
|
||||
extern "C" void entropy_to_hdf5(hid_t group)
|
||||
extern "C" void write_eigenvalue_hdf5(hid_t group)
|
||||
{
|
||||
write_dataset(group, "n_inactive", settings::n_inactive);
|
||||
write_dataset(group, "generations_per_batch", settings::gen_per_batch);
|
||||
write_dataset(group, "k_generation", simulation::k_generation);
|
||||
if (settings::entropy_on) {
|
||||
write_dataset(group, "entropy", entropy);
|
||||
}
|
||||
write_dataset(group, "k_col_abs", simulation::k_col_abs);
|
||||
write_dataset(group, "k_col_tra", simulation::k_col_tra);
|
||||
write_dataset(group, "k_abs_tra", simulation::k_abs_tra);
|
||||
std::array<double, 2> k_combined;
|
||||
openmc_get_keff(k_combined.data());
|
||||
write_dataset(group, "k_combined", k_combined);
|
||||
}
|
||||
|
||||
extern "C" void entropy_from_hdf5(hid_t group)
|
||||
extern "C" void read_eigenvalue_hdf5(hid_t group)
|
||||
{
|
||||
read_dataset(group, "generations_per_batch", settings::gen_per_batch);
|
||||
int n = simulation::restart_batch*settings::gen_per_batch;
|
||||
simulation::k_generation.resize(n);
|
||||
read_dataset(group, "k_generation", simulation::k_generation);
|
||||
if (settings::entropy_on) {
|
||||
read_dataset(group, "entropy", entropy);
|
||||
}
|
||||
read_dataset(group, "k_col_abs", simulation::k_col_abs);
|
||||
read_dataset(group, "k_col_tra", simulation::k_col_tra);
|
||||
read_dataset(group, "k_abs_tra", simulation::k_abs_tra);
|
||||
}
|
||||
|
||||
extern "C" double entropy_c(int i)
|
||||
|
|
|
|||
|
|
@ -130,14 +130,20 @@ contains
|
|||
character(*) :: message
|
||||
integer, optional :: error_code ! error code
|
||||
|
||||
integer :: code ! error code
|
||||
integer(C_INT) :: code ! error code
|
||||
integer :: i_start ! starting position
|
||||
integer :: i_end ! ending position
|
||||
integer :: line_wrap ! length of line
|
||||
integer :: length ! length of message
|
||||
integer :: indent ! length of indentation
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
integer :: mpi_err
|
||||
interface
|
||||
subroutine abort_mpi(code) bind(C)
|
||||
import C_INT
|
||||
integer(C_INT), value :: code
|
||||
end subroutine
|
||||
end interface
|
||||
#endif
|
||||
|
||||
|
||||
|
|
@ -190,7 +196,7 @@ contains
|
|||
|
||||
#ifdef OPENMC_MPI
|
||||
! Abort MPI
|
||||
call MPI_ABORT(mpi_intracomm, code, mpi_err)
|
||||
call abort_mpi(code)
|
||||
#endif
|
||||
|
||||
! Abort program
|
||||
|
|
|
|||
14
src/error.cpp
Normal file
14
src/error.cpp
Normal file
|
|
@ -0,0 +1,14 @@
|
|||
#include "openmc/error.h"
|
||||
|
||||
#include "openmc/message_passing.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
void abort_mpi(int code)
|
||||
{
|
||||
MPI_Abort(mpi::intracomm, code);
|
||||
}
|
||||
#endif
|
||||
|
||||
} // namespace openmc
|
||||
|
|
@ -2,9 +2,13 @@
|
|||
|
||||
#include "openmc/message_passing.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
void openmc_free_bank()
|
||||
{
|
||||
#ifdef OPENMC_MPI
|
||||
MPI_Type_free(&openmc::mpi::bank);
|
||||
MPI_Type_free(&mpi::bank);
|
||||
#endif
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -91,7 +91,7 @@ find_cell(Particle* p, int search_surf) {
|
|||
if (cells[i_cell]->contains(r, u, surf)) {
|
||||
p->coord[p->n_coord-1].cell = i_cell;
|
||||
|
||||
if (settings::verbosity >= 10 || openmc_trace) {
|
||||
if (settings::verbosity >= 10 || simulation::trace) {
|
||||
std::stringstream msg;
|
||||
msg << " Entering cell " << cells[i_cell]->id_;
|
||||
write_message(msg, 1);
|
||||
|
|
@ -255,7 +255,7 @@ cross_lattice(Particle* p, int lattice_translation[3])
|
|||
{
|
||||
Lattice& lat {*lattices[p->coord[p->n_coord-1].lattice-1]};
|
||||
|
||||
if (settings::verbosity >= 10 || openmc_trace) {
|
||||
if (settings::verbosity >= 10 || simulation::trace) {
|
||||
std::stringstream msg;
|
||||
msg << " Crossing lattice " << lat.id_ << ". Current position ("
|
||||
<< p->coord[p->n_coord-1].lattice_x << ","
|
||||
|
|
|
|||
|
|
@ -13,9 +13,6 @@ module hdf5_interface
|
|||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use error, only: fatal_error
|
||||
#ifdef PHDF5
|
||||
use message_passing, only: mpi_intracomm, MPI_INFO_NULL
|
||||
#endif
|
||||
use string, only: to_c_string, to_f_string
|
||||
|
||||
implicit none
|
||||
|
|
|
|||
|
|
@ -49,44 +49,17 @@ contains
|
|||
! setting up timers, etc.
|
||||
!===============================================================================
|
||||
|
||||
function openmc_init_f(intracomm) result(err) bind(C)
|
||||
integer, intent(in), optional :: intracomm ! MPI intracommunicator
|
||||
function openmc_init_f() result(err) bind(C)
|
||||
integer(C_INT) :: err
|
||||
|
||||
#ifdef _OPENMP
|
||||
character(MAX_WORD_LEN) :: envvar
|
||||
#endif
|
||||
|
||||
! Copy the communicator to a new variable. This is done to avoid changing
|
||||
! the signature of this subroutine. If MPI is being used but no communicator
|
||||
! was passed, assume MPI_COMM_WORLD.
|
||||
#ifdef OPENMC_MPI
|
||||
#ifdef OPENMC_MPIF08
|
||||
type(MPI_Comm), intent(in) :: comm ! MPI intracommunicator
|
||||
if (present(intracomm)) then
|
||||
comm % MPI_VAL = intracomm
|
||||
else
|
||||
comm = MPI_COMM_WORLD
|
||||
end if
|
||||
#else
|
||||
integer :: comm
|
||||
if (present(intracomm)) then
|
||||
comm = intracomm
|
||||
else
|
||||
comm = MPI_COMM_WORLD
|
||||
end if
|
||||
#endif
|
||||
#endif
|
||||
|
||||
! Start total and initialization timer
|
||||
call time_total%start()
|
||||
call time_initialize%start()
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
! Setup MPI
|
||||
call initialize_mpi(comm)
|
||||
#endif
|
||||
|
||||
#ifdef _OPENMP
|
||||
! Change schedule of main parallel-do loop if OMP_SCHEDULE is set
|
||||
call get_environment_variable("OMP_SCHEDULE", envvar)
|
||||
|
|
@ -114,57 +87,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
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -15,9 +15,11 @@
|
|||
#include "openmc/hdf5_interface.h"
|
||||
#include "openmc/message_passing.h"
|
||||
#include "openmc/settings.h"
|
||||
#include "openmc/simulation.h"
|
||||
#include "openmc/string_utils.h"
|
||||
|
||||
// data/functions from Fortran side
|
||||
extern "C" void openmc_init_f();
|
||||
extern "C" void print_usage();
|
||||
extern "C" void print_version();
|
||||
|
||||
|
|
@ -46,12 +48,7 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
|
|||
if (err) return err;
|
||||
|
||||
// Continue with rest of initialization
|
||||
#ifdef OPENMC_MPI
|
||||
MPI_Fint fcomm = MPI_Comm_c2f(comm);
|
||||
openmc_init_f(&fcomm);
|
||||
#else
|
||||
openmc_init_f(nullptr);
|
||||
#endif
|
||||
openmc_init_f();
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -79,16 +76,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
|
||||
|
|
@ -172,13 +171,13 @@ parse_command_line(int argc, char* argv[])
|
|||
|
||||
#ifdef _OPENMP
|
||||
// Read and set number of OpenMP threads
|
||||
openmc_n_threads = std::stoi(argv[i]);
|
||||
if (openmc_n_threads < 1) {
|
||||
simulation::n_threads = std::stoi(argv[i]);
|
||||
if (simulation::n_threads < 1) {
|
||||
std::string msg {"Number of threads must be positive."};
|
||||
strcpy(openmc_err_msg, msg.c_str());
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
omp_set_num_threads(openmc_n_threads);
|
||||
omp_set_num_threads(simulation::n_threads);
|
||||
#else
|
||||
if (openmc_master)
|
||||
warning("Ignoring number of threads specified on command line.");
|
||||
|
|
|
|||
|
|
@ -223,7 +223,7 @@ contains
|
|||
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
! Preallocate space for keff and entropy by generation
|
||||
call k_generation % reserve(n_max_batches*gen_per_batch)
|
||||
call k_generation_reserve(n_max_batches*gen_per_batch)
|
||||
end if
|
||||
|
||||
! Particle tracks
|
||||
|
|
|
|||
|
|
@ -2,28 +2,12 @@ module message_passing
|
|||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
#ifdef OPENMC_MPIF08
|
||||
use mpi_f08
|
||||
#else
|
||||
use mpi
|
||||
#endif
|
||||
#endif
|
||||
|
||||
! The defaults set here for the number of processors, rank, and master and
|
||||
! mpi_enabled flag are for when MPI is not being used at all, i.e. a serial
|
||||
! run. In this case, these variables are still used at times.
|
||||
! The defaults set here for the number of processors, rank, and master and are
|
||||
! for when MPI is not being used at all, i.e. a serial run. In this case, these
|
||||
! variables are still used at times.
|
||||
|
||||
integer(C_INT), bind(C, name='openmc_n_procs') :: n_procs = 1 ! number of processes
|
||||
integer(C_INT), bind(C, name='openmc_rank') :: rank = 0 ! rank of process
|
||||
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
|
||||
|
||||
end module message_passing
|
||||
|
|
|
|||
|
|
@ -13,4 +13,23 @@ MPI_Datatype bank;
|
|||
#endif
|
||||
|
||||
} // namespace mpi
|
||||
|
||||
//==============================================================================
|
||||
// Fortran compatibility functions
|
||||
//==============================================================================
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
extern "C" void
|
||||
send_int(void* buffer, int count, int dest, int tag)
|
||||
{
|
||||
MPI_Send(buffer, count, MPI_INTEGER, dest, tag, mpi::intracomm);
|
||||
}
|
||||
|
||||
extern "C" void
|
||||
recv_int(void* buffer, int count, int source, int tag)
|
||||
{
|
||||
MPI_Recv(buffer, count, MPI_INTEGER, source, tag, mpi::intracomm, MPI_STATUS_IGNORE);
|
||||
}
|
||||
#endif
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -339,7 +339,7 @@ contains
|
|||
! write out information about batch and generation
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
|
||||
trim(to_str(current_batch)) // "/" // trim(to_str(current_gen))
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation % data(i)
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation(i)
|
||||
|
||||
! write out entropy info
|
||||
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
|
||||
|
|
@ -373,7 +373,7 @@ contains
|
|||
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
|
||||
trim(to_str(current_batch)) // "/" // trim(to_str(gen_per_batch))
|
||||
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') &
|
||||
k_generation % data(i)
|
||||
k_generation(i)
|
||||
|
||||
! write out entropy info
|
||||
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
|
||||
|
|
@ -511,9 +511,9 @@ contains
|
|||
end if
|
||||
write(ou,100) " Time in active batches", time_active % elapsed
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
write(ou,100) " Time synchronizing fission bank", time_bank % elapsed
|
||||
write(ou,100) " Sampling source sites", time_bank_sample % elapsed
|
||||
write(ou,100) " SEND/RECV source sites", time_bank_sendrecv % elapsed
|
||||
write(ou,100) " Time synchronizing fission bank", time_bank_elapsed()
|
||||
write(ou,100) " Sampling source sites", time_bank_sample_elapsed()
|
||||
write(ou,100) " SEND/RECV source sites", time_bank_sendrecv_elapsed()
|
||||
end if
|
||||
write(ou,100) " Time accumulating tallies", time_tallies % elapsed
|
||||
if (cmfd_run) write(ou,100) " Time in CMFD", time_cmfd % elapsed
|
||||
|
|
|
|||
|
|
@ -127,15 +127,15 @@ Particle::mark_as_lost(const char* message)
|
|||
// Increment number of lost particles
|
||||
alive = false;
|
||||
#pragma omp atomic
|
||||
openmc_n_lost_particles += 1;
|
||||
simulation::n_lost_particles += 1;
|
||||
|
||||
// Count the total number of simulated particles (on this processor)
|
||||
auto n = openmc_current_batch * settings::gen_per_batch * openmc_work;
|
||||
auto n = simulation::current_batch * settings::gen_per_batch * simulation::work;
|
||||
|
||||
// Abort the simulation if the maximum number of lost particles has been
|
||||
// reached
|
||||
if (openmc_n_lost_particles >= MAX_LOST_PARTICLES &&
|
||||
openmc_n_lost_particles >= REL_MAX_LOST_PARTICLES*n) {
|
||||
if (simulation::n_lost_particles >= MAX_LOST_PARTICLES &&
|
||||
simulation::n_lost_particles >= REL_MAX_LOST_PARTICLES*n) {
|
||||
fatal_error("Maximum number of lost particles has been reached.");
|
||||
}
|
||||
}
|
||||
|
|
@ -148,7 +148,7 @@ Particle::write_restart() const
|
|||
|
||||
// Set up file name
|
||||
std::stringstream filename;
|
||||
filename << settings::path_output << "particle_" << openmc_current_batch
|
||||
filename << settings::path_output << "particle_" << simulation::current_batch
|
||||
<< '_' << id << ".h5";
|
||||
|
||||
#pragma omp critical (WriteParticleRestart)
|
||||
|
|
@ -165,9 +165,9 @@ Particle::write_restart() const
|
|||
#endif
|
||||
|
||||
// Write data to file
|
||||
write_dataset(file_id, "current_batch", openmc_current_batch);
|
||||
write_dataset(file_id, "current_batch", simulation::current_batch);
|
||||
write_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
|
||||
write_dataset(file_id, "current_generation", openmc_current_gen);
|
||||
write_dataset(file_id, "current_generation", simulation::current_gen);
|
||||
write_dataset(file_id, "n_particles", settings::n_particles);
|
||||
switch (settings::run_mode) {
|
||||
case RUN_MODE_FIXEDSOURCE:
|
||||
|
|
@ -188,7 +188,7 @@ Particle::write_restart() const
|
|||
int64_t n;
|
||||
openmc_source_bank(&src, &n);
|
||||
|
||||
int64_t i = openmc_current_work;
|
||||
int64_t i = simulation::current_work;
|
||||
write_dataset(file_id, "weight", src[i-1].wgt);
|
||||
write_dataset(file_id, "energy", src[i-1].E);
|
||||
hsize_t dims[] {3};
|
||||
|
|
|
|||
|
|
@ -16,6 +16,7 @@
|
|||
#include "openmc/random_lcg.h"
|
||||
#include "openmc/settings.h"
|
||||
#include "openmc/simulation.h"
|
||||
#include "openmc/tallies/tally.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -30,7 +31,7 @@ collision_mg(Particle* p, const double* energy_bin_avg,
|
|||
sample_reaction(p, energy_bin_avg, material_xs);
|
||||
|
||||
// Display information about collision
|
||||
if ((settings::verbosity >= 10) || (openmc_trace)) {
|
||||
if ((settings::verbosity >= 10) || (simulation::trace)) {
|
||||
std::stringstream msg;
|
||||
msg << " Energy Group = " << p->g;
|
||||
write_message(msg, 1);
|
||||
|
|
@ -115,7 +116,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
|
|||
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
|
||||
|
||||
// Determine the expected number of neutrons produced
|
||||
double nu_t = p->wgt / openmc_keff * weight * material_xs->nu_fission /
|
||||
double nu_t = p->wgt / simulation::keff * weight * material_xs->nu_fission /
|
||||
material_xs->total;
|
||||
|
||||
// Sample the number of neutrons produced
|
||||
|
|
|
|||
|
|
@ -17,6 +17,7 @@
|
|||
#include "openmc/mesh.h"
|
||||
#include "openmc/output.h"
|
||||
#include "openmc/random_lcg.h"
|
||||
#include "openmc/simulation.h"
|
||||
#include "openmc/source.h"
|
||||
#include "openmc/string_utils.h"
|
||||
#include "openmc/xml_interface.h"
|
||||
|
|
@ -57,7 +58,7 @@ bool urr_ptables_on {true};
|
|||
bool write_all_tracks {false};
|
||||
bool write_initial_source {false};
|
||||
bool dagmc {false};
|
||||
|
||||
|
||||
std::string path_cross_sections;
|
||||
std::string path_input;
|
||||
std::string path_multipole;
|
||||
|
|
@ -218,7 +219,7 @@ void read_settings_xml()
|
|||
fatal_error("DAGMC mode unsupported for this build of OpenMC");
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
// To this point, we haven't displayed any output since we didn't know what
|
||||
// the verbosity is. Now that we checked for it, show the title if necessary
|
||||
if (openmc_master) {
|
||||
|
|
@ -393,14 +394,14 @@ void read_settings_xml()
|
|||
// Number of OpenMP threads
|
||||
if (check_for_node(root, "threads")) {
|
||||
#ifdef _OPENMP
|
||||
if (openmc_n_threads == 0) {
|
||||
openmc_n_threads = std::stoi(get_node_value(root, "threads"));
|
||||
if (openmc_n_threads < 1) {
|
||||
if (simulation::n_threads == 0) {
|
||||
simulation::n_threads = std::stoi(get_node_value(root, "threads"));
|
||||
if (simulation::n_threads < 1) {
|
||||
std::stringstream msg;
|
||||
msg << "Invalid number of threads: " << openmc_n_threads;
|
||||
msg << "Invalid number of threads: " << simulation::n_threads;
|
||||
fatal_error(msg);
|
||||
}
|
||||
omp_set_num_threads(openmc_n_threads);
|
||||
omp_set_num_threads(simulation::n_threads);
|
||||
}
|
||||
#else
|
||||
if (openmc_master) warning("OpenMC was not compiled with OpenMP support; "
|
||||
|
|
@ -414,7 +415,7 @@ void read_settings_xml()
|
|||
omp_set_num_threads(1);
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
// ==========================================================================
|
||||
// EXTERNAL SOURCE
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
k_sum
|
||||
#ifdef _OPENMP
|
||||
use eigenvalue, only: join_bank_from_threads
|
||||
#endif
|
||||
|
|
@ -57,6 +57,9 @@ module simulation
|
|||
subroutine initialize_source() bind(C)
|
||||
end subroutine
|
||||
|
||||
subroutine initialize_generation() bind(C)
|
||||
end subroutine
|
||||
|
||||
function sample_external_source() result(site) bind(C)
|
||||
import Bank
|
||||
type(Bank) :: site
|
||||
|
|
@ -222,30 +225,6 @@ contains
|
|||
|
||||
end subroutine initialize_batch
|
||||
|
||||
!===============================================================================
|
||||
! INITIALIZE_GENERATION
|
||||
!===============================================================================
|
||||
|
||||
subroutine initialize_generation()
|
||||
|
||||
interface
|
||||
subroutine ufs_count_sites() bind(C)
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
! Reset number of fission bank sites
|
||||
n_bank = 0
|
||||
|
||||
! Count source sites if using uniform fission source weighting
|
||||
if (ufs) call ufs_count_sites()
|
||||
|
||||
! Store current value of tracklength k
|
||||
keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH)
|
||||
end if
|
||||
|
||||
end subroutine initialize_generation
|
||||
|
||||
!===============================================================================
|
||||
! FINALIZE_GENERATION
|
||||
!===============================================================================
|
||||
|
|
@ -258,6 +237,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
|
||||
|
|
@ -291,9 +273,7 @@ contains
|
|||
#endif
|
||||
|
||||
! Distribute fission bank across processors evenly
|
||||
call time_bank % start()
|
||||
call synchronize_bank()
|
||||
call time_bank % stop()
|
||||
|
||||
! Calculate shannon entropy
|
||||
if (entropy_on) call shannon_entropy()
|
||||
|
|
@ -325,11 +305,13 @@ contains
|
|||
subroutine finalize_batch()
|
||||
|
||||
integer(C_INT) :: err
|
||||
#ifdef OPENMC_MPI
|
||||
integer :: mpi_err ! MPI error code
|
||||
#endif
|
||||
character(MAX_FILE_LEN) :: filename
|
||||
|
||||
interface
|
||||
subroutine broadcast_triggers() bind(C)
|
||||
end subroutine broadcast_triggers
|
||||
end interface
|
||||
|
||||
! Reduce tallies onto master process and accumulate
|
||||
call time_tallies % start()
|
||||
call accumulate_tallies()
|
||||
|
|
@ -351,8 +333,7 @@ contains
|
|||
! Check_triggers
|
||||
if (master) call check_triggers()
|
||||
#ifdef OPENMC_MPI
|
||||
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
|
||||
mpi_intracomm, mpi_err)
|
||||
call broadcast_triggers()
|
||||
#endif
|
||||
if (satisfy_triggers .or. &
|
||||
(trigger_on .and. current_batch == n_max_batches)) then
|
||||
|
|
@ -402,9 +383,6 @@ contains
|
|||
! Set up tally procedure pointers
|
||||
call init_tally_routines()
|
||||
|
||||
! Determine how much work each processor should do
|
||||
call calculate_work()
|
||||
|
||||
! Allocate source bank, and for eigenvalue simulations also allocate the
|
||||
! fission bank
|
||||
call allocate_banks()
|
||||
|
|
@ -436,7 +414,7 @@ contains
|
|||
! Reset global variables -- this is done before loading state point (as that
|
||||
! will potentially populate k_generation and entropy)
|
||||
current_batch = 0
|
||||
call k_generation % clear()
|
||||
call k_generation_clear()
|
||||
call entropy_clear()
|
||||
need_depletion_rx = .false.
|
||||
|
||||
|
|
@ -473,21 +451,13 @@ contains
|
|||
integer(C_INT) :: err
|
||||
|
||||
integer :: i ! loop index
|
||||
#ifdef OPENMC_MPI
|
||||
integer :: n ! size of arrays
|
||||
integer :: mpi_err ! MPI error code
|
||||
integer :: count_per_filter ! number of result values for one filter bin
|
||||
real(8) :: tempr(3) ! temporary array for communication
|
||||
#ifdef OPENMC_MPIF08
|
||||
type(MPI_Datatype) :: result_block
|
||||
#else
|
||||
integer :: result_block
|
||||
#endif
|
||||
#endif
|
||||
|
||||
interface
|
||||
subroutine print_overlap_check() bind(C)
|
||||
end subroutine print_overlap_check
|
||||
|
||||
subroutine broadcast_results() bind(C)
|
||||
end subroutine broadcast_results
|
||||
end interface
|
||||
|
||||
err = 0
|
||||
|
|
@ -515,37 +485,7 @@ contains
|
|||
total_gen = total_gen + current_batch*gen_per_batch
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
! Broadcast tally results so that each process has access to results
|
||||
if (allocated(tallies)) then
|
||||
do i = 1, size(tallies)
|
||||
associate (results => tallies(i) % obj % results)
|
||||
! Create a new datatype that consists of all values for a given filter
|
||||
! bin and then use that to broadcast. This is done to minimize the
|
||||
! chance of the 'count' argument of MPI_BCAST exceeding 2**31
|
||||
n = size(results, 3)
|
||||
count_per_filter = size(results, 1) * size(results, 2)
|
||||
call MPI_TYPE_CONTIGUOUS(count_per_filter, MPI_DOUBLE, &
|
||||
result_block, mpi_err)
|
||||
call MPI_TYPE_COMMIT(result_block, mpi_err)
|
||||
call MPI_BCAST(results, n, result_block, 0, mpi_intracomm, mpi_err)
|
||||
call MPI_TYPE_FREE(result_block, mpi_err)
|
||||
end associate
|
||||
end do
|
||||
end if
|
||||
|
||||
! Also broadcast global tally results
|
||||
n = size(global_tallies)
|
||||
call MPI_BCAST(global_tallies, n, MPI_DOUBLE, 0, mpi_intracomm, mpi_err)
|
||||
|
||||
! These guys are needed so that non-master processes can calculate the
|
||||
! combined estimate of k-effective
|
||||
tempr(1) = k_col_abs
|
||||
tempr(2) = k_col_tra
|
||||
tempr(3) = k_abs_tra
|
||||
call MPI_BCAST(tempr, 3, MPI_REAL8, 0, mpi_intracomm, mpi_err)
|
||||
k_col_abs = tempr(1)
|
||||
k_col_tra = tempr(2)
|
||||
k_abs_tra = tempr(3)
|
||||
call broadcast_results()
|
||||
#endif
|
||||
|
||||
! Write tally results to tallies.out
|
||||
|
|
@ -573,46 +513,6 @@ contains
|
|||
|
||||
end function openmc_simulation_finalize
|
||||
|
||||
!===============================================================================
|
||||
! CALCULATE_WORK determines how many particles each processor should simulate
|
||||
!===============================================================================
|
||||
|
||||
subroutine calculate_work()
|
||||
|
||||
integer :: i ! loop index
|
||||
integer :: remainder ! Number of processors with one extra particle
|
||||
integer(8) :: i_bank ! Running count of number of particles
|
||||
integer(8) :: min_work ! Minimum number of particles on each proc
|
||||
integer(8) :: work_i ! Number of particles on rank i
|
||||
|
||||
if (.not. allocated(work_index)) allocate(work_index(0:n_procs))
|
||||
|
||||
! Determine minimum amount of particles to simulate on each processor
|
||||
min_work = n_particles/n_procs
|
||||
|
||||
! Determine number of processors that have one extra particle
|
||||
remainder = int(mod(n_particles, int(n_procs,8)), 4)
|
||||
|
||||
i_bank = 0
|
||||
work_index(0) = 0
|
||||
do i = 0, n_procs - 1
|
||||
! Number of particles for rank i
|
||||
if (i < remainder) then
|
||||
work_i = min_work + 1
|
||||
else
|
||||
work_i = min_work
|
||||
end if
|
||||
|
||||
! Set number of particles
|
||||
if (rank == i) work = work_i
|
||||
|
||||
! Set index into source bank for rank i
|
||||
i_bank = i_bank + work_i
|
||||
work_index(i+1) = i_bank
|
||||
end do
|
||||
|
||||
end subroutine calculate_work
|
||||
|
||||
!===============================================================================
|
||||
! ALLOCATE_BANKS allocates memory for the fission and source banks
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -1,8 +1,12 @@
|
|||
#include "openmc/simulation.h"
|
||||
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/eigenvalue.h"
|
||||
#include "openmc/message_passing.h"
|
||||
#include "openmc/settings.h"
|
||||
#include "openmc/tallies/tally.h"
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
// OPENMC_RUN encompasses all the main logic where iterations are performed
|
||||
// over the batches, generations, and histories in a fixed source or k-eigenvalue
|
||||
|
|
@ -27,10 +31,39 @@ namespace openmc {
|
|||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
namespace simulation {
|
||||
|
||||
int current_batch;
|
||||
int current_gen;
|
||||
int64_t current_work;
|
||||
double keff {1.0};
|
||||
double keff_std;
|
||||
double k_col_abs {0.0};
|
||||
double k_col_tra {0.0};
|
||||
double k_abs_tra {0.0};
|
||||
double log_spacing;
|
||||
int n_lost_particles {0};
|
||||
bool need_depletion_rx {false};
|
||||
int restart_batch;
|
||||
bool satisfy_triggers {false};
|
||||
bool simulation_initialized {false};
|
||||
int total_gen {0};
|
||||
int64_t work;
|
||||
|
||||
std::vector<double> k_generation;
|
||||
std::vector<int64_t> work_index;
|
||||
|
||||
// Threadprivate variables
|
||||
bool trace; //!< flag to show debug information
|
||||
#ifdef _OPENMP
|
||||
int n_threads {-1}; //!< number of OpenMP threads
|
||||
int thread_id; //!< ID of a given thread
|
||||
#endif
|
||||
|
||||
} // namespace simulation
|
||||
|
||||
//==============================================================================
|
||||
// Functions
|
||||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
void openmc_simulation_init_c()
|
||||
|
|
@ -39,6 +72,26 @@ void openmc_simulation_init_c()
|
|||
calculate_work();
|
||||
}
|
||||
|
||||
void initialize_generation()
|
||||
{
|
||||
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
|
||||
// Reset number of fission bank sites
|
||||
n_bank = 0;
|
||||
|
||||
// Count source sites if using uniform fission source weighting
|
||||
if (settings::ufs_on) ufs_count_sites();
|
||||
|
||||
// Store current value of tracklength k
|
||||
keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE);
|
||||
}
|
||||
}
|
||||
|
||||
int overall_generation()
|
||||
{
|
||||
using namespace simulation;
|
||||
return settings::gen_per_batch*(current_batch - 1) + current_gen;
|
||||
}
|
||||
|
||||
void calculate_work()
|
||||
{
|
||||
// Determine minimum amount of particles to simulate on each processor
|
||||
|
|
@ -48,19 +101,67 @@ void calculate_work()
|
|||
int64_t remainder = settings::n_particles % mpi::n_procs;
|
||||
|
||||
int64_t i_bank = 0;
|
||||
work_index.reserve(mpi::n_procs);
|
||||
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;
|
||||
|
||||
// Set number of particles
|
||||
if (mpi::rank == i) openmc_work = work_i;
|
||||
if (mpi::rank == i) simulation::work = work_i;
|
||||
|
||||
// Set index into source bank for rank i
|
||||
i_bank += work_i;
|
||||
work_index.push_back(i_bank);
|
||||
simulation::work_index[i + 1] = i_bank;
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
void broadcast_results() {
|
||||
// Broadcast tally results so that each process has access to results
|
||||
for (int i = 1; i <= n_tallies; ++i) {
|
||||
// Create a new datatype that consists of all values for a given filter
|
||||
// bin and then use that to broadcast. This is done to minimize the
|
||||
// chance of the 'count' argument of MPI_BCAST exceeding 2**31
|
||||
auto results = tally_results(i);
|
||||
|
||||
auto shape = results.shape();
|
||||
int count_per_filter = shape[1] * shape[2];
|
||||
MPI_Datatype result_block;
|
||||
MPI_Type_contiguous(count_per_filter, MPI_DOUBLE, &result_block);
|
||||
MPI_Type_commit(&result_block);
|
||||
MPI_Bcast(results.data(), shape[0], result_block, 0, mpi::intracomm);
|
||||
MPI_Type_free(&result_block);
|
||||
}
|
||||
|
||||
// Also broadcast global tally results
|
||||
auto gt = global_tallies();
|
||||
MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm);
|
||||
|
||||
// These guys are needed so that non-master processes can calculate the
|
||||
// combined estimate of k-effective
|
||||
double temp[] {simulation::k_col_abs, simulation::k_col_tra,
|
||||
simulation::k_abs_tra};
|
||||
MPI_Bcast(temp, 3, MPI_DOUBLE, 0, mpi::intracomm);
|
||||
simulation::k_col_abs = temp[0];
|
||||
simulation::k_col_tra = temp[1];
|
||||
simulation::k_abs_tra = temp[2];
|
||||
}
|
||||
|
||||
void broadcast_triggers()
|
||||
{
|
||||
MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
|
||||
}
|
||||
#endif
|
||||
|
||||
//==============================================================================
|
||||
// Fortran compatibility
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double k_generation(int i) { return simulation::k_generation.at(i - 1); }
|
||||
extern "C" int k_generation_size() { return simulation::k_generation.size(); }
|
||||
extern "C" void k_generation_clear() { simulation::k_generation.clear(); }
|
||||
extern "C" void k_generation_reserve(int i) { simulation::k_generation.reserve(i); }
|
||||
extern "C" int64_t work_index(int rank) { return simulation::work_index[rank]; }
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -13,84 +13,98 @@ module simulation_header
|
|||
! GEOMETRY-RELATED VARIABLES
|
||||
|
||||
! Number of lost particles
|
||||
integer(C_INT), bind(C, name='openmc_n_lost_particles') :: n_lost_particles = 0
|
||||
integer(C_INT), bind(C) :: n_lost_particles
|
||||
|
||||
real(8) :: log_spacing ! spacing on logarithmic grid
|
||||
real(C_DOUBLE), bind(C) :: log_spacing ! spacing on logarithmic grid
|
||||
|
||||
! ============================================================================
|
||||
! SIMULATION VARIABLES
|
||||
|
||||
integer(C_INT), bind(C, name='openmc_current_batch') :: current_batch ! current batch
|
||||
integer(C_INT), bind(C, name='openmc_current_gen') :: current_gen ! current generation within a batch
|
||||
integer(C_INT), bind(C, name='openmc_total_gen') :: total_gen = 0 ! total number of generations simulated
|
||||
logical(C_BOOL), bind(C, name='openmc_simulation_initialized') :: &
|
||||
simulation_initialized = .false.
|
||||
logical :: need_depletion_rx ! need to calculate depletion reaction rx?
|
||||
integer(C_INT), bind(C) :: current_batch ! current batch
|
||||
integer(C_INT), bind(C) :: current_gen ! current generation within a batch
|
||||
integer(C_INT), bind(C) :: total_gen ! total number of generations simulated
|
||||
logical(C_BOOL), bind(C) :: simulation_initialized
|
||||
logical(C_BOOL), bind(C) :: need_depletion_rx ! need to calculate depletion reaction rx?
|
||||
|
||||
! ============================================================================
|
||||
! TALLY PRECISION TRIGGER VARIABLES
|
||||
|
||||
logical :: satisfy_triggers = .false. ! whether triggers are satisfied
|
||||
logical(C_BOOL), bind(C) :: satisfy_triggers ! whether triggers are satisfied
|
||||
|
||||
integer(C_INT64_T), bind(C, name='openmc_work') :: work ! number of particles per processor
|
||||
integer(C_INT64_T), allocatable :: work_index(:) ! starting index in source bank for each process
|
||||
integer(C_INT64_T), bind(C, name='openmc_current_work') :: current_work ! index in source bank of current history simulated
|
||||
integer(C_INT64_T), bind(C) :: work ! number of particles per processor
|
||||
integer(C_INT64_T), bind(C) :: current_work ! index in source bank of current history simulated
|
||||
|
||||
! ============================================================================
|
||||
! K-EIGENVALUE SIMULATION VARIABLES
|
||||
|
||||
! Temporary k-effective values
|
||||
type(VectorReal) :: k_generation ! single-generation estimates of k
|
||||
real(C_DOUBLE), bind(C, name='openmc_keff') :: keff = ONE ! average k over active batches
|
||||
real(C_DOUBLE), bind(C, name='openmc_keff_std') :: keff_std ! standard deviation of average k
|
||||
real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption
|
||||
real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength
|
||||
real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength
|
||||
real(C_DOUBLE), bind(C) :: keff ! average k over active batches
|
||||
real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k
|
||||
real(C_DOUBLE), bind(C) :: k_col_abs ! sum over batches of k_collision * k_absorption
|
||||
real(C_DOUBLE), bind(C) :: k_col_tra ! sum over batches of k_collision * k_tracklength
|
||||
real(C_DOUBLE), bind(C) :: k_abs_tra ! sum over batches of k_absorption * k_tracklength
|
||||
|
||||
! ============================================================================
|
||||
! PARALLEL PROCESSING VARIABLES
|
||||
|
||||
#ifdef _OPENMP
|
||||
integer(C_INT), bind(C, name='openmc_n_threads') :: n_threads = NONE ! number of OpenMP threads
|
||||
integer :: thread_id ! ID of a given thread
|
||||
integer(C_INT), bind(C) :: n_threads ! number of OpenMP threads
|
||||
integer(C_INT), bind(C) :: thread_id ! ID of a given thread
|
||||
#endif
|
||||
|
||||
! ============================================================================
|
||||
! MISCELLANEOUS VARIABLES
|
||||
|
||||
integer :: restart_batch
|
||||
integer(C_INT), bind(C) :: restart_batch
|
||||
|
||||
logical(C_BOOL), bind(C, name='openmc_trace') :: trace
|
||||
logical(C_BOOL), bind(C) :: trace
|
||||
|
||||
!$omp threadprivate(trace, thread_id, current_work)
|
||||
|
||||
interface
|
||||
subroutine entropy_clear() bind(C)
|
||||
end subroutine
|
||||
|
||||
pure function overall_generation() result(gen) bind(C)
|
||||
import C_INT
|
||||
integer(C_INT) :: gen
|
||||
end function overall_generation
|
||||
|
||||
function k_generation(i) result(k) bind(C)
|
||||
import C_DOUBLE, C_INT
|
||||
integer(C_INT), value :: i
|
||||
real(C_DOUBLE) :: k
|
||||
end function
|
||||
|
||||
function k_generation_size() result(sz) bind(C)
|
||||
import C_INT
|
||||
integer(C_INT) :: sz
|
||||
end function
|
||||
|
||||
subroutine k_generation_clear() bind(C)
|
||||
end subroutine
|
||||
|
||||
subroutine k_generation_reserve(i) bind(C)
|
||||
import C_INT
|
||||
integer(C_INT), value :: i
|
||||
end subroutine
|
||||
|
||||
function work_index(rank) result(i) bind(C)
|
||||
import C_INT, C_INT64_T
|
||||
integer(C_INT), value :: rank
|
||||
integer(C_INT64_T) :: i
|
||||
end function
|
||||
end interface
|
||||
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! OVERALL_GENERATION determines the overall generation number
|
||||
!===============================================================================
|
||||
|
||||
pure function overall_generation() result(gen) bind(C)
|
||||
integer(C_INT) :: gen
|
||||
gen = gen_per_batch*(current_batch - 1) + current_gen
|
||||
end function overall_generation
|
||||
|
||||
!===============================================================================
|
||||
! FREE_MEMORY_SIMULATION deallocates global arrays defined in this module
|
||||
!===============================================================================
|
||||
|
||||
subroutine free_memory_simulation()
|
||||
|
||||
if (allocated(work_index)) deallocate(work_index)
|
||||
|
||||
call k_generation % clear()
|
||||
call k_generation % shrink_to_fit()
|
||||
call k_generation_clear()
|
||||
call entropy_clear()
|
||||
end subroutine free_memory_simulation
|
||||
|
||||
|
|
|
|||
|
|
@ -264,17 +264,17 @@ void initialize_source()
|
|||
}
|
||||
|
||||
// Read in the source bank
|
||||
read_source_bank(file_id, work_index.data(), source_bank);
|
||||
read_source_bank(file_id, source_bank);
|
||||
|
||||
// Close file
|
||||
file_close(file_id);
|
||||
|
||||
} else {
|
||||
// Generation source sites from specified distribution in user input
|
||||
for (int64_t i = 0; i < openmc_work; ++i) {
|
||||
for (int64_t i = 0; i < simulation::work; ++i) {
|
||||
// initialize random number seed
|
||||
int64_t id = openmc_total_gen*settings::n_particles +
|
||||
work_index[openmc::mpi::rank] + i + 1;
|
||||
int64_t id = simulation::total_gen*settings::n_particles +
|
||||
simulation::work_index[mpi::rank] + i + 1;
|
||||
set_particle_seed(id);
|
||||
|
||||
// sample external source distribution
|
||||
|
|
@ -287,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, work_index.data(), source_bank);
|
||||
write_source_bank(file_id, source_bank);
|
||||
file_close(file_id);
|
||||
}
|
||||
}
|
||||
|
|
@ -364,10 +364,10 @@ extern "C" void fill_source_bank_fixedsource()
|
|||
int64_t n;
|
||||
openmc_source_bank(&source_bank, &n);
|
||||
|
||||
for (int64_t i = 0; i < openmc_work; ++i) {
|
||||
for (int64_t i = 0; i < simulation::work; ++i) {
|
||||
// initialize random number seed
|
||||
int64_t id = (openmc_total_gen + overall_generation()) *
|
||||
settings::n_particles + work_index[openmc::mpi::rank] + i + 1;
|
||||
int64_t id = (simulation::total_gen + overall_generation()) *
|
||||
settings::n_particles + simulation::work_index[mpi::rank] + i + 1;
|
||||
set_particle_seed(id);
|
||||
|
||||
// sample external source distribution
|
||||
|
|
|
|||
|
|
@ -36,17 +36,15 @@ module state_point
|
|||
implicit none
|
||||
|
||||
interface
|
||||
subroutine write_source_bank(group_id, work_index, bank_) bind(C)
|
||||
subroutine write_source_bank(group_id, bank_) bind(C)
|
||||
import HID_T, C_INT64_T, Bank
|
||||
integer(HID_T), value :: group_id
|
||||
integer(C_INT64_T), intent(in) :: work_index(*)
|
||||
type(Bank), intent(in) :: bank_(*)
|
||||
end subroutine write_source_bank
|
||||
|
||||
subroutine read_source_bank(group_id, work_index, bank_) bind(C)
|
||||
subroutine read_source_bank(group_id, bank_) bind(C)
|
||||
import HID_T, C_INT64_T, Bank
|
||||
integer(HID_T), value :: group_id
|
||||
integer(C_INT64_T), intent(in) :: work_index(*)
|
||||
type(Bank), intent(out) :: bank_(*)
|
||||
end subroutine read_source_bank
|
||||
end interface
|
||||
|
|
@ -63,15 +61,13 @@ contains
|
|||
integer(C_INT) :: err
|
||||
|
||||
logical :: write_source_
|
||||
integer :: i, j, k
|
||||
integer :: i, j
|
||||
integer :: i_xs
|
||||
integer, allocatable :: id_array(:)
|
||||
integer(HID_T) :: file_id
|
||||
integer(HID_T) :: cmfd_group, tallies_group, tally_group, &
|
||||
filters_group, filter_group, derivs_group, &
|
||||
deriv_group, runtime_group
|
||||
integer(C_INT) :: ignored_err
|
||||
real(C_DOUBLE) :: k_combined(2)
|
||||
character(MAX_WORD_LEN), allocatable :: str_array(:)
|
||||
character(C_CHAR), pointer :: string(:)
|
||||
character(len=:, kind=C_CHAR), allocatable :: filename_
|
||||
|
|
@ -83,10 +79,14 @@ contains
|
|||
import HID_T
|
||||
integer(HID_T), value :: group
|
||||
end subroutine
|
||||
subroutine entropy_to_hdf5(group) bind(C)
|
||||
subroutine write_eigenvalue_hdf5(group) bind(C)
|
||||
import HID_T
|
||||
integer(HID_T), value :: group
|
||||
end subroutine
|
||||
subroutine write_tally_results_nr(file_id) bind(C)
|
||||
import HID_T
|
||||
integer(HID_T), value :: file_id
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
err = 0
|
||||
|
|
@ -169,16 +169,7 @@ contains
|
|||
|
||||
! Write out information for eigenvalue run
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
call write_dataset(file_id, "n_inactive", n_inactive)
|
||||
call write_dataset(file_id, "generations_per_batch", gen_per_batch)
|
||||
k = k_generation % size()
|
||||
call write_dataset(file_id, "k_generation", k_generation % data(1:k))
|
||||
call entropy_to_hdf5(file_id)
|
||||
call write_dataset(file_id, "k_col_abs", k_col_abs)
|
||||
call write_dataset(file_id, "k_col_tra", k_col_tra)
|
||||
call write_dataset(file_id, "k_abs_tra", k_abs_tra)
|
||||
ignored_err = openmc_get_keff(k_combined)
|
||||
call write_dataset(file_id, "k_combined", k_combined)
|
||||
call write_eigenvalue_hdf5(file_id)
|
||||
|
||||
! Write out CMFD info
|
||||
if (cmfd_on) then
|
||||
|
|
@ -415,11 +406,11 @@ contains
|
|||
time_active % get_value())
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
call write_dataset(runtime_group, "synchronizing fission bank", &
|
||||
time_bank % get_value())
|
||||
time_bank_elapsed())
|
||||
call write_dataset(runtime_group, "sampling source sites", &
|
||||
time_bank_sample % get_value())
|
||||
time_bank_sample_elapsed())
|
||||
call write_dataset(runtime_group, "SEND-RECV source sites", &
|
||||
time_bank_sendrecv % get_value())
|
||||
time_bank_sendrecv_elapsed())
|
||||
end if
|
||||
call write_dataset(runtime_group, "accumulating tallies", &
|
||||
time_tallies % get_value())
|
||||
|
|
@ -447,7 +438,7 @@ contains
|
|||
if (master .or. parallel) then
|
||||
file_id = file_open(filename_, 'a', parallel=.true.)
|
||||
end if
|
||||
call write_source_bank(file_id, work_index, source_bank)
|
||||
call write_source_bank(file_id, source_bank)
|
||||
if (master .or. parallel) call file_close(file_id)
|
||||
end if
|
||||
end function openmc_statepoint_write
|
||||
|
|
@ -485,134 +476,11 @@ contains
|
|||
file_id = file_open(filename_, 'w', parallel=.true.)
|
||||
call write_attribute(file_id, "filetype", 'source')
|
||||
end if
|
||||
call write_source_bank(file_id, work_index, source_bank)
|
||||
call write_source_bank(file_id, source_bank)
|
||||
if (master .or. parallel) call file_close(file_id)
|
||||
|
||||
end subroutine write_source_point
|
||||
|
||||
!===============================================================================
|
||||
! WRITE_TALLY_RESULTS_NR
|
||||
!===============================================================================
|
||||
|
||||
subroutine write_tally_results_nr(file_id)
|
||||
integer(HID_T), intent(in) :: file_id
|
||||
|
||||
integer :: i ! loop index
|
||||
integer :: n ! number of filter bins
|
||||
integer :: m ! number of score bins
|
||||
integer :: n_bins ! total number of bins
|
||||
integer(HID_T) :: tallies_group, tally_group
|
||||
real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results
|
||||
real(8), target :: global_temp(3,N_GLOBAL_TALLIES)
|
||||
#ifdef OPENMC_MPI
|
||||
integer :: mpi_err ! MPI error code
|
||||
real(8) :: dummy ! temporary receive buffer for non-root reduces
|
||||
#endif
|
||||
type(TallyObject) :: dummy_tally
|
||||
|
||||
! ==========================================================================
|
||||
! COLLECT AND WRITE GLOBAL TALLIES
|
||||
|
||||
if (master) then
|
||||
! Write number of realizations
|
||||
call write_dataset(file_id, "n_realizations", n_realizations)
|
||||
|
||||
! Write number of global tallies
|
||||
call write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES)
|
||||
|
||||
tallies_group = open_group(file_id, "tallies")
|
||||
end if
|
||||
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
! Reduce global tallies
|
||||
n_bins = size(global_tallies)
|
||||
call MPI_REDUCE(global_tallies, global_temp, n_bins, MPI_REAL8, MPI_SUM, &
|
||||
0, mpi_intracomm, mpi_err)
|
||||
#endif
|
||||
|
||||
if (master) then
|
||||
! Transfer values to value on master
|
||||
if (current_batch == n_max_batches .or. satisfy_triggers) then
|
||||
global_tallies(:,:) = global_temp(:,:)
|
||||
end if
|
||||
|
||||
! Write out global tallies sum and sum_sq
|
||||
call write_dataset(file_id, "global_tallies", global_temp)
|
||||
end if
|
||||
|
||||
if (active_tallies % size() > 0) then
|
||||
! Indicate that tallies are on
|
||||
if (master) then
|
||||
call write_attribute(file_id, "tallies_present", 1)
|
||||
end if
|
||||
|
||||
! Write all tally results
|
||||
TALLY_RESULTS: do i = 1, n_tallies
|
||||
associate (t => tallies(i) % obj)
|
||||
! Determine size of tally results array
|
||||
m = size(t % results, 2)
|
||||
n = size(t % results, 3)
|
||||
n_bins = m*n*2
|
||||
|
||||
! Allocate array for storing sums and sums of squares, but
|
||||
! contiguously in memory for each
|
||||
allocate(tally_temp(2,m,n))
|
||||
tally_temp(1,:,:) = t % results(RESULT_SUM,:,:)
|
||||
tally_temp(2,:,:) = t % results(RESULT_SUM_SQ,:,:)
|
||||
|
||||
if (master) then
|
||||
tally_group = open_group(tallies_group, "tally " // &
|
||||
trim(to_str(t % id)))
|
||||
|
||||
! The MPI_IN_PLACE specifier allows the master to copy values into
|
||||
! a receive buffer without having a temporary variable
|
||||
#ifdef OPENMC_MPI
|
||||
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, &
|
||||
MPI_SUM, 0, mpi_intracomm, mpi_err)
|
||||
#endif
|
||||
|
||||
! At the end of the simulation, store the results back in the
|
||||
! regular TallyResults array
|
||||
if (current_batch == n_max_batches .or. satisfy_triggers) then
|
||||
t % results(RESULT_SUM,:,:) = tally_temp(1,:,:)
|
||||
t % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:)
|
||||
end if
|
||||
|
||||
! Put in temporary tally result
|
||||
allocate(dummy_tally % results(3,m,n))
|
||||
dummy_tally % results(RESULT_SUM,:,:) = tally_temp(1,:,:)
|
||||
dummy_tally % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:)
|
||||
|
||||
! Write reduced tally results to file
|
||||
call dummy_tally % write_results_hdf5(tally_group)
|
||||
|
||||
! Deallocate temporary tally result
|
||||
deallocate(dummy_tally % results)
|
||||
else
|
||||
! Receive buffer not significant at other processors
|
||||
#ifdef OPENMC_MPI
|
||||
call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, &
|
||||
0, mpi_intracomm, mpi_err)
|
||||
#endif
|
||||
end if
|
||||
|
||||
! Deallocate temporary copy of tally results
|
||||
deallocate(tally_temp)
|
||||
|
||||
if (master) call close_group(tally_group)
|
||||
end associate
|
||||
end do TALLY_RESULTS
|
||||
|
||||
if (master) call close_group(tallies_group)
|
||||
else
|
||||
! Indicate that tallies are off
|
||||
if (master) call write_dataset(file_id, "tallies_present", 0)
|
||||
end if
|
||||
|
||||
|
||||
end subroutine write_tally_results_nr
|
||||
|
||||
!===============================================================================
|
||||
! LOAD_STATE_POINT
|
||||
!===============================================================================
|
||||
|
|
@ -632,7 +500,9 @@ contains
|
|||
character(MAX_WORD_LEN) :: word
|
||||
|
||||
interface
|
||||
subroutine entropy_from_hdf5() bind(C)
|
||||
subroutine read_eigenvalue_hdf5(group) bind(C)
|
||||
import HID_T
|
||||
integer(HID_T), value :: group
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
|
|
@ -711,16 +581,7 @@ contains
|
|||
! Read information specific to eigenvalue run
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
call read_dataset(int_array(1), file_id, "n_inactive")
|
||||
call read_dataset(gen_per_batch, file_id, "generations_per_batch")
|
||||
|
||||
n = restart_batch*gen_per_batch
|
||||
call k_generation % resize(n)
|
||||
call read_dataset(k_generation % data(1:n), file_id, "k_generation")
|
||||
|
||||
call entropy_from_hdf5()
|
||||
call read_dataset(k_col_abs, file_id, "k_col_abs")
|
||||
call read_dataset(k_col_tra, file_id, "k_col_tra")
|
||||
call read_dataset(k_abs_tra, file_id, "k_abs_tra")
|
||||
call read_eigenvalue_hdf5(file_id)
|
||||
|
||||
! Take maximum of statepoint n_inactive and input n_inactive
|
||||
n_inactive = max(n_inactive, int_array(1))
|
||||
|
|
@ -753,13 +614,14 @@ contains
|
|||
! of active cycle or inactive cycle
|
||||
if (restart_batch > n_inactive) then
|
||||
do i = n_inactive + 1, restart_batch
|
||||
k_sum(1) = k_sum(1) + k_generation % data(i)
|
||||
k_sum(2) = k_sum(2) + k_generation % data(i)**2
|
||||
k_sum(1) = k_sum(1) + k_generation(i)
|
||||
k_sum(2) = k_sum(2) + k_generation(i)**2
|
||||
end do
|
||||
n = gen_per_batch*n_realizations
|
||||
keff = k_sum(1) / n
|
||||
else
|
||||
keff = k_generation % data(n)
|
||||
n = k_generation_size()
|
||||
keff = k_generation(n)
|
||||
end if
|
||||
current_batch = restart_batch
|
||||
|
||||
|
|
@ -820,7 +682,7 @@ contains
|
|||
end if
|
||||
|
||||
! Write out source
|
||||
call read_source_bank(file_id, work_index, source_bank)
|
||||
call read_source_bank(file_id, source_bank)
|
||||
|
||||
end if
|
||||
|
||||
|
|
|
|||
|
|
@ -1,16 +1,20 @@
|
|||
#include "openmc/state_point.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
#include "mpi.h"
|
||||
#endif
|
||||
#include "xtensor/xbuilder.hpp" // for empty_like
|
||||
#include "xtensor/xview.hpp"
|
||||
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/error.h"
|
||||
#include "openmc/hdf5_interface.h"
|
||||
#include "openmc/message_passing.h"
|
||||
#include "openmc/settings.h"
|
||||
#include "openmc/simulation.h"
|
||||
#include "openmc/tallies/tally.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -34,7 +38,7 @@ hid_t h5banktype() {
|
|||
|
||||
|
||||
void
|
||||
write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
|
||||
write_source_bank(hid_t group_id, Bank* source_bank)
|
||||
{
|
||||
hid_t banktype = h5banktype();
|
||||
|
||||
|
|
@ -46,11 +50,11 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
|
|||
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
|
||||
|
||||
// Create another data space but for each proc individually
|
||||
hsize_t count[] {static_cast<hsize_t>(openmc_work)};
|
||||
hsize_t count[] {static_cast<hsize_t>(simulation::work)};
|
||||
hid_t memspace = H5Screate_simple(1, count, nullptr);
|
||||
|
||||
// Select hyperslab for this dataspace
|
||||
hsize_t start[] {static_cast<hsize_t>(work_index[openmc::mpi::rank])};
|
||||
hsize_t start[] {static_cast<hsize_t>(simulation::work_index[mpi::rank])};
|
||||
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
|
||||
|
||||
// Set up the property list for parallel writing
|
||||
|
|
@ -77,24 +81,25 @@ write_source_bank(hid_t group_id, int64_t* work_index, 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 + openmc_work};
|
||||
std::vector<Bank> temp_source {source_bank, source_bank + simulation::work};
|
||||
#endif
|
||||
|
||||
for (int i = 0; i < openmc::mpi::n_procs; ++i) {
|
||||
for (int i = 0; i < mpi::n_procs; ++i) {
|
||||
// Create memory space
|
||||
hsize_t count[] {static_cast<hsize_t>(work_index[i+1] - work_index[i])};
|
||||
hsize_t count[] {static_cast<hsize_t>(simulation::work_index[i+1] -
|
||||
simulation::work_index[i])};
|
||||
hid_t memspace = H5Screate_simple(1, count, nullptr);
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
// Receive source sites from other processes
|
||||
if (i > 0)
|
||||
MPI_Recv(source_bank, count[0], openmc::mpi::bank, i, i,
|
||||
openmc::mpi::intracomm, MPI_STATUS_IGNORE);
|
||||
MPI_Recv(source_bank, count[0], mpi::bank, i, i,
|
||||
mpi::intracomm, MPI_STATUS_IGNORE);
|
||||
#endif
|
||||
|
||||
// Select hyperslab for this dataspace
|
||||
dspace = H5Dget_space(dset);
|
||||
hsize_t start[] {static_cast<hsize_t>(work_index[i])};
|
||||
hsize_t start[] {static_cast<hsize_t>(simulation::work_index[i])};
|
||||
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
|
||||
|
||||
// Write data to hyperslab
|
||||
|
|
@ -113,8 +118,8 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
|
|||
#endif
|
||||
} else {
|
||||
#ifdef OPENMC_MPI
|
||||
MPI_Send(source_bank, openmc_work, openmc::mpi::bank, 0, openmc::mpi::rank,
|
||||
openmc::mpi::intracomm);
|
||||
MPI_Send(source_bank, simulation::work, mpi::bank, 0, mpi::rank,
|
||||
mpi::intracomm);
|
||||
#endif
|
||||
}
|
||||
#endif
|
||||
|
|
@ -123,7 +128,7 @@ write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
|
|||
}
|
||||
|
||||
|
||||
void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
|
||||
void read_source_bank(hid_t group_id, Bank* source_bank)
|
||||
{
|
||||
hid_t banktype = h5banktype();
|
||||
|
||||
|
|
@ -131,20 +136,20 @@ void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
|
|||
hid_t dset = H5Dopen(group_id, "source_bank", H5P_DEFAULT);
|
||||
|
||||
// Create another data space but for each proc individually
|
||||
hsize_t dims[] {static_cast<hsize_t>(openmc_work)};
|
||||
hsize_t dims[] {static_cast<hsize_t>(simulation::work)};
|
||||
hid_t memspace = H5Screate_simple(1, dims, nullptr);
|
||||
|
||||
// Make sure source bank is big enough
|
||||
hid_t dspace = H5Dget_space(dset);
|
||||
hsize_t dims_all[1];
|
||||
H5Sget_simple_extent_dims(dspace, dims_all, nullptr);
|
||||
if (work_index[openmc::mpi::n_procs] > dims_all[0]) {
|
||||
if (simulation::work_index[mpi::n_procs] > dims_all[0]) {
|
||||
fatal_error("Number of source sites in source file is less "
|
||||
"than number of source particles per generation.");
|
||||
}
|
||||
|
||||
// Select hyperslab for each process
|
||||
hsize_t start[] {static_cast<hsize_t>(work_index[openmc::mpi::rank])};
|
||||
hsize_t start[] {static_cast<hsize_t>(simulation::work_index[mpi::rank])};
|
||||
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, dims, nullptr);
|
||||
|
||||
#ifdef PHDF5
|
||||
|
|
@ -164,4 +169,110 @@ void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
|
|||
H5Tclose(banktype);
|
||||
}
|
||||
|
||||
void write_tally_results_nr(hid_t file_id)
|
||||
{
|
||||
// ==========================================================================
|
||||
// COLLECT AND WRITE GLOBAL TALLIES
|
||||
|
||||
hid_t tallies_group;
|
||||
if (mpi::master) {
|
||||
// Write number of realizations
|
||||
write_dataset(file_id, "n_realizations", n_realizations);
|
||||
|
||||
// Write number of global tallies
|
||||
write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES);
|
||||
|
||||
tallies_group = open_group(file_id, "tallies");
|
||||
}
|
||||
|
||||
// Get pointer to global tallies
|
||||
auto gt = global_tallies();
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
// Reduce global tallies
|
||||
xt::xtensor<double, 2> gt_reduced = xt::empty_like(gt);
|
||||
MPI_Reduce(gt.data(), gt_reduced.data(), gt.size(), MPI_DOUBLE,
|
||||
MPI_SUM, 0, mpi::intracomm);
|
||||
|
||||
// Transfer values to value on master
|
||||
if (mpi::master) {
|
||||
if (simulation::current_batch == settings::n_max_batches ||
|
||||
simulation::satisfy_triggers) {
|
||||
std::copy(gt_reduced.begin(), gt_reduced.end(), gt.begin());
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
// Write out global tallies sum and sum_sq
|
||||
if (mpi::master) {
|
||||
write_dataset(file_id, "global_tallies", gt);
|
||||
}
|
||||
|
||||
for (int i = 1; i <= n_tallies; ++i) {
|
||||
// Skip any tallies that are not active
|
||||
bool active;
|
||||
openmc_tally_get_active(i, &active);
|
||||
if (!active) continue;
|
||||
|
||||
if (mpi::master && !object_exists(file_id, "tallies_present")) {
|
||||
write_attribute(file_id, "tallies_present", 1);
|
||||
}
|
||||
|
||||
// Get view of accumulated tally values
|
||||
auto results = tally_results(i);
|
||||
auto values_view = xt::view(results, xt::all(), xt::all(),
|
||||
xt::range(RESULT_SUM, RESULT_SUM_SQ + 1));
|
||||
|
||||
// Make copy of tally values in contiguous array
|
||||
xt::xtensor<double, 2> values = values_view;
|
||||
|
||||
if (mpi::master) {
|
||||
// Open group for tally
|
||||
int id;
|
||||
openmc_tally_get_id(i, &id);
|
||||
std::string groupname {"tally " + std::to_string(id)};
|
||||
hid_t tally_group = open_group(tallies_group, groupname.c_str());
|
||||
|
||||
// The MPI_IN_PLACE specifier allows the master to copy values into
|
||||
// a receive buffer without having a temporary variable
|
||||
#ifdef OPENMC_MPI
|
||||
MPI_Reduce(MPI_IN_PLACE, values.data(), values.size(), MPI_DOUBLE,
|
||||
MPI_SUM, 0, mpi::intracomm);
|
||||
#endif
|
||||
|
||||
// At the end of the simulation, store the results back in the
|
||||
// regular TallyResults array
|
||||
if (simulation::current_batch == settings::n_max_batches ||
|
||||
simulation::satisfy_triggers) {
|
||||
values_view = values;
|
||||
}
|
||||
|
||||
// Put in temporary tally result
|
||||
xt::xtensor<double, 3> results_copy = xt::zeros_like(results);
|
||||
auto copy_view = xt::view(results_copy, xt::all(), xt::all(),
|
||||
xt::range(RESULT_SUM, RESULT_SUM_SQ + 1));
|
||||
copy_view = values;
|
||||
|
||||
// Write reduced tally results to file
|
||||
auto shape = results_copy.shape();
|
||||
write_tally_results(tally_group, shape[0], shape[1], results_copy.data());
|
||||
|
||||
close_group(tally_group);
|
||||
} else {
|
||||
// Receive buffer not significant at other processors
|
||||
#ifdef OPENMC_MPI
|
||||
MPI_Reduce(values.data(), nullptr, values.size(), MPI_REAL8, MPI_SUM,
|
||||
0, mpi::intracomm);
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
if (mpi::master) {
|
||||
if (!object_exists(file_id, "tallies_present")) {
|
||||
// Indicate that tallies are off
|
||||
write_dataset(file_id, "tallies_present", 0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -3778,6 +3778,11 @@ contains
|
|||
real(C_DOUBLE) :: val
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
interface
|
||||
subroutine reduce_tally_results() bind(C)
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
! Combine tally results onto master process
|
||||
if (reduce_tallies) call reduce_tally_results()
|
||||
#endif
|
||||
|
|
@ -3821,67 +3826,6 @@ contains
|
|||
|
||||
end subroutine accumulate_tallies
|
||||
|
||||
!===============================================================================
|
||||
! REDUCE_TALLY_RESULTS collects all the results from tallies onto one processor
|
||||
!===============================================================================
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
subroutine reduce_tally_results()
|
||||
|
||||
integer :: i
|
||||
integer :: n ! number of filter bins
|
||||
integer :: m ! number of score bins
|
||||
integer :: n_bins ! total number of bins
|
||||
integer :: mpi_err ! MPI error code
|
||||
real(C_DOUBLE), allocatable :: tally_temp(:,:) ! contiguous array of results
|
||||
real(C_DOUBLE), allocatable :: tally_temp2(:,:) ! reduced contiguous results
|
||||
real(C_DOUBLE) :: temp(N_GLOBAL_TALLIES), temp2(N_GLOBAL_TALLIES)
|
||||
|
||||
do i = 1, active_tallies % size()
|
||||
associate (t => tallies(active_tallies % data(i)) % obj)
|
||||
|
||||
m = size(t % results, 2)
|
||||
n = size(t % results, 3)
|
||||
n_bins = m*n
|
||||
|
||||
allocate(tally_temp(m,n), tally_temp2(m,n))
|
||||
|
||||
! Reduce contiguous set of tally results
|
||||
tally_temp = t % results(RESULT_VALUE,:,:)
|
||||
call MPI_REDUCE(tally_temp, tally_temp2, n_bins, MPI_DOUBLE, &
|
||||
MPI_SUM, 0, mpi_intracomm, mpi_err)
|
||||
|
||||
if (master) then
|
||||
! Transfer values to value on master
|
||||
t % results(RESULT_VALUE,:,:) = tally_temp2
|
||||
else
|
||||
! Reset value on other processors
|
||||
t % results(RESULT_VALUE,:,:) = ZERO
|
||||
end if
|
||||
|
||||
deallocate(tally_temp, tally_temp2)
|
||||
end associate
|
||||
end do
|
||||
|
||||
! Reduce global tallies onto master
|
||||
temp = global_tallies(RESULT_VALUE, :)
|
||||
call MPI_REDUCE(temp, temp2, N_GLOBAL_TALLIES, MPI_DOUBLE, MPI_SUM, &
|
||||
0, mpi_intracomm, mpi_err)
|
||||
if (master) then
|
||||
global_tallies(RESULT_VALUE, :) = temp2
|
||||
else
|
||||
global_tallies(RESULT_VALUE, :) = ZERO
|
||||
end if
|
||||
|
||||
! We also need to determine the total starting weight of particles from the
|
||||
! last realization
|
||||
temp(1) = total_weight
|
||||
call MPI_REDUCE(temp, total_weight, 1, MPI_REAL8, MPI_SUM, &
|
||||
0, mpi_intracomm, mpi_err)
|
||||
|
||||
end subroutine reduce_tally_results
|
||||
#endif
|
||||
|
||||
!===============================================================================
|
||||
! SETUP_ACTIVE_TALLIES
|
||||
!===============================================================================
|
||||
|
|
|
|||
102
src/tallies/tally.cpp
Normal file
102
src/tallies/tally.cpp
Normal file
|
|
@ -0,0 +1,102 @@
|
|||
#include "openmc/tallies/tally.h"
|
||||
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/message_passing.h"
|
||||
|
||||
#include "xtensor/xadapt.hpp"
|
||||
#include "xtensor/xbuilder.hpp" // for empty_like
|
||||
#include "xtensor/xview.hpp"
|
||||
|
||||
#include <array>
|
||||
#include <cstddef>
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
adaptor_type<2> global_tallies()
|
||||
{
|
||||
// Get pointer to global tallies
|
||||
double* buffer;
|
||||
openmc_global_tallies(&buffer);
|
||||
|
||||
// Adapt into xtensor
|
||||
std::array<size_t, 2> shape = {N_GLOBAL_TALLIES, 3};
|
||||
std::size_t size {3*N_GLOBAL_TALLIES};
|
||||
|
||||
return xt::adapt(buffer, size, xt::no_ownership(), shape);
|
||||
}
|
||||
|
||||
adaptor_type<3> tally_results(int idx)
|
||||
{
|
||||
// Get pointer to tally results
|
||||
double* results;
|
||||
std::array<std::size_t, 3> shape;
|
||||
openmc_tally_results(idx, &results, shape.data());
|
||||
|
||||
// Adapt array into xtensor with no ownership
|
||||
std::size_t size {shape[0] * shape[1] * shape[2]};
|
||||
return xt::adapt(results, size, xt::no_ownership(), shape);
|
||||
}
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
void reduce_tally_results()
|
||||
{
|
||||
for (int i = 1; i <= n_tallies; ++i) {
|
||||
// Skip any tallies that are not active
|
||||
bool active;
|
||||
openmc_tally_get_active(i, &active);
|
||||
if (!active) continue;
|
||||
|
||||
// Get view of accumulated tally values
|
||||
auto results = tally_results(i);
|
||||
auto values_view = xt::view(results, xt::all(), xt::all(), RESULT_VALUE);
|
||||
|
||||
// Make copy of tally values in contiguous array
|
||||
xt::xtensor<double, 2> values = values_view;
|
||||
xt::xtensor<double, 2> values_reduced = xt::empty_like(values);
|
||||
|
||||
// Reduce contiguous set of tally results
|
||||
MPI_Reduce(values.data(), values_reduced.data(), values.size(),
|
||||
MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm);
|
||||
|
||||
// Transfer values on master and reset on other ranks
|
||||
if (mpi::master) {
|
||||
values_view = values_reduced;
|
||||
} else {
|
||||
values_view = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
// Get view of global tally values
|
||||
auto gt = global_tallies();
|
||||
auto gt_values_view = xt::view(gt, xt::all(), RESULT_VALUE);
|
||||
|
||||
// Make copy of values in contiguous array
|
||||
xt::xtensor<double, 1> gt_values = gt_values_view;
|
||||
xt::xtensor<double, 1> gt_values_reduced = xt::empty_like(gt_values);
|
||||
|
||||
// Reduce contiguous data
|
||||
MPI_Reduce(gt_values.data(), gt_values_reduced.data(), N_GLOBAL_TALLIES,
|
||||
MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm);
|
||||
|
||||
// Transfer values on master and reset on other ranks
|
||||
if (mpi::master) {
|
||||
gt_values_view = gt_values_reduced;
|
||||
} else {
|
||||
gt_values_view = 0.0;
|
||||
}
|
||||
|
||||
// We also need to determine the total starting weight of particles from the
|
||||
// last realization
|
||||
double weight_reduced;
|
||||
MPI_Reduce(&total_weight, &weight_reduced, 1, MPI_DOUBLE, MPI_SUM,
|
||||
0, mpi::intracomm);
|
||||
if (mpi::master) total_weight = weight_reduced;
|
||||
}
|
||||
#endif
|
||||
|
||||
} // namespace openmc
|
||||
|
|
@ -153,7 +153,7 @@ module tally_header
|
|||
|
||||
! Normalization for statistics
|
||||
integer(C_INT32_T), public, bind(C) :: n_realizations = 0 ! # of independent realizations
|
||||
real(8), public :: total_weight ! total starting particle weight in realization
|
||||
real(C_DOUBLE), public, bind(C) :: total_weight ! total starting particle weight in realization
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -681,14 +681,19 @@ contains
|
|||
! allows a user to obtain in-memory tally results from Python directly.
|
||||
integer(C_INT32_T), intent(in), value :: index
|
||||
type(C_PTR), intent(out) :: ptr
|
||||
integer(C_INT), intent(out) :: shape_(3)
|
||||
integer(C_SIZE_T), intent(out) :: shape_(3)
|
||||
integer(C_INT) :: err
|
||||
|
||||
if (index >= 1 .and. index <= size(tallies)) then
|
||||
associate (t => tallies(index) % obj)
|
||||
if (allocated(t % results)) then
|
||||
ptr = C_LOC(t % results(1,1,1))
|
||||
shape_(:) = shape(t % results)
|
||||
|
||||
! Note that shape is reversed since it is assumed to be used from
|
||||
! C/C++ code
|
||||
shape_(1) = size(t % results, 3)
|
||||
shape_(2) = size(t % results, 2)
|
||||
shape_(3) = size(t % results, 1)
|
||||
err = 0
|
||||
else
|
||||
err = E_ALLOCATE
|
||||
|
|
|
|||
60
src/timer.cpp
Normal file
60
src/timer.cpp
Normal file
|
|
@ -0,0 +1,60 @@
|
|||
#include "openmc/timer.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// Timer implementation
|
||||
//==============================================================================
|
||||
|
||||
void Timer::start ()
|
||||
{
|
||||
running_ = true;
|
||||
start_ = clock::now();
|
||||
}
|
||||
|
||||
void Timer::stop()
|
||||
{
|
||||
elapsed_ = elapsed();
|
||||
running_ = false;
|
||||
}
|
||||
|
||||
void Timer::reset()
|
||||
{
|
||||
running_ = false;
|
||||
elapsed_ = 0.0;
|
||||
}
|
||||
|
||||
double Timer::elapsed()
|
||||
{
|
||||
if (running_) {
|
||||
std::chrono::duration<double> diff = clock::now() - start_;
|
||||
return elapsed_ += diff.count();
|
||||
} else {
|
||||
return elapsed_;
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
Timer time_bank;
|
||||
Timer time_bank_sample;
|
||||
Timer time_bank_sendrecv;
|
||||
|
||||
//==============================================================================
|
||||
// Fortran compatibility
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double time_bank_elapsed() { return time_bank.elapsed(); }
|
||||
extern "C" double time_bank_sample_elapsed() { return time_bank_sample.elapsed(); }
|
||||
extern "C" double time_bank_sendrecv_elapsed() { return time_bank_sendrecv.elapsed(); }
|
||||
|
||||
extern "C" void reset_timers()
|
||||
{
|
||||
time_bank.reset();
|
||||
time_bank_sample.reset();
|
||||
time_bank_sendrecv.reset();
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
@ -1,9 +1,28 @@
|
|||
module timer_header
|
||||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use constants, only: ZERO
|
||||
|
||||
implicit none
|
||||
|
||||
interface
|
||||
function time_bank_elapsed() result(t) bind(C)
|
||||
import C_DOUBLE
|
||||
real(C_DOUBLE) :: t
|
||||
end function
|
||||
function time_bank_sample_elapsed() result(t) bind(C)
|
||||
import C_DOUBLE
|
||||
real(C_DOUBLE) :: t
|
||||
end function
|
||||
function time_bank_sendrecv_elapsed() result(t) bind(C)
|
||||
import C_DOUBLE
|
||||
real(C_DOUBLE) :: t
|
||||
end function
|
||||
subroutine reset_timers() bind(C)
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
!===============================================================================
|
||||
! TIMER represents a timer that can be started and stopped to measure how long
|
||||
! different routines run. The intrinsic routine system_clock is used to measure
|
||||
|
|
@ -29,9 +48,6 @@ module timer_header
|
|||
type(Timer) :: time_initialize ! timer for initialization
|
||||
type(Timer) :: time_read_xs ! timer for reading cross sections
|
||||
type(Timer) :: time_unionize ! timer for material xs-energy grid union
|
||||
type(Timer) :: time_bank ! timer for fission bank synchronization
|
||||
type(Timer) :: time_bank_sample ! timer for fission bank sampling
|
||||
type(Timer) :: time_bank_sendrecv ! timer for fission bank SEND/RECV
|
||||
type(Timer) :: time_tallies ! timer for accumulate tallies
|
||||
type(Timer) :: time_inactive ! timer for inactive batches
|
||||
type(Timer) :: time_active ! timer for active batches
|
||||
|
|
|
|||
|
|
@ -322,7 +322,9 @@ contains
|
|||
real(8) :: norm ! "norm" of surface normal
|
||||
real(8) :: xyz(3) ! Saved global coordinate
|
||||
integer :: i_surface ! index in surfaces
|
||||
#ifdef DAGMC
|
||||
integer :: i_cell ! index of new cell
|
||||
#endif
|
||||
logical :: rotational ! if rotational periodic BC applied
|
||||
logical :: found ! particle found in universe?
|
||||
class(Surface), pointer :: surf
|
||||
|
|
|
|||
|
|
@ -138,16 +138,35 @@ contains
|
|||
integer :: min_samples ! minimum number of samples per process
|
||||
integer :: remainder ! leftover samples from uneven divide
|
||||
#ifdef OPENMC_MPI
|
||||
integer :: mpi_err ! MPI error code
|
||||
integer :: m ! index over materials
|
||||
integer :: n ! number of materials
|
||||
integer, allocatable :: data(:) ! array used to send number of hits
|
||||
integer(C_INT) :: n ! number of materials
|
||||
integer(C_INT), allocatable :: data(:) ! array used to send number of hits
|
||||
#endif
|
||||
real(8) :: f ! fraction of hits
|
||||
real(8) :: var_f ! variance of fraction of hits
|
||||
real(8) :: volume_sample ! total volume of sampled region
|
||||
real(8) :: atoms(2, size(nuclides))
|
||||
|
||||
#ifdef OPENMC_MPI
|
||||
interface
|
||||
subroutine send_int(buffer, count, dest, tag) bind(C)
|
||||
import C_INT
|
||||
integer(C_INT), intent(in) :: buffer
|
||||
integer(C_INT), value :: count
|
||||
integer(C_INT), value :: dest
|
||||
integer(C_INT), value :: tag
|
||||
end subroutine
|
||||
|
||||
subroutine recv_int(buffer, count, source, tag) bind(C)
|
||||
import C_INT
|
||||
integer(C_INT), intent(out) :: buffer
|
||||
integer(C_INT), value :: count
|
||||
integer(C_INT), value :: source
|
||||
integer(C_INT), value :: tag
|
||||
end subroutine
|
||||
end interface
|
||||
#endif
|
||||
|
||||
! Divide work over MPI processes
|
||||
min_samples = this % samples / n_procs
|
||||
remainder = mod(this % samples, n_procs)
|
||||
|
|
@ -285,12 +304,10 @@ contains
|
|||
if (master) then
|
||||
#ifdef OPENMC_MPI
|
||||
do j = 1, n_procs - 1
|
||||
call MPI_RECV(n, 1, MPI_INTEGER, j, 0, mpi_intracomm, &
|
||||
MPI_STATUS_IGNORE, mpi_err)
|
||||
call recv_int(n, 1, j, 0)
|
||||
|
||||
allocate(data(2*n))
|
||||
call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, mpi_intracomm, &
|
||||
MPI_STATUS_IGNORE, mpi_err)
|
||||
call recv_int(data(1), 2*n, j, 1)
|
||||
do k = 0, n - 1
|
||||
do m = 1, master_indices(i_domain) % size()
|
||||
if (data(2*k + 1) == master_indices(i_domain) % data(m)) then
|
||||
|
|
@ -353,8 +370,8 @@ contains
|
|||
data(2*k + 2) = master_hits(i_domain) % data(k + 1)
|
||||
end do
|
||||
|
||||
call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, mpi_intracomm, mpi_err)
|
||||
call MPI_SEND(data, 2*n, MPI_INTEGER, 0, 1, mpi_intracomm, mpi_err)
|
||||
call send_int(n, 1, 0, 0)
|
||||
call send_int(data(1), 2*n, 0, 1)
|
||||
deallocate(data)
|
||||
#endif
|
||||
end if
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue