crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge remote-tracking branch 'upstream/develop' into cpp_tallies

Browse source
This commit is contained in:
Sterling Harper 2018-10-17 19:56:59 -04:00
commit 9b5e402a11
72 changed files with 1985 additions and 1923 deletions
Download patch file Download diff file Expand all files Collapse all files

13
CMakeLists.txt
View file

@ -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
@ -345,7 +344,6 @@ add_library(libopenmc SHARED
src/photon_physics.F90
src/physics_common.F90
src/physics.F90
src/physics_mg.F90
src/plot.F90
src/plot_header.F90
src/progress_header.F90
@ -408,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
@ -423,6 +422,8 @@ add_library(libopenmc SHARED
src/nuclide.cpp
src/output.cpp
src/particle.cpp
src/physics_common.cpp
src/physics_mg.cpp
src/plot.cpp
src/position.cpp
src/pugixml/pugixml_c.cpp
@ -441,10 +442,13 @@ add_library(libopenmc SHARED
src/string_functions.cpp
src/summary.cpp
src/surface.cpp
src/tallies/tally.cpp
src/tallies/tally_filter.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 +475,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

862
examples/jupyter/pandas-dataframes.ipynb
View file

File diff suppressed because one or more lines are too long

10
include/openmc/capi.h
View file

@ -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,15 +147,12 @@ extern "C" {
extern int32_t n_surfaces;
extern int32_t n_tallies;
extern int32_t n_universes;
extern bool openmc_simulation_initialized;
// 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
}

15
include/openmc/constants.h
View file

@ -315,9 +315,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
@ -408,10 +408,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};

7
include/openmc/eigenvalue.h
View file

@ -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();

5
include/openmc/error.h
View file

@ -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

4
include/openmc/finalize.h
View file

@ -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

11
include/openmc/hdf5_interface.h
View file

@ -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,

1
include/openmc/material.h
View file

@ -26,6 +26,7 @@ class Material
public:
int32_t id_; //!< Unique ID
double volume_ {-1.0}; //!< Volume in [cm^3]
bool fissionable {false}; //!< Does this material contain fissionable nuclides
//! \brief Default temperature for cells containing this material.
//!

8
include/openmc/math_functions.h
View file

@ -22,7 +22,7 @@ namespace openmc {
//! @return The requested percentile
//==============================================================================
extern "C" double normal_percentile_c(double p);
extern "C" double normal_percentile(double p);
//==============================================================================
//! Calculate the percentile of the Student's t distribution with a specified
@ -58,7 +58,7 @@ extern "C" void calc_pn_c(int n, double x, double pnx[]);
//! evaluated at x
//==============================================================================
extern "C" double evaluate_legendre_c(int n, const double data[], double x);
extern "C" double evaluate_legendre(int n, const double data[], double x);
//==============================================================================
//! Calculate the n-th order real spherical harmonics for a given angle (in
@ -144,7 +144,7 @@ Direction rotate_angle(Direction u, double mu, double* phi);
//! @result The sampled outgoing energy
//==============================================================================
extern "C" double maxwell_spectrum_c(double T);
extern "C" double maxwell_spectrum(double T);
//==============================================================================
//! Samples an energy from a Watt energy-dependent fission distribution.
@ -159,7 +159,7 @@ extern "C" double maxwell_spectrum_c(double T);
//! @result The sampled outgoing energy
//==============================================================================
extern "C" double watt_spectrum_c(double a, double b);
extern "C" double watt_spectrum(double a, double b);
//==============================================================================
//! Doppler broadens the windowed multipole curvefit.

2
include/openmc/message_passing.h
View file

@ -2,7 +2,7 @@
#define OPENMC_MESSAGE_PASSING_H
#ifdef OPENMC_MPI
#include "mpi.h"
#include <mpi.h>
#endif
namespace openmc {

10
include/openmc/mgxs_interface.h
View file

@ -18,6 +18,9 @@ extern std::vector<Mgxs> nuclides_MG;
extern std::vector<Mgxs> macro_xs;
extern "C" int num_energy_groups;
//TODO: When more of the Fortran is converted (input_xml, tallies, etc, also
// bring over energy_bin_avg, energy_bins, etc, as vectors)
//==============================================================================
// Mgxs data loading interface methods
//==============================================================================
@ -44,13 +47,6 @@ extern "C" void
calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs);
extern "C" void
sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
double uvw[3]);
extern "C" void
sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout);
extern "C" double
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);

16
include/openmc/physics_common.h Normal file
View file

@ -0,0 +1,16 @@
//! \file physics_common.h
//! A collection of physics methods common to MG, CE, photon, etc.
#ifndef OPENMC_PHYSICS_COMMON_H
#define OPENMC_PHYSICS_COMMON_H
#include "openmc/particle.h"
namespace openmc {
//! \brief Performs the russian roulette operation for a particle
extern "C" void
russian_roulette(Particle* p);
} // namespace openmc
#endif // OPENMC_PHYSICS_COMMON_H

60
include/openmc/physics_mg.h Normal file
View file

