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

Merge pull request #1154 from paulromano/cpp-output

Browse source 
Translate most output functions to C++
This commit is contained in:
Sterling Harper 2019-02-12 11:23:25 -05:00 • committed by GitHub
commit e86761814b
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
36 changed files with 977 additions and 1613 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -311,7 +311,6 @@ add_library(libopenmc SHARED
src/dict_header.F90
src/eigenvalue.F90
src/endf.F90
src/endf_header.F90
src/error.F90
src/geometry.F90
src/geometry_header.F90

2
include/openmc/capi.h
View file

@ -100,7 +100,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[], bool* write_source);
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);

1
include/openmc/constants.h
View file

@ -28,6 +28,7 @@ constexpr std::array<int, 3> VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELE
constexpr int HDF5_VERSION[] {2, 0};
// Version numbers for binary files
constexpr std::array<int, 2> VERSION_STATEPOINT {17, 0};
constexpr std::array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
constexpr std::array<int, 2> VERSION_TRACK {2, 0};
constexpr std::array<int, 2> VERSION_SUMMARY {6, 0};

9
include/openmc/eigenvalue.h
View file

@ -9,6 +9,7 @@
#include <vector>
#include "xtensor/xtensor.hpp"
#include <hdf5.h>
#include "openmc/particle.h"
@ -78,6 +79,14 @@ void ufs_count_sites();
//! Get UFS weight corresponding to particle's location
extern "C" double ufs_get_weight(const Particle* p);
//! Write data related to k-eigenvalue to statepoint
//! \param[in] group HDF5 group
void write_eigenvalue_hdf5(hid_t group);
//! Read data related to k-eigenvalue from statepoint
//! \param[in] group HDF5 group
void read_eigenvalue_hdf5(hid_t group);
} // namespace openmc
#endif // OPENMC_EIGENVALUE_H

35
include/openmc/error.h
View file

@ -9,11 +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,
int level);
inline void
set_errmsg(const char* message)
{
@ -32,17 +27,7 @@ set_errmsg(const std::stringstream& message)
std::strcpy(openmc_err_msg, message.str().c_str());
}
inline
void fatal_error(const char* message)
{
fatal_error_from_c(message, std::strlen(message));
}
inline
void fatal_error(const std::string& message)
{
fatal_error_from_c(message.c_str(), message.length());
}
[[noreturn]] void fatal_error(const std::string& message, int err=-1);
inline
void fatal_error(const std::stringstream& message)
@ -51,28 +36,20 @@ void fatal_error(const std::stringstream& message)
}
inline
void warning(const std::string& message)
void fatal_error(const char* message)
{
warning_from_c(message.c_str(), message.length());
fatal_error({message, std::strlen(message)});
}
void warning(const std::string& message);
inline
void warning(const std::stringstream& message)
{
warning(message.str());
}
inline
void write_message(const char* message, int level)
{
write_message_from_c(message, std::strlen(message), level);
}
inline
void write_message(const std::string& message, int level)
{
write_message_from_c(message.c_str(), message.length(), level);
}
void write_message(const std::string& message, int level);
inline
void write_message(const std::stringstream& message, int level)

13
include/openmc/hdf5_interface.h
View file

@ -229,6 +229,19 @@ void read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)
read_dataset(obj_id, name, H5TypeMap<T>::type_id, &buffer, indep);
}
// overload for std::string
inline void
read_dataset(hid_t obj_id, const char* name, std::string& str, bool indep=false)
{
// Create buffer to read data into
auto n = attribute_typesize(obj_id, name);
char buffer[n];
// Read attribute and set string
read_string(obj_id, name, n, buffer, indep);
str = std::string{buffer, n};
}
template <typename T>
void read_dataset(hid_t dset, std::vector<T>& vec, bool indep=false)
{

2
include/openmc/math_functions.h
View file

@ -34,7 +34,7 @@ extern "C" double normal_percentile(double p);
//! \return The requested percentile
//==============================================================================
extern "C" double t_percentile_c(double p, int df);
extern "C" double t_percentile(double p, int df);
//==============================================================================
//! Calculate the n-th order Legendre polynomials at the value of x.

2
include/openmc/mesh.h
View file

@ -128,7 +128,7 @@ extern "C" void read_meshes(pugi::xml_node* root);
//! Write mesh data to an HDF5 group
//! \param[in] group HDF5 group
extern "C" void meshes_to_hdf5(hid_t group);
void meshes_to_hdf5(hid_t group);
} // namespace openmc

45
include/openmc/output.h
View file

@ -10,42 +10,51 @@
namespace openmc {
//==============================================================================
//! \brief Display the main title banner as well as information about the
//! program developers, version, and date/time which the problem was run.
void title();
//! Display a header block.
//!
//
//! \param msg The main text of the header
//! \param level The lowest verbosity level at which this header is printed
//==============================================================================
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 the attributes of a particle.
//==============================================================================
extern "C" void print_particle(Particle* p);
//==============================================================================
//! Display plot information.
//==============================================================================
void print_plot();
//==============================================================================
//! Display information regarding cell overlap checking.
//==============================================================================
void print_overlap_check();
extern "C" void title();
//! Display information about command line usage of OpenMC
void print_usage();
//! Display current version and copright/license information
void print_version();
//! Display header listing what physical values will displayed
void print_columns();
//! Display information about a generation of neutrons
void print_generation();
//! \brief Display last batch's tallied value of the neutron multiplication
//! factor as well as the average value if we're in active batches
void print_batch_keff();
//! Display time elapsed for various stages of a run
void print_runtime();
//! Display results for global tallies including k-effective estimators
void print_results();
} // namespace openmc
#endif // OPENMC_OUTPUT_H

5
include/openmc/settings.h
View file

