diff --git a/CMakeLists.txt b/CMakeLists.txt index 61a53b193..0e91ad9d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 38653d32d..66d718ff5 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -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); diff --git a/include/openmc/constants.h b/include/openmc/constants.h index b55ea061f..916f7693b 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -28,6 +28,7 @@ constexpr std::array VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELE constexpr int HDF5_VERSION[] {2, 0}; // Version numbers for binary files +constexpr std::array VERSION_STATEPOINT {17, 0}; constexpr std::array VERSION_PARTICLE_RESTART {2, 0}; constexpr std::array VERSION_TRACK {2, 0}; constexpr std::array VERSION_SUMMARY {6, 0}; diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index f7e93ec2b..d29f404b9 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -9,6 +9,7 @@ #include #include "xtensor/xtensor.hpp" +#include #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 diff --git a/include/openmc/error.h b/include/openmc/error.h index b02c0a9f8..ee7168a4e 100644 --- a/include/openmc/error.h +++ b/include/openmc/error.h @@ -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) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 6881cac77..398f03b09 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -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::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 void read_dataset(hid_t dset, std::vector& vec, bool indep=false) { diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index 87710fb6b..e07374df3 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -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. diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 7fcd58a06..1ca4d57b4 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -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 diff --git a/include/openmc/output.h b/include/openmc/output.h index 6c8226c65..e2387dd27 100644 --- a/include/openmc/output.h +++ b/include/openmc/output.h @@ -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 diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 79db8032d..1e2b85c15 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -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 diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 0c1db86f3..e44c57f1a 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -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 diff --git a/openmc/capi/math.py b/openmc/capi/math.py index e5a4b12c6..1baa9cf64 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -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 diff --git a/src/algorithm.F90 b/src/algorithm.F90 index 101b11b9a..d9f9e2e22 100644 --- a/src/algorithm.F90 +++ b/src/algorithm.F90 @@ -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. !=============================================================================== diff --git a/src/constants.F90 b/src/constants.F90 index 5e6b8629b..e68f5a6c3 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 1acef1262..22f9d392a 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -134,7 +134,7 @@ void load_dagmc_geometry() // initialize cell objects model::n_cells = model::DAG->num_entities(3); - + moab::EntityHandle graveyard = 0; for (int i = 0; i < model::n_cells; i++) { moab::EntityHandle vol_handle = model::DAG->entity_by_index(3, i+1); @@ -175,10 +175,22 @@ void load_dagmc_geometry() if (model::DAG->is_implicit_complement(vol_handle)) { if (model::DAG->has_prop(vol_handle, "mat")) { // if the implicit complement has been assigned a material, use it - std::string comp_mat = DMD.volume_material_property_data_eh[vol_handle]; - // Note: material numbers are set by UWUW - int mat_number = uwuw.material_library[comp_mat].metadata["mat_number"].asInt(); - c->material_.push_back(mat_number); + if (using_uwuw) { + std::string comp_mat = DMD.volume_material_property_data_eh[vol_handle]; + // Note: material numbers are set by UWUW + int mat_number = uwuw.material_library[comp_mat].metadata["mat_number"].asInt(); + c->material_.push_back(mat_number); + } else { + std::string mat_value; + rval= model::DAG->prop_value(vol_handle, "mat", mat_value); + MB_CHK_ERR_CONT(rval); + // remove _comp + std::string _comp = "_comp"; + size_t _comp_pos = mat_value.find(_comp); + if (_comp_pos != std::string::npos) { mat_value.erase(_comp_pos, _comp.length()); } + // assign IC material by id + c->material_.push_back(std::stoi(mat_value)); + } } else { // if no material is found, the implicit complement is void c->material_.push_back(MATERIAL_VOID); @@ -199,6 +211,11 @@ void load_dagmc_geometry() std::string cmp_str = mat_value; to_lower(cmp_str); + + if (cmp_str.find("graveyard") != std::string::npos) { + graveyard = vol_handle; + } + // material void checks if (cmp_str.find("void") != std::string::npos || cmp_str.find("vacuum") != std::string::npos || @@ -231,6 +248,11 @@ void load_dagmc_geometry() model::overlap_check_count.resize(model::cells.size(), 0); } + if (!graveyard) { + warning("No graveyard volume found in the DagMC model." + "This may result in lost particles and rapid simulation failure."); + } + /// Surfaces \\\ // initialize surface objects @@ -271,6 +293,17 @@ void load_dagmc_geometry() s->bc_ = BC_TRANSMIT; } + // graveyard check + moab::Range parent_vols; + rval = model::DAG->moab_instance()->get_parent_meshsets(surf_handle, parent_vols); + MB_CHK_ERR_CONT(rval); + + // if this surface belongs to the graveyard + if (graveyard && parent_vols.find(graveyard) != parent_vols.end()) { + // set BC to vacuum + s->bc_ = BC_VACUUM; + } + // add to global array and map model::surfaces[i] = s; model::surface_map[s->id_] = s->id_; diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index b3b8b850c..0f52649d8 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -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; diff --git a/src/endf_header.F90 b/src/endf_header.F90 deleted file mode 100644 index d45fbe68f..000000000 --- a/src/endf_header.F90 +++ /dev/null @@ -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 diff --git a/src/error.F90 b/src/error.F90 index ad07987eb..ca94c6ce7 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -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 diff --git a/src/error.cpp b/src/error.cpp index d17219f31..38001b72a 100644 --- a/src/error.cpp +++ b/src/error.cpp @@ -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 // for isatty +#endif + +#include // for exit +#include // for setw +#include 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 diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 1ac6ba423..8c18a4030 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 45ed2047e..b51d81b5b 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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) diff --git a/src/initialize.cpp b/src/initialize.cpp index 69986fdc0..e454917d0 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -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(); diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 40795200d..09a40660b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/math.F90 b/src/math.F90 index 5607a9556..81856c93c 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -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 diff --git a/src/math_functions.cpp b/src/math_functions.cpp index fe2c53230..919c305c4 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -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) { diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 057b7c113..ac02a9c71 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -893,6 +893,8 @@ void check_data_version(hid_t file_id) //============================================================================== extern "C" void extend_nuclides(); +extern "C" void nuclide_from_hdf5(hid_t group, const Nuclide* ptr, + const double* temps, int n, int n_nuclide); extern "C" int openmc_load_nuclide(const char* name) { @@ -917,6 +919,11 @@ extern "C" int openmc_load_nuclide(const char* name) int i_nuclide = data::nuclides.size(); data::nuclides.push_back(std::make_unique( group, temperature, i_nuclide)); + + // Read from Fortran too + nuclide_from_hdf5(group, data::nuclides.back().get(), + &temperature.front(), temperature.size(), i_nuclide + 1); + close_group(group); file_close(file_id); diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 4308bc487..f84716964 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 diff --git a/src/output.F90 b/src/output.F90 index aa66ae017..d026f7f87 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -20,7 +20,6 @@ module output use tally_header use tally_derivative_header use tally_filter - use timer_header implicit none @@ -29,12 +28,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 @@ -46,73 +39,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 !=============================================================================== @@ -167,299 +93,4 @@ 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 & - &" - 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 - -!=============================================================================== -! MEAN_STDEV computes the sample mean and standard deviation of the mean of a -! single tally score -!=============================================================================== - - pure function mean_stdev(result_, n) result(x) - real(8), intent(in) :: result_(2) ! sum and sum-of-squares - integer, intent(in) :: n ! number of realizations - real(8) :: x(2) ! mean and standard deviation - - x(1) = result_(1) / n - if (n > 1) then - x(2) = sqrt((result_(2) / n - x(1)*x(1))/(n - 1)) - else - x(2) = ZERO - end if - end function mean_stdev - end module output diff --git a/src/output.cpp b/src/output.cpp index efc4ba519..dbfb2ab7b 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -2,17 +2,22 @@ #include // for std::transform #include // for strlen -#include -#include // for setw, setprecision -#include // for left +#include // for time, localtime +#include // for setw, setprecision, put_time +#include // for fixed, scientific, left #include #include #include #include +#include // for pair + +#include +#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" @@ -23,16 +28,73 @@ #include "openmc/plot.h" #include "openmc/reaction.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/surface.h" #include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" #include "openmc/tallies/tally_scoring.h" +#include "openmc/timer.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'; +} + +//============================================================================== + std::string header(const char* msg) { // Determine how many times to repeat the '=' character. @@ -62,20 +124,16 @@ header(const char* msg, int level) { // Print header based on verbosity level. if (settings::verbosity >= level) - std::cout << out << "\n\n"; + std::cout << '\n' << out << "\n\n"; } //============================================================================== 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(); } @@ -252,19 +310,295 @@ print_overlap_check() //============================================================================== -std::pair -mean_stdev(double sum, double sum_sq, int n) +void print_usage() { - double mean, std_dev; - mean = sum / n; - if (n > 1) { - std_dev = std::sqrt((sum_sq / n - mean*mean) / (n - 1)); - } else { - std_dev = 0; + 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"; } - return {mean, std_dev}; } +//============================================================================== + +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 " + "\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 +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(>(K_COLLISION, 0), n); + std::cout << " k-effective (Collision) = " + << mean << " +/- " << t_n1 * stdev << '\n'; + std::tie(mean, stdev) = mean_stdev(>(K_TRACKLENGTH, 0), n); + std::cout << " k-effective (Track-length) = " + << mean << " +/- " << t_n1 * stdev << '\n'; + std::tie(mean, stdev) = mean_stdev(>(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(>(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); +} + +//============================================================================== + const std::unordered_map score_names = { {SCORE_FLUX, "Flux"}, {SCORE_TOTAL, "Total Reaction Rate"}, @@ -308,13 +642,13 @@ write_tallies() double t_value = 1; if (settings::confidence_intervals) { auto alpha = 1 - CONFIDENCE_LEVEL; - t_value = t_percentile_c(1 - alpha*0.5, n_realizations - 1); + t_value = t_percentile(1 - alpha*0.5, n_realizations - 1); } // Write header block. std::string tally_header("TALLY " + std::to_string(tally.id_)); if (!tally.name_.empty()) tally_header += ": " + tally.name_; - tallies_out << "\n" << header(tally_header) << "\n\n"; + tallies_out << header(tally_header) << "\n\n"; // Write derivative information. if (tally.deriv_ != C_NONE) { @@ -385,9 +719,7 @@ write_tallies() double mean, stdev; //TODO: off-by-one std::tie(mean, stdev) = mean_stdev( - results(filter_index-1, score_index, RESULT_SUM), - results(filter_index-1, score_index, RESULT_SUM_SQ), - n_realizations); + &results(filter_index-1, score_index, 0), n_realizations); tallies_out << std::string(indent+1, ' ') << std::left << std::setw(36) << score_name << " " << mean << " +/- " << t_value * stdev << "\n"; diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index d64717cf4..84dac0a1b 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -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 diff --git a/src/settings.F90 b/src/settings.F90 index 37cab546e..a9bbb5d31 100644 --- a/src/settings.F90 +++ b/src/settings.F90 @@ -86,14 +86,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 diff --git a/src/settings.cpp b/src/settings.cpp index 06410a424..47f5ff072 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -39,9 +39,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}; @@ -63,7 +63,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; diff --git a/src/simulation.cpp b/src/simulation.cpp index 89ad54c4f..7c7bbc4e6 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -30,12 +30,6 @@ namespace openmc { extern "C" void accumulate_tallies(); extern "C" void allocate_tally_results(); 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(); diff --git a/src/state_point.F90 b/src/state_point.F90 index 8ca3464c1..c96a395a6 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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,327 +20,186 @@ 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 - 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 - type(TallyDerivative), pointer :: deriv + type(TallyDerivative), pointer :: deriv - 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' - end if - - ! Determine whether or not to write the source bank - if (present(write_source)) then - write_source_ = write_source - else - write_source_ = .true. - 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 (n_tally_derivs() > 0) then - derivs_group = create_group(tallies_group, "derivatives") - do i = 0, n_tally_derivs() - 1 - !associate(deriv => tally_derivs(i)) - deriv => tally_deriv_c(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) - - ! Write number of 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) - - ! 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()))) - - ! 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") + ! Write information for derivatives. + if (n_tally_derivs() > 0) then + derivs_group = create_group(tallies_group, "derivatives") + do i = 0, n_tally_derivs() - 1 + deriv => tally_deriv_c(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 write_dataset(tally_group, "n_realizations", & - tally % n_realizations) + call close_group(deriv_group) + end do + call close_group(derivs_group) + end if - call write_dataset(tally_group, "n_filters", tally % n_filters()) - if (tally % n_filters() > 0) then - ! Write IDs of filters - allocate(id_array(tally % n_filters())) - do j = 1, tally % n_filters() - id_array(j) = filters(tally % filter(j) + 1) % obj % id - end do - call write_dataset(tally_group, "filters", id_array) - deallocate(id_array) - end if + ! Write number of filters + filters_group = create_group(tallies_group, "filters") + call write_attribute(filters_group, "n_filters", n_filters) - ! 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)+1) % name, '.') - if (i_xs > 0) then - str_array(j) = nuclides(tally % nuclide_bins(j)+1) % name(1 : i_xs-1) - else - str_array(j) = nuclides(tally % nuclide_bins(j)+1) % name - end if + 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) + + ! Write number of 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) + + ! 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()))) + + ! 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) + + call write_dataset(tally_group, "n_filters", tally % n_filters()) + if (tally % n_filters() > 0) then + ! Write IDs of filters + allocate(id_array(tally % n_filters())) + do j = 1, tally % n_filters() + id_array(j) = filters(tally % filter(j) + 1) % 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)+1) % name, '.') + if (i_xs > 0) then + str_array(j) = nuclides(tally % nuclide_bins(j)+1) % name(1 : i_xs-1) else - call get_name_c(tally % nuclide_bins(j)+1, 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)+1) % name end if else - str_array(j) = 'total' + call get_name_c(tally % nuclide_bins(j)+1, 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() /= C_NONE) then - deriv => tally_deriv_c(tally % deriv()) - call write_dataset(tally_group, "derivative", & - 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(tally % n_score_bins())) - do j = 1, tally % n_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() /= C_NONE) then + deriv => tally_deriv_c(tally % deriv()) + call write_dataset(tally_group, "derivative", deriv % id) + end if - deallocate(str_array) + ! Write scores. + call write_dataset(tally_group, "n_score_bins", tally % n_score_bins()) + allocate(str_array(tally % n_score_bins())) + do j = 1, tally % n_score_bins() + str_array(j) = reaction_name(tally % score_bins(j)) + end do + call write_dataset(tally_group, "score_bins", str_array) + deallocate(str_array) - call close_group(tally_group) - end associate - end do TALLY_METADATA - - 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) @@ -350,8 +208,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) @@ -362,241 +218,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 diff --git a/src/state_point.cpp b/src/state_point.cpp index aa1fc0e65..107ecdaaf 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -1,6 +1,7 @@ #include "openmc/state_point.h" #include +#include // for int64_t #include // for setfill, setw #include #include @@ -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 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}; @@ -318,18 +598,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 diff --git a/src/string.F90 b/src/string.F90 index 7e1ca41ba..f0419e2c1 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -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 diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 6d2c83ab2..96b2fff06 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -2,18 +2,15 @@ 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 - use output, only: header use particle_header, only: LocalCoord, Particle use settings use simulation_header diff --git a/src/timer_header.F90 b/src/timer_header.F90 index f1c89242e..3bbc20762 100644 --- a/src/timer_header.F90 +++ b/src/timer_header.F90 @@ -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 !=============================================================================== diff --git a/tests/regression_tests/dagmc/dagmc.h5m b/tests/regression_tests/dagmc/dagmc.h5m index 49c28c4f5..a55e061dc 100644 Binary files a/tests/regression_tests/dagmc/dagmc.h5m and b/tests/regression_tests/dagmc/dagmc.h5m differ diff --git a/tests/regression_tests/deplete/example_geometry.py b/tests/regression_tests/deplete/example_geometry.py index ca10c1f72..e54c087f6 100644 --- a/tests/regression_tests/deplete/example_geometry.py +++ b/tests/regression_tests/deplete/example_geometry.py @@ -47,10 +47,8 @@ def generate_initial_number_density(): fuel_dict['U238'] = 2.21371e22 fuel_dict['O16'] = 4.62954e22 fuel_dict['O17'] = 1.127684e20 - fuel_dict['I135'] = 1.0e10 fuel_dict['Xe135'] = 1.0e10 fuel_dict['Xe136'] = 1.0e10 - fuel_dict['Cs135'] = 1.0e10 fuel_dict['Gd156'] = 1.0e10 fuel_dict['Gd157'] = 1.0e10 # fuel_dict['O18'] = 9.51352e19 # Does not exist in ENDF71, merged into 17 @@ -62,10 +60,8 @@ def generate_initial_number_density(): fuel_gd_dict['Gd156'] = 3.95517E+10 fuel_gd_dict['Gd157'] = 1.08156e20 fuel_gd_dict['O16'] = 4.64035e22 - fuel_dict['I135'] = 1.0e10 fuel_dict['Xe136'] = 1.0e10 fuel_dict['Xe135'] = 1.0e10 - fuel_dict['Cs135'] = 1.0e10 # There are a whole bunch of 1e-10 stuff here. # Concentration to be used for cladding diff --git a/tests/regression_tests/deplete/test_reference.h5 b/tests/regression_tests/deplete/test_reference.h5 index 8750924b1..2f1ef375a 100644 Binary files a/tests/regression_tests/deplete/test_reference.h5 and b/tests/regression_tests/deplete/test_reference.h5 differ diff --git a/tests/regression_tests/uwuw/dagmc.h5m b/tests/regression_tests/uwuw/dagmc.h5m index 8fca1fff1..4fbba3bb6 100644 Binary files a/tests/regression_tests/uwuw/dagmc.h5m and b/tests/regression_tests/uwuw/dagmc.h5m differ diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 27faf01d5..89cf5d708 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -391,7 +391,11 @@ def test_restart(capi_init): def test_load_nuclide(capi_init): + # load multiple nuclides openmc.capi.load_nuclide('H3') + assert 'H3' in openmc.capi.nuclides openmc.capi.load_nuclide('Pu239') + assert 'Pu239' in openmc.capi.nuclides + # load non-existent nuclide with pytest.raises(exc.DataError): openmc.capi.load_nuclide('Pu3') diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index 994ee531a..5a066be6e 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -20,8 +20,7 @@ def test_get_atoms(res): t, n = res.get_atoms("1", "Xe135") t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0] - n_ref = [6.674732823364922e+08, 4.022377919003244e+14, - 3.469325382312551e+14, 3.639463751297369e+14] + n_ref = [6.67473282e+08, 3.76986925e+14, 3.68587383e+14, 3.91338675e+14] np.testing.assert_allclose(t, t_ref) np.testing.assert_allclose(n, n_ref) @@ -32,13 +31,11 @@ def test_get_reaction_rate(res): t, r = res.get_reaction_rate("1", "Xe135", "(n,gamma)") t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0] - n_ref = [6.674732823364922e+08, 4.022377919003244e+14, - 3.469325382312551e+14, 3.639463751297369e+14] - xs_ref = np.array([3.272630849911638e-05, 2.663333206501429e-05, - 3.378999816491878e-05, 3.277013860196171e-05]) + n_ref = [6.67473282e+08, 3.76986925e+14, 3.68587383e+14, 3.91338675e+14] + xs_ref = [3.32282266e-05, 2.76207120e-05, 4.10986677e-05, 3.72453665e-05] np.testing.assert_allclose(t, t_ref) - np.testing.assert_allclose(r, n_ref * xs_ref) + np.testing.assert_allclose(r, np.array(n_ref) * xs_ref) def test_get_eigenvalue(res): @@ -46,8 +43,7 @@ def test_get_eigenvalue(res): t, k = res.get_eigenvalue() t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0] - k_ref = [1.2162584685809341, 1.1021840248590875, 1.223040353281542, - 1.2360979537312164] + k_ref = [1.16984322, 1.19097427, 1.03012572, 1.20045627] np.testing.assert_allclose(t, t_ref) np.testing.assert_allclose(k, k_ref)