@ -0,0 +1,60 @@
//! \file physics_mg.h
//! Methods needed to perform the collision physics for multi-group mode
#ifndef OPENMC_PHYSICS_MG_H
#define OPENMC_PHYSICS_MG_H
#include "openmc/capi.h"
#include "openmc/particle.h"
#include "openmc/nuclide.h"
namespace openmc {
//TODO: Remove energy_bin_avg and material_xs parameters when they reside on
// the C-side this should happen after materials, physics, input, and tallies
// are brought over
//! \brief samples particle behavior after a collision event.
//! \param p Particle to operate on
//! \param energy_bin_avg Average energy within each energy bin
//! \param material_xs The cross section cache for the current material
extern "C" void
collision_mg(Particle* p, const double* energy_bin_avg,
const MaterialMacroXS* material_xs);
//! \brief samples a reaction type.
//!
//! Note that there is special logic when suvival biasing is turned on since
//! fission and disappearance are treated implicitly.
//! \param p Particle to operate on
//! \param energy_bin_avg Average energy within each energy bin
//! \param material_xs The cross section cache for the current material
void
sample_reaction(Particle* p, const double* energy_bin_avg,
const MaterialMacroXS* material_xs);
//! \brief Samples the scattering event
//! \param p Particle to operate on
//! \param energy_bin_avg Average energy within each energy bin
void
scatter(Particle* p, const double* energy_bin_avg);
//! \brief Determines the average total, prompt and delayed neutrons produced
//! from fission and creates the appropriate bank sites.
//! \param p Particle to operate on
//! \param bank_array The particle bank to populate
//! \param size_bank Number of particles currently in the bank
//! \param bank_array_size Allocated size of the bank
//! \param material_xs The cross section cache for the current material
void
create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
int64_t bank_array_size, const MaterialMacroXS* material_xs);
//! \brief Handles an absorption event
//! \param p Particle to operate on
//! \param material_xs The cross section cache for the current material
void
absorption(Particle* p, const MaterialMacroXS* material_xs);
} // namespace openmc
#endif // OPENMC_PHYSICS_MG_H

10
include/openmc/search.h
View file

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

53
include/openmc/simulation.h
View file

@ -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

8
include/openmc/state_point.h
View file

@ -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
View 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
View 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

6
openmc/capi/core.py
View file

@ -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)

23
openmc/capi/math.py
View file

@ -12,8 +12,8 @@ _dll.t_percentile_c.argtypes = [c_double, c_int]
_dll.calc_pn_c.restype = None
_dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)]
_dll.evaluate_legendre_c.restype = c_double
_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double]
_dll.evaluate_legendre.restype = c_double
_dll.evaluate_legendre.argtypes = [c_int, POINTER(c_double), c_double]
_dll.calc_rn_c.restype = None
_dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
@ -27,11 +27,11 @@ _dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
_dll.rotate_angle_c.restype = None
_dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
POINTER(c_double)]
_dll.maxwell_spectrum_c.restype = c_double
_dll.maxwell_spectrum_c.argtypes = [c_double]
_dll.maxwell_spectrum.restype = c_double
_dll.maxwell_spectrum.argtypes = [c_double]
_dll.watt_spectrum_c.restype = c_double
_dll.watt_spectrum_c.argtypes = [c_double, c_double]
_dll.watt_spectrum.restype = c_double
_dll.watt_spectrum.argtypes = [c_double, c_double]
_dll.broaden_wmp_polynomials_c.restype = None
_dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int,
@ -100,9 +100,8 @@ def evaluate_legendre(data, x):
"""
data_arr = np.array(data, dtype=np.float64)
return _dll.evaluate_legendre_c(len(data),
data_arr.ctypes.data_as(POINTER(c_double)),
x)
return _dll.evaluate_legendre(len(data),
data_arr.ctypes.data_as(POINTER(c_double)), x)
def calc_rn(n, uvw):
@ -182,7 +181,7 @@ def calc_zn_rad(n, rho):
zn_rad = np.zeros(num_bins, dtype=np.float64)
_dll.calc_zn_rad_c(n, rho, zn_rad)
return zn_rad
def rotate_angle(uvw0, mu, phi=None):
""" Rotates direction cosines through a polar angle whose cosine is
@ -231,7 +230,7 @@ def maxwell_spectrum(T):
"""
return _dll.maxwell_spectrum_c(T)
return _dll.maxwell_spectrum(T)
def watt_spectrum(a, b):
@ -251,7 +250,7 @@ def watt_spectrum(a, b):
"""
return _dll.watt_spectrum_c(a, b)
return _dll.watt_spectrum(a, b)
def broaden_wmp_polynomials(E, dopp, n):

8
openmc/capi/tally.py
View file

@ -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):

2
openmc/mgxs/mgxs.py
View file

@ -2774,7 +2774,7 @@ class TransportXS(MGXS):
p1_tally = p1_tally.get_slice(filters=[openmc.LegendreFilter],
filter_bins=[('P1',)],
squeeze=True)
p1_tally.scores = ['scatter-1']
p1_tally._scores = ['scatter-1']
# Compute total cross section
total_xs = self.tallies['total'] / self.tallies['flux (tracklength)']

17
openmc/mgxs_library.py
View file

@ -1840,10 +1840,10 @@ class XSdata(object):
mu = np.linspace(-1, 1, xsdata.num_orders)
# Evaluate the legendre on the mu grid
for imu in range(len(mu)):
new_data[..., imu] = \
np.sum((l + 0.5) * eval_legendre(l, mu[imu]) *
orig_data[..., l]
for l in range(self.num_orders))
for l in range(self.num_orders):
new_data[..., imu] += (
(l + 0.5) * eval_legendre(l, mu[imu]) *
orig_data[..., l])
elif target_format == 'histogram':
# This code uses the vectorized integration capabilities
@ -1857,11 +1857,10 @@ class XSdata(object):
mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU)
table_fine = np.zeros(new_data.shape[:-1] + (_NMU,))
for imu in range(len(mu_fine)):
table_fine[..., imu] = \
np.sum((l + 0.5) *
eval_legendre(l, mu_fine[imu]) *
orig_data[..., l]
for l in range(self.num_orders))
for l in range(self.num_orders):
table_fine[..., imu] += ((l + 0.5)
* eval_legendre(l, mu_fine[imu]) *
orig_data[..., l])
new_data[..., h_bin] = simps(table_fine, mu_fine)
elif self.scatter_format == 'tabular':

