diff --git a/CMakeLists.txt b/CMakeLists.txt index bb5c45c427..55f4d7f4ee 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -381,6 +381,7 @@ add_library(libopenmc SHARED src/tallies/tally_header.F90 src/tallies/trigger.F90 src/tallies/trigger_header.F90 + src/bank.cpp src/dagmc.cpp src/cell.cpp src/cmfd_execute.cpp @@ -408,6 +409,8 @@ add_library(libopenmc SHARED src/nuclide.cpp src/output.cpp src/particle.cpp + src/photon.cpp + src/physics.cpp src/physics_common.cpp src/physics_mg.cpp src/plot.cpp diff --git a/include/openmc/bank.h b/include/openmc/bank.h new file mode 100644 index 0000000000..b5dcb4d911 --- /dev/null +++ b/include/openmc/bank.h @@ -0,0 +1,47 @@ +#ifndef OPENMC_BANK_H +#define OPENMC_BANK_H + +#include +#include + +namespace openmc { + +struct Bank { + double wgt; + double xyz[3]; + double uvw[3]; + double E; + int delayed_group; + int particle; +}; + +} // namespace openmc + +// Without an explicit instantiation of vector, the Intel compiler +// will complain about the threadprivate directive on filter_matches. Note that +// this has to happen *outside* of the openmc namespace +extern template class std::vector; + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +namespace simulation { + +extern "C" int64_t n_bank; + +extern std::vector source_bank; +extern std::vector fission_bank; +#ifdef _OPENMP +extern std::vector master_fission_bank; +#endif + +#pragma omp threadprivate(fission_bank, n_bank) + +} // namespace simulation + +} // namespace openmc + +#endif // OPENMC_BANK_H diff --git a/include/openmc/capi.h b/include/openmc/capi.h index fdea48e9e1..50ae72000f 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -6,9 +6,11 @@ #include #ifdef __cplusplus +#include "openmc/bank.h" extern "C" { -#endif - + int openmc_fission_bank(openmc::Bank** ptr, int64_t* n); + int openmc_source_bank(openmc::Bank** ptr, int64_t* n); +#else struct Bank { double wgt; double xyz[3]; @@ -18,6 +20,10 @@ extern "C" { int particle; }; + int openmc_fission_bank(struct Bank** ptr, int64_t* n); + int openmc_source_bank(struct Bank** ptr, int64_t* n); +#endif + int openmc_calculate_volumes(); int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n); int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n); @@ -38,7 +44,6 @@ extern "C" { int openmc_filter_set_type(int32_t index, const char* type); int openmc_finalize(); int openmc_find_cell(double* xyz, int32_t* index, int32_t* instance); - int openmc_fission_bank(struct Bank** ptr, int64_t* n); int openmc_get_cell_index(int32_t id, int32_t* index); int openmc_get_filter_index(int32_t id, int32_t* index); void openmc_get_filter_next_id(int32_t* id); @@ -85,7 +90,6 @@ extern "C" { void openmc_set_seed(int64_t new_seed); int openmc_simulation_finalize(); int openmc_simulation_init(); - int openmc_source_bank(struct Bank** ptr, int64_t* n); 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/eigenvalue.h b/include/openmc/eigenvalue.h index 6d918eb403..f7e93ec2b7 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -25,9 +25,6 @@ extern std::array k_sum; //!< Used to reduce sum and sum_sq extern std::vector entropy; //!< Shannon entropy at each generation extern xt::xtensor source_frac; //!< Source fraction for UFS -extern "C" int64_t n_bank; -#pragma omp threadprivate(n_bank) - } // namespace simulation //============================================================================== @@ -35,15 +32,22 @@ extern "C" int64_t n_bank; //============================================================================== //! Collect/normalize the tracklength keff from each process -extern "C" void calculate_generation_keff(); +void calculate_generation_keff(); //! Calculate mean/standard deviation of keff during active generations //! //! This function sets the global variables keff and keff_std which represent //! the mean and standard deviation of the mean of k-effective over active //! generations. It also broadcasts the value from the master process. -extern "C" void calculate_average_keff(); +void calculate_average_keff(); +#ifdef _OPENMP +//! Join threadprivate fission banks into a single fission bank +//! +//! Note that this operation is necessarily sequential to preserve the order of +//! the bank when using varying numbers of threads. +void join_bank_from_threads(); +#endif //! Calculates a minimum variance estimate of k-effective //! @@ -60,16 +64,16 @@ extern "C" void calculate_average_keff(); extern "C" int openmc_get_keff(double* k_combined); //! Sample/redistribute source sites from accumulated fission sites -extern "C" void synchronize_bank(); +void synchronize_bank(); //! Calculates the Shannon entropy of the fission source distribution to assess //! source convergence -extern "C" void shannon_entropy(); +void shannon_entropy(); //! Determines the source fraction in each UFS mesh cell and reweights the //! source bank so that the sum of the weights is equal to n_particles. The //! 'source_frac' variable is used later to bias the production of fission sites -extern "C" void ufs_count_sites(); +void ufs_count_sites(); //! Get UFS weight corresponding to particle's location extern "C" double ufs_get_weight(const Particle* p); diff --git a/include/openmc/endf.h b/include/openmc/endf.h index c781bd0b84..a7030d0355 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -22,6 +22,17 @@ Interpolation int2interp(int i); //! \return Whether corresponding reaction is a fission reaction bool is_fission(int MT); +//! Determine if a given MT number is that of a disappearance reaction, i.e., a +//! reaction with no neutron in the exit channel +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is a disappearance reaction +bool is_disappearance(int MT); + +//! Determine if a given MT number is that of an inelastic scattering reaction +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is an inelastic scattering reaction +bool is_inelastic_scatter(int MT); + //============================================================================== //! Abstract one-dimensional function //============================================================================== diff --git a/include/openmc/material.h b/include/openmc/material.h index fcab023a52..ad9b9c9f5a 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -43,5 +43,11 @@ public: explicit Material(pugi::xml_node material_node); }; +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" bool material_isotropic(int i_material, int i_nuc_mat); + } // namespace openmc #endif // OPENMC_MATERIAL_H diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index db30a1a50a..9a99f71e82 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -5,23 +5,75 @@ #define OPENMC_NUCLIDE_H #include +#include // for unique_ptr +#include + +#include #include "openmc/constants.h" +#include "openmc/endf.h" +#include "openmc/reaction.h" +#include "openmc/reaction_product.h" namespace openmc { //============================================================================== -// Global variables +// Constants //============================================================================== -namespace data { +constexpr double CACHE_INVALID {-1.0}; -// Minimum/maximum transport energy for each particle type. Order corresponds to -// that of the ParticleType enum -extern std::array energy_min; -extern std::array energy_max; +//=============================================================================== +// Data for a nuclide +//=============================================================================== -} // namespace data +class Nuclide { +public: + // Types, aliases + using EmissionMode = ReactionProduct::EmissionMode; + struct EnergyGrid { + std::vector grid_index; + std::vector energy; + }; + + // Constructors + Nuclide(hid_t group, const double* temperature, int n); + + // Methods + double nu(double E, EmissionMode mode, int group=0) const; + void calculate_elastic_xs(int i_nuclide) const; + + //! Determines the microscopic 0K elastic cross section at a trial relative + //! energy used in resonance scattering + double elastic_xs_0K(double E) const; + + // Data members + std::string name_; //! Name of nuclide, e.g. "U235" + int Z_; //! Atomic number + int A_; //! Mass number + int metastable_; //! Metastable state + double awr_; //! Atomic weight ratio + std::vector kTs_; //! temperatures in eV (k*T) + std::vector grid_; //! Energy grid at each temperature + + bool fissionable_ {false}; //! Whether nuclide is fissionable + bool has_partial_fission_ {false}; //! has partial fission reactions? + std::vector fission_rx_; //! Fission reactions + int n_precursor_ {0}; //! Number of delayed neutron precursors + std::unique_ptr total_nu_; //! Total neutron yield + + // Resonance scattering information + bool resonant_ {false}; + std::vector energy_0K_; + std::vector elastic_0K_; + std::vector xs_cdf_; + + std::vector> reactions_; //! Reactions + std::vector index_inelastic_scatter_; + +private: + void create_derived(); +}; //=============================================================================== //! Cached microscopic cross sections for a particular nuclide at the current @@ -66,19 +118,53 @@ struct NuclideMicroXS { // particle is traveling through //=============================================================================== - struct MaterialMacroXS { - double total; //!< macroscopic total xs - double absorption; //!< macroscopic absorption xs - double fission; //!< macroscopic fission xs - double nu_fission; //!< macroscopic production xs - double photon_prod; //!< macroscopic photon production xs +struct MaterialMacroXS { + double total; //!< macroscopic total xs + double absorption; //!< macroscopic absorption xs + double fission; //!< macroscopic fission xs + double nu_fission; //!< macroscopic production xs + double photon_prod; //!< macroscopic photon production xs + + // Photon cross sections + double coherent; //!< macroscopic coherent xs + double incoherent; //!< macroscopic incoherent xs + double photoelectric; //!< macroscopic photoelectric xs + double pair_production; //!< macroscopic pair production xs +}; + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +// Minimum/maximum transport energy for each particle type. Order corresponds to +// that of the ParticleType enum +extern std::array energy_min; +extern std::array energy_max; + +extern std::vector> nuclides; + +} // namespace data + +namespace simulation { + +// Cross section caches +extern NuclideMicroXS* micro_xs; +extern "C" MaterialMacroXS material_xs; +#pragma omp threadprivate(micro_xs, material_xs) + +} // namespace simulation + +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" void set_micro_xs(); +extern "C" bool nuclide_wmp_present(int i_nuclide); +extern "C" double nuclide_wmp_emin(int i_nuclide); +extern "C" double nuclide_wmp_emax(int i_nuclide); - // Photon cross sections - double coherent; //!< macroscopic coherent xs - double incoherent; //!< macroscopic incoherent xs - double photoelectric; //!< macroscopic photoelectric xs - double pair_production; //!< macroscopic pair production xs - }; } // namespace openmc diff --git a/include/openmc/photon.h b/include/openmc/photon.h new file mode 100644 index 0000000000..3f36decdb3 --- /dev/null +++ b/include/openmc/photon.h @@ -0,0 +1,37 @@ +#ifndef OPENMC_PHOTON_H +#define OPENMC_PHOTON_H + +#include + +#include +#include + +namespace openmc { + +//============================================================================== +//! Photon interaction data for a single element +//============================================================================== + +class PhotonInteraction { +public: + // Constructors + PhotonInteraction(hid_t group); + + // Data members + std::string name_; //! Name of element, e.g. "Zr" + int Z_; //! Atomic number +}; + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +extern std::vector elements; + +} // namespace data + +} // namespace openmc + +#endif // OPENMC_PHOTON_H diff --git a/include/openmc/physics.h b/include/openmc/physics.h new file mode 100644 index 0000000000..58bc7e9906 --- /dev/null +++ b/include/openmc/physics.h @@ -0,0 +1,89 @@ +#ifndef OPENMC_PHYSICS_H +#define OPENMC_PHYSICS_H + +#include "openmc/bank.h" +#include "openmc/nuclide.h" +#include "openmc/particle.h" +#include "openmc/position.h" +#include "openmc/reaction.h" + +namespace openmc { + +//============================================================================== +// Non-member functions +//============================================================================== + +//! Sample a nuclide and reaction and then calls the appropriate routine +extern "C" void collision(Particle* p); + +//! Samples an incident neutron reaction +void sample_neutron_reaction(Particle* p); + +//! Samples an element based on the macroscopic cross sections for each nuclide +//! within a material and then samples a reaction for that element and calls the +//! appropriate routine to process the physics. +extern "C" void sample_photon_reaction(Particle* p); + +//! Terminates the particle and either deposits all energy locally +//! (electron_treatment = ELECTRON_LED) or creates secondary bremsstrahlung +//! photons from electron deflections with charged particles (electron_treatment +//! = ELECTRON_TTB). +void sample_electron_reaction(Particle* p); + +//! Terminates the particle and either deposits all energy locally +//! (electron_treatment = ELECTRON_LED) or creates secondary bremsstrahlung +//! photons from electron deflections with charged particles (electron_treatment +//! = ELECTRON_TTB). Two annihilation photons of energy MASS_ELECTRON_EV (0.511 +//! MeV) are created and travel in opposite directions. +void sample_positron_reaction(Particle* p); + +void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat); + +//! Determine the average total, prompt, and delayed neutrons produced from +//! fission and creates appropriate bank sites. +void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, + Bank* bank_array, int64_t* bank_size, int64_t bank_capacity); + +// void sample_element(Particle* p); + +Reaction* sample_fission(int i_nuclide, double E); + +void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product); + +void absorption(Particle* p, int i_nuclide); + +void scatter(Particle*, int i_nuclide, int i_nuc_mat); + +//! Treats the elastic scattering of a neutron with a target. +void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, + double* uvw, double* mu_lab, double* wgt); + +extern "C" void sab_scatter(int i_nuclide, int i_sab, double* E, + double* uvw, double* mu); + +//! samples the target velocity. The constant cross section free gas model is +//! the default method. Methods for correctly accounting for the energy +//! dependence of cross sections in treating resonance elastic scattering such +//! as the DBRC, WCM, and a new, accelerated scheme are also implemented here. +Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, + Direction v_neut, double xs_eff, double kT, double* wgt); + +//! samples a target velocity based on the free gas scattering formulation, used +//! by most Monte Carlo codes, in which cross section is assumed to be constant +//! in energy. Excellent documentation for this method can be found in +//! FRA-TM-123. +Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT); + +void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site); + +//! handles all reactions with a single secondary neutron (other than fission), +//! i.e. level scattering, (n,np), (n,na), etc. +void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p); + +void sample_secondary_photons(Particle* p, int i_nuclide); + +extern "C" void thick_target_bremsstrahlung(Particle* p, double* E_lost); + +} // namespace openmc + +#endif // OPENMC_PHYSICS_H diff --git a/include/openmc/position.h b/include/openmc/position.h index 62c7f804fc..0deb298dd2 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -25,6 +25,8 @@ struct Position { Position& operator-=(double); Position& operator*=(Position); Position& operator*=(double); + Position& operator/=(Position); + Position& operator/=(double); const double& operator[](int i) const { switch (i) { @@ -76,6 +78,10 @@ inline Position operator*(Position a, Position b) { return a *= b; } inline Position operator*(Position a, double b) { return a *= b; } inline Position operator*(double a, Position b) { return b *= a; } +inline Position operator/(Position a, Position b) { return a /= b; } +inline Position operator/(Position a, double b) { return a /= b; } +inline Position operator/(double a, Position b) { return b /= a; } + inline bool operator==(Position a, Position b) {return a.x == b.x && a.y == b.y && a.z == b.z;} diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 1015a4ad27..7f20e7a956 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -4,6 +4,7 @@ #ifndef OPENMC_REACTION_H #define OPENMC_REACTION_H +#include #include #include "hdf5.h" @@ -39,13 +40,17 @@ public: std::vector products_; //!< Reaction products }; +//============================================================================== +// Non-member functions +//============================================================================== + +std::string reaction_name(int mt); + //============================================================================== // Fortran compatibility functions //============================================================================== extern "C" { - Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n); - void reaction_delete(Reaction* rx); int reaction_mt(Reaction* rx); double reaction_q_value(Reaction* rx); bool reaction_scatter_in_cm(Reaction* rx); diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 2002e95fed..cba0a2b713 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -56,6 +56,9 @@ extern "C" int thread_id; //!< ID of a given thread // Functions //============================================================================== +//! Allocate space for source and fission banks +void allocate_banks(); + //! Determine number of particles to transport per process void calculate_work(); diff --git a/include/openmc/source.h b/include/openmc/source.h index 45f4896b46..d9f5cbefa0 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -9,9 +9,9 @@ #include "pugixml.hpp" +#include "openmc/bank.h" #include "openmc/distribution_multi.h" #include "openmc/distribution_spatial.h" -#include "openmc/capi.h" #include "openmc/particle.h" namespace openmc { @@ -62,7 +62,7 @@ extern "C" void initialize_source(); //! Sample a site from all external source distributions in proportion to their //! source strength //! \return Sampled source site -extern "C" Bank sample_external_source(); +Bank sample_external_source(); //! Fill source bank at end of generation for fixed source simulations void fill_source_bank_fixedsource(); diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 6b019e7b8f..0c1db86f3f 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -10,8 +10,8 @@ namespace openmc { void write_source_point(const char* filename); -extern "C" void write_source_bank(hid_t group_id, Bank* source_bank); -extern "C" void read_source_bank(hid_t group_id, Bank* source_bank); +extern "C" void write_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(); diff --git a/openmc/capi/__init__.py b/openmc/capi/__init__.py index 95ac77a64c..1015df2166 100644 --- a/openmc/capi/__init__.py +++ b/openmc/capi/__init__.py @@ -38,8 +38,11 @@ else: from unittest.mock import Mock _dll = Mock() -dagmc_enabled = bool(c_bool.in_dll(_dll, "dagmc_enabled")) - + +def _dagmc_enabled(): + return c_bool.in_dll(_dll, "dagmc_enabled").value + + from .error import * from .core import * from .nuclide import * diff --git a/src/api.F90 b/src/api.F90 index efcfa16a6a..d7fedcbe30 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -82,6 +82,9 @@ contains subroutine free_memory_settings() bind(C) end subroutine free_memory_settings + + subroutine free_memory_bank() bind(C) + end subroutine free_memory_bank end interface call free_memory_geometry() diff --git a/src/bank.cpp b/src/bank.cpp new file mode 100644 index 0000000000..89c493babd --- /dev/null +++ b/src/bank.cpp @@ -0,0 +1,107 @@ +#include "openmc/bank.h" + +#include "openmc/capi.h" +#include "openmc/error.h" + +#include + +// Explicit template instantiation definition +template class std::vector; + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +namespace simulation { + +int64_t n_bank; + +std::vector source_bank; +std::vector fission_bank; +#ifdef _OPENMP +std::vector master_fission_bank; +#endif + +} // namespace simulation + +//============================================================================== +// C API +//============================================================================== + +extern "C" int openmc_source_bank(Bank** ptr, int64_t* n) +{ + if (simulation::source_bank.size() == 0) { + set_errmsg("Source bank has not been allocated."); + return OPENMC_E_ALLOCATE; + } else { + *ptr = simulation::source_bank.data(); + *n = simulation::source_bank.size(); + return 0; + } +} + +extern "C" int openmc_fission_bank(Bank** ptr, int64_t* n) +{ + if (simulation::fission_bank.size() == 0) { + set_errmsg("Fission bank has not been allocated."); + return OPENMC_E_ALLOCATE; + } else { + *ptr = simulation::fission_bank.data(); + *n = simulation::fission_bank.size(); + return 0; + } +} + +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" void free_memory_bank() +{ + simulation::source_bank.clear(); +#pragma omp parallel + { + simulation::fission_bank.clear(); + } +#ifdef _OPENMP + simulation::master_fission_bank.clear(); +#endif +} + +extern "C" int fission_bank_delayed_group(int64_t i) { + return simulation::fission_bank[i-1].delayed_group; +} + +extern "C" double fission_bank_E(int64_t i) { + return simulation::fission_bank[i-1].E; +} + +extern "C" double fission_bank_wgt(int64_t i) { + return simulation::fission_bank[i-1].wgt; +} + +extern "C" void source_bank_xyz(int64_t i, double* xyz) +{ + xyz[0] = simulation::source_bank[i-1].xyz[0]; + xyz[1] = simulation::source_bank[i-1].xyz[1]; + xyz[2] = simulation::source_bank[i-1].xyz[2]; +} + +extern "C" double source_bank_E(int64_t i) +{ + return simulation::source_bank[i-1].E; +} + +extern "C" double source_bank_wgt(int64_t i) +{ + return simulation::source_bank[i-1].wgt; +} + +extern "C" void source_bank_set_wgt(int64_t i, double wgt) +{ + simulation::source_bank[i-1].wgt = wgt; +} + +} // namespace openmc diff --git a/src/bank_header.F90 b/src/bank_header.F90 index c39bf8b784..992d357fb6 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -2,8 +2,6 @@ module bank_header use, intrinsic :: ISO_C_BINDING - use error, only: E_ALLOCATE, set_errmsg - implicit none !=============================================================================== @@ -21,70 +19,58 @@ module bank_header integer(C_INT) :: particle ! particle type (neutron, photon, etc.) end type Bank - ! Source and fission bank - type(Bank), allocatable, target :: source_bank(:) - type(Bank), allocatable, target :: fission_bank(:) -#ifdef _OPENMP - type(Bank), allocatable, target :: master_fission_bank(:) -#endif - integer(C_INT64_T), bind(C) :: n_bank ! # of sites in fission bank +!$omp threadprivate(n_bank) -!$omp threadprivate(fission_bank, n_bank) + interface + function openmc_fission_bank(ptr, n) result(err) bind(C) + import C_PTR, C_INT64_T, C_INT + type(C_PTR), intent(out) :: ptr + integer(C_INT64_T), intent(out) :: n + integer(C_INT) :: err + end function -contains + function fission_bank_delayed_group(i) result(g) bind(C) + import C_INT64_T, C_INT + integer(C_INT64_T), value :: i + integer(C_INT) :: g + end function -!=============================================================================== -! FREE_MEMORY_BANK deallocates global arrays defined in this module -!=============================================================================== + function fission_bank_E(i) result(E) bind(C, name='fission_bank_E') + import C_INT64_T, C_DOUBLE + integer(C_INT64_T), value :: i + real(C_DOUBLE) :: E + end function - subroutine free_memory_bank() + function fission_bank_wgt(i) result(wgt) bind(C) + import C_INT64_T, C_DOUBLE + integer(C_INT64_T), value :: i + real(C_DOUBLE) :: wgt + end function - ! Deallocate fission and source bank and entropy -!$omp parallel - if (allocated(fission_bank)) deallocate(fission_bank) -!$omp end parallel -#ifdef _OPENMP - if (allocated(master_fission_bank)) deallocate(master_fission_bank) -#endif - if (allocated(source_bank)) deallocate(source_bank) + subroutine source_bank_xyz(i, xyz) bind(C) + import C_INT64_T, C_DOUBLE + integer(C_INT64_T), value :: i + real(C_DOUBLE), intent(in) :: xyz(*) + end subroutine - end subroutine free_memory_bank + function source_bank_E(i) result(E) bind(C, name='source_bank_E') + import C_INT64_T, C_DOUBLE + integer(C_INT64_T), value :: i + real(C_DOUBLE) :: E + end function -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== + function source_bank_wgt(i) result(wgt) bind(C) + import C_INT64_T, C_DOUBLE + integer(C_INT64_T), value :: i + real(C_DOUBLE) :: wgt + end function - function openmc_source_bank(ptr, n) result(err) bind(C) - ! Return a pointer to the source bank - type(C_PTR), intent(out) :: ptr - integer(C_INT64_T), intent(out) :: n - integer(C_INT) :: err - - if (.not. allocated(source_bank)) then - err = E_ALLOCATE - call set_errmsg("Source bank has not been allocated.") - else - err = 0 - ptr = C_LOC(source_bank) - n = size(source_bank) - end if - end function openmc_source_bank - - function openmc_fission_bank(ptr, n) result(err) bind(C) - ! Return a pointer to the source bank - type(C_PTR), intent(out) :: ptr - integer(C_INT64_T), intent(out) :: n - integer(C_INT) :: err - - if (.not. allocated(fission_bank)) then - err = E_ALLOCATE - call set_errmsg("Fission bank has not been allocated.") - else - err = 0 - ptr = C_LOC(fission_bank) - n = size(fission_bank) - end if - end function openmc_fission_bank + subroutine source_bank_set_wgt(i, wgt) bind(C) + import C_INT64_T, C_DOUBLE + integer(C_INT64_T), value :: i + real(C_DOUBLE), value :: wgt + end subroutine + end interface end module bank_header diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 681d08c5c1..4b97b5e327 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -218,7 +218,7 @@ contains subroutine cmfd_reweight(new_weights) use algorithm, only: binary_search - use bank_header, only: source_bank + use bank_header use constants, only: ZERO, ONE use error, only: warning, fatal_error use message_passing @@ -230,13 +230,14 @@ contains integer :: ny ! maximum number of cells in y direction integer :: nz ! maximum number of cells in z direction integer(C_INT) :: ng ! maximum number of energy groups - integer :: i ! iteration counter + integer(8) :: i ! iteration counter integer :: g ! index for group integer :: ijk(3) ! spatial bin location integer :: e_bin ! energy bin of source particle integer :: mesh_bin ! mesh bin of soruce particle integer :: n_groups ! number of energy groups real(8) :: norm ! normalization factor + real(C_DOUBLE) :: xyz(3) logical(C_BOOL) :: outside ! any source sites outside mesh logical :: in_mesh ! source site is inside mesh @@ -309,21 +310,22 @@ contains end if ! begin loop over source bank - do i = 1, int(work,4) + do i = 1, work ! Determine spatial bin - call cmfd_mesh % get_indices(source_bank(i) % xyz, ijk, in_mesh) + call source_bank_xyz(i, xyz) + call cmfd_mesh % get_indices(xyz, ijk, in_mesh) ! Determine energy bin n_groups = size(cmfd % egrid) - 1 - if (source_bank(i) % E < cmfd % egrid(1)) then + if (source_bank_E(i) < cmfd % egrid(1)) then e_bin = 1 if (master) call warning('Source pt below energy grid') - elseif (source_bank(i) % E > cmfd % egrid(n_groups + 1)) then + elseif (source_bank_E(i) > cmfd % egrid(n_groups + 1)) then e_bin = n_groups if (master) call warning('Source pt above energy grid') else - e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank(i) % E) + e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank_E(i)) end if ! Reverese energy bin (lowest grp is highest energy bin) @@ -335,8 +337,8 @@ contains end if ! Reweight particle - source_bank(i) % wgt = source_bank(i) % wgt * & - cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3)) + call source_bank_set_wgt(i, source_bank_wgt(i) * & + cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3))) end do end subroutine cmfd_reweight diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp index 0af584c1b0..28561b00a7 100644 --- a/src/cmfd_execute.cpp +++ b/src/cmfd_execute.cpp @@ -5,6 +5,7 @@ #include "xtensor/xarray.hpp" #include "xtensor/xio.hpp" +#include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/mesh.h" #include "openmc/message_passing.h" @@ -18,14 +19,10 @@ extern "C" void cmfd_populate_sourcecounts(int n_energy, const double* energies, double* source_counts, bool* outside) { - // Get pointer to source bank - Bank* source_bank; - int64_t n; - openmc_source_bank(&source_bank, &n); - // Get source counts in each mesh bin / energy bin auto& m = model::meshes.at(settings::index_cmfd_mesh); - xt::xarray counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside); + xt::xarray counts = m->count_sites(simulation::work, + simulation::source_bank.data(), n_energy, energies, outside); // Copy data from the xarray into the source counts array std::copy(counts.begin(), counts.end(), source_counts); diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index bdbbd1c2bc..4de3ae8a67 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -2,9 +2,6 @@ module eigenvalue use, intrinsic :: ISO_C_BINDING - use bank_header - use simulation_header - implicit none interface @@ -15,49 +12,4 @@ module eigenvalue end function end interface -contains - -#ifdef _OPENMP -!=============================================================================== -! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission -! bank that can be sampled. Note that this operation is necessarily sequential -! to preserve the order of the bank when using varying numbers of threads. -!=============================================================================== - - subroutine join_bank_from_threads() bind(C) - - integer(8) :: total ! total number of fission bank sites - integer :: i ! loop index for threads - - ! Initialize the total number of fission bank sites - total = 0 - -!$omp parallel - - ! Copy thread fission bank sites to one shared copy -!$omp do ordered schedule(static) - do i = 1, n_threads -!$omp ordered - master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank) - total = total + n_bank -!$omp end ordered - end do -!$omp end do - - ! Make sure all threads have made it to this point -!$omp barrier - - ! Now copy the shared fission bank sites back to the master thread's copy. - if (thread_id == 0) then - n_bank = total - fission_bank(1:n_bank) = master_fission_bank(1:n_bank) - else - n_bank = 0 - end if - -!$omp end parallel - - end subroutine join_bank_from_threads -#endif - end module eigenvalue diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 64e0c4e420..b3b8b850c1 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -5,6 +5,7 @@ #include "xtensor/xtensor.hpp" #include "xtensor/xview.hpp" +#include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/error.h" @@ -69,13 +70,6 @@ void synchronize_bank() { simulation::time_bank.start(); - // Get pointers to source/fission bank - Bank* source_bank; - Bank* fission_bank; - int64_t n; - openmc_source_bank(&source_bank, &n); - openmc_fission_bank(&fission_bank, &n); - // In order to properly understand the fission bank algorithm, you need to // think of the fission and source bank as being one global array divided // over multiple processors. At the start, each processor has a random amount @@ -147,14 +141,14 @@ void synchronize_bank() // and the remaining 100 would be randomly sampled. if (total < settings::n_particles) { for (int64_t j = 1; j <= settings::n_particles / total; ++j) { - temp_sites[index_temp] = fission_bank[i]; + temp_sites[index_temp] = simulation::fission_bank[i]; ++index_temp; } } // Randomly sample sites needed if (prn() < p_sample) { - temp_sites[index_temp] = fission_bank[i]; + temp_sites[index_temp] = simulation::fission_bank[i]; ++index_temp; } } @@ -195,7 +189,8 @@ void synchronize_bank() // fission bank sites_needed = settings::n_particles - finish; for (int i = 0; i < sites_needed; ++i) { - temp_sites[index_temp] = fission_bank[simulation::n_bank - sites_needed + i]; + int i_bank = simulation::n_bank - sites_needed + i; + temp_sites[index_temp] = simulation::fission_bank[i_bank]; ++index_temp; } } @@ -274,7 +269,7 @@ void synchronize_bank() // asynchronous receive for the source sites requests.emplace_back(); - MPI_Irecv(&source_bank[index_local], static_cast(n), mpi::bank, + MPI_Irecv(&simulation::source_bank[index_local], static_cast(n), mpi::bank, neighbor, neighbor, mpi::intracomm, &requests.back()); } else { @@ -283,7 +278,7 @@ void synchronize_bank() index_temp = start - bank_position[mpi::rank]; std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n], - &source_bank[index_local]); + &simulation::source_bank[index_local]); } // Increment all indices @@ -300,7 +295,8 @@ void synchronize_bank() MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE); #else - std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, source_bank); + std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, + simulation::source_bank.begin()); #endif simulation::time_bank_sendrecv.stop(); @@ -347,6 +343,46 @@ void calculate_average_keff() } } +#ifdef _OPENMP +void join_bank_from_threads() +{ + // Initialize the total number of fission bank sites + int64_t total = 0; + +#pragma omp parallel + { + // Copy thread fission bank sites to one shared copy +#pragma omp for ordered schedule(static) + for (int i = 0; i < simulation::n_threads; ++i) { +#pragma omp ordered + { + std::copy( + &simulation::fission_bank[0], + &simulation::fission_bank[0] + simulation::n_bank, + &simulation::master_fission_bank[total] + ); + total += simulation::n_bank; + } + } + + // Make sure all threads have made it to this point +#pragma omp barrier + + // Now copy the shared fission bank sites back to the master thread's copy. + if (simulation::thread_id == 0) { + simulation::n_bank = total; + std::copy( + &simulation::master_fission_bank[0], + &simulation::master_fission_bank[0] + simulation::n_bank, + &simulation::fission_bank[0] + ); + } else { + simulation::n_bank = 0; + } + } +} +#endif + int openmc_get_keff(double* k_combined) { k_combined[0] = 0.0; @@ -499,15 +535,10 @@ void shannon_entropy() // Get pointer to entropy mesh auto& m = model::meshes[settings::index_entropy_mesh]; - // Get pointer to fission bank - Bank* fission_bank; - int64_t n; - openmc_fission_bank(&fission_bank, &n); - // Get source weight in each mesh bin bool sites_outside; - xt::xtensor p = m->count_sites( - simulation::n_bank, fission_bank, 0, nullptr, &sites_outside); + xt::xtensor p = m->count_sites(simulation::n_bank, + simulation::fission_bank.data(), 0, nullptr, &sites_outside); // display warning message if there were sites outside entropy box if (sites_outside) { @@ -544,15 +575,10 @@ void ufs_count_sites() s = m->volume_frac_; } else { - // Get pointer to source bank - Bank* source_bank; - int64_t n; - openmc_source_bank(&source_bank, &n); - // count number of source sites in each ufs mesh cell bool sites_outside; - simulation::source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr, - &sites_outside); + simulation::source_frac = m->count_sites(simulation::work, + simulation::source_bank.data(), 0, nullptr, &sites_outside); // Check for sites outside of the mesh if (mpi::master && sites_outside) { @@ -572,7 +598,7 @@ void ufs_count_sites() // Since the total starting weight is not equal to n_particles, we need to // renormalize the weight of the source sites for (int i = 0; i < simulation::work; ++i) { - source_bank[i].wgt *= settings::n_particles / total; + simulation::source_bank[i].wgt *= settings::n_particles / total; } } } diff --git a/src/endf.cpp b/src/endf.cpp index 8ffb75345d..db1332591e 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -45,6 +45,37 @@ bool is_fission(int mt) return mt == 18 || mt == 19 || mt == 20 || mt == 21 || mt == 38; } +bool is_disappearance(int mt) +{ + if (mt >= N_DISAPPEAR && mt <= N_DA) { + return true; + } else if (mt >= N_P0 && mt <= N_AC) { + return true; + } else if (mt == N_TA || mt == N_DT || mt == N_P3HE || mt == N_D3HE + || mt == N_3HEA || mt == N_3P) { + return true; + } else { + return false; + } +} + +bool is_inelastic_scatter(int mt) +{ + if (mt < 100) { + if (is_fission(mt)) { + return false; + } else { + return mt >= MISC && mt != 27; + } + } else if (mt <= 200) { + return !is_disappearance(mt); + } else if (mt >= N_2N0 && mt <= N_2NC) { + return true; + } else { + return false; + } +} + //============================================================================== // Polynomial implementation //============================================================================== diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1da1bf9df6..dab02f52bc 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2113,9 +2113,6 @@ contains call close_group(group_id) call file_close(file_id) - ! Assign resonant scattering data - if (res_scat_on) call nuclides(i_nuclide) % assign_0K_elastic_scattering() - ! Determine if minimum/maximum energy for this nuclide is greater/less ! than the previous if (size(nuclides(i_nuclide) % grid) >= 1) then @@ -2297,7 +2294,6 @@ contains logical :: file_exists ! Does multipole library exist? character(7) :: readable ! Is multipole library readable? character(MAX_FILE_LEN) :: filename ! Path to multipole xs library - character(kind=C_CHAR), pointer :: string(:) integer(HID_T) :: file_id integer(HID_T) :: group_id diff --git a/src/material_header.F90 b/src/material_header.F90 index 3bcf96db8b..e041f52f4a 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -1041,4 +1041,21 @@ contains end subroutine bremsstrahlung_init +!=============================================================================== +! Fortran compatibility +!=============================================================================== + + function material_isotropic(i_material, i_nuc_mat) result(iso) bind(C) + integer(C_INT), value :: i_material + integer(C_INT), value :: i_nuc_mat + logical(C_BOOL) :: iso + + iso = .false. + associate (mat => materials(i_material)) + if (mat % has_isotropic_nuclides) then + iso = mat % p0(i_nuc_mat) + end if + end associate + end function + end module material_header diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 18e7e5a7e1..ea1e213d57 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1,5 +1,17 @@ #include "openmc/nuclide.h" +#include "openmc/container_util.h" +#include "openmc/endf.h" +#include "openmc/error.h" +#include "openmc/hdf5_interface.h" +#include "openmc/message_passing.h" +#include "openmc/search.h" +#include "openmc/settings.h" +#include "openmc/string_utils.h" + +#include // for sort +#include // for to_string, stoi + namespace openmc { //============================================================================== @@ -7,12 +19,348 @@ namespace openmc { //============================================================================== namespace data { - std::array energy_min {0.0, 0.0}; std::array energy_max {INFTY, INFTY}; - +std::vector> nuclides; } // namespace data +namespace simulation { +NuclideMicroXS* micro_xs; +MaterialMacroXS material_xs; +} // namespace simulation + +//============================================================================== +// Nuclide implementation +//============================================================================== + +Nuclide::Nuclide(hid_t group, const double* temperature, int n) +{ + // Get name of nuclide from group, removing leading '/' + name_ = object_name(group).substr(1); + + read_attribute(group, "Z", Z_); + read_attribute(group, "A", A_); + read_attribute(group, "metastable", metastable_); + read_attribute(group, "atomic_weight_ratio", awr_); + + // Determine temperatures available + hid_t kT_group = open_group(group, "kTs"); + auto dset_names = dataset_names(kT_group); + std::vector temps_available; + for (const auto& name : dset_names) { + double T; + read_dataset(kT_group, name.c_str(), T); + temps_available.push_back(T / K_BOLTZMANN); + } + std::sort(temps_available.begin(), temps_available.end()); + + // If only one temperature is available, revert to nearest temperature + if (temps_available.size() == 1 && settings::temperature_method == TEMPERATURE_INTERPOLATION) { + if (mpi::master) { + warning("Cross sections for " + name_ + " are only available at one " + "temperature. Reverting to nearest temperature method."); + } + settings::temperature_method = TEMPERATURE_NEAREST; + } + + // Determine actual temperatures to read -- start by checking whether a + // temperature range was given, in which case all temperatures in the range + // are loaded irrespective of what temperatures actually appear in the model + std::vector temps_to_read; + double T_min = n > 0 ? settings::temperature_range[0] : 0.0; + double T_max = n > 0 ? settings::temperature_range[1] : INFTY; + if (T_max > 0.0) { + for (auto T : temps_available) { + if (T_min <= T && T <= T_max) { + temps_to_read.push_back(std::round(T)); + } + } + } + + switch (settings::temperature_method) { + case TEMPERATURE_NEAREST: + // Find nearest temperatures + for (int i = 0; i < n; ++i) { + double T_desired = temperature[i]; + + // Determine closest temperature + double min_delta_T = INFTY; + double T_actual; + for (auto T : temps_available) { + double delta_T = std::abs(T - T_desired); + if (delta_T < min_delta_T) { + T_actual = T; + min_delta_T = delta_T; + } + } + + if (std::abs(T_actual - T_desired) < settings::temperature_tolerance) { + if (!contains(temps_to_read, std::round(T_actual))) { + temps_to_read.push_back(std::round(T_actual)); + + // Write warning for resonance scattering data if 0K is not available + if (std::abs(T_actual - T_desired) > 0 && T_desired == 0 && mpi::master) { + warning(name_ + " does not contain 0K data needed for resonance " + "scattering options selected. Using data at " + std::to_string(T_actual) + + " K instead."); + } + } + } else { + fatal_error("Nuclear data library does not contain cross sections for " + + name_ + " at or near " + std::to_string(T_desired) + " K."); + } + } + break; + + case TEMPERATURE_INTERPOLATION: + // If temperature interpolation or multipole is selected, get a list of + // bounding temperatures for each actual temperature present in the model + for (int i = 0; i < n; ++i) { + double T_desired = temperature[i]; + + bool found_pair = false; + for (int j = 0; j < temps_available.size() - 1; ++j) { + if (temps_available[j] <= T_desired && T_desired < temps_available[j + 1]) { + int T_j = std::round(temps_available[j]); + int T_j1 = std::round(temps_available[j+1]); + if (!contains(temps_to_read, T_j)) { + temps_to_read.push_back(T_j); + } + if (!contains(temps_to_read, T_j1)) { + temps_to_read.push_back(T_j1); + } + found_pair = true; + } + } + + if (!found_pair) { + fatal_error("Nuclear data library does not contain cross sections for " + + name_ +" at temperatures that bound " + std::to_string(T_desired) + " K."); + } + } + break; + } + + // Sort temperatures to read + std::sort(temps_to_read.begin(), temps_to_read.end()); + + hid_t energy_group = open_group(group, "energy"); + for (const auto& T : temps_to_read) { + std::string dset {std::to_string(T) + "K"}; + + // Determine exact kT values + double kT; + read_dataset(kT_group, dset.c_str(), kT); + kTs_.push_back(kT); + + // Read energy grid + grid_.emplace_back(); + read_dataset(energy_group, dset.c_str(), grid_.back().energy); + } + close_group(kT_group); + + // Check for 0K energy grid + if (object_exists(energy_group, "0K")) { + read_dataset(energy_group, "0K", energy_0K_); + } + close_group(energy_group); + + // Read reactions + hid_t rxs_group = open_group(group, "reactions"); + for (auto name : group_names(rxs_group)) { + if (starts_with(name, "reaction_")) { + hid_t rx_group = open_group(rxs_group, name.c_str()); + reactions_.push_back(std::make_unique(rx_group, temps_to_read)); + + // Check for 0K elastic scattering + const auto& rx = reactions_.back(); + if (rx->mt_ == ELASTIC) { + if (object_exists(rx_group, "0K")) { + hid_t temp_group = open_group(rx_group, "0K"); + read_dataset(temp_group, "xs", elastic_0K_); + close_group(temp_group); + } + } + close_group(rx_group); + + // Determine reaction indices for inelastic scattering reactions + if (is_inelastic_scatter(rx->mt_) && !rx->redundant_) { + index_inelastic_scatter_.push_back(reactions_.size() - 1); + } + } + } + close_group(rxs_group); + + // Check for nu-total + if (object_exists(group, "total_nu")) { + // Read total nu data + hid_t nu_group = open_group(group, "total_nu"); + hid_t nu_dset = open_dataset(nu_group, "yield"); + std::string func_type; + read_attribute(nu_dset, "type", func_type); + if (func_type == "Tabulated1D") { + total_nu_ = std::make_unique(nu_dset); + } else if (func_type == "Polynomial") { + total_nu_ = std::make_unique(nu_dset); + } + close_dataset(nu_dset); + close_group(nu_group); + } + + this->create_derived(); +} + +void Nuclide::create_derived() +{ + for (int i = 0; i < reactions_.size(); ++i) { + const auto& rx {reactions_[i]}; + + for (int t = 0; t < kTs_.size(); ++t) { + // Skip redundant reactions + if (rx->redundant_) continue; + + if (is_fission(rx->mt_)) { + fissionable_ = true; + + // Keep track of fission reactions + if (t == 0) { + fission_rx_.push_back(rx.get()); + if (rx->mt_ == N_F) has_partial_fission_ = true; + } + } + } + } + + // Determine number of delayed neutron precursors + if (fissionable_) { + for (const auto& product : fission_rx_[0]->products_) { + if (product.emission_mode_ == EmissionMode::delayed) { + ++n_precursor_; + } + } + } + + if (settings::res_scat_on) { + // Determine if this nuclide should be treated as a resonant scatterer + if (!settings::res_scat_nuclides.empty()) { + // If resonant nuclides were specified, check the list explicitly + for (const auto& name : settings::res_scat_nuclides) { + if (name_ == name) { + resonant_ = true; + + // Make sure nuclide has 0K data + if (energy_0K_.empty()) { + fatal_error("Cannot treat " + name_ + " as a resonant scatterer " + "because 0 K elastic scattering data is not present."); + } + break; + } + } + } else { + // Otherwise, assume that any that have 0 K elastic scattering data are + // resonant + resonant_ = !energy_0K_.empty(); + } + + if (resonant_) { + // Build CDF for 0K elastic scattering + double xs_cdf_sum = 0.0; + xs_cdf_.resize(energy_0K_.size()); + xs_cdf_[0] = 0.0; + + const auto& E = energy_0K_; + auto& xs = elastic_0K_; + for (int i = 0; i < E.size() - 1; ++i) { + // Negative cross sections result in a CDF that is not monotonically + // increasing. Set all negative xs values to zero. + if (xs[i] < 0.0) xs[i] = 0.0; + + // build xs cdf + xs_cdf_sum += (std::sqrt(E[i])*xs[i] + std::sqrt(E[i+1])*xs[i+1]) + / 2.0 * (E[i+1] - E[i]); + xs_cdf_[i] = xs_cdf_sum; + } + } + } +} + +double Nuclide::nu(double E, EmissionMode mode, int group) const +{ + if (!fissionable_) return 0.0; + + switch (mode) { + case EmissionMode::prompt: + return (*fission_rx_[0]->products_[0].yield_)(E); + case EmissionMode::delayed: + if (n_precursor_ > 0) { + auto rx = fission_rx_[0]; + if (group >= 1 && group < rx->products_.size()) { + // If delayed group specified, determine yield immediately + return (*rx->products_[group].yield_)(E); + } else { + double nu {0.0}; + + for (int i = 1; i < rx->products_.size(); ++i) { + // Skip any non-neutron products + const auto& product = rx->products_[i]; + if (product.particle_ != ParticleType::neutron) continue; + + // Evaluate yield + if (product.emission_mode_ == EmissionMode::delayed) { + nu += (*product.yield_)(E); + } + } + return nu; + } + } else { + return 0.0; + } + case EmissionMode::total: + if (total_nu_) { + return (*total_nu_)(E); + } else { + return (*fission_rx_[0]->products_[0].yield_)(E); + } + } +} + +void Nuclide::calculate_elastic_xs(int i_nuclide) const +{ + // Get temperature index, grid index, and interpolation factor + auto& micro = simulation::micro_xs[i_nuclide]; + int i_temp = micro.index_temp - 1; + int i_grid = micro.index_grid - 1; + double f = micro.interp_factor; + + if (i_temp >= 0) { + const auto& xs = reactions_[0]->xs_[i_temp].value; + micro.elastic = (1.0 - f)*xs[i_grid] + f*xs[i_grid + 1]; + } +} + +double Nuclide::elastic_xs_0K(double E) const +{ + // Determine index on nuclide energy grid + int i_grid; + if (E < energy_0K_.front()) { + i_grid = 0; + } else if (E > energy_0K_.back()) { + i_grid = energy_0K_.size() - 2; + } else { + i_grid = lower_bound_index(energy_0K_.begin(), energy_0K_.end(), E); + } + + // check for rare case where two energy points are the same + if (energy_0K_[i_grid] == energy_0K_[i_grid+1]) ++i_grid; + + // calculate interpolation factor + double f = (E - energy_0K_[i_grid]) / + (energy_0K_[i_grid + 1] - energy_0K_[i_grid]); + + // Calculate microscopic nuclide elastic cross section + return (1.0 - f)*elastic_0K_[i_grid] + f*elastic_0K_[i_grid + 1]; +} + //============================================================================== // Fortran compatibility functions //============================================================================== @@ -24,4 +372,27 @@ set_particle_energy_bounds(int particle, double E_min, double E_max) data::energy_max[particle - 1] = E_max; } +extern "C" Nuclide* nuclide_from_hdf5_c(hid_t group, const double* temperature, int n) +{ + data::nuclides.push_back(std::make_unique(group, temperature, n)); + return data::nuclides.back().get(); +} + +extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx) +{ + return nuc->reactions_[i_rx-1].get(); +} + +extern "C" void nuclides_clear() { data::nuclides.clear(); } + + +extern "C" NuclideMicroXS* micro_xs_ptr(); +void set_micro_xs() +{ +#pragma omp parallel + { + simulation::micro_xs = micro_xs_ptr(); + } +} + } // namespace openmc diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 921dcd3bd2..becb865240 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -71,15 +71,7 @@ module nuclide_header ! Microscopic cross sections type(SumXS), allocatable :: xs(:) - ! Resonance scattering info - logical :: resonant = .false. ! resonant scatterer? - real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs - real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section - real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section - ! Fission information - logical :: has_partial_fission = .false. ! nuclide has partial fission reactions? - integer :: n_fission = 0 ! # of fission reactions integer :: n_precursor = 0 ! # of delayed neutron precursors integer, allocatable :: index_fission(:) ! indices in reactions class(Function1D), allocatable :: total_nu @@ -95,7 +87,6 @@ module nuclide_header ! Reactions type(Reaction), allocatable :: reactions(:) - integer, allocatable :: index_inelastic_scatter(:) ! Array that maps MT values to index in reactions; used at tally-time. Note ! that ENDF-102 does not assign any MT values above 891. @@ -105,8 +96,9 @@ module nuclide_header class(Function1D), allocatable :: fission_q_prompt ! fragments and prompt neutrons, gammas class(Function1D), allocatable :: fission_q_recov ! fragments, neutrons, gammas, betas + type(C_PTR) :: ptr + contains - procedure :: assign_0K_elastic_scattering procedure :: clear => nuclide_clear procedure :: from_hdf5 => nuclide_from_hdf5 procedure :: init_grid => nuclide_init_grid @@ -184,7 +176,7 @@ module nuclide_header ! Cross section caches type(NuclideMicroXS), allocatable, target :: micro_xs(:) ! Cache for each nuclide - type(MaterialMacroXS) :: material_xs ! Cache for current material + type(MaterialMacroXS), bind(C) :: material_xs ! Cache for current material !$omp threadprivate(micro_xs, material_xs) ! Minimum/maximum energies @@ -231,78 +223,10 @@ contains b = library_present_c(type, to_c_string(name)) end function -!=============================================================================== -! ASSIGN_0K_ELASTIC_SCATTERING -!=============================================================================== - - subroutine assign_0K_elastic_scattering(this) - class(Nuclide), intent(inout) :: this - - integer :: i - real(8) :: xs_cdf_sum - - interface - function res_scat_nuclides_empty() result(empty) bind(C) - import C_BOOL - logical(C_BOOL) :: empty - end function - - function res_scat_nuclides_size() result(n) bind(C) - import C_INT - integer(C_INT) :: n - end function - - function res_scat_nuclides_cmp(i, name) result(b) bind(C) - import C_INT, C_CHAR, C_BOOL - integer(C_INT), value :: i - character(kind=C_CHAR), intent(in) :: name(*) - logical(C_BOOL) :: b - end function - end interface - - this % resonant = .false. - if (.not. res_scat_nuclides_empty()) then - ! If resonant nuclides were specified, check the list explicitly - do i = 1, res_scat_nuclides_size() - if (res_scat_nuclides_cmp(i, to_c_string(this % name))) then - this % resonant = .true. - - ! Make sure nuclide has 0K data - if (.not. allocated(this % energy_0K)) then - call fatal_error("Cannot treat " // trim(this % name) // " as a & - &resonant scatterer because 0 K elastic scattering data is & - ¬ present.") - end if - - exit - end if - end do - else - ! Otherwise, assume that any that have 0 K elastic scattering data are - ! resonant - this % resonant = allocated(this % energy_0K) - end if - - if (this % resonant) then - ! Build CDF for 0K elastic scattering - xs_cdf_sum = ZERO - allocate(this % xs_cdf(0:size(this % energy_0K))) - this % xs_cdf(0) = ZERO - - associate (E => this % energy_0K, xs => this % elastic_0K) - do i = 1, size(E) - 1 - ! Negative cross sections result in a CDF that is not monotonically - ! increasing. Set all negative xs values to zero. - if (xs(i) < ZERO) xs(i) = ZERO - - ! build xs cdf - xs_cdf_sum = xs_cdf_sum + (sqrt(E(i))*xs(i) + sqrt(E(i+1))*xs(i+1))& - / TWO * (E(i+1) - E(i)) - this % xs_cdf(i) = xs_cdf_sum - end do - end associate - end if - end subroutine assign_0K_elastic_scattering + function micro_xs_ptr() result(ptr) bind(C) + type(C_PTR) :: ptr + ptr = C_LOC(micro_xs(1)) + end function !=============================================================================== ! NUCLIDE_CLEAR resets and deallocates data in Nuclide @@ -310,19 +234,6 @@ contains subroutine nuclide_clear(this) class(Nuclide), intent(inout) :: this ! The Nuclide object to clear - integer :: i - - interface - subroutine reaction_delete(rx) bind(C) - import C_PTR - type(C_PTR), value :: rx - end subroutine reaction_delete - end interface - - do i = 1, size(this % reactions) - call reaction_delete(this % reactions(i) % ptr) - end do - deallocate(this % reactions) if (associated(this % multipole)) deallocate(this % multipole) @@ -332,7 +243,7 @@ contains minmax, master, i_nuclide) class(Nuclide), intent(inout) :: this integer(HID_T), intent(in) :: group_id - type(VectorReal), intent(in) :: temperature ! list of desired temperatures + type(VectorReal), intent(in), target :: temperature ! list of desired temperatures integer, intent(inout) :: method real(8), intent(in) :: tolerance real(8), intent(in) :: minmax(2) ! range of temperatures @@ -347,7 +258,6 @@ contains integer(HID_T) :: kT_group integer(HID_T) :: rxs_group integer(HID_T) :: rx_group - integer(HID_T) :: xs, temp_group integer(HID_T) :: total_nu integer(HID_T) :: fer_group ! fission_energy_release group integer(HID_T) :: fer_dset @@ -361,7 +271,24 @@ contains real(8) :: temp_actual type(VectorInt) :: MTs type(VectorInt) :: temps_to_read - type(VectorInt) :: index_inelastic_scatter + + interface + function nuclide_from_hdf5_c(group, temperature, n) result(ptr) bind(C) + import HID_T, C_DOUBLE, C_INT, C_PTR + integer(HID_T), value :: group + type(C_PTR), value :: temperature + integer(C_INT), value :: n + type(C_PTR) :: ptr + end function + end interface + + ! Read data on C++ side + if (temperature % size() > 0) then + this % ptr = nuclide_from_hdf5_c(group_id, C_LOC(temperature % data(1)), & + temperature % size()) + else + this % ptr = nuclide_from_hdf5_c(group_id, C_NULL_PTR, 0) + end if ! Get name of nuclide from group this % name = get_name(group_id) @@ -487,15 +414,6 @@ contains call read_dataset(this % grid(i) % energy, energy_dset) call close_dataset(energy_dset) end do - - ! Check for 0K energy grid - if (object_exists(energy_group, '0K')) then - energy_dset = open_dataset(energy_group, '0K') - call get_shape(energy_dset, dims) - allocate(this % energy_0K(int(dims(1), 4))) - call read_dataset(this % energy_0K, energy_dset) - call close_dataset(energy_dset) - end if call close_group(energy_group) ! Get MT values based on group names @@ -513,36 +431,13 @@ contains rx_group = open_group(rxs_group, 'reaction_' // trim(& zero_padded(MTs % data(i), 3))) - call this % reactions(i) % from_hdf5(rx_group, temps_to_read) - - ! Check for 0K elastic scattering - if (this % reactions(i) % MT == 2) then - if (object_exists(rx_group, '0K')) then - temp_group = open_group(rx_group, '0K') - xs = open_dataset(temp_group, 'xs') - call get_shape(xs, dims) - allocate(this % elastic_0K(int(dims(1), 4))) - call read_dataset(this % elastic_0K, xs) - call close_dataset(xs) - call close_group(temp_group) - end if - end if - - ! Add the reaction index to the scattering array if this is an inelastic - ! scatter reaction - if (is_inelastic_scatter(MTs % data(i))) then - call index_inelastic_scatter % push_back(i) - end if + ! Set pointer for each reaction + call this % reactions(i) % init(this % ptr, i) call close_group(rx_group) end do call close_group(rxs_group) - ! Recast to a regular array to save space - allocate(this % index_inelastic_scatter(index_inelastic_scatter % size())) - this % index_inelastic_scatter = & - index_inelastic_scatter % data(1: index_inelastic_scatter % size()) - ! Read unresolved resonance probability tables if present if (object_exists(group_id, 'urr')) then this % urr_present = .true. @@ -718,7 +613,6 @@ contains allocate(this % index_fission(1)) elseif (rx % MT == N_F) then allocate(this % index_fission(PARTIAL_FISSION_MAX)) - this % has_partial_fission = .true. end if end if @@ -738,7 +632,6 @@ contains if (t == 1) then i_fission = i_fission + 1 this % index_fission(i_fission) = i - this % n_fission = this % n_fission + 1 end if end if ! fission end do ! temperature @@ -1353,45 +1246,6 @@ contains end if end subroutine multipole_deriv_eval -!=============================================================================== -! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section -! for a given nuclide at the trial relative energy used in resonance scattering -!=============================================================================== - - pure function elastic_xs_0K(E, nuc) result(xs_out) - real(8), intent(in) :: E ! trial energy - type(Nuclide), intent(in) :: nuc ! target nuclide at temperature - real(8) :: xs_out ! 0K xs at trial energy - - integer :: i_grid ! index on nuclide energy grid - integer :: n_grid - real(8) :: f ! interp factor on nuclide energy grid - - ! Determine index on nuclide energy grid - n_grid = size(nuc % energy_0K) - if (E < nuc % energy_0K(1)) then - i_grid = 1 - elseif (E > nuc % energy_0K(n_grid)) then - i_grid = n_grid - 1 - else - i_grid = binary_search(nuc % energy_0K, n_grid, E) - end if - - ! check for rare case where two energy points are the same - if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then - i_grid = i_grid + 1 - end if - - ! calculate interpolation factor - f = (E - nuc % energy_0K(i_grid)) & - & / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid)) - - ! Calculate microscopic nuclide elastic cross section - xs_out = (ONE - f) * nuc % elastic_0K(i_grid) & - & + f * nuc % elastic_0K(i_grid + 1) - - end function elastic_xs_0K - !=============================================================================== ! CALCULATE_URR_XS determines cross sections in the unresolved resonance range ! from probability tables @@ -1573,6 +1427,9 @@ contains interface subroutine library_clear() bind(C) end subroutine + + subroutine nuclides_clear() bind(C) + end subroutine end interface ! Deallocate cross section data, listings, and cache @@ -1582,6 +1439,7 @@ contains call nuclides(i) % clear() end do deallocate(nuclides) + call nuclides_clear() end if n_nuclides = 0 @@ -1666,9 +1524,6 @@ contains call nuclide_dict % set(to_lower(name_), n) n_nuclides = n - ! Assign resonant scattering data - if (res_scat_on) call nuclides(n) % assign_0K_elastic_scattering() - ! Initialize nuclide grid call nuclides(n) % init_grid(energy_min(NEUTRON), & energy_max(NEUTRON), n_log_bins) @@ -1706,4 +1561,23 @@ contains end if end function openmc_nuclide_name + function nuclide_wmp_present(i_nuclide) result(b) bind(C) + integer(C_INT), value :: i_nuclide + logical(C_BOOL) :: b + b = nuclides(i_nuclide + 1) % mp_present + end function + + function nuclide_wmp_emin(i_nuclide) result(E) bind(C) + integer(C_INT), value :: i_nuclide + real(C_DOUBLE) :: E + E = nuclides(i_nuclide + 1) % multipole % E_min + end function + + function nuclide_wmp_emax(i_nuclide) result(E) bind(C) + integer(C_INT), value :: i_nuclide + real(C_DOUBLE) :: E + E = nuclides(i_nuclide + 1) % multipole % E_max + end function + + end module nuclide_header diff --git a/src/particle.cpp b/src/particle.cpp index f87ff30666..c8a7538014 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -3,6 +3,7 @@ #include #include +#include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/error.h" @@ -184,17 +185,12 @@ Particle::write_restart() const write_dataset(file_id, "id", id); write_dataset(file_id, "type", type); - // Get pointer to source bank - Bank* src; - int64_t n; - openmc_source_bank(&src, &n); - int64_t i = simulation::current_work; - write_dataset(file_id, "weight", src[i-1].wgt); - write_dataset(file_id, "energy", src[i-1].E); + write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt); + write_dataset(file_id, "energy", simulation::source_bank[i-1].E); hsize_t dims[] {3}; - write_double(file_id, 1, dims, "xyz", src[i-1].xyz, false); - write_double(file_id, 1, dims, "uvw", src[i-1].uvw, false); + write_double(file_id, 1, dims, "xyz", simulation::source_bank[i-1].xyz, false); + write_double(file_id, 1, dims, "uvw", simulation::source_bank[i-1].uvw, false); // Close file file_close(file_id); diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 0190c6e979..78a558d346 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -33,12 +33,20 @@ contains integer :: previous_run_mode type(Particle) :: p + interface + subroutine set_micro_xs() bind(C) + end subroutine + end interface + err = 0 ! Set verbosity high verbosity = 10 + !$omp parallel allocate(micro_xs(n_nuclides)) + !$omp end parallel + call set_micro_xs() ! Initialize the particle to be tracked call particle_initialize(p) diff --git a/src/photon.cpp b/src/photon.cpp new file mode 100644 index 0000000000..b36f29a5c3 --- /dev/null +++ b/src/photon.cpp @@ -0,0 +1,38 @@ +#include "openmc/photon.h" + +#include "openmc/hdf5_interface.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +std::vector elements; + +} // namespace data + +//============================================================================== +// PhotonInteraction implementation +//============================================================================== + +PhotonInteraction::PhotonInteraction(hid_t group) +{ + // Get name of nuclide from group, removing leading '/' + name_ = object_name(group).substr(1); + + read_attribute(group, "Z", Z_); +} + +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" void photon_from_hdf5_c(hid_t group) +{ + data::elements.emplace_back(group); +} + +} // namespace openmc diff --git a/src/photon_header.F90 b/src/photon_header.F90 index a284c7cc02..9da958a6f4 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -130,6 +130,16 @@ contains real(8), allocatable :: matrix(:,:) real(8), allocatable :: dcs(:,:) + interface + subroutine photon_from_hdf5_c(group) bind(C) + import HID_T + integer(HID_T), value :: group + end subroutine + end interface + + ! Read element data on C++ side + call photon_from_hdf5_c(group_id) + ! Get name of nuclide from group this % name = get_name(group_id) diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index ff631b5cf1..1a9e4edeac 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -540,9 +540,9 @@ contains ! THICK_TARGET_BREMSSTRAHLUNG !=============================================================================== - subroutine thick_target_bremsstrahlung(p, E_lost) + subroutine thick_target_bremsstrahlung(p, E_lost) bind(C) type(Particle), intent(inout) :: p - real(8), intent(inout) :: E_lost + real(C_DOUBLE), intent(out) :: E_lost integer :: i, j integer :: i_e, i_w diff --git a/src/physics.F90 b/src/physics.F90 index 29a8064951..0576e2a887 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,8 +1,6 @@ module physics - use algorithm, only: binary_search use constants - use endf, only: reaction_name use error, only: fatal_error, warning, write_message use material_header, only: Material, materials use math @@ -14,147 +12,30 @@ module physics atomic_relaxation, pair_production, & thick_target_bremsstrahlung use physics_common - use random_lcg, only: prn, advance_prn_seed, prn_set_stream - use reaction_header, only: Reaction + use random_lcg, only: prn use sab_header, only: sab_tables use settings use simulation_header - use string, only: to_str use tally_header implicit none + interface + subroutine collision(p) bind(C) + import Particle + type(Particle), intent(inout) :: p + end subroutine + end interface + contains -!=============================================================================== -! COLLISION samples a nuclide and reaction and then calls the appropriate -! routine for that reaction -!=============================================================================== - - subroutine collision(p) - - type(Particle), intent(inout) :: p - - ! Add to collision counter for particle - p % n_collision = p % n_collision + 1 - - ! Sample reaction for the material the particle is in - if (p % type == NEUTRON) then - call sample_neutron_reaction(p) - else if (p % type == PHOTON) then - call sample_photon_reaction(p) - else if (p % type == ELECTRON) then - call sample_electron_reaction(p) - else if (p % type == POSITRON) then - call sample_positron_reaction(p) - end if - - ! Kill particle if energy falls below cutoff - if (p % E < energy_cutoff(p % type)) then - p % alive = .false. - p % wgt = ZERO - p % last_wgt = ZERO - end if - - ! Display information about collision - if (verbosity >= 10 .or. trace) then - if (p % type == NEUTRON) then - call write_message(" " // trim(reaction_name(p % event_MT)) & - &// " with " // trim(adjustl(nuclides(p % event_nuclide) % name)) & - &// ". Energy = " // trim(to_str(p % E)) // " eV.") - else - call write_message(" " // trim(reaction_name(p % event_MT)) & - &// " with " // trim(adjustl(elements(p % event_nuclide) % name)) & - &// ". Energy = " // trim(to_str(p % E)) // " eV.") - end if - end if - - end subroutine collision - -!=============================================================================== -! SAMPLE_NEUTRON_REACTION samples a nuclide based on the macroscopic cross -! sections for each nuclide within a material and then samples a reaction for -! that nuclide and calls the appropriate routine to process the physics. Note -! that there is special logic when suvival biasing is turned on since fission -! and disappearance are treated implicitly. -!=============================================================================== - - subroutine sample_neutron_reaction(p) - - type(Particle), intent(inout) :: p - - integer :: i_nuclide ! index in nuclides array - integer :: i_nuc_mat ! index in material's nuclides array - integer :: i_reaction ! index in nuc % reactions array - type(Nuclide), pointer :: nuc - - call sample_nuclide(p, 'total ', i_nuclide, i_nuc_mat) - - ! Get pointer to table - nuc => nuclides(i_nuclide) - - ! Save which nuclide particle had collision with - p % event_nuclide = i_nuclide - - ! Create fission bank sites. Note that while a fission reaction is sampled, - ! it never actually "happens", i.e. the weight of the particle does not - ! change when sampling fission sites. The following block handles all - ! absorption (including fission) - - if (nuc % fissionable) then - if (run_mode == MODE_EIGENVALUE) then - call sample_fission(i_nuclide, p % E, i_reaction) - call create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank) - elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then - call sample_fission(i_nuclide, p % E, i_reaction) - call create_fission_sites(p, i_nuclide, i_reaction, & - p % secondary_bank, p % n_secondary) - end if - end if - - ! Create secondary photons - if (photon_transport) then - call prn_set_stream(STREAM_PHOTON) - call sample_secondary_photons(p, i_nuclide) - call prn_set_stream(STREAM_TRACKING) - end if - - ! If survival biasing is being used, the following subroutine adjusts the - ! weight of the particle. Otherwise, it checks to see if absorption occurs - - if (micro_xs(i_nuclide) % absorption > ZERO) then - call absorption(p, i_nuclide) - else - p % absorb_wgt = ZERO - end if - if (.not. p % alive) return - - ! Sample a scattering reaction and determine the secondary energy of the - ! exiting neutron - call scatter(p, i_nuclide, i_nuc_mat) - - ! Advance URR seed stream 'N' times after energy changes - if (p % E /= p % last_E) then - call prn_set_stream(STREAM_URR_PTABLE) - call advance_prn_seed(size(nuclides, kind=8)) - call prn_set_stream(STREAM_TRACKING) - end if - - ! Play russian roulette if survival biasing is turned on - if (survival_biasing) then - call russian_roulette(p) - if (.not. p % alive) return - end if - - end subroutine sample_neutron_reaction - !=============================================================================== ! SAMPLE_PHOTON_REACTION samples an element based on the macroscopic cross ! sections for each nuclide within a material and then samples a reaction for ! that element and calls the appropriate routine to process the physics. !=============================================================================== - subroutine sample_photon_reaction(p) + subroutine sample_photon_reaction(p) bind(C) type(Particle), intent(inout) :: p integer :: i_shell ! index in subshells @@ -325,130 +206,6 @@ contains end subroutine sample_photon_reaction -!=============================================================================== -! SAMPLE_ELECTRON_REACTION terminates the particle and either deposits all -! energy locally (electron_treatment = ELECTRON_LED) or creates secondary -! bremsstrahlung photons from electron deflections with charged particles -! (electron_treatment = ELECTRON_TTB). -!=============================================================================== - - subroutine sample_electron_reaction(p) - type(Particle), intent(inout) :: p - - real(8) :: E_lost ! energy lost to bremsstrahlung photons - - ! TODO: create reaction types - - if (electron_treatment == ELECTRON_TTB) then - call thick_target_bremsstrahlung(p, E_lost) - end if - - p % E = ZERO - p % alive = .false. - - end subroutine sample_electron_reaction - -!=============================================================================== -! SAMPLE_POSITRON_REACTION terminates the particle and either deposits all -! energy locally (electron_treatment = ELECTRON_LED) or creates secondary -! bremsstrahlung photons from electron deflections with charged particles -! (electron_treatment = ELECTRON_TTB). Two annihilation photons of energy -! MASS_ELECTRON_EV (0.511 MeV) are created and travel in opposite directions. -!=============================================================================== - - subroutine sample_positron_reaction(p) - type(Particle), intent(inout) :: p - - real(8) :: mu ! scattering cosine - real(8) :: phi ! azimuthal angle - real(8) :: uvw(3) ! new direction - - real(8) :: E_lost ! energy lost to bremsstrahlung photons - - ! TODO: create reaction types - - if (electron_treatment == ELECTRON_TTB) then - call thick_target_bremsstrahlung(p, E_lost) - end if - - ! Sample angle isotropically - mu = TWO*prn() - ONE - phi = TWO*PI*prn() - uvw(1) = mu - uvw(2) = sqrt(ONE - mu*mu)*cos(phi) - uvw(3) = sqrt(ONE - mu*mu)*sin(phi) - - ! Create annihilation photon pair traveling in opposite directions - call particle_create_secondary(p, uvw, MASS_ELECTRON_EV, PHOTON, .true._C_BOOL) - call particle_create_secondary(p, -uvw, MASS_ELECTRON_EV, PHOTON, .true._C_BOOL) - - p % E = ZERO - p % alive = .false. - - end subroutine sample_positron_reaction - -!=============================================================================== -! SAMPLE_NUCLIDE -!=============================================================================== - - subroutine sample_nuclide(p, base, i_nuclide, i_nuc_mat) - - type(Particle), intent(in) :: p - character(7), intent(in) :: base ! which reaction to sample based on - integer, intent(out) :: i_nuclide - integer, intent(out) :: i_nuc_mat - - real(8) :: prob - real(8) :: cutoff - real(8) :: atom_density ! atom density of nuclide in atom/b-cm - real(8) :: sigma ! microscopic total xs for nuclide - type(Material), pointer :: mat - - ! Get pointer to current material - mat => materials(p % material) - - ! Sample cumulative distribution function - select case (base) - case ('total') - cutoff = prn() * material_xs % total - case ('scatter') - cutoff = prn() * (material_xs % total - material_xs % absorption) - case ('fission') - cutoff = prn() * material_xs % fission - end select - - i_nuc_mat = 0 - prob = ZERO - do while (prob < cutoff) - i_nuc_mat = i_nuc_mat + 1 - - ! Check to make sure that a nuclide was sampled - if (i_nuc_mat > mat % n_nuclides) then - call particle_write_restart(p) - call fatal_error("Did not sample any nuclide during collision.") - end if - - ! Find atom density - i_nuclide = mat % nuclide(i_nuc_mat) - atom_density = mat % atom_density(i_nuc_mat) - - ! Determine microscopic cross section - select case (base) - case ('total') - sigma = atom_density * micro_xs(i_nuclide) % total - case ('scatter') - sigma = atom_density * (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) - case ('fission') - sigma = atom_density * micro_xs(i_nuclide) % fission - end select - - ! Increment probability to compare to cutoff - prob = prob + sigma - end do - - end subroutine sample_nuclide - !=============================================================================== ! SAMPLE_ELEMENT !=============================================================================== @@ -492,402 +249,23 @@ contains end function sample_element -!=============================================================================== -! SAMPLE_FISSION -!=============================================================================== - - subroutine sample_fission(i_nuclide, E, i_reaction) - integer, intent(in) :: i_nuclide ! index in nuclides array - real(8), intent(in) :: E ! incident neutron energy - integer, intent(out) :: i_reaction ! index in nuc % reactions array - - integer :: i - integer :: i_grid - integer :: i_temp - integer :: threshold - real(8) :: f - real(8) :: prob - real(8) :: cutoff - type(Nuclide), pointer :: nuc - - ! Get pointer to nuclide - nuc => nuclides(i_nuclide) - - ! If we're in the URR, by default use the first fission reaction. We also - ! default to the first reaction if we know that there are no partial fission - ! reactions - if (micro_xs(i_nuclide) % use_ptable .or. & - .not. nuc % has_partial_fission) then - i_reaction = nuc % index_fission(1) - return - end if - - ! Check to see if we are in a windowed multipole range. WMP only supports - ! the first fission reaction. - if (nuc % mp_present) then - if (E >= nuc % multipole % E_min .and. & - E <= nuc % multipole % E_max) then - i_reaction = nuc % index_fission(1) - return - end if - end if - - ! Get grid index and interpolatoin factor and sample fission cdf - i_temp = micro_xs(i_nuclide) % index_temp - i_grid = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - cutoff = prn() * micro_xs(i_nuclide) % fission - prob = ZERO - - ! Loop through each partial fission reaction type - - FISSION_REACTION_LOOP: do i = 1, nuc % n_fission - i_reaction = nuc % index_fission(i) - - associate (rx => nuc % reactions(i_reaction)) - ! if energy is below threshold for this reaction, skip it - threshold = rx % xs_threshold(i_temp) - if (i_grid < threshold) cycle - - ! add to cumulative probability - prob = prob + ((ONE - f) * rx % xs(i_temp, i_grid - threshold + 1) & - + f*(rx % xs(i_temp, i_grid - threshold + 2))) - end associate - - ! Create fission bank sites if fission occurs - if (prob > cutoff) exit FISSION_REACTION_LOOP - end do FISSION_REACTION_LOOP - - end subroutine sample_fission - -!=============================================================================== -! SAMPLE_PHOTON_PRODUCT -!=============================================================================== - - subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product) - integer, intent(in) :: i_nuclide ! index in nuclides array - real(8), intent(in) :: E ! energy of neutron - integer, intent(out) :: i_reaction ! index in nuc % reactions array - integer, intent(out) :: i_product ! index in reaction % products array - - integer :: i_grid - integer :: i_temp - integer :: threshold - integer :: last_valid_reaction - integer :: last_valid_product - real(8) :: f - real(8) :: prob - real(8) :: cutoff - real(8) :: yield - - ! Get pointer to nuclide - associate (nuc => nuclides(i_nuclide)) - - ! Get grid index and interpolation factor and sample photon production cdf - i_temp = micro_xs(i_nuclide) % index_temp - i_grid = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - cutoff = prn() * micro_xs(i_nuclide) % photon_prod - prob = ZERO - - ! Loop through each reaction type - REACTION_LOOP: do i_reaction = 1, size(nuc % reactions) - associate (rx => nuc % reactions(i_reaction)) - threshold = rx % xs_threshold(i_temp) - - ! if energy is below threshold for this reaction, skip it - if (i_grid < threshold) cycle - - do i_product = 1, rx % products_size() - if (rx % product_particle(i_product) == PHOTON) then - ! add to cumulative probability - yield = rx % product_yield(i_product, E) - prob = prob + ((ONE - f) * rx % xs(i_temp, i_grid - threshold + 1) & - + f*(rx % xs(i_temp, i_grid - threshold + 2))) * yield - - if (prob > cutoff) return - last_valid_reaction = i_reaction - last_valid_product = i_product - end if - end do - end associate - end do REACTION_LOOP - end associate - - i_reaction = last_valid_reaction - i_product = last_valid_product - - end subroutine sample_photon_product - -!=============================================================================== -! ABSORPTION -!=============================================================================== - - subroutine absorption(p, i_nuclide) - type(Particle), intent(inout) :: p - integer, intent(in) :: i_nuclide - - if (survival_biasing) then - ! Determine weight absorbed in survival biasing - p % absorb_wgt = p % wgt * micro_xs(i_nuclide) % absorption / & - micro_xs(i_nuclide) % total - - ! Adjust weight of particle by probability of absorption - p % wgt = p % wgt - p % absorb_wgt - p % last_wgt = p % wgt - - ! Score implicit absorption estimate of keff - if (run_mode == MODE_EIGENVALUE) then - global_tally_absorption = global_tally_absorption + p % absorb_wgt * & - micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption - end if - else - ! See if disappearance reaction happens - if (micro_xs(i_nuclide) % absorption > & - prn() * micro_xs(i_nuclide) % total) then - ! Score absorption estimate of keff - if (run_mode == MODE_EIGENVALUE) then - global_tally_absorption = global_tally_absorption + p % wgt * & - micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption - end if - - p % alive = .false. - p % event = EVENT_ABSORB - p % event_MT = N_DISAPPEAR - end if - end if - - end subroutine absorption - -!=============================================================================== -! SCATTER -!=============================================================================== - - subroutine scatter(p, i_nuclide, i_nuc_mat) - type(Particle), intent(inout) :: p - integer, intent(in) :: i_nuclide - integer, intent(in) :: i_nuc_mat - - integer :: i - integer :: j - integer :: i_temp - integer :: i_grid - integer :: threshold - real(8) :: f - real(8) :: prob - real(8) :: cutoff - real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering - real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering - real(8) :: phi ! azimuthal angle for iso-in-lab scattering - real(8) :: kT ! temperature in eV - logical :: sampled ! whether or not a reaction type has been sampled - type(Nuclide), pointer :: nuc - - ! copy incoming direction - uvw_old(:) = p % coord(1) % uvw - - ! Get pointer to nuclide and grid index/interpolation factor - nuc => nuclides(i_nuclide) - i_temp = micro_xs(i_nuclide) % index_temp - i_grid = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - - ! For tallying purposes, this routine might be called directly. In that - ! case, we need to sample a reaction via the cutoff variable - cutoff = prn() * (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) - sampled = .false. - - ! Calculate elastic cross section if it wasn't precalculated - if (micro_xs(i_nuclide) % elastic == CACHE_INVALID) then - call nuc % calculate_elastic_xs(micro_xs(i_nuclide)) - end if - - prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal - if (prob > cutoff) then - ! ======================================================================= - ! NON-S(A,B) ELASTIC SCATTERING - - ! Determine temperature - if (nuc % mp_present) then - kT = p % sqrtkT**2 - else - kT = nuc % kTs(micro_xs(i_nuclide) % index_temp) - end if - - ! Perform collision physics for elastic scattering - call elastic_scatter(i_nuclide, nuc % reactions(1), kT, p % E, & - p % coord(1) % uvw, p % mu, p % wgt) - - p % event_MT = ELASTIC - sampled = .true. - end if - - prob = micro_xs(i_nuclide) % elastic - if (prob > cutoff .and. .not. sampled) then - ! ======================================================================= - ! S(A,B) SCATTERING - - call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, p % E, & - p % coord(1) % uvw, p % mu) - - p % event_MT = ELASTIC - sampled = .true. - end if - - if (.not. sampled) then - ! ======================================================================= - ! INELASTIC SCATTERING - - j = 0 - do while (prob < cutoff) - j = j + 1 - i = nuc % index_inelastic_scatter(j) - - ! Check to make sure inelastic scattering reaction sampled - if (i > size(nuc % reactions)) then - call particle_write_restart(p) - call fatal_error("Did not sample any reaction for nuclide " & - &// trim(nuc % name)) - end if - - associate (rx => nuc % reactions(i)) - ! if energy is below threshold for this reaction, skip it - threshold = rx % xs_threshold(i_temp) - if (i_grid < threshold) cycle - - ! add to cumulative probability - prob = prob + ((ONE - f)*rx % xs(i_temp, i_grid - threshold + 1) & - + f*(rx % xs(i_temp, i_grid - threshold + 2))) - end associate - end do - - ! Perform collision physics for inelastic scattering - call inelastic_scatter(nuc, nuc%reactions(i), p) - p % event_MT = nuc % reactions(i) % MT - - end if - - ! Set event component - p % event = EVENT_SCATTER - - ! Sample new outgoing angle for isotropic-in-lab scattering - associate (mat => materials(p % material)) - if (mat % has_isotropic_nuclides) then - if (materials(p % material) % p0(i_nuc_mat)) then - ! Sample isotropic-in-lab outgoing direction - uvw_new(1) = TWO * prn() - ONE - phi = TWO * PI * prn() - uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) - uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) - p % mu = dot_product(uvw_old, uvw_new) - - ! Change direction of particle - p % coord(1) % uvw = uvw_new - end if - end if - end associate - - end subroutine scatter - -!=============================================================================== -! ELASTIC_SCATTER treats the elastic scattering of a neutron with a -! target. -!=============================================================================== - - subroutine elastic_scatter(i_nuclide, rxn, kT, E, uvw, mu_lab, wgt) - integer, intent(in) :: i_nuclide - type(Reaction), intent(in) :: rxn - real(8), intent(in) :: kT ! temperature in eV - real(8), intent(inout) :: E - real(8), intent(inout) :: uvw(3) - real(8), intent(out) :: mu_lab - real(8), intent(inout) :: wgt - - real(8) :: awr ! atomic weight ratio of target - real(8) :: mu_cm ! cosine of polar angle in center-of-mass - real(8) :: vel ! magnitude of velocity - real(8) :: v_n(3) ! velocity of neutron - real(8) :: v_cm(3) ! velocity of center-of-mass - real(8) :: v_t(3) ! velocity of target nucleus - real(8) :: uvw_cm(3) ! directional cosines in center-of-mass - type(Nuclide), pointer :: nuc - - ! get pointer to nuclide - nuc => nuclides(i_nuclide) - - vel = sqrt(E) - awr = nuc % awr - - ! Neutron velocity in LAB - v_n = vel * uvw - - ! Sample velocity of target nucleus - if (.not. micro_xs(i_nuclide) % use_ptable) then - call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, & - micro_xs(i_nuclide) % elastic, kT) - else - v_t = ZERO - end if - - ! Velocity of center-of-mass - v_cm = (v_n + awr*v_t)/(awr + ONE) - - ! Transform to CM frame - v_n = v_n - v_cm - - ! Find speed of neutron in CM - vel = sqrt(dot_product(v_n, v_n)) - - ! Sample scattering angle - mu_cm = rxn % sample_elastic_mu(E) - - ! Determine direction cosines in CM - uvw_cm = v_n/vel - - ! Rotate neutron velocity vector to new angle -- note that the speed of the - ! neutron in CM does not change in elastic scattering. However, the speed - ! will change when we convert back to LAB - v_n = vel * rotate_angle(uvw_cm, mu_cm) - - ! Transform back to LAB frame - v_n = v_n + v_cm - - E = dot_product(v_n, v_n) - vel = sqrt(E) - - ! compute cosine of scattering angle in LAB frame by taking dot product of - ! neutron's pre- and post-collision angle - mu_lab = dot_product(uvw, v_n) / vel - - ! Set energy and direction of particle in LAB frame - uvw = v_n / vel - - ! Because of floating-point roundoff, it may be possible for mu_lab to be - ! outside of the range [-1,1). In these cases, we just set mu_lab to exactly - ! -1 or 1 - - if (abs(mu_lab) > ONE) mu_lab = sign(ONE,mu_lab) - - end subroutine elastic_scatter - !=============================================================================== ! SAB_SCATTER performs thermal scattering of a particle with a bound scatterer ! according to a specified S(a,b) table. !=============================================================================== - subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu) - integer, intent(in) :: i_nuclide ! index in micro_xs - integer, intent(in) :: i_sab ! index in sab_tables - real(8), intent(inout) :: E ! incoming/outgoing energy - real(8), intent(inout) :: uvw(3) ! directional cosines - real(8), intent(out) :: mu ! scattering cosine + subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu) bind(C) + integer(C_INT), value :: i_nuclide ! index in micro_xs + integer(C_INT), value :: i_sab ! index in sab_tables + real(C_DOUBLE), intent(inout) :: E ! incoming/outgoing energy + real(C_DOUBLE), intent(inout) :: uvw(3) ! directional cosines + real(C_DOUBLE), intent(out) :: mu ! scattering cosine real(C_DOUBLE) :: E_out type(C_PTR) :: ptr ! Sample from C++ side - ptr = C_LOC(micro_xs(i_nuclide)) + ptr = C_LOC(micro_xs(i_nuclide + 1)) call sab_tables(i_sab) % sample(ptr, E, E_out, mu) ! Set energy to outgoing, change direction of particle @@ -895,602 +273,4 @@ contains uvw = rotate_angle(uvw, mu) end subroutine sab_scatter -!=============================================================================== -! SAMPLE_TARGET_VELOCITY samples the target velocity. The constant cross section -! free gas model is the default method. Methods for correctly accounting -! for the energy dependence of cross sections in treating resonance elastic -! scattering such as the DBRC, WCM, and a new, accelerated scheme are also -! implemented here. -!=============================================================================== - - subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff, kT) - type(Nuclide), intent(in) :: nuc ! target nuclide at temperature T - real(8), intent(out) :: v_target(3) ! target velocity - real(8), intent(in) :: E ! particle energy - real(8), intent(in) :: uvw(3) ! direction cosines - real(8), intent(in) :: v_neut(3) ! neutron velocity - real(8), intent(inout) :: wgt ! particle weight - real(8), intent(in) :: xs_eff ! effective elastic xs at temperature T - real(8), intent(in) :: kT ! equilibrium temperature of target in eV - - real(8) :: awr ! target/neutron mass ratio - real(8) :: E_rel ! trial relative energy - real(8) :: xs_0K ! 0K xs at E_rel - real(8) :: wcf ! weight correction factor - real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1) - real(8) :: E_low ! lowest practical relative energy - real(8) :: E_up ! highest practical relative energy - real(8) :: E_t ! trial target energy - real(8) :: xs_max ! max 0K xs over practical relative energies - real(8) :: xs_low ! 0K xs at lowest practical relative energy - real(8) :: xs_up ! 0K xs at highest practical relative energy - real(8) :: m ! slope for interpolation - real(8) :: R ! rejection criterion for DBRC / target speed - real(8) :: cdf_low ! xs cdf at lowest practical relative energy - real(8) :: cdf_up ! xs cdf at highest practical relative energy - real(8) :: cdf_rel ! trial xs cdf value - real(8) :: mu ! cosine between neutron and target velocities - - integer :: i_E_low ! 0K index to lowest practical relative energy - integer :: i_E_up ! 0K index to highest practical relative energy - integer :: i_E_rel ! index to trial relative energy - integer :: n_grid ! number of energies on 0K grid - - integer :: sampling_method ! method of target velocity sampling - - awr = nuc % awr - - ! check if nuclide is a resonant scatterer - if (nuc % resonant) then - - ! sampling method to use - sampling_method = res_scat_method - - ! upper resonance scattering energy bound (target is at rest above this E) - if (E > res_scat_energy_max) then - v_target = ZERO - return - - ! lower resonance scattering energy bound (should be no resonances below) - else if (E < res_scat_energy_min) then - sampling_method = RES_SCAT_CXS - end if - - ! otherwise, use free gas model - else - if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then - v_target = ZERO - return - else - sampling_method = RES_SCAT_CXS - end if - end if - - ! use appropriate target velocity sampling method - select case (sampling_method) - case (RES_SCAT_CXS) - - ! sample target velocity with the constant cross section (cxs) approx. - call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - - case (RES_SCAT_WCM) - - ! sample target velocity with the constant cross section (cxs) approx. - call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - - ! adjust weight as prescribed by the weight correction method (wcm) - E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - xs_0K = elastic_xs_0K(E_rel, nuc) - wcf = xs_0K / xs_eff - wgt = wcf * wgt - - case (RES_SCAT_DBRC, RES_SCAT_ARES) - E_red = sqrt(awr * E / kT) - E_low = max(ZERO, E_red - FOUR)**2 * kT / awr - E_up = (E_red + FOUR)**2 * kT / awr - - ! find lower and upper energy bound indices - ! lower index - n_grid = size(nuc % energy_0K) - if (E_low < nuc % energy_0K(1)) then - i_E_low = 1 - elseif (E_low > nuc % energy_0K(n_grid)) then - i_E_low = n_grid - 1 - else - i_E_low = binary_search(nuc % energy_0K, n_grid, E_low) - end if - - ! upper index - if (E_up < nuc % energy_0K(1)) then - i_E_up = 1 - elseif (E_up > nuc % energy_0K(n_grid)) then - i_E_up = n_grid - 1 - else - i_E_up = binary_search(nuc % energy_0K, n_grid, E_up) - end if - - if (i_E_up == i_E_low) then - ! Handle degenerate case -- if the upper/lower bounds occur for the same - ! index, then using cxs is probably a good approximation - call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - - else - if (sampling_method == RES_SCAT_DBRC) then - ! interpolate xs since we're not exactly at the energy indices - xs_low = nuc % elastic_0K(i_E_low) - m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & - / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) - xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low)) - xs_up = nuc % elastic_0K(i_E_up) - m = (nuc % elastic_0K(i_E_up + 1) - xs_up) & - / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) - xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up)) - - ! get max 0K xs value over range of practical relative energies - xs_max = max(xs_low, & - maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up)), xs_up) - - DBRC_REJECT_LOOP: do - TARGET_ENERGY_LOOP: do - ! sample target velocity with the constant cross section (cxs) approx. - call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - if (E_rel < E_up) exit TARGET_ENERGY_LOOP - end do TARGET_ENERGY_LOOP - - ! perform Doppler broadening rejection correction (dbrc) - xs_0K = elastic_xs_0K(E_rel, nuc) - R = xs_0K / xs_max - if (prn() < R) exit DBRC_REJECT_LOOP - end do DBRC_REJECT_LOOP - - elseif (sampling_method == RES_SCAT_ARES) then - ! interpolate xs CDF since we're not exactly at the energy indices - ! cdf value at lower bound attainable energy - m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) & - / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) - cdf_low = nuc % xs_cdf(i_E_low - 1) & - + m * (E_low - nuc % energy_0K(i_E_low)) - if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO - - ! cdf value at upper bound attainable energy - m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) & - / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) - cdf_up = nuc % xs_cdf(i_E_up - 1) & - + m * (E_up - nuc % energy_0K(i_E_up)) - - ARES_REJECT_LOOP: do - - ! directly sample Maxwellian - E_t = -kT * log(prn()) - - ! sample a relative energy using the xs cdf - cdf_rel = cdf_low + prn() * (cdf_up - cdf_low) - i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), & - i_E_up - i_E_low + 2, cdf_rel) - E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1) - m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) & - - nuc % xs_cdf(i_E_low + i_E_rel - 2)) & - / (nuc % energy_0K(i_E_low + i_E_rel) & - - nuc % energy_0K(i_E_low + i_E_rel - 1)) - E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m - - ! perform rejection sampling on cosine between - ! neutron and target velocities - mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t)) - - if (abs(mu) < ONE) then - ! set and accept target velocity - E_t = E_t / awr - v_target = sqrt(E_t) * rotate_angle(uvw, mu) - exit ARES_REJECT_LOOP - end if - end do ARES_REJECT_LOOP - end if - end if - end select - - end subroutine sample_target_velocity - -!=============================================================================== -! SAMPLE_CXS_TARGET_VELOCITY samples a target velocity based on the free gas -! scattering formulation, used by most Monte Carlo codes, in which cross section -! is assumed to be constant in energy. Excellent documentation for this method -! can be found in FRA-TM-123. -!=============================================================================== - - subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - type(Nuclide), intent(in) :: nuc ! target nuclide at temperature - real(8), intent(out) :: v_target(3) - real(8), intent(in) :: E - real(8), intent(in) :: uvw(3) - real(8), intent(in) :: kT ! equilibrium temperature of target in eV - - real(8) :: awr ! target/neutron mass ratio - real(8) :: alpha ! probability of sampling f2 over f1 - real(8) :: mu ! cosine of angle between neutron and target vel - real(8) :: r1, r2 ! pseudo-random numbers - real(8) :: c ! cosine used in maxwell sampling - real(8) :: accept_prob ! probability of accepting combination of vt and mu - real(8) :: beta_vn ! beta * speed of neutron - real(8) :: beta_vt ! beta * speed of target - real(8) :: beta_vt_sq ! (beta * speed of target)^2 - real(8) :: vt ! speed of target - - awr = nuc % awr - - beta_vn = sqrt(awr * E / kT) - alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) - - do - ! Sample two random numbers - r1 = prn() - r2 = prn() - - if (prn() < alpha) then - ! With probability alpha, we sample the distribution p(y) = - ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte - ! Carlo sampler - - beta_vt_sq = -log(r1*r2) - - else - ! With probability 1-alpha, we sample the distribution p(y) = y^2 * - ! e^(-y^2). This can be done with sampling scheme C61 from the Monte - ! Carlo sampler - - c = cos(PI/TWO * prn()) - beta_vt_sq = -log(r1) - log(r2)*c*c - end if - - ! Determine beta * vt - beta_vt = sqrt(beta_vt_sq) - - ! Sample cosine of angle between neutron and target velocity - mu = TWO*prn() - ONE - - ! Determine rejection probability - accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & - /(beta_vn + beta_vt) - - ! Perform rejection sampling on vt and mu - if (prn() < accept_prob) exit - end do - - ! Determine speed of target nucleus - vt = sqrt(beta_vt_sq*kT/awr) - - ! Determine velocity vector of target nucleus based on neutron's velocity - ! and the sampled angle between them - v_target = vt * rotate_angle(uvw, mu) - - end subroutine sample_cxs_target_velocity - -!=============================================================================== -! CREATE_FISSION_SITES determines the average total, prompt, and delayed -! neutrons produced from fission and creates appropriate bank sites. -!=============================================================================== - - subroutine create_fission_sites(p, i_nuclide, i_reaction, bank_array, size_bank) - type(Particle), intent(inout) :: p - integer, intent(in) :: i_nuclide - integer, intent(in) :: i_reaction - type(Bank), intent(inout) :: bank_array(:) - integer(8), intent(inout) :: size_bank - - integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born - integer :: i ! loop index - integer :: nu ! actual number of neutrons produced - real(8) :: nu_t ! total nu - real(8) :: weight ! weight adjustment for ufs method - type(Nuclide), pointer :: nuc - - interface - function ufs_get_weight(p) result(weight) bind(C) - import Particle, C_DOUBLE - type(Particle), intent(in) :: p - real(C_DOUBLE) :: WEIGHT - end function - end interface - - ! Get pointers - nuc => nuclides(i_nuclide) - - ! TODO: Heat generation from fission - - ! If uniform fission source weighting is turned on, we increase of decrease - ! the expected number of fission sites produced - - if (ufs) then - weight = ufs_get_weight(p) - else - weight = ONE - end if - - ! Determine expected number of neutrons produced - nu_t = p % wgt / keff * weight * micro_xs(i_nuclide) % nu_fission / & - micro_xs(i_nuclide) % total - - ! Sample number of neutrons produced - if (prn() > nu_t - int(nu_t)) then - nu = int(nu_t) - else - nu = int(nu_t) + 1 - end if - - ! Check for bank size getting hit. For fixed source calculations, this is a - ! fatal error. For eigenvalue calculations, it just means that k-effective - ! was too high for a single batch. - if (size_bank + nu > size(bank_array)) then - if (run_mode == MODE_FIXEDSOURCE) then - call fatal_error("Secondary particle bank size limit reached. If you & - &are running a subcritical multiplication problem, k-effective & - &may be too close to one.") - else - if (master) call warning("Maximum number of sites in fission bank & - &reached. This can result in irreproducible results using different & - &numbers of processes/threads.") - end if - end if - - ! Bank source neutrons - if (nu == 0 .or. size_bank == size(bank_array)) return - - ! Initialize counter of delayed neutrons encountered for each delayed group - ! to zero. - nu_d(:) = 0 - - p % fission = .true. ! Fission neutrons will be banked - do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4) - ! Bank source neutrons by copying particle data - bank_array(i) % xyz = p % coord(1) % xyz - - ! Set particle as neutron - bank_array(i) % particle = NEUTRON - - ! Set weight of fission bank site - bank_array(i) % wgt = ONE/weight - - ! Sample delayed group and angle/energy for fission reaction - call sample_fission_neutron(nuc, nuc % reactions(i_reaction), & - p % E, bank_array(i)) - - ! Set delayed group on particle too - p % delayed_group = bank_array(i) % delayed_group - - ! Increment the number of neutrons born delayed - if (p % delayed_group > 0) then - nu_d(p % delayed_group) = nu_d(p % delayed_group) + 1 - end if - end do - - ! increment number of bank sites - size_bank = min(size_bank + nu, int(size(bank_array),8)) - - ! Store total and delayed weight banked for analog fission tallies - p % n_bank = nu - p % wgt_bank = nu/weight - p % n_delayed_bank(:) = nu_d(:) - - end subroutine create_fission_sites - -!=============================================================================== -! SAMPLE_FISSION_NEUTRON -!=============================================================================== - - subroutine sample_fission_neutron(nuc, rxn, E_in, site) - type(Nuclide), intent(in) :: nuc - type(Reaction), intent(in) :: rxn - real(8), intent(in) :: E_in - type(Bank), intent(inout) :: site - - integer :: group ! index on nu energy grid / precursor group - integer :: n_sample ! number of resamples - real(8) :: nu_t ! total nu - real(8) :: nu_d ! delayed nu - real(8) :: beta ! delayed neutron fraction - real(8) :: xi ! random number - real(8) :: yield ! delayed neutron precursor yield - real(8) :: prob ! cumulative probability - real(8) :: mu ! cosine of scattering angle - real(8) :: phi ! azimuthal angle - - ! Sample cosine of angle -- fission neutrons are always emitted - ! isotropically. Sometimes in ACE data, fission reactions actually have - ! an angular distribution listed, but for those that do, it's simply just - ! a uniform distribution in mu - mu = TWO * prn() - ONE - - ! Sample azimuthal angle uniformly in [0,2*pi) - phi = TWO*PI*prn() - site % uvw(1) = mu - site % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) - site % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) - - ! Determine total nu, delayed nu, and delayed neutron fraction - nu_t = nuc % nu(E_in, EMISSION_TOTAL) - nu_d = nuc % nu(E_in, EMISSION_DELAYED) - beta = nu_d / nu_t - - if (prn() < beta) then - ! ==================================================================== - ! DELAYED NEUTRON SAMPLED - - ! sampled delayed precursor group - xi = prn()*nu_d - prob = ZERO - do group = 1, nuc % n_precursor - - ! determine delayed neutron precursor yield for group j - yield = rxn % product_yield(1 + group, E_in) - - ! Check if this group is sampled - prob = prob + yield - if (xi < prob) exit - end do - - ! if the sum of the probabilities is slightly less than one and the - ! random number is greater, j will be greater than nuc % - ! n_precursor -- check for this condition - group = min(group, nuc % n_precursor) - - ! set the delayed group for the particle born from fission - site % delayed_group = group - - n_sample = 0 - do - ! sample from energy/angle distribution -- note that mu has already been - ! sampled above and doesn't need to be resampled - call rxn % product_sample(1 + group, E_in, site % E, mu) - - ! resample if energy is greater than maximum neutron energy - if (site % E < energy_max(NEUTRON)) exit - - ! check for large number of resamples - n_sample = n_sample + 1 - if (n_sample == MAX_SAMPLE) then - ! call particle_write_restart(p) - call fatal_error("Resampled energy distribution maximum number of " & - // "times for nuclide " // nuc % name) - end if - end do - - else - ! ==================================================================== - ! PROMPT NEUTRON SAMPLED - - ! set the delayed group for the particle born from fission to 0 - site % delayed_group = 0 - - ! sample from prompt neutron energy distribution - n_sample = 0 - do - call rxn % product_sample(1, E_in, site % E, mu) - - ! resample if energy is greater than maximum neutron energy - if (site % E < energy_max(NEUTRON)) exit - - ! check for large number of resamples - n_sample = n_sample + 1 - if (n_sample == MAX_SAMPLE) then - ! call particle_write_restart(p) - call fatal_error("Resampled energy distribution maximum number of " & - // "times for nuclide " // nuc % name) - end if - end do - end if - - end subroutine sample_fission_neutron - -!=============================================================================== -! INELASTIC_SCATTER handles all reactions with a single secondary neutron (other -! than fission), i.e. level scattering, (n,np), (n,na), etc. -!=============================================================================== - - subroutine inelastic_scatter(nuc, rxn, p) - type(Nuclide), intent(in) :: nuc - type(Reaction), intent(in) :: rxn - type(Particle), intent(inout) :: p - - integer :: i ! loop index - real(8) :: E ! energy in lab (incoming/outgoing) - real(8) :: mu ! cosine of scattering angle in lab - real(8) :: A ! atomic weight ratio of nuclide - real(8) :: E_in ! incoming energy - real(8) :: E_cm ! outgoing energy in center-of-mass - real(8) :: yield ! neutron yield - - ! copy energy of neutron - E_in = p % E - - ! sample outgoing energy and scattering cosine - call rxn % product_sample(1, E_in, E, mu) - - ! if scattering system is in center-of-mass, transfer cosine of scattering - ! angle and outgoing energy from CM to LAB - if (rxn % scatter_in_cm) then - E_cm = E - - ! determine outgoing energy in lab - A = nuc%awr - E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & - / ((A+ONE)*(A+ONE)) - - ! determine outgoing angle in lab - mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) - end if - - ! Because of floating-point roundoff, it may be possible for mu to be - ! outside of the range [-1,1). In these cases, we just set mu to exactly -1 - ! or 1 - if (abs(mu) > ONE) mu = sign(ONE,mu) - - ! Set outgoing energy and scattering angle - p % E = E - p % mu = mu - - ! change direction of particle - p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) - - ! evaluate yield - yield = rxn % product_yield(1, E_in) - if (mod(yield, ONE) == ZERO) then - ! If yield is integral, create exactly that many secondary particles - do i = 1, nint(yield) - 1 - call particle_create_secondary(p, p % coord(1) % uvw, p % E, & - NEUTRON, run_CE=.true._C_BOOL) - end do - else - ! Otherwise, change weight of particle based on yield - p % wgt = yield * p % wgt - end if - - end subroutine inelastic_scatter - -!=============================================================================== -! SAMPLE_SECONDARY_PHOTONS -!=============================================================================== - - subroutine sample_secondary_photons(p, i_nuclide) - type(Particle), intent(inout) :: p - integer, intent(in) :: i_nuclide - - integer :: i_reaction ! index in nuc % reactions array - integer :: i_product ! index in nuc % reactions % products array - - real(8) :: nu_t - real(8) :: mu - real(8) :: E - real(8) :: uvw(3) - integer :: nu - integer :: i - - ! Sample the number of photons produced - nu_t = p % wgt * micro_xs(i_nuclide) % photon_prod / & - micro_xs(i_nuclide) % total - if (prn() > nu_t - int(nu_t)) then - nu = int(nu_t) - else - nu = int(nu_t) + 1 - end if - - ! Sample each secondary photon - do i = 1, nu - - ! Sample the reaction and product - call sample_photon_product(i_nuclide, p % E, i_reaction, i_product) - - ! Sample the outgoing energy and angle - call nuclides(i_nuclide) % reactions(i_reaction) % & - product_sample(i_product, p % E, E, mu) - - ! Sample the new direction - uvw = rotate_angle(p % coord(1) % uvw, mu) - - ! Create the secondary photon - call particle_create_secondary(p, uvw, E, PHOTON, run_CE=.true._C_BOOL) - end do - - end subroutine sample_secondary_photons - end module physics diff --git a/src/physics.cpp b/src/physics.cpp new file mode 100644 index 0000000000..43fc2871b4 --- /dev/null +++ b/src/physics.cpp @@ -0,0 +1,1181 @@ +#include "openmc/physics.h" + +#include "openmc/bank.h" +#include "openmc/constants.h" +#include "openmc/eigenvalue.h" +#include "openmc/error.h" +#include "openmc/material.h" +#include "openmc/math_functions.h" +#include "openmc/message_passing.h" +#include "openmc/nuclide.h" +#include "openmc/photon.h" +#include "openmc/physics_common.h" +#include "openmc/random_lcg.h" +#include "openmc/reaction.h" +#include "openmc/secondary_uncorrelated.h" +#include "openmc/search.h" +#include "openmc/settings.h" +#include "openmc/simulation.h" +#include "openmc/thermal.h" +#include "openmc/tallies/tally.h" + +#include // for max, min, max_element +#include // for sqrt, exp, log, abs, copysign +#include + +namespace openmc { + +//============================================================================== +// Non-member functions +//============================================================================== + +void collision(Particle* p) +{ + // Add to collision counter for particle + ++(p->n_collision); + + // Sample reaction for the material the particle is in + switch (static_cast(p->type)) { + case ParticleType::neutron: + sample_neutron_reaction(p); + break; + case ParticleType::photon: + sample_photon_reaction(p); + break; + case ParticleType::electron: + sample_electron_reaction(p); + break; + case ParticleType::positron: + sample_positron_reaction(p); + break; + } + + // Kill particle if energy falls below cutoff + if (p->E < settings::energy_cutoff[p->type - 1]) { + p->alive = false; + p->wgt = 0.0; + p->last_wgt = 0.0; + } + + // Display information about collision + if (settings::verbosity >= 10 || simulation::trace) { + std::stringstream msg; + if (static_cast(p->type) == ParticleType::neutron) { + msg << " " << reaction_name(p->event_MT) << " with " << + data::nuclides[p->event_nuclide-1]->name_ << ". Energy = " << p->E << " eV."; + } else { + msg << " " << reaction_name(p->event_MT) << " with " << + data::elements[p->event_nuclide-1].name_ << ". Energy = " << p->E << " eV."; + } + write_message(msg, 1); + } +} + +void sample_neutron_reaction(Particle* p) +{ + int i_nuclide; + int i_nuc_mat; + sample_nuclide(p, SCORE_TOTAL, &i_nuclide, &i_nuc_mat); + + // Save which nuclide particle had collision with + // TODO: off-by-one + p->event_nuclide = i_nuclide + 1; + + // Create fission bank sites. Note that while a fission reaction is sampled, + // it never actually "happens", i.e. the weight of the particle does not + // change when sampling fission sites. The following block handles all + // absorption (including fission) + + const auto& nuc {data::nuclides[i_nuclide]}; + + if (nuc->fissionable_) { + Reaction* rx = sample_fission(i_nuclide, p->E); + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + create_fission_sites(p, i_nuclide, rx, simulation::fission_bank.data(), + &simulation::n_bank, simulation::fission_bank.size()); + } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE && + settings::create_fission_neutrons) { + create_fission_sites(p, i_nuclide, rx, p->secondary_bank, + &p->n_secondary, MAX_SECONDARY); + } + } + + // Create secondary photons + if (settings::photon_transport) { + prn_set_stream(STREAM_PHOTON); + sample_secondary_photons(p, i_nuclide); + prn_set_stream(STREAM_TRACKING); + } + + // If survival biasing is being used, the following subroutine adjusts the + // weight of the particle. Otherwise, it checks to see if absorption occurs + + if (simulation::micro_xs[i_nuclide].absorption > 0.0) { + absorption(p, i_nuclide); + } else { + p->absorb_wgt = 0.0; + } + if (!p->alive) return; + + // Sample a scattering reaction and determine the secondary energy of the + // exiting neutron + scatter(p, i_nuclide, i_nuc_mat); + + // Advance URR seed stream 'N' times after energy changes + if (p->E != p->last_E) { + prn_set_stream(STREAM_URR_PTABLE); + advance_prn_seed(data::nuclides.size()); + prn_set_stream(STREAM_TRACKING); + } + + // Play russian roulette if survival biasing is turned on + if (settings::survival_biasing) { + russian_roulette(p); + if (!p->alive) return; + } +} + +void +create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_array, + int64_t* size_bank, int64_t bank_capacity) +{ + // TODO: Heat generation from fission + + // If uniform fission source weighting is turned on, we increase or decrease + // the expected number of fission sites produced + double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; + + // Determine the expected number of neutrons produced + double nu_t = p->wgt / simulation::keff * weight * simulation::micro_xs[ + i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].total; + + // Sample the number of neutrons produced + int nu = static_cast(nu_t); + if (prn() <= (nu_t - nu)) ++nu; + + // Check for the bank size getting hit. For fixed source calculations, this + // is a fatal error; for eigenvalue calculations, it just means that k-eff + // was too high for a single batch. + if (*size_bank + nu > bank_capacity) { + if (settings::run_mode == RUN_MODE_FIXEDSOURCE) { + throw std::runtime_error{"Secondary particle bank size limit reached." + " If you are running a subcritical multiplication problem," + " k-effective may be too close to one."}; + } else { + if (mpi::master) { + warning("Maximum number of sites in fission bank reached. This can" + " result in irreproducible results using different numbers of" + " processes/threads."); + } + } + } + + // Begin banking the source neutrons + // First, if our bank is full then don't continue + if (nu == 0 || *size_bank == bank_capacity) return; + + // Initialize the counter of delayed neutrons encountered for each delayed + // group. + double nu_d[MAX_DELAYED_GROUPS] = {0.}; + + p->fission = true; + for (size_t i = *size_bank; i < std::min(*size_bank + nu, bank_capacity); ++i) { + // Bank source neutrons by copying the particle data + bank_array[i].xyz[0] = p->coord[0].xyz[0]; + bank_array[i].xyz[1] = p->coord[0].xyz[1]; + bank_array[i].xyz[2] = p->coord[0].xyz[2]; + + // Set that the bank particle is a neutron + bank_array[i].particle = static_cast(ParticleType::neutron); + + // Set the weight of the fission bank site + bank_array[i].wgt = 1. / weight; + + // Sample delayed group and angle/energy for fission reaction + sample_fission_neutron(i_nuclide, rx, p->E, &bank_array[i]); + + // Set the delayed group on the particle as well + p->delayed_group = bank_array[i].delayed_group; + + // Increment the number of neutrons born delayed + if (p->delayed_group > 0) { + nu_d[p->delayed_group-1]++; + } + } + + // Increment number of bank sites + *size_bank = std::min(*size_bank + nu, bank_capacity); + + // Store the total weight banked for analog fission tallies + p->n_bank = nu; + p->wgt_bank = nu / weight; + for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { + p->n_delayed_bank[d] = nu_d[d]; + } +} + +// TODO: Finish converting photon physics functions + +// void sample_photon_reaction(Particle* p) +// { +// // Kill photon if below energy cutoff -- an extra check is made here because +// // photons with energy below the cutoff may have been produced by neutrons +// // reactions or atomic relaxation +// if (p->E < energy_cutoff(PHOTON)) { +// p->E = 0.0 +// p->alive = false +// return +// } + +// // Sample element within material +// i_element = sample_element(p) +// p->event_nuclide = i_element + +// // Calculate photon energy over electron rest mass equivalent +// alpha = p->E/MASS_ELECTRON_EV + +// // For tallying purposes, this routine might be called directly. In that +// // case, we need to sample a reaction via the cutoff variable +// prob = 0.0 +// cutoff = prn() * micro_photon_xs(i_element) % total + +// associate (elm => elements(i_element)) +// // Coherent (Rayleigh) scattering +// prob = prob + micro_photon_xs(i_element) % coherent +// if (prob > cutoff) { +// rayleigh_scatter(elm, alpha, mu) +// p->coord(1) % uvw = rotate_angle(p->coord(1) % uvw, mu) +// p->event_MT = COHERENT +// return +// } + +// // Incoherent (Compton) scattering +// prob = prob + micro_photon_xs(i_element) % incoherent +// if (prob > cutoff) { +// compton_scatter(elm, alpha, alpha_out, mu, i_shell, true) + +// // Determine binding energy of shell. The binding energy is 0.0 if +// // doppler broadening is not used. +// if (i_shell == 0) { +// e_b = 0.0 +// } else { +// e_b = elm % binding_energy(i_shell) +// } + +// // Create Compton electron +// E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b +// mu_electron = (alpha - alpha_out*mu) & +// / std::sqrt(alpha**2 + alpha_out**2 - 2.0*alpha*alpha_out*mu) +// phi = 2.0*PI*prn() +// uvw = rotate_angle(p->coord(1) % uvw, mu_electron, phi) +// particle_create_secondary(p, uvw, E_electron, ELECTRON, true) + +// // TODO: Compton subshell data does not match atomic relaxation data +// // Allow electrons to fill orbital and produce auger electrons +// // and fluorescent photons +// if (i_shell > 0) { +// atomic_relaxation(p, elm, i_shell) +// } + +// phi = phi + PI +// p->E = alpha_out*MASS_ELECTRON_EV +// p->coord(1) % uvw = rotate_angle(p->coord(1) % uvw, mu, phi) +// p->event_MT = INCOHERENT +// return +// } + +// // Photoelectric effect +// prob_after = prob + micro_photon_xs(i_element) % photoelectric +// if (prob_after > cutoff) { +// do i_shell = 1, size(elm % shells) +// // Get grid index and interpolation factor +// i_grid = micro_photon_xs(i_element) % index_grid +// f = micro_photon_xs(i_element) % interp_factor + +// // Check threshold of reaction +// i_start = elm % shells(i_shell) % threshold +// if (i_grid <= i_start) cycle + +// // Evaluation subshell photoionization cross section +// xs = std::exp(elm % shells(i_shell) % cross_section(i_grid - i_start) + & +// f*(elm % shells(i_shell) % cross_section(i_grid + 1 - i_start) - & +// elm % shells(i_shell) % cross_section(i_grid - i_start))) + +// prob = prob + xs +// if (prob > cutoff) { +// E_electron = p->E - elm % shells(i_shell) % binding_energy + +// // Sample mu using non-relativistic Sauter distribution. +// // See Eqns 3.19 and 3.20 in "Implementing a photon physics +// // model in Serpent 2" by Toni Kaltiaisenaho +// SAMPLE_MU: do +// r = prn() +// if (FOUR * (1.0 - r) * r >= prn()) { +// rel_vel = std::sqrt(E_electron * (E_electron + 2.0 * MASS_ELECTRON_EV))& +// / (E_electron + MASS_ELECTRON_EV) +// mu = (2.0 * r + rel_vel - 1.0) / & +// (2.0 * rel_vel * r - rel_vel + 1.0) +// exit SAMPLE_MU +// } +// end do SAMPLE_MU + +// phi = 2.0*PI*prn() +// uvw(1) = mu +// uvw(2) = std::sqrt(1.0 - mu*mu)*std::cos(phi) +// uvw(3) = std::sqrt(1.0 - mu*mu)*std::sin(phi) + +// // Create secondary electron +// particle_create_secondary(p, uvw, E_electron, ELECTRON, & +// run_CE=true) + +// // Allow electrons to fill orbital and produce auger electrons +// // and fluorescent photons +// atomic_relaxation(p, elm, i_shell) +// p->event_MT = 533 + elm % shells(i_shell) % index_subshell +// p->alive = false +// p->E = 0.0 + +// return +// } +// end do +// } +// prob = prob_after + +// // Pair production +// prob = prob + micro_photon_xs(i_element) % pair_production +// if (prob > cutoff) { +// pair_production(elm, alpha, E_electron, E_positron, mu_electron, & +// mu_positron) + +// // Create secondary electron +// uvw = rotate_angle(p->coord(1) % uvw, mu_electron) +// particle_create_secondary(p, uvw, E_electron, ELECTRON, true) + +// // Create secondary positron +// uvw = rotate_angle(p->coord(1) % uvw, mu_positron) +// particle_create_secondary(p, uvw, E_positron, POSITRON, true) + +// p->event_MT = PAIR_PROD +// p->alive = false +// p->E = 0.0 +// } + +// end associate +// } + +void sample_electron_reaction(Particle* p) +{ + // TODO: create reaction types + + if (settings::electron_treatment == ELECTRON_TTB) { + double E_lost; + thick_target_bremsstrahlung(p, &E_lost); + } + + p->E = 0.0; + p->alive = false; +} + +void sample_positron_reaction(Particle* p) +{ + // TODO: create reaction types + + if (settings::electron_treatment == ELECTRON_TTB) { + double E_lost; + thick_target_bremsstrahlung(p, &E_lost); + } + + // Sample angle isotropically + double mu = 2.0*prn() - 1.0; + double phi = 2.0*PI*prn(); + std::array uvw; + uvw[0] = mu; + uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi); + uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi); + + // Create annihilation photon pair traveling in opposite directions + int photon = static_cast(ParticleType::photon); + p->create_secondary(uvw.data(), MASS_ELECTRON_EV, photon, true); + + uvw[0] = -uvw[0]; + uvw[1] = -uvw[1]; + uvw[2] = -uvw[2]; + p->create_secondary(uvw.data(), MASS_ELECTRON_EV, photon, true); + + p->E = 0.0; + p->alive = false; +} + +void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat) +{ + // Sample cumulative distribution function + double cutoff; + switch (mt) { + case SCORE_TOTAL: + cutoff = prn() * simulation::material_xs.total; + break; + case SCORE_SCATTER: + cutoff = prn() * (simulation::material_xs.total - + simulation::material_xs.absorption); + break; + case SCORE_FISSION: + cutoff = prn() * simulation::material_xs.fission; + break; + } + + // Get pointers to nuclide/density arrays + int* nuclides; + double* densities; + int n; + openmc_material_get_densities(p->material, &nuclides, &densities, &n); + + *i_nuc_mat = 0; + double prob = 0.0; + while (prob < cutoff) { + // Check to make sure that a nuclide was sampled + if (*i_nuc_mat > n) { + p->write_restart(); + fatal_error("Did not sample any nuclide during collision."); + } + + // Find atom density + // TODO: off-by-one + *i_nuclide = nuclides[*i_nuc_mat] - 1; + double atom_density = densities[*i_nuc_mat]; + + // Determine microscopic cross section + double sigma; + switch (mt) { + case SCORE_TOTAL: + sigma = atom_density * simulation::micro_xs[*i_nuclide].total; + break; + case SCORE_SCATTER: + sigma = atom_density * (simulation::micro_xs[*i_nuclide].total - + simulation::micro_xs[*i_nuclide].absorption); + break; + case SCORE_FISSION: + sigma = atom_density * simulation::micro_xs[*i_nuclide].fission; + break; + } + + // Increment probability to compare to cutoff + prob += sigma; + + ++(*i_nuc_mat); + } +} + +// TODO: Finish converting photon physics functions + +// void sample_element(Particle* p) +// { +// associate (mat => materials(p->material)) +// // Sample cumulative distribution function +// cutoff = prn() * simulation::material_xs.total + +// i = 0 +// prob = 0.0 +// do while (prob < cutoff) +// i = i + 1 + +// // Check to make sure that a nuclide was sampled +// if (i > mat % n_nuclides) { +// particle_write_restart(p) +// fatal_error("Did not sample any element during collision.") +// } + +// // Find atom density +// i_element = mat % element(i) +// atom_density = mat % atom_density(i) + +// // Determine microscopic cross section +// sigma = atom_density * micro_photon_xs(i_element) % total + +// // Increment probability to compare to cutoff +// prob = prob + sigma +// end do +// end associate +// } + +Reaction* sample_fission(int i_nuclide, double E) +{ + // Get pointer to nuclide + const auto& nuc {data::nuclides[i_nuclide]}; + + // If we're in the URR, by default use the first fission reaction. We also + // default to the first reaction if we know that there are no partial fission + // reactions + if (simulation::micro_xs[i_nuclide].use_ptable || !nuc->has_partial_fission_) { + return nuc->fission_rx_[0]; + } + + // Check to see if we are in a windowed multipole range. WMP only supports + // the first fission reaction. + if (nuclide_wmp_present(i_nuclide)) { + if (E >= nuclide_wmp_emin(i_nuclide) && E <= nuclide_wmp_emax(i_nuclide)) { + return nuc->fission_rx_[0]; + } + } + + // Get grid index and interpolatoin factor and sample fission cdf + int i_temp = simulation::micro_xs[i_nuclide].index_temp - 1; + int i_grid = simulation::micro_xs[i_nuclide].index_grid; + double f = simulation::micro_xs[i_nuclide].interp_factor; + double cutoff = prn() * simulation::micro_xs[i_nuclide].fission; + double prob = 0.0; + + // Loop through each partial fission reaction type + for (auto& rx : nuc->fission_rx_) { + // if energy is below threshold for this reaction, skip it + int threshold = rx->xs_[i_temp].threshold; + if (i_grid < threshold) continue; + + // add to cumulative probability + prob += (1.0 - f) * rx->xs_[i_temp].value[i_grid - threshold] + + f*rx->xs_[i_temp].value[i_grid - threshold + 1]; + + // Create fission bank sites if fission occurs + if (prob > cutoff) return rx; + } +} + +void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product) +{ + // Get grid index and interpolation factor and sample photon production cdf + // TODO: off-by-one + int i_temp = simulation::micro_xs[i_nuclide].index_temp - 1; + int i_grid = simulation::micro_xs[i_nuclide].index_grid; + double f = simulation::micro_xs[i_nuclide].interp_factor; + double cutoff = prn() * simulation::micro_xs[i_nuclide].photon_prod; + double prob = 0.0; + + // Loop through each reaction type + const auto& nuc {data::nuclides[i_nuclide]}; + for (int i = 0; i < nuc->reactions_.size(); ++i) { + const auto& rx = nuc->reactions_[i]; + int threshold = rx->xs_[i_temp].threshold; + + // if energy is below threshold for this reaction, skip it + if (i_grid < threshold) continue; + + // Evaluate neutron cross section + double xs = ((1.0 - f) * rx->xs_[i_temp].value[i_grid - threshold] + + f*(rx->xs_[i_temp].value[i_grid - threshold + 1])); + + for (int j = 0; j < rx->products_.size(); ++j) { + if (rx->products_[j].particle_ == ParticleType::photon) { + // add to cumulative probability + prob += (*rx->products_[j].yield_)(E) * xs; + + *i_rx = i; + *i_product = j; + if (prob > cutoff) return; + } + } + } +} + +void absorption(Particle* p, int i_nuclide) +{ + if (settings::survival_biasing) { + // Determine weight absorbed in survival biasing + p->absorb_wgt = p->wgt * simulation::micro_xs[i_nuclide].absorption / + simulation::micro_xs[i_nuclide].total; + + // Adjust weight of particle by probability of absorption + p->wgt -= p->absorb_wgt; + p->last_wgt = p->wgt; + + // Score implicit absorption estimate of keff + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + global_tally_absorption += p->absorb_wgt * simulation::micro_xs[ + i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption; + } + } else { + // See if disappearance reaction happens + if (simulation::micro_xs[i_nuclide].absorption > + prn() * simulation::micro_xs[i_nuclide].total) { + // Score absorption estimate of keff + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + global_tally_absorption += p->wgt * simulation::micro_xs[ + i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption; + } + + p->alive = false; + p->event = EVENT_ABSORB; + p->event_MT = N_DISAPPEAR; + } + } +} + +void scatter(Particle* p, int i_nuclide, int i_nuc_mat) +{ + // copy incoming direction + Direction u_old {p->coord[0].uvw}; + + // Get pointer to nuclide and grid index/interpolation factor + const auto& nuc {data::nuclides[i_nuclide]}; + const auto& micro {simulation::micro_xs[i_nuclide]}; + int i_temp = micro.index_temp - 1; + int i_grid = micro.index_grid - 1; + double f = micro.interp_factor; + + // For tallying purposes, this routine might be called directly. In that + // case, we need to sample a reaction via the cutoff variable + double cutoff = prn() * (micro.total - micro.absorption); + bool sampled = false; + + // Calculate elastic cross section if it wasn't precalculated + if (micro.elastic == CACHE_INVALID) { + nuc->calculate_elastic_xs(i_nuclide); + } + + double prob = micro.elastic - micro.thermal; + if (prob > cutoff) { + // ======================================================================= + // NON-S(A,B) ELASTIC SCATTERING + + // Determine temperature + double kT = nuclide_wmp_present(i_nuclide) ? + p->sqrtkT*p->sqrtkT : nuc->kTs_[i_temp]; + + // Perform collision physics for elastic scattering + elastic_scatter(i_nuclide, nuc->reactions_[0].get(), kT, + &p->E, p->coord[0].uvw, &p->mu, &p->wgt); + + p->event_MT = ELASTIC; + sampled = true; + } + + prob = micro.elastic; + if (prob > cutoff && !sampled) { + // ======================================================================= + // S(A,B) SCATTERING + + sab_scatter(i_nuclide, micro.index_sab, &p->E, p->coord[0].uvw, &p->mu); + + p->event_MT = ELASTIC; + sampled = true; + } + + if (!sampled) { + // ======================================================================= + // INELASTIC SCATTERING + + int j = 0; + int i; + while (prob < cutoff) { + i = nuc->index_inelastic_scatter_[j]; + ++j; + + // Check to make sure inelastic scattering reaction sampled + if (i >= nuc->reactions_.size()) { + p->write_restart(); + fatal_error("Did not sample any reaction for nuclide " + nuc->name_); + } + + // if energy is below threshold for this reaction, skip it + const auto& xs {nuc->reactions_[i]->xs_[i_temp]}; + int threshold = xs.threshold - 1; + if (i_grid < threshold) continue; + + // add to cumulative probability + prob += (1.0 - f)*xs.value[i_grid - threshold] + + f*xs.value[i_grid - threshold + 1]; + } + + // Perform collision physics for inelastic scattering + const auto& rx {nuc->reactions_[i]}; + inelastic_scatter(nuc.get(), rx.get(), p); + p->event_MT = rx->mt_; + } + + // Set event component + p->event = EVENT_SCATTER; + + // Sample new outgoing angle for isotropic-in-lab scattering + if (material_isotropic(p->material, i_nuc_mat)) { + // Sample isotropic-in-lab outgoing direction + double mu = 2.0*prn() - 1.0; + double phi = 2.0*PI*prn(); + Direction u_new; + u_new.x = mu; + u_new.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); + u_new.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); + + p->mu = u_old.dot(u_new); + + // Change direction of particle + p->coord[0].uvw[0] = u_new.x; + p->coord[0].uvw[1] = u_new.y; + p->coord[0].uvw[2] = u_new.z; + } +} + +void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, + double* uvw, double* mu_lab, double* wgt) +{ + // get pointer to nuclide + const auto& nuc {data::nuclides[i_nuclide]}; + + double vel = std::sqrt(*E); + double awr = nuc->awr_; + + // Neutron velocity in LAB + Direction u {uvw}; + Direction v_n = vel*u; + + // Sample velocity of target nucleus + Direction v_t {}; + if (!simulation::micro_xs[i_nuclide].use_ptable) { + v_t = sample_target_velocity(nuc.get(), *E, u, v_n, + simulation::micro_xs[i_nuclide].elastic, kT, wgt); + } + + // Velocity of center-of-mass + Direction v_cm = (v_n + awr*v_t)/(awr + 1.0); + + // Transform to CM frame + v_n -= v_cm; + + // Find speed of neutron in CM + vel = v_n.norm(); + + // Sample scattering angle, checking if it is an ncorrelated angle-energy + // distribution + double mu_cm; + auto& d = rx->products_[0].distribution_[0]; + auto d_ = dynamic_cast(d.get()); + if (d_) { + mu_cm = d_->angle().sample(*E); + } else { + mu_cm = 2.0*prn() - 1.0; + } + + // Determine direction cosines in CM + Direction u_cm = v_n/vel; + + // Rotate neutron velocity vector to new angle -- note that the speed of the + // neutron in CM does not change in elastic scattering. However, the speed + // will change when we convert back to LAB + v_n = vel * rotate_angle(u_cm, mu_cm, nullptr); + + // Transform back to LAB frame + v_n += v_cm; + + *E = v_n.dot(v_n); + vel = std::sqrt(*E); + + // compute cosine of scattering angle in LAB frame by taking dot product of + // neutron's pre- and post-collision angle + *mu_lab = u.dot(v_n) / vel; + + // Set energy and direction of particle in LAB frame + u = v_n / vel; + uvw[0] = u.x; + uvw[1] = u.y; + uvw[2] = u.z; + + // Because of floating-point roundoff, it may be possible for mu_lab to be + // outside of the range [-1,1). In these cases, we just set mu_lab to exactly + // -1 or 1 + if (std::abs(*mu_lab) > 1.0) *mu_lab = std::copysign(1.0, *mu_lab); +} + +// void sab_scatter(int i_nuclide, int i_sab, double* E, Direction* u, double* mu) +// { +// // Sample from C++ side +// ptr = C_LOC(micro_xs(i_nuclide)) +// sab_tables(i_sab) % sample(ptr, E, E_out, mu) + +// // Set energy to outgoing, change direction of particle +// E = E_out +// uvw = rotate_angle(uvw, mu) +// } + +Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, + Direction v_neut, double xs_eff, double kT, double* wgt) +{ + // check if nuclide is a resonant scatterer + int sampling_method; + if (nuc->resonant_) { + + // sampling method to use + sampling_method = settings::res_scat_method; + + // upper resonance scattering energy bound (target is at rest above this E) + if (E > settings::res_scat_energy_max) { + return {}; + + // lower resonance scattering energy bound (should be no resonances below) + } else if (E < settings::res_scat_energy_min) { + sampling_method = RES_SCAT_CXS; + } + + // otherwise, use free gas model + } else { + if (E >= FREE_GAS_THRESHOLD * kT && nuc->awr_ > 1.0) { + return {}; + } else { + sampling_method = RES_SCAT_CXS; + } + } + + // use appropriate target velocity sampling method + switch (sampling_method) { + case RES_SCAT_CXS: + + // sample target velocity with the constant cross section (cxs) approx. + return sample_cxs_target_velocity(nuc->awr_, E, u, kT); + + case RES_SCAT_WCM: { + // sample target velocity with the constant cross section (cxs) approx. + Direction v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT); + + // adjust weight as prescribed by the weight correction method (wcm) + Direction v_rel = v_neut - v_target; + double E_rel = v_rel.dot(v_rel); + double xs_0K = nuc->elastic_xs_0K(E_rel); + *wgt *= xs_0K / xs_eff; + return v_target; + } + + case RES_SCAT_DBRC: + case RES_SCAT_ARES: { + double E_red = std::sqrt(nuc->awr_ * E / kT); + double E_low = std::pow(std::max(0.0, E_red - 4.0), 2) * kT / nuc->awr_; + double E_up = (E_red + 4.0)*(E_red + 4.0) * kT / nuc->awr_; + + // find lower and upper energy bound indices + // lower index + int i_E_low; + if (E_low < nuc->energy_0K_.front()) { + i_E_low = 0; + } else if (E_low > nuc->energy_0K_.back()) { + i_E_low = nuc->energy_0K_.size() - 2; + } else { + i_E_low = lower_bound_index(nuc->energy_0K_.begin(), + nuc->energy_0K_.end(), E_low); + } + + // upper index + int i_E_up; + if (E_up < nuc->energy_0K_.front()) { + i_E_up = 0; + } else if (E_up > nuc->energy_0K_.back()) { + i_E_up = nuc->energy_0K_.size() - 2; + } else { + i_E_up = lower_bound_index(nuc->energy_0K_.begin(), + nuc->energy_0K_.end(), E_up); + } + + if (i_E_up == i_E_low) { + // Handle degenerate case -- if the upper/lower bounds occur for the same + // index, then using cxs is probably a good approximation + return sample_cxs_target_velocity(nuc->awr_, E, u, kT); + } + + if (sampling_method == RES_SCAT_DBRC) { + // interpolate xs since we're not exactly at the energy indices + double xs_low = nuc->elastic_0K_[i_E_low]; + double m = (nuc->elastic_0K_[i_E_low + 1] - xs_low) + / (nuc->energy_0K_[i_E_low + 1] - nuc->energy_0K_[i_E_low]); + xs_low += m * (E_low - nuc->energy_0K_[i_E_low]); + double xs_up = nuc->elastic_0K_[i_E_up]; + m = (nuc->elastic_0K_[i_E_up + 1] - xs_up) + / (nuc->energy_0K_[i_E_up + 1] - nuc->energy_0K_[i_E_up]); + xs_up += m * (E_up - nuc->energy_0K_[i_E_up]); + + // get max 0K xs value over range of practical relative energies + double xs_max = *std::max_element(&nuc->elastic_0K_[i_E_low + 1], + &nuc->elastic_0K_[i_E_up + 1]); + xs_max = std::max({xs_low, xs_max, xs_up}); + + while (true) { + double E_rel; + Direction v_target; + while (true) { + // sample target velocity with the constant cross section (cxs) approx. + v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT); + Direction v_rel = v_neut - v_target; + E_rel = v_rel.dot(v_rel); + if (E_rel < E_up) break; + } + + // perform Doppler broadening rejection correction (dbrc) + double xs_0K = nuc->elastic_xs_0K(E_rel); + double R = xs_0K / xs_max; + if (prn() < R) return v_target; + } + + } else if (sampling_method == RES_SCAT_ARES) { + // interpolate xs CDF since we're not exactly at the energy indices + // cdf value at lower bound attainable energy + double m = (nuc->xs_cdf_[i_E_low] - nuc->xs_cdf_[i_E_low - 1]) + / (nuc->energy_0K_[i_E_low + 1] - nuc->energy_0K_[i_E_low]); + double cdf_low = nuc->xs_cdf_[i_E_low - 1] + + m * (E_low - nuc->energy_0K_[i_E_low]); + if (E_low <= nuc->energy_0K_.front()) cdf_low = 0.0; + + // cdf value at upper bound attainable energy + m = (nuc->xs_cdf_[i_E_up] - nuc->xs_cdf_[i_E_up - 1]) + / (nuc->energy_0K_[i_E_up + 1] - nuc->energy_0K_[i_E_up]); + double cdf_up = nuc->xs_cdf_[i_E_up - 1] + + m*(E_up - nuc->energy_0K_[i_E_up]); + + while (true) { + // directly sample Maxwellian + double E_t = -kT * std::log(prn()); + + // sample a relative energy using the xs cdf + double cdf_rel = cdf_low + prn()*(cdf_up - cdf_low); + int i_E_rel = lower_bound_index(&nuc->xs_cdf_[i_E_low-1], + &nuc->xs_cdf_[i_E_up+1], cdf_rel); + double E_rel = nuc->energy_0K_[i_E_low + i_E_rel]; + double m = (nuc->xs_cdf_[i_E_low + i_E_rel] + - nuc->xs_cdf_[i_E_low + i_E_rel - 1]) + / (nuc->energy_0K_[i_E_low + i_E_rel + 1] + - nuc->energy_0K_[i_E_low + i_E_rel]); + E_rel += (cdf_rel - nuc->xs_cdf_[i_E_low + i_E_rel - 1]) / m; + + // perform rejection sampling on cosine between + // neutron and target velocities + double mu = (E_t + nuc->awr_ * (E - E_rel)) / + (2.0 * std::sqrt(nuc->awr_ * E * E_t)); + + if (std::abs(mu) < 1.0) { + // set and accept target velocity + E_t /= nuc->awr_; + return std::sqrt(E_t) * rotate_angle(u, mu, nullptr); + } + } + } + } // case RES_SCAT_ARES, RES_SCAT_DBRC + } // switch (sampling_method) +} + +Direction +sample_cxs_target_velocity(double awr, double E, Direction u, double kT) +{ + double beta_vn = std::sqrt(awr * E / kT); + double alpha = 1.0/(1.0 + std::sqrt(PI)*beta_vn/2.0); + + double beta_vt_sq; + double mu; + while (true) { + // Sample two random numbers + double r1 = prn(); + double r2 = prn(); + + if (prn() < alpha) { + // With probability alpha, we sample the distribution p(y) = + // y*e^(-y). This can be done with sampling scheme C45 frmo the Monte + // Carlo sampler + + beta_vt_sq = -std::log(r1*r2); + + } else { + // With probability 1-alpha, we sample the distribution p(y) = y^2 * + // e^(-y^2). This can be done with sampling scheme C61 from the Monte + // Carlo sampler + + double c = std::cos(PI/2.0 * prn()); + beta_vt_sq = -std::log(r1) - std::log(r2)*c*c; + } + + // Determine beta * vt + double beta_vt = std::sqrt(beta_vt_sq); + + // Sample cosine of angle between neutron and target velocity + mu = 2.0*prn() - 1.0; + + // Determine rejection probability + double accept_prob = std::sqrt(beta_vn*beta_vn + beta_vt_sq - + 2*beta_vn*beta_vt*mu) / (beta_vn + beta_vt); + + // Perform rejection sampling on vt and mu + if (prn() < accept_prob) break; + } + + // Determine speed of target nucleus + double vt = std::sqrt(beta_vt_sq*kT/awr); + + // Determine velocity vector of target nucleus based on neutron's velocity + // and the sampled angle between them + return vt * rotate_angle(u, mu, nullptr); +} + +void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site) +{ + // Sample cosine of angle -- fission neutrons are always emitted + // isotropically. Sometimes in ACE data, fission reactions actually have + // an angular distribution listed, but for those that do, it's simply just + // a uniform distribution in mu + double mu = 2.0 * prn() - 1.0; + + // Sample azimuthal angle uniformly in [0,2*pi) + double phi = 2.0*PI*prn(); + site->uvw[0] = mu; + site->uvw[1] = std::sqrt(1.0 - mu*mu) * std::cos(phi); + site->uvw[2] = std::sqrt(1.0 - mu*mu) * std::sin(phi); + + // Determine total nu, delayed nu, and delayed neutron fraction + const auto& nuc {data::nuclides[i_nuclide]}; + double nu_t = nuc->nu(E_in, Nuclide::EmissionMode::total); + double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed); + double beta = nu_d / nu_t; + + if (prn() < beta) { + // ==================================================================== + // DELAYED NEUTRON SAMPLED + + // sampled delayed precursor group + double xi = prn()*nu_d; + double prob = 0.0; + int group; + for (group = 1; group < nuc->n_precursor_; ++group) { + // determine delayed neutron precursor yield for group j + double yield = (*rx->products_[group].yield_)(E_in); + + // Check if this group is sampled + prob += yield; + if (xi < prob) break; + } + + // if the sum of the probabilities is slightly less than one and the + // random number is greater, j will be greater than nuc % + // n_precursor -- check for this condition + group = std::min(group, nuc->n_precursor_); + + // set the delayed group for the particle born from fission + site->delayed_group = group; + + int n_sample = 0; + while (true) { + // sample from energy/angle distribution -- note that mu has already been + // sampled above and doesn't need to be resampled + rx->products_[group].sample(E_in, site->E, mu); + + // resample if energy is greater than maximum neutron energy + constexpr int neutron = static_cast(ParticleType::neutron); + if (site->E < data::energy_max[neutron]) break; + + // check for large number of resamples + ++n_sample; + if (n_sample == MAX_SAMPLE) { + // particle_write_restart(p) + fatal_error("Resampled energy distribution maximum number of times " + "for nuclide " + nuc->name_); + } + } + + } else { + // ==================================================================== + // PROMPT NEUTRON SAMPLED + + // set the delayed group for the particle born from fission to 0 + site->delayed_group = 0; + + // sample from prompt neutron energy distribution + int n_sample = 0; + while (true) { + rx->products_[0].sample(E_in, site->E, mu); + + // resample if energy is greater than maximum neutron energy + constexpr int neutron = static_cast(ParticleType::neutron); + if (site->E < data::energy_max[neutron]) break; + + // check for large number of resamples + ++n_sample; + if (n_sample == MAX_SAMPLE) { + // particle_write_restart(p) + fatal_error("Resampled energy distribution maximum number of times " + "for nuclide " + nuc->name_); + } + } + } +} + +void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) +{ + // copy energy of neutron + double E_in = p->E; + + // sample outgoing energy and scattering cosine + double E; + double mu; + rx->products_[0].sample(E_in, E, mu); + + // if scattering system is in center-of-mass, transfer cosine of scattering + // angle and outgoing energy from CM to LAB + if (rx->scatter_in_cm_) { + double E_cm = E; + + // determine outgoing energy in lab + double A = nuc->awr_; + E = E_cm + (E_in + 2.0*mu*(A + 1.0) * std::sqrt(E_in*E_cm)) + / ((A + 1.0)*(A + 1.0)); + + // determine outgoing angle in lab + mu = mu*std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E); + } + + // Because of floating-point roundoff, it may be possible for mu to be + // outside of the range [-1,1). In these cases, we just set mu to exactly -1 + // or 1 + if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu); + + // Set outgoing energy and scattering angle + p->E = E; + p->mu = mu; + + // change direction of particle + rotate_angle_c(p->coord[0].uvw, mu, nullptr); + + // evaluate yield + double yield = (*rx->products_[0].yield_)(E_in); + if (std::floor(yield) == yield) { + // If yield is integral, create exactly that many secondary particles + for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { + int neutron = static_cast(ParticleType::neutron); + p->create_secondary(p->coord[0].uvw, p->E, neutron, true); + } + } else { + // Otherwise, change weight of particle based on yield + p->wgt *= yield; + } +} + +void sample_secondary_photons(Particle* p, int i_nuclide) +{ + // Sample the number of photons produced + double y_t = p->wgt * simulation::micro_xs[i_nuclide].photon_prod / + simulation::micro_xs[i_nuclide].total; + int y = static_cast(y_t); + if (prn() <= y_t - y) ++y; + + // Sample each secondary photon + for (int i = 0; i < y; ++i) { + // Sample the reaction and product + int i_rx; + int i_product; + sample_photon_product(i_nuclide, p->E, &i_rx, &i_product); + + // Sample the outgoing energy and angle + auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; + double E; + double mu; + rx->products_[i_product].sample(p->E, E, mu); + + // Sample the new direction + double uvw[3]; + std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + rotate_angle_c(uvw, mu, nullptr); + + // Create the secondary photon + int photon = static_cast(ParticleType::photon); + p->create_secondary(uvw, E, photon, true); + } +} + +} // namespace openmc diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index a4795f82a2..202e5abe30 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -5,6 +5,7 @@ #include "xtensor/xarray.hpp" +#include "openmc/bank.h" #include "openmc/constants.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" @@ -47,12 +48,8 @@ sample_reaction(Particle* p, const MaterialMacroXS* material_xs) if (model::materials[p->material - 1]->fissionable) { if (settings::run_mode == RUN_MODE_EIGENVALUE) { - Bank* result_bank; - int64_t result_bank_size; - // Get pointer to fission bank from Fortran side - openmc_fission_bank(&result_bank, &result_bank_size); - create_fission_sites(p, result_bank, &simulation::n_bank, result_bank_size, - material_xs); + create_fission_sites(p, simulation::fission_bank.data(), &simulation::n_bank, + simulation::fission_bank.size(), material_xs); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { create_fission_sites(p, p->secondary_bank, &(p->n_secondary), diff --git a/src/position.cpp b/src/position.cpp index 18ec33e803..ecf5be631f 100644 --- a/src/position.cpp +++ b/src/position.cpp @@ -60,4 +60,22 @@ Position::operator*=(double v) return *this; } +Position& +Position::operator/=(Position other) +{ + x /= other.x; + y /= other.y; + z /= other.z; + return *this; +} + +Position& +Position::operator/=(double v) +{ + x /= v; + y /= v; + z /= v; + return *this; +} + } // namespace openmc diff --git a/src/reaction.cpp b/src/reaction.cpp index 94dc7f4377..918d61fea0 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -3,6 +3,7 @@ #include #include // for move +#include "openmc/constants.h" #include "openmc/hdf5_interface.h" #include "openmc/endf.h" #include "openmc/random_lcg.h" @@ -10,6 +11,10 @@ namespace openmc { +//============================================================================== +// Reaction implementation +//============================================================================== + Reaction::Reaction(hid_t group, const std::vector& temperatures) { read_attribute(group, "Q_value", q_value_); @@ -77,16 +82,198 @@ Reaction::Reaction(hid_t group, const std::vector& temperatures) } //============================================================================== -// Fortran compatibility functions +// Non-member functions //============================================================================== -Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n) +std::string reaction_name(int mt) { - std::vector temps {temperatures, temperatures + n}; - return new Reaction{group, temps}; + if (mt == SCORE_FLUX) { + return "flux"; + } else if (mt == SCORE_TOTAL) { + return "total"; + } else if (mt == SCORE_SCATTER) { + return "scatter"; + } else if (mt == SCORE_NU_SCATTER) { + return "nu-scatter"; + } else if (mt == SCORE_ABSORPTION) { + return "absorption"; + } else if (mt == SCORE_FISSION) { + return "fission"; + } else if (mt == SCORE_NU_FISSION) { + return "nu-fission"; + } else if (mt == SCORE_DECAY_RATE) { + return "decay-rate"; + } else if (mt == SCORE_DELAYED_NU_FISSION) { + return "delayed-nu-fission"; + } else if (mt == SCORE_PROMPT_NU_FISSION) { + return "prompt-nu-fission"; + } else if (mt == SCORE_KAPPA_FISSION) { + return "kappa-fission"; + } else if (mt == SCORE_CURRENT) { + return "current"; + } else if (mt == SCORE_EVENTS) { + return "events"; + } else if (mt == SCORE_INVERSE_VELOCITY) { + return "inverse-velocity"; + } else if (mt == SCORE_FISS_Q_PROMPT) { + return "fission-q-prompt"; + } else if (mt == SCORE_FISS_Q_RECOV) { + return "fission-q-recoverable"; + + // Normal ENDF-based reactions + } else if (mt == TOTAL_XS) { + return "(n,total)"; + } else if (mt == ELASTIC) { + return "(n,elastic)"; + } else if (mt == N_LEVEL) { + return "(n,level)"; + } else if (mt == N_2ND) { + return "(n,2nd)"; + } else if (mt == N_2N) { + return "(n,2n)"; + } else if (mt == N_3N) { + return "(n,3n)"; + } else if (mt == N_FISSION) { + return "(n,fission)"; + } else if (mt == N_F) { + return "(n,f)"; + } else if (mt == N_NF) { + return "(n,nf)"; + } else if (mt == N_2NF) { + return "(n,2nf)"; + } else if (mt == N_NA) { + return "(n,na)"; + } else if (mt == N_N3A) { + return "(n,n3a)"; + } else if (mt == N_2NA) { + return "(n,2na)"; + } else if (mt == N_3NA) { + return "(n,3na)"; + } else if (mt == N_NP) { + return "(n,np)"; + } else if (mt == N_N2A) { + return "(n,n2a)"; + } else if (mt == N_2N2A) { + return "(n,2n2a)"; + } else if (mt == N_ND) { + return "(n,nd)"; + } else if (mt == N_NT) { + return "(n,nt)"; + } else if (mt == N_N3HE) { + return "(n,nHe-3)"; + } else if (mt == N_ND2A) { + return "(n,nd2a)"; + } else if (mt == N_NT2A) { + return "(n,nt2a)"; + } else if (mt == N_4N) { + return "(n,4n)"; + } else if (mt == N_3NF) { + return "(n,3nf)"; + } else if (mt == N_2NP) { + return "(n,2np)"; + } else if (mt == N_3NP) { + return "(n,3np)"; + } else if (mt == N_N2P) { + return "(n,n2p)"; + } else if (mt == N_NPA) { + return "(n,npa)"; + } else if (N_N1 <= mt && mt <= N_N40) { + return "(n,n" + std::to_string(mt-50) + ")"; + } else if (mt == N_NC) { + return "(n,nc)"; + } else if (mt == N_DISAPPEAR) { + return "(n,disappear)"; + } else if (mt == N_GAMMA) { + return "(n,gamma)"; + } else if (mt == N_P) { + return "(n,p)"; + } else if (mt == N_D) { + return "(n,d)"; + } else if (mt == N_T) { + return "(n,t)"; + } else if (mt == N_3HE) { + return "(n,3He)"; + } else if (mt == N_A) { + return "(n,a)"; + } else if (mt == N_2A) { + return "(n,2a)"; + } else if (mt == N_3A) { + return "(n,3a)"; + } else if (mt == N_2P) { + return "(n,2p)"; + } else if (mt == N_PA) { + return "(n,pa)"; + } else if (mt == N_T2A) { + return "(n,t2a)"; + } else if (mt == N_D2A) { + return "(n,d2a)"; + } else if (mt == N_PD) { + return "(n,pd)"; + } else if (mt == N_PT) { + return "(n,pt)"; + } else if (mt == N_DA) { + return "(n,da)"; + } else if (mt == 201) { + return "(n,Xn)"; + } else if (mt == 202) { + return "(n,Xgamma)"; + } else if (mt == 203) { + return "(n,Xp)"; + } else if (mt == 204) { + return "(n,Xd)"; + } else if (mt == 205) { + return "(n,Xt)"; + } else if (mt == 206) { + return "(n,X3He)"; + } else if (mt == 207) { + return "(n,Xa)"; + } else if (mt == 444) { + return "(damage)"; + } else if (mt == COHERENT) { + return "coherent scatter"; + } else if (mt == INCOHERENT) { + return "incoherent scatter"; + } else if (mt == PAIR_PROD_ELEC) { + return "pair production, electron"; + } else if (mt == PAIR_PROD) { + return "pair production"; + } else if (mt == PAIR_PROD_NUC) { + return "pair production, nuclear"; + } else if (mt == PHOTOELECTRIC) { + return "photoelectric"; + } else if (534 <= mt && mt <= 572) { + std::stringstream name; + name << "photoelectric, " << SUBSHELLS[mt - 534] << " subshell"; + return name.str(); + } else if (600 <= mt && mt <= 648) { + return "(n,p" + std::to_string(mt-600) + ")"; + } else if (mt == 649) { + return "(n,pc)"; + } else if (650 <= mt && mt <= 698) { + return "(n,d" + std::to_string(mt-650) + ")"; + } else if (mt == 699) { + return "(n,dc)"; + } else if (700 <= mt && mt <= 748) { + return "(n,t" + std::to_string(mt-700) + ")"; + } else if (mt == 749) { + return "(n,tc)"; + } else if (750 <= mt && mt <= 798) { + return "(n,3He" + std::to_string(mt-750) + ")"; + } else if (mt == 799) { + return "(n,3Hec)"; + } else if (800 <= mt && mt <= 848) { + return "(n,a" + std::to_string(mt-800) + ")"; + } else if (mt == 849) { + return "(n,ac)"; + } else { + return "MT=" + std::to_string(mt); + } } -void reaction_delete(Reaction* rx) { delete rx; } + +//============================================================================== +// Fortran compatibility functions +//============================================================================== int reaction_mt(Reaction* rx) { return rx->mt_; } diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index 0762846efc..9d57db5ccf 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -22,7 +22,7 @@ module reaction_header logical(C_BOOL) :: scatter_in_cm ! scattering system in center-of-mass? logical(C_BOOL) :: redundant ! redundant reaction? contains - procedure :: from_hdf5 + procedure :: init procedure :: mt_ procedure :: q_value_ procedure :: scatter_in_cm_ @@ -40,11 +40,10 @@ module reaction_header end type Reaction interface - function reaction_from_hdf5(group, temperatures, n) result(ptr) bind(C) - import C_PTR, HID_T, C_INT - integer(HID_T), value :: group - integer(C_INT), intent(in) :: temperatures - integer(C_INT), value :: n + function nuclide_reaction(nuc_ptr, i_rx) result(ptr) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: nuc_ptr + integer(C_INT), value :: i_rx type(C_PTR) :: ptr end function @@ -148,27 +147,17 @@ module reaction_header contains - subroutine from_hdf5(this, group_id, temperatures) + subroutine init(this, nuc_ptr, i_rx) class(Reaction), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - type(VectorInt), intent(in) :: temperatures + type(C_PTR), intent(in) :: nuc_ptr + integer(C_INT), intent(in) :: i_rx - integer(C_INT) :: dummy - integer(C_INT) :: n - - n = temperatures % size() - if (n > 0) then - this % ptr = reaction_from_hdf5(group_id, temperatures % data(1), n) - else - ! In this case, temperatures % data(1) doesn't exist, so we just pass a - ! dummy value - this % ptr = reaction_from_hdf5(group_id, dummy, n) - end if + this % ptr = nuclide_reaction(nuc_ptr, i_rx) this % MT = reaction_mt(this % ptr) this % Q_value = reaction_q_value(this % ptr) this % scatter_in_cm = reaction_scatter_in_cm(this % ptr) this % redundant = reaction_redundant(this % ptr) - end subroutine from_hdf5 + end subroutine function mt_(this) result(mt) class(Reaction), intent(in) :: this diff --git a/src/settings.cpp b/src/settings.cpp index 36657c6c62..303a35ddd9 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -806,18 +806,6 @@ void read_settings_xml() //============================================================================== extern "C" { - bool res_scat_nuclides_empty() { - return settings::res_scat_nuclides.empty(); - } - - int res_scat_nuclides_size() { - return settings::res_scat_nuclides.size(); - } - - bool res_scat_nuclides_cmp(int i, const char* name) { - return settings::res_scat_nuclides[i - 1] == name; - } - const char* path_cross_sections_c() { return settings::path_cross_sections.c_str(); } diff --git a/src/simulation.F90 b/src/simulation.F90 index 3eabbfecb4..f00b58ed1c 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -2,20 +2,9 @@ module simulation use, intrinsic :: ISO_C_BINDING -#ifdef _OPENMP - use omp_lib -#endif - - use bank_header, only: source_bank - use constants, only: ZERO - use error, only: fatal_error use material_header, only: n_materials, materials - use message_passing use nuclide_header, only: micro_xs, n_nuclides use photon_header, only: micro_photon_xs, n_elements - use settings - use simulation_header - use tally_header use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer implicit none @@ -69,56 +58,4 @@ contains end subroutine -!=============================================================================== -! ALLOCATE_BANKS allocates memory for the fission and source banks -!=============================================================================== - - subroutine allocate_banks() bind(C) - - integer :: alloc_err ! allocation error code - - ! Allocate source bank - if (allocated(source_bank)) deallocate(source_bank) - allocate(source_bank(work), STAT=alloc_err) - - ! Check for allocation errors - if (alloc_err /= 0) then - call fatal_error("Failed to allocate source bank.") - end if - - if (run_mode == MODE_EIGENVALUE) then - -#ifdef _OPENMP - ! If OpenMP is being used, each thread needs its own private fission - ! bank. Since the private fission banks need to be combined at the end of - ! a generation, there is also a 'master_fission_bank' that is used to - ! collect the sites from each thread. - - n_threads = omp_get_max_threads() - -!$omp parallel - thread_id = omp_get_thread_num() - - if (allocated(fission_bank)) deallocate(fission_bank) - if (thread_id == 0) then - allocate(fission_bank(3*work)) - else - allocate(fission_bank(3*work/n_threads)) - end if -!$omp end parallel - if (allocated(master_fission_bank)) deallocate(master_fission_bank) - allocate(master_fission_bank(3*work), STAT=alloc_err) -#else - if (allocated(fission_bank)) deallocate(fission_bank) - allocate(fission_bank(3*work), STAT=alloc_err) -#endif - - ! Check for allocation errors - if (alloc_err /= 0) then - call fatal_error("Failed to allocate fission bank.") - end if - end if - - end subroutine allocate_banks - end module simulation diff --git a/src/simulation.cpp b/src/simulation.cpp index 50bf12521c..8459c9c71f 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,10 +1,12 @@ #include "openmc/simulation.h" +#include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" #include "openmc/message_passing.h" +#include "openmc/nuclide.h" #include "openmc/output.h" #include "openmc/particle.h" #include "openmc/random_lcg.h" @@ -15,6 +17,8 @@ #include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" +#include + #include #include @@ -24,14 +28,12 @@ namespace openmc { extern "C" bool cmfd_on; extern "C" void accumulate_tallies(); -extern "C" void allocate_banks(); extern "C" void allocate_tally_results(); extern "C" void check_triggers(); extern "C" void cmfd_init_batch(); extern "C" void cmfd_tally_init(); extern "C" void execute_cmfd(); extern "C" void init_tally_routines(); -extern "C" void join_bank_from_threads(); extern "C" void load_state_point(); extern "C" void print_batch_keff(); extern "C" void print_columns(); @@ -99,6 +101,7 @@ int openmc_simulation_init() // Call Fortran initialization simulation_init_f(); + set_micro_xs(); // Reset global variables -- this is done before loading state point (as that // will potentially populate k_generation and entropy) @@ -280,6 +283,36 @@ int thread_id; //!< ID of a given thread // Non-member functions //============================================================================== +void allocate_banks() +{ + // Allocate source bank + simulation::source_bank.resize(simulation::work); + + if (settings::run_mode == RUN_MODE_EIGENVALUE) { +#ifdef _OPENMP + // If OpenMP is being used, each thread needs its own private fission + // bank. Since the private fission banks need to be combined at the end of + // a generation, there is also a 'master_fission_bank' that is used to + // collect the sites from each thread. + + simulation::n_threads = omp_get_max_threads(); + +#pragma omp parallel + { + simulation::thread_id = omp_get_thread_num(); + if (simulation::thread_id == 0) { + simulation::fission_bank.resize(3*simulation::work); + } else { + simulation::fission_bank.resize(3*simulation::work / simulation::n_threads); + } + } + simulation::master_fission_bank.resize(3*simulation::work); +#else + simulation::fission_bank.resize(3*simulation::work); +#endif + } +} + void initialize_batch() { // Increment current batch @@ -453,13 +486,8 @@ void finalize_generation() void initialize_history(Particle* p, int64_t index_source) { - // Get pointer to source bank - Bank* source_bank; - int64_t n; - openmc_source_bank(&source_bank, &n); - // set defaults - p->from_source(&source_bank[index_source - 1]); + p->from_source(&simulation::source_bank[index_source - 1]); // set identifier for particle p->id = simulation::work_index[mpi::rank] + index_source; diff --git a/src/source.cpp b/src/source.cpp index 6d071c0c7a..bfa1198231 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -5,6 +5,7 @@ #include "xtensor/xadapt.hpp" +#include "openmc/bank.h" #include "openmc/cell.h" #include "openmc/error.h" #include "openmc/file_utils.h" @@ -242,11 +243,6 @@ 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 (settings::path_source != "") { // Read the source from a binary file instead of sampling from some // assumed source distribution @@ -268,7 +264,7 @@ void initialize_source() } // Read in the source bank - read_source_bank(file_id, source_bank); + read_source_bank(file_id); // Close file file_close(file_id); @@ -282,7 +278,7 @@ void initialize_source() set_particle_seed(id); // sample external source distribution - source_bank[i] = sample_external_source(); + simulation::source_bank[i] = sample_external_source(); } } @@ -291,7 +287,7 @@ void initialize_source() write_message("Writing out initial source...", 5); std::string filename = settings::path_output + "initial_source.h5"; hid_t file_id = file_open(filename, 'w', true); - write_source_bank(file_id, source_bank); + write_source_bank(file_id); file_close(file_id); } } @@ -354,11 +350,6 @@ extern "C" double total_source_strength() void fill_source_bank_fixedsource() { if (settings::path_source.empty()) { - // Get pointer to source bank - Bank* source_bank; - int64_t n; - openmc_source_bank(&source_bank, &n); - for (int64_t i = 0; i < simulation::work; ++i) { // initialize random number seed int64_t id = (simulation::total_gen + overall_generation()) * @@ -366,7 +357,7 @@ void fill_source_bank_fixedsource() set_particle_seed(id); // sample external source distribution - source_bank[i] = sample_external_source(); + simulation::source_bank[i] = sample_external_source(); } } } diff --git a/src/state_point.F90 b/src/state_point.F90 index 2692856dc0..9172f1fba0 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -35,16 +35,14 @@ module state_point implicit none interface - subroutine write_source_bank(group_id, bank_) bind(C) - import HID_T, C_INT64_T, Bank + subroutine write_source_bank(group_id) bind(C) + import HID_T integer(HID_T), value :: group_id - type(Bank), intent(in) :: bank_(*) end subroutine write_source_bank - subroutine read_source_bank(group_id, bank_) bind(C) - import HID_T, C_INT64_T, Bank + subroutine read_source_bank(group_id) bind(C) + import HID_T integer(HID_T), value :: group_id - type(Bank), intent(out) :: bank_(*) end subroutine read_source_bank end interface @@ -437,7 +435,7 @@ contains if (master .or. parallel) then file_id = file_open(filename_, 'a', parallel=.true.) end if - call write_source_bank(file_id, source_bank) + call write_source_bank(file_id) if (master .or. parallel) call file_close(file_id) end if end function openmc_statepoint_write @@ -634,8 +632,8 @@ contains file_id = file_open(path_source_point, 'r', parallel=.true.) end if - ! Write out source - call read_source_bank(file_id, source_bank) + ! Read source + call read_source_bank(file_id) end if diff --git a/src/state_point.cpp b/src/state_point.cpp index fdf700c444..a808e3ef59 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -8,6 +8,7 @@ #include "xtensor/xbuilder.hpp" // for empty_like #include "xtensor/xview.hpp" +#include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/eigenvalue.h" @@ -70,16 +71,13 @@ write_source_point(const char* filename) } // Get pointer to source bank and write to file - Bank* source_bank; - int64_t n; - openmc_source_bank(&source_bank, &n); - write_source_bank(file_id, source_bank); + write_source_bank(file_id); if (mpi::master || parallel) file_close(file_id); } void -write_source_bank(hid_t group_id, Bank* source_bank) +write_source_bank(hid_t group_id) { hid_t banktype = h5banktype(); @@ -103,7 +101,7 @@ write_source_bank(hid_t group_id, Bank* source_bank) H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE); // Write data to file in parallel - H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank); + H5Dwrite(dset, banktype, memspace, dspace, plist, simulation::source_bank.data()); // Free resources H5Sclose(dspace); @@ -122,7 +120,8 @@ write_source_bank(hid_t group_id, Bank* source_bank) // Save source bank sites since the souce_bank array is overwritten below #ifdef OPENMC_MPI - std::vector temp_source {source_bank, source_bank + simulation::work}; + std::vector temp_source {simulation::source_bank.begin(), + simulation::source_bank.begin() + simulation::work}; #endif for (int i = 0; i < mpi::n_procs; ++i) { @@ -134,7 +133,7 @@ write_source_bank(hid_t group_id, Bank* source_bank) #ifdef OPENMC_MPI // Receive source sites from other processes if (i > 0) - MPI_Recv(source_bank, count[0], mpi::bank, i, i, + MPI_Recv(simulation::source_bank.data(), count[0], mpi::bank, i, i, mpi::intracomm, MPI_STATUS_IGNORE); #endif @@ -144,7 +143,8 @@ write_source_bank(hid_t group_id, Bank* source_bank) H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr); // Write data to hyperslab - H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank); + H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, + simulation::source_bank.data()); H5Sclose(memspace); H5Sclose(dspace); @@ -155,12 +155,12 @@ write_source_bank(hid_t group_id, Bank* source_bank) #ifdef OPENMC_MPI // Restore state of source bank - std::copy(temp_source.begin(), temp_source.end(), source_bank); + std::copy(temp_source.begin(), temp_source.end(), simulation::source_bank.begin()); #endif } else { #ifdef OPENMC_MPI - MPI_Send(source_bank, simulation::work, mpi::bank, 0, mpi::rank, - mpi::intracomm); + MPI_Send(simulation::source_bank.data(), simulation::work, mpi::bank, + 0, mpi::rank, mpi::intracomm); #endif } #endif @@ -169,7 +169,7 @@ write_source_bank(hid_t group_id, Bank* source_bank) } -void read_source_bank(hid_t group_id, Bank* source_bank) +void read_source_bank(hid_t group_id) { hid_t banktype = h5banktype(); @@ -197,10 +197,10 @@ void read_source_bank(hid_t group_id, Bank* source_bank) // Read data in parallel hid_t plist = H5Pcreate(H5P_DATASET_XFER); H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE); - H5Dread(dset, banktype, memspace, dspace, plist, source_bank); + H5Dread(dset, banktype, memspace, dspace, plist, simulation::source_bank.data()); H5Pclose(plist); #else - H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank); + H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, simulation::source_bank.data()); #endif // Close all ids diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index b82837ea1e..1889d0dcba 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -3,6 +3,7 @@ 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 @@ -687,7 +688,7 @@ contains do k = 1, p % n_bank ! get the delayed group - g = fission_bank(n_bank - p % n_bank + k) % delayed_group + g = fission_bank_delayed_group(n_bank - p % n_bank + k) ! Case for tallying delayed emissions if (g /= 0) then @@ -697,8 +698,8 @@ contains reactions(nuclides(p % event_nuclide) % index_fission(1))) ! determine score based on bank site weight and keff. - score = score + keff * fission_bank(n_bank - p % n_bank + k) & - % wgt * rxn % product_decay_rate(1 + g) * flux + score = score + keff * fission_bank_wgt(n_bank - p % n_bank + k) & + * rxn % product_decay_rate(1 + g) * flux end associate ! if the delayed group filter is present, tally to corresponding @@ -1843,7 +1844,7 @@ contains do k = 1, p % n_bank ! get the delayed group - g = fission_bank(n_bank - p % n_bank + k) % delayed_group + g = fission_bank_delayed_group(n_bank - p % n_bank + k) ! Case for tallying delayed emissions if (g /= 0) then @@ -1851,13 +1852,13 @@ contains ! determine score based on bank site weight and keff. if (i_nuclide > 0) then score = score + keff * atom_density * & - fission_bank(n_bank - p % n_bank + k) % wgt * & + fission_bank_wgt(n_bank - p % n_bank + k) * & get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * & get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / & get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux else score = score + keff * & - fission_bank(n_bank - p % n_bank + k) % wgt * & + fission_bank_wgt(n_bank - p % n_bank + k) * & get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux end if @@ -2375,10 +2376,10 @@ contains do k = 1, p % n_bank ! get the delayed group - g = fission_bank(n_bank - p % n_bank + k) % delayed_group + g = fission_bank_delayed_group(n_bank - p % n_bank + k) ! determine score based on bank site weight and keff - score = keff * fission_bank(n_bank - p % n_bank + k) % wgt + score = keff * fission_bank_wgt(n_bank - p % n_bank + k) ! Add derivative information for differential tallies. Note that the ! i_nuclide and atom_density arguments do not matter since this is an @@ -2390,7 +2391,7 @@ contains if (.not. run_CE .and. eo_filt % matches_transport_groups) then ! determine outgoing energy group from fission bank - g_out = int(fission_bank(n_bank - p % n_bank + k) % E) + g_out = int(fission_bank_E(n_bank - p % n_bank + k)) ! modify the value so that g_out = 1 corresponds to the highest ! energy bin @@ -2403,10 +2404,9 @@ contains ! determine outgoing energy from fission bank if (run_CE) then - E_out = fission_bank(n_bank - p % n_bank + k) % E + E_out = fission_bank_E(n_bank - p % n_bank + k) else - E_out = energy_bin_avg(int(fission_bank(n_bank - p % n_bank + k) & - % E)) + E_out = energy_bin_avg(int(fission_bank_E(n_bank - p % n_bank + k))) end if ! If this outgoing energy falls within the energyout filter's range, diff --git a/tests/regression_tests/dagmc/test.py b/tests/regression_tests/dagmc/test.py index 6db936d12b..9766c9732f 100644 --- a/tests/regression_tests/dagmc/test.py +++ b/tests/regression_tests/dagmc/test.py @@ -6,7 +6,7 @@ import pytest from tests.testing_harness import PyAPITestHarness pytestmark = pytest.mark.skipif( - not openmc.capi.dagmc_enabled, + not openmc.capi._dagmc_enabled(), reason="DAGMC CAD geometry is not enabled.") def test_dagmc(): @@ -22,7 +22,7 @@ def test_dagmc(): model.settings.source = source model.settings.dagmc = True - + # tally tally = openmc.Tally() tally.scores = ['total']