crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge remote-tracking branch 'upstream/develop' into cmfd-capi

Browse source
This commit is contained in:
Shikhar Kumar 2018-11-12 16:50:45 -05:00
commit c33e6b3dd9
220 changed files with 8685 additions and 8047 deletions
Download patch file Download diff file Expand all files Collapse all files

1
.gitignore vendored
View file

@ -69,6 +69,7 @@ scripts/G4EMLOW*/
# Images
*.ppm
*.voxel
*.vti
# PyCharm project configuration files
.idea

52
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
@ -132,7 +131,7 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel)
# Intel compiler options
list(APPEND f90flags -fpp -std08 -assume byterecl -traceback)
if(debug)
list(APPEND f90flags -g -warn -ftrapuv -fp-stack-check
list(APPEND f90flags -g -warn -fp-stack-check
"-check all" -fpe0 -O0)
list(APPEND ldflags -g)
endif()
@ -345,10 +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
src/pugixml/pugixml_f.F90
src/random_lcg.F90
src/reaction_header.F90
@ -375,26 +370,14 @@ add_library(libopenmc SHARED
src/tallies/tally_derivative_header.F90
src/tallies/tally_filter.F90
src/tallies/tally_filter_header.F90
src/tallies/tally_filter_azimuthal.F90
src/tallies/tally_filter_cell.F90
src/tallies/tally_filter_cellborn.F90
src/tallies/tally_filter_cellfrom.F90
src/tallies/tally_filter_delayedgroup.F90
src/tallies/tally_filter_distribcell.F90
src/tallies/tally_filter_energy.F90
src/tallies/tally_filter_energyfunc.F90
src/tallies/tally_filter_legendre.F90
src/tallies/tally_filter_material.F90
src/tallies/tally_filter_mesh.F90
src/tallies/tally_filter_meshsurface.F90
src/tallies/tally_filter_mu.F90
src/tallies/tally_filter_particle.F90
src/tallies/tally_filter_polar.F90
src/tallies/tally_filter_sph_harm.F90
src/tallies/tally_filter_sptl_legendre.F90
src/tallies/tally_filter_surface.F90
src/tallies/tally_filter_universe.F90
src/tallies/tally_filter_zernike.F90
src/tallies/tally_header.F90
src/tallies/trigger.F90
src/tallies/trigger_header.F90
@ -408,6 +391,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,8 +407,11 @@ 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/progress_bar.cpp
src/pugixml/pugixml_c.cpp
src/random_lcg.cpp
src/reaction.cpp
@ -441,9 +428,33 @@ add_library(libopenmc SHARED
src/string_functions.cpp
src/summary.cpp
src/surface.cpp
src/tallies/filter.cpp
src/tallies/filter_azimuthal.cpp
src/tallies/filter_cellborn.cpp
src/tallies/filter_cellfrom.cpp
src/tallies/filter_cell.cpp
src/tallies/filter_delayedgroup.cpp
src/tallies/filter_distribcell.cpp
src/tallies/filter_energyfunc.cpp
src/tallies/filter_energy.cpp
src/tallies/filter_legendre.cpp
src/tallies/filter_material.cpp
src/tallies/filter_mesh.cpp
src/tallies/filter_meshsurface.cpp
src/tallies/filter_mu.cpp
src/tallies/filter_particle.cpp
src/tallies/filter_polar.cpp
src/tallies/filter_sph_harm.cpp
src/tallies/filter_sptl_legendre.cpp
src/tallies/filter_surface.cpp
src/tallies/filter_universe.cpp
src/tallies/filter_zernike.cpp
src/tallies/tally.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
@ -470,9 +481,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

11
docs/source/io_formats/data_wmp.rst
View file

@ -4,10 +4,15 @@
Windowed Multipole Library Format
=================================
**/version** (*char[]*)
The format version of the file. The current version is "v1.0"
**/**
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file
- **version** (*int[2]*) -- Major and minor version of the data
**/<nuclide name>/**
:Datasets:
**/nuclide/**
- **broaden_poly** (*int[]*)
If 1, Doppler broaden curve fit for window with corresponding index.
If 0, do not.

9
docs/source/pythonapi/mgxs.rst
View file

@ -5,6 +5,15 @@
Energy Groups
-------------
Module Variables
++++++++++++++++
.. autodata:: openmc.mgxs.GROUP_STRUCTURES
:annotation:
Classes
+++++++
.. autosummary::
:toctree: generated
:nosignatures:

4
examples/jupyter/mgxs-part-ii.ipynb
View file

@ -2048,8 +2048,8 @@
"o16 = df[df['nuclide'] == 'O16']['mean']\n",
"\n",
"# Cast DataFrames as NumPy arrays\n",
"h1 = h1.as_matrix()\n",
"o16 = o16.as_matrix()\n",
"h1 = h1.values\n",
"o16 = o16.values\n",
"\n",
"# Reshape arrays to 2D matrix for plotting\n",
"h1.shape = (fine_groups.num_groups, fine_groups.num_groups)\n",

2
examples/jupyter/nuclear-data.ipynb
View file

@ -1466,7 +1466,7 @@
},
"outputs": [],
"source": [
"url = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/092238.h5'\n",
"url = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/092238.h5'\n",
"filename, headers = urllib.request.urlretrieve(url, '092238.h5')"
]
},

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

File diff suppressed because one or more lines are too long

13
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[]);
@ -94,7 +95,7 @@ extern "C" {
int openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]);
int openmc_sphharm_filter_set_order(int32_t index, int order);
int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
int openmc_statepoint_write(const char filename[]);
int openmc_statepoint_write(const char filename[], bool* write_source);
int openmc_tally_allocate(int32_t index, const char* type);
int openmc_tally_get_active(int32_t index, bool* active);
int openmc_tally_get_estimator(int32_t index, int32_t* estimator);
@ -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,10 +135,7 @@ 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;
extern int32_t n_materials;
extern int n_nuclides;
@ -148,15 +146,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
}

17
include/openmc/constants.h
View file

@ -62,6 +62,7 @@ constexpr int MAX_SAMPLE {100000};
// Maximum number of words in a single line, length of line, and length of
// single word
constexpr int MAX_LINE_LEN {250};
constexpr int MAX_WORD_LEN {150};
// Maximum number of external source spatial resamples to encounter before an
@ -314,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
@ -407,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};
@ -418,6 +420,7 @@ constexpr int DIFF_NUCLIDE_DENSITY {2};
constexpr int DIFF_TEMPERATURE {3};
constexpr int C_NONE {-1};
constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE
// Interpolation rules
enum class Interpolation {

17
include/openmc/container_util.h Normal file
View file

@ -0,0 +1,17 @@
#ifndef OPENMC_CONTAINER_UTIL_H
#define OPENMC_CONTAINER_UTIL_H
#include <algorithm> // for find
#include <iterator> // for begin, end
namespace openmc {
template<class C, class T>
inline bool contains(const C& v, const T& x)
{
return std::end(v) != std::find(std::begin(v), std::end(v), x);
}
}
#endif // OPENMC_CONTAINER_UTIL_H

34
include/openmc/eigenvalue.h
View file

@ -1,6 +1,10 @@
//! \file eigenvalue.h
//! \brief Data/functions related to k-eigenvalue calculations
#ifndef OPENMC_EIGENVALUE_H
#define OPENMC_EIGENVALUE_H
#include <array>
#include <cstdint> // for int64_t
#include <vector>
@ -14,6 +18,8 @@ namespace openmc {
// Global variables
//==============================================================================
extern double keff_generation; //!< Single-generation k on each processor
extern std::array<double, 2> k_sum; //!< Used to reduce sum and sum_sq
extern std::vector<double> entropy; //!< Shannon entropy at each generation
extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
@ -24,6 +30,34 @@ extern "C" int64_t n_bank;
// Non-member functions
//==============================================================================
//! Collect/normalize the tracklength keff from each process
extern "C" void calculate_generation_keff();
//! Calculate mean/standard deviation of keff during active generations
//!
//! This function sets the global variables keff and keff_std which represent
//! the mean and standard deviation of the mean of k-effective over active
//! generations. It also broadcasts the value from the master process.
extern "C" void calculate_average_keff();
//! Calculates a minimum variance estimate of k-effective
//!
//! The minimum variance estimate is based on a linear combination of the
//! collision, absorption, and tracklength estimates. The theory behind this can
//! be found in M. Halperin, "Almost linearly-optimum combination of unbiased
//! estimates," J. Am. Stat. Assoc., 56, 36-43 (1961),
//! doi:10.1080/01621459.1961.10482088. The implementation here follows that
//! described in T. Urbatsch et al., "Estimation and interpretation of keff
//! confidence intervals in MCNP," Nucl. Technol., 111, 169-182 (1995).
//!
//! \param[out] k_combined Estimate of k-effective and its standard deviation
//! \return Error status
extern "C" int openmc_get_keff(double* k_combined);
//! 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,8 @@
#ifndef OPENMC_FINALIZE_H
#define OPENMC_FINALIZE_H
extern "C" void openmc_free_bank();
namespace openmc {
} // namespace openmc
#endif // OPENMC_FINALIZE_H

49
include/openmc/geometry_aux.h
View file

@ -5,6 +5,7 @@
#define OPENMC_GEOMETRY_AUX_H
#include <cstdint>
#include <string>
namespace openmc {
@ -26,7 +27,7 @@ extern "C" void assign_temperatures();
//!
//! This function looks for a universe that is not listed in a Cell::fill or in
//! a Lattice.
//! @return The index of the root universe.
//! \return The index of the root universe.
//==============================================================================
extern "C" int32_t find_root_universe();
@ -41,14 +42,14 @@ extern "C" void neighbor_lists();
//! Populate all data structures needed for distribcells.
//==============================================================================
extern "C" void prepare_distribcell(int32_t* filter_cell_list, int n);
extern "C" void prepare_distribcell();
//==============================================================================
//! Recursively search through the geometry and count cell instances.
//!
//! This function will update the Cell::n_instances value for each cell in the
//! geometry.
//! @param univ_indx The index of the universe to begin searching from (probably
//! \param univ_indx The index of the universe to begin searching from (probably
//! the root universe).
//==============================================================================
@ -56,51 +57,33 @@ extern "C" void count_cell_instances(int32_t univ_indx);
//==============================================================================
//! Recursively search through universes and count universe instances.
//! @param search_univ The index of the universe to begin searching from.
//! @param target_univ_id The ID of the universe to be counted.
//! @return The number of instances of target_univ_id in the geometry tree under
//! \param search_univ The index of the universe to begin searching from.
//! \param target_univ_id The ID of the universe to be counted.
//! \return The number of instances of target_univ_id in the geometry tree under
//! search_univ.
//==============================================================================
extern "C" int
count_universe_instances(int32_t search_univ, int32_t target_univ_id);
//==============================================================================
//! Find the length necessary for a string to contain a distribcell path.
//! @param target_cell The index of the Cell in the global Cell array.
//! @param map The index of the distribcell mapping corresponding to the target
//! cell.
//! @param target_offset An instance number for a distributed cell.
//! @param root_univ The index of the root Universe in the global Universe
//! array.
//! @return The size of a character array needed to fit the distribcell path.
//==============================================================================
extern "C" int
distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset,
int32_t root_univ);
//==============================================================================
//! Build a character array representing the path to a distribcell instance.
//! @param target_cell The index of the Cell in the global Cell array.
//! @param map The index of the distribcell mapping corresponding to the target
//! \param target_cell The index of the Cell in the global Cell array.
//! \param map The index of the distribcell mapping corresponding to the target
//! cell.
//! @param target_offset An instance number for a distributed cell.
//! @param root_univ The index of the root Universe in the global Universe
//! array.
//! @param[out] path The unique traversal through the geometry tree that leads
//! to the desired instance of the target cell.
//! \param target_offset An instance number for a distributed cell.
//! \return The unique traversal through the geometry tree that leads to the
//! desired instance of the target cell.
//==============================================================================
extern "C" void
distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset,
int32_t root_univ, char *path);
std::string
distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset);
//==============================================================================
//! Determine the maximum number of nested coordinate levels in the geometry.
//! @param univ The index of the universe to begin seraching from (probably the
//! \param univ The index of the universe to begin seraching from (probably the
//! root universe).
//! @return The number of coordinate levels.
//! \return The number of coordinate levels.
//==============================================================================
extern "C" int maximum_levels(int32_t univ);

27
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"
@ -51,7 +52,7 @@ void write_string(hid_t group_id, const char* name, const std::string& buffer,
std::vector<hsize_t> attribute_shape(hid_t obj_id, const char* name);
std::vector<std::string> dataset_names(hid_t group_id);
void ensure_exists(hid_t group_id, const char* name);
void ensure_exists(hid_t obj_id, const char* name, bool attribute=false);
std::vector<std::string> group_names(hid_t group_id);
std::vector<hsize_t> object_shape(hid_t obj_id);
std::string object_name(hid_t obj_id);
@ -201,7 +202,13 @@ read_attribute(hid_t obj_id, const char* name, std::vector<std::string>& vec)
read_attr_string(obj_id, name, n, buffer[0]);
for (int i = 0; i < m; ++i) {
vec.emplace_back(&buffer[i][0], std::min(strlen(buffer[i]), n));
// Determine proper length of string -- strlen doesn't work because
// buffer[i] might not have any null characters
std::size_t k = 0;
for (; k < n; ++k) if (buffer[i][k] == '\0') break;
// Create string based on (char*, size_t) constructor
vec.emplace_back(&buffer[i][0], k);
}
}
@ -333,7 +340,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 +368,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,
@ -375,5 +386,11 @@ write_dataset(hid_t obj_id, const char* name, Position r)
write_dataset(obj_id, name, buffer);
}
inline void
write_dataset(hid_t obj_id, const char* name, std::string buffer)
{
write_string(obj_id, name, buffer.c_str(), false);
}
} // namespace openmc
#endif // OPENMC_HDF5_INTERFACE_H

3
include/openmc/initialize.h
View file

@ -5,9 +5,6 @@
#include "mpi.h"
#endif
extern "C" void print_usage();
extern "C" void print_version();
namespace openmc {
int parse_command_line(int argc, char* argv[]);

3
include/openmc/material.h
View file

@ -24,8 +24,9 @@ extern std::unordered_map<int32_t, int32_t> material_map;
class Material
{
public:
int32_t id; //!< Unique ID
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.
//!

12
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
@ -90,7 +90,7 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
//! evaluated at rho and phi.
//==============================================================================
extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
extern "C" void calc_zn(int n, double rho, double phi, double zn[]);
//==============================================================================
//! Calculate only the even radial components of n-th order modified Zernike
@ -113,7 +113,7 @@ extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
//! evaluated at rho and phi when m = 0.
//==============================================================================
extern "C" void calc_zn_rad_c(int n, double rho, double zn_rad[]);
extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]);
//==============================================================================
//! Rotate the direction cosines through a polar angle whose cosine is mu and
@ -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 {

13
include/openmc/mgxs_interface.h
View file

@ -7,6 +7,7 @@
#include "hdf5_interface.h"
#include "mgxs.h"
#include <vector>
namespace openmc {
@ -17,6 +18,9 @@ namespace openmc {
extern std::vector<Mgxs> nuclides_MG;
extern std::vector<Mgxs> macro_xs;
extern "C" int num_energy_groups;
extern std::vector<double> energy_bins;
extern std::vector<double> energy_bin_avg;
extern std::vector<double> rev_energy_bins;
//==============================================================================
// Mgxs data loading interface methods
@ -36,6 +40,8 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
int n_temps, const double temps[], const double atom_densities[],
double tolerance, int& method);
extern "C" void read_mg_cross_sections_header_c(hid_t file_id);
//==============================================================================
// Mgxs tracking/transport/tallying interface methods
//==============================================================================
@ -44,13 +50,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);

18
include/openmc/output.h
View file

@ -4,6 +4,7 @@
#ifndef OPENMC_OUTPUT_H
#define OPENMC_OUTPUT_H
#include <string>
namespace openmc {
@ -16,11 +17,26 @@ namespace openmc {
void header(const char* msg, int level);
//==============================================================================
//! Retrieve a time stamp
//!
//! \return current time stamp (format: "yyyy-mm-dd hh:mm:ss")
//==============================================================================
std::string time_stamp();
//==============================================================================
//! Display plot information
//==============================================================================
extern "C" void print_plot();
//==============================================================================
//! Display information regarding cell overlap checking.
//==============================================================================
extern "C" void print_overlap_check();
void print_overlap_check();
extern "C" void title();

7
include/openmc/particle.h
View file

@ -146,9 +146,7 @@ extern "C" {
//! site may have been produced from an external source, from fission, or
//! simply as a secondary particle.
//! \param src Source site data
//! \param run_CE Whether continuous-energy data is being used
//! \param energy_bin_avg An array of energy group bin averages
void from_source(const Bank* src, bool run_CE, const double* energy_bin_avg);
void from_source(const Bank* src);
//! mark a particle as lost and create a particle restart file
//! \param message A warning message to display
@ -174,8 +172,7 @@ extern "C" {
void particle_create_secondary(Particle* p, const double* uvw, double E,
int type, bool run_CE);
void particle_initialize(Particle* p);
void particle_from_source(Particle* p, const Bank* src, bool run_CE,
const double* energy_bin_avg);
void particle_from_source(Particle* p, const Bank* src);
void particle_mark_as_lost(Particle* p, const char* message);
void particle_write_restart(Particle* p);

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

172
include/openmc/plot.h
View file

@ -1,15 +1,181 @@
#ifndef OPENMC_PLOT_H
#define OPENMC_PLOT_H
#include <unordered_map>
#include <sstream>
#include "xtensor/xarray.hpp"
#include "hdf5.h"
#include "openmc/position.h"
#include "openmc/constants.h"
#include "openmc/particle.h"
#include "openmc/xml_interface.h"
namespace openmc {
extern "C" void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
//===============================================================================
// Global variables
//===============================================================================
extern int PLOT_LEVEL_LOWEST; //!< lower bound on plot universe level
extern std::unordered_map<int, int> plot_map; //!< map of plot ids to index
extern "C" int32_t n_plots; //!< number of plots in openmc run
class Plot;
extern std::vector<Plot> plots; //!< Plot instance container
//===============================================================================
// RGBColor holds color information for plotted objects
//===============================================================================
struct RGBColor {
//Constructors
RGBColor() : red(0), green(0), blue(0) { };
RGBColor(const int v[3]) : red(v[0]), green(v[1]), blue(v[2]) { };
RGBColor(int r, int g, int b) : red(r), green(g), blue(b) { };
RGBColor(const std::vector<int> &v) {
if (v.size() != 3) {
throw std::out_of_range("Incorrect vector size for RGBColor.");
}
red = v[0];
green = v[1];
blue = v[2];
}
// Members
uint8_t red, green, blue;
};
typedef xt::xtensor<RGBColor, 2> ImageData;
enum class PlotType {
slice = 1,
voxel = 2
};
enum class PlotBasis {
xy = 1,
xz = 2,
yz = 3
};
enum class PlotColorBy {
cells = 1,
mats = 2
};
//===============================================================================
// Plot class
//===============================================================================
class Plot
{
public:
// Constructor
Plot(pugi::xml_node plot);
// Methods
private:
void set_id(pugi::xml_node plot_node);
void set_type(pugi::xml_node plot_node);
void set_output_path(pugi::xml_node plot_node);
void set_bg_color(pugi::xml_node plot_node);
void set_basis(pugi::xml_node plot_node);
void set_origin(pugi::xml_node plot_node);
void set_width(pugi::xml_node plot_node);
void set_universe(pugi::xml_node plot_node);
void set_default_colors(pugi::xml_node plot_node);
void set_user_colors(pugi::xml_node plot_node);
void set_meshlines(pugi::xml_node plot_node);
void set_mask(pugi::xml_node plot_node);
// Members
public:
int id_; //!< Plot ID
PlotType type_; //!< Plot type (Slice/Voxel)
PlotColorBy color_by_; //!< Plot coloring (cell/material)
Position origin_; //!< Plot origin in geometry
Position width_; //!< Plot width in geometry
PlotBasis basis_; //!< Plot basis (XY/XZ/YZ)
std::array<int, 3> pixels_; //!< Plot size in pixels
int meshlines_width_; //!< Width of lines added to the plot
int level_; //!< Plot universe level
int index_meshlines_mesh_; //!< Index of the mesh to draw on the plot
RGBColor meshlines_color_; //!< Color of meshlines on the plot
RGBColor not_found_; //!< Plot background color
std::vector<RGBColor> colors_; //!< Plot colors
std::string path_plot_; //!< Plot output filename
};
//===============================================================================
// Non-member functions
//===============================================================================
//! Add mesh lines to image data of a plot object
//! \param[in] plot object
//! \param[out] image data associated with the plot object
void draw_mesh_lines(Plot pl, ImageData& data);
//! Write a ppm image to file using a plot object's image data
//! \param[in] plot object
//! \param[out] image data associated with the plot object
void output_ppm(Plot pl,
const ImageData& data);
//! Get the rgb color for a given particle position in a plot
//! \param[in] particle with position for current pixel
//! \param[in] plot object
//! \param[out] rgb color
//! \param[out] cell or material id for particle position
void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id);
//! Initialize a voxel file
//! \param[in] id of an open hdf5 file
//! \param[in] dimensions of the voxel file (dx, dy, dz)
//! \param[out] dataspace pointer to voxel data
//! \param[out] dataset pointer to voxesl data
//! \param[out] pointer to memory space of voxel data
void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
hid_t* dset, hid_t* memspace);
extern "C" void voxel_write_slice(int x, hid_t dspace, hid_t dset,
//! Write a section of the voxel data to hdf5
//! \param[in] voxel slice
//! \param[out] dataspace pointer to voxel data
//! \param[out] dataset pointer to voxesl data
//! \param[out] pointer to data to write
void voxel_write_slice(int x, hid_t dspace, hid_t dset,
hid_t memspace, void* buf);
extern "C" void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
//! Close voxel file entities
//! \param[in] data space to close
//! \param[in] dataset to close
//! \param[in] memory space to close
void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
//===============================================================================
// External functions
//===============================================================================
//! Read plots from plots.xml node
//! \param[in] plot node of plots.xml
extern "C" void read_plots(pugi::xml_node* plot_node);
//! Create a ppm image for a plot object
//! \param[in] plot object
void create_ppm(Plot pl);
//! Create an hdf5 voxel file for a plot object
//! \param[in] plot object
void create_voxel(Plot pl);
//! Create a randomly generated RGB color
//! \return RGBColor with random value
RGBColor random_color();
} // namespace openmc
#endif // OPENMC_PLOT_H

22
include/openmc/progress_bar.h Normal file
View file

@ -0,0 +1,22 @@
#ifndef OPENMC_PROGRESSBAR_H
#define OPENMC_PROGRESSBAR_H
#include <string>
class ProgressBar {
public:
// Constructor
ProgressBar();
void set_value(double val);
private:
std::string bar;
char bar_old[72] = "???% | |";
};
#endif // OPENMC_PROGRESSBAR_H

1
include/openmc/random_lcg.h
View file

@ -17,6 +17,7 @@ extern "C" const int STREAM_SOURCE;
extern "C" const int STREAM_URR_PTABLE;
extern "C" const int STREAM_VOLUME;
extern "C" const int STREAM_PHOTON;
constexpr int64_t DEFAULT_SEED = 1;
//==============================================================================
//! Generate a pseudo-random number using a linear congruential generator.

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

13
include/openmc/settings.h
View file

@ -7,6 +7,8 @@
#include <array>
#include <cstdint>
#include <string>
#include <unordered_set>
#include <vector>
#include "pugixml.hpp"
@ -59,30 +61,33 @@ extern std::string path_statepoint; //!< path to a statepoint file
extern "C" int32_t index_entropy_mesh; //!< Index of entropy mesh in global mesh array
extern "C" int32_t index_ufs_mesh; //!< Index of UFS mesh in global mesh array
extern "C" int32_t index_cmfd_mesh; //!< Index of CMFD mesh in global mesh array
extern "C" int32_t n_batches; //!< number of (inactive+active) batches
extern "C" int32_t n_inactive; //!< number of inactive batches
extern "C" int32_t gen_per_batch; //!< number of generations per batch
extern "C" int64_t n_particles; //!< number of particles per generation
extern "C" int electron_treatment; //!< how to treat secondary electrons
extern "C" double energy_cutoff[4]; //!< Energy cutoff in [eV] for each particle type
extern "C" std::array<double, 4> energy_cutoff; //!< Energy cutoff in [eV] for each particle type
extern "C" int legendre_to_tabular_points; //!< number of points to convert Legendres
extern "C" int max_order; //!< Maximum Legendre order for multigroup data
extern "C" int n_log_bins; //!< number of bins for logarithmic energy grid
extern "C" int n_max_batches; //!< Maximum number of batches
extern "C" int res_scat_method; //!< resonance upscattering method
extern "C" double res_scat_energy_min; //!< Min energy in [eV] for res. upscattering
extern "C" double res_scat_energy_max; //!< Max energy in [eV] for res. upscattering
extern "C" int run_mode; //!< Run mode (eigenvalue, fixed src, etc.)
extern std::unordered_set<int> sourcepoint_batch; //!< Batches when source should be written
extern std::unordered_set<int> statepoint_batch; //!< Batches when state should be written
extern "C" int temperature_method; //!< method for choosing temperatures
extern "C" double temperature_tolerance; //!< Tolerance in [K] on choosing temperatures
extern "C" double temperature_default; //!< Default T in [K]
extern "C" double temperature_range[2]; //!< Min/max T in [K] over which to load xs
extern "C" std::array<double, 2> temperature_range; //!< Min/max T in [K] over which to load xs
extern "C" int trace_batch; //!< Batch to trace particle on
extern "C" int trace_gen; //!< Generation to trace particle on
extern "C" int64_t trace_particle; //!< Particle ID to enable trace on
extern std::vector<std::array<int, 3>> track_identifiers; //!< Particle numbers for writing tracks
extern "C" int trigger_batch_interval; //!< Batch interval for triggers
extern "C" int verbosity; //!< How verbose to make output
extern "C" double weight_cutoff; //!< Weight cutoff for Russian roulette

72
include/openmc/simulation.h
View file

@ -4,36 +4,86 @@
#ifndef OPENMC_SIMULATION_H
#define OPENMC_SIMULATION_H
#include "openmc/particle.h"
#include <cstdint>
#include <vector>
namespace openmc {
constexpr int STATUS_EXIT_NORMAL {0};
constexpr int STATUS_EXIT_MAX_BATCH {1};
constexpr int STATUS_EXIT_ON_TRIGGER {2};
//==============================================================================
// 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" bool initialized; //!< has simulation been initialized?
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" 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
//==============================================================================
//! Initialize simulation
extern "C" void openmc_simulation_init_c();
//! Determine number of particles to transport per process
void calculate_work();
//! Initialize a batch
void initialize_batch();
//! Initialize a fission generation
void initialize_generation();
void initialize_history(Particle* p, int64_t index_source);
//! Finalize a batch
//!
//! Handles synchronization and accumulation of tallies, calculation of Shannon
//! entropy, getting single-batch estimate of keff, and turning on tallies when
//! appropriate
void finalize_batch();
//! Finalize a fission generation
void finalize_generation();
//! Determine overall generation number
extern "C" int overall_generation();
#ifdef OPENMC_MPI
void broadcast_results();
#endif
} // namespace openmc
#endif // OPENMC_SIMULATION_H

3
include/openmc/source.h
View file

@ -58,6 +58,9 @@ extern "C" void initialize_source();
//! \return Sampled source site
extern "C" Bank sample_external_source();
//! Fill source bank at end of generation for fixed source simulations
void fill_source_bank_fixedsource();
} // namespace openmc
#endif // OPENMC_SOURCE_H

9
include/openmc/state_point.h
View file

@ -9,10 +9,11 @@
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);
void write_source_point(const char* filename);
extern "C" void write_source_bank(hid_t group_id, Bank* source_bank);
extern "C" void read_source_bank(hid_t group_id, Bank* source_bank);
extern "C" void write_tally_results_nr(hid_t file_id);
extern "C" void restart_set_keff();
} // namespace openmc
#endif // OPENMC_STATE_POINT_H

2
include/openmc/string_functions.h
View file

@ -14,5 +14,7 @@ char* strtrim(char* c_str);
void to_lower(std::string& str);
int word_count(std::string const& str);
} // namespace openmc
#endif // STRING_FUNCTIONS_H

99
include/openmc/tallies/filter.h Normal file
View file

@ -0,0 +1,99 @@
#ifndef OPENMC_TALLIES_FILTER_H
#define OPENMC_TALLIES_FILTER_H
#include <cstdint>
#include <memory>
#include <string>
#include <vector>
#include "openmc/hdf5_interface.h"
#include "openmc/particle.h"
#include "pugixml.hpp"
namespace openmc {
//==============================================================================
//! Stores bins and weights for filtered tally events.
//==============================================================================
class FilterMatch
{
public:
std::vector<int> bins_;
std::vector<double> weights_;
};
} // namespace openmc
// Without an explicit instantiation of vector<FilterMatch>, the Intel compiler
// will complain about the threadprivate directive on filter_matches. Note that
// this has to happen *outside* of the openmc namespace
template class std::vector<openmc::FilterMatch>;
namespace openmc {
//==============================================================================
//! Modifies tally score events.
//==============================================================================
class Filter
{
public:
virtual ~Filter() = default;
virtual std::string type() const = 0;
//! Uses an XML input to fill the filter's data fields.
virtual void from_xml(pugi::xml_node node) = 0;
//! Matches a tally event to a set of filter bins and weights.
//!
//! \param[out] match will contain the matching bins and corresponding
//! weights; note that there may be zero matching bins
virtual void
get_all_bins(const Particle* p, int estimator, FilterMatch& match) const = 0;
//! Writes data describing this filter to an HDF5 statepoint group.
virtual void
to_statepoint(hid_t filter_group) const
{
write_dataset(filter_group, "type", type());
write_dataset(filter_group, "n_bins", n_bins_);
}
//! Return a string describing a filter bin for the tallies.out file.
//
//! For example, an `EnergyFilter` might return the string
//! "Incoming Energy [0.625E-6, 20.0)".
virtual std::string text_label(int bin) const = 0;
virtual void initialize() {}
int n_bins_;
};
//==============================================================================
// Global variables
//==============================================================================
extern "C" int32_t n_filters;
extern std::vector<FilterMatch> filter_matches;
#pragma omp threadprivate(filter_matches)
extern std::vector<std::unique_ptr<Filter>> tally_filters;
//==============================================================================
extern "C" void free_memory_tally_c();
//==============================================================================
// Filter-related Fortran functions that will be called from C++
extern "C" int verify_filter(int32_t index);
extern "C" Filter* filter_from_f(int32_t index);
extern "C" void filter_update_n_bins(int32_t index);
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_H

35
include/openmc/tallies/filter_azimuthal.h Normal file
View file

@ -0,0 +1,35 @@
#ifndef OPENMC_TALLIES_FILTER_AZIMUTHAL_H
#define OPENMC_TALLIES_FILTER_AZIMUTHAL_H
#include <string>
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Bins the incident neutron azimuthal angle (relative to the global xy-plane).
//==============================================================================
class AzimuthalFilter : public Filter
{
public:
~AzimuthalFilter() = default;
std::string type() const override {return "azimuthal";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
std::vector<double> bins_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_AZIMUTHAL_H

42
include/openmc/tallies/filter_cell.h Normal file
View file

@ -0,0 +1,42 @@
#ifndef OPENMC_TALLIES_FILTER_CELL_H
#define OPENMC_TALLIES_FILTER_CELL_H
#include <cstdint>
#include <unordered_map>
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Specifies which geometric cells tally events reside in.
//==============================================================================
class CellFilter : public Filter
{
public:
~CellFilter() = default;
std::string type() const override {return "cell";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
//! The indices of the cells binned by this filter.
std::vector<int32_t> cells_;
//! A map from cell indices to filter bin indices.
std::unordered_map<int32_t, int> map_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_CELL_H

26
include/openmc/tallies/filter_cellborn.h Normal file
View file

@ -0,0 +1,26 @@
#ifndef OPENMC_TALLIES_FILTER_CELLBORN_H
#define OPENMC_TALLIES_FILTER_CELLBORN_H
#include <string>
#include "openmc/tallies/filter_cell.h"
namespace openmc {
//==============================================================================
//! Specifies which cell the particle was born in.
//==============================================================================
class CellbornFilter : public CellFilter
{
public:
std::string type() const override {return "cellborn";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
std::string text_label(int bin) const override;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_CELLBORN_H

26
include/openmc/tallies/filter_cellfrom.h Normal file
View file

@ -0,0 +1,26 @@
#ifndef OPENMC_TALLIES_FILTER_CELLFROM_H
#define OPENMC_TALLIES_FILTER_CELLFROM_H
#include <string>
#include "openmc/tallies/filter_cell.h"
namespace openmc {
//==============================================================================
//! Specifies which geometric cells particles exit when crossing a surface.
//==============================================================================
class CellFromFilter : public CellFilter
{
public:
std::string type() const override {return "cellfrom";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
std::string text_label(int bin) const override;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_CELLFROM_H

37
include/openmc/tallies/filter_delayedgroup.h Normal file
View file

@ -0,0 +1,37 @@
#ifndef OPENMC_TALLIES_FILTER_DELAYEDGROUP_H
#define OPENMC_TALLIES_FILTER_DELAYEDGROUP_H
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Bins outgoing fission neutrons in their delayed groups.
//!
//! The get_all_bins functionality is not actually used. The bins are manually
//! iterated over in the scoring subroutines.
//==============================================================================
class DelayedGroupFilter : public Filter
{
public:
~DelayedGroupFilter() = default;
std::string type() const override {return "delayedgroup";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
std::vector<int> groups_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_DELAYEDGROUP_H

36
include/openmc/tallies/filter_distribcell.h Normal file
View file

@ -0,0 +1,36 @@
#ifndef OPENMC_TALLIES_FILTER_DISTRIBCELL_H
#define OPENMC_TALLIES_FILTER_DISTRIBCELL_H
#include <string>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Specifies which distributed geometric cells tally events reside in.
//==============================================================================
class DistribcellFilter : public Filter
{
public:
~DistribcellFilter() = default;
std::string type() const override {return "distribcell";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
int32_t cell_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_DISTRIBCELL_H

55
include/openmc/tallies/filter_energy.h Normal file
View file

@ -0,0 +1,55 @@
#ifndef OPENMC_TALLIES_FILTER_ENERGY_H
#define OPENMC_TALLIES_FILTER_ENERGY_H
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Bins the incident neutron energy.
//==============================================================================
class EnergyFilter : public Filter
{
public:
~EnergyFilter() = default;
std::string type() const override {return "energy";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
std::vector<double> bins_;
//! True if transport group number can be used directly to get bin number
bool matches_transport_groups_ {false};
};
//==============================================================================
//! Bins the outgoing neutron energy.
//!
//! Only scattering events use the get_all_bins functionality. Nu-fission
//! tallies manually iterate over the filter bins.
//==============================================================================
class EnergyoutFilter : public EnergyFilter
{
public:
std::string type() const override {return "energyout";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
std::string text_label(int bin) const override;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_ENERGY_H

45
include/openmc/tallies/filter_energyfunc.h Normal file
View file

@ -0,0 +1,45 @@
#ifndef OPENMC_TALLIES_FILTER_ENERGYFUNC_H
#define OPENMC_TALLIES_FILTER_ENERGYFUNC_H
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Multiplies tally scores by an arbitrary function of incident energy
//! described by a piecewise linear-linear interpolation.
//==============================================================================
class EnergyFunctionFilter : public Filter
{
public:
EnergyFunctionFilter()
: Filter {}
{
n_bins_ = 1;
}
~EnergyFunctionFilter() = default;
std::string type() const override {return "energyfunction";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
//! Incident neutron energy interpolation grid.
std::vector<double> energy_;
//! Interpolant values.
std::vector<double> y_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_ENERGYFUNC_H

34
include/openmc/tallies/filter_legendre.h Normal file
View file

@ -0,0 +1,34 @@
#ifndef OPENMC_TALLIES_FILTER_LEGENDRE_H
#define OPENMC_TALLIES_FILTER_LEGENDRE_H
#include <string>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Gives Legendre moments of the change in scattering angle
//==============================================================================
class LegendreFilter : public Filter
{
public:
~LegendreFilter() = default;
std::string type() const override {return "legendre";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
int order_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_LEGENDRE_H

42
include/openmc/tallies/filter_material.h Normal file
View file

@ -0,0 +1,42 @@
#ifndef OPENMC_TALLIES_FILTER_MATERIAL_H
#define OPENMC_TALLIES_FILTER_MATERIAL_H
#include <cstdint>
#include <unordered_map>
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Specifies which material tally events reside in.
//==============================================================================
class MaterialFilter : public Filter
{
public:
~MaterialFilter() = default;
std::string type() const override {return "material";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
//! The indices of the materials binned by this filter.
std::vector<int32_t> materials_;
//! A map from material indices to filter bin indices.
std::unordered_map<int32_t, int> map_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_MATERIAL_H

41
include/openmc/tallies/filter_mesh.h Normal file
View file

@ -0,0 +1,41 @@
#ifndef OPENMC_TALLIES_FILTER_MESH_H
#define OPENMC_TALLIES_FILTER_MESH_H
#include <cstdint>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Indexes the location of particle events to a regular mesh. For tracklength
//! tallies, it will produce multiple valid bins and the bin weight will
//! correspond to the fraction of the track length that lies in that bin.
//==============================================================================
class MeshFilter : public Filter
{
public:
~MeshFilter() = default;
std::string type() const override {return "mesh";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
virtual int32_t mesh() const {return mesh_;}
virtual void set_mesh(int32_t mesh);
protected:
int32_t mesh_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_MESH_H

22
include/openmc/tallies/filter_meshsurface.h Normal file
View file

@ -0,0 +1,22 @@
#ifndef OPENMC_TALLIES_FILTER_MESHSURFACE_H
#define OPENMC_TALLIES_FILTER_MESHSURFACE_H
#include "openmc/tallies/filter_mesh.h"
namespace openmc {
class MeshSurfaceFilter : public MeshFilter
{
public:
std::string type() const override {return "meshsurface";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
std::string text_label(int bin) const override;
void set_mesh(int32_t mesh) override;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_MESHSURFACE_H

35
include/openmc/tallies/filter_mu.h Normal file
View file

@ -0,0 +1,35 @@
#ifndef OPENMC_TALLIES_FILTER_MU_H
#define OPENMC_TALLIES_FILTER_MU_H
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Bins the incoming-outgoing direction cosine. This is only used for scatter
//! reactions.
//==============================================================================
class MuFilter : public Filter
{
public:
~MuFilter() = default;
std::string type() const override {return "mu";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
std::vector<double> bins_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_MU_H

34
include/openmc/tallies/filter_particle.h Normal file
View file

@ -0,0 +1,34 @@
#ifndef OPENMC_TALLIES_FILTER_PARTICLE_H
#define OPENMC_TALLIES_FILTER_PARTICLE_H
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Bins by type of particle (e.g. neutron, photon).
//==============================================================================
class ParticleFilter : public Filter
{
public:
~ParticleFilter() = default;
std::string type() const override {return "particle";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
std::vector<int> particles_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_PARTICLE_H

35
include/openmc/tallies/filter_polar.h Normal file
View file

@ -0,0 +1,35 @@
#ifndef OPENMC_TALLIES_FILTER_POLAR_H
#define OPENMC_TALLIES_FILTER_POLAR_H
#include <cmath>
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Bins the incident neutron polar angle (relative to the global z-axis).
//==============================================================================
class PolarFilter : public Filter
{
public:
~PolarFilter() = default;
std::string type() const override {return "polar";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
std::vector<double> bins_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_POLAR_H

42
include/openmc/tallies/filter_sph_harm.h Normal file
View file

@ -0,0 +1,42 @@
#ifndef OPENMC_TALLIES_FILTER_SPH_HARM_H
#define OPENMC_TALLIES_FILTER_SPH_HARM_H
#include <string>
#include "openmc/tallies/filter.h"
namespace openmc {
//TODO: those integer values are not needed when Fortran interop is removed
enum class SphericalHarmonicsCosine {
scatter = 1, particle = 2
};
//==============================================================================
//! Gives spherical harmonics expansion moments of a tally score
//==============================================================================
class SphericalHarmonicsFilter : public Filter
{
public:
~SphericalHarmonicsFilter() = default;
std::string type() const override {return "sphericalharmonics";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
int order_;
//! The type of angle that this filter measures when binning events.
SphericalHarmonicsCosine cosine_ {SphericalHarmonicsCosine::particle};
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_SPH_HARM_H

48
include/openmc/tallies/filter_sptl_legendre.h Normal file
View file

@ -0,0 +1,48 @@
#ifndef OPENMC_TALLIES_FILTER_SPTL_LEGENDRE_H
#define OPENMC_TALLIES_FILTER_SPTL_LEGENDRE_H
#include <string>
#include "openmc/tallies/filter.h"
namespace openmc {
//TODO: those integer values are not needed when Fortran interop is removed
enum class LegendreAxis {
x = 1, y = 2, z = 3
};
//==============================================================================
//! Gives Legendre moments of the particle's normalized position along an axis
//==============================================================================
class SpatialLegendreFilter : public Filter
{
public:
~SpatialLegendreFilter() = default;
std::string type() const override {return "spatiallegendre";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
int order_;
//! The Cartesian coordinate axis that the Legendre expansion is applied to.
LegendreAxis axis_;
//! The minimum coordinate along the reference axis that the expansion covers.
double min_;
//! The maximum coordinate along the reference axis that the expansion covers.
double max_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_SPTL_LEGENDRE_H

42
include/openmc/tallies/filter_surface.h Normal file
View file

@ -0,0 +1,42 @@
#ifndef OPENMC_TALLIES_FILTER_SURFACE_H
#define OPENMC_TALLIES_FILTER_SURFACE_H
#include <cstdint>
#include <unordered_map>
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Specifies which surface particles are crossing
//==============================================================================
class SurfaceFilter : public Filter
{
public:
~SurfaceFilter() = default;
std::string type() const override {return "surface";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
//! The indices of the surfaces binned by this filter.
std::vector<int32_t> surfaces_;
//! A map from surface indices to filter bin indices.
std::unordered_map<int32_t, int> map_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_SURFACE_H

42
include/openmc/tallies/filter_universe.h Normal file
View file

@ -0,0 +1,42 @@
#ifndef OPENMC_TALLIES_FILTER_UNIVERSE_H
#define OPENMC_TALLIES_FILTER_UNIVERSE_H
#include <cstdint>
#include <unordered_map>
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Specifies which geometric universes tally events reside in.
//==============================================================================
class UniverseFilter : public Filter
{
public:
~UniverseFilter() = default;
std::string type() const override {return "universe";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
//! The indices of the universes binned by this filter.
std::vector<int32_t> universes_;
//! A map from universe indices to filter bin indices.
std::unordered_map<int32_t, int> map_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_UNIVERSE_H

65
include/openmc/tallies/filter_zernike.h Normal file
View file

@ -0,0 +1,65 @@
#ifndef OPENMC_TALLIES_FILTER_ZERNIKE_H
#define OPENMC_TALLIES_FILTER_ZERNIKE_H
#include <string>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Gives Zernike polynomial moments of a particle's position
//==============================================================================
class ZernikeFilter : public Filter
{
public:
std::string type() const override {return "zernike";}
~ZernikeFilter() = default;
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
int order() const {return order_;}
virtual void set_order(int order);
//! Cartesian x coordinate for the origin of this expansion.
double x_;
//! Cartesian y coordinate for the origin of this expansion.
double y_;
//! Maximum radius from the origin covered by this expansion.
double r_;
protected:
int order_;
};
//==============================================================================
//! Gives even order radial Zernike polynomial moments of a particle's position
//==============================================================================
class ZernikeRadialFilter : public ZernikeFilter
{
public:
std::string type() const override {return "zernikeradial";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
std::string text_label(int bin) const override;
void set_order(int order) override;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_ZERNIKE_H

49
include/openmc/tallies/tally.h Normal file
View file

@ -0,0 +1,49 @@
#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;
extern "C" double global_tally_collision;
extern "C" double global_tally_tracklength;
extern "C" double global_tally_leakage;
#pragma omp threadprivate(global_tally_absorption, global_tally_collision, \
global_tally_tracklength, global_tally_leakage)
//==============================================================================
// 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

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

@ -0,0 +1,60 @@
#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_ {0.0}; //!< elasped time in [s]
};
//==============================================================================
// Global variables
//==============================================================================
extern Timer time_active;
extern Timer time_bank;
extern Timer time_bank_sample;
extern Timer time_bank_sendrecv;
extern Timer time_finalize;
extern Timer time_inactive;
extern Timer time_initialize;
extern Timer time_tallies;
extern Timer time_total;
extern Timer time_transport;
//==============================================================================
// Non-member functions
//==============================================================================
void reset_timers();
} // namespace openmc
#endif // OPENMC_TIMER_H

1
include/openmc/xml_interface.h
View file

@ -10,7 +10,6 @@
#include "xtensor/xarray.hpp"
#include "xtensor/xadapt.hpp"
namespace openmc {
inline bool

2
openmc/arithmetic.py
View file

@ -1,6 +1,6 @@
import sys
import copy
from collections import Iterable
from collections.abc import Iterable
import numpy as np
import pandas as pd

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)

6
openmc/capi/filter.py
View file

@ -238,13 +238,15 @@ class MaterialFilter(Filter):
materials = POINTER(c_int32)()
n = c_int32()
_dll.openmc_material_filter_get_bins(self._index, materials, n)
return [Material(index=materials[i]) for i in range(n.value)]
#TODO: fix this off-by-one when materials become 0-indexed
return [Material(index=materials[i]+1) for i in range(n.value)]
@bins.setter
def bins(self, materials):
# Get material indices as int32_t[]
n = len(materials)
bins = (c_int32*n)(*(m._index for m in materials))
#TODO: fix this off-by-one when materials become 0-indexed
bins = (c_int32*n)(*(m._index-1 for m in materials))
_dll.openmc_material_filter_set_bins(self._index, n, bins)

35
openmc/capi/math.py
View file

@ -12,26 +12,26 @@ _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)]
_dll.calc_zn_c.restype = None
_dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)]
_dll.calc_zn.restype = None
_dll.calc_zn.argtypes = [c_int, c_double, c_double, ndpointer(c_double)]
_dll.calc_zn_rad_c.restype = None
_dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
_dll.calc_zn_rad.restype = None
_dll.calc_zn_rad.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):
@ -154,7 +153,7 @@ def calc_zn(n, rho, phi):
num_bins = ((n + 1) * (n + 2)) // 2
zn = np.zeros(num_bins, dtype=np.float64)
_dll.calc_zn_c(n, rho, phi, zn)
_dll.calc_zn(n, rho, phi, zn)
return zn
@ -180,9 +179,9 @@ def calc_zn_rad(n, rho):
num_bins = n // 2 + 1
zn_rad = np.zeros(num_bins, dtype=np.float64)
_dll.calc_zn_rad_c(n, rho, zn_rad)
_dll.calc_zn_rad(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):

4
openmc/data/__init__.py
View file

@ -4,7 +4,9 @@ HDF5_VERSION_MINOR = 0
HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR)
# Version of WMP nuclear data format
WMP_VERSION = 'v1.0'
WMP_VERSION_MAJOR = 1
WMP_VERSION_MINOR = 1
WMP_VERSION = (WMP_VERSION_MAJOR, WMP_VERSION_MINOR)
from .data import *

58
openmc/data/multipole.py
View file

@ -4,7 +4,7 @@ from math import exp, erf, pi, sqrt
import h5py
import numpy as np
from . import WMP_VERSION
from . import WMP_VERSION, WMP_VERSION_MAJOR
from .data import K_BOLTZMANN
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
@ -133,6 +133,8 @@ class WindowedMultipole(EqualityMixin):
Parameters
----------
name : str
Name of the nuclide using the GND naming convention
Attributes
----------
@ -171,7 +173,8 @@ class WindowedMultipole(EqualityMixin):
a/E + b/sqrt(E) + c + d sqrt(E) + ...
"""
def __init__(self):
def __init__(self, name):
self.name = name
self.spacing = None
self.sqrtAWR = None
self.E_min = None
@ -181,6 +184,10 @@ class WindowedMultipole(EqualityMixin):
self.broaden_poly = None
self.curvefit = None
@property
def name(self):
return self._name
@property
def fit_order(self):
return self.curvefit.shape[1] - 1
@ -221,6 +228,11 @@ class WindowedMultipole(EqualityMixin):
def curvefit(self):
return self._curvefit
@name.setter
def name(self, name):
cv.check_type('name', name, str)
self._name = name
@spacing.setter
def spacing(self, spacing):
if spacing is not None:
@ -322,17 +334,25 @@ class WindowedMultipole(EqualityMixin):
group = group_or_filename
else:
h5file = h5py.File(group_or_filename, 'r')
try:
version = h5file['version'].value.decode()
except AttributeError:
version = h5file['version'].value[0].decode()
if version != WMP_VERSION:
raise ValueError('The given WMP data uses version '
+ version + ' whereas your installation of the OpenMC '
'Python API expects version ' + WMP_VERSION)
group = h5file['nuclide']
out = cls()
# Make sure version matches
if 'version' in h5file.attrs:
major, minor = h5file.attrs['version']
if major != WMP_VERSION_MAJOR:
raise IOError(
'WMP data format uses version {}. {} whereas your '
'installation of the OpenMC Python API expects version '
'{}.x.'.format(major, minor, WMP_VERSION_MAJOR))
else:
raise IOError(
'WMP data does not indicate a version. Your installation of '
'the OpenMC Python API expects version {}.x data.'
.format(WMP_VERSION_MAJOR))
group = list(h5file.values())[0]
name = group.name[1:]
out = cls(name)
# Read scalars.
@ -478,13 +498,16 @@ class WindowedMultipole(EqualityMixin):
fun = np.vectorize(lambda x: self._evaluate(x, T))
return fun(E)
def export_to_hdf5(self, path, libver='earliest'):
def export_to_hdf5(self, path, mode='a', libver='earliest'):
"""Export windowed multipole data to an HDF5 file.
Parameters
----------
path : str
Path to write HDF5 file to
mode : {'r', r+', 'w', 'x', 'a'}
Mode that is used to open the HDF5 file. This is the second argument
to the :class:`h5py.File` constructor.
libver : {'earliest', 'latest'}
Compatibility mode for the HDF5 file. 'latest' will produce files
that are less backwards compatible but have performance benefits.
@ -492,12 +515,11 @@ class WindowedMultipole(EqualityMixin):
"""
# Open file and write version.
with h5py.File(path, 'w', libver=libver) as f:
f.create_dataset('version', (1, ), dtype='S10')
f['version'][:] = WMP_VERSION.encode('ASCII')
with h5py.File(path, mode, libver=libver) as f:
f.attrs['filetype'] = np.string_('data_wmp')
f.attrs['version'] = np.array(WMP_VERSION)
# Make a nuclide group.
g = f.create_group('nuclide')
g = f.create_group(self.name)
# Write scalars.
g.create_dataset('spacing', data=np.array(self.spacing))

5
openmc/data/photon.py
View file

@ -1,4 +1,5 @@
from collections import OrderedDict, Mapping, Callable
from collections import OrderedDict
from collections.abc import Mapping, Callable
from copy import deepcopy
from io import StringIO
from numbers import Integral, Real
@ -740,7 +741,7 @@ class IncidentPhoton(EqualityMixin):
shell_values.insert(0, None)
df = relax.transitions[shell].replace(
shell_values, range(len(shell_values)))
sub_group.create_dataset('transitions', data=df.as_matrix())
sub_group.create_dataset('transitions', data=df.values)
# Determine threshold
threshold = rx.xs.x[0]

2
openmc/data/reaction.py
View file

@ -924,7 +924,7 @@ class Reaction(EqualityMixin):
rx = cls(mt)
rx.q_value = group.attrs['Q_value']
rx.center_of_mass = bool(group.attrs['center_of_mass'])
rx.redundant = bool(group.attrs['redundant'])
rx.redundant = bool(group.attrs.get('redundant', False))
# Read cross section at each temperature
for T, Tgroup in group.items():

2
openmc/data/resonance_covariance.py
View file

@ -1,4 +1,4 @@
from collections import MutableSequence
from collections.abc import MutableSequence
import warnings
import io
import copy

3
openmc/filter.py
View file

@ -1,5 +1,6 @@
from abc import ABCMeta
from collections import Iterable, OrderedDict
from collections import OrderedDict
from collections.abc import Iterable
import copy
import hashlib
from itertools import product

2
openmc/mesh.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections.abc import Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys

40
openmc/mgxs/__init__.py
View file

@ -1,4 +1,44 @@
import numpy as np
from openmc.mgxs.groups import EnergyGroups
from openmc.mgxs.library import Library
from openmc.mgxs.mgxs import *
from openmc.mgxs.mdgxs import *
GROUP_STRUCTURES = {}
"""Dictionary of commonly used energy group structures, including "CASMO-X" (where X
is 2, 4, 8, 16, 25, 40 or 70) from the CASMO_ lattice physics code.
.. _CASMO: https://www.studsvik.com/SharepointFiles/CASMO-5%20Development%20and%20Applications.pdf
"""
GROUP_STRUCTURES['CASMO-2'] = np.array([
0., 6.25e-1, 2.e7])
GROUP_STRUCTURES['CASMO-4'] = np.array([
0., 6.25e-1, 5.53e3, 8.21e5, 2.e7])
GROUP_STRUCTURES['CASMO-8'] = np.array([
0., 5.8e-2, 1.4e-1, 2.8e-1, 6.25e-1, 4., 5.53e3, 8.21e5, 2.e7])
GROUP_STRUCTURES['CASMO-16'] = np.array([
0., 3.e-2, 5.8e-2, 1.4e-1, 2.8e-1, 3.5e-1, 6.25e-1, 8.5e-1,
9.72e-1, 1.02, 1.097, 1.15, 1.3, 4., 5.53e3, 8.21e5, 2.e7])
GROUP_STRUCTURES['CASMO-25'] = np.array([
0., 3.e-2, 5.8e-2, 1.4e-1, 2.8e-1, 3.5e-1, 6.25e-1, 9.72e-1, 1.02, 1.097,
1.15, 1.855, 4., 9.877, 1.5968e1, 1.4873e2, 5.53e3, 9.118e3, 1.11e5, 5.e5,
8.21e5, 1.353e6, 2.231e6, 3.679e6, 6.0655e6, 2.e7])
GROUP_STRUCTURES['CASMO-40'] = np.array([
0., 1.5e-2, 3.e-2, 4.2e-2, 5.8e-2, 8.e-2, 1.e-1, 1.4e-1,
1.8e-1, 2.2e-1, 2.8e-1, 3.5e-1, 6.25e-1, 8.5e-1, 9.5e-1,
9.72e-1, 1.02, 1.097, 1.15, 1.3, 1.5, 1.855, 2.1, 2.6, 3.3, 4.,
9.877, 1.5968e1, 2.77e1, 4.8052e1, 1.4873e2, 5.53e3, 9.118e3,
1.11e5, 5.e5, 8.21e5, 1.353e6, 2.231e6, 3.679e6, 6.0655e6, 2.e7])
GROUP_STRUCTURES['CASMO-70'] = np.array([
0., 5.e-3, 1.e-2, 1.5e-2, 2.e-2, 2.5e-2, 3.e-2, 3.5e-2, 4.2e-2,
5.e-2, 5.8e-2, 6.7e-2, 8.e-2, 1.e-1, 1.4e-1, 1.8e-1, 2.2e-1,
2.5e-1, 2.8e-1, 3.e-1, 3.2e-1, 3.5e-1, 4.e-1, 5.e-1, 6.25e-1,
7.8e-1, 8.5e-1, 9.1e-1, 9.5e-1, 9.72e-1, 9.96e-1, 1.02, 1.045,
1.071, 1.097, 1.123, 1.15, 1.3, 1.5, 1.855, 2.1, 2.6, 3.3, 4.,
9.877, 1.5968e1, 2.77e1, 4.8052e1, 7.5501e1, 1.4873e2,
3.6726e2, 9.069e2, 1.4251e3, 2.2395e3, 3.5191e3, 5.53e3,
9.118e3, 1.503e4, 2.478e4, 4.085e4, 6.734e4, 1.11e5, 1.83e5,
3.025e5, 5.e5, 8.21e5, 1.353e6, 2.231e6, 3.679e6, 6.0655e6, 2.e7])

2
openmc/mgxs/groups.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections.abc import Iterable
from numbers import Real
import copy
import sys

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':

3
openmc/model/model.py
View file

@ -181,7 +181,8 @@ class Model(object):
"""Export model to XML files."""
self.settings.export_to_xml()
self.geometry.export_to_xml()
if not self.settings.dagmc:
self.geometry.export_to_xml()
# If a materials collection was specified, export it. Otherwise, look
# for all materials in the geometry and use that to automatically build

2
openmc/polynomial.py
View file

@ -1,6 +1,5 @@
import numpy as np
import openmc
import openmc.capi as capi
from collections.abc import Iterable
@ -74,6 +73,7 @@ class ZernikeRadial(Polynomial):
return self._order
def __call__(self, r):
import openmc.capi as capi
if isinstance(r, Iterable):
return [np.sum(self._norm_coef * capi.calc_zn_rad(self.order, r_i))
for r_i in r]

16
openmc/settings.py
View file

@ -226,8 +226,8 @@ class Settings(object):
self._create_fission_neutrons = None
self._log_grid_bins = None
self._dagmc = None
self._dagmc = False
@property
def run_mode(self):
return self._run_mode
@ -375,7 +375,7 @@ class Settings(object):
@property
def dagmc(self):
return self._dagmc
@run_mode.setter
def run_mode(self, run_mode):
cv.check_value('run mode', run_mode, _RUN_MODES)
@ -523,7 +523,7 @@ class Settings(object):
def dagmc(self, dagmc):
cv.check_type('dagmc geometry', dagmc, bool)
self._dagmc = dagmc
@ptables.setter
def ptables(self, ptables):
cv.check_type('probability tables', ptables, bool)
@ -959,10 +959,10 @@ class Settings(object):
elem.text = str(self._log_grid_bins)
def _create_dagmc_subelement(self, root):
if self._dagmc is not None:
if self._dagmc:
elem = ET.SubElement(root, "dagmc")
element.text = str(self._dagmc).lower()
elem.text = str(self._dagmc).lower()
def export_to_xml(self, path='settings.xml'):
"""Export simulation settings to an XML file.
@ -1011,7 +1011,7 @@ class Settings(object):
self._create_create_fission_neutrons_subelement(root_element)
self._create_log_grid_bins_subelement(root_element)
self._create_dagmc_subelement(root_element)
# Clean the indentation in the file to be user-readable
clean_indentation(root_element)

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

3
openmc/universe.py
View file

@ -1,4 +1,5 @@
from collections import OrderedDict, Iterable
from collections import OrderedDict
from collections.abc import Iterable
from copy import copy, deepcopy
from numbers import Integral, Real
import random

6
scripts/openmc-get-multipole-data
View file

@ -29,9 +29,9 @@ parser.add_argument('-b', '--batch', action='store_true',
args = parser.parse_args()
baseUrl = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/'
files = ['WMP_Library_v1.0.tar.gz']
checksums = ['22cb675734cfccb278dffd40dcfbf26a']
baseUrl = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/'
files = ['WMP_Library_v1.1.tar.gz']
checksums = ['8523895928dd6ba63fba803e3a45d4f3']
block_size = 16384
# ==============================================================================

6
scripts/openmc-plot-mesh-tally
View file

@ -11,8 +11,10 @@ import tkinter.font as font
import tkinter.messagebox as messagebox
import tkinter.ttk as ttk
import matplotlib
matplotlib.use("TkAgg")
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
from matplotlib.backends.backend_tkagg import NavigationToolbar2TkAgg
from matplotlib.backends.backend_tkagg import NavigationToolbar2Tk
from matplotlib.figure import Figure
import matplotlib.pyplot as plt
import numpy as np
@ -56,7 +58,7 @@ class MeshPlotter(tk.Frame):
self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=1)
# Create the navigation toolbar, tied to the canvas
self.mpl_toolbar = NavigationToolbar2TkAgg(self.canvas, figureFrame)
self.mpl_toolbar = NavigationToolbar2Tk(self.canvas, figureFrame)
self.mpl_toolbar.update()
self.canvas._tkcanvas.pack(side=tk.TOP, fill=tk.BOTH, expand=1)

85
scripts/openmc-voxel-to-silovtk
View file

@ -1,85 +0,0 @@
#!/usr/bin/env python3
import struct
import sys
from argparse import ArgumentParser
import numpy as np
import h5py
def main():
# Process command line arguments
parser = ArgumentParser()
parser.add_argument('voxel_file', help='Path to voxel file')
parser.add_argument('-o', '--output', action='store',
default='plot', help='Path to output SILO or VTK file.')
parser.add_argument('-s', '--silo', action='store_true',
default=False, help='Flag to convert to SILO instead of VTK.')
args = parser.parse_args()
# Read data from voxel file
fh = h5py.File(args.voxel_file, 'r')
dimension = fh.attrs['num_voxels']
width = fh.attrs['voxel_width']
lower_left = fh.attrs['lower_left']
voxel_data = fh['data'].value
nx, ny, nz = dimension
upper_right = lower_left + width*dimension
if not args.silo:
import vtk
grid = vtk.vtkImageData()
grid.SetDimensions(nx+1, ny+1, nz+1)
grid.SetOrigin(*lower_left)
grid.SetSpacing(*width)
data = vtk.vtkDoubleArray()
data.SetName("id")
data.SetNumberOfTuples(nx*ny*nz)
for x in range(nx):
sys.stdout.write(" {}%\r".format(int(x/nx*100)))
sys.stdout.flush()
for y in range(ny):
for z in range(nz):
i = z*nx*ny + y*nx + x
data.SetValue(i, voxel_data[x, y, z])
grid.GetCellData().AddArray(data)
writer = vtk.vtkXMLImageDataWriter()
if vtk.vtkVersion.GetVTKMajorVersion() > 5:
writer.SetInputData(grid)
else:
writer.SetInput(grid)
if not args.output.endswith(".vti"):
args.output += ".vti"
writer.SetFileName(args.output)
writer.Write()
else:
import silomesh
if not args.output.endswith(".silo"):
args.output += ".silo"
silomesh.init_silo(args.output)
meshparams = list(map(int, dimension)) + list(map(float, lower_left)) + \
list(map(float, upper_right))
silomesh.init_mesh('plot', *meshparams)
silomesh.init_var("id")
for x in range(nx):
sys.stdout.write(" {}%\r".format(int(x/nx*100)))
sys.stdout.flush()
for y in range(ny):
for z in range(nz):
silomesh.set_value(float(voxel_data[x, y, z]),
x + 1, y + 1, z + 1)
print()
silomesh.finalize_var()
silomesh.finalize_mesh()
silomesh.finalize_silo()
if __name__ == '__main__':
main()

58
scripts/openmc-voxel-to-vtk Executable file
View file

@ -0,0 +1,58 @@
#!/usr/bin/env python3
import struct
import sys
from argparse import ArgumentParser
import numpy as np
import h5py
import vtk
def main():
# Process command line arguments
parser = ArgumentParser()
parser.add_argument('voxel_file', help='Path to voxel file')
parser.add_argument('-o', '--output', action='store',
default='plot', help='Path to output VTK file.')
args = parser.parse_args()
# Read data from voxel file
fh = h5py.File(args.voxel_file, 'r')
dimension = fh.attrs['num_voxels']
width = fh.attrs['voxel_width']
lower_left = fh.attrs['lower_left']
voxel_data = fh['data'].value
nx, ny, nz = dimension
upper_right = lower_left + width*dimension
grid = vtk.vtkImageData()
grid.SetDimensions(nx+1, ny+1, nz+1)
grid.SetOrigin(*lower_left)
grid.SetSpacing(*width)
data = vtk.vtkDoubleArray()
data.SetName("id")
data.SetNumberOfTuples(nx*ny*nz)
for x in range(nx):
sys.stdout.write(" {}%\r".format(int(x/nx*100)))
sys.stdout.flush()
for y in range(ny):
for z in range(nz):
i = z*nx*ny + y*nx + x
data.SetValue(i, voxel_data[x, y, z])
grid.GetCellData().AddArray(data)
writer = vtk.vtkXMLImageDataWriter()
if vtk.vtkVersion.GetVTKMajorVersion() > 5:
writer.SetInputData(grid)
else:
writer.SetInput(grid)
if not args.output.endswith(".vti"):
args.output += ".vti"
writer.SetFileName(args.output)
writer.Write()
if __name__ == '__main__':
main()

2
setup.py
View file

@ -64,7 +64,7 @@ kwargs = {
# Optional dependencies
'extras_require': {
'test': ['pytest', 'pytest-cov'],
'vtk': ['vtk', 'silomesh'],
'vtk': ['vtk'],
},
}

252
src/api.F90
View file

@ -2,191 +2,20 @@ module openmc_api
use, intrinsic :: ISO_C_BINDING
use bank_header, only: openmc_source_bank
use constants, only: K_BOLTZMANN
use eigenvalue, only: k_sum, openmc_get_keff
use constants
use error
use geometry, only: find_cell
use geometry_header
use hdf5_interface
use material_header
use math
use message_passing
use nuclide_header
use initialize, only: openmc_init_f
use particle_header
use plot, only: openmc_plot_geometry
use random_lcg, only: openmc_get_seed, openmc_set_seed
use settings
use simulation_header
use state_point, only: openmc_statepoint_write
use tally_header
use tally_filter_header
use tally_filter
use tally, only: openmc_tally_allocate
use simulation
use string, only: to_f_string
use timer_header
use volume_calc, only: openmc_calculate_volumes
use string, only: to_str
#ifdef DAGMC
use dagmc_header, only: free_memory_dagmc
#endif
implicit none
private
public :: openmc_calculate_volumes
public :: openmc_cell_filter_get_bins
public :: openmc_cell_get_id
public :: openmc_cell_set_id
public :: openmc_energy_filter_get_bins
public :: openmc_energy_filter_set_bins
public :: openmc_extend_filters
public :: openmc_extend_cells
public :: openmc_extend_materials
public :: openmc_extend_tallies
public :: openmc_filter_get_id
public :: openmc_filter_get_type
public :: openmc_filter_set_id
public :: openmc_filter_set_type
public :: openmc_finalize
public :: openmc_find_cell
public :: openmc_get_cell_index
public :: openmc_get_keff
public :: openmc_get_filter_index
public :: openmc_get_filter_next_id
public :: openmc_get_material_index
public :: openmc_get_nuclide_index
public :: openmc_get_seed
public :: openmc_get_tally_index
public :: openmc_get_tally_next_id
public :: openmc_global_tallies
public :: openmc_hard_reset
public :: openmc_init_f
public :: openmc_load_nuclide
public :: openmc_material_add_nuclide
public :: openmc_material_get_id
public :: openmc_material_get_densities
public :: openmc_material_set_density
public :: openmc_material_set_densities
public :: openmc_material_set_id
public :: openmc_material_filter_get_bins
public :: openmc_material_filter_set_bins
public :: openmc_mesh_filter_set_mesh
public :: openmc_meshsurface_filter_set_mesh
public :: openmc_next_batch
public :: openmc_nuclide_name
public :: openmc_plot_geometry
public :: openmc_reset
public :: openmc_set_seed
public :: openmc_simulation_finalize
public :: openmc_simulation_init
public :: openmc_source_bank
public :: openmc_tally_allocate
public :: openmc_tally_get_estimator
public :: openmc_tally_get_id
public :: openmc_tally_get_filters
public :: openmc_tally_get_n_realizations
public :: openmc_tally_get_nuclides
public :: openmc_tally_get_scores
public :: openmc_tally_get_type
public :: openmc_tally_results
public :: openmc_tally_set_estimator
public :: openmc_tally_set_filters
public :: openmc_tally_set_id
public :: openmc_tally_set_nuclides
public :: openmc_tally_set_scores
public :: openmc_tally_set_type
contains
!===============================================================================
! OPENMC_FINALIZE frees up memory by deallocating arrays and resetting global
! variables
!===============================================================================
function openmc_finalize() result(err) bind(C)
integer(C_INT) :: err
interface
subroutine openmc_free_bank() bind(C)
end subroutine openmc_free_bank
end interface
! Clear results
err = openmc_reset()
! Reset global variables
assume_separate = .false.
check_overlaps = .false.
confidence_intervals = .false.
create_fission_neutrons = .true.
electron_treatment = ELECTRON_LED
energy_cutoff(:) = [ZERO, 1000.0_8, ZERO, ZERO]
energy_max(:) = [INFINITY, INFINITY]
energy_min(:) = [ZERO, ZERO]
entropy_on = .false.
gen_per_batch = 1
index_entropy_mesh = -1
index_ufs_mesh = -1
keff = ONE
legendre_to_tabular = .true.
legendre_to_tabular_points = C_NONE
n_batch_interval = 1
n_lost_particles = 0
n_particles = -1
n_source_points = 0
n_state_points = 0
n_tallies = 0
output_summary = .true.
output_tallies = .true.
particle_restart_run = .false.
photon_transport = .false.
pred_batches = .false.
reduce_tallies = .true.
res_scat_on = .false.
res_scat_method = RES_SCAT_ARES
res_scat_energy_min = 0.01_8
res_scat_energy_max = 1000.0_8
restart_run = .false.
root_universe = -1
run_CE = .true.
run_mode = -1
dagmc = .false.
satisfy_triggers = .false.
call openmc_set_seed(DEFAULT_SEED)
source_latest = .false.
source_separate = .false.
source_write = .true.
survival_biasing = .false.
temperature_default = 293.6_8
temperature_method = TEMPERATURE_NEAREST
temperature_multipole = .false.
temperature_range = [ZERO, ZERO]
temperature_tolerance = 10.0_8
total_gen = 0
trigger_on = .false.
ufs = .false.
urr_ptables_on = .true.
verbosity = 7
weight_cutoff = 0.25_8
weight_survive = ONE
write_all_tracks = .false.
write_initial_source = .false.
! Deallocate arrays
call free_memory()
err = 0
#ifdef OPENMC_MPI
! Free all MPI types
call MPI_TYPE_FREE(MPI_BANK, err)
call openmc_free_bank()
#endif
end function openmc_finalize
!===============================================================================
! OPENMC_FIND_CELL determines what cell contains a given point in space
!===============================================================================
@ -222,87 +51,25 @@ contains
end function openmc_find_cell
!===============================================================================
! OPENMC_HARD_RESET reset tallies and timers as well as the pseudorandom
! generator state
!===============================================================================
function openmc_hard_reset() result(err) bind(C)
integer(C_INT) :: err
! Reset all tallies and timers
err = openmc_reset()
! Reset total generations and keff guess
keff = ONE
total_gen = 0
! Reset the random number generator state
call openmc_set_seed(DEFAULT_SEED)
end function openmc_hard_reset
!===============================================================================
! OPENMC_RESET resets tallies and timers
!===============================================================================
function openmc_reset() result(err) bind(C)
integer(C_INT) :: err
integer :: i
if (allocated(tallies)) then
do i = 1, size(tallies)
associate (t => tallies(i) % obj)
t % n_realizations = 0
if (allocated(t % results)) then
t % results(:, :, :) = ZERO
end if
end associate
end do
end if
! Reset global tallies
n_realizations = 0
if (allocated(global_tallies)) then
global_tallies(:, :) = ZERO
end if
k_col_abs = ZERO
k_col_tra = ZERO
k_abs_tra = ZERO
k_sum(:) = ZERO
! Reset timers
call time_total % reset()
call time_total % reset()
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()
err = 0
end function openmc_reset
!===============================================================================
! FREE_MEMORY deallocates and clears all global allocatable arrays in the
! program
!===============================================================================
subroutine free_memory()
subroutine free_memory() bind(C)
use bank_header
use cmfd_header
use geometry_header
use material_header
use photon_header
use plot_header
use sab_header
use settings
use simulation_header
use surface_header
use tally_derivative_header
use tally_filter_header
use tally_header
use trigger_header
use volume_header
@ -317,7 +84,6 @@ contains
call free_memory_geometry()
call free_memory_surfaces()
call free_memory_material()
call free_memory_plot()
call free_memory_volume()
call free_memory_simulation()
call free_memory_nuclide()

35
src/cell.cpp
View file

@ -476,7 +476,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
std::vector<int32_t> mat_ids;
for (auto i_mat : material_) {
if (i_mat != MATERIAL_VOID) {
mat_ids.push_back(materials[i_mat]->id);
mat_ids.push_back(materials[i_mat]->id_);
} else {
mat_ids.push_back(MATERIAL_VOID);
}
@ -649,6 +649,19 @@ read_cells(pugi::xml_node* node)
cells.push_back(new CSGCell(cell_node));
}
// Fill the cell map.
for (int i = 0; i < cells.size(); i++) {
int32_t id = cells[i]->id_;
auto search = cell_map.find(id);
if (search == cell_map.end()) {
cell_map[id] = i;
} else {
std::stringstream err_msg;
err_msg << "Two or more cells use the same unique ID: " << id;
fatal_error(err_msg);
}
}
// Populate the Universe vector and map.
for (int i = 0; i < cells.size(); i++) {
int32_t uid = cells[i]->universe_;
@ -766,7 +779,19 @@ extern "C" {
int32_t cell_id(Cell* c) {return c->id_;}
void cell_set_id(Cell* c, int32_t id) {c->id_ = id;}
void
cell_set_id(Cell* c, int32_t id)
{
c->id_ = id;
// Find the index of this cell and update the cell map.
for (int i = 0; i < cells.size(); i++) {
if (cells[i] == c) {
cell_map[id] = i;
break;
}
}
}
int cell_type(Cell* c) {return c->type_;}
@ -789,10 +814,6 @@ extern "C" {
int32_t cell_fill(Cell* c) {return c->fill_;}
int32_t cell_n_instances(Cell* c) {return c->n_instances_;}
int cell_distribcell_index(Cell* c) {return c->distribcell_index_;}
int cell_material_size(Cell* c) {return c->material_.size();}
//TODO: off-by-one
@ -807,8 +828,6 @@ extern "C" {
double cell_sqrtkT(Cell* c, int i) {return c->sqrtkT_[i];}
int32_t cell_offset(Cell* c, int map) {return c->offset_[map];}
void extend_cells_c(int32_t n)
{
cells.reserve(cells.size() + n);

52
src/cmfd_data.F90
View file

@ -99,7 +99,7 @@ contains
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m = meshes(filt % mesh)
m = meshes(filt % mesh())
end select
! Set mesh widths
@ -154,27 +154,27 @@ contains
! Reset all bins to 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins % clear()
call filter_matches(t % filter(l)) % bins % push_back(1)
call filter_matches(t % filter(l)) % bins_clear()
call filter_matches(t % filter(l)) % bins_push_back(1)
end do
! Set ijk as mesh indices
ijk = (/ i, j, k /)
! Get bin number for mesh indices
filter_matches(i_filter_mesh) % bins % data(1) = &
m % get_bin_from_indices(ijk)
call filter_matches(i_filter_mesh) &
% bins_set_data(1, m % get_bin_from_indices(ijk))
! Apply energy in filter
if (energy_filters) then
filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
end if
! Calculate score index from bins
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
% bins_data(1) - 1) * t % stride(l)
end do
! Get flux
@ -198,30 +198,30 @@ contains
! Reset all bins to 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins % clear()
call filter_matches(t % filter(l)) % bins % push_back(1)
call filter_matches(t % filter(l)) % bins_clear()
call filter_matches(t % filter(l)) % bins_push_back(1)
end do
! Set ijk as mesh indices
ijk = (/ i, j, k /)
! Get bin number for mesh indices
filter_matches(i_filter_mesh) % bins % data(1) = &
m % get_bin_from_indices(ijk)
call filter_matches(i_filter_mesh) &
% bins_set_data(1, m % get_bin_from_indices(ijk))
if (energy_filters) then
! Apply energy in filter
filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
! Set energy out bin
filter_matches(i_filter_eout) % bins % data(1) = ng - g + 1
call filter_matches(i_filter_eout) % bins_set_data(1, ng - g + 1)
end if
! Calculate score index from bins
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
% bins_data(1) - 1) * t % stride(l)
end do
! Get scattering
@ -244,23 +244,23 @@ contains
! Initialize and filter for energy
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins % clear()
call filter_matches(t % filter(l)) % bins % push_back(1)
call filter_matches(t % filter(l)) % bins_clear()
call filter_matches(t % filter(l)) % bins_push_back(1)
end do
! Set the bin for this mesh cell
i_mesh = m % get_bin_from_indices([ i, j, k ])
filter_matches(i_filter_mesh) % bins % data(1) = 12*(i_mesh - 1) + 1
call filter_matches(i_filter_mesh) % bins_set_data(1, 12*(i_mesh - 1) + 1)
! Set the energy bin if needed
if (energy_filters) then
filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
end if
score_index = 0
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
% bins_data(1) - 1) * t % stride(l)
end do
! Left surface
@ -303,30 +303,30 @@ contains
! Reset all bins to 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins % clear()
call filter_matches(t % filter(l)) % bins % push_back(1)
call filter_matches(t % filter(l)) % bins_clear()
call filter_matches(t % filter(l)) % bins_push_back(1)
end do
! Set ijk as mesh indices
ijk = (/ i, j, k /)
! Get bin number for mesh indices
filter_matches(i_filter_mesh) % bins % data(1) = &
m % get_bin_from_indices(ijk)
call filter_matches(i_filter_mesh) &
% bins_set_data(1, m % get_bin_from_indices(ijk))
! Apply energy in filter
if (energy_filters) then
filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
end if
! Apply Legendre filter
filter_matches(i_filter_legendre) % bins % data(1) = 2
call filter_matches(i_filter_legendre) % bins_set_data(1, 2)
! Calculate score index from bins
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
% bins_data(1) - 1) * t % stride(l)
end do
! Get p1 scatter rr and convert to p1 scatter xs

27
src/cmfd_execute.F90
View file

@ -13,13 +13,23 @@ 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
!==============================================================================
! EXECUTE_CMFD runs the CMFD calculation
!==============================================================================
subroutine execute_cmfd()
subroutine execute_cmfd() bind(C)
use cmfd_data, only: set_up_cmfd
use cmfd_solver, only: cmfd_solver_execute
@ -63,7 +73,7 @@ contains
! CMFD_INIT_BATCH handles cmfd options at the start of every batch
!==============================================================================
subroutine cmfd_init_batch()
subroutine cmfd_init_batch() bind(C)
! Check to activate CMFD diffusion and possible feedback
! this guarantees that when cmfd begins at least one batch of tallies are
@ -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
@ -372,7 +375,7 @@ contains
! CMFD_TALLY_INIT
!===============================================================================
subroutine cmfd_tally_init()
subroutine cmfd_tally_init() bind(C)
integer :: i
if (cmfd_run) then
do i = 1, size(cmfd_tallies)

17
src/cmfd_execute.cpp
View file

@ -7,10 +7,12 @@
#include "openmc/capi.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/simulation.h"
#include "openmc/settings.h"
namespace openmc {
extern "C" int index_cmfd_mesh;
extern "C" void
cmfd_populate_sourcecounts(int n_energy, const double* energies,
@ -22,11 +24,20 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies,
openmc_source_bank(&source_bank, &n);
// 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);
auto& m = meshes.at(settings::index_cmfd_mesh);
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

6
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(:,:)
@ -143,7 +143,7 @@ module cmfd_header
logical, public :: cmfd_run_adjoint = .false.
! CMFD run logicals
logical, public :: cmfd_on = .false.
logical(C_BOOL), public, bind(C) :: cmfd_on = .false.
! CMFD display info
character(len=25), public :: cmfd_display = 'balance'

4
src/constants.F90
View file

@ -17,6 +17,9 @@ module constants
! HDF5 data format
integer, parameter :: HDF5_VERSION(2) = [1, 0]
! WMP data format
integer, parameter :: WMP_VERSION(2) = [1, 1]
! Version numbers for binary files
integer, parameter :: VERSION_STATEPOINT(2) = [17, 0]
integer, parameter :: VERSION_PARTICLE_RESTART(2) = [2, 0]
@ -25,7 +28,6 @@ module constants
integer, parameter :: VERSION_VOLUME(2) = [1, 0]
integer, parameter :: VERSION_VOXEL(2) = [1, 0]
integer, parameter :: VERSION_MGXS_LIBRARY(2) = [1, 0]
character(10), parameter :: VERSION_MULTIPOLE = "v1.0"
! ============================================================================
! ADJUSTABLE PARAMETERS

2
src/dagmc.cpp
View file

@ -54,7 +54,7 @@ void load_dagmc_geometry()
c->fill_ = C_NONE; // no fill, single universe
cells.push_back(c);
cell_map[c->id_] = c->id_;
cell_map[c->id_] = i;
// Populate the Universe vector and dict
auto it = universe_map.find(dagmc_univ_id);

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;

536
src/eigenvalue.F90
View file

@ -2,539 +2,21 @@ module eigenvalue
use, intrinsic :: ISO_C_BINDING
use algorithm, only: binary_search
use constants, only: ZERO
use error, only: fatal_error, warning
use math, only: t_percentile
use message_passing
use random_lcg, only: prn, set_particle_seed, advance_prn_seed
use settings
use bank_header
use simulation_header
use string, only: to_str
use tally_header
use timer_header
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
function openmc_get_keff(k_combined) result(err) bind(C)
import C_INT, C_DOUBLE
real(C_DOUBLE), intent(out) :: k_combined(2)
integer(C_INT) :: err
end function
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
! processors
!===============================================================================
subroutine calculate_average_keff()
integer :: i ! overall generation within simulation
integer :: n ! number of active generations
real(8) :: alpha ! significance level for CI
real(8) :: t_value ! t-value for confidence intervals
! Determine overall generation and number of active generations
i = overall_generation()
if (current_batch > n_inactive) then
n = gen_per_batch*n_realizations + current_gen
else
n = 0
end if
if (n <= 0) then
! For inactive generations, use current generation k as estimate for next
! generation
keff = k_generation % data(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
! Determine mean
keff = k_sum(1) / n
if (n > 1) then
if (confidence_intervals) then
! Calculate t-value for confidence intervals
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, n - 1)
else
t_value = ONE
end if
! Standard deviation of the sample mean of k
keff_std = t_value * sqrt((k_sum(2)/n - keff**2) / (n - 1))
end if
end if
end subroutine calculate_average_keff
!===============================================================================
! OPENMC_GET_KEFF calculates a minimum variance estimate of k-effective based on
! a linear combination of the collision, absorption, and tracklength
! estimates. The theory behind this can be found in M. Halperin, "Almost
! linearly-optimum combination of unbiased estimates," J. Am. Stat. Assoc., 56,
! 36-43 (1961), doi:10.1080/01621459.1961.10482088. The implementation here
! follows that described in T. Urbatsch et al., "Estimation and interpretation
! of keff confidence intervals in MCNP," Nucl. Technol., 111, 169-182 (1995).
!===============================================================================
function openmc_get_keff(k_combined) result(err) bind(C)
real(C_DOUBLE), intent(out) :: k_combined(2)
integer(C_INT) :: err
integer :: l ! loop index
integer :: i, j, k ! indices referring to collision, absorption, or track
real(8) :: n ! number of realizations
real(8) :: kv(3) ! vector of k-effective estimates
real(8) :: cov(3,3) ! sample covariance matrix
real(8) :: f ! weighting factor
real(8) :: g ! sum of weighting factors
real(8) :: S(3) ! sums used for variance calculation
k_combined = ZERO
! Make sure we have at least four realizations. Notice that at the end,
! there is a N-3 term in a denominator.
if (n_realizations <= 3) then
err = -1
return
end if
! Initialize variables
n = real(n_realizations, 8)
! Copy estimates of k-effective and its variance (not variance of the mean)
kv(1) = global_tallies(RESULT_SUM, K_COLLISION) / n
kv(2) = global_tallies(RESULT_SUM, K_ABSORPTION) / n
kv(3) = global_tallies(RESULT_SUM, K_TRACKLENGTH) / n
cov(1, 1) = (global_tallies(RESULT_SUM_SQ, K_COLLISION) - &
n * kv(1) * kv(1)) / (n - ONE)
cov(2, 2) = (global_tallies(RESULT_SUM_SQ, K_ABSORPTION) - &
n * kv(2) * kv(2)) / (n - ONE)
cov(3, 3) = (global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) - &
n * kv(3) * kv(3)) / (n - ONE)
! Calculate covariances based on sums with Bessel's correction
cov(1, 2) = (k_col_abs - n * kv(1) * kv(2)) / (n - ONE)
cov(1, 3) = (k_col_tra - n * kv(1) * kv(3)) / (n - ONE)
cov(2, 3) = (k_abs_tra - n * kv(2) * kv(3)) / (n - ONE)
cov(2, 1) = cov(1, 2)
cov(3, 1) = cov(1, 3)
cov(3, 2) = cov(2, 3)
! Check to see if two estimators are the same; this is guaranteed to happen
! in MG-mode with survival biasing when the collision and absorption
! estimators are the same, but can theoretically happen at anytime.
! If it does, the standard estimators will produce floating-point
! exceptions and an expression specifically derived for the combination of
! two estimators (vice three) should be used instead.
! First we will identify if there are any matching estimators
if ((abs(kv(1) - kv(2)) / kv(1) < FP_REL_PRECISION) .and. &
(abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) then
! 1 and 2 match, so only use 1 and 3 in our comparisons
i = 1
j = 3
else if ((abs(kv(1) - kv(3)) / kv(1) < FP_REL_PRECISION) .and. &
(abs(cov(1, 1) - cov(3, 3)) / cov(1, 1) < FP_REL_PRECISION)) then
! 1 and 3 match, so only use 1 and 2 in our comparisons
i = 1
j = 2
else if ((abs(kv(2) - kv(3)) / kv(2) < FP_REL_PRECISION) .and. &
(abs(cov(2, 2) - cov(3, 3)) / cov(2, 2) < FP_REL_PRECISION)) then
! 2 and 3 match, so only use 1 and 2 in our comparisons
i = 1
j = 2
else
! No two estimators match, so set i to 0 and this will be the indicator
! to use all three estimators.
i = 0
end if
if (i == 0) then
! Use three estimators as derived in the paper by Urbatsch
! Initialize variables
g = ZERO
S = ZERO
do l = 1, 3
! Permutations of estimates
if (l == 1) then
! i = collision, j = absorption, k = tracklength
i = 1
j = 2
k = 3
elseif (l == 2) then
! i = absortion, j = tracklength, k = collision
i = 2
j = 3
k = 1
elseif (l == 3) then
! i = tracklength, j = collision, k = absorption
i = 3
j = 1
k = 2
end if
! Calculate weighting
f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) + &
cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k))
! Add to S sums for variance of combined estimate
S(1) = S(1) + f * cov(1, l)
S(2) = S(2) + (cov(j, j) + cov(k, k) - TWO * cov(j, k)) * kv(l) * kv(l)
S(3) = S(3) + (cov(k, k) + cov(i, j) - cov(j, k) - &
cov(i, k)) * kv(l) * kv(j)
! Add to sum for combined k-effective
k_combined(1) = k_combined(1) + f * kv(l)
g = g + f
end do
! Complete calculations of S sums
S = (n - ONE) * S
S(1) = (n - ONE)**2 * S(1)
! Calculate combined estimate of k-effective
k_combined(1) = k_combined(1) / g
! Calculate standard deviation of combined estimate
g = (n - ONE)**2 * g
k_combined(2) = sqrt(S(1) / &
(g * n * (n - THREE)) * (ONE + n * ((S(2) - TWO * S(3)) / g)))
else
! Use only two estimators
! These equations are derived analogously to that done in the paper by
! Urbatsch, but are simpler than for the three estimators case since the
! block matrices of the three estimator equations reduces to scalars here
! Store the commonly used term
f = kv(i) - kv(j)
g = cov(i, i) + cov(j, j) - TWO * cov(i, j)
! Calculate combined estimate of k-effective
k_combined(1) = kv(i) - (cov(i, i) - cov(i, j)) / g * f
! Calculate standard deviation of combined estimate
k_combined(2) = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) * &
(g + n * f * f) / (n * (n - TWO) * g * g)
k_combined(2) = sqrt(k_combined(2))
end if
err = 0
end function openmc_get_keff
#ifdef _OPENMP
!===============================================================================
! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission
@ -542,7 +24,7 @@ contains
! to preserve the order of the bank when using varying numbers of threads.
!===============================================================================
subroutine join_bank_from_threads()
subroutine join_bank_from_threads() bind(C)
integer(8) :: total ! total number of fission bank sites
integer :: i ! loop index for threads

496
src/eigenvalue.cpp
View file

@ -1,16 +1,28 @@
#include "openmc/eigenvalue.h"
#include "xtensor/xbuilder.hpp"
#include "xtensor/xmath.hpp"
#include "xtensor/xtensor.hpp"
#include "xtensor/xview.hpp"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.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 <array>
#include <cmath> // for sqrt, abs, pow
#include <string>
namespace openmc {
@ -18,6 +30,8 @@ namespace openmc {
// Global variables
//==============================================================================
double keff_generation;
std::array<double, 2> k_sum;
std::vector<double> entropy;
xt::xtensor<double, 1> source_frac;
@ -25,6 +39,457 @@ 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;
std::vector<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.data(), temp_sites.data() + settings::n_particles, source_bank);
#endif
time_bank_sendrecv.stop();
time_bank.stop();
}
void calculate_average_keff()
{
// Determine overall generation and number of active generations
int i = overall_generation() - 1;
int n;
if (simulation::current_batch > settings::n_inactive) {
n = settings::gen_per_batch*n_realizations + simulation::current_gen;
} else {
n = 0;
}
if (n <= 0) {
// For inactive generations, use current generation k as estimate for next
// generation
simulation::keff = simulation::k_generation[i];
} else {
// Sample mean of keff
k_sum[0] += simulation::k_generation[i];
k_sum[1] += std::pow(simulation::k_generation[i], 2);
// Determine mean
simulation::keff = k_sum[0] / n;
if (n > 1) {
double t_value;
if (settings::confidence_intervals) {
// Calculate t-value for confidence intervals
double alpha = 1.0 - CONFIDENCE_LEVEL;
t_value = t_percentile_c(1.0 - alpha/2.0, n - 1);
} else {
t_value = 1.0;
}
// Standard deviation of the sample mean of k
simulation::keff_std = t_value * std::sqrt((k_sum[1]/n -
std::pow(simulation::keff, 2)) / (n - 1));
}
}
}
int openmc_get_keff(double* k_combined)
{
k_combined[0] = 0.0;
k_combined[1] = 0.0;
// Make sure we have at least four realizations. Notice that at the end,
// there is a N-3 term in a denominator.
if (n_realizations <= 3) {
return -1;
}
// Initialize variables
int64_t n = n_realizations;
// Copy estimates of k-effective and its variance (not variance of the mean)
auto gt = global_tallies();
std::array<double, 3> kv {};
xt::xtensor<double, 2> cov = xt::zeros<double>({3, 3});
kv[0] = gt(K_COLLISION, RESULT_SUM) / n;
kv[1] = gt(K_ABSORPTION, RESULT_SUM) / n;
kv[2] = gt(K_TRACKLENGTH, RESULT_SUM) / n;
cov(0, 0) = (gt(K_COLLISION, RESULT_SUM_SQ) - n*kv[0]*kv[0]) / (n - 1);
cov(1, 1) = (gt(K_ABSORPTION, RESULT_SUM_SQ) - n*kv[1]*kv[1]) / (n - 1);
cov(2, 2) = (gt(K_TRACKLENGTH, RESULT_SUM_SQ) - n*kv[2]*kv[2]) / (n - 1);
// Calculate covariances based on sums with Bessel's correction
cov(0, 1) = (simulation::k_col_abs - n * kv[0] * kv[1]) / (n - 1);
cov(0, 2) = (simulation::k_col_tra - n * kv[0] * kv[2]) / (n - 1);
cov(1, 2) = (simulation::k_abs_tra - n * kv[1] * kv[2]) / (n - 1);
cov(1, 0) = cov(0, 1);
cov(2, 0) = cov(0, 2);
cov(2, 1) = cov(1, 2);
// Check to see if two estimators are the same; this is guaranteed to happen
// in MG-mode with survival biasing when the collision and absorption
// estimators are the same, but can theoretically happen at anytime.
// If it does, the standard estimators will produce floating-point
// exceptions and an expression specifically derived for the combination of
// two estimators (vice three) should be used instead.
// First we will identify if there are any matching estimators
int i, j, k;
if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) &&
(std::abs(cov(0, 0) - cov(1, 1)) / cov(0, 0) < FP_REL_PRECISION)) {
// 0 and 1 match, so only use 0 and 2 in our comparisons
i = 0;
j = 2;
} else if ((std::abs(kv[0] - kv[2]) / kv[0] < FP_REL_PRECISION) &&
(std::abs(cov(0, 0) - cov(2, 2)) / cov(0, 0) < FP_REL_PRECISION)) {
// 0 and 2 match, so only use 0 and 1 in our comparisons
i = 0;
j = 1;
} else if ((std::abs(kv[1] - kv[2]) / kv[1] < FP_REL_PRECISION) &&
(std::abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) {
// 1 and 2 match, so only use 0 and 1 in our comparisons
i = 0;
j = 1;
} else {
// No two estimators match, so set i to -1 and this will be the indicator
// to use all three estimators.
i = -1;
}
if (i == -1) {
// Use three estimators as derived in the paper by Urbatsch
// Initialize variables
double g = 0.0;
std::array<double, 3> S {};
for (int l = 0; l < 3; ++l) {
// Permutations of estimates
switch (l) {
case 0:
// i = collision, j = absorption, k = tracklength
i = 0;
j = 1;
k = 2;
break;
case 1:
// i = absortion, j = tracklength, k = collision
i = 1;
j = 2;
k = 0;
break;
case 2:
// i = tracklength, j = collision, k = absorption
i = 2;
j = 0;
k = 1;
break;
}
// Calculate weighting
double f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) +
cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k));
// Add to S sums for variance of combined estimate
S[0] += f * cov(0, l);
S[1] += (cov(j, j) + cov(k, k) - 2.0 * cov(j, k)) * kv[l] * kv[l];
S[2] += (cov(k, k) + cov(i, j) - cov(j, k) - cov(i, k)) * kv[l] * kv[j];
// Add to sum for combined k-effective
k_combined[0] += f * kv[l];
g += f;
}
// Complete calculations of S sums
for (auto& S_i : S) {
S_i *= (n - 1);
}
S[0] *= (n - 1)*(n - 1);
// Calculate combined estimate of k-effective
k_combined[0] /= g;
// Calculate standard deviation of combined estimate
g *= (n - 1)*(n - 1);
k_combined[1] = std::sqrt(S[0] / (g*n*(n - 3)) *
(1 + n*((S[1] - 2*S[2]) / g)));
} else {
// Use only two estimators
// These equations are derived analogously to that done in the paper by
// Urbatsch, but are simpler than for the three estimators case since the
// block matrices of the three estimator equations reduces to scalars here
// Store the commonly used term
double f = kv[i] - kv[j];
double g = cov(i, i) + cov(j, j) - 2.0*cov(i, j);
// Calculate combined estimate of k-effective
k_combined[0] = kv[i] - (cov(i, i) - cov(i, j)) / g * f;
// Calculate standard deviation of combined estimate
k_combined[1] = (cov(i, i)*cov(j, j) - cov(i, j)*cov(i, j)) *
(g + n*f*f) / (n*(n - 2)*g*g);
k_combined[1] = std::sqrt(k_combined[1]);
}
return 0;
}
void shannon_entropy()
{
// Get pointer to entropy mesh
@ -66,7 +531,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 +547,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 +567,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,20 +592,40 @@ 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);
}
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" double entropy_c(int i)
{
return entropy.at(i - 1);
@ -151,5 +636,6 @@ extern "C" void entropy_clear()
entropy.clear();
}
extern "C" void k_sum_reset() { k_sum.fill(0.0); }
} // namespace openmc

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

131
src/finalize.cpp
View file

@ -1,10 +1,137 @@
#include "openmc/finalize.h"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/geometry.h"
#include "openmc/message_passing.h"
#include "openmc/nuclide.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/timer.h"
#include "openmc/tallies/tally.h"
void openmc_free_bank()
using namespace openmc;
// Functions defined in Fortran
extern "C" void free_memory();
extern "C" void reset_timers_f();
int openmc_finalize()
{
// Clear results
openmc_reset();
// Reset global variables
settings::assume_separate = false;
settings::check_overlaps = false;
settings::confidence_intervals = false;
settings::create_fission_neutrons = true;
settings::electron_treatment = ELECTRON_LED;
settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0};
settings::entropy_on = false;
settings::gen_per_batch = 1;
settings::index_entropy_mesh = -1;
settings::index_ufs_mesh = -1;
settings::legendre_to_tabular = true;
settings::legendre_to_tabular_points = -1;
settings::n_particles = -1;
settings::output_summary = true;
settings::output_tallies = true;
settings::particle_restart_run = false;
settings::photon_transport = false;
settings::reduce_tallies = true;
settings::res_scat_on = false;
settings::res_scat_method = RES_SCAT_ARES;
settings::res_scat_energy_min = 0.01;
settings::res_scat_energy_max = 1000.0;
settings::restart_run = false;
settings::run_CE = true;
settings::run_mode = -1;
settings::dagmc = false;
settings::source_latest = false;
settings::source_separate = false;
settings::source_write = true;
settings::survival_biasing = false;
settings::temperature_default = 293.6;
settings::temperature_method = TEMPERATURE_NEAREST;
settings::temperature_multipole = false;
settings::temperature_range = {0.0, 0.0};
settings::temperature_tolerance = 10.0;
settings::trigger_on = false;
settings::trigger_predict = false;
settings::trigger_batch_interval = 1;
settings::ufs_on = false;
settings::urr_ptables_on = true;
settings::verbosity = 7;
settings::weight_cutoff = 0.25;
settings::weight_survive = 1.0;
settings::write_all_tracks = false;
settings::write_initial_source = false;
simulation::keff = 1.0;
simulation::n_lost_particles = 0;
simulation::satisfy_triggers = false;
simulation::total_gen = 0;
energy_max = {INFTY, INFTY};
energy_min = {0.0, 0.0};
n_tallies = 0;
openmc_root_universe = -1;
openmc_set_seed(DEFAULT_SEED);
// Deallocate arrays
free_memory();
// Free all MPI types
#ifdef OPENMC_MPI
MPI_Type_free(&openmc::mpi::bank);
MPI_Type_free(&mpi::bank);
#endif
return 0;
}
int openmc_reset()
{
for (int i = 1; i <= n_tallies; ++i) {
openmc_tally_reset(i);
}
// Reset global tallies (can't really use global_tallies() right now because
// it doesn't have any information about whether the underlying buffer was
// allocated)
n_realizations = 0;
double* buffer = nullptr;
openmc_global_tallies(&buffer);
if (buffer) {
for (int i = 0; i < 3*N_GLOBAL_TALLIES; ++i) {
buffer[i] = 0.0;
}
}
simulation::k_col_abs = 0.0;
simulation::k_col_tra = 0.0;
simulation::k_abs_tra = 0.0;
k_sum = {0.0, 0.0};
// Reset timers
reset_timers();
reset_timers_f();
return 0;
}
int openmc_hard_reset()
{
// Reset all tallies and timers
openmc_reset();
// Reset total generations and keff guess
simulation::keff = 1.0;
simulation::total_gen = 0;
// Reset the random number generator state
openmc_set_seed(DEFAULT_SEED);
return 0;
}

6
src/geometry.cpp
View file

@ -16,8 +16,6 @@ namespace openmc {
std::vector<int64_t> overlap_check_count;
constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE
//==============================================================================
extern "C" bool
@ -91,7 +89,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 +253,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 << ","

39
src/geometry_aux.cpp
View file

@ -12,6 +12,8 @@
#include "openmc/material.h"
#include "openmc/settings.h"
#include "openmc/surface.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/filter_distribcell.h"
namespace openmc {
@ -181,12 +183,15 @@ neighbor_lists()
//==============================================================================
void
prepare_distribcell(int32_t* filter_cell_list, int n)
prepare_distribcell()
{
// Read the list of cells contained in distribcell filters from Fortran.
// Find all cells listed in a DistribcellFilter.
std::unordered_set<int32_t> distribcells;
for (int i = 0; i < n; i++) {
distribcells.insert(filter_cell_list[i]);
for (auto& filt : tally_filters) {
auto* distrib_filt = dynamic_cast<DistribcellFilter*>(filt.get());
if (distrib_filt) {
distribcells.insert(distrib_filt->cell_);
}
}
// Find all cells with distributed materials or temperatures. Make sure that
@ -393,29 +398,11 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
}
}
//==============================================================================
int
distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset,
int32_t root_univ)
std::string
distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset)
{
Universe& root = *universes[root_univ];
std::string path_ {distribcell_path_inner(target_cell, map, target_offset,
root, 0)};
return path_.size() + 1;
}
//==============================================================================
void
distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset,
int32_t root_univ, char* path)
{
Universe& root = *universes[root_univ];
std::string path_ {distribcell_path_inner(target_cell, map, target_offset,
root, 0)};
path_.copy(path, path_.size());
path[path_.size()] = '\0';
auto& root_univ = *universes[openmc_root_universe];
return distribcell_path_inner(target_cell, map, target_offset, root_univ, 0);
}
//==============================================================================

Some files were not shown because too many files have changed in this diff Show more