4
openmc/tallies.py
View file

@ -376,7 +376,7 @@ class Tally(IDManagerMixin):
def scores(self, scores):
cv.check_type('tally scores', scores, MutableSequence)
for i, score in enumerate(scores[:-1]):
for i, score in enumerate(scores):
# If the score is already in the Tally, raise an error
if score in scores[i+1:]:
msg = 'Unable to add a duplicate score "{0}" to Tally ID="{1}" ' \
@ -390,7 +390,7 @@ class Tally(IDManagerMixin):
for deprecated in ['scatter-', 'nu-scatter-', 'scatter-p',
'nu-scatter-p', 'scatter-y', 'nu-scatter-y',
'flux-y', 'total-y']:
if score.startswith(deprecated):
if score.strip().startswith(deprecated):
msg = score.strip() + ' is no longer supported.'
raise ValueError(msg)
scores[i] = score.strip()

5
src/api.F90
View file

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

21
src/cmfd_execute.F90
View file

@ -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

13
src/cmfd_execute.cpp
View file

@ -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

4
src/cmfd_header.F90
View file

@ -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(:,:)

4
src/distribution.cpp
View file

@ -90,7 +90,7 @@ Maxwell::Maxwell(pugi::xml_node node)
double Maxwell::sample() const
{
return maxwell_spectrum_c(theta_);
return maxwell_spectrum(theta_);
}
//==============================================================================
@ -110,7 +110,7 @@ Watt::Watt(pugi::xml_node node)
double Watt::sample() const
{
return watt_spectrum_c(a_, b_);
return watt_spectrum(a_, b_);
}
//==============================================================================

4
src/distribution_energy.cpp
View file