@ -25,12 +25,12 @@ namespace settings {
// Boolean flags
extern "C" bool assume_separate; //!< assume tallies are spatially separate?
extern "C" bool check_overlaps; //!< check overlaps in geometry?
extern "C" bool cmfd_run; //!< use CMFD?
extern "C" bool confidence_intervals; //!< use confidence intervals for results?
extern "C" bool create_fission_neutrons; //!< create fission neutrons (fixed source)?
extern "C" bool dagmc; //!< indicator of DAGMC geometry
extern "C" bool entropy_on; //!< calculate Shannon entropy?
extern "C" bool legendre_to_tabular; //!< convert Legendre distributions to tabular?
extern "C" bool output_summary; //!< write summary.h5?
extern bool output_summary; //!< write summary.h5?
extern "C" bool output_tallies; //!< write tallies.out?
extern "C" bool particle_restart_run; //!< particle restart run?
extern "C" bool photon_transport; //!< photon transport turned on?
@ -49,7 +49,6 @@ extern bool ufs_on; //!< uniform fission site method on?
extern bool urr_ptables_on; //!< use unresolved resonance prob. tables?
extern "C" bool write_all_tracks; //!< write track files for every particle?
extern "C" bool write_initial_source; //!< write out initial source file?
extern "C" bool dagmc; //!< indicator of DAGMC geometry
// Paths to various files
extern std::string path_cross_sections; //!< path to cross_sections.xml

9
include/openmc/state_point.h
View file

@ -9,11 +9,12 @@
namespace openmc {
void load_state_point();
void write_source_point(const char* filename);
extern "C" void write_source_bank(hid_t group_id);
extern "C" void read_source_bank(hid_t group_id);
extern "C" void write_tally_results_nr(hid_t file_id);
extern "C" void restart_set_keff();
void write_source_bank(hid_t group_id);
void read_source_bank(hid_t group_id);
void write_tally_results_nr(hid_t file_id);
void restart_set_keff();
} // namespace openmc
#endif // OPENMC_STATE_POINT_H

8
openmc/capi/math.py
View file

@ -6,8 +6,8 @@ from numpy.ctypeslib import ndpointer
from . import _dll
_dll.t_percentile_c.restype = c_double
_dll.t_percentile_c.argtypes = [c_double, c_int]
_dll.t_percentile.restype = c_double
_dll.t_percentile.argtypes = [c_double, c_int]
_dll.calc_pn_c.restype = None
_dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)]
@ -58,7 +58,7 @@ def t_percentile(p, df):
"""
return _dll.t_percentile_c(p, df)
return _dll.t_percentile(p, df)
def calc_pn(n, x):
@ -261,7 +261,7 @@ def normal_variate(mean_value, std_dev):
Parameters
----------
mean_value : float
Mean of the Normal distribution
Mean of the Normal distribution
std_dev : float
Standard deviation of the normal distribution

134
src/algorithm.F90
View file

@ -5,12 +5,6 @@ module algorithm
implicit none
integer, parameter :: MAX_ITERATION = 64
interface binary_search
module procedure binary_search_real, binary_search_int4, binary_search_int8
end interface binary_search
interface sort
module procedure sort_int, sort_real, sort_vector_int, sort_vector_real
end interface sort
@ -21,134 +15,6 @@ module algorithm
contains
!===============================================================================
! BINARY_SEARCH performs a binary search of an array to find where a specific
! value lies in the array. This is used extensively for energy grid searching
!===============================================================================
pure function binary_search_real(array, n, val) result(array_index)
integer, intent(in) :: n
real(8), intent(in) :: array(n)
real(8), intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
array_index = -1
return
end if
n_iteration = 0
do while (R - L > 1)
! Find values at midpoint
array_index = L + (R - L)/2
if (val >= array(array_index)) then
L = array_index
else
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
array_index = -2
return
end if
end do
array_index = L
end function binary_search_real
pure function binary_search_int4(array, n, val) result(array_index)
integer, intent(in) :: n
integer, intent(in) :: array(n)
integer, intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
array_index = -1
return
end if
n_iteration = 0
do while (R - L > 1)
! Find values at midpoint
array_index = L + (R - L)/2
if (val >= array(array_index)) then
L = array_index
else
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
array_index = -2
return
end if
end do
array_index = L
end function binary_search_int4
pure function binary_search_int8(array, n, val) result(array_index)
integer, intent(in) :: n
integer(8), intent(in) :: array(n)
integer(8), intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
array_index = -1
return
end if
n_iteration = 0
do while (R - L > 1)
! Find values at midpoint
array_index = L + (R - L)/2
if (val >= array(array_index)) then
L = array_index
else
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
array_index = -2
return
end if
end do
array_index = L
end function binary_search_int8
!===============================================================================
! SORT sorts an array in place using an insertion sort.
!===============================================================================

51
src/constants.F90
View file

@ -14,20 +14,9 @@ module constants
integer, parameter :: &
VERSION(3) = [VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE]
! HDF5 data format
integer, parameter :: HDF5_VERSION(2) = [2, 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]
integer, parameter :: VERSION_TRACK(2) = [2, 0]
integer, parameter :: VERSION_SUMMARY(2) = [6, 0]
integer, parameter :: VERSION_VOLUME(2) = [1, 0]
integer, parameter :: VERSION_VOXEL(2) = [1, 0]
integer, parameter :: VERSION_MGXS_LIBRARY(2) = [1, 0]
! ============================================================================
! ADJUSTABLE PARAMETERS
@ -41,21 +30,14 @@ module constants
! Used for surface current tallies
real(8), parameter :: TINY_BIT = 1e-8_8
! User for precision in geometry
real(C_DOUBLE), bind(C, name='FP_PRECISION') :: FP_PRECISION = 1e-14_8
real(C_DOUBLE), bind(C, name='FP_REL_PRECISION') :: FP_REL_PRECISION = 1e-5_8
real(C_DOUBLE), bind(C, name='FP_COINCIDENT') :: FP_COINCIDENT = 1e-12_8
! Maximum number of collisions/crossings
integer, parameter :: MAX_EVENTS = 1000000
integer, parameter :: MAX_SAMPLE = 100000
! Maximum number of secondary particles created
integer, parameter :: MAX_SECONDARY = 1000
! Maximum number of words in a single line, length of line, and length of
! single word
integer, parameter :: MAX_WORDS = 500
integer, parameter :: MAX_LINE_LEN = 250
integer, parameter :: MAX_WORD_LEN = 150
integer, parameter :: MAX_FILE_LEN = 255
@ -70,9 +52,7 @@ module constants
PI = 3.1415926535898_8, & ! pi
MASS_NEUTRON = 1.00866491588_8, & ! mass of a neutron in amu
MASS_NEUTRON_EV = 939.5654133e6_8, & ! mass of a neutron in eV/c^2
AMU = 1.660539040e-27_8, & ! 1 amu in kg
C_LIGHT = 2.99792458e8_8, & ! speed of light in m/s
N_AVOGADRO = 0.6022140857_8, & ! Avogadro's number in 10^24/mol
K_BOLTZMANN = 8.6173303e-5_8, & ! Boltzmann constant in eV/K
INFINITY = huge(0.0_8), & ! positive infinity
ZERO = 0.0_8, &
@ -112,14 +92,6 @@ module constants
! ============================================================================
! CROSS SECTION RELATED CONSTANTS
! Interpolation flag
integer, parameter :: &
HISTOGRAM = 1, & ! y is constant in x
LINEAR_LINEAR = 2, & ! y is linear in x
LINEAR_LOG = 3, & ! y is linear in ln(x)
LOG_LINEAR = 4, & ! ln(y) is linear in x
LOG_LOG = 5 ! ln(y) is linear in ln(x)
! Particle type
integer, parameter :: &
NEUTRON = 0, &
@ -127,17 +99,6 @@ module constants
ELECTRON = 2, &
POSITRON = 3
! Angular distribution type
integer, parameter :: &
ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution (CE)
ANGLE_32_EQUI = 2, & ! 32 equiprobable bins (CE)
ANGLE_TABULAR = 3, & ! Tabular angular distribution (CE or MG)
ANGLE_LEGENDRE = 4, & ! Legendre angular distribution (MG)
ANGLE_HISTOGRAM = 5 ! Histogram angular distribution (MG)
! Number of mu bins to use when converting Legendres to tabular type
integer, parameter :: DEFAULT_NMU = 33
! Reaction types
integer, parameter :: &
TOTAL_XS = 1, ELASTIC = 2, N_NONELASTIC = 3, N_LEVEL = 4, MISC = 5, &
@ -173,24 +134,12 @@ module constants
MGXS_ISOTROPIC = 1, & ! Isotropically Weighted Data
MGXS_ANGLE = 2 ! Data by Angular Bins
! Flag to denote this was a macroscopic data object
real(8), parameter :: &
MACROSCOPIC_AWR = -TWO
! Secondary particle emission type
integer, parameter :: &
EMISSION_PROMPT = 1, & ! Prompt emission of secondary particle
EMISSION_DELAYED = 2, & ! Delayed emission of secondary particle
EMISSION_TOTAL = 3 ! Yield represents total emission (prompt + delayed)
! Library types
integer, parameter :: &
LIBRARY_NEUTRON = 1, &
LIBRARY_THERMAL = 2, &
LIBRARY_PHOTON = 3, &
LIBRARY_MULTIGROUP = 4, &
LIBRARY_WMP = 5
! Maximum number of partial fission reactions
integer, parameter :: PARTIAL_FISSION_MAX = 4

6
src/eigenvalue.cpp
View file

@ -331,7 +331,7 @@ void calculate_average_keff()
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);
t_value = t_percentile(1.0 - alpha/2.0, n - 1);
} else {
t_value = 1.0;
}
@ -622,7 +622,7 @@ double ufs_get_weight(const Particle* p)
}
}
extern "C" void write_eigenvalue_hdf5(hid_t group)
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);
@ -638,7 +638,7 @@ extern "C" void write_eigenvalue_hdf5(hid_t group)
write_dataset(group, "k_combined", k_combined);
}
extern "C" void read_eigenvalue_hdf5(hid_t group)
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;

191
src/endf_header.F90
View file

@ -1,191 +0,0 @@
module endf_header
use algorithm, only: binary_search
use constants, only: ZERO, HISTOGRAM, LINEAR_LINEAR, LINEAR_LOG, &
LOG_LINEAR, LOG_LOG
use hdf5_interface
implicit none
type, abstract :: Function1D
contains
procedure(function1d_evaluate_), deferred :: evaluate
procedure(function1d_from_hdf5_), deferred :: from_hdf5
end type Function1D
abstract interface
pure function function1d_evaluate_(this, x) result(y)
import Function1D
class(Function1D), intent(in) :: this
real(8), intent(in) :: x
real(8) :: y
end function function1d_evaluate_
subroutine function1d_from_hdf5_(this, dset_id)
import Function1D, HID_T
class(Function1D), intent(inout) :: this
integer(HID_T), intent(in) :: dset_id
end subroutine function1d_from_hdf5_
end interface
!===============================================================================
! POLYNOMIAL represents a one-dimensional function expressed as a polynomial
!===============================================================================
type, extends(Function1D) :: Polynomial
real(8), allocatable :: coef(:) ! coefficients
contains
procedure :: from_hdf5 => polynomial_from_hdf5
procedure :: evaluate => polynomial_evaluate
end type Polynomial
!===============================================================================
! TABULATED1D represents a one-dimensional interpolable function
!===============================================================================
type, extends(Function1D) :: Tabulated1D
integer :: n_regions = 0 ! # of interpolation regions
integer, allocatable :: nbt(:) ! values separating interpolation regions
integer, allocatable :: int(:) ! interpolation scheme
integer :: n_pairs ! # of pairs of (x,y) values
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
contains
procedure :: from_hdf5 => tabulated1d_from_hdf5
procedure :: evaluate => tabulated1d_evaluate
end type Tabulated1D
contains
!===============================================================================
! Polynomial implementation
!===============================================================================
subroutine polynomial_from_hdf5(this, dset_id)
class(Polynomial), intent(inout) :: this
integer(HID_T), intent(in) :: dset_id
integer(HSIZE_T) :: dims(1)
call get_shape(dset_id, dims)
allocate(this % coef(dims(1)))
call read_dataset(this % coef, dset_id)
end subroutine polynomial_from_hdf5
pure function polynomial_evaluate(this, x) result(y)
class(Polynomial), intent(in) :: this
real(8), intent(in) :: x
real(8) :: y
integer :: i
! Use Horner's rule to evaluate polynomial. Note that coefficients are
! ordered in increasing powers of x.
y = ZERO
do i = size(this % coef), 1, -1
y = y*x + this % coef(i)
end do
end function polynomial_evaluate
!===============================================================================
! Tabulated1D implementation
!===============================================================================
subroutine tabulated1d_from_hdf5(this, dset_id)
class(Tabulated1D), intent(inout) :: this
integer(HID_T), intent(in) :: dset_id
real(8), allocatable :: xy(:,:)
integer(HSIZE_T) :: dims(2)
call read_attribute(this % nbt, dset_id, 'breakpoints')
call read_attribute(this % int, dset_id, 'interpolation')
this % n_regions = size(this % nbt)
call get_shape(dset_id, dims)
this % n_pairs = int(dims(1), 4)
allocate(this % x(this % n_pairs))
allocate(this % y(this % n_pairs))
allocate(xy(dims(1), dims(2)))
call read_dataset(xy, dset_id)
this % x(:) = xy(:,1)
this % y(:) = xy(:,2)
end subroutine tabulated1d_from_hdf5
pure function tabulated1d_evaluate(this, x) result(y)
class(Tabulated1D), intent(in) :: this
real(8), intent(in) :: x ! x value to find y at
real(8) :: y ! y(x)
integer :: i ! bin in which to interpolate
integer :: j ! index for interpolation region
integer :: n_regions ! number of interpolation regions
integer :: n_pairs ! number of tabulated values
integer :: interp ! ENDF interpolation scheme
real(8) :: r ! interpolation factor
real(8) :: x0, x1 ! bounding x values
real(8) :: y0, y1 ! bounding y values
! determine number of interpolation regions and pairs
n_regions = this % n_regions
n_pairs = this % n_pairs
! find which bin the abscissa is in -- if the abscissa is outside the
! tabulated range, the first or last point is chosen, i.e. no interpolation
! is done outside the energy range
if (x < this % x(1)) then
y = this % y(1)
return
elseif (x > this % x(n_pairs)) then
y = this % y(n_pairs)
return
else
i = binary_search(this % x, n_pairs, x)
end if
! determine interpolation scheme
if (n_regions == 0) then
interp = LINEAR_LINEAR
elseif (n_regions == 1) then
interp = this % int(1)
elseif (n_regions > 1) then
do j = 1, n_regions
if (i < this % nbt(j)) then
interp = this % int(j)
exit
end if
end do
end if
! handle special case of histogram interpolation
if (interp == HISTOGRAM) then
y = this % y(i)
return
end if
! determine bounding values
x0 = this % x(i)
x1 = this % x(i + 1)
y0 = this % y(i)
y1 = this % y(i + 1)
! determine interpolation factor and interpolated value
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function tabulated1d_evaluate
end module endf_header

26
src/error.F90
View file

@ -111,14 +111,6 @@ contains
end subroutine warning
subroutine warning_from_c(message, message_len) bind(C)
integer(C_INT), intent(in), value :: message_len
character(kind=C_CHAR), intent(in) :: message(message_len)
character(message_len+1) :: message_out
write(message_out, *) message
call warning(message_out)
end subroutine
!===============================================================================
! FATAL_ERROR alerts the user that an error has been encountered and displays a
! message about the particular problem. Errors are considered 'fatal' and hence
@ -208,14 +200,6 @@ contains
end subroutine fatal_error
subroutine fatal_error_from_c(message, message_len) bind(C)
integer(C_INT), intent(in), value :: message_len
character(kind=C_CHAR), intent(in) :: message(message_len)
character(message_len+1) :: message_out
write(message_out, *) message
call fatal_error(message_out)
end subroutine
!===============================================================================
! WRITE_MESSAGE displays an informational message to the log file and the
! standard output stream.
@ -271,14 +255,4 @@ contains
end subroutine write_message
subroutine write_message_from_c(message, message_len, level) bind(C)
integer(C_INT), intent(in), value :: message_len
character(kind=C_CHAR), intent(in) :: message(message_len)
integer(C_INT), intent(in), value :: level
character(message_len+1) :: message_out
! Using * in the internal write adds an extra space at the beginning
write(message_out, *) message
call write_message(message_out(2:), level)
end subroutine write_message_from_c
end module error

101
src/error.cpp
View file

@ -1,14 +1,113 @@
#include "openmc/error.h"
#include "openmc/message_passing.h"
#include "openmc/settings.h"
#if defined(__unix__) || defined(__unix) || (defined(__APPLE__) && defined(__MACH__))
#include <unistd.h> // for isatty
#endif
#include <cstdlib> // for exit
#include <iomanip> // for setw
#include <iostream>
namespace openmc {
#ifdef OPENMC_MPI
void abort_mpi(int code)
{
MPI_Abort(mpi::intracomm, code);
MPI_Abort(mpi::intracomm, code);
}
#endif
void output(const std::string& message, std::ostream& out, int indent)
{
// Set line wrapping and indentation
int line_wrap = 80;
// Determine length of message
int length = message.size();
int i_start = 0;
int line_len = line_wrap - indent + 1;
while (i_start < length) {
if (length - i_start < line_len) {
// Remainder of message will fit on line
std::cerr << message.substr(i_start) << '\n';
break;
} else {
// Determine last space in current line
std::string s = message.substr(i_start, line_len);
auto pos = s.find_last_of(' ');
// Write up to last space, or whole line if no space is present
std::cerr << s.substr(0, pos) << '\n' << std::setw(indent) << " ";
// Advance starting position
i_start += (pos == std::string::npos) ? line_len : pos + 1;
}
}
}
void warning(const std::string& message)
{
#ifdef _POSIX_VERSION
// Make output yellow if user is in a terminal
if (isatty(STDERR_FILENO)) {
std::cerr << "\033[0;33m";
}
#endif
// Write warning
std::cerr << " WARNING: ";
output(message, std::cerr, 10);
#ifdef _POSIX_VERSION
// Reset color for terminal
if (isatty(STDERR_FILENO)) {
std::cerr << "\033[0m";
}
#endif
}
void write_message(const std::string& message, int level)
{
// Only allow master to print to screen
if (!mpi::master) return;
if (level <= settings::verbosity) {
std::cout << " ";
output(message, std::cout, 1);
}
}
void fatal_error(const std::string& message, int err)
{
#ifdef _POSIX_VERSION
// Make output red if user is in a terminal
if (isatty(STDERR_FILENO)) {
std::cerr << "\033[0;31m";
}
#endif
// Write error message
std::cerr << " ERROR: ";
output(message, std::cerr, 8);
#ifdef _POSIX_VERSION
// Reset color for terminal
if (isatty(STDERR_FILENO)) {
std::cerr << "\033[0m";
}
#endif
#ifdef OPENMC_MPI
MPI_Abort(mpi::intracomm, err);
#endif
// Abort the program
std::exit(err);
}
} // namespace openmc

4
src/geometry_header.F90
View file

@ -2,12 +2,8 @@ module geometry_header
use, intrinsic :: ISO_C_BINDING
use algorithm, only: find
use constants, only: K_BOLTZMANN, MATERIAL_VOID
use dict_header, only: DictIntInt
use nuclide_header
use stl_vector, only: VectorReal
use string, only: to_lower
implicit none

16
src/initialize.F90
View file

@ -20,14 +20,6 @@ module initialize
import C_PTR
type(C_PTR) :: ptr
end function
function path_statepoint_c() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
function path_sourcepoint_c() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
end interface
contains
@ -55,14 +47,6 @@ contains
else
path_input = ''
end if
if (.not. is_null(path_statepoint_c())) then
call c_f_pointer(path_statepoint_c(), string, [255])
path_state_point = to_f_string(string)
end if
if (.not. is_null(path_sourcepoint_c())) then
call c_f_pointer(path_sourcepoint_c(), string, [255])
path_source_point = to_f_string(string)
end if
if (.not. is_null(path_particle_restart_c())) then
call c_f_pointer(path_particle_restart_c(), string, [255])
path_particle_restart = to_f_string(string)

2
src/initialize.cpp
View file

@ -30,8 +30,6 @@
#include "openmc/timer.h"
// data/functions from Fortran side
extern "C" void print_usage();
extern "C" void print_version();
extern "C" void read_command_line();
extern "C" void read_geometry_xml();
extern "C" void read_materials_xml();

6
src/input_xml.F90
View file

@ -17,15 +17,13 @@ module input_xml
use message_passing
use mgxs_interface
use nuclide_header
use output, only: title, header
use photon_header
use random_lcg, only: prn
use surface_header
use settings
use stl_vector, only: VectorInt, VectorReal, VectorChar
use string, only: to_lower, to_str, str_to_int, str_to_real, &
starts_with, ends_with, split_string, &
zero_padded, to_c_string
use string, only: to_lower, to_str, str_to_int, &
starts_with, ends_with, to_c_string
use tally
use tally_header, only: openmc_extend_tallies
use tally_derivative_header

50
src/math.F90
View file

@ -5,64 +5,16 @@ module math
implicit none
private
public :: t_percentile
public :: calc_pn
public :: calc_rn
public :: rotate_angle
interface
pure function t_percentile(p, df) bind(C, name='t_percentile_c') &
result(t)
pure function t_percentile(p, df) bind(C) result(t)
use ISO_C_BINDING
implicit none
real(C_DOUBLE), value, intent(in) :: p
integer(C_INT), value, intent(in) :: df
real(C_DOUBLE) :: t
end function t_percentile
pure subroutine calc_pn(n, x, pnx) bind(C, name='calc_pn_c')
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n
real(C_DOUBLE), value, intent(in) :: x
real(C_DOUBLE), intent(out) :: pnx(n + 1)
end subroutine calc_pn
pure subroutine calc_rn(n, uvw, rn) bind(C, name='calc_rn_c')
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n
real(C_DOUBLE), intent(in) :: uvw(3)
real(C_DOUBLE), intent(out) :: rn(2 * n + 1)
end subroutine calc_rn
subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c')
use ISO_C_BINDING
implicit none
real(C_DOUBLE), intent(inout) :: uvw(3)
real(C_DOUBLE), value, intent(in) :: mu
real(C_DOUBLE), optional, intent(in) :: phi
end subroutine rotate_angle_c_intfc
end interface
contains
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done
! with direct sampling rather than rejection as is done in MCNP and SERPENT.
!===============================================================================
function rotate_angle(uvw0, mu, phi) result(uvw)
real(C_DOUBLE), intent(in) :: uvw0(3) ! directional cosine
real(C_DOUBLE), intent(in) :: mu ! cosine of angle in lab or CM
real(C_DOUBLE), intent(in), optional :: phi ! azimuthal angle
real(C_DOUBLE) :: uvw(3) ! rotated directional cosine
uvw = uvw0
call rotate_angle_c_intfc(uvw, mu, phi)
end function rotate_angle
end module math

2
src/math_functions.cpp
View file

@ -60,7 +60,7 @@ double normal_percentile(double p) {
}
double t_percentile_c(double p, int df){
double t_percentile(double p, int df){
double t;
if (df == 1) {

3
src/nuclide_header.F90
View file

@ -3,10 +3,9 @@ module nuclide_header
use, intrinsic :: ISO_FORTRAN_ENV
use, intrinsic :: ISO_C_BINDING
use algorithm, only: sort, find, binary_search
use algorithm, only: sort, find
use constants
use endf, only: is_fission, is_disappearance
use endf_header, only: Function1D, Polynomial, Tabulated1D
use error
use hdf5_interface
use message_passing

351
src/output.F90
View file

@ -22,7 +22,6 @@ module output
use tally_filter
use tally_filter_mesh, only: MeshFilter
use tally_filter_header, only: TallyFilterMatch
use timer_header
implicit none
@ -31,12 +30,6 @@ module output
integer :: eu = ERROR_UNIT
interface
function entropy(i) result(h) bind(C, name='entropy_c')
import C_INT, C_DOUBLE
integer(C_INT), value :: i
real(C_DOUBLE) :: h
end function
subroutine print_particle(p) bind(C)
import Particle
type(Particle), intent(in) :: p
@ -45,73 +38,6 @@ module output
contains
!===============================================================================
! TITLE prints the main title banner as well as information about the program
! developers, version, and date/time which the problem was run.
!===============================================================================
subroutine title() bind(C)
#ifdef _OPENMP
use omp_lib
#endif
write(UNIT=OUTPUT_UNIT, FMT='(/23(A/))') &
' %%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%', &
' ############### %%%%%%%%%%%%%%%%%%%%%%%%', &
' ################## %%%%%%%%%%%%%%%%%%%%%%%', &
' ################### %%%%%%%%%%%%%%%%%%%%%%%', &
' #################### %%%%%%%%%%%%%%%%%%%%%%', &
' ##################### %%%%%%%%%%%%%%%%%%%%%', &
' ###################### %%%%%%%%%%%%%%%%%%%%', &
' ####################### %%%%%%%%%%%%%%%%%%', &
' ####################### %%%%%%%%%%%%%%%%%', &
' ###################### %%%%%%%%%%%%%%%%%', &
' #################### %%%%%%%%%%%%%%%%%', &
' ################# %%%%%%%%%%%%%%%%%', &
' ############### %%%%%%%%%%%%%%%%', &
' ############ %%%%%%%%%%%%%%%', &
' ######## %%%%%%%%%%%%%%', &
' %%%%%%%%%%%'
! Write version information
write(UNIT=OUTPUT_UNIT, FMT=*) &
' | The OpenMC Monte Carlo Code'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' Copyright | 2011-2018 MIT and OpenMC contributors'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' License | http://openmc.readthedocs.io/en/latest/license.html'
write(UNIT=OUTPUT_UNIT, FMT='(11X,"Version | ",I1,".",I2,".",I1)') &
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(10X,"Git SHA1 | ",A)') GIT_SHA1
#endif
! Write the date and time
write(UNIT=OUTPUT_UNIT, FMT='(9X,"Date/Time | ",A)') time_stamp()
#ifdef OPENMC_MPI
! Write number of processors
write(UNIT=OUTPUT_UNIT, FMT='(5X,"MPI Processes | ",A)') &
trim(to_str(n_procs))
#endif
#ifdef _OPENMP
! Write number of OpenMP threads
write(UNIT=OUTPUT_UNIT, FMT='(4X,"OpenMP Threads | ",A)') &
trim(to_str(omp_get_max_threads()))
#endif
write(UNIT=OUTPUT_UNIT, FMT=*)
end subroutine title
!===============================================================================
! TIME_STAMP returns the current date and time in a formatted string
!===============================================================================
@ -166,283 +92,6 @@ contains
end subroutine header
!===============================================================================
! PRINT_VERSION shows the current version as well as copright and license
! information
!===============================================================================
subroutine print_version() bind(C)
if (master) then
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I2,".",I1)') &
"OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,A)') "Git SHA1: ", GIT_SHA1
#endif
write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2018 &
&Massachusetts Institute of Technology and OpenMC contributors"
write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at &
&<http://openmc.readthedocs.io/en/latest/license.html>"
end if
end subroutine print_version
!===============================================================================
! PRINT_USAGE displays information about command line usage of OpenMC
!===============================================================================
subroutine print_usage() bind(C)
if (master) then
write(OUTPUT_UNIT,*) 'Usage: openmc [options] [directory]'
write(OUTPUT_UNIT,*)
write(OUTPUT_UNIT,*) 'Options:'
write(OUTPUT_UNIT,*) ' -c, --volume Run in stochastic volume calculation mode'
write(OUTPUT_UNIT,*) ' -g, --geometry-debug Run with geometry debugging on'
write(OUTPUT_UNIT,*) ' -n, --particles Number of particles per generation'
write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode'
write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run from a state point'
write(OUTPUT_UNIT,*) ' or a particle restart file'
write(OUTPUT_UNIT,*) ' -s, --threads Number of OpenMP threads'
write(OUTPUT_UNIT,*) ' -t, --track Write tracks for all particles'
write(OUTPUT_UNIT,*) ' -v, --version Show version information'
write(OUTPUT_UNIT,*) ' -h, --help Show this message'
end if
end subroutine print_usage
!===============================================================================
! PRINT_COLUMNS displays a header listing what physical values will displayed
! below them
!===============================================================================
subroutine print_columns() bind(C)
write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "Bat./Gen."
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') " k "
if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "Entropy "
write(UNIT=ou, FMT='(A20,3X)', ADVANCE='NO') " Average k "
write(UNIT=ou, FMT=*)
write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "========="
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
write(UNIT=ou, FMT='(A20,3X)', ADVANCE='NO') "===================="
write(UNIT=ou, FMT=*)
end subroutine print_columns
!===============================================================================
! PRINT_GENERATION displays information for a generation of neutrons.
!===============================================================================
subroutine print_generation() bind(C)
integer :: i ! overall generation
integer :: n ! number of active generations
! 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
! write out information about batch and generation
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
trim(to_str(current_batch)) // "/" // trim(to_str(current_gen))
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation(i)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(i)
if (n > 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
keff, keff_std
end if
! next line
write(UNIT=OUTPUT_UNIT, FMT=*)
end subroutine print_generation
!===============================================================================
! PRINT_BATCH_KEFF displays the last batch's tallied value of the neutron
! multiplication factor as well as the average value if we're in active batches
!===============================================================================
subroutine print_batch_keff() bind(C)
integer :: i ! overall generation
integer :: n ! number of active generations
! Determine overall generation and number of active generations
i = current_batch*gen_per_batch
n = n_realizations*gen_per_batch
! write out information batch and option independent output
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
trim(to_str(current_batch)) // "/" // trim(to_str(gen_per_batch))
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') &
k_generation(i)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(i)
! write out accumulated k-effective if after first active batch
if (n > 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
keff, keff_std
else
write(UNIT=OUTPUT_UNIT, FMT='(23X)', ADVANCE='NO')
end if
! next line
write(UNIT=OUTPUT_UNIT, FMT=*)
end subroutine print_batch_keff
!===============================================================================
! PRINT_RUNTIME displays the total time elapsed for the entire run, for
! initialization, for computation, and for intergeneration synchronization.
!===============================================================================
subroutine print_runtime() bind(C)
integer :: n_active
real(8) :: speed_inactive ! # of neutrons/second in inactive batches
real(8) :: speed_active ! # of neutrons/second in active batches
character(15) :: string
! display header block
call header("Timing Statistics", 6)
! display time elapsed for various sections
write(ou,100) "Total time for initialization", time_initialize_elapsed()
write(ou,100) " Reading cross sections", time_read_xs_elapsed()
write(ou,100) "Total time in simulation", time_inactive_elapsed() + &
time_active_elapsed()
write(ou,100) " Time in transport only", time_transport_elapsed()
if (run_mode == MODE_EIGENVALUE) then
write(ou,100) " Time in inactive batches", time_inactive_elapsed()
end if
write(ou,100) " Time in active batches", time_active_elapsed()
if (run_mode == MODE_EIGENVALUE) then
write(ou,100) " Time synchronizing fission bank", time_bank_elapsed()
write(ou,100) " Sampling source sites", time_bank_sample_elapsed()
write(ou,100) " SEND/RECV source sites", time_bank_sendrecv_elapsed()
end if
write(ou,100) " Time accumulating tallies", time_tallies_elapsed()
write(ou,100) "Total time for finalization", time_finalize_elapsed()
write(ou,100) "Total time elapsed", time_total_elapsed()
! Calculate particle rate in active/inactive batches
n_active = current_batch - n_inactive
if (restart_run) then
if (restart_batch < n_inactive) then
speed_inactive = real(n_particles * (n_inactive - restart_batch) * &
gen_per_batch) / time_inactive_elapsed()
speed_active = real(n_particles * n_active * gen_per_batch) / &
time_active_elapsed()
else
speed_inactive = ZERO
speed_active = real(n_particles * (n_batches - restart_batch) * &
gen_per_batch) / time_active_elapsed()
end if
else
if (n_inactive > 0) then
speed_inactive = real(n_particles * n_inactive * gen_per_batch) / &
time_inactive_elapsed()
end if
speed_active = real(n_particles * n_active * gen_per_batch) / &
time_active_elapsed()
end if
! display calculation rate
if (.not. (restart_run .and. (restart_batch >= n_inactive)) &
.and. n_inactive > 0) then
string = to_str(speed_inactive)
write(ou,101) "Calculation Rate (inactive)", trim(string)
end if
string = to_str(speed_active)
write(ou,101) "Calculation Rate (active)", trim(string)
! format for write statements
100 format (1X,A,T36,"= ",ES11.4," seconds")
101 format (1X,A,T36,"= ",A," particles/second")
end subroutine print_runtime
!===============================================================================
! PRINT_RESULTS displays various estimates of k-effective as well as the global
! leakage rate.
!===============================================================================
subroutine print_results() bind(C)
integer :: n ! number of realizations
real(8) :: alpha ! significance level for CI
real(8) :: t_n1 ! t-value with N-1 degrees of freedom
real(8) :: t_n3 ! t-value with N-3 degrees of freedom
real(8) :: x(2) ! mean and standard deviation
real(C_DOUBLE) :: k_combined(2)
integer(C_INT) :: err
! display header block for results
call header("Results", 4)
n = n_realizations
if (confidence_intervals) then
! Calculate t-value for confidence intervals
alpha = ONE - CONFIDENCE_LEVEL
t_n1 = t_percentile(ONE - alpha/TWO, n - 1)
t_n3 = t_percentile(ONE - alpha/TWO, n - 3)
else
t_n1 = ONE
t_n3 = ONE
end if
! write global tallies
if (n > 1) then
associate (r => global_tallies(RESULT_SUM:RESULT_SUM_SQ, :))
if (run_mode == MODE_EIGENVALUE) then
x(:) = mean_stdev(r(:, K_COLLISION), n)
write(ou,102) "k-effective (Collision)", x(1), t_n1 * x(2)
x(:) = mean_stdev(r(:, K_TRACKLENGTH), n)
write(ou,102) "k-effective (Track-length)", x(1), t_n1 * x(2)
x(:) = mean_stdev(r(:, K_ABSORPTION), n)
write(ou,102) "k-effective (Absorption)", x(1), t_n1 * x(2)
if (n > 3) then
err = openmc_get_keff(k_combined)
write(ou,102) "Combined k-effective", k_combined(1), &
t_n3 * k_combined(2)
end if
end if
x(:) = mean_stdev(r(:, LEAKAGE), n)
write(ou,102) "Leakage Fraction", x(1), t_n1 * x(2)
end associate
else
if (master) call warning("Could not compute uncertainties -- only one &
&active batch simulated!")
if (run_mode == MODE_EIGENVALUE) then
write(ou,103) "k-effective (Collision)", global_tallies(RESULT_SUM, K_COLLISION) / n
write(ou,103) "k-effective (Track-length)", global_tallies(RESULT_SUM, K_TRACKLENGTH) / n
write(ou,103) "k-effective (Absorption)", global_tallies(RESULT_SUM, K_ABSORPTION) / n
end if
write(ou,103) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE) / n
end if
write(ou,*)
102 format (1X,A,T30,"= ",F8.5," +/- ",F8.5)
103 format (1X,A,T30,"= ",F8.5)
end subroutine print_results
!===============================================================================
! WRITE_TALLIES creates an output file and writes out the mean values of all
! tallies and their standard deviations

368
src/output.cpp
View file

@ -2,25 +2,91 @@
#include <algorithm> // for std::transform
#include <cstring> // for strlen
#include <iomanip> // for setw
#include <ctime> // for time, localtime
#include <iomanip> // for setw, setprecision, put_time
#include <ios> // for fixed, scientific, left
#include <iostream>
#include <sstream>
#include <ctime>
#include <utility> // for pair
#include <omp.h>
#include "xtensor/xview.hpp"
#include "openmc/capi.h"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/lattice.h"
#include "openmc/math_functions.h"
#include "openmc/message_passing.h"
#include "openmc/plot.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/surface.h"
#include "openmc/timer.h"
#include "openmc/tallies/tally.h"
namespace openmc {
//==============================================================================
void title()
{
std::cout <<
" %%%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%%%%%%%%%%%%%%\n" <<
" ############### %%%%%%%%%%%%%%%%%%%%%%%%\n" <<
" ################## %%%%%%%%%%%%%%%%%%%%%%%\n" <<
" ################### %%%%%%%%%%%%%%%%%%%%%%%\n" <<
" #################### %%%%%%%%%%%%%%%%%%%%%%\n" <<
" ##################### %%%%%%%%%%%%%%%%%%%%%\n" <<
" ###################### %%%%%%%%%%%%%%%%%%%%\n" <<
" ####################### %%%%%%%%%%%%%%%%%%\n" <<
" ####################### %%%%%%%%%%%%%%%%%\n" <<
" ###################### %%%%%%%%%%%%%%%%%\n" <<
" #################### %%%%%%%%%%%%%%%%%\n" <<
" ################# %%%%%%%%%%%%%%%%%\n" <<
" ############### %%%%%%%%%%%%%%%%\n" <<
" ############ %%%%%%%%%%%%%%%\n" <<
" ######## %%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%\n";
// Write version information
std::cout <<
" | The OpenMC Monte Carlo Code\n" <<
" Copyright | 2011-2019 MIT and OpenMC contributors\n" <<
" License | http://openmc.readthedocs.io/en/latest/license.html\n" <<
" Version | " << VERSION_MAJOR << '.' << VERSION_MINOR << '.'
<< VERSION_RELEASE << '\n';
#ifdef GIT_SHA1
std::cout << " Git SHA1 | " << GIT_SHA1 << '\n';
#endif
// Write the date and time
std::cout << " Date/Time | " << time_stamp() << '\n';
#ifdef OPENMC_MPI
// Write number of processors
std::cout << " MPI Processes | " << mpi::n_procs << '\n';
#endif
#ifdef _OPENMP
// Write number of OpenMP threads
std::cout << " OpenMC Threads | " << omp_get_max_threads() << '\n';
#endif
std::cout << '\n';
}
//==============================================================================
void
header(const char* msg, int level) {
// Determine how many times to repeat the '=' character.
@ -41,7 +107,7 @@ header(const char* msg, int level) {
// Print header based on verbosity level.
if (settings::verbosity >= level) {
std::cout << out.str() << "\n\n";
std::cout << '\n' << out.str() << "\n\n";
}
}
@ -49,13 +115,9 @@ header(const char* msg, int level) {
std::string time_stamp()
{
int base_year = 1990;
std::stringstream ts;
std::time_t t = std::time(0); // get time now
std::tm* now = std::localtime(&t);
ts << now->tm_year + base_year << "-" << now->tm_mon
<< "-" << now->tm_mday << " " << now->tm_hour
<< ":" << now->tm_min << ":" << now->tm_sec;
std::time_t t = std::time(nullptr); // get time now
ts << std::put_time(std::localtime(&t), "%Y-%m-%d %H:%M:%S");
return ts.str();
}
@ -230,4 +292,292 @@ print_overlap_check()
}
}
//==============================================================================
void print_usage()
{
if (mpi::master) {
std::cout <<
"Usage: openmc [options] [directory]\n\n"
"Options:\n"
" -c, --volume Run in stochastic volume calculation mode\n"
" -g, --geometry-debug Run with geometry debugging on\n"
" -n, --particles Number of particles per generation\n"
" -p, --plot Run in plotting mode\n"
" -r, --restart Restart a previous run from a state point\n"
" or a particle restart file\n"
" -s, --threads Number of OpenMP threads\n"
" -t, --track Write tracks for all particles\n"
" -v, --version Show version information\n"
" -h, --help Show this message\n";
}
}
//==============================================================================
void print_version()
{
if (mpi::master) {
std::cout << "OpenMC version " << VERSION_MAJOR << '.' << VERSION_MINOR
<< '.' << VERSION_RELEASE << '\n';
#ifdef GIT_SHA1
std::cout << "Git SHA1: " << GIT_SHA1 << '\n';
#endif
std::cout << "Copyright (c) 2011-2019 Massachusetts Institute of "
"Technology and OpenMC contributors\nMIT/X license at "
"<http://openmc.readthedocs.io/en/latest/license.html>\n";
}
}
//==============================================================================
void print_columns()
{
if (settings::entropy_on) {
std::cout <<
" Bat./Gen. k Entropy Average k \n"
" ========= ======== ======== ====================\n";
} else {
std::cout <<
" Bat./Gen. k Average k\n"
" ========= ======== ====================\n";
}
}
//==============================================================================
void print_generation()
{
// Save state of cout
auto f {std::cout.flags()};
// Determine overall generation and number of active generations
int i = overall_generation() - 1;
int n = simulation::current_batch > settings::n_inactive ?
settings::gen_per_batch*n_realizations + simulation::current_gen : 0;
// Set format for values
std::cout << std::fixed << std::setprecision(5);
// write out information batch and option independent output
std::cout << " " << std::setw(9) << std::to_string(simulation::current_batch)
+ "/" + std::to_string(simulation::current_gen) << " " << std::setw(8)
<< simulation::k_generation[i];
// write out entropy info
if (settings::entropy_on) {
std::cout << " " << std::setw(8) << simulation::entropy[i];
}
if (n > 1) {
std::cout << " " << std::setw(8) << simulation::keff << " +/-"
<< std::setw(8) << simulation::keff_std;
}
std::cout << '\n';
// Restore state of cout
std::cout.flags(f);
}
//==============================================================================
void print_batch_keff()
{
// Save state of cout
auto f {std::cout.flags()};
// Determine overall generation and number of active generations
int i = simulation::current_batch*settings::gen_per_batch - 1;
int n = n_realizations*settings::gen_per_batch;
// Set format for values
std::cout << std::fixed << std::setprecision(5);
// write out information batch and option independent output
std::cout << " " << std::setw(9) << std::to_string(simulation::current_batch)
+ "/" + std::to_string(settings::gen_per_batch) << " " << std::setw(8)
<< simulation::k_generation[i];
// write out entropy info
if (settings::entropy_on) {
std::cout << " " << std::setw(8) << simulation::entropy[i];
}
if (n > 1) {
std::cout << " " << std::setw(8) << simulation::keff << " +/-"
<< std::setw(8) << simulation::keff_std;
}
std::cout << '\n';
// Restore state of cout
std::cout.flags(f);
}
//==============================================================================
void show_time(const char* label, double secs, int indent_level=0)
{
std::cout << std::string(2*indent_level, ' ');
int width = 33 - indent_level*2;
std::cout << " " << std::setw(width) << std::left << label << " = "
<< std::setw(10) << std::right << secs << " seconds\n";
}
void show_rate(const char* label, double particles_per_sec)
{
std::cout << " " << std::setw(33) << std::left << label << " = " <<
particles_per_sec << " particles/second\n";
}
void print_runtime()
{
using namespace simulation;
// display header block
header("Timing Statistics", 6);
// Save state of cout
auto f {std::cout.flags()};
// display time elapsed for various sections
std::cout << std::scientific << std::setprecision(4);
show_time("Total time for initialization", time_initialize.elapsed());
show_time("Reading cross sections", time_read_xs.elapsed(), 1);
show_time("Total time in simulation", time_inactive.elapsed() +
time_active.elapsed());
show_time("Time in transport only", time_transport.elapsed(), 1);
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
show_time("Time in inactive batches", time_inactive.elapsed(), 1);
}
show_time("Time in active batches", time_active.elapsed(), 1);
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
show_time("Time synchronizing fission bank", time_bank.elapsed(), 1);
show_time("Sampling source sites", time_bank_sample.elapsed(), 2);
show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2);
}
show_time("Time accumulating tallies", time_tallies.elapsed(), 1);
show_time("Total time for finalization", time_finalize.elapsed());
show_time("Total time elapsed", time_total.elapsed());
// Restore state of cout
std::cout.flags(f);
// Calculate particle rate in active/inactive batches
int n_active = simulation::current_batch - settings::n_inactive;
double speed_inactive;
double speed_active;
if (settings::restart_run) {
if (simulation::restart_batch < settings::n_inactive) {
speed_inactive = (settings::n_particles * (settings::n_inactive
- simulation::restart_batch) * settings::gen_per_batch)
/ time_inactive.elapsed();
speed_active = (settings::n_particles * n_active
* settings::gen_per_batch) / time_active.elapsed();
} else {
speed_inactive = 0.0;
speed_active = (settings::n_particles * (settings::n_batches
- simulation::restart_batch) * settings::gen_per_batch)
/ time_active.elapsed();
}
} else {
if (settings::n_inactive > 0) {
speed_inactive = (settings::n_particles * settings::n_inactive
* settings::gen_per_batch) / time_inactive.elapsed();
}
speed_active = (settings::n_particles * n_active * settings::gen_per_batch)
/ time_active.elapsed();
}
// display calculation rate
std::cout << std::setprecision(6) << std::showpoint;
if (!(settings::restart_run && (simulation::restart_batch >= settings::n_inactive))
&& settings::n_inactive > 0) {
show_rate("Calculation Rate (inactive)", speed_inactive);
}
show_rate("Calculation Rate (active)", speed_active);
// Restore state of cout
std::cout.flags(f);
}
//==============================================================================
std::pair<double, double> mean_stdev(const double* x, int n)
{
double mean = x[RESULT_SUM] / n;
double stdev = n > 1 ? std::sqrt((x[RESULT_SUM_SQ]/n
- mean*mean)/(n - 1)) : 0.0;
return {mean, stdev};
}
//==============================================================================
void print_results()
{
// Save state of cout
auto f {std::cout.flags()};
// display header block for results
header("Results", 4);
// Calculate t-value for confidence intervals
int n = n_realizations;
double alpha, t_n1, t_n3;
if (settings::confidence_intervals) {
alpha = 1.0 - CONFIDENCE_LEVEL;
t_n1 = t_percentile(1.0 - alpha/2.0, n - 1);
t_n3 = t_percentile(1.0 - alpha/2.0, n - 3);
} else {
t_n1 = 1.0;
t_n3 = 1.0;
}
// Set formatting for floats
std::cout << std::fixed << std::setprecision(5);
// write global tallies
auto gt = global_tallies();
double mean, stdev;
if (n > 1) {
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
std::tie(mean, stdev) = mean_stdev(&gt(K_COLLISION, 0), n);
std::cout << " k-effective (Collision) = "
<< mean << " +/- " << t_n1 * stdev << '\n';
std::tie(mean, stdev) = mean_stdev(&gt(K_TRACKLENGTH, 0), n);
std::cout << " k-effective (Track-length) = "
<< mean << " +/- " << t_n1 * stdev << '\n';
std::tie(mean, stdev) = mean_stdev(&gt(K_ABSORPTION, 0), n);
std::cout << " k-effective (Absorption) = "
<< mean << " +/- " << t_n1 * stdev << '\n';
if (n > 3) {
double k_combined[2];
openmc_get_keff(k_combined);
std::cout << " Combined k-effective = "
<< k_combined[0] << " +/- " << t_n3 * k_combined[1] << '\n';
}
}
std::tie(mean, stdev) = mean_stdev(&gt(LEAKAGE, 0), n);
std::cout << " Leakage Fraction = "
<< mean << " +/- " << t_n1 * stdev << '\n';
} else {
if (mpi::master) warning("Could not compute uncertainties -- only one "
"active batch simulated!");
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
std::cout << " k-effective (Collision) = "
<< gt(K_COLLISION, RESULT_SUM) / n << '\n';
std::cout << " k-effective (Track-length) = "
<< gt(K_TRACKLENGTH, RESULT_SUM) / n << '\n';
std::cout << " k-effective (Absorption) = "
<< gt(K_ABSORPTION, RESULT_SUM) / n << '\n';
}
std::cout << " Leakage Fraction = "
<< gt(LEAKAGE, RESULT_SUM) / n << '\n';
}
std::cout << '\n';
// Restore state of cout
std::cout.flags(f);
}
} // namespace openmc

25
src/random_lcg.F90
View file

@ -11,41 +11,16 @@ module random_lcg
real(C_DOUBLE) :: pseudo_rn
end function prn
function future_prn(n) result(pseudo_rn) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT64_T), value :: n
real(C_DOUBLE) :: pseudo_rn
end function future_prn
subroutine set_particle_seed(id) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT64_T), value :: id
end subroutine set_particle_seed
subroutine advance_prn_seed(n) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT64_T), value :: n
end subroutine advance_prn_seed
subroutine prn_set_stream(n) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value :: n
end subroutine prn_set_stream
function openmc_get_seed() result(seed) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT64_T) :: seed
end function openmc_get_seed
subroutine openmc_set_seed(new_seed) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT64_T), value :: new_seed
end subroutine openmc_set_seed
end interface
end module random_lcg

5
src/settings.F90
View file

@ -88,14 +88,9 @@ module settings
character(MAX_FILE_LEN) :: path_input ! Path to input file
character(MAX_FILE_LEN) :: path_cross_sections = '' ! Path to cross_sections.xml
character(MAX_FILE_LEN) :: path_state_point ! Path to binary state point
character(MAX_FILE_LEN) :: path_source_point ! Path to binary source point
character(MAX_FILE_LEN) :: path_particle_restart ! Path to particle restart
character(MAX_FILE_LEN) :: path_output = '' ! Path to output directory
! Various output options
logical(C_BOOL), bind(C) :: output_summary
! No reduction at end of batch
logical(C_BOOL), bind(C) :: reduce_tallies

3
src/settings.cpp
View file

@ -38,9 +38,9 @@ namespace settings {
// Default values for boolean flags
bool assume_separate {false};
bool check_overlaps {false};
bool cmfd_run {false};
bool confidence_intervals {false};
bool create_fission_neutrons {true};
bool dagmc {false};
bool entropy_on {false};
bool legendre_to_tabular {true};
bool output_summary {true};
@ -62,7 +62,6 @@ bool ufs_on {false};
bool urr_ptables_on {true};
bool write_all_tracks {false};
bool write_initial_source {false};
bool dagmc {false};
std::string path_cross_sections;
std::string path_input;

6
src/simulation.cpp
View file

@ -30,12 +30,6 @@ extern "C" void accumulate_tallies();
extern "C" void allocate_tally_results();
extern "C" void check_triggers();
extern "C" void init_tally_routines();
extern "C" void load_state_point();
extern "C" void print_batch_keff();
extern "C" void print_columns();
extern "C" void print_generation();
extern "C" void print_results();
extern "C" void print_runtime();
extern "C" void setup_active_tallies();
extern "C" void simulation_init_f();
extern "C" void simulation_finalize_f();

645
src/state_point.F90
View file

@ -13,7 +13,6 @@ module state_point
use, intrinsic :: ISO_C_BINDING
use bank_header, only: Bank
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message
@ -21,324 +20,187 @@ module state_point
use message_passing
use mgxs_interface
use nuclide_header, only: nuclides
use output, only: time_stamp
use random_lcg, only: openmc_get_seed, openmc_set_seed
use settings
use simulation_header
use string, only: to_str, count_digits, zero_padded, to_f_string
use string, only: to_str
use tally_header
use tally_filter_header
use tally_derivative_header, only: tally_derivs
use timer_header
implicit none
interface
subroutine write_source_bank(group_id) bind(C)
import HID_T
integer(HID_T), value :: group_id
end subroutine write_source_bank
subroutine read_source_bank(group_id) bind(C)
import HID_T
integer(HID_T), value :: group_id
end subroutine read_source_bank
end interface
contains
!===============================================================================
! OPENMC_STATEPOINT_WRITE writes an HDF5 statepoint file to disk
!===============================================================================
function openmc_statepoint_write(filename, write_source) result(err) bind(C)
type(C_PTR), value :: filename
logical(C_BOOL), intent(in), optional :: write_source
integer(C_INT) :: err
subroutine statepoint_write_f(file_id) bind(C)
integer(HID_T), value :: file_id
logical :: write_source_
integer :: i, j
integer :: i_xs
integer, allocatable :: id_array(:)
integer(HID_T) :: file_id
integer(HID_T) :: tallies_group, tally_group, &
filters_group, filter_group, derivs_group, &
deriv_group, runtime_group
deriv_group
character(MAX_WORD_LEN), allocatable :: str_array(:)
character(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: filename_
character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name
logical :: parallel
interface
subroutine meshes_to_hdf5(group) bind(C)
import HID_T
integer(HID_T), value :: group
end subroutine
subroutine write_eigenvalue_hdf5(group) bind(C)
import HID_T
integer(HID_T), value :: group
end subroutine
subroutine write_tally_results_nr(file_id) bind(C)
import HID_T
integer(HID_T), value :: file_id
end subroutine
end interface
! Open tallies group
tallies_group = open_group(file_id, "tallies")
err = 0
! Set the filename
if (c_associated(filename)) then
call c_f_pointer(filename, string, [MAX_FILE_LEN])
filename_ = to_f_string(string)
else
! Set filename for state point
filename_ = trim(path_output) // 'statepoint.' // &
& zero_padded(current_batch, count_digits(n_max_batches))
filename_ = trim(filename_) // '.h5'
! Write information for derivatives.
if (size(tally_derivs) > 0) then
derivs_group = create_group(tallies_group, "derivatives")
do i = 1, size(tally_derivs)
associate(deriv => tally_derivs(i))
deriv_group = create_group(derivs_group, "derivative " &
// trim(to_str(deriv % id)))
select case (deriv % variable)
case (DIFF_DENSITY)
call write_dataset(deriv_group, "independent variable", "density")
call write_dataset(deriv_group, "material", deriv % diff_material)
case (DIFF_NUCLIDE_DENSITY)
call write_dataset(deriv_group, "independent variable", &
"nuclide_density")
call write_dataset(deriv_group, "material", deriv % diff_material)
call write_dataset(deriv_group, "nuclide", &
nuclides(deriv % diff_nuclide) % name)
case (DIFF_TEMPERATURE)
call write_dataset(deriv_group, "independent variable", &
"temperature")
call write_dataset(deriv_group, "material", deriv % diff_material)
case default
call fatal_error("Independent variable for derivative " &
// trim(to_str(deriv % id)) // " not defined in &
&state_point.F90.")
end select
call close_group(deriv_group)
end associate
end do
call close_group(derivs_group)
end if
! Determine whether or not to write the source bank
if (present(write_source)) then
write_source_ = write_source
else
write_source_ = .true.
! Write number of filters
filters_group = create_group(tallies_group, "filters")
call write_attribute(filters_group, "n_filters", n_filters)
if (n_filters > 0) then
! Write IDs of filters
allocate(id_array(n_filters))
do i = 1, n_filters
id_array(i) = filters(i) % obj % id
end do
call write_attribute(filters_group, "ids", id_array)
deallocate(id_array)
! Write filter information
FILTER_LOOP: do i = 1, n_filters
filter_group = create_group(filters_group, "filter " // &
trim(to_str(filters(i) % obj % id)))
call filters(i) % obj % to_statepoint(filter_group)
call close_group(filter_group)
end do FILTER_LOOP
end if
! Write message
call write_message("Creating state point " // trim(filename_) // "...", 5)
if (master) then
! Create statepoint file
file_id = file_open(filename_, 'w')
! Write file type
call write_attribute(file_id, "filetype", "statepoint")
! Write revision number for state point file
call write_attribute(file_id, "version", VERSION_STATEPOINT)
! Write OpenMC version
call write_attribute(file_id, "openmc_version", VERSION)
#ifdef GIT_SHA1
call write_attribute(file_id, "git_sha1", GIT_SHA1)
#endif
! Write current date and time
call write_attribute(file_id, "date_and_time", time_stamp())
! Write path to input
call write_attribute(file_id, "path", path_input)
! Write out random number seed
call write_dataset(file_id, "seed", openmc_get_seed())
! Write run information
if (run_CE) then
call write_dataset(file_id, "energy_mode", "continuous-energy")
else
call write_dataset(file_id, "energy_mode", "multi-group")
end if
select case(run_mode)
case (MODE_FIXEDSOURCE)
call write_dataset(file_id, "run_mode", "fixed source")
case (MODE_EIGENVALUE)
call write_dataset(file_id, "run_mode", "eigenvalue")
end select
if (photon_transport) then
call write_attribute(file_id, "photon_transport", 1)
else
call write_attribute(file_id, "photon_transport", 0)
end if
call write_dataset(file_id, "n_particles", n_particles)
call write_dataset(file_id, "n_batches", n_batches)
! Write out current batch number
call write_dataset(file_id, "current_batch", current_batch)
! Indicate whether source bank is stored in statepoint
if (write_source_) then
call write_attribute(file_id, "source_present", 1)
else
call write_attribute(file_id, "source_present", 0)
end if
! Write out information for eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call write_eigenvalue_hdf5(file_id)
end if
tallies_group = create_group(file_id, "tallies")
! Write meshes
call meshes_to_hdf5(tallies_group)
! Write information for derivatives.
if (size(tally_derivs) > 0) then
derivs_group = create_group(tallies_group, "derivatives")
do i = 1, size(tally_derivs)
associate(deriv => tally_derivs(i))
deriv_group = create_group(derivs_group, "derivative " &
// trim(to_str(deriv % id)))
select case (deriv % variable)
case (DIFF_DENSITY)
call write_dataset(deriv_group, "independent variable", "density")
call write_dataset(deriv_group, "material", deriv % diff_material)
case (DIFF_NUCLIDE_DENSITY)
call write_dataset(deriv_group, "independent variable", &
"nuclide_density")
call write_dataset(deriv_group, "material", deriv % diff_material)
call write_dataset(deriv_group, "nuclide", &
nuclides(deriv % diff_nuclide) % name)
case (DIFF_TEMPERATURE)
call write_dataset(deriv_group, "independent variable", &
"temperature")
call write_dataset(deriv_group, "material", deriv % diff_material)
case default
call fatal_error("Independent variable for derivative " &
// trim(to_str(deriv % id)) // " not defined in &
&state_point.F90.")
end select
call close_group(deriv_group)
end associate
end do
call close_group(derivs_group)
end if
! Write number of filters
filters_group = create_group(tallies_group, "filters")
call write_attribute(filters_group, "n_filters", n_filters)
if (n_filters > 0) then
! Write IDs of filters
allocate(id_array(n_filters))
do i = 1, n_filters
id_array(i) = filters(i) % obj % id
end do
call write_attribute(filters_group, "ids", id_array)
deallocate(id_array)
! Write filter information
FILTER_LOOP: do i = 1, n_filters
filter_group = create_group(filters_group, "filter " // &
trim(to_str(filters(i) % obj % id)))
call filters(i) % obj % to_statepoint(filter_group)
call close_group(filter_group)
end do FILTER_LOOP
end if
call close_group(filters_group)
call close_group(filters_group)
! Write number of tallies
call write_attribute(tallies_group, "n_tallies", n_tallies)
call write_attribute(tallies_group, "n_tallies", n_tallies)
if (n_tallies > 0) then
! Write array of tally IDs
allocate(id_array(n_tallies))
do i = 1, n_tallies
id_array(i) = tallies(i) % obj % id
end do
call write_attribute(tallies_group, "ids", id_array)
deallocate(id_array)
if (n_tallies > 0) then
! Write array of tally IDs
allocate(id_array(n_tallies))
do i = 1, n_tallies
id_array(i) = tallies(i) % obj % id
end do
call write_attribute(tallies_group, "ids", id_array)
deallocate(id_array)
! Write all tally information except results
TALLY_METADATA: do i = 1, n_tallies
! Write all tally information except results
TALLY_METADATA: do i = 1, n_tallies
! Get pointer to tally
associate (tally => tallies(i) % obj)
tally_group = create_group(tallies_group, "tally " // &
trim(to_str(tally % id)))
! Get pointer to tally
associate (tally => tallies(i) % obj)
tally_group = create_group(tallies_group, "tally " // &
trim(to_str(tally % id)))
! Write the name for this tally
call write_dataset(tally_group, "name", tally % name)
! Write the name for this tally
call write_dataset(tally_group, "name", tally % name)
select case(tally % estimator)
case (ESTIMATOR_ANALOG)
call write_dataset(tally_group, "estimator", "analog")
case (ESTIMATOR_TRACKLENGTH)
call write_dataset(tally_group, "estimator", "tracklength")
case (ESTIMATOR_COLLISION)
call write_dataset(tally_group, "estimator", "collision")
end select
call write_dataset(tally_group, "n_realizations", &
tally % n_realizations)
select case(tally % estimator)
case (ESTIMATOR_ANALOG)
call write_dataset(tally_group, "estimator", "analog")
case (ESTIMATOR_TRACKLENGTH)
call write_dataset(tally_group, "estimator", "tracklength")
case (ESTIMATOR_COLLISION)
call write_dataset(tally_group, "estimator", "collision")
end select
call write_dataset(tally_group, "n_realizations", &
tally % n_realizations)
call write_dataset(tally_group, "n_filters", size(tally % filter))
if (size(tally % filter) > 0) then
! Write IDs of filters
allocate(id_array(size(tally % filter)))
do j = 1, size(tally % filter)
id_array(j) = filters(tally % filter(j)) % obj % id
end do
call write_dataset(tally_group, "filters", id_array)
deallocate(id_array)
end if
call write_dataset(tally_group, "n_filters", size(tally % filter))
if (size(tally % filter) > 0) then
! Write IDs of filters
allocate(id_array(size(tally % filter)))
do j = 1, size(tally % filter)
id_array(j) = filters(tally % filter(j)) % obj % id
end do
call write_dataset(tally_group, "filters", id_array)
deallocate(id_array)
end if
! Set up nuclide bin array and then write
allocate(str_array(tally % n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
if (tally % nuclide_bins(j) > 0) then
if (run_CE) then
i_xs = index(nuclides(tally % nuclide_bins(j)) % name, '.')
if (i_xs > 0) then
str_array(j) = nuclides(tally % nuclide_bins(j)) % name(1 : i_xs-1)
else
str_array(j) = nuclides(tally % nuclide_bins(j)) % name
end if
! Set up nuclide bin array and then write
allocate(str_array(tally % n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
if (tally % nuclide_bins(j) > 0) then
if (run_CE) then
i_xs = index(nuclides(tally % nuclide_bins(j)) % name, '.')
if (i_xs > 0) then
str_array(j) = nuclides(tally % nuclide_bins(j)) % name(1 : i_xs-1)
else
call get_name_c(tally % nuclide_bins(j), len(temp_name), &
temp_name)
i_xs = index(temp_name, '.')
if (i_xs > 0) then
str_array(j) = trim(temp_name(1 : i_xs-1))
else
str_array(j) = trim(temp_name)
end if
str_array(j) = nuclides(tally % nuclide_bins(j)) % name
end if
else
str_array(j) = 'total'
call get_name_c(tally % nuclide_bins(j), len(temp_name), &
temp_name)
i_xs = index(temp_name, '.')
if (i_xs > 0) then
str_array(j) = trim(temp_name(1 : i_xs-1))
else
str_array(j) = trim(temp_name)
end if
end if
end do NUCLIDE_LOOP
call write_dataset(tally_group, "nuclides", str_array)
deallocate(str_array)
! Write derivative information.
if (tally % deriv /= NONE) then
call write_dataset(tally_group, "derivative", &
tally_derivs(tally % deriv) % id)
else
str_array(j) = 'total'
end if
end do NUCLIDE_LOOP
call write_dataset(tally_group, "nuclides", str_array)
deallocate(str_array)
! Write scores.
call write_dataset(tally_group, "n_score_bins", tally % n_score_bins)
allocate(str_array(size(tally % score_bins)))
do j = 1, size(tally % score_bins)
str_array(j) = reaction_name(tally % score_bins(j))
end do
call write_dataset(tally_group, "score_bins", str_array)
! Write derivative information.
if (tally % deriv /= NONE) then
call write_dataset(tally_group, "derivative", &
tally_derivs(tally % deriv) % id)
end if
deallocate(str_array)
! Write scores.
call write_dataset(tally_group, "n_score_bins", tally % n_score_bins)
allocate(str_array(size(tally % score_bins)))
do j = 1, size(tally % score_bins)
str_array(j) = reaction_name(tally % score_bins(j))
end do
call write_dataset(tally_group, "score_bins", str_array)
call close_group(tally_group)
end associate
end do TALLY_METADATA
deallocate(str_array)
end if
call close_group(tallies_group)
call close_group(tally_group)
end associate
end do TALLY_METADATA
end if
! Check for the no-tally-reduction method
if (.not. reduce_tallies) then
! If using the no-tally-reduction method, we need to collect tally
! results before writing them to the state point file.
call write_tally_results_nr(file_id)
elseif (master) then
! Write number of global realizations
call write_dataset(file_id, "n_realizations", n_realizations)
if (reduce_tallies) then
! Write global tallies
call write_dataset(file_id, "global_tallies", global_tallies)
@ -347,8 +209,6 @@ contains
! Indicate that tallies are on
call write_attribute(file_id, "tallies_present", 1)
tallies_group = open_group(file_id, "tallies")
! Write all tally results
TALLY_RESULTS: do i = 1, n_tallies
associate (tally => tallies(i) % obj)
@ -359,241 +219,52 @@ contains
call close_group(tally_group)
end associate
end do TALLY_RESULTS
call close_group(tallies_group)
else
! Indicate tallies are off
call write_attribute(file_id, "tallies_present", 0)
end if
! Write out the runtime metrics.
runtime_group = create_group(file_id, "runtime")
call write_dataset(runtime_group, "total initialization", &
time_initialize_elapsed())
call write_dataset(runtime_group, "reading cross sections", &
time_read_xs_elapsed())
call write_dataset(runtime_group, "simulation", &
time_inactive_elapsed() + time_active_elapsed())
call write_dataset(runtime_group, "transport", &
time_transport_elapsed())
if (run_mode == MODE_EIGENVALUE) then
call write_dataset(runtime_group, "inactive batches", &
time_inactive_elapsed())
end if
call write_dataset(runtime_group, "active batches", &
time_active_elapsed())
if (run_mode == MODE_EIGENVALUE) then
call write_dataset(runtime_group, "synchronizing fission bank", &
time_bank_elapsed())
call write_dataset(runtime_group, "sampling source sites", &
time_bank_sample_elapsed())
call write_dataset(runtime_group, "SEND-RECV source sites", &
time_bank_sendrecv_elapsed())
end if
call write_dataset(runtime_group, "accumulating tallies", &
time_tallies_elapsed())
call write_dataset(runtime_group, "total", time_total_elapsed())
call close_group(runtime_group)
call file_close(file_id)
end if
#ifdef PHDF5
parallel = .true.
#else
parallel = .false.
#endif
! Write the source bank if desired
if (write_source_) then
if (master .or. parallel) then
file_id = file_open(filename_, 'a', parallel=.true.)
end if
call write_source_bank(file_id)
if (master .or. parallel) call file_close(file_id)
end if
end function openmc_statepoint_write
call close_group(tallies_group)
end subroutine
!===============================================================================
! LOAD_STATE_POINT
!===============================================================================
subroutine load_state_point() bind(C)
subroutine load_state_point_f(file_id) bind(C)
integer(HID_T), value :: file_id
integer :: i
integer :: int_array(3)
integer, allocatable :: array(:)
integer(C_INT64_T) :: seed
integer(HID_T) :: file_id
integer :: temp
integer(HID_T) :: tallies_group
integer(HID_T) :: tally_group
logical :: source_present
character(MAX_WORD_LEN) :: word
interface
subroutine read_eigenvalue_hdf5(group) bind(C)
import HID_T
integer(HID_T), value :: group
end subroutine
! Read global tally data
call read_dataset(global_tallies, file_id, "global_tallies")
subroutine restart_set_keff() bind(C)
end subroutine
end interface
! Check if tally results are present
call read_attribute(temp, file_id, "tallies_present")
! Write message
call write_message("Loading state point " // trim(path_state_point) &
// "...", 5)
! Read in sum and sum squared
if (temp == 1) then
tallies_group = open_group(file_id, "tallies")
! Open file for reading
file_id = file_open(path_state_point, 'r', parallel=.true.)
TALLY_RESULTS: do i = 1, n_tallies
associate (t => tallies(i) % obj)
! Read sum, sum_sq, and N for each bin
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(t % id)))
call t % read_results_hdf5(tally_group)
call read_dataset(t % n_realizations, tally_group, &
"n_realizations")
call close_group(tally_group)
end associate
end do TALLY_RESULTS
! Read filetype
call read_attribute(word, file_id, "filetype")
if (word /= 'statepoint') then
call fatal_error("OpenMC tried to restart from a non-statepoint file.")
call close_group(tallies_group)
end if
! Read revision number for state point file and make sure it matches with
! current version
call read_attribute(array, file_id, "version")
if (any(array /= VERSION_STATEPOINT)) then
call fatal_error("State point version does not match current version &
&in OpenMC.")
end if
! Read and overwrite random number seed
call read_dataset(seed, file_id, "seed")
call openmc_set_seed(seed)
! It is not impossible for a state point to be generated from a CE run but
! to be loaded in to an MG run (or vice versa), check to prevent that.
call read_dataset(word, file_id, "energy_mode")
if (word == "multi-group" .and. run_CE) then
call fatal_error("State point file is from multi-group run but &
& current run is continous-energy!")
else if (word == "continuous-energy" .and. .not. run_CE) then
call fatal_error("State point file is from continuous-energy run but &
& current run is multi-group!")
end if
! Read and overwrite run information except number of batches
call read_dataset(word, file_id, "run_mode")
select case(word)
case ('fixed source')
run_mode = MODE_FIXEDSOURCE
case ('eigenvalue')
run_mode = MODE_EIGENVALUE
end select
call read_attribute(int_array(1), file_id, "photon_transport")
if (int_array(1) == 1) then
photon_transport = .true.
else
photon_transport = .false.
end if
call read_dataset(n_particles, file_id, "n_particles")
call read_dataset(int_array(1), file_id, "n_batches")
! Take maximum of statepoint n_batches and input n_batches
n_batches = max(n_batches, int_array(1))
! Read batch number to restart at
call read_dataset(restart_batch, file_id, "current_batch")
! Check for source in statepoint if needed
call read_attribute(int_array(1), file_id, "source_present")
if (int_array(1) == 1) then
source_present = .true.
else
source_present = .false.
end if
if (restart_batch > n_batches) then
call fatal_error("The number batches specified in settings.xml is fewer &
& than the number of batches in the given statepoint file.")
end if
! Read information specific to eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call read_dataset(int_array(1), file_id, "n_inactive")
call read_eigenvalue_hdf5(file_id)
! Take maximum of statepoint n_inactive and input n_inactive
n_inactive = max(n_inactive, int_array(1))
end if
! Read number of realizations for global tallies
call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.)
! Set k_sum, keff, and current_batch based on whether restart file is part
! of active cycle or inactive cycle
call restart_set_keff()
current_batch = restart_batch
! Check to make sure source bank is present
if (path_source_point == path_state_point .and. .not. source_present) then
call fatal_error("Source bank must be contained in statepoint restart &
&file")
end if
! Read tallies to master. If we are using Parallel HDF5, all processes
! need to be included in the HDF5 calls.
#ifdef PHDF5
if (.true.) then
#else
if (master) then
#endif
! Read global tally data
call read_dataset(global_tallies, file_id, "global_tallies")
! Check if tally results are present
call read_attribute(int_array(1), file_id, "tallies_present")
! Read in sum and sum squared
if (int_array(1) == 1) then
tallies_group = open_group(file_id, "tallies")
TALLY_RESULTS: do i = 1, n_tallies
associate (t => tallies(i) % obj)
! Read sum, sum_sq, and N for each bin
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(t % id)))
call t % read_results_hdf5(tally_group)
call read_dataset(t % n_realizations, tally_group, &
"n_realizations")
call close_group(tally_group)
end associate
end do TALLY_RESULTS
call close_group(tallies_group)
end if
end if
! Read source if in eigenvalue mode
if (run_mode == MODE_EIGENVALUE) then
! Check if source was written out separately
if (.not. source_present) then
! Close statepoint file
call file_close(file_id)
! Write message
call write_message("Loading source file " // trim(path_source_point) &
// "...", 5)
! Open source file
file_id = file_open(path_source_point, 'r', parallel=.true.)
end if
! Read source
call read_source_bank(file_id)
end if
! Close file
call file_close(file_id)
end subroutine load_state_point
end subroutine load_state_point_f
end module state_point

294
src/state_point.cpp
View file

@ -1,6 +1,7 @@
#include "openmc/state_point.h"
#include <algorithm>
#include <cstdint> // for int64_t
#include <iomanip> // for setfill, setw
#include <string>
#include <vector>
@ -14,13 +15,292 @@
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/output.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/timer.h"
namespace openmc {
extern "C" void statepoint_write_f(hid_t file_id);
extern "C" void load_state_point_f(hid_t file_id);
extern "C" int
openmc_statepoint_write(const char* filename, bool* write_source)
{
// Set the filename
std::string filename_;
if (filename) {
filename_ = filename;
} else {
// Determine width for zero padding
int w = std::to_string(settings::n_max_batches).size();
// Set filename for state point
std::stringstream ss;
ss << settings::path_output << "statepoint." << std::setfill('0')
<< std::setw(w) << simulation::current_batch << ".h5";
filename_ = ss.str();
}
// Determine whether or not to write the source bank
bool write_source_ = write_source ? *write_source : true;
// Write message
write_message("Creating state point " + filename_ + "...", 5);
hid_t file_id;
if (mpi::master) {
// Create statepoint file
file_id = file_open(filename_, 'w');
// Write file type
write_attribute(file_id, "filetype", "statepoint");
// Write revision number for state point file
write_attribute(file_id, "version", VERSION_STATEPOINT);
// Write OpenMC version
write_attribute(file_id, "openmc_version", VERSION);
#ifdef GIT_SHA1
write_attribute(file_id, "git_sha1", GIT_SHA1);
#endif
// Write current date and time
write_attribute(file_id, "date_and_time", time_stamp());
// Write path to input
write_attribute(file_id, "path", settings::path_input);
// Write out random number seed
write_dataset(file_id, "seed", openmc_get_seed());
// Write run information
write_dataset(file_id, "energy_mode", settings::run_CE ?
"continuous-energy" : "multi-group");
switch (settings::run_mode) {
case RUN_MODE_FIXEDSOURCE:
write_dataset(file_id, "run_mode", "fixed source");
break;
case RUN_MODE_EIGENVALUE:
write_dataset(file_id, "run_mode", "eigenvalue");
break;
}
write_attribute(file_id, "photon_transport", settings::photon_transport);
write_dataset(file_id, "n_particles", settings::n_particles);
write_dataset(file_id, "n_batches", settings::n_batches);
// Write out current batch number
write_dataset(file_id, "current_batch", simulation::current_batch);
// Indicate whether source bank is stored in statepoint
write_attribute(file_id, "source_present", write_source_);
// Write out information for eigenvalue run
if (settings::run_mode == RUN_MODE_EIGENVALUE) write_eigenvalue_hdf5(file_id);
// Write tally-related data
hid_t tallies_group = create_group(file_id, "tallies");
meshes_to_hdf5(tallies_group);
statepoint_write_f(file_id);
close_group(tallies_group);
}
// Check for the no-tally-reduction method
if (!settings::reduce_tallies) {
// If using the no-tally-reduction method, we need to collect tally
// results before writing them to the state point file.
write_tally_results_nr(file_id);
} else if (mpi::master) {
// Write number of global realizations
write_dataset(file_id, "n_realizations", n_realizations);
// Write out the runtime metrics.
using namespace simulation;
hid_t runtime_group = create_group(file_id, "runtime");
write_dataset(runtime_group, "total initialization", time_initialize.elapsed());
write_dataset(runtime_group, "reading cross sections", time_read_xs.elapsed());
write_dataset(runtime_group, "simulation", time_inactive.elapsed()
+ time_active.elapsed());
write_dataset(runtime_group, "transport", time_transport.elapsed());
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
write_dataset(runtime_group, "inactive batches", time_inactive.elapsed());
}
write_dataset(runtime_group, "active batches", time_active.elapsed());
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
write_dataset(runtime_group, "synchronizing fission bank", time_bank.elapsed());
write_dataset(runtime_group, "sampling source sites", time_bank_sample.elapsed());
write_dataset(runtime_group, "SEND-RECV source sites", time_bank_sendrecv.elapsed());
}
write_dataset(runtime_group, "accumulating tallies", time_tallies.elapsed());
write_dataset(runtime_group, "total", time_total.elapsed());
close_group(runtime_group);
file_close(file_id);
}
#ifdef PHDF5
bool parallel = true;
#else
bool parallel = false;
#endif
// Write the source bank if desired
if (write_source_) {
if (mpi::master || parallel) file_id = file_open(filename_, 'a', true);
write_source_bank(file_id);
if (mpi::master || parallel) file_close(file_id);
}
return 0;
}
void restart_set_keff()
{
if (simulation::restart_batch > settings::n_inactive) {
for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) {
simulation::k_sum[0] += simulation::k_generation[i];
simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
}
int n = settings::gen_per_batch*n_realizations;
simulation::keff = simulation::k_sum[0] / n;
} else {
simulation::keff = simulation::k_generation.back();
}
}
void load_state_point()
{
// Write message
write_message("Loading state point " + settings::path_statepoint + "...", 5);
// Open file for reading
hid_t file_id = file_open(settings::path_statepoint.c_str(), 'r', true);
// Read filetype
std::string word;
read_attribute(file_id, "filetype", word);
if (word != "statepoint") {
fatal_error("OpenMC tried to restart from a non-statepoint file.");
}
// Read revision number for state point file and make sure it matches with
// current version
std::array<int, 2> array;
read_attribute(file_id, "version", array);
if (array != VERSION_STATEPOINT) {
fatal_error("State point version does not match current version in OpenMC.");
}
// Read and overwrite random number seed
int64_t seed;
read_dataset(file_id, "seed", seed);
openmc_set_seed(seed);
// It is not impossible for a state point to be generated from a CE run but
// to be loaded in to an MG run (or vice versa), check to prevent that.
read_dataset(file_id, "energy_mode", word);
if (word == "multi-group" && settings::run_CE) {
fatal_error("State point file is from multigroup run but current run is "
"continous energy.");
} else if (word == "continuous-energy" && !settings::run_CE) {
fatal_error("State point file is from continuous-energy run but current "
"run is multigroup!");
}
// Read and overwrite run information except number of batches
read_dataset(file_id, "run_mode", word);
if (word == "fixed source") {
settings::run_mode = RUN_MODE_FIXEDSOURCE;
} else if (word == "eigenvalue") {
settings::run_mode = RUN_MODE_EIGENVALUE;
}
read_attribute(file_id, "photon_transport", settings::photon_transport);
read_dataset(file_id, "n_particles", settings::n_particles);
int temp;
read_dataset(file_id, "n_batches", temp);
// Take maximum of statepoint n_batches and input n_batches
settings::n_batches = std::max(settings::n_batches, temp);
// Read batch number to restart at
read_dataset(file_id, "current_batch", simulation::restart_batch);
// Check for source in statepoint if needed
bool source_present;
read_attribute(file_id, "source_present", source_present);
if (simulation::restart_batch > settings::n_batches) {
fatal_error("The number batches specified in settings.xml is fewer "
" than the number of batches in the given statepoint file.");
}
// Read information specific to eigenvalue run
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
read_dataset(file_id, "n_inactive", temp);
read_eigenvalue_hdf5(file_id);
// Take maximum of statepoint n_inactive and input n_inactive
settings::n_inactive = std::max(settings::n_inactive, temp);
}
// Read number of realizations for global tallies
read_dataset(file_id, "n_realizations", n_realizations);
// Set k_sum, keff, and current_batch based on whether restart file is part
// of active cycle or inactive cycle
restart_set_keff();
simulation::current_batch = simulation::restart_batch;
// Check to make sure source bank is present
if (settings::path_sourcepoint == settings::path_statepoint &&
!source_present) {
fatal_error("Source bank must be contained in statepoint restart file");
}
// Read tallies to master. If we are using Parallel HDF5, all processes
// need to be included in the HDF5 calls.
#ifdef PHDF5
if (true) {
#else
if (mpi::master) {
#endif
// Read tally results
load_state_point_f(file_id);
}
// Read source if in eigenvalue mode
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
// Check if source was written out separately
if (!source_present) {
// Close statepoint file
file_close(file_id);
// Write message
write_message("Loading source file " + settings::path_sourcepoint
+ "...", 5);
// Open source file
file_id = file_open(settings::path_source.c_str(), 'r', true);
}
// Read source
read_source_bank(file_id);
}
// Close file
file_close(file_id);
}
hid_t h5banktype() {
// Create type for array of 3 reals
hsize_t dims[] {3};
@ -316,18 +596,4 @@ void write_tally_results_nr(hid_t file_id)
}
}
void restart_set_keff()
{
if (simulation::restart_batch > settings::n_inactive) {
for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) {
simulation::k_sum[0] += simulation::k_generation[i];
simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
}
int n = settings::gen_per_batch*n_realizations;
simulation::keff = simulation::k_sum[0] / n;
} else {
simulation::keff = simulation::k_generation.back();
}
}
} // namespace openmc

113
src/string.F90
View file

@ -2,7 +2,7 @@ module string
use, intrinsic :: ISO_C_BINDING
use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL
use constants, only: ERROR_INT, ERROR_REAL
use error, only: fatal_error, warning
use stl_vector, only: VectorInt
@ -14,81 +14,6 @@ module string
contains
!===============================================================================
! SPLIT_STRING takes a sentence and splits it into separate words much like the
! Python string.split() method.
!
! Arguments:
! string = input line
! words = array of words
! n = total number of words
!===============================================================================
subroutine split_string(string, words, n)
character(*), intent(in) :: string
character(*), intent(out) :: words(MAX_WORDS)
integer, intent(out) :: n
character(1) :: chr ! current character
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: i_end ! ending index of word
i_start = 0
i_end = 0
n = 0
do i = 1, len_trim(string)
chr = string(i:i)
! Note that ACHAR(9) is a horizontal tab
if ((i_start == 0) .and. (chr /= ' ') .and. (chr /= achar(9))) then
i_start = i
end if
if (i_start > 0) then
if ((chr == ' ') .or. (chr == achar(9))) i_end = i - 1
if (i == len_trim(string)) i_end = i
if (i_end > 0) then
n = n + 1
if (i_end - i_start + 1 > len(words(n))) then
call warning("The word '" // string(i_start:i_end) &
// "' is longer than the space allocated for it.")
end if
words(n) = string(i_start:i_end)
! reset indices
i_start = 0
i_end = 0
end if
end if
end do
end subroutine split_string
!===============================================================================
! CONCATENATE takes an array of words and concatenates them together in one
! string with a single space between words
!
! Arguments:
! words = array of words
! n_words = total number of words
! string = concatenated string
!===============================================================================
pure function concatenate(words, n_words) result(string)
integer, intent(in) :: n_words
character(*), intent(in) :: words(n_words)
character(MAX_LINE_LEN) :: string
integer :: i ! index
string = words(1)
if (n_words == 1) return
do i = 2, n_words
string = trim(string) // ' ' // words(i)
end do
end function concatenate
!===============================================================================
! TO_LOWER converts a string to all lower case characters
!===============================================================================
@ -247,25 +172,6 @@ contains
end function ends_with
!===============================================================================
! COUNT_DIGITS returns the number of digits needed to represent the input
! integer.
!===============================================================================
pure function count_digits(num) result(n_digits)
integer, intent(in) :: num
integer :: n_digits
n_digits = 1
do while (num / 10**(n_digits) /= 0 .and. abs(num / 10 **(n_digits-1)) /= 1&
&.and. n_digits /= 10)
! Note that 10 digits is the maximum needed to represent an integer(4) so
! the loop automatically exits when n_digits = 10.
n_digits = n_digits + 1
end do
end function count_digits
!===============================================================================
! INT4_TO_STR converts an integer(4) to a string.
!===============================================================================
@ -319,23 +225,6 @@ contains
end function str_to_int
!===============================================================================
! STR_TO_REAL converts an arbitrary string to a real(8)
!===============================================================================
pure function str_to_real(string) result(num)
character(*), intent(in) :: string
real(8) :: num
integer :: ioError
! Read string
read(UNIT=string, FMT=*, IOSTAT=ioError) num
if (ioError > 0) num = ERROR_REAL
end function str_to_real
!===============================================================================
! REAL_TO_STR converts a real(8) to a string based on how large the value is and
! how many significant digits are desired. By default, six significants digits

2
src/tallies/tally.F90
View file

@ -2,14 +2,12 @@ module tally
use, intrinsic :: ISO_C_BINDING
use algorithm, only: binary_search
use bank_header
use constants
use dict_header, only: EMPTY
use error, only: fatal_error
use geometry_header
use material_header
use math, only: t_percentile
use message_passing
use mgxs_interface
use nuclide_header

54
src/timer_header.F90
View file

@ -4,57 +4,6 @@ module timer_header
use constants, only: ZERO
implicit none
interface
function time_active_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_bank_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_bank_sample_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_bank_sendrecv_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_finalize_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_inactive_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_initialize_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_read_xs_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_tallies_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_total_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_transport_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
subroutine reset_timers() bind(C)
end subroutine
end interface
!===============================================================================
! TIMER represents a timer that can be started and stopped to measure how long
! different routines run. The intrinsic routine system_clock is used to measure
@ -73,9 +22,6 @@ module timer_header
procedure :: reset => timer_reset
end type Timer
! ============================================================================
! TIMING VARIABLES
contains
!===============================================================================