From b92677e70b6bec782dafd3ffebaa9708f756f0ad Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Jul 2018 06:43:35 -0500 Subject: [PATCH] Convert source distributions to C++ --- CMakeLists.txt | 6 +- include/openmc/capi.h | 3 +- include/openmc/distribution.h | 3 + include/openmc/distribution_multi.h | 2 + include/openmc/distribution_spatial.h | 8 +- include/openmc/file_utils.h | 17 ++ include/openmc/mgxs_interface.h | 1 + include/openmc/nuclide.h | 5 + include/openmc/particle.h | 2 +- include/openmc/settings.h | 3 + include/openmc/simulation.h | 11 + include/openmc/source.h | 54 ++++ src/api.F90 | 9 +- src/distribution_multivariate.F90 | 257 ----------------- src/distribution_univariate.F90 | 373 ------------------------- src/error.F90 | 3 +- src/initialize.cpp | 1 - src/input_xml.F90 | 59 ++-- src/material_header.F90 | 16 ++ src/mgxs_interface.F90 | 13 +- src/nuclide.cpp | 16 ++ src/settings.F90 | 7 +- src/settings.cpp | 22 ++ src/simulation.F90 | 19 +- src/simulation.cpp | 39 +++ src/simulation_header.F90 | 2 +- src/source.F90 | 141 ---------- src/source.cpp | 355 +++++++++++++++++++++++ src/source_header.F90 | 388 -------------------------- src/tallies/tally_header.F90 | 10 +- src/xml_interface.F90 | 11 + 31 files changed, 629 insertions(+), 1227 deletions(-) create mode 100644 include/openmc/file_utils.h create mode 100644 include/openmc/source.h delete mode 100644 src/distribution_multivariate.F90 delete mode 100644 src/distribution_univariate.F90 create mode 100644 src/nuclide.cpp delete mode 100644 src/source.F90 create mode 100644 src/source.cpp delete mode 100644 src/source_header.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 23a6176a8..fd3c58631 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -299,8 +299,6 @@ add_library(libopenmc SHARED src/cmfd_solver.F90 src/constants.F90 src/dict_header.F90 - src/distribution_multivariate.F90 - src/distribution_univariate.F90 src/eigenvalue.F90 src/endf.F90 src/endf_header.F90 @@ -341,8 +339,6 @@ add_library(libopenmc SHARED src/settings.F90 src/simulation_header.F90 src/simulation.F90 - src/source.F90 - src/source_header.F90 src/state_point.F90 src/stl_vector.F90 src/string.F90 @@ -400,6 +396,7 @@ add_library(libopenmc SHARED src/message_passing.cpp src/mgxs.cpp src/mgxs_interface.cpp + src/nuclide.cpp src/particle.cpp src/plot.cpp src/position.cpp @@ -414,6 +411,7 @@ add_library(libopenmc SHARED src/scattdata.cpp src/settings.cpp src/simulation.cpp + src/source.cpp src/state_point.cpp src/string_functions.cpp src/surface.cpp diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 8fd89334d..bc4371592 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -30,7 +30,6 @@ extern "C" { int openmc_extend_filters(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end); - int openmc_extend_sources(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_filter_get_id(int32_t index, int32_t* id); int openmc_filter_get_type(int32_t index, char* type); @@ -57,6 +56,7 @@ extern "C" { int openmc_material_add_nuclide(int32_t index, const char name[], double density); int openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n); int openmc_material_get_id(int32_t index, int32_t* id); + int openmc_material_get_fissionable(int32_t index, bool* fissionable); int openmc_material_get_volume(int32_t index, double* volume); int openmc_material_set_density(int32_t index, double density); int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density); @@ -84,7 +84,6 @@ extern "C" { int openmc_simulation_finalize(); int openmc_simulation_init(); int openmc_source_bank(struct Bank** ptr, int64_t* n); - int openmc_source_set_strength(int32_t index, double strength); int openmc_spatial_legendre_filter_get_order(int32_t index, int* order); int openmc_spatial_legendre_filter_get_params(int32_t index, int* axis, double* min, double* max); int openmc_spatial_legendre_filter_set_order(int32_t index, int order); diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index ed71a1ee2..b8ebcabf4 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -36,6 +36,9 @@ public: //! Sample a value from the distribution //! \return Sampled value double sample() const; + + const std::vector& x() const { return x_; } + const std::vector& p() const { return p_; } private: std::vector x_; //!< Possible outcomes std::vector p_; //!< Probability of each outcome diff --git a/include/openmc/distribution_multi.h b/include/openmc/distribution_multi.h index 9a5f9e93b..487f52831 100644 --- a/include/openmc/distribution_multi.h +++ b/include/openmc/distribution_multi.h @@ -73,6 +73,8 @@ public: Direction sample() const; }; +using UPtrAngle = std::unique_ptr; + } // namespace openmc #endif // DISTRIBUTION_MULTI_H diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 753496211..e94005f0f 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -1,4 +1,4 @@ -#ifndef OPENMC_DISTRIBTUION_SPATIAL_H +#ifndef OPENMC_DISTRIBUTION_SPATIAL_H #define OPENMC_DISTRIBUTION_SPATIAL_H #include "pugixml.hpp" @@ -48,10 +48,11 @@ public: //! Sample a position from the distribution //! \return Sampled position Position sample() const; + bool only_fissionable() const { return only_fissionable_; } private: Position lower_left_; //!< Lower-left coordinates of box Position upper_right_; //!< Upper-right coordinates of box - bool only_fissionable {false}; //!< Only accept sites in fissionable region? + bool only_fissionable_ {false}; //!< Only accept sites in fissionable region? }; //============================================================================== @@ -61,6 +62,7 @@ private: class SpatialPoint : public SpatialDistribution { public: SpatialPoint() : r_{} { }; + SpatialPoint(Position r) : r_{r} { }; explicit SpatialPoint(pugi::xml_node node); //! Sample a position from the distribution @@ -70,6 +72,8 @@ private: Position r_; //!< Single position at which sites are generated }; +using UPtrSpace = std::unique_ptr; + } // namespace openmc #endif // OPENMC_DISTRIBUTION_SPATIAL_H diff --git a/include/openmc/file_utils.h b/include/openmc/file_utils.h new file mode 100644 index 000000000..0a24cb547 --- /dev/null +++ b/include/openmc/file_utils.h @@ -0,0 +1,17 @@ +#ifndef OPENMC_FILE_UTILS_H +#define OPENMC_FILE_UTILS_H + +#include // for ifstream +#include + +namespace openmc { + +inline bool file_exists(const std::string& filename) +{ + std::ifstream s {filename}; + return s.good(); +} + +} + +#endif // OPENMC_FILE_UTILS_H \ No newline at end of file diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h index e5f56e997..02a43fece 100644 --- a/include/openmc/mgxs_interface.h +++ b/include/openmc/mgxs_interface.h @@ -16,6 +16,7 @@ namespace openmc { extern std::vector nuclides_MG; extern std::vector macro_xs; +extern "C" int num_energy_groups; //============================================================================== // Mgxs data loading interface methods diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 60769139f..8a19791ba 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -1,10 +1,15 @@ #ifndef OPENMC_NUCLIDE_H #define OPENMC_NUCLIDE_H +#include + #include "openmc/constants.h" namespace openmc { +extern std::array energy_min; +extern std::array energy_max; + //=============================================================================== //! Cached microscopic cross sections for a particular nuclide at the current //! energy diff --git a/include/openmc/particle.h b/include/openmc/particle.h index f351d1186..d186cdfd4 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -33,7 +33,7 @@ constexpr double REL_MAX_LOST_PARTICLES {1.0e-6}; //! Particle types enum class ParticleType { - neutron, photon, electron, positron + neutron = 1, photon = 2, electron = 3, positron = 4 }; extern "C" { diff --git a/include/openmc/settings.h b/include/openmc/settings.h index ea4dfebdc..49e0ab69e 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -18,8 +18,11 @@ namespace openmc { // Defined on Fortran side extern "C" bool openmc_check_overlaps; extern "C" bool openmc_particle_restart_run; +extern "C" bool openmc_photon_transport; extern "C" bool openmc_restart_run; +extern "C" bool openmc_run_CE; extern "C" bool openmc_write_all_tracks; +extern "C" bool openmc_write_initial_source; // Defined in .cpp // TODO: Make strings instead of char* once Fortran is gone diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 14bb56f27..a1350c936 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -2,12 +2,23 @@ #define OPENMC_SIMULATION_H #include +#include + +namespace openmc { extern "C" int openmc_current_batch; extern "C" int openmc_current_gen; extern "C" int64_t openmc_current_work; extern "C" int openmc_n_lost_particles; +extern "C" int openmc_total_gen; + +extern std::vector work_index; #pragma omp threadprivate(openmc_current_work) +extern "C" void openmc_simulation_init_c(); +void calculate_work(); + +} // namespace openmc + #endif // OPENMC_SIMULATION_H diff --git a/include/openmc/source.h b/include/openmc/source.h new file mode 100644 index 000000000..ac5f69a2f --- /dev/null +++ b/include/openmc/source.h @@ -0,0 +1,54 @@ +#ifndef OPENMC_SOURCE_H +#define OPENMC_SOURCE_H + +#include +#include + +#include "pugixml.hpp" + +#include "openmc/distribution_multi.h" +#include "openmc/distribution_spatial.h" +#include "openmc/capi.h" +#include "openmc/particle.h" + +namespace openmc { + +//============================================================================== +//! External source distribution +//============================================================================== + +class SourceDistribution { +public: + SourceDistribution(UPtrSpace space, UPtrAngle angle, UPtrDist energy); + explicit SourceDistribution(pugi::xml_node node); + + Bank sample() const; + double strength() const { return strength_; } +private: + ParticleType particle_ {ParticleType::neutron}; + double strength_ {1.0}; + UPtrSpace space_; + UPtrAngle angle_; + UPtrDist energy_; +}; + +//============================================================================== +// Global variables +//============================================================================== + +extern std::vector external_sources; + +//============================================================================== +// Functions +//============================================================================== + +//! Initialize source bank from file/distribution +extern "C" void initialize_source(); + +//! Sample a site from all external source distributions in proportion to their +//! source strength +extern "C" Bank sample_external_source(); + +} // namespace openmc + +#endif // OPENMC_SOURCE_H diff --git a/src/api.F90 b/src/api.F90 index 76be2f446..bfa1314c6 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -20,7 +20,6 @@ module openmc_api use random_lcg, only: openmc_get_seed, openmc_set_seed use settings use simulation_header - use source_header, only: openmc_extend_sources, openmc_source_set_strength use state_point, only: openmc_statepoint_write use tally_header use tally_filter_header @@ -44,7 +43,6 @@ module openmc_api public :: openmc_extend_filters public :: openmc_extend_cells public :: openmc_extend_materials - public :: openmc_extend_sources public :: openmc_extend_tallies public :: openmc_filter_get_id public :: openmc_filter_get_type @@ -83,7 +81,6 @@ module openmc_api public :: openmc_simulation_finalize public :: openmc_simulation_init public :: openmc_source_bank - public :: openmc_source_set_strength public :: openmc_tally_allocate public :: openmc_tally_get_estimator public :: openmc_tally_get_id @@ -301,12 +298,16 @@ contains use plot_header use sab_header use settings - use source_header use surface_header use tally_derivative_header use trigger_header use volume_header + interface + subroutine free_memory_source() bind(C) + end subroutine + end interface + call free_memory_geometry() call free_memory_surfaces() call free_memory_material() diff --git a/src/distribution_multivariate.F90 b/src/distribution_multivariate.F90 deleted file mode 100644 index 19db1c297..000000000 --- a/src/distribution_multivariate.F90 +++ /dev/null @@ -1,257 +0,0 @@ -module distribution_multivariate - - use constants, only: ONE, TWO, PI - use distribution_univariate - use error, only: fatal_error - use random_lcg, only: prn - use math, only: rotate_angle - use xml_interface - - implicit none - -!=============================================================================== -! UNITSPHEREDISTRIBUTION type defines a probability density function for points -! on the unit sphere. Extensions of this type are used to sample angular -! distributions for starting sources -!=============================================================================== - - type, abstract :: UnitSphereDistribution - real(8) :: reference_uvw(3) - contains - procedure(unitsphere_distribution_sample_), deferred :: sample - end type UnitSphereDistribution - - abstract interface - function unitsphere_distribution_sample_(this) result(uvw) - import UnitSphereDistribution - class(UnitSphereDistribution), intent(in) :: this - real(8) :: uvw(3) - end function unitsphere_distribution_sample_ - end interface - -!=============================================================================== -! Derived classes of UnitSphereDistribution -!=============================================================================== - - ! Explicit distribution of polar and azimuthal angles - type, extends(UnitSphereDistribution) :: PolarAzimuthal - class(Distribution), allocatable :: mu - class(Distribution), allocatable :: phi - contains - procedure :: sample => polar_azimuthal_sample - end type PolarAzimuthal - - ! Uniform distribution on the unit sphere - type, extends(UnitSphereDistribution) :: Isotropic - contains - procedure :: sample => isotropic_sample - end type Isotropic - - ! Monodirectional distribution - type, extends(UnitSphereDistribution) :: Monodirectional - contains - procedure :: sample => monodirectional_sample - end type Monodirectional - -!=============================================================================== -! SPATIALDISTRIBUTION type defines a probability density function for arbitrary -! points in Euclidean space. -!=============================================================================== - - type, abstract :: SpatialDistribution - contains - procedure(spatial_distribution_from_xml_), deferred :: from_xml - procedure(spatial_distribution_sample_), deferred :: sample - end type SpatialDistribution - - abstract interface - subroutine spatial_distribution_from_xml_(this, node) - import SpatialDistribution, XMLNode - class(SpatialDistribution), intent(inout) :: this - type(XMLNode), intent(in) :: node - end subroutine spatial_distribution_from_xml_ - - function spatial_distribution_sample_(this) result(xyz) - import SpatialDistribution - class(SpatialDistribution), intent(in) :: this - real(8) :: xyz(3) - end function spatial_distribution_sample_ - end interface - - type, extends(SpatialDistribution) :: CartesianIndependent - class(Distribution), allocatable :: x - class(Distribution), allocatable :: y - class(Distribution), allocatable :: z - contains - procedure :: from_xml => cartesian_independent_from_xml - procedure :: sample => cartesian_independent_sample - end type CartesianIndependent - - type, extends(SpatialDistribution) :: SpatialBox - real(8) :: lower_left(3) - real(8) :: upper_right(3) - logical :: only_fissionable = .false. - contains - procedure :: from_xml => spatial_box_from_xml - procedure :: sample => spatial_box_sample - end type SpatialBox - - type, extends(SpatialDistribution) :: SpatialPoint - real(8) :: xyz(3) - contains - procedure :: from_xml => spatial_point_from_xml - procedure :: sample => spatial_point_sample - end type SpatialPoint - -contains - - function polar_azimuthal_sample(this) result(uvw) - class(PolarAzimuthal), intent(in) :: this - real(8) :: uvw(3) - - real(8) :: mu ! cosine of polar angle - real(8) :: phi ! azimuthal angle - - ! Sample cosine of polar angle - mu = this % mu % sample() - if (mu == ONE) then - uvw(:) = this % reference_uvw - else - ! Sample azimuthal angle - phi = this % phi % sample() - uvw = rotate_angle(this % reference_uvw, mu, phi) - end if - end function polar_azimuthal_sample - - function isotropic_sample(this) result(uvw) - class(Isotropic), intent(in) :: this - real(8) :: uvw(3) - - real(8) :: phi - real(8) :: mu - - phi = TWO*PI*prn() - mu = TWO*prn() - ONE - uvw(1) = mu - uvw(2) = sqrt(ONE - mu*mu) * cos(phi) - uvw(3) = sqrt(ONE - mu*mu) * sin(phi) - end function isotropic_sample - - function monodirectional_sample(this) result(uvw) - class(Monodirectional), intent(in) :: this - real(8) :: uvw(3) - - uvw(:) = this % reference_uvw - end function monodirectional_sample - - subroutine cartesian_independent_from_xml(this, node) - class(CartesianIndependent), intent(inout) :: this - type(XMLNode), intent(in) :: node - - type(XMLNode) :: node_dist - - ! Read distribution for x coordinate - if (check_for_node(node, "x")) then - node_dist = node % child("x") - call distribution_from_xml(this % x, node_dist) - else - allocate(Discrete :: this % x) - select type (dist => this % x) - type is (Discrete) - allocate(dist % x(1), dist % p(1)) - dist % x(1) = ZERO - dist % p(1) = ONE - end select - end if - - ! Read distribution for y coordinate - if (check_for_node(node, "y")) then - node_dist = node % child("y") - call distribution_from_xml(this % y, node_dist) - else - allocate(Discrete :: this % y) - select type (dist => this % y) - type is (Discrete) - allocate(dist % x(1), dist % p(1)) - dist % x(1) = ZERO - dist % p(1) = ONE - end select - end if - - if (check_for_node(node, "z")) then - node_dist = node % child("z") - call distribution_from_xml(this % z, node_dist) - else - allocate(Discrete :: this % z) - select type (dist => this % z) - type is (Discrete) - allocate(dist % x(1), dist % p(1)) - dist % x(1) = ZERO - dist % p(1) = ONE - end select - end if - - end subroutine cartesian_independent_from_xml - - function cartesian_independent_sample(this) result(xyz) - class(CartesianIndependent), intent(in) :: this - real(8) :: xyz(3) - - xyz(1) = this % x % sample() - xyz(2) = this % y % sample() - xyz(3) = this % z % sample() - end function cartesian_independent_sample - - subroutine spatial_box_from_xml(this, node) - class(SpatialBox), intent(inout) :: this - type(XMLNode), intent(in) :: node - - real(8), allocatable :: temp_real(:) - - ! Make sure correct number of parameters are given - if (node_word_count(node, "parameters") /= 6) then - call fatal_error('Box/fission spatial source must have & - &six parameters specified.') - end if - - ! Read lower-right/upper-left coordinates - allocate(temp_real(6)) - call get_node_array(node, "parameters", temp_real) - this % lower_left(:) = temp_real(1:3) - this % upper_right(:) = temp_real(4:6) - deallocate(temp_real) - end subroutine spatial_box_from_xml - - function spatial_box_sample(this) result(xyz) - class(SpatialBox), intent(in) :: this - real(8) :: xyz(3) - - integer :: i - real(8) :: r(3) - - r = [ (prn(), i = 1,3) ] - xyz(:) = this % lower_left + r*(this % upper_right - this % lower_left) - end function spatial_box_sample - - subroutine spatial_point_from_xml(this, node) - class(SpatialPoint), intent(inout) :: this - type(XMLNode), intent(in) :: node - - ! Make sure correct number of parameters are given - if (node_word_count(node, "parameters") /= 3) then - call fatal_error('Point spatial source must have & - &three parameters specified.') - end if - - ! Read location of point source - call get_node_array(node, "parameters", this % xyz) - end subroutine spatial_point_from_xml - - function spatial_point_sample(this) result(xyz) - class(SpatialPoint), intent(in) :: this - real(8) :: xyz(3) - - xyz(:) = this % xyz - end function spatial_point_sample - -end module distribution_multivariate diff --git a/src/distribution_univariate.F90 b/src/distribution_univariate.F90 deleted file mode 100644 index 3842db95c..000000000 --- a/src/distribution_univariate.F90 +++ /dev/null @@ -1,373 +0,0 @@ -module distribution_univariate - - use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, & - MAX_LINE_LEN, MAX_WORD_LEN - use error, only: fatal_error - use random_lcg, only: prn - use math, only: maxwell_spectrum, watt_spectrum - use pugixml - use string, only: to_lower - use xml_interface - - implicit none - -!=============================================================================== -! DISTRIBUTION type defines a probability density function -!=============================================================================== - - type, abstract :: Distribution - contains - procedure(distribution_sample_), deferred :: sample - end type Distribution - - type DistributionContainer - class(Distribution), allocatable :: obj - end type DistributionContainer - - abstract interface - function distribution_sample_(this) result(x) - import Distribution - class(Distribution), intent(in) :: this - real(8) :: x - end function distribution_sample_ - end interface - -!=============================================================================== -! Derived classes of Distribution -!=============================================================================== - - ! Discrete distribution - type, extends(Distribution) :: Discrete - real(8), allocatable :: x(:) - real(8), allocatable :: p(:) - contains - procedure :: sample => discrete_sample - procedure :: initialize => discrete_initialize - end type Discrete - - ! Uniform distribution over the interval [a,b] - type, extends(Distribution) :: Uniform - real(8) :: a - real(8) :: b - contains - procedure :: sample => uniform_sample - end type Uniform - - ! Maxwellian distribution of form c*E*exp(-E/a) - type, extends(Distribution) :: Maxwell - real(8) :: theta - contains - procedure :: sample => maxwell_sample - end type Maxwell - - ! Watt fission spectrum with form c*exp(-E/a)*sinh(sqrt(b*E)) - type, extends(Distribution) :: Watt - real(8) :: a - real(8) :: b - contains - procedure :: sample => watt_sample - end type Watt - - ! Histogram or linear-linear interpolated tabular distribution - type, extends(Distribution) :: Tabular - integer :: interpolation - real(8), allocatable :: x(:) ! tabulated independent variable - real(8), allocatable :: p(:) ! tabulated probability density - real(8), allocatable :: c(:) ! cumulative distribution at tabulated values - contains - procedure :: sample => tabular_sample - procedure :: initialize => tabular_initialize - end type Tabular - - type, extends(Distribution) :: Equiprobable - real(8), allocatable :: x(:) - contains - procedure :: sample => equiprobable_sample - end type Equiprobable - -contains - - function discrete_sample(this) result(x) - class(Discrete), intent(in) :: this - real(8) :: x - - integer :: i ! loop counter - integer :: n ! size of distribution - real(8) :: c ! cumulative frequency - real(8) :: xi ! sampled CDF value - - n = size(this%x) - if (n > 1) then - xi = prn() - c = ZERO - do i = 1, size(this%x) - c = c + this%p(i) - if (xi < c) exit - end do - x = this%x(i) - else - x = this%x(1) - end if - end function discrete_sample - - subroutine discrete_initialize(this, x, p) - class(Discrete), intent(inout) :: this - real(8), intent(in) :: x(:) - real(8), intent(in) :: p(:) - - integer :: n - - ! Check length of x, p arrays - if (size(x) /= size(p)) then - call fatal_error('Tabulated probabilities not of same length as & - &independent variable.') - end if - - ! Copy probability density function - n = size(x) - allocate(this%x(n), this%p(n)) - this%x(:) = x(:) - this%p(:) = p(:) - - ! Normalize density function - this%p(:) = this%p(:)/sum(this%p) - end subroutine - - function uniform_sample(this) result(x) - class(Uniform), intent(in) :: this - real(8) :: x - - x = this%a + prn()*(this%b - this%a) - end function uniform_sample - - function maxwell_sample(this) result(x) - class(Maxwell), intent(in) :: this - real(8) :: x - - x = maxwell_spectrum(this%theta) - end function maxwell_sample - - function watt_sample(this) result(x) - class(Watt), intent(in) :: this - real(8) :: x - - x = watt_spectrum(this%a, this%b) - end function watt_sample - - function tabular_sample(this) result(x) - class(Tabular), intent(in) :: this - real(8) :: x - - integer :: i - real(8) :: c ! sampled cumulative frequency - real(8) :: m ! slope of PDF - real(8) :: x_i, x_i1 ! i-th and (i+1)th x values - real(8) :: c_i, c_i1 ! i-th and (i+1)th cumulative distribution values - real(8) :: p_i, p_i1 ! i-th and (i+1)th probability density values - - ! Sample value of CDF - c = prn() - - ! Find first CDF bin which is above the sampled value - c_i = this%c(1) - do i = 1, size(this%c) - 1 - c_i1 = this%c(i + 1) - if (c <= c_i1) exit - c_i = c_i1 - end do - - ! Determine bounding PDF values - x_i = this%x(i) - p_i = this%p(i) - - if (this%interpolation == HISTOGRAM) then - ! Histogram interpolation - if (p_i > ZERO) then - x = x_i + (c - c_i)/p_i - else - x = x_i - end if - else - ! Linear-linear interpolation - x_i1 = this%x(i + 1) - p_i1 = this%p(i + 1) - - m = (p_i1 - p_i)/(x_i1 - x_i) - if (m == ZERO) then - x = x_i + (c - c_i)/p_i - else - x = x_i + (sqrt(max(ZERO, p_i*p_i + 2*m*(c - c_i))) - p_i)/m - end if - end if - end function tabular_sample - - subroutine tabular_initialize(this, x, p, interp) - class(Tabular), intent(inout) :: this - real(8), intent(in) :: x(:) - real(8), intent(in) :: p(:) - integer, intent(in) :: interp - - integer :: i - integer :: n - - ! Check interpolation parameter - if (interp /= HISTOGRAM .and. interp /= LINEAR_LINEAR) then - call fatal_error('Only histogram and linear-linear interpolation for tabular & - &distribution is supported.') - end if - - ! Check length of x, p arrays - if (size(x) /= size(p)) then - call fatal_error('Tabulated probabilities not of same length as & - &independent variable.') - end if - - ! Copy probability density function and interpolation parameter - n = size(x) - allocate(this%x(n), this%p(n), this%c(n)) - this%interpolation = interp - this%x(:) = x(:) - this%p(:) = p(:) - - ! Calculate cumulative distribution function - this%c(1) = ZERO - do i = 2, n - if (this%interpolation == HISTOGRAM) then - this%c(i) = this%c(i-1) + this%p(i-1)*(this%x(i) - this%x(i-1)) - elseif (this%interpolation == LINEAR_LINEAR) then - this%c(i) = this%c(i-1) + HALF*(this%p(i-1) + this%p(i)) * & - (this%x(i) - this%x(i-1)) - end if - end do - - ! Normalize density and distribution functions - this%p(:) = this%p(:)/this%c(n) - this%c(:) = this%c(:)/this%c(n) - end subroutine tabular_initialize - - function equiprobable_sample(this) result(x) - class(Equiprobable), intent(in) :: this - real(8) :: x - - integer :: i - integer :: n - real(8) :: r - real(8) :: xl, xr - - n = size(this%x) - - r = prn() - i = 1 + int((n - 1)*r) - - xl = this%x(i) - xr = this%x(i+1) - x = xl + ((n - 1)*r - i + ONE) * (xr - xl) - end function equiprobable_sample - - subroutine distribution_from_xml(dist, node_dist) - class(Distribution), allocatable, intent(inout) :: dist - type(XMLNode), intent(in) :: node_dist - - character(MAX_WORD_LEN) :: type - character(MAX_LINE_LEN) :: temp_str - integer :: n - integer :: temp_int - real(8), allocatable :: temp_real(:) - - if (check_for_node(node_dist, "type")) then - ! Determine type of distribution - call get_node_value(node_dist, "type", type) - - ! Determine number of parameters specified - if (check_for_node(node_dist, "parameters")) then - n = node_word_count(node_dist, "parameters") - else - n = 0 - end if - - ! Allocate extension of Distribution - select case (to_lower(type)) - case ('uniform') - allocate(Uniform :: dist) - if (n /= 2) then - call fatal_error('Uniform distribution must have two & - ¶meters specified.') - end if - - case ('maxwell') - allocate(Maxwell :: dist) - if (n /= 1) then - call fatal_error('Maxwell energy distribution must have one & - ¶meter specified.') - end if - - case ('watt') - allocate(Watt :: dist) - if (n /= 2) then - call fatal_error('Watt energy distribution must have two & - ¶meters specified.') - end if - - case ('discrete') - allocate(Discrete :: dist) - - case ('tabular') - allocate(Tabular :: dist) - - case default - call fatal_error('Invalid distribution type: ' // trim(type) // '.') - - end select - - ! Read parameters and interpolation for distribution - select type (dist) - type is (Uniform) - allocate(temp_real(2)) - call get_node_array(node_dist, "parameters", temp_real) - dist%a = temp_real(1) - dist%b = temp_real(2) - deallocate(temp_real) - - type is (Maxwell) - call get_node_value(node_dist, "parameters", dist%theta) - - type is (Watt) - allocate(temp_real(2)) - call get_node_array(node_dist, "parameters", temp_real) - dist%a = temp_real(1) - dist%b = temp_real(2) - deallocate(temp_real) - - type is (Discrete) - allocate(temp_real(n)) - call get_node_array(node_dist, "parameters", temp_real) - call dist%initialize(temp_real(1:n/2), temp_real(n/2+1:n)) - deallocate(temp_real) - - type is (Tabular) - ! Read interpolation - if (check_for_node(node_dist, "interpolation")) then - call get_node_value(node_dist, "interpolation", temp_str) - select case(to_lower(temp_str)) - case ('histogram') - temp_int = HISTOGRAM - case ('linear-linear') - temp_int = LINEAR_LINEAR - case default - call fatal_error("Unknown interpolation type for distribution: " & - // trim(temp_str)) - end select - else - temp_int = HISTOGRAM - end if - - ! Read and initialize tabular distribution - allocate(temp_real(n)) - call get_node_array(node_dist, "parameters", temp_real) - call dist%initialize(temp_real(1:n/2), temp_real(n/2+1:n), temp_int) - deallocate(temp_real) - end select - end if - end subroutine distribution_from_xml - -end module distribution_univariate diff --git a/src/error.F90 b/src/error.F90 index 82b142b51..c00006c8d 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -270,8 +270,9 @@ contains 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, level) + call write_message(message_out(2:), level) end subroutine write_message_from_c end module error diff --git a/src/initialize.cpp b/src/initialize.cpp index 59e13512e..21ab93b72 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -2,7 +2,6 @@ #include #include -#include #include #include diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a2b1a1fd9..1d7979041 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -7,8 +7,6 @@ module input_xml use cmfd_header, only: cmfd_mesh use constants use dict_header, only: DictIntInt, DictCharInt, DictEntryCI - use distribution_multivariate - use distribution_univariate use endf, only: reaction_name use error, only: fatal_error, warning, write_message, openmc_err_msg use geometry, only: neighbor_lists @@ -28,7 +26,6 @@ module input_xml use surface_header use set_header, only: SetChar use settings - use source_header use stl_vector, only: VectorInt, VectorReal, VectorChar use string, only: to_lower, to_str, str_to_int, str_to_real, & starts_with, ends_with, tokenize, split_string, & @@ -102,6 +99,13 @@ module input_xml integer(C_INT32_T), intent(in), value :: univ integer(C_INT) :: n end function maximum_levels + + subroutine set_particle_energy_bounds(particle, E_min, E_max) bind(C) + import C_INT, C_DOUBLE + integer(C_INT), value :: particle + real(C_DOUBLE), value :: E_min + real(C_DOUBLE), value :: E_max + end subroutine end interface contains @@ -212,7 +216,6 @@ contains integer, allocatable :: temp_int_array(:) integer :: n_tracks logical :: file_exists - character(MAX_WORD_LEN) :: type character(MAX_LINE_LEN) :: filename type(XMLDocument) :: doc type(XMLNode) :: root @@ -227,7 +230,6 @@ contains type(XMLNode) :: node_vol type(XMLNode) :: node_tab_leg type(XMLNode), allocatable :: node_mesh_list(:) - type(XMLNode), allocatable :: node_source_list(:) type(XMLNode), allocatable :: node_vol_list(:) ! Check if settings.xml exists @@ -456,46 +458,13 @@ contains ! ========================================================================== ! EXTERNAL SOURCE - ! Get point to list of elements and make sure there is at least one - call get_node_list(root, "source", node_source_list) - n = size(node_source_list) - - if (n == 0) then - ! Default source is isotropic point source at origin with Watt spectrum - allocate(external_source(1)) - external_source % particle = NEUTRON - external_source % strength = ONE - - allocate(SpatialPoint :: external_source(1) % space) - select type (space => external_source(1) % space) - type is (SpatialPoint) - space % xyz(:) = [ZERO, ZERO, ZERO] - end select - - allocate(Isotropic :: external_source(1) % angle) - external_source(1) % angle % reference_uvw(:) = [ZERO, ZERO, ONE] - - allocate(Watt :: external_source(1) % energy) - select type(energy => external_source(1) % energy) - type is (Watt) - energy % a = 0.988e6_8 - energy % b = 2.249e-6_8 - end select - else - ! Allocate array for sources - allocate(external_source(n)) - end if + ! Handled on C++ side ! Check if we want to write out source if (check_for_node(root, "write_initial_source")) then call get_node_value(root, "write_initial_source", write_initial_source) end if - ! Read each source - do i = 1, n - call external_source(i) % from_xml(node_source_list(i), path_source) - end do - ! Survival biasing if (check_for_node(root, "survival_biasing")) then call get_node_value(root, "survival_biasing", survival_biasing) @@ -1004,7 +973,7 @@ contains subroutine read_geometry_xml() integer :: i, j, k - integer :: n, n_mats, n_rlats, n_hlats + integer :: n, n_rlats, n_hlats integer :: id integer :: univ_id integer :: n_cells_in_univ @@ -1012,7 +981,6 @@ contains logical :: file_exists logical :: boundary_exists character(MAX_LINE_LEN) :: filename - character(MAX_WORD_LEN), allocatable :: sarray(:) character(:), allocatable :: region_spec type(Cell), pointer :: c class(Lattice), pointer :: lat @@ -1473,7 +1441,6 @@ contains character(3) :: element ! name of element, e.g. Zr character(MAX_WORD_LEN) :: units ! units on density character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml - character(MAX_LINE_LEN) :: temp_str ! temporary string when reading character(MAX_WORD_LEN), allocatable :: sarray(:) real(8) :: val ! value entered for density real(8) :: temp_dble ! temporary double prec. real @@ -3479,6 +3446,8 @@ contains ! Get the minimum and maximum energies energy_min(NEUTRON) = energy_bins(num_energy_groups + 1) energy_max(NEUTRON) = energy_bins(1) + call set_particle_energy_bounds(NEUTRON, energy_min(NEUTRON), & + energy_max(NEUTRON)) ! Get the datasets present in the library call get_groups(file_id, names) @@ -3639,6 +3608,8 @@ contains nuclides(i_nuclide) % grid(1) % energy(1)) energy_max(NEUTRON) = min(energy_max(NEUTRON), nuclides(i_nuclide) % & grid(1) % energy(size(nuclides(i_nuclide) % grid(1) % energy))) + call set_particle_energy_bounds(NEUTRON, energy_min(NEUTRON), & + energy_max(NEUTRON)) end if ! Add name and alias to dictionary @@ -3672,6 +3643,8 @@ contains energy_max(PHOTON) = min(energy_max(PHOTON), & exp(elements(i_element) % energy(size(elements(i_element) & % energy)))) + call set_particle_energy_bounds(PHOTON, energy_min(PHOTON), & + energy_max(PHOTON)) end if ! Add element to set @@ -3713,6 +3686,8 @@ contains if (size(ttb_e_grid) >= 1) then energy_min(PHOTON) = max(energy_min(PHOTON), ttb_e_grid(1)) energy_max(PHOTON) = min(energy_max(PHOTON), ttb_e_grid(size(ttb_e_grid))) + call set_particle_energy_bounds(PHOTON, energy_min(PHOTON), & + energy_max(PHOTON)) end if ! Take logarithm of energies since they are log-log interpolated diff --git a/src/material_header.F90 b/src/material_header.F90 index 1f48168af..71e6ec5f1 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -678,6 +678,22 @@ contains end function openmc_material_get_id + function openmc_material_get_fissionable(index, fissionable) result(err) bind(C) + ! returns whether a material is fissionable + integer(C_INT32_T), value :: index + logical(C_BOOL), intent(out) :: fissionable + integer(C_INT) :: err + + if (index >= 1 .and. index <= size(materials)) then + fissionable = materials(index) % fissionable + err = 0 + else + err = E_OUT_OF_BOUNDS + call set_errmsg("Index in materials array is out of bounds.") + end if + end function openmc_material_get_fissionable + + function openmc_material_set_id(index, id) result(err) bind(C) ! Set the ID of a material integer(C_INT32_T), value, intent(in) :: index diff --git a/src/mgxs_interface.F90 b/src/mgxs_interface.F90 index 2a579e930..9dccc9230 100644 --- a/src/mgxs_interface.F90 +++ b/src/mgxs_interface.F90 @@ -147,7 +147,7 @@ module mgxs_interface end interface ! Number of energy groups - integer(C_INT) :: num_energy_groups + integer(C_INT), bind(C) :: num_energy_groups ! Number of delayed groups integer(C_INT) :: num_delayed_groups @@ -159,6 +159,13 @@ module mgxs_interface real(8), allocatable :: energy_bin_avg(:) ! Energy group structure with increasing energy - real(8), allocatable :: rev_energy_bins(:) + real(C_DOUBLE), allocatable, target :: rev_energy_bins(:) -end module mgxs_interface \ No newline at end of file +contains + + function rev_energy_bins_ptr() result(ptr) bind(C) + type(C_PTR) :: ptr + ptr = C_LOC(rev_energy_bins(1)) + end function + +end module mgxs_interface diff --git a/src/nuclide.cpp b/src/nuclide.cpp new file mode 100644 index 000000000..f7f5399ba --- /dev/null +++ b/src/nuclide.cpp @@ -0,0 +1,16 @@ +#include "openmc/nuclide.h" + +namespace openmc { + +// Order corresponds to ParticleType enum +std::array energy_min {0.0, 0.0}; +std::array energy_max {INFTY, INFTY}; + +extern "C" void +set_particle_energy_bounds(int particle, double E_min, double E_max) +{ + energy_min[particle - 1] = E_min; + energy_max[particle - 1] = E_max; +} + +} // namespace openmc diff --git a/src/settings.F90 b/src/settings.F90 index 972eba1d5..dfe4afbe3 100644 --- a/src/settings.F90 +++ b/src/settings.F90 @@ -9,7 +9,7 @@ module settings ! ============================================================================ ! ENERGY TREATMENT RELATED VARIABLES - logical(C_BOOL) :: run_CE = .true. ! Run in CE mode? + logical(C_BOOL), bind(C, name='openmc_run_CE') :: run_CE = .true. ! Run in CE mode? ! ============================================================================ ! CONTINUOUS-ENERGY CROSS SECTION RELATED VARIABLES @@ -26,7 +26,7 @@ module settings integer :: n_log_bins ! number of bins for logarithmic grid - logical :: photon_transport = .false. + logical(C_BOOL), bind(C, name='openmc_photon_transport') :: photon_transport = .false. integer :: electron_treatment = ELECTRON_TTB ! ============================================================================ @@ -104,7 +104,8 @@ module settings particle_restart_run = .false. ! Write out initial source - logical :: write_initial_source = .false. + logical(C_BOOL), bind(C, name='openmc_write_initial_source') :: & + write_initial_source = .false. ! Whether create fission neutrons or not. Only applied for MODE_FIXEDSOURCE logical :: create_fission_neutrons = .true. diff --git a/src/settings.cpp b/src/settings.cpp index b23013c31..aafa0fb23 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -2,7 +2,11 @@ #include "openmc/capi.h" #include "openmc/constants.h" +#include "openmc/distribution.h" +#include "openmc/distribution_multi.h" +#include "openmc/distribution_spatial.h" #include "openmc/error.h" +#include "openmc/source.h" #include "openmc/string_utils.h" #include "openmc/xml_interface.h" @@ -98,6 +102,24 @@ void read_settings(pugi::xml_node* root) temperature_range[0] = range[0]; temperature_range[1] = range[1]; } + + // ========================================================================== + // EXTERNAL SOURCE + + // Get point to list of elements and make sure there is at least one + for (pugi::xml_node node : root->children("source")) { + external_sources.emplace_back(node); + } + + // If no source specified, default to isotropic point source at origin with Watt spectrum + if (external_sources.empty()) { + SourceDistribution source { + UPtrSpace{new SpatialPoint({0.0, 0.0, 0.0})}, + UPtrAngle{new Isotropic()}, + UPtrDist{new Watt(0.988, 2.249e-6)} + }; + external_sources.push_back(std::move(source)); + } } } // namespace openmc diff --git a/src/simulation.F90 b/src/simulation.F90 index df6bf13ba..b5d24180c 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -30,7 +30,6 @@ module simulation use random_lcg, only: set_particle_seed use settings use simulation_header - use source, only: initialize_source, sample_external_source use state_point, only: openmc_statepoint_write, write_source_point, load_state_point use string, only: to_str use tally, only: accumulate_tallies, setup_active_tallies, & @@ -52,6 +51,19 @@ module simulation integer(C_INT), parameter :: STATUS_EXIT_MAX_BATCH = 1 integer(C_INT), parameter :: STATUS_EXIT_ON_TRIGGER = 2 + interface + subroutine openmc_simulation_init_c() bind(C) + end subroutine + + subroutine initialize_source() bind(C) + end subroutine + + function sample_external_source() result(site) bind(C) + import Bank + type(Bank) :: site + end function + end interface + contains !=============================================================================== @@ -298,7 +310,7 @@ contains do i = 1, work call set_particle_seed((total_gen + overall_generation()) * & n_particles + work_index(rank) + i) - call sample_external_source(source_bank(i)) + source_bank(i) = sample_external_source() end do end if end if @@ -425,6 +437,9 @@ contains ! Skip if simulation has already been initialized if (simulation_initialized) return + ! Call initialization on C++ side + call openmc_simulation_init_c() + ! Set up tally procedure pointers call init_tally_routines() diff --git a/src/simulation.cpp b/src/simulation.cpp index e1dd0ac93..19cc78d19 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,4 +1,7 @@ +#include "openmc/simulation.h" + #include "openmc/capi.h" +#include "openmc/message_passing.h" // OPENMC_RUN encompasses all the main logic where iterations are performed // over the batches, generations, and histories in a fixed source or k-eigenvalue @@ -16,3 +19,39 @@ int openmc_run() { openmc_simulation_finalize(); return err; } + +namespace openmc { + +std::vector work_index; + +void calculate_work() +{ + // Determine minimum amount of particles to simulate on each processor + int64_t min_work = n_particles/mpi::n_procs; + + // Determine number of processors that have one extra particle + int64_t remainder = n_particles % mpi::n_procs; + + int64_t i_bank = 0; + work_index.reserve(mpi::n_procs); + work_index.push_back(0); + for (int i = 0; i < mpi::n_procs; ++i) { + // Number of particles for rank i + int64_t work_i = i < remainder ? min_work + 1 : min_work; + + // Set number of particles + if (mpi::rank == i) openmc_work = work_i; + + // Set index into source bank for rank i + i_bank += work_i; + work_index.push_back(i_bank); + } +} + +void openmc_simulation_init_c() +{ + // Determine how much work each process should do + calculate_work(); +} + +} // namespace openmc diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 60bed6425..7addf9b2c 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -22,7 +22,7 @@ module simulation_header integer(C_INT), bind(C, name='openmc_current_batch') :: current_batch ! current batch integer(C_INT), bind(C, name='openmc_current_gen') :: current_gen ! current generation within a batch - integer :: total_gen = 0 ! total number of generations simulated + integer(C_INT), bind(C, name='openmc_total_gen') :: total_gen = 0 ! total number of generations simulated logical(C_BOOL), bind(C, name='openmc_simulation_initialized') :: & simulation_initialized = .false. logical :: need_depletion_rx ! need to calculate depletion reaction rx? diff --git a/src/source.F90 b/src/source.F90 deleted file mode 100644 index a6b6a92ac..000000000 --- a/src/source.F90 +++ /dev/null @@ -1,141 +0,0 @@ -module source - -#ifdef OPENMC_MPI - use message_passing -#endif - - use algorithm, only: binary_search - use bank_header, only: Bank, source_bank - use constants - use distribution_univariate, only: Discrete - use distribution_multivariate, only: SpatialBox - use error, only: fatal_error - use geometry, only: find_cell - use hdf5_interface - use math - use message_passing, only: rank - use mgxs_interface, only: rev_energy_bins, num_energy_groups - use output, only: write_message - use particle_header, only: Particle - use random_lcg, only: prn, set_particle_seed, prn_set_stream - use settings - use simulation_header - use source_header, only: external_source - use string, only: to_str - use state_point, only: read_source_bank, write_source_bank - - implicit none - -contains - -!=============================================================================== -! INITIALIZE_SOURCE initializes particles in the source bank -!=============================================================================== - - subroutine initialize_source() - - integer(8) :: i ! loop index over bank sites - integer(8) :: id ! particle id - integer(HID_T) :: file_id - character(MAX_WORD_LEN) :: filetype - character(MAX_FILE_LEN) :: filename - type(Bank), pointer :: src ! source bank site - - call write_message("Initializing source particles...", 5) - - if (path_source /= '') then - ! Read the source from a binary file instead of sampling from some - ! assumed source distribution - - call write_message('Reading source file from ' // trim(path_source) & - &// '...', 6) - - ! Open the binary file - file_id = file_open(path_source, 'r', parallel=.true.) - - ! Read the file type - call read_attribute(filetype, file_id, "filetype") - - ! Check to make sure this is a source file - if (filetype /= 'source' .and. filetype /= 'statepoint') then - call fatal_error("Specified starting source file not a source file & - &type.") - end if - - ! Read in the source bank - call read_source_bank(file_id, work_index, source_bank) - - ! Close file - call file_close(file_id) - - else - ! Generation source sites from specified distribution in user input - do i = 1, work - ! Get pointer to source bank site - src => source_bank(i) - - ! initialize random number seed - id = total_gen*n_particles + work_index(rank) + i - call set_particle_seed(id) - - ! sample external source distribution - call sample_external_source(src) - end do - end if - - ! Write out initial source - if (write_initial_source) then - call write_message('Writing out initial source...', 5) - filename = trim(path_output) // 'initial_source.h5' - file_id = file_open(filename, 'w', parallel=.true.) - call write_source_bank(file_id, work_index, source_bank) - call file_close(file_id) - end if - - end subroutine initialize_source - -!=============================================================================== -! SAMPLE_EXTERNAL_SOURCE samples the user-specified external source and stores -! the position, angle, and energy in a Bank type. -!=============================================================================== - - subroutine sample_external_source(site) - type(Bank), intent(inout) :: site ! source site - - integer :: i ! dummy loop index - integer :: n_source ! number of source distributions - real(8) :: c ! cumulative frequency - real(8) :: xi - - ! Set the random number generator to the source stream. - call prn_set_stream(STREAM_SOURCE) - - ! Sample from among multiple source distributions - n_source = size(external_source) - if (n_source > 1) then - xi = prn()*sum(external_source(:) % strength) - c = ZERO - do i = 1, n_source - c = c + external_source(i) % strength - if (xi < c) exit - end do - else - i = 1 - end if - - ! Sample source site from i-th source distribution - site = external_source(i) % sample() - - ! If running in MG, convert site % E to group - if (.not. run_CE) then - site % E = real(binary_search(rev_energy_bins, num_energy_groups + 1, & - site % E), 8) - site % E = num_energy_groups + 1 - site % E - end if - - ! Set the random number generator back to the tracking stream. - call prn_set_stream(STREAM_TRACKING) - - end subroutine sample_external_source - -end module source diff --git a/src/source.cpp b/src/source.cpp new file mode 100644 index 000000000..e01c9795e --- /dev/null +++ b/src/source.cpp @@ -0,0 +1,355 @@ +#include "openmc/source.h" + +#include // for move +#include // for stringstream + +#include "xtensor/xadapt.hpp" + +#include "openmc/cell.h" +#include "openmc/error.h" +#include "openmc/file_utils.h" +#include "openmc/hdf5_interface.h" +#include "openmc/material.h" +#include "openmc/message_passing.h" +#include "openmc/mgxs_interface.h" +#include "openmc/nuclide.h" +#include "openmc/capi.h" +#include "openmc/random_lcg.h" +#include "openmc/search.h" +#include "openmc/settings.h" +#include "openmc/simulation.h" +#include "openmc/state_point.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +std::vector external_sources; + +//============================================================================== +// SourceDistribution implementation +//============================================================================== + +SourceDistribution::SourceDistribution(UPtrSpace space, UPtrAngle angle, UPtrDist energy) + : space_{std::move(space)}, angle_{std::move(angle)}, energy_{std::move(energy)} { } + +SourceDistribution::SourceDistribution(pugi::xml_node node) +{ + // Check for particle type + if (check_for_node(node, "particle")) { + auto temp_str = get_node_value(node, "particle", true, true); + if (temp_str == "neutron") { + particle_ = ParticleType::neutron; + } else if (temp_str == "photon") { + particle_ = ParticleType::photon; + openmc_photon_transport = true; + } else { + fatal_error(std::string("Unknown source particle type: ") + temp_str); + } + } + + // Check for source strength + if (check_for_node(node, "strength")) { + strength_ = std::stod(get_node_value(node, "strength")); + } + + // Check for external source file + if (check_for_node(node, "file")) { + // Copy path of source file + path_source = get_node_value(node, "file", false, true); + + // Check if source file exists + if (!file_exists(path_source)) { + std::stringstream msg; + msg << "Source file '" << path_source << "' does not exist."; + fatal_error(msg); + } + + } else { + + // Spatial distribution for external source + if (check_for_node(node, "space")) { + // Get pointer to spatial distribution + pugi::xml_node node_space = node.child("space"); + + // Check for type of spatial distribution and read + std::string type; + if (check_for_node(node_space, "type")) + type = get_node_value(node_space, "type", true, true); + if (type == "cartesian") { + space_ = UPtrSpace{new CartesianIndependent(node_space)}; + } else if (type == "box") { + space_ = UPtrSpace{new SpatialBox(node_space)}; + } else if (type == "fission") { + space_ = UPtrSpace{new SpatialBox(node_space, true)}; + } else if (type == "point") { + space_ = UPtrSpace{new SpatialPoint(node_space)}; + } else { + std::stringstream msg; + msg << "Invalid spatial distribution for external source: " << type; + fatal_error(msg); + } + + } else { + // If no spatial distribution specified, make it a point source + space_ = UPtrSpace{new SpatialPoint()}; + } + + // Determine external source angular distribution + if (check_for_node(node, "angle")) { + // Get pointer to angular distribution + pugi::xml_node node_angle = node.child("angle"); + + // Check for type of angular distribution + std::string type; + if (check_for_node(node_angle, "type")) + type = get_node_value(node_angle, "type", true, true); + if (type == "isotropic") { + angle_ = UPtrAngle{new Isotropic()}; + } else if (type == "monodirectional") { + angle_ = UPtrAngle{new Monodirectional(node_angle)}; + } else if (type == "mu-phi") { + angle_ = UPtrAngle{new PolarAzimuthal(node_angle)}; + } else { + std::stringstream msg; + msg << "Invalid angular distribution for external source: " << type; + fatal_error(msg); + } + + } else { + angle_ = UPtrAngle{new Isotropic()}; + } + + // Determine external source energy distribution + if (check_for_node(node, "energy")) { + pugi::xml_node node_dist = node.child("energy"); + energy_ = distribution_from_xml(node_dist); + } else { + // Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1 + energy_ = UPtrDist{new Watt(0.988e6, 2.249e-6)}; + } + } +} + + +Bank SourceDistribution::sample() const +{ + Bank site; + + // Set weight to one by default + site.wgt = 1.0; + + // Repeat sampling source location until a good site has been found + bool found = false; + int n_reject = 0; + static int n_accept = 0; + while (!found) { + // Set particle type + site.particle = static_cast(particle_); + + // Sample spatial distribution + Position r = space_->sample(); + site.xyz[0] = r.x; + site.xyz[1] = r.y; + site.xyz[2] = r.z; + + // Now search to see if location exists in geometry + int32_t cell_index, instance; + int err = openmc_find_cell(site.xyz, &cell_index, &instance); + found = (err != OPENMC_E_GEOMETRY); + + // Check if spatial site is in fissionable material + if (found) { + auto space_box = dynamic_cast(space_.get()); + if (space_box) { + if (space_box->only_fissionable()) { + // Determine material + auto c = global_cells[cell_index - 1]; + int32_t mat_index = c->material[instance]; + auto m = global_materials[mat_index]; + + if (mat_index == MATERIAL_VOID) { + found = false; + } else { + bool fissionable; + openmc_material_get_fissionable(mat_index + 1, &fissionable); + if (!fissionable) found = false; + } + } + } + } + + // Check for rejection + if (!found) { + ++n_reject; + if (n_reject >= EXTSRC_REJECT_THRESHOLD && + static_cast(n_accept)/n_reject <= EXTSRC_REJECT_FRACTION) { + fatal_error("More than 95% of external source sites sampled were " + "rejected. Please check your external source definition."); + } + } + } + + // Increment number of accepted samples + ++n_accept; + + // Sample angle + Direction u = angle_->sample(); + site.uvw[0] = u.x; + site.uvw[1] = u.y; + site.uvw[2] = u.z; + + // Check for monoenergetic source above maximum particle energy + auto p = static_cast(particle_); + auto energy_ptr = dynamic_cast(energy_.get()); + if (energy_ptr) { + auto energies = xt::adapt(energy_ptr->x()); + if (xt::any(energies > energy_max[p-1])) { + fatal_error("Source energy above range of energies of at least " + "one cross section table"); + } else if (xt::any(energies < energy_min[p-1])) { + fatal_error("Source energy below range of energies of at least " + "one cross section table"); + } + } + + while (true) { + // Sample energy spectrum + site.E = energy_->sample(); + + // Resample if energy falls outside minimum or maximum particle energy + if (site.E < energy_max[p-1] && site.E > energy_min[p-1]) break; + } + + // Set delayed group + site.delayed_group = 0; + + return site; +} + +//============================================================================== +// Non-member functions +//============================================================================== + +void initialize_source() +{ + write_message("Initializing source particles...", 5); + + // Get pointer to source bank + Bank* source_bank; + int64_t n; + openmc_source_bank(&source_bank, &n); + + if (path_source != "") { + // Read the source from a binary file instead of sampling from some + // assumed source distribution + + std::stringstream msg; + msg << "Reading source file from " << path_source << "..."; + write_message(msg, 6); + + // Open the binary file + hid_t file_id = file_open(path_source, 'r', true); + + // Read the file type + std::string filetype; + read_attribute(file_id, "filetype", filetype); + + // Check to make sure this is a source file + if (filetype != "source" && filetype != "statepoint") { + fatal_error("Specified starting source file not a source file type."); + } + + // Read in the source bank + read_source_bank(file_id, work_index.data(), source_bank); + + // Close file + file_close(file_id); + + } else { + // Generation source sites from specified distribution in user input + for (int64_t i = 0; i < openmc_work; ++i) { + // initialize random number seed + int64_t id = openmc_total_gen*n_particles + work_index[openmc::mpi::rank] + i + 1; + set_particle_seed(id); + + // sample external source distribution + source_bank[i] = sample_external_source(); + } + } + + // Write out initial source + if (openmc_write_initial_source) { + write_message("Writing out initial source...", 5); + std::string filename = path_output + "initial_source.h5"; + hid_t file_id = file_open(filename, 'w', true); + write_source_bank(file_id, work_index.data(), source_bank); + file_close(file_id); + } +} + +extern "C" double* rev_energy_bins_ptr(); + +Bank sample_external_source() +{ + // Set the random number generator to the source stream. + prn_set_stream(STREAM_SOURCE); + + // Determine total source strength + double total_strength = 0.0; + for (auto& s : external_sources) + total_strength += s.strength(); + + // Sample from among multiple source distributions + auto n_source = external_sources.size(); + int i = 0; + if (n_source > 1) { + double xi = prn()*total_strength; + double c = 0.0; + for (i = 0; i < external_sources.size(); ++i) { + c += external_sources[i].strength(); + if (xi < c) break; + } + } + + // Sample source site from i-th source distribution + Bank site {external_sources[i].sample()}; + + // If running in MG, convert site % E to group + if (!openmc_run_CE) { + // Get pointer to rev_energy_bins array on Fortran side + double* rev_energy_bins = rev_energy_bins_ptr(); + + int n = num_energy_groups + 1; + site.E = lower_bound_index(rev_energy_bins, rev_energy_bins + n, site.E); + site.E = num_energy_groups - site.E; + } + + // Set the random number generator back to the tracking stream. + prn_set_stream(STREAM_TRACKING); + + return site; +} + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" void free_memory_source() +{ + external_sources.clear(); +} + +extern "C" double total_source_strength() +{ + double strength = 0.0; + for (const auto& s : external_sources) { + strength += s.strength(); + } + return strength; +} + +} // namespace openmc diff --git a/src/source_header.F90 b/src/source_header.F90 deleted file mode 100644 index 163ae1c19..000000000 --- a/src/source_header.F90 +++ /dev/null @@ -1,388 +0,0 @@ -module source_header - - use, intrinsic :: ISO_C_BINDING - - use bank_header, only: Bank - use constants - use distribution_univariate - use distribution_multivariate - use error - use geometry, only: find_cell - use material_header, only: materials - use nuclide_header, only: energy_min, energy_max - use particle_header - use settings, only: photon_transport - use string, only: to_lower - use xml_interface - - implicit none - private - public :: free_memory_source - public :: openmc_extend_sources - public :: openmc_source_set_strength - - integer :: n_accept = 0 ! Number of samples accepted - integer :: n_reject = 0 ! Number of samples rejected - -!=============================================================================== -! SOURCEDISTRIBUTION describes an external source of particles for a -! fixed-source problem or for the starting source in a k eigenvalue problem -!=============================================================================== - - type, public :: SourceDistribution - integer :: particle ! particle type - real(8) :: strength = ONE ! source strength - class(SpatialDistribution), allocatable :: space ! spatial distribution - class(UnitSphereDistribution), allocatable :: angle ! angle distribution - class(Distribution), allocatable :: energy ! energy distribution - contains - procedure :: from_xml - procedure :: sample - end type SourceDistribution - - ! Number of external source distributions - integer(C_INT32_T), public, bind(C) :: n_sources = 0 - - ! External source distributions - type(SourceDistribution), public, allocatable :: external_source(:) - -contains - - subroutine from_xml(this, node, path_source) - class(SourceDistribution), intent(inout) :: this - type(XMLNode), intent(in) :: node - character(MAX_FILE_LEN), intent(out) :: path_source - - integer :: n - logical :: file_exists - character(MAX_WORD_LEN) :: type, temp_str - type(XMLNode) :: node_space - type(XMLNode) :: node_angle - type(XMLNode) :: node_dist - - ! Check for particle type - if (check_for_node(node, "particle")) then - call get_node_value(node, "particle", temp_str) - select case (to_lower(temp_str)) - case ('neutron') - this % particle = NEUTRON - case ('photon') - this % particle = PHOTON - photon_transport = .true. - case default - call fatal_error('Unknown source particle type: ' // trim(temp_str)) - end select - else - this % particle = NEUTRON - end if - - ! Check for source strength - if (check_for_node(node, "strength")) then - call get_node_value(node, "strength", this % strength) - end if - - ! Check for external source file - if (check_for_node(node, "file")) then - ! Copy path of source file - call get_node_value(node, "file", path_source) - - ! Check if source file exists - inquire(FILE=path_source, EXIST=file_exists) - if (.not. file_exists) then - call fatal_error("Source file '" // trim(path_source) & - // "' does not exist!") - end if - - else - - ! Spatial distribution for external source - if (check_for_node(node, "space")) then - - ! Get pointer to spatial distribution - node_space = node % child("space") - - ! Check for type of spatial distribution - type = '' - if (check_for_node(node_space, "type")) & - call get_node_value(node_space, "type", type) - select case (to_lower(type)) - case ('cartesian') - allocate(CartesianIndependent :: this % space) - - case ('box') - allocate(SpatialBox :: this % space) - - case ('fission') - allocate(SpatialBox :: this % space) - select type(space => this % space) - type is (SpatialBox) - space % only_fissionable = .true. - end select - - case ('point') - allocate(SpatialPoint :: this % space) - - case default - call fatal_error("Invalid spatial distribution for external source: "& - // trim(type)) - end select - - ! Read spatial distribution from XML - call this % space % from_xml(node_space) - - else - ! If no spatial distribution specified, make it a point source - allocate(SpatialPoint :: this % space) - select type (space => this % space) - type is (SpatialPoint) - space % xyz(:) = [ZERO, ZERO, ZERO] - end select - end if - - ! Determine external source angular distribution - if (check_for_node(node, "angle")) then - - ! Get pointer to angular distribution - node_angle = node % child("angle") - - ! Check for type of angular distribution - type = '' - if (check_for_node(node_angle, "type")) & - call get_node_value(node_angle, "type", type) - select case (to_lower(type)) - case ('isotropic') - allocate(Isotropic :: this % angle) - - case ('monodirectional') - allocate(Monodirectional :: this % angle) - - case ('mu-phi') - allocate(PolarAzimuthal :: this % angle) - - case default - call fatal_error("Invalid angular distribution for external source: "& - // trim(type)) - end select - - ! Read reference directional unit vector - if (check_for_node(node_angle, "reference_uvw")) then - n = node_word_count(node_angle, "reference_uvw") - if (n /= 3) then - call fatal_error('Angular distribution reference direction must have & - &three parameters specified.') - end if - call get_node_array(node_angle, "reference_uvw", & - this % angle % reference_uvw) - else - ! By default, set reference unit vector to be positive z-direction - this % angle % reference_uvw(:) = [ZERO, ZERO, ONE] - end if - - ! Read parameters for angle distribution - select type (angle => this % angle) - type is (Monodirectional) - call get_node_array(node_angle, "reference_uvw", & - this % angle % reference_uvw) - - type is (PolarAzimuthal) - if (check_for_node(node_angle, "mu")) then - node_dist = node_angle % child("mu") - call distribution_from_xml(angle % mu, node_dist) - else - allocate(Uniform :: angle%mu) - select type (mu => angle%mu) - type is (Uniform) - mu % a = -ONE - mu % b = ONE - end select - end if - - if (check_for_node(node_angle, "phi")) then - node_dist = node_angle % child("phi") - call distribution_from_xml(angle % phi, node_dist) - else - allocate(Uniform :: angle%phi) - select type (phi => angle%phi) - type is (Uniform) - phi % a = ZERO - phi % b = TWO*PI - end select - end if - end select - - else - ! Set default angular distribution isotropic - allocate(Isotropic :: this % angle) - this % angle % reference_uvw(:) = [ZERO, ZERO, ONE] - end if - - ! Determine external source energy distribution - if (check_for_node(node, "energy")) then - node_dist = node % child("energy") - call distribution_from_xml(this % energy, node_dist) - else - ! Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1 - allocate(Watt :: this % energy) - select type(energy => this % energy) - type is (Watt) - energy % a = 0.988e6_8 - energy % b = 2.249e-6_8 - end select - end if - end if - - end subroutine from_xml - - function sample(this) result(site) - class(SourceDistribution), intent(in) :: this - type(Bank) :: site - - logical :: found ! Does the source particle exist within geometry? - type(Particle) :: p ! Temporary particle for using find_cell - - ! Set weight to one by default - site % wgt = ONE - - ! Repeat sampling source location until a good site has been found - found = .false. - do while (.not. found) - ! Set particle defaults - call particle_initialize(p) - - ! Set particle type - site % particle = this % particle - - ! Sample spatial distribution - site % xyz(:) = this % space % sample() - - ! Fill p with needed data - p % coord(1) % xyz(:) = site % xyz - p % coord(1) % uvw(:) = [ ONE, ZERO, ZERO ] - - ! Now search to see if location exists in geometry - call find_cell(p, found) - - ! Check if spatial site is in fissionable material - if (found) then - select type (space => this % space) - type is (SpatialBox) - if (space % only_fissionable) then - if (p % material == MATERIAL_VOID) then - found = .false. - elseif (.not. materials(p % material) % fissionable) then - found = .false. - end if - end if - end select - end if - - ! Check for rejection - if (.not. found) then - n_reject = n_reject + 1 - if (n_reject >= EXTSRC_REJECT_THRESHOLD .and. & - real(n_accept, 8)/n_reject <= EXTSRC_REJECT_FRACTION) then - call fatal_error("More than 95% of external source sites sampled & - &were rejected. Please check your external source definition.") - end if - end if - end do - - ! Increment number of accepted samples - n_accept = n_accept + 1 - - call particle_clear(p) - - ! Sample angle - site % uvw(:) = this % angle % sample() - - ! Check for monoenergetic source above maximum particle energy - select type (energy => this % energy) - type is (Discrete) - if (any(energy % x > energy_max(this % particle))) then - call fatal_error("Source energy above range of energies of at least & - &one cross section table") - else if (any(energy % x < energy_min(this % particle))) then - call fatal_error("Source energy below range of energies of at least & - &one cross section table") - end if - end select - - do - ! Sample energy spectrum - site % E = this % energy % sample() - - ! Resample if energy falls outside minimum or maximum particle energy - if (site % E < energy_max(this % particle) .and. & - site % E > energy_min(this % particle)) exit - end do - - ! Set delayed group - site % delayed_group = 0 - - end function sample - -!=============================================================================== -! FREE_MEMORY_SOURCE deallocates global arrays defined in this module -!=============================================================================== - - subroutine free_memory_source() - n_sources = 0 - if (allocated(external_source)) deallocate(external_source) - end subroutine free_memory_source - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_extend_sources(n, index_start, index_end) result(err) bind(C) - ! Extend the external_source array by n elements - integer(C_INT32_T), value, intent(in) :: n - integer(C_INT32_T), optional, intent(out) :: index_start - integer(C_INT32_T), optional, intent(out) :: index_end - integer(C_INT) :: err - - type(SourceDistribution), allocatable :: temp(:) ! temporary array - - if (n_sources == 0) then - ! Allocate external_source array - allocate(external_source(n)) - else - ! Allocate external_source array with increased size - allocate(temp(n_sources + n)) - - ! Copy original source array to temporary array - temp(1:n_sources) = external_source - - ! Move allocation from temporary array - call move_alloc(FROM=temp, TO=external_source) - end if - - ! Return indices in external_source array - if (present(index_start)) index_start = n_sources + 1 - if (present(index_end)) index_end = n_sources + n - n_sources = n_sources + n - - err = 0 - end function openmc_extend_sources - - - function openmc_source_set_strength(index, strength) result(err) bind(C) - integer(C_INT32_T), value, intent(in) :: index - real(C_DOUBLE), value, intent(in) :: strength - integer(C_INT) :: err - - if (index >= 1 .and. index <= n_sources) then - if (strength > ZERO) then - external_source(index) % strength = strength - err = 0 - else - err = E_INVALID_ARGUMENT - call set_errmsg("Source strength must be positive.") - end if - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in external source array is out of bounds.") - end if - end function openmc_source_set_strength - -end module source_header diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 32f094f89..367214a93 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -9,7 +9,6 @@ module tally_header use message_passing, only: n_procs use nuclide_header, only: nuclide_dict use settings, only: reduce_tallies, run_mode - use source_header, only: external_source use stl_vector, only: VectorInt use string, only: to_lower, to_f_string, str_to_int, to_str use tally_filter_header, only: TallyFilterContainer, filters, n_filters @@ -169,6 +168,13 @@ contains real(C_DOUBLE) :: val real(C_DOUBLE) :: total_source + interface + function total_source_strength() result(strength) bind(C) + import C_DOUBLE + real(C_DOUBLE) :: strength + end function + end interface + ! Increment number of realizations if (reduce_tallies) then this % n_realizations = this % n_realizations + 1 @@ -178,7 +184,7 @@ contains ! Calculate total source strength for normalization if (run_mode == MODE_FIXEDSOURCE) then - total_source = sum(external_source(:) % strength) + total_source = total_source_strength() else total_source = ONE end if diff --git a/src/xml_interface.F90 b/src/xml_interface.F90 index 0e08f75d1..7aa86556b 100644 --- a/src/xml_interface.F90 +++ b/src/xml_interface.F90 @@ -18,6 +18,7 @@ module xml_interface interface get_node_value module procedure get_node_value_bool + module procedure get_node_value_cbool module procedure get_node_value_integer module procedure get_node_value_long module procedure get_node_value_double @@ -118,6 +119,16 @@ contains end if end subroutine get_node_value_bool + subroutine get_node_value_cbool(node, name, val) + type(XMLNode), intent(in) :: node + character(*), intent(in) :: name + logical(kind=C_BOOL), intent(out) :: val + + logical :: val_ + call get_node_value_bool(node, name, val_) + val = val_ + end subroutine get_node_value_cbool + !=============================================================================== ! GET_NODE_VALUE_INTEGER returns a integer value from an attribute or node !===============================================================================