@ -279,7 +279,7 @@ double MaxwellEnergy::sample(double E) const
while (true) {
// Sample maxwell fission spectrum
double E_out = maxwell_spectrum_c(theta);
double E_out = maxwell_spectrum(theta);
// Accept energy based on restriction energy
if (E_out <= E - u_) return E_out;
@ -343,7 +343,7 @@ double WattEnergy::sample(double E) const
while (true) {
// Sample energy-dependent Watt fission spectrum
double E_out = watt_spectrum_c(a, b);
double E_out = watt_spectrum(a, b);
// Accept energy based on restriction energy
if (E_out <= E - u_) return E_out;

314
src/eigenvalue.F90
View file

@ -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

299
src/eigenvalue.cpp
View file

@ -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)

12
src/error.F90
View file

@ -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
View 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

6
src/finalize.cpp
View file

@ -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

4
src/geometry.cpp
View file

@ -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 << ","

3
src/hdf5_interface.F90
View file

@ -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

80
src/initialize.F90
View file

@ -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
!===============================================================================

37
src/initialize.cpp
View file

@ -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.");

4
src/input_xml.F90
View file

@ -219,7 +219,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
@ -2975,7 +2975,7 @@ contains
! Check if material is fissionable
if (nuclides(materials(i) % nuclide(j)) % fissionable) then
materials(i) % fissionable = .true.
call materials(i) % set_fissionable(.true.)
end if
end do

7
src/material.cpp
View file

@ -122,6 +122,13 @@ extern "C" {
material_map[id] = index - 1;
}
bool material_fissionable(Material* mat) {return mat->fissionable;}
void material_set_fissionable(Material* mat, bool fissionable)
{
mat->fissionable = fissionable;
}
void extend_materials_c(int32_t n)
{
materials.reserve(materials.size() + n);

31
src/material_header.F90
View file

@ -51,6 +51,20 @@ module material_header
integer(C_INT32_T), intent(in), value :: index
end subroutine material_set_id_c
function material_fissionable_c(mat_ptr) &
bind(C, name='material_fissionable') result(fissionable)
import C_PTR, C_BOOL
type(C_PTR), intent(in), value :: mat_ptr
logical(C_BOOL) :: fissionable
end function material_fissionable_c
subroutine material_set_fissionable_c(mat_ptr, fissionable) &
bind(C, name='material_set_fissionable')
import C_PTR, C_BOOL
type(C_PTR), intent(in), value :: mat_ptr
logical(C_BOOL), intent(in), value :: fissionable
end subroutine material_set_fissionable_c
subroutine extend_materials_c(n) bind(C)
import C_INT32_T
integer(C_INT32_T), intent(in), value :: n
@ -95,7 +109,6 @@ module material_header
character(20), allocatable :: sab_names(:) ! name of S(a,b) table
! Does this material contain fissionable nuclides? Is it depletable?
logical :: fissionable = .false.
logical :: depletable = .false.
! enforce isotropic scattering in lab for specific nuclides
@ -105,6 +118,8 @@ module material_header
contains
procedure :: id => material_id
procedure :: set_id => material_set_id
procedure :: fissionable => material_fissionable
procedure :: set_fissionable => material_set_fissionable
procedure :: set_density => material_set_density
procedure :: init_nuclide_index => material_init_nuclide_index
procedure :: assign_sab_tables => material_assign_sab_tables
@ -139,6 +154,18 @@ contains
call material_set_id_c(this % ptr, id, index)
end subroutine material_set_id
function material_fissionable(this) result(fissionable)
class(Material), intent(in) :: this
logical(C_BOOL) :: fissionable
fissionable = material_fissionable_c(this % ptr)
end function material_fissionable
subroutine material_set_fissionable(this, fissionable)
class(Material),intent(in) :: this
logical, intent(in) :: fissionable
call material_set_fissionable_c(this % ptr, logical(fissionable, C_BOOL))
end subroutine material_set_fissionable
function material_set_density(this, density) result(err)
class(Material), intent(inout) :: this
real(8), intent(in) :: density
@ -685,7 +712,7 @@ contains
integer(C_INT) :: err
if (index >= 1 .and. index <= size(materials)) then
fissionable = materials(index) % fissionable
fissionable = materials(index) % fissionable()
err = 0
else
err = E_OUT_OF_BOUNDS

44
src/math.F90
View file

@ -12,10 +12,7 @@ module math
public :: calc_rn
public :: calc_zn
public :: calc_zn_rad
public :: evaluate_legendre
public :: rotate_angle
public :: maxwell_spectrum
public :: watt_spectrum
public :: faddeeva
public :: w_derivative
public :: broaden_wmp_polynomials
@ -42,16 +39,6 @@ module math
real(C_DOUBLE), intent(out) :: pnx(n + 1)
end subroutine calc_pn
pure function evaluate_legendre_c_intfc(n, data, x) &
bind(C, name='evaluate_legendre_c') result(val)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n
real(C_DOUBLE), intent(in) :: data(n)
real(C_DOUBLE), value, intent(in) :: x
real(C_DOUBLE) :: val
end function evaluate_legendre_c_intfc
pure subroutine calc_rn(n, uvw, rn) bind(C, name='calc_rn_c')
use ISO_C_BINDING
implicit none
@ -85,23 +72,6 @@ module math
real(C_DOUBLE), optional, intent(in) :: phi
end subroutine rotate_angle_c_intfc
function maxwell_spectrum(T) bind(C, name='maxwell_spectrum_c') &
result(E_out)
use ISO_C_BINDING
implicit none
real(C_DOUBLE), value, intent(in) :: T
real(C_DOUBLE) :: E_out
end function maxwell_spectrum
function watt_spectrum(a, b) bind(C, name='watt_spectrum_c') &
result(E_out)
use ISO_C_BINDING
implicit none
real(C_DOUBLE), value, intent(in) :: a
real(C_DOUBLE), value, intent(in) :: b
real(C_DOUBLE) :: E_out
end function watt_spectrum
subroutine broaden_wmp_polynomials(E, dopp, n, factors) &
bind(C, name='broaden_wmp_polynomials_c')
use ISO_C_BINDING
@ -157,20 +127,6 @@ module math
contains
!===============================================================================
! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
! and the value of x
!===============================================================================
pure function evaluate_legendre(data, x) result(val) bind(C)
real(C_DOUBLE), intent(in) :: data(:)
real(C_DOUBLE), intent(in) :: x
real(C_DOUBLE) :: val
val = evaluate_legendre_c_intfc(size(data) - 1, data, x)
end function evaluate_legendre
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done

12
src/math_functions.cpp
View file

@ -6,7 +6,7 @@ namespace openmc {
// Mathematical methods
//==============================================================================
double normal_percentile_c(double p) {
double normal_percentile(double p) {
constexpr double p_low = 0.02425;
constexpr double a[6] = {-3.969683028665376e1, 2.209460984245205e2,
-2.759285104469687e2, 1.383577518672690e2,
@ -79,7 +79,7 @@ double t_percentile_c(double p, int df){
// 16 (4), pp. 1123-1132 (1987).
double n = df;
double k = 1. / (n - 2.);
double z = normal_percentile_c(p);
double z = normal_percentile(p);
double z2 = z * z;
t = std::sqrt(n * k) * (z + (z2 - 3.) * z * k / 4. + ((5. * z2 - 56.) * z2 +
75.) * z * k * k / 96. + (((z2 - 27.) * 3. * z2 + 417.) * z2 - 315.) *
@ -103,7 +103,7 @@ void calc_pn_c(int n, double x, double pnx[]) {
}
double evaluate_legendre_c(int n, const double data[], double x) {
double evaluate_legendre(int n, const double data[], double x) {
double pnx[n + 1];
double val = 0.0;
calc_pn_c(n, x, pnx);
@ -658,7 +658,7 @@ Direction rotate_angle(Direction u, double mu, double* phi)
}
double maxwell_spectrum_c(double T) {
double maxwell_spectrum(double T) {
// Set the random numbers
double r1 = prn();
double r2 = prn();
@ -674,8 +674,8 @@ double maxwell_spectrum_c(double T) {
}
double watt_spectrum_c(double a, double b) {
double w = maxwell_spectrum_c(a);
double watt_spectrum(double a, double b) {
double w = maxwell_spectrum(a);
double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w);
return E_out;

22
src/message_passing.F90
View file

@ -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

19
src/message_passing.cpp
View file

@ -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

3
src/mgxs_data.F90
View file

@ -90,7 +90,8 @@ contains
end if
end do NUCLIDE_LOOP
mat % fissionable = query_fissionable_c(mat % n_nuclides, mat % nuclide)
call mat % set_fissionable( &
logical(query_fissionable_c(mat % n_nuclides, mat % nuclide)))
end do MATERIAL_LOOP

22
src/mgxs_interface.F90
View file

@ -62,26 +62,6 @@ module mgxs_interface
real(C_DOUBLE), intent(inout) :: nu_fiss_xs
end subroutine calculate_xs_c
subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: i_mat
integer(C_INT), value, intent(in) :: gin
integer(C_INT), intent(inout) :: gout
real(C_DOUBLE), intent(inout) :: mu
real(C_DOUBLE), intent(inout) :: wgt
real(C_DOUBLE), intent(inout) :: uvw(1:3)
end subroutine sample_scatter_c
subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: i_mat
integer(C_INT), value, intent(in) :: gin
integer(C_INT), intent(inout) :: dg
integer(C_INT), intent(inout) :: gout
end subroutine sample_fission_energy_c
subroutine get_name_c(index, name_len, name) bind(C)
use ISO_C_BINDING
implicit none
@ -156,7 +136,7 @@ module mgxs_interface
real(8), allocatable :: energy_bins(:)
! Midpoint of the energy group structure
real(8), allocatable :: energy_bin_avg(:)
real(C_DOUBLE), allocatable :: energy_bin_avg(:)
! Energy group structure with increasing energy
real(C_DOUBLE), allocatable, target :: rev_energy_bins(:)

30
src/mgxs_interface.cpp
View file

@ -97,36 +97,6 @@ calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
//==============================================================================
void
sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
double uvw[3])
{
int gout_c = gout - 1;
macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
// adjust return value for fortran indexing
gout = gout_c + 1;
// Rotate the angle
rotate_angle_c(uvw, mu, nullptr);
}
//==============================================================================
void
sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout)
{
int dg_c = 0;
int gout_c = 0;
macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
// adjust return values for fortran indexing
dg = dg_c + 1;
gout = gout_c + 1;
}
//==============================================================================
double
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
{

10
src/output.F90
View file

@ -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

16
src/particle.cpp
View file

@ -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};

31
src/physics_common.F90
View file

@ -1,33 +1,20 @@
module physics_common
use constants
use, intrinsic :: ISO_C_BINDING
use particle_header, only: Particle
use random_lcg, only: prn
use settings, only: weight_cutoff, weight_survive
implicit none
contains
!===============================================================================
! RUSSIAN_ROULETTE
! RUSSIAN_ROULETTE FROM C
!===============================================================================
subroutine russian_roulette(p)
type(Particle), intent(inout) :: p
if (p % wgt < weight_cutoff) then
if (prn() < p % wgt / weight_survive) then
p % wgt = weight_survive
p % last_wgt = p % wgt
else
p % wgt = ZERO
p % last_wgt = ZERO
p % alive = .false.
end if
end if
end subroutine russian_roulette
interface
subroutine russian_roulette(p) bind(C)
import Particle
type(Particle), intent(inout) :: p
end subroutine russian_roulette
end interface
end module physics_common

26
src/physics_common.cpp Normal file
View file

@ -0,0 +1,26 @@
#include "openmc/physics_common.h"
#include "openmc/settings.h"
#include "openmc/random_lcg.h"
namespace openmc {
//==============================================================================
// RUSSIAN_ROULETTE
//==============================================================================
void russian_roulette(Particle* p)
{
if (p->wgt < settings::weight_cutoff) {
if (prn() < p->wgt / settings::weight_survive) {
p->wgt = settings::weight_survive;
p->last_wgt = p->wgt;
} else {
p->wgt = 0.;
p->last_wgt = 0.;
p->alive = false;
}
}
}
} //namespace openmc

275
src/physics_mg.F90
View file

@ -1,275 +0,0 @@
module physics_mg
! This module contains the multi-group specific physics routines so as to not
! hinder performance of the CE versions with multiple if-thens.
use bank_header
use constants
use error, only: fatal_error, warning, write_message
use material_header, only: Material, materials
use math, only: rotate_angle
use mgxs_interface
use message_passing
use nuclide_header, only: material_xs
use particle_header
use physics_common
use random_lcg, only: prn
use settings
use simulation_header
use string, only: to_str
use tally_header
implicit none
contains
!===============================================================================
! COLLISION_MG samples a nuclide and reaction and then calls the appropriate
! routine for that reaction
!===============================================================================
subroutine collision_mg(p)
type(Particle), intent(inout) :: p
! Add to collision counter for particle
p % n_collision = p % n_collision + 1
! Sample nuclide/reaction for the material the particle is in
call sample_reaction(p)
! Display information about collision
if (verbosity >= 10 .or. trace) then
call write_message(" " // "Energy Group = " // trim(to_str(p % g)))
end if
end subroutine collision_mg
!===============================================================================
! SAMPLE_REACTION samples a nuclide based on the macroscopic cross sections for
! each nuclide within a material and then samples a reaction for that nuclide
! and calls the appropriate routine to process the physics. Note that there is
! special logic when suvival biasing is turned on since fission and
! disappearance are treated implicitly.
!===============================================================================
subroutine sample_reaction(p)
type(Particle), intent(inout) :: p
type(Material), pointer :: mat
mat => materials(p % material)
! Create fission bank sites. Note that while a fission reaction is sampled,
! it never actually "happens", i.e. the weight of the particle does not
! change when sampling fission sites. The following block handles all
! absorption (including fission)
if (mat % fissionable) then
if (run_mode == MODE_EIGENVALUE) then
call create_fission_sites(p, fission_bank, n_bank)
elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then
call create_fission_sites(p, p % secondary_bank, p % n_secondary)
end if
end if
! If survival biasing is being used, the following subroutine adjusts the
! weight of the particle. Otherwise, it checks to see if absorption occurs
if (material_xs % absorption > ZERO) then
call absorption(p)
else
p % absorb_wgt = ZERO
end if
if (.not. p % alive) return
! Sample a scattering reaction and determine the secondary energy of the
! exiting neutron
call scatter(p)
! Play russian roulette if survival biasing is turned on
if (survival_biasing) then
call russian_roulette(p)
if (.not. p % alive) return
end if
end subroutine sample_reaction
!===============================================================================
! ABSORPTION
!===============================================================================
subroutine absorption(p)
type(Particle), intent(inout) :: p
if (survival_biasing) then
! Determine weight absorbed in survival biasing
p % absorb_wgt = (p % wgt * &
material_xs % absorption / material_xs % total)
! Adjust weight of particle by probability of absorption
p % wgt = p % wgt - p % absorb_wgt
p % last_wgt = p % wgt
! Score implicit absorption estimate of keff
!$omp atomic
global_tallies(RESULT_VALUE, K_ABSORPTION) = &
global_tallies(RESULT_VALUE, K_ABSORPTION) + p % absorb_wgt * &
material_xs % nu_fission / material_xs % absorption
else
! See if disappearance reaction happens
if (material_xs % absorption > prn() * material_xs % total) then
! Score absorption estimate of keff
!$omp atomic
global_tallies(RESULT_VALUE, K_ABSORPTION) = &
global_tallies(RESULT_VALUE, K_ABSORPTION) + p % wgt * &
material_xs % nu_fission / material_xs % absorption
p % alive = .false.
p % event = EVENT_ABSORB
end if
end if
end subroutine absorption
!===============================================================================
! SCATTER
!===============================================================================
subroutine scatter(p)
type(Particle), intent(inout) :: p
call sample_scatter_c(p % material, p % last_g, p % g, p % mu, &
p % wgt, p % coord(1) % uvw)
! Update energy value for downstream compatability (in tallying)
p % E = energy_bin_avg(p % g)
! Set event component
p % event = EVENT_SCATTER
end subroutine scatter
!===============================================================================
! CREATE_FISSION_SITES determines the average total, prompt, and delayed
! neutrons produced from fission and creates appropriate bank sites.
!===============================================================================
subroutine create_fission_sites(p, bank_array, size_bank)
type(Particle), intent(inout) :: p
type(Bank), intent(inout) :: bank_array(:)
integer(8), intent(inout) :: size_bank
integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born
integer :: i ! loop index
integer :: dg ! delayed group
integer :: gout ! group out
integer :: nu ! actual number of neutrons produced
real(8) :: nu_t ! total nu
real(8) :: mu ! fission neutron angular cosine
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
interface
function ufs_get_weight(p) result(weight) bind(C)
import Particle, C_DOUBLE
type(Particle), intent(in) :: p
real(C_DOUBLE) :: WEIGHT
end function
end interface
! TODO: Heat generation from fission
! If uniform fission source weighting is turned on, we increase of decrease
! the expected number of fission sites produced
if (ufs) then
weight = ufs_get_weight(p)
else
weight = ONE
end if
! Determine expected number of neutrons produced
nu_t = p % wgt / keff * weight * &
material_xs % nu_fission / material_xs % total
! Sample number of neutrons produced
if (prn() > nu_t - int(nu_t)) then
nu = int(nu_t)
else
nu = int(nu_t) + 1
end if
! Check for bank size getting hit. For fixed source calculations, this is a
! fatal error. For eigenvalue calculations, it just means that k-effective
! was too high for a single batch.
if (size_bank + nu > size(bank_array)) then
if (run_mode == MODE_FIXEDSOURCE) then
call fatal_error("Secondary particle bank size limit reached. If you &
&are running a subcritical multiplication problem, k-effective &
&may be too close to one.")
else
if (master) call warning("Maximum number of sites in fission bank &
&reached. This can result in irreproducible results using different &
&numbers of processes/threads.")
end if
end if
! Bank source neutrons
if (nu == 0 .or. size_bank == size(bank_array)) return
! Initialize counter of delayed neutrons encountered for each delayed group
! to zero.
nu_d(:) = 0
p % fission = .true. ! Fission neutrons will be banked
do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4)
! Bank source neutrons by copying particle data
bank_array(i) % xyz = p % coord(1) % xyz
! Set particle as neutron
bank_array(i) % particle = NEUTRON
! Set weight of fission bank site
bank_array(i) % wgt = ONE/weight
! Sample cosine of angle -- fission neutrons are treated as being emitted
! isotropically.
mu = TWO * prn() - ONE
! Sample azimuthal angle uniformly in [0,2*pi)
phi = TWO * PI * prn()
bank_array(i) % uvw(1) = mu
bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
call sample_fission_energy_c(p % material, p % g, dg, gout)
bank_array(i) % E = real(gout, 8)
bank_array(i) % delayed_group = dg
! Set delayed group on particle too
p % delayed_group = dg
! Increment the number of neutrons born delayed
if (p % delayed_group > 0) then
nu_d(p % delayed_group) = nu_d(p % delayed_group) + 1
end if
end do
! increment number of bank sites
size_bank = min(size_bank + nu, int(size(bank_array),8))
! Store total weight banked for analog fission tallies
p % n_bank = nu
p % wgt_bank = nu/weight
p % n_delayed_bank(:) = nu_d(:)
end subroutine create_fission_sites
end module physics_mg

234
src/physics_mg.cpp Normal file
View file

@ -0,0 +1,234 @@
#include "openmc/physics_mg.h"
#include <stdexcept>
#include <sstream>
#include "xtensor/xarray.hpp"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/math_functions.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/physics_common.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/tally.h"
namespace openmc {
void
collision_mg(Particle* p, const double* energy_bin_avg,
const MaterialMacroXS* material_xs)
{
// Add to the collision counter for the particle
p->n_collision++;
// Sample the reaction type
sample_reaction(p, energy_bin_avg, material_xs);
// Display information about collision
if ((settings::verbosity >= 10) || (simulation::trace)) {
std::stringstream msg;
msg << " Energy Group = " << p->g;
write_message(msg, 1);
}
}
void
sample_reaction(Particle* p, const double* energy_bin_avg,
const MaterialMacroXS* material_xs)
{
// Create fission bank sites. Note that while a fission reaction is sampled,
// it never actually "happens", i.e. the weight of the particle does not
// change when sampling fission sites. The following block handles all
// absorption (including fission)
if (materials[p->material - 1]->fissionable) {
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
Bank* result_bank;
int64_t result_bank_size;
// Get pointer to fission bank from Fortran side
openmc_fission_bank(&result_bank, &result_bank_size);
create_fission_sites(p, result_bank, &n_bank, result_bank_size,
material_xs);
} else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
(settings::create_fission_neutrons)) {
create_fission_sites(p, p->secondary_bank, &(p->n_secondary),
MAX_SECONDARY, material_xs);
}
}
// If survival biasing is being used, the following subroutine adjusts the
// weight of the particle. Otherwise, it checks to see if absorption occurs.
if (material_xs->absorption > 0.) {
absorption(p, material_xs);
} else {
p->absorb_wgt = 0.;
}
if (!p->alive) return;
// Sample a scattering event to determine the energy of the exiting neutron
scatter(p, energy_bin_avg);
// Play Russian roulette if survival biasing is turned on
if (settings::survival_biasing) {
russian_roulette(p);
if (!p->alive) return;
}
}
void
scatter(Particle* p, const double* energy_bin_avg)
{
// Adjust indices for Fortran to C++ indexing
// TODO: Remove when no longer needed
int gin = p->last_g - 1;
int gout = p->g - 1;
int i_mat = p->material - 1;
macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt);
// Adjust return value for fortran indexing
// TODO: Remove when no longer needed
p->g = gout + 1;
// Rotate the angle
rotate_angle_c(p->coord[0].uvw, p->mu, nullptr);
// Update energy value for downstream compatability (in tallying)
p->E = energy_bin_avg[gout];
// Set event component
p->event = EVENT_SCATTER;
}
void
create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
int64_t bank_array_size, const MaterialMacroXS* material_xs)
{
// TODO: Heat generation from fission
// If uniform fission source weighting is turned on, we increase or decrease
// the expected number of fission sites produced
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
// Determine the expected number of neutrons produced
double nu_t = p->wgt / simulation::keff * weight * material_xs->nu_fission /
material_xs->total;
// Sample the number of neutrons produced
int nu = static_cast<int>(nu_t);
if (prn() <= (nu_t - int(nu_t))) {
nu++;
}
// Check for the bank size getting hit. For fixed source calculations, this
// is a fatal error; for eigenvalue calculations, it just means that k-eff
// was too high for a single batch.
if (*size_bank + nu > bank_array_size) {
if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
throw std::runtime_error{"Secondary particle bank size limit reached."
" If you are running a subcritical multiplication problem,"
" k-effective may be too close to one."};
} else {
if (mpi::master) {
std::stringstream msg;
msg << "Maximum number of sites in fission bank reached. This can"
" result in irreproducible results using different numbers of"
" processes/threads.";
warning(msg);
}
}
}
// Begin banking the source neutrons
// First, if our bank is full then don't continue
if ((nu == 0) || (*size_bank == bank_array_size)) return;
// Initialize the counter of delayed neutrons encountered for each delayed
// group.
double nu_d[MAX_DELAYED_GROUPS] = {0.};
p->fission = true;
for (size_t i = static_cast<size_t>(*size_bank);
i < static_cast<size_t>(std::min(*size_bank + nu, bank_array_size)); i++) {
// Bank source neutrons by copying the particle data
bank_array[i].xyz[0] = p->coord[0].xyz[0];
bank_array[i].xyz[1] = p->coord[0].xyz[1];
bank_array[i].xyz[2] = p->coord[0].xyz[2];
// Set that the bank particle is a neutron
bank_array[i].particle = static_cast<int>(ParticleType::neutron);
// Set the weight of the fission bank site
bank_array[i].wgt = 1. / weight;
// Sample the cosine of the angle, assuming fission neutrons are emitted
// isotropically
double mu = 2. * prn() - 1.;
// Sample the azimuthal angle uniformly in [0, 2.pi)
double phi = 2. * PI * prn();
bank_array[i].uvw[0] = mu;
bank_array[i].uvw[1] = std::sqrt(1. - mu * mu) * std::cos(phi);
bank_array[i].uvw[2] = std::sqrt(1. - mu * mu) * std::sin(phi);
// Sample secondary energy distribution for the fission reaction and set
// the energy in the fission bank
int dg;
int gout;
macro_xs[p->material - 1].sample_fission_energy(p->g - 1, dg, gout);
bank_array[i].E = static_cast<double>(gout + 1);
bank_array[i].delayed_group = dg + 1;
// Set the delayed group on the particle as well
p->delayed_group = dg + 1;
// Increment the number of neutrons born delayed
if (p->delayed_group > 0) {
nu_d[dg]++;
}
}
// Increment number of bank sites
*size_bank = std::min(*size_bank + nu, bank_array_size);
// Store the total weight banked for analog fission tallies
p->n_bank = nu;
p->wgt_bank = nu / weight;
for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
p->n_delayed_bank[d] = nu_d[d];
}
}
void
absorption(Particle* p, const MaterialMacroXS* material_xs)
{
if (settings::survival_biasing) {
// Determine weight absorbed in survival biasing
p->absorb_wgt = p->wgt * material_xs->absorption / material_xs->total;
// Adjust weight of particle by the probability of absorption
p->wgt -= p->absorb_wgt;
p->last_wgt = p->wgt;
// Score implicit absorpion estimate of keff
#pragma omp atomic
global_tally_absorption += p->absorb_wgt * material_xs->nu_fission /
material_xs->absorption;
} else {
if (material_xs->absorption > prn() * material_xs->total) {
#pragma omp atomic
global_tally_absorption += p->wgt * material_xs->nu_fission /
material_xs->absorption;
p->alive = false;
p->event = EVENT_ABSORB;
}
}
}
} //namespace openmc

12
src/scattdata.cpp
View file

@ -316,8 +316,8 @@ ScattDataLegendre::update_max_val()
}
// Calculate probability
double f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
dist[gin][i_gout].data(), mu);
double f = evaluate_legendre(dist[gin][i_gout].size() - 1,
dist[gin][i_gout].data(), mu);
// if this is a new maximum, store it
if (f > max_val[gin][i_gout]) max_val[gin][i_gout] = f;
@ -339,8 +339,8 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu)
f = 0.;
} else {
int i_gout = gout - gmin[gin];
f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
dist[gin][i_gout].data(), mu);
f = evaluate_legendre(dist[gin][i_gout].size() - 1,
dist[gin][i_gout].data(), mu);
}
return f;
}
@ -881,8 +881,8 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
tab.fmu[gin][i_gout].resize(n_mu);
for (int imu = 0; imu < n_mu; imu++) {
tab.fmu[gin][i_gout][imu] =
evaluate_legendre_c(leg.dist[gin][i_gout].size() - 1,
leg.dist[gin][i_gout].data(), tab.mu[imu]);
evaluate_legendre(leg.dist[gin][i_gout].size() - 1,
leg.dist[gin][i_gout].data(), tab.mu[imu]);
}
// Ensure positivity

4
src/secondary_nbody.cpp
View file

@ -32,13 +32,13 @@ void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu) const
double E_max = (Ap - 1.0)/Ap * (A_/(A_ + 1.0)*E_in + Q_);
// x is essentially a Maxwellian distribution
double x = maxwell_spectrum_c(1.0);
double x = maxwell_spectrum(1.0);
double y;
double r1, r2, r3, r4, r5, r6;
switch (n_bodies_) {
case 3:
y = maxwell_spectrum_c(1.0);
y = maxwell_spectrum(1.0);
break;
case 4:
r1 = prn();

17
src/settings.cpp
View file

@ -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

136
src/simulation.F90
View file

@ -11,7 +11,7 @@ module simulation
use cmfd_header, only: cmfd_on
use constants, only: ZERO
use eigenvalue, only: calculate_average_keff, calculate_generation_keff, &
synchronize_bank, keff_generation, k_sum
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()
@ -435,7 +413,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.
@ -472,17 +450,6 @@ 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 openmc_simulation_finalize_c() bind(C)
@ -490,6 +457,9 @@ contains
subroutine print_overlap_check() bind(C)
end subroutine print_overlap_check
subroutine broadcast_results() bind(C)
end subroutine broadcast_results
end interface
err = 0
@ -513,37 +483,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
!===============================================================================

139
src/simulation.cpp
View file

@ -1,10 +1,14 @@
#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 "openmc/tallies/tally_filter.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
// calculation.
@ -29,10 +33,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()
@ -47,27 +80,17 @@ void openmc_simulation_init_c()
}
}
void calculate_work()
void initialize_generation()
{
// Determine minimum amount of particles to simulate on each processor
int64_t min_work = settings::n_particles / mpi::n_procs;
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
// Reset number of fission bank sites
n_bank = 0;
// Determine number of processors that have one extra particle
int64_t remainder = settings::n_particles % mpi::n_procs;
// Count source sites if using uniform fission source weighting
if (settings::ufs_on) ufs_count_sites();
int64_t i_bank = 0;
work_index.reserve(mpi::n_procs);
work_index.push_back(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;
// Set index into source bank for rank i
i_bank += work_i;
work_index.push_back(i_bank);
// Store current value of tracklength k
keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE);
}
}
@ -80,4 +103,82 @@ openmc_simulation_finalize_c()
}
}
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
int64_t min_work = settings::n_particles / mpi::n_procs;
// Determine number of processors that have one extra particle
int64_t remainder = settings::n_particles % mpi::n_procs;
int64_t i_bank = 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) simulation::work = work_i;
// Set index into source bank for rank i
i_bank += work_i;
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

86
src/simulation_header.F90
View file

@ -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

16
src/source.cpp
View file

@ -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

184
src/state_point.F90
View file

@ -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

147
src/state_point.cpp
View file

@ -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

66
src/tallies/tally.F90
View file

@ -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
View 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

13
src/tallies/tally_header.F90
View file

@ -137,7 +137,7 @@ module tally_header
! invalidation. Thus, we use threadprivate variables to accumulate global
! tallies and then reduce at the end of a generation.
real(C_DOUBLE), public :: global_tally_collision = ZERO
real(C_DOUBLE), public :: global_tally_absorption = ZERO
real(C_DOUBLE), public, bind(C) :: global_tally_absorption = ZERO
real(C_DOUBLE), public :: global_tally_tracklength = ZERO
real(C_DOUBLE), public :: global_tally_leakage = ZERO
!$omp threadprivate(global_tally_collision, global_tally_absorption, &
@ -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
@ -688,14 +688,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
View 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

22
src/timer_header.F90
View file

@ -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

15
src/tracking.F90
View file

@ -17,7 +17,6 @@ module tracking
use nuclide_header
use particle_header
use physics, only: collision
use physics_mg, only: collision_mg
use random_lcg, only: prn, prn_set_stream
use settings
use simulation_header
@ -33,6 +32,16 @@ module tracking
implicit none
interface
subroutine collision_mg(p, energy_bin_avg, material_xs) bind(C)
import Particle, C_DOUBLE, MaterialMacroXS
type(Particle), intent(inout) :: p
real(C_DOUBLE), intent(in) :: energy_bin_avg(*)
type(MaterialMacroXS), intent(in) :: material_xs
end subroutine collision_mg
end interface
contains
!===============================================================================
@ -226,7 +235,7 @@ contains
if (run_CE) then
call collision(p)
else
call collision_mg(p)
call collision_mg(p, energy_bin_avg, material_xs)
end if
! Score collision estimator tallies -- this is done after a collision
@ -313,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

35
src/volume_calc.F90
View file

@ -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