diff --git a/CMakeLists.txt b/CMakeLists.txt index c9bb9f3c80..9b9f5a3c95 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -354,7 +354,6 @@ add_library(libopenmc SHARED src/timer_header.F90 src/tracking.F90 src/track_output.F90 - src/urr_header.F90 src/vector_header.F90 src/volume_calc.F90 src/volume_header.F90 @@ -374,6 +373,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_solver.cpp @@ -401,12 +401,15 @@ 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 src/position.cpp src/progress_bar.cpp src/pugixml/pugixml_c.cpp + src/urr.cpp src/random_lcg.cpp src/reaction.cpp src/reaction_product.cpp diff --git a/docs/source/methods/random_numbers.rst b/docs/source/methods/random_numbers.rst index 138a6c75e4..ce118556e4 100644 --- a/docs/source/methods/random_numbers.rst +++ b/docs/source/methods/random_numbers.rst @@ -58,7 +58,7 @@ relationship: \xi_{i+k} = g^k \xi_i + c \frac{g^k - 1}{g - 1} \mod M -Note that :eq:`lcg-skipahead` has the same general form as \eqref{eq:lcg}, so +Note that equation :eq:`lcg-skipahead` has the same general form as equation :eq:`lcg`, so the idea is to determine the new multiplicative and additive constants in :math:`O(\log_2 N)` operations. diff --git a/examples/jupyter/mdgxs-part-i.ipynb b/examples/jupyter/mdgxs-part-i.ipynb index 6d84ca3c21..4ca079b318 100644 --- a/examples/jupyter/mdgxs-part-i.ipynb +++ b/examples/jupyter/mdgxs-part-i.ipynb @@ -325,7 +325,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings_file.source = openmc.source.Source(space=uniform_dist)\n", + "settings_file.source = openmc.Source(space=uniform_dist)\n", "\n", "# Export to \"settings.xml\"\n", "settings_file.export_to_xml()" diff --git a/examples/jupyter/mdgxs-part-ii.ipynb b/examples/jupyter/mdgxs-part-ii.ipynb index 14736f7be5..4c2cc40486 100644 --- a/examples/jupyter/mdgxs-part-ii.ipynb +++ b/examples/jupyter/mdgxs-part-ii.ipynb @@ -320,7 +320,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-10.71, -10.71, -10, 10.71, 10.71, 10.]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings.source = openmc.source.Source(space=uniform_dist)\n", + "settings.source = openmc.Source(space=uniform_dist)\n", "\n", "# Export to \"settings.xml\"\n", "settings.export_to_xml()" diff --git a/examples/jupyter/mg-mode-part-i.ipynb b/examples/jupyter/mg-mode-part-i.ipynb index 7fcb3d079a..8e52787879 100644 --- a/examples/jupyter/mg-mode-part-i.ipynb +++ b/examples/jupyter/mg-mode-part-i.ipynb @@ -148,7 +148,7 @@ "source": [ "To build the actual 2-D model, we will first begin by creating the `materials.xml` file.\n", "\n", - "First we need to define materials that will be used in the problem. In other notebooks, either nuclides or elements were added to materials at the equivalent stage. We can do that in multi-group mode as well. However, multi-group cross-sections are sometimes provided as macroscopic cross-sections; the C5G7 benchmark data are macroscopic. In this case, we can instead use the `Material.add_macroscopic` method to specific a macroscopic object. Unlike for nuclides and elements, we do not need provide information on atom/weight percents as no number densities are needed.\n", + "First we need to define materials that will be used in the problem. In other notebooks, either nuclides or elements were added to materials at the equivalent stage. We can do that in multi-group mode as well. However, multi-group cross-sections are sometimes provided as macroscopic cross-sections; the C5G7 benchmark data are macroscopic. In this case, we can instead use the `Material.add_macroscopic` method to specify a macroscopic object. Unlike for nuclides and elements, we do not need provide information on atom/weight percents as no number densities are needed.\n", "\n", "When assigning macroscopic objects to a material, the density can still be scaled by setting the density to a value that is not 1.0. This would be useful, for example, when slightly perturbing the density of water due to a small change in temperature (while of course ignoring any resultant spectral shift). The density of a macroscopic dataset is set to 1.0 in the `openmc.Material` object by default when a macroscopic dataset is used; so we will show its use the first time and then afterwards it will not be required.\n", "\n", @@ -514,7 +514,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-32.13, -10.71, -1e50, 10.71, 32.13, 1e50]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings_file.source = openmc.source.Source(space=uniform_dist)\n", + "settings_file.source = openmc.Source(space=uniform_dist)\n", "\n", "# Tell OpenMC we want to run in eigenvalue mode\n", "settings_file.run_mode = 'eigenvalue'\n", diff --git a/examples/jupyter/mg-mode-part-ii.ipynb b/examples/jupyter/mg-mode-part-ii.ipynb index 6380c32b08..e2df6fd3db 100644 --- a/examples/jupyter/mg-mode-part-ii.ipynb +++ b/examples/jupyter/mg-mode-part-ii.ipynb @@ -351,7 +351,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-10.71, -10.71, -10, 10.71, 10.71, 10.]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings_file.source = openmc.source.Source(space=uniform_dist)\n", + "settings_file.source = openmc.Source(space=uniform_dist)\n", "\n", "# Export to \"settings.xml\"\n", "settings_file.export_to_xml()" diff --git a/examples/jupyter/mg-mode-part-iii.ipynb b/examples/jupyter/mg-mode-part-iii.ipynb index 2f37070c14..7487b842c6 100644 --- a/examples/jupyter/mg-mode-part-iii.ipynb +++ b/examples/jupyter/mg-mode-part-iii.ipynb @@ -464,7 +464,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [pin_pitch, pin_pitch, 10, length - pin_pitch, length - pin_pitch, 10]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings_file.source = openmc.source.Source(space=uniform_dist)\n", + "settings_file.source = openmc.Source(space=uniform_dist)\n", "\n", "# Export to \"settings.xml\"\n", "settings_file.export_to_xml()" diff --git a/examples/jupyter/mgxs-part-i.ipynb b/examples/jupyter/mgxs-part-i.ipynb index d09aeaa464..90beeaae26 100644 --- a/examples/jupyter/mgxs-part-i.ipynb +++ b/examples/jupyter/mgxs-part-i.ipynb @@ -310,7 +310,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings_file.source = openmc.source.Source(space=uniform_dist)\n", + "settings_file.source = openmc.Source(space=uniform_dist)\n", "\n", "# Export to \"settings.xml\"\n", "settings_file.export_to_xml()" diff --git a/examples/jupyter/mgxs-part-ii.ipynb b/examples/jupyter/mgxs-part-ii.ipynb index 1058b71349..797a357f01 100644 --- a/examples/jupyter/mgxs-part-ii.ipynb +++ b/examples/jupyter/mgxs-part-ii.ipynb @@ -256,7 +256,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings_file.source = openmc.source.Source(space=uniform_dist)\n", + "settings_file.source = openmc.Source(space=uniform_dist)\n", "\n", "# Activate tally precision triggers\n", "settings_file.trigger_active = True\n", diff --git a/examples/jupyter/mgxs-part-iii.ipynb b/examples/jupyter/mgxs-part-iii.ipynb index 4171486a81..d3e90d0946 100644 --- a/examples/jupyter/mgxs-part-iii.ipynb +++ b/examples/jupyter/mgxs-part-iii.ipynb @@ -335,7 +335,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-10.71, -10.71, -10, 10.71, 10.71, 10.]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings_file.source = openmc.source.Source(space=uniform_dist)\n", + "settings_file.source = openmc.Source(space=uniform_dist)\n", "\n", "# Export to \"settings.xml\"\n", "settings_file.export_to_xml()" diff --git a/examples/jupyter/pandas-dataframes.ipynb b/examples/jupyter/pandas-dataframes.ipynb index 140a0d6847..b8b10311c4 100644 --- a/examples/jupyter/pandas-dataframes.ipynb +++ b/examples/jupyter/pandas-dataframes.ipynb @@ -236,7 +236,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-10.71, -10.71, -10, 10.71, 10.71, 10.]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings.source = openmc.source.Source(space=uniform_dist)\n", + "settings.source = openmc.Source(space=uniform_dist)\n", "\n", "# Export to \"settings.xml\"\n", "settings.export_to_xml()" diff --git a/examples/jupyter/post-processing.ipynb b/examples/jupyter/post-processing.ipynb index 4c0f4ef004..fd25e810d7 100644 --- a/examples/jupyter/post-processing.ipynb +++ b/examples/jupyter/post-processing.ipynb @@ -239,7 +239,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings_file.source = openmc.source.Source(space=uniform_dist)\n", + "settings_file.source = openmc.Source(space=uniform_dist)\n", "\n", "# Export to \"settings.xml\"\n", "settings_file.export_to_xml()" diff --git a/examples/jupyter/tally-arithmetic.ipynb b/examples/jupyter/tally-arithmetic.ipynb index 5e904759cf..2bccf73f2a 100644 --- a/examples/jupyter/tally-arithmetic.ipynb +++ b/examples/jupyter/tally-arithmetic.ipynb @@ -240,7 +240,7 @@ "# Create an initial uniform spatial source distribution over fissionable zones\n", "bounds = [-0.63, -0.63, -100., 0.63, 0.63, 100.]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", - "settings_file.source = openmc.source.Source(space=uniform_dist)\n", + "settings_file.source = openmc.Source(space=uniform_dist)\n", "\n", "# Export to \"settings.xml\"\n", "settings_file.export_to_xml()" 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 ca2b0f3ab9..aad72a1386 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/constants.h b/include/openmc/constants.h index c3c9babac8..96fc19b2bd 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -257,12 +257,12 @@ constexpr int LIBRARY_PHOTON {3}; constexpr int LIBRARY_MULTIGROUP {4}; // Probability table parameters -constexpr int URR_CUM_PROB {1}; -constexpr int URR_TOTAL {2}; -constexpr int URR_ELASTIC {3}; -constexpr int URR_FISSION {4}; -constexpr int URR_N_GAMMA {5}; -constexpr int URR_HEATING {6}; +constexpr int URR_CUM_PROB {0}; +constexpr int URR_TOTAL {1}; +constexpr int URR_ELASTIC {2}; +constexpr int URR_FISSION {3}; +constexpr int URR_N_GAMMA {4}; +constexpr int URR_HEATING {5}; // Maximum number of partial fission reactions constexpr int PARTIAL_FISSION_MAX {4}; @@ -424,7 +424,7 @@ constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE // Interpolation rules enum class Interpolation { - histogram, lin_lin, lin_log, log_lin, log_log + histogram = 1, lin_lin = 2, lin_log = 3, log_lin = 4, log_log = 5 }; // Run modes 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/hdf5_interface.h b/include/openmc/hdf5_interface.h index 313b841e20..5b5d2a9991 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -272,6 +272,33 @@ void read_dataset(hid_t obj_id, const char* name, xt::xarray& arr, bool indep } +template +void read_dataset(hid_t obj_id, const char* name, xt::xtensor& arr, bool indep=false) +{ + // Open dataset and read array + hid_t dset = open_dataset(obj_id, name); + + // Get shape of dataset + std::vector hsize_t_shape = object_shape(dset); + close_dataset(dset); + + // cast from hsize_t to size_t + std::vector shape(hsize_t_shape.size()); + for (int i = 0; i < shape.size(); i++) { + shape[i] = static_cast(hsize_t_shape[i]); + } + + // Allocate new xarray to read data into + xt::xarray xarr(shape); + + // Read data from the dataset + read_dataset(obj_id, name, xarr); + + // Copy into xtensor + arr = xarr; + +} + template void read_dataset_as_shape(hid_t obj_id, const char* name, xt::xtensor& arr, bool indep=false) 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..455dab5fc7 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -5,23 +5,24 @@ #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" +#include "urr.h" namespace openmc { //============================================================================== -// Global variables +// Constants //============================================================================== -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; - -} // namespace data +constexpr double CACHE_INVALID {-1.0}; //=============================================================================== //! Cached microscopic cross sections for a particular nuclide at the current @@ -66,20 +67,117 @@ 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 + // 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 +}; + +//=============================================================================== +// Data for a nuclide +//=============================================================================== + +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, int i_nuclide); + + // Methods + double nu(double E, EmissionMode mode, int group=0) const; + void calculate_elastic_xs() const; + + //! Determines the microscopic 0K elastic cross section at a trial relative + //! energy used in resonance scattering + double elastic_xs_0K(double E) const; + + //! \brief Determines cross sections in the unresolved resonance range + //! from probability tables. + void calculate_urr_xs(int i_temp, double E); + + // 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 + int i_nuclide_; //! Index in the nuclides array + + 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_; + + // Unresolved resonance range information + bool urr_present_ {false}; + int urr_inelastic_ {C_NONE}; + std::vector urr_data_; + + std::vector> reactions_; //! Reactions + std::vector index_inelastic_scatter_; + +private: + void create_derived(); +}; + +//============================================================================== +// 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); +extern "C" void nuclide_calculate_urr_xs(bool use_mp, int i_nuclide, + int i_temp, double E); + } // namespace openmc #endif // OPENMC_NUCLIDE_H 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/physics_mg.h b/include/openmc/physics_mg.h index e5648c4a52..8febad19a5 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -10,24 +10,18 @@ namespace openmc { -//TODO: Remove material_xs parameters when they reside on -// the C-side this should happen after materials, physics, input, and tallies -// are brought over - //! \brief samples particle behavior after a collision event. //! \param p Particle to operate on -//! \param material_xs The cross section cache for the current material extern "C" void -collision_mg(Particle* p, const MaterialMacroXS* material_xs); +collision_mg(Particle* p); //! \brief samples a reaction type. //! //! Note that there is special logic when suvival biasing is turned on since //! fission and disappearance are treated implicitly. //! \param p Particle to operate on -//! \param material_xs The cross section cache for the current material void -sample_reaction(Particle* p, const MaterialMacroXS* material_xs); +sample_reaction(Particle* p); //! \brief Samples the scattering event //! \param p Particle to operate on @@ -40,16 +34,14 @@ scatter(Particle* p); //! \param bank_array The particle bank to populate //! \param size_bank Number of particles currently in the bank //! \param bank_array_size Allocated size of the bank -//! \param material_xs The cross section cache for the current material void create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, - int64_t bank_array_size, const MaterialMacroXS* material_xs); + int64_t bank_array_size); //! \brief Handles an absorption event //! \param p Particle to operate on -//! \param material_xs The cross section cache for the current material void -absorption(Particle* p, const MaterialMacroXS* material_xs); +absorption(Particle* p); } // namespace openmc #endif // OPENMC_PHYSICS_MG_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/include/openmc/urr.h b/include/openmc/urr.h new file mode 100644 index 0000000000..6764254425 --- /dev/null +++ b/include/openmc/urr.h @@ -0,0 +1,33 @@ +//! \brief UrrData information for the unresolved resonance treatment + +#ifndef OPENMC_URR_H +#define OPENMC_URR_H + +#include "xtensor/xtensor.hpp" + +#include "openmc/constants.h" +#include "openmc/hdf5_interface.h" + +namespace openmc { + +//============================================================================== +//! UrrData contains probability tables for the unresolved resonance range. +//============================================================================== + +class UrrData{ +public: + Interpolation interp_; //!< interpolation type + int inelastic_flag_; //!< inelastic competition flag + int absorption_flag_; //!< other absorption flag + bool multiply_smooth_; //!< multiply by smooth cross section? + int n_energy_; //!< number of energy points + xt::xtensor energy_; //!< incident energies + xt::xtensor prob_; //!< Actual probability tables + + //! \brief Load the URR data from the provided HDF5 group + explicit UrrData(hid_t group_id); +}; + +} // namespace openmc + +#endif // OPENMC_URR_H 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/openmc/capi/core.py b/openmc/capi/core.py index d5acded941..ef606d6352 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -136,7 +136,7 @@ def find_material(xyz): _dll.openmc_find_cell((c_double*3)(*xyz), index, instance) mats = openmc.capi.Cell(index=index.value).fill - if isinstance(mats, openmc.capi.Material): + if isinstance(mats, (openmc.capi.Material, type(None))): return mats else: return mats[instance.value] diff --git a/openmc/mgxs/__init__.py b/openmc/mgxs/__init__.py index 86e9dce9ef..b3ce9438c3 100644 --- a/openmc/mgxs/__init__.py +++ b/openmc/mgxs/__init__.py @@ -7,10 +7,14 @@ from openmc.mgxs.mdgxs import * GROUP_STRUCTURES = {} """Dictionary of commonly used energy group structures, including "CASMO-X" (where X -is 2, 4, 8, 16, 25, 40 or 70) from the CASMO_ lattice physics code. +is 2, 4, 8, 16, 25, 40 or 70) from the CASMO_ lattice physics code and other commonly +used activation_ energy group structures "VITAMIN-J-175", "TRIPOLI-315", "CCFE-709_" +and "UKAEA-1102_" .. _CASMO: https://www.studsvik.com/SharepointFiles/CASMO-5%20Development%20and%20Applications.pdf - +.. _activation: https://fispact.ukaea.uk/wiki/Keyword:GETXS +.. _CCFE-709: https://fispact.ukaea.uk/wiki/CCFE-709_group_structure +.. _UKAEA-1102: https://fispact.ukaea.uk/wiki/UKAEA-1102_group_structure """ GROUP_STRUCTURES['CASMO-2'] = np.array([ @@ -42,3 +46,445 @@ GROUP_STRUCTURES['CASMO-70'] = np.array([ 3.6726e2, 9.069e2, 1.4251e3, 2.2395e3, 3.5191e3, 5.53e3, 9.118e3, 1.503e4, 2.478e4, 4.085e4, 6.734e4, 1.11e5, 1.83e5, 3.025e5, 5.e5, 8.21e5, 1.353e6, 2.231e6, 3.679e6, 6.0655e6, 2.e7]) +GROUP_STRUCTURES['VITAMIN-J-175'] = np.array([ + 1.0000e-5, 1.0000e-1, 4.1399e-1, 5.3158e-1, 6.8256e-1, + 8.7643e-1, 1.1253, 1.4450, 1.8554, 2.3824, 3.0590, + 3.9279, 5.0435, 6.4759, 8.3153, 1.0677e1, 1.3710e1, + 1.7604e1, 2.2603e1, 2.9023e1, 3.7266e1, 4.7851e1, 6.1442e1, + 7.8893e1, 1.0130e2, 1.3007e2, 1.6702e2, 2.1445e2, 2.7536e2, + 3.5358e2, 4.5400e2, 5.8295e2, 7.4852e2, 9.6112e2, 1.2341e3, + 1.5846e3, 2.0347e3, 2.2487e3, 2.4852e3, 2.6126e3, 2.7465e3, + 3.0354e3, 3.3546e3, 3.7074e3, 4.3074e3, 5.5308e3, 7.1017e3, + 9.1188e3, 1.0595e4, 1.1709e4, 1.5034e4, 1.9304e4, 2.1875e4, + 2.3579e4, 2.4176e4, 2.4788e4, 2.6058e4, 2.7000e4, 2.8501e4, + 3.1828e4, 3.4307e4, 4.0868e4, 4.6309e4, 5.2475e4, 5.6562e4, + 6.7380e4, 7.2024e4, 7.9499e4, 8.2503e4, 8.6517e4, 9.8036e4, + 1.1109e5, 1.1679e5, 1.2277e5, 1.2907e5, 1.3569e5, 1.4264e5, + 1.4996e5, 1.5764e5, 1.6573e5, 1.7422e5, 1.8316e5, 1.9255e5, + 2.0242e5, 2.1280e5, 2.2371e5, 2.3518e5, 2.4724e5, 2.7324e5, + 2.8725e5, 2.9452e5, 2.9721e5, 2.9849e5, 3.0197e5, 3.3373e5, + 3.6883e5, 3.8774e5, 4.0762e5, 4.5049e5, 4.9787e5, 5.2340e5, + 5.5023e5, 5.7844e5, 6.0810e5, 6.3928e5, 6.7206e5, 7.0651e5, + 7.4274e5, 7.8082e5, 8.2085e5, 8.6294e5, 9.0718e5, 9.6167e5, + 1.0026e6, 1.1080e6, 1.1648e6, 1.2246e6, 1.2874e6, 1.3534e6, + 1.4227e6, 1.4957e6, 1.5724e6, 1.6530e6, 1.7377e6, 1.8268e6, + 1.9205e6, 2.0190e6, 2.1225e6, 2.2313e6, 2.3069e6, 2.3457e6, + 2.3653e6, 2.3851e6, 2.4660e6, 2.5924e6, 2.7253e6, 2.8650e6, + 3.0119e6, 3.1664e6, 3.3287e6, 3.6788e6, 4.0657e6, 4.4933e6, + 4.7237e6, 4.9658e6, 5.2205e6, 5.4881e6, 5.7695e6, 6.0653e6, + 6.3763e6, 6.5924e6, 6.7032e6, 7.0469e6, 7.4082e6, 7.7880e6, + 8.1873e6, 8.6071e6, 9.0484e6, 9.5123e6, 1.0000e7, 1.0513e7, + 1.1052e7, 1.1618e7, 1.2214e7, 1.2523e7, 1.2840e7, 1.3499e7, + 1.3840e7, 1.4191e7, 1.4550e7, 1.4918e7, 1.5683e7, 1.6487e7, + 1.6905e7, 1.7332e7, 1.9640e7]) +GROUP_STRUCTURES['TRIPOLI-315,'] = np.array([ + 1.0e-5, 1.1e-4, 3.000e-3, 5.500e-3, 1.000e-2, 1.500e-2, 2.000e-2, 3.000e-2, + 3.200e-2, 3.238e-2, 4.300e-2, 5.900e-2, 7.700e-2, 9.500e-2, 1.000e-1, + 1.150e-1, 1.340e-1, 1.600e-1, 1.890e-1, 2.200e-1, 2.480e-1, 2.825e-1, + 3.145e-1, 3.520e-1, 3.910e-1, 4.140e-1, 4.330e-1, 4.850e-1, 5.316e-1, + 5.400e-1, 6.250e-1, 6.826e-1, 7.050e-1, 7.900e-1, 8.600e-1, 8.764e-1, + 9.300e-1, 9.860e-1, 1.010, 1.035, 1.070, 1.080, 1.090, + 1.110, 1.125, 1.170, 1.235, 1.305, 1.370, 1.440, + 1.445, 1.510, 1.590, 1.670, 1.755, 1.840, 1.855, + 1.930, 2.020, 2.130, 2.360, 2.372, 2.768, 3.059, + 3.381, 3.928, 4.129, 4.470, 4.670, 5.043, 5.623, + 6.160, 6.476, 7.079, 7.524, 7.943, 8.315, 8.913, + 9.190, 1.000e1, 1.068e1, 1.122e1, 1.259e1, 1.371e1, 1.523e1, + 1.674e1, 1.760e1, 1.903e1, 2.045e1, 2.260e1, 2.498e1, 2.792e1, + 2.920e1, 3.051e1, 3.389e1, 3.727e1, 3.981e1, 4.552e1, 4.785e1, + 5.012e1, 5.559e1, 6.144e1, 6.310e1, 6.790e1, 7.079e1, 7.889e1, + 8.528e1, 9.166e1, 1.013e2, 1.122e2, 1.301e2, 1.367e2, 1.585e2, + 1.670e2, 1.778e2, 2.040e2, 2.145e2, 2.430e2, 2.754e2, 3.043e2, + 3.536e2, 3.981e2, 4.540e2, 5.145e2, 5.830e2, 6.310e2, 6.773e2, + 7.079e2, 7.485e2, 8.482e2, 9.611e2, 1.010e3, 1.117e3, 1.234e3, + 1.364e3, 1.507e3, 1.585e3, 1.796e3, 2.035e3, 2.113e3, 2.249e3, + 2.371e3, 2.485e3, 2.613e3, 2.661e3, 2.747e3, 2.818e3, 3.035e3, + 3.162e3, 3.355e3, 3.548e3, 3.707e3, 3.981e3, 4.307e3, 4.643e3, + 5.004e3, 5.531e3, 6.267e3, 7.102e3, 7.466e3, 8.251e3, 9.119e3, + 1.008e4, 1.114e4, 1.171e4, 1.273e4, 1.383e4, 1.503e4, 1.585e4, + 1.662e4, 1.778e4, 1.931e4, 1.995e4, 2.054e4, 2.113e4, 2.187e4, + 2.239e4, 2.304e4, 2.358e4, 2.418e4, 2.441e4, 2.479e4, 2.512e4, + 2.585e4, 2.606e4, 2.661e4, 2.700e4, 2.738e4, 2.818e4, 2.850e4, + 2.901e4, 2.985e4, 3.073e4, 3.162e4, 3.183e4, 3.431e4, 3.698e4, + 4.087e4, 4.359e4, 4.631e4, 4.939e4, 5.248e4, 5.517e4, 5.656e4, + 6.173e4, 6.738e4, 7.200e4, 7.499e4, 7.950e4, 8.230e4, 8.250e4, + 8.652e4, 9.804e4, 1.111e5, 1.168e5, 1.228e5, 1.291e5, 1.357e5, + 1.426e5, 1.500e5, 1.576e5, 1.657e5, 1.742e5, 1.832e5, 1.925e5, + 2.024e5, 2.128e5, 2.237e5, 2.352e5, 2.472e5, 2.732e5, 2.873e5, + 2.945e5, 2.972e5, 2.985e5, 3.020e5, 3.337e5, 3.688e5, 3.877e5, + 4.076e5, 4.505e5, 5.234e5, 5.502e5, 5.784e5, 6.081e5, 6.393e5, + 6.721e5, 7.065e5, 7.427e5, 7.808e5, 8.209e5, 8.629e5, 9.072e5, + 9.616e5, 1.003e6, 1.108e6, 1.165e6, 1.225e6, 1.287e6, 1.353e6, + 1.423e6, 1.496e6, 1.572e6, 1.653e6, 1.738e6, 1.827e6, 1.921e6, + 2.019e6, 2.122e6, 2.231e6, 2.307e6, 2.346e6, 2.365e6, 2.385e6, + 2.466e6, 2.592e6, 2.725e6, 2.865e6, 3.012e6, 3.166e6, 3.329e6, + 3.679e6, 4.066e6, 4.493e6, 4.724e6, 4.966e6, 5.220e6, 5.488e6, + 5.769e6, 6.065e6, 6.376e6, 6.592e6, 6.703e6, 7.047e6, 7.408e6, + 7.788e6, 8.187e6, 8.607e6, 9.048e6, 9.512e6, 1.000e7, 1.051e7, + 1.105e7, 1.162e7, 1.221e7, 1.284e7, 1.350e7, 1.384e7, 1.419e7, + 1.455e7, 1.492e7, 1.568e7, 1.649e7, 1.691e7, 1.733e7, 1.964e7]) +GROUP_STRUCTURES['CCFE-709'] = np.array([ + 1.e-5, 1.0471e-5, 1.0965e-5, 1.1482e-5, 1.2023e-5, + 1.2589e-5, 1.3183e-5, 1.3804e-5, 1.4454e-5, 1.5136e-5, + 1.5849e-5, 1.6596e-5, 1.7378e-5, 1.8197e-5, 1.9055e-5, + 1.9953e-5, 2.0893e-5, 2.1878e-5, 2.2909e-5, 2.3988e-5, + 2.5119e-5, 2.6303e-5, 2.7542e-5, 2.8840e-5, 3.0200e-5, + 3.1623e-5, 3.3113e-5, 3.4674e-5, 3.6308e-5, 3.8019e-5, + 3.9811e-5, 4.1687e-5, 4.3652e-5, 4.5709e-5, 4.7863e-5, + 5.0119e-5, 5.2481e-5, 5.4954e-5, 5.7544e-5, 6.0256e-5, + 6.3096e-5, 6.6069e-5, 6.9183e-5, 7.2444e-5, 7.5858e-5, + 7.9433e-5, 8.3176e-5, 8.7096e-5, 9.1201e-5, 9.5499e-5, + 1.0000e-4, 1.0471e-4, 1.0965e-4, 1.1482e-4, 1.2023e-4, + 1.2589e-4, 1.3183e-4, 1.3804e-4, 1.4454e-4, 1.5136e-4, + 1.5849e-4, 1.6596e-4, 1.7378e-4, 1.8197e-4, 1.9055e-4, + 1.9953e-4, 2.0893e-4, 2.1878e-4, 2.2909e-4, 2.3988e-4, + 2.5119e-4, 2.6303e-4, 2.7542e-4, 2.8840e-4, 3.0200e-4, + 3.1623e-4, 3.3113e-4, 3.4674e-4, 3.6308e-4, 3.8019e-4, + 3.9811e-4, 4.1687e-4, 4.3652e-4, 4.5709e-4, 4.7863e-4, + 5.0119e-4, 5.2481e-4, 5.4954e-4, 5.7544e-4, 6.0256e-4, + 6.3096e-4, 6.6069e-4, 6.9183e-4, 7.2444e-4, 7.5858e-4, + 7.9433e-4, 8.3176e-4, 8.7096e-4, 9.1201e-4, 9.5499e-4, + 1.0000e-3, 1.0471e-3, 1.0965e-3, 1.1482e-3, 1.2023e-3, + 1.2589e-3, 1.3183e-3, 1.3804e-3, 1.4454e-3, 1.5136e-3, + 1.5849e-3, 1.6596e-3, 1.7378e-3, 1.8197e-3, 1.9055e-3, + 1.9953e-3, 2.0893e-3, 2.1878e-3, 2.2909e-3, 2.3988e-3, + 2.5119e-3, 2.6303e-3, 2.7542e-3, 2.8840e-3, 3.0200e-3, + 3.1623e-3, 3.3113e-3, 3.4674e-3, 3.6308e-3, 3.8019e-3, + 3.9811e-3, 4.1687e-3, 4.3652e-3, 4.5709e-3, 4.7863e-3, + 5.0119e-3, 5.2481e-3, 5.4954e-3, 5.7544e-3, 6.0256e-3, + 6.3096e-3, 6.6069e-3, 6.9183e-3, 7.2444e-3, 7.5858e-3, + 7.9433e-3, 8.3176e-3, 8.7096e-3, 9.1201e-3, 9.5499e-3, + 1.0000e-2, 1.0471e-2, 1.0965e-2, 1.1482e-2, 1.2023e-2, + 1.2589e-2, 1.3183e-2, 1.3804e-2, 1.4454e-2, 1.5136e-2, + 1.5849e-2, 1.6596e-2, 1.7378e-2, 1.8197e-2, 1.9055e-2, + 1.9953e-2, 2.0893e-2, 2.1878e-2, 2.2909e-2, 2.3988e-2, + 2.5119e-2, 2.6303e-2, 2.7542e-2, 2.8840e-2, 3.0200e-2, + 3.1623e-2, 3.3113e-2, 3.4674e-2, 3.6308e-2, 3.8019e-2, + 3.9811e-2, 4.1687e-2, 4.3652e-2, 4.5709e-2, 4.7863e-2, + 5.0119e-2, 5.2481e-2, 5.4954e-2, 5.7544e-2, 6.0256e-2, + 6.3096e-2, 6.6069e-2, 6.9183e-2, 7.2444e-2, 7.5858e-2, + 7.9433e-2, 8.3176e-2, 8.7096e-2, 9.1201e-2, 9.5499e-2, + 1.0000e-1, 1.0471e-1, 1.0965e-1, 1.1482e-1, 1.2023e-1, + 1.2589e-1, 1.3183e-1, 1.3804e-1, 1.4454e-1, 1.5136e-1, + 1.5849e-1, 1.6596e-1, 1.7378e-1, 1.8197e-1, 1.9055e-1, + 1.9953e-1, 2.0893e-1, 2.1878e-1, 2.2909e-1, 2.3988e-1, + 2.5119e-1, 2.6303e-1, 2.7542e-1, 2.8840e-1, 3.0200e-1, + 3.1623e-1, 3.3113e-1, 3.4674e-1, 3.6308e-1, 3.8019e-1, + 3.9811e-1, 4.1687e-1, 4.3652e-1, 4.5709e-1, 4.7863e-1, + 5.0119e-1, 5.2481e-1, 5.4954e-1, 5.7544e-1, 6.0256e-1, + 6.3096e-1, 6.6069e-1, 6.9183e-1, 7.2444e-1, 7.5858e-1, + 7.9433e-1, 8.3176e-1, 8.7096e-1, 9.1201e-1, 9.5499e-1, + 1.0000e0, 1.0471e0, 1.0965e0, 1.1482e0, 1.2023e0, + 1.2589e0, 1.3183e0, 1.3804e0, 1.4454e0, 1.5136e0, + 1.5849e0, 1.6596e0, 1.7378e0, 1.8197e0, 1.9055e0, + 1.9953e0, 2.0893e0, 2.1878e0, 2.2909e0, 2.3988e0, + 2.5119e0, 2.6303e0, 2.7542e0, 2.8840e0, 3.0200e0, + 3.1623e0, 3.3113e0, 3.4674e0, 3.6308e0, 3.8019e0, + 3.9811e0, 4.1687e0, 4.3652e0, 4.5709e0, 4.7863e0, + 5.0119e0, 5.2481e0, 5.4954e0, 5.7544e0, 6.0256e0, + 6.3096e0, 6.6069e0, 6.9183e0, 7.2444e0, 7.5858e0, + 7.9433e0, 8.3176e0, 8.7096e0, 9.1201e0, 9.5499e0, + 1.0000e1, 1.0471e1, 1.0965e1, 1.1482e1, 1.2023e1, + 1.2589e1, 1.3183e1, 1.3804e1, 1.4454e1, 1.5136e1, + 1.5849e1, 1.6596e1, 1.7378e1, 1.8197e1, 1.9055e1, + 1.9953e1, 2.0893e1, 2.1878e1, 2.2909e1, 2.3988e1, + 2.5119e1, 2.6303e1, 2.7542e1, 2.8840e1, 3.0200e1, + 3.1623e1, 3.3113e1, 3.4674e1, 3.6308e1, 3.8019e1, + 3.9811e1, 4.1687e1, 4.3652e1, 4.5709e1, 4.7863e1, + 5.0119e1, 5.2481e1, 5.4954e1, 5.7544e1, 6.0256e1, + 6.3096e1, 6.6069e1, 6.9183e1, 7.2444e1, 7.5858e1, + 7.9433e1, 8.3176e1, 8.7096e1, 9.1201e1, 9.5499e1, + 1.0000e2, 1.0471e2, 1.0965e2, 1.1482e2, 1.2023e2, + 1.2589e2, 1.3183e2, 1.3804e2, 1.4454e2, 1.5136e2, + 1.5849e2, 1.6596e2, 1.7378e2, 1.8197e2, 1.9055e2, + 1.9953e2, 2.0893e2, 2.1878e2, 2.2909e2, 2.3988e2, + 2.5119e2, 2.6303e2, 2.7542e2, 2.8840e2, 3.0200e2, + 3.1623e2, 3.3113e2, 3.4674e2, 3.6308e2, 3.8019e2, + 3.9811e2, 4.1687e2, 4.3652e2, 4.5709e2, 4.7863e2, + 5.0119e2, 5.2481e2, 5.4954e2, 5.7544e2, 6.0256e2, + 6.3096e2, 6.6069e2, 6.9183e2, 7.2444e2, 7.5858e2, + 7.9433e2, 8.3176e2, 8.7096e2, 9.1201e2, 9.5499e2, + 1.0000e3, 1.0471e3, 1.0965e3, 1.1482e3, 1.2023e3, + 1.2589e3, 1.3183e3, 1.3804e3, 1.4454e3, 1.5136e3, + 1.5849e3, 1.6596e3, 1.7378e3, 1.8197e3, 1.9055e3, + 1.9953e3, 2.0893e3, 2.1878e3, 2.2909e3, 2.3988e3, + 2.5119e3, 2.6303e3, 2.7542e3, 2.8840e3, 3.0200e3, + 3.1623e3, 3.3113e3, 3.4674e3, 3.6308e3, 3.8019e3, + 3.9811e3, 4.1687e3, 4.3652e3, 4.5709e3, 4.7863e3, + 5.0119e3, 5.2481e3, 5.4954e3, 5.7544e3, 6.0256e3, + 6.3096e3, 6.6069e3, 6.9183e3, 7.2444e3, 7.5858e3, + 7.9433e3, 8.3176e3, 8.7096e3, 9.1201e3, 9.5499e3, + 1.0000e4, 1.0471e4, 1.0965e4, 1.1482e4, 1.2023e4, + 1.2589e4, 1.3183e4, 1.3804e4, 1.4454e4, 1.5136e4, + 1.5849e4, 1.6596e4, 1.7378e4, 1.8197e4, 1.9055e4, + 1.9953e4, 2.0893e4, 2.1878e4, 2.2909e4, 2.3988e4, + 2.5119e4, 2.6303e4, 2.7542e4, 2.8840e4, 3.0200e4, + 3.1623e4, 3.3113e4, 3.4674e4, 3.6308e4, 3.8019e4, + 3.9811e4, 4.1687e4, 4.3652e4, 4.5709e4, 4.7863e4, + 5.0119e4, 5.2481e4, 5.4954e4, 5.7544e4, 6.0256e4, + 6.3096e4, 6.6069e4, 6.9183e4, 7.2444e4, 7.5858e4, + 7.9433e4, 8.3176e4, 8.7096e4, 9.1201e4, 9.5499e4, + 1.0000e5, 1.0471e5, 1.0965e5, 1.1482e5, 1.2023e5, + 1.2589e5, 1.3183e5, 1.3804e5, 1.4454e5, 1.5136e5, + 1.5849e5, 1.6596e5, 1.7378e5, 1.8197e5, 1.9055e5, + 1.9953e5, 2.0893e5, 2.1878e5, 2.2909e5, 2.3988e5, + 2.5119e5, 2.6303e5, 2.7542e5, 2.8840e5, 3.0200e5, + 3.1623e5, 3.3113e5, 3.4674e5, 3.6308e5, 3.8019e5, + 3.9811e5, 4.1687e5, 4.3652e5, 4.5709e5, 4.7863e5, + 5.0119e5, 5.2481e5, 5.4954e5, 5.7544e5, 6.0256e5, + 6.3096e5, 6.6069e5, 6.9183e5, 7.2444e5, 7.5858e5, + 7.9433e5, 8.3176e5, 8.7096e5, 9.1201e5, 9.5499e5, + 1.0000e6, 1.0471e6, 1.0965e6, 1.1482e6, 1.2023e6, + 1.2589e6, 1.3183e6, 1.3804e6, 1.4454e6, 1.5136e6, + 1.5849e6, 1.6596e6, 1.7378e6, 1.8197e6, 1.9055e6, + 1.9953e6, 2.0893e6, 2.1878e6, 2.2909e6, 2.3988e6, + 2.5119e6, 2.6303e6, 2.7542e6, 2.8840e6, 3.0200e6, + 3.1623e6, 3.3113e6, 3.4674e6, 3.6308e6, 3.8019e6, + 3.9811e6, 4.1687e6, 4.3652e6, 4.5709e6, 4.7863e6, + 5.0119e6, 5.2481e6, 5.4954e6, 5.7544e6, 6.0256e6, + 6.3096e6, 6.6069e6, 6.9183e6, 7.2444e6, 7.5858e6, + 7.9433e6, 8.3176e6, 8.7096e6, 9.1201e6, 9.5499e6, + 1.0000e7, 1.0200e7, 1.0400e7, 1.0600e7, 1.0800e7, + 1.1000e7, 1.1200e7, 1.1400e7, 1.1600e7, 1.1800e7, + 1.2000e7, 1.2200e7, 1.2400e7, 1.2600e7, 1.2800e7, + 1.3000e7, 1.3200e7, 1.3400e7, 1.3600e7, 1.3800e7, + 1.4000e7, 1.4200e7, 1.4400e7, 1.4600e7, 1.4800e7, + 1.5000e7, 1.5200e7, 1.5400e7, 1.5600e7, 1.5800e7, + 1.6000e7, 1.6200e7, 1.6400e7, 1.6600e7, 1.6800e7, + 1.7000e7, 1.7200e7, 1.7400e7, 1.7600e7, 1.7800e7, + 1.8000e7, 1.8200e7, 1.8400e7, 1.8600e7, 1.8800e7, + 1.9000e7, 1.9200e7, 1.9400e7, 1.9600e7, 1.9800e7, + 2.0000e7, 2.1000e7, 2.2000e7, 2.3000e7, 2.4000e7, + 2.5000e7, 2.6000e7, 2.7000e7, 2.8000e7, 2.9000e7, + 3.0000e7, 3.2000e7, 3.4000e7, 3.6000e7, 3.8000e7, + 4.0000e7, 4.2000e7, 4.4000e7, 4.6000e7, 4.8000e7, + 5.0000e7, 5.2000e7, 5.4000e7, 5.6000e7, 5.8000e7, + 6.0000e7, 6.5000e7, 7.0000e7, 7.5000e7, 8.0000e7, + 9.0000e7, 1.0000e8, 1.1000e8, 1.2000e8, 1.3000e8, + 1.4000e8, 1.5000e8, 1.6000e8, 1.8000e8, 2.0000e8, + 2.4000e8, 2.8000e8, 3.2000e8, 3.6000e8, 4.0000e8, + 4.4000e8, 4.8000e8, 5.2000e8, 5.6000e8, 6.0000e8, + 6.4000e8, 6.8000e8, 7.2000e8, 7.6000e8, 8.0000e8, + 8.4000e8, 8.8000e8, 9.2000e8, 9.6000e8, 1.0000e9,]) +GROUP_STRUCTURES['UKAEA-1102'] = np.array([ + 1.0000e-5, 1.0471e-5, 1.0965e-5, 1.1482e-5, 1.2023e-5, + 1.2589e-5, 1.3183e-5, 1.3804e-5, 1.4454e-5, 1.5136e-5, + 1.5849e-5, 1.6596e-5, 1.7378e-5, 1.8197e-5, 1.9055e-5, + 1.9953e-5, 2.0893e-5, 2.1878e-5, 2.2909e-5, 2.3988e-5, + 2.5119e-5, 2.6303e-5, 2.7542e-5, 2.8840e-5, 3.0200e-5, + 3.1623e-5, 3.3113e-5, 3.4674e-5, 3.6308e-5, 3.8019e-5, + 3.9811e-5, 4.1687e-5, 4.3652e-5, 4.5709e-5, 4.7863e-5, + 5.0119e-5, 5.2481e-5, 5.4954e-5, 5.7544e-5, 6.0256e-5, + 6.3096e-5, 6.6069e-5, 6.9183e-5, 7.2444e-5, 7.5858e-5, + 7.9433e-5, 8.3176e-5, 8.7096e-5, 9.1201e-5, 9.5499e-5, + 1.0000e-4, 1.0471e-4, 1.0965e-4, 1.1482e-4, 1.2023e-4, + 1.2589e-4, 1.3183e-4, 1.3804e-4, 1.4454e-4, 1.5136e-4, + 1.5849e-4, 1.6596e-4, 1.7378e-4, 1.8197e-4, 1.9055e-4, + 1.9953e-4, 2.0893e-4, 2.1878e-4, 2.2909e-4, 2.3988e-4, + 2.5119e-4, 2.6303e-4, 2.7542e-4, 2.8840e-4, 3.0200e-4, + 3.1623e-4, 3.3113e-4, 3.4674e-4, 3.6308e-4, 3.8019e-4, + 3.9811e-4, 4.1687e-4, 4.3652e-4, 4.5709e-4, 4.7863e-4, + 5.0119e-4, 5.2481e-4, 5.4954e-4, 5.7544e-4, 6.0256e-4, + 6.3096e-4, 6.6069e-4, 6.9183e-4, 7.2444e-4, 7.5858e-4, + 7.9433e-4, 8.3176e-4, 8.7096e-4, 9.1201e-4, 9.5499e-4, + 1.0000e-3, 1.0471e-3, 1.0965e-3, 1.1482e-3, 1.2023e-3, + 1.2589e-3, 1.3183e-3, 1.3804e-3, 1.4454e-3, 1.5136e-3, + 1.5849e-3, 1.6596e-3, 1.7378e-3, 1.8197e-3, 1.9055e-3, + 1.9953e-3, 2.0893e-3, 2.1878e-3, 2.2909e-3, 2.3988e-3, + 2.5119e-3, 2.6303e-3, 2.7542e-3, 2.8840e-3, 3.0200e-3, + 3.1623e-3, 3.3113e-3, 3.4674e-3, 3.6308e-3, 3.8019e-3, + 3.9811e-3, 4.1687e-3, 4.3652e-3, 4.5709e-3, 4.7863e-3, + 5.0119e-3, 5.2481e-3, 5.4954e-3, 5.7544e-3, 6.0256e-3, + 6.3096e-3, 6.6069e-3, 6.9183e-3, 7.2444e-3, 7.5858e-3, + 7.9433e-3, 8.3176e-3, 8.7096e-3, 9.1201e-3, 9.5499e-3, + 1.0000e-2, 1.0471e-2, 1.0965e-2, 1.1482e-2, 1.2023e-2, + 1.2589e-2, 1.3183e-2, 1.3804e-2, 1.4454e-2, 1.5136e-2, + 1.5849e-2, 1.6596e-2, 1.7378e-2, 1.8197e-2, 1.9055e-2, + 1.9953e-2, 2.0893e-2, 2.1878e-2, 2.2909e-2, 2.3988e-2, + 2.5119e-2, 2.6303e-2, 2.7542e-2, 2.8840e-2, 3.0200e-2, + 3.1623e-2, 3.3113e-2, 3.4674e-2, 3.6308e-2, 3.8019e-2, + 3.9811e-2, 4.1687e-2, 4.3652e-2, 4.5709e-2, 4.7863e-2, + 5.0119e-2, 5.2481e-2, 5.4954e-2, 5.7544e-2, 6.0256e-2, + 6.3096e-2, 6.6069e-2, 6.9183e-2, 7.2444e-2, 7.5858e-2, + 7.9433e-2, 8.3176e-2, 8.7096e-2, 9.1201e-2, 9.5499e-2, + 1.0000e-1, 1.0471e-1, 1.0965e-1, 1.1482e-1, 1.2023e-1, + 1.2589e-1, 1.3183e-1, 1.3804e-1, 1.4454e-1, 1.5136e-1, + 1.5849e-1, 1.6596e-1, 1.7378e-1, 1.8197e-1, 1.9055e-1, + 1.9953e-1, 2.0893e-1, 2.1878e-1, 2.2909e-1, 2.3988e-1, + 2.5119e-1, 2.6303e-1, 2.7542e-1, 2.8840e-1, 3.0200e-1, + 3.1623e-1, 3.3113e-1, 3.4674e-1, 3.6308e-1, 3.8019e-1, + 3.9811e-1, 4.1687e-1, 4.3652e-1, 4.5709e-1, 4.7863e-1, + 5.0119e-1, 5.2481e-1, 5.5000e-1, 5.7500e-1, 6.0000e-1, + 6.2500e-1, 6.5000e-1, 6.7500e-1, 7.0000e-1, 7.2500e-1, + 7.5000e-1, 7.7500e-1, 8.0000e-1, 8.2500e-1, 8.5000e-1, + 8.7500e-1, 9.0000e-1, 9.2500e-1, 9.5000e-1, 9.7500e-1, + 1.0000, 1.0250, 1.0500, 1.0750, 1.1000, + 1.1250, 1.1500, 1.1750, 1.2000, 1.2250, + 1.2500, 1.2750, 1.3000, 1.3250, 1.3500, + 1.3750, 1.4000, 1.4250, 1.4500, 1.4750, + 1.5000, 1.5250, 1.5500, 1.5750, 1.6000, + 1.6250, 1.6500, 1.6750, 1.7000, 1.7250, + 1.7500, 1.7750, 1.8000, 1.8250, 1.8500, + 1.8750, 1.9000, 1.9250, 1.9500, 1.9750, + 2.0000, 2.0250, 2.0500, 2.0750, 2.1000, + 2.1250, 2.1500, 2.1750, 2.2000, 2.2250, + 2.2500, 2.2750, 2.3000, 2.3250, 2.3500, + 2.3750, 2.4000, 2.4250, 2.4500, 2.4750, + 2.5000, 2.5250, 2.5500, 2.5750, 2.6000, + 2.6250, 2.6500, 2.6750, 2.7000, 2.7250, + 2.7500, 2.7750, 2.8000, 2.8250, 2.8500, + 2.8750, 2.9000, 2.9250, 2.9500, 2.9750, + 3.0000, 3.0250, 3.0500, 3.0750, 3.1000, + 3.1250, 3.1500, 3.1750, 3.2000, 3.2250, + 3.2500, 3.2750, 3.3000, 3.3250, 3.3500, + 3.3750, 3.4000, 3.4250, 3.4500, 3.4750, + 3.5000, 3.5250, 3.5500, 3.5750, 3.6000, + 3.6250, 3.6500, 3.6750, 3.7000, 3.7250, + 3.7500, 3.7750, 3.8000, 3.8250, 3.8500, + 3.8750, 3.9000, 3.9250, 3.9500, 3.9750, + 4.0000, 4.0250, 4.0500, 4.0750, 4.1000, + 4.1250, 4.1500, 4.1750, 4.2000, 4.2250, + 4.2500, 4.2750, 4.3000, 4.3250, 4.3500, + 4.3750, 4.4000, 4.4250, 4.4500, 4.4750, + 4.5000, 4.5250, 4.5500, 4.5750, 4.6000, + 4.6250, 4.6500, 4.6750, 4.7000, 4.7250, + 4.7500, 4.7750, 4.8000, 4.8250, 4.8500, + 4.8750, 4.9000, 4.9250, 4.9500, 4.9750, + 5.0000, 5.0250, 5.0500, 5.0750, 5.1000, + 5.1250, 5.1500, 5.1750, 5.2000, 5.2250, + 5.2500, 5.2750, 5.3000, 5.3250, 5.3500, + 5.3750, 5.4000, 5.4250, 5.4500, 5.4750, + 5.5000, 5.5250, 5.5500, 5.5750, 5.6000, + 5.6250, 5.6500, 5.6750, 5.7000, 5.7250, + 5.7500, 5.7750, 5.8000, 5.8250, 5.8500, + 5.8750, 5.9000, 5.9250, 5.9500, 5.9750, + 6.0000, 6.0250, 6.0500, 6.0750, 6.1000, + 6.1250, 6.1500, 6.1750, 6.2000, 6.2250, + 6.2500, 6.2750, 6.3000, 6.3250, 6.3500, + 6.3750, 6.4000, 6.4250, 6.4500, 6.4750, + 6.5000, 6.5250, 6.5500, 6.5750, 6.6000, + 6.6250, 6.6500, 6.6750, 6.7000, 6.7250, + 6.7500, 6.7750, 6.8000, 6.8250, 6.8500, + 6.8750, 6.9000, 6.9250, 6.9500, 6.9750, + 7.0000, 7.0250, 7.0500, 7.0750, 7.1000, + 7.1250, 7.1500, 7.1750, 7.2000, 7.2250, + 7.2500, 7.2750, 7.3000, 7.3250, 7.3500, + 7.3750, 7.4000, 7.4250, 7.4500, 7.4750, + 7.5000, 7.5250, 7.5500, 7.5750, 7.6000, + 7.6250, 7.6500, 7.6750, 7.7000, 7.7250, + 7.7500, 7.7750, 7.8000, 7.8250, 7.8500, + 7.8750, 7.9000, 7.9250, 7.9500, 7.9750, + 8.0000, 8.0250, 8.0500, 8.0750, 8.1000, + 8.1250, 8.1500, 8.1750, 8.2000, 8.2250, + 8.2500, 8.2750, 8.3000, 8.3250, 8.3500, + 8.3750, 8.4000, 8.4250, 8.4500, 8.4750, + 8.5000, 8.5250, 8.5500, 8.5750, 8.6000, + 8.6250, 8.6500, 8.6750, 8.7000, 8.7250, + 8.7500, 8.7750, 8.8000, 8.8250, 8.8500, + 8.8750, 8.9000, 8.9250, 8.9500, 8.9750, + 9.0000, 9.0250, 9.0500, 9.0750, 9.1000, + 9.1250, 9.1500, 9.1750, 9.2000, 9.2250, + 9.2500, 9.2750, 9.3000, 9.3250, 9.3500, + 9.3750, 9.4000, 9.4250, 9.4500, 9.4750, + 9.5000, 9.5250, 9.5500, 9.5750, 9.6000, + 9.6250, 9.6500, 9.6750, 9.7000, 9.7250, + 9.7500, 9.7750, 9.8000, 9.8250, 9.8500, + 9.8750, 9.9000, 9.9250, 9.9500, 9.9750, + 1.0000e1, 1.0471e1, 1.0965e1, 1.1482e1, 1.2023e1, + 1.2589e1, 1.3183e1, 1.3804e1, 1.4454e1, 1.5136e1, + 1.5849e1, 1.6596e1, 1.7378e1, 1.8197e1, 1.9055e1, + 1.9953e1, 2.0893e1, 2.1878e1, 2.2909e1, 2.3988e1, + 2.5119e1, 2.6303e1, 2.7542e1, 2.8840e1, 3.0200e1, + 3.1623e1, 3.3113e1, 3.4674e1, 3.6308e1, 3.8019e1, + 3.9811e1, 4.1687e1, 4.3652e1, 4.5709e1, 4.7863e1, + 5.0119e1, 5.2481e1, 5.4954e1, 5.7544e1, 6.0256e1, + 6.3096e1, 6.6069e1, 6.9183e1, 7.2444e1, 7.5858e1, + 7.9433e1, 8.3176e1, 8.7096e1, 9.1201e1, 9.5499e1, + 1.0000e2, 1.0471e2, 1.0965e2, 1.1482e2, 1.2023e2, + 1.2589e2, 1.3183e2, 1.3804e2, 1.4454e2, 1.5136e2, + 1.5849e2, 1.6596e2, 1.7378e2, 1.8197e2, 1.9055e2, + 1.9953e2, 2.0893e2, 2.1878e2, 2.2909e2, 2.3988e2, + 2.5119e2, 2.6303e2, 2.7542e2, 2.8840e2, 3.0200e2, + 3.1623e2, 3.3113e2, 3.4674e2, 3.6308e2, 3.8019e2, + 3.9811e2, 4.1687e2, 4.3652e2, 4.5709e2, 4.7863e2, + 5.0119e2, 5.2481e2, 5.4954e2, 5.7544e2, 6.0256e2, + 6.3096e2, 6.6069e2, 6.9183e2, 7.2444e2, 7.5858e2, + 7.9433e2, 8.3176e2, 8.7096e2, 9.1201e2, 9.5499e2, + 1.0000e3, 1.0471e3, 1.0965e3, 1.1482e3, 1.2023e3, + 1.2589e3, 1.3183e3, 1.3804e3, 1.4454e3, 1.5136e3, + 1.5849e3, 1.6596e3, 1.7378e3, 1.8197e3, 1.9055e3, + 1.9953e3, 2.0893e3, 2.1878e3, 2.2909e3, 2.3988e3, + 2.5119e3, 2.6303e3, 2.7542e3, 2.8840e3, 3.0200e3, + 3.1623e3, 3.3113e3, 3.4674e3, 3.6308e3, 3.8019e3, + 3.9811e3, 4.1687e3, 4.3652e3, 4.5709e3, 4.7863e3, + 5.0119e3, 5.2481e3, 5.4954e3, 5.7544e3, 6.0256e3, + 6.3096e3, 6.6069e3, 6.9183e3, 7.2444e3, 7.5858e3, + 7.9433e3, 8.3176e3, 8.7096e3, 9.1201e3, 9.5499e3, + 1.0000e4, 1.0471e4, 1.0965e4, 1.1482e4, 1.2023e4, + 1.2589e4, 1.3183e4, 1.3804e4, 1.4454e4, 1.5136e4, + 1.5849e4, 1.6596e4, 1.7378e4, 1.8197e4, 1.9055e4, + 1.9953e4, 2.0893e4, 2.1878e4, 2.2909e4, 2.3988e4, + 2.5119e4, 2.6303e4, 2.7542e4, 2.8840e4, 3.0200e4, + 3.1623e4, 3.3113e4, 3.4674e4, 3.6308e4, 3.8019e4, + 3.9811e4, 4.1687e4, 4.3652e4, 4.5709e4, 4.7863e4, + 5.0119e4, 5.2481e4, 5.4954e4, 5.7544e4, 6.0256e4, + 6.3096e4, 6.6069e4, 6.9183e4, 7.2444e4, 7.5858e4, + 7.9433e4, 8.3176e4, 8.7096e4, 9.1201e4, 9.5499e4, + 1.0000e5, 1.0471e5, 1.0965e5, 1.1482e5, 1.2023e5, + 1.2589e5, 1.3183e5, 1.3804e5, 1.4454e5, 1.5136e5, + 1.5849e5, 1.6596e5, 1.7378e5, 1.8197e5, 1.9055e5, + 1.9953e5, 2.0893e5, 2.1878e5, 2.2909e5, 2.3988e5, + 2.5119e5, 2.6303e5, 2.7542e5, 2.8840e5, 3.0200e5, + 3.1623e5, 3.3113e5, 3.4674e5, 3.6308e5, 3.8019e5, + 3.9811e5, 4.1687e5, 4.3652e5, 4.5709e5, 4.7863e5, + 5.0119e5, 5.2481e5, 5.4954e5, 5.7544e5, 6.0256e5, + 6.3096e5, 6.6069e5, 6.9183e5, 7.2444e5, 7.5858e5, + 7.9433e5, 8.3176e5, 8.7096e5, 9.1201e5, 9.5499e5, + 1.0000e6, 1.0471e6, 1.0965e6, 1.1482e6, 1.2023e6, + 1.2589e6, 1.3183e6, 1.3804e6, 1.4454e6, 1.5136e6, + 1.5849e6, 1.6596e6, 1.7378e6, 1.8197e6, 1.9055e6, + 1.9953e6, 2.0893e6, 2.1878e6, 2.2909e6, 2.3988e6, + 2.5119e6, 2.6303e6, 2.7542e6, 2.8840e6, 3.0200e6, + 3.1623e6, 3.3113e6, 3.4674e6, 3.6308e6, 3.8019e6, + 3.9811e6, 4.1687e6, 4.3652e6, 4.5709e6, 4.7863e6, + 5.0000e6, 5.2000e6, 5.4000e6, 5.6000e6, 5.8000e6, + 6.0000e6, 6.2000e6, 6.4000e6, 6.6000e6, 6.8000e6, + 7.0000e6, 7.2000e6, 7.4000e6, 7.6000e6, 7.8000e6, + 8.0000e6, 8.2000e6, 8.4000e6, 8.6000e6, 8.8000e6, + 9.0000e6, 9.2000e6, 9.4000e6, 9.6000e6, 9.8000e6, + 1.0000e7, 1.0200e7, 1.0400e7, 1.0600e7, 1.0800e7, + 1.1000e7, 1.1200e7, 1.1400e7, 1.1600e7, 1.1800e7, + 1.2000e7, 1.2200e7, 1.2400e7, 1.2600e7, 1.2800e7, + 1.3000e7, 1.3200e7, 1.3400e7, 1.3600e7, 1.3800e7, + 1.4000e7, 1.4200e7, 1.4400e7, 1.4600e7, 1.4800e7, + 1.5000e7, 1.5200e7, 1.5400e7, 1.5600e7, 1.5800e7, + 1.6000e7, 1.6200e7, 1.6400e7, 1.6600e7, 1.6800e7, + 1.7000e7, 1.7200e7, 1.7400e7, 1.7600e7, 1.7800e7, + 1.8000e7, 1.8200e7, 1.8400e7, 1.8600e7, 1.8800e7, + 1.9000e7, 1.9200e7, 1.9400e7, 1.9600e7, 1.9800e7, + 2.0000e7, 2.0200e7, 2.0400e7, 2.0600e7, 2.0800e7, + 2.1000e7, 2.1200e7, 2.1400e7, 2.1600e7, 2.1800e7, + 2.2000e7, 2.2200e7, 2.2400e7, 2.2600e7, 2.2800e7, + 2.3000e7, 2.3200e7, 2.3400e7, 2.3600e7, 2.3800e7, + 2.4000e7, 2.4200e7, 2.4400e7, 2.4600e7, 2.4800e7, + 2.5000e7, 2.5200e7, 2.5400e7, 2.5600e7, 2.5800e7, + 2.6000e7, 2.6200e7, 2.6400e7, 2.6600e7, 2.6800e7, + 2.7000e7, 2.7200e7, 2.7400e7, 2.7600e7, 2.7800e7, + 2.8000e7, 2.8200e7, 2.8400e7, 2.8600e7, 2.8800e7, + 2.9000e7, 2.9200e7, 2.9400e7, 2.9600e7, 2.9800e7, + 3.0000e7, 3.0200e7, 3.1623e7, 3.3113e7, 3.4674e7, + 3.6308e7, 3.8019e7, 3.9811e7, 4.1687e7, 4.3652e7, + 4.5709e7, 4.7863e7, 5.0119e7, 5.2481e7, 5.4954e7, + 5.7544e7, 6.0256e7, 6.3096e7, 6.6069e7, 6.9183e7, + 7.2444e7, 7.5858e7, 7.9433e7, 8.3176e7, 8.7096e7, + 9.1201e7, 9.5499e7, 1.0000e8, 1.0471e8, 1.0965e8, + 1.1482e8, 1.2023e8, 1.2589e8, 1.3183e8, 1.3804e8, + 1.4454e8, 1.5136e8, 1.5849e8, 1.6596e8, 1.7378e8, + 1.8197e8, 1.9055e8, 1.9953e8, 2.0893e8, 2.1878e8, + 2.2909e8, 2.3988e8, 2.5119e8, 2.6303e8, 2.7542e8, + 2.8840e8, 3.0200e8, 3.1623e8, 3.3113e8, 3.4674e8, + 3.6308e8, 3.8019e8, 3.9811e8, 4.1687e8, 4.3652e8, + 4.5709e8, 4.7863e8, 5.0119e8, 5.2481e8, 5.4954e8, + 5.7544e8, 6.0256e8, 6.3096e8, 6.6069e8, 6.9183e8, + 7.2444e8, 7.5858e8, 7.9433e8, 8.3176e8, 8.7096e8, + 9.1201e8, 9.5499e8, 1.e9]) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index b077d04d1d..11b2b24ed7 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -3900,7 +3900,7 @@ class ScatterMatrixXS(MatrixMGXS): if self.formulation == 'simple': if self.scatter_format == 'legendre': - # If using P0 correction subtract P2 scatter from the diag. + # If using P0 correction subtract P1 scatter from the diag. if self.correction == 'P0' and self.legendre_order == 0: scatter_p0 = self.tallies[self.rxn_type].get_slice( filters=[openmc.LegendreFilter], diff --git a/scripts/casl_chain.py b/scripts/casl_chain.py new file mode 100755 index 0000000000..86389e3ced --- /dev/null +++ b/scripts/casl_chain.py @@ -0,0 +1,277 @@ +# This dictionary contains the 255-nuclides, simplified burnup chain used in +# CASL-ORIGEN, which can be found in Appendix A of Kang Seog Kim, "Specification +# for the VERA Depletion Benchmark Suite", CASL-U-2015-1014-000, Rev. 0, +# ORNL/TM-2016/53, 2016. +# +# Note 32 of the 255 nuclides appeare twice as they are both activation +# nuclides (category 1) and fission product nuclides (category 3). + +CASL_CHAIN = { + # Nuclide: (Stable, CAT, IFPY, Special yield treatment) + # Stable: True if nuclide has no decay reactions + # CAT: Category of nuclides + # 1-Activation nuclides + # 2-Heavy metal nuclides + # 3-Fission product nuclides + # IFPY: Indicator of fission product yield + # 0-Non FPY + # 1-Direct FPY (-1 indicates (stable+metastable) direct FPY) + # 2-Cumulative FPY + # 3-Special treatment with weight fractions + # Special yield: (nuclide_i/weight_i/IFPY_i) + 'B10': (True, 1, 0, None), + 'B11': (True, 1, 0, None), + 'O16': (True, 1, 0, None), + 'Ag107': (True, 1, 0, None), + 'Ag109': (True, 1, 0, None), # redundant as FP + 'Ag110': (False, 1, 0, None), # redundant as FP + 'Cd110': (True, 1, 0, None), # redundant as FP + 'Cd111': (True, 1, 0, None), # redundant as FP + 'Cd112': (True, 1, 0, None), + 'Cd113': (True, 1, 0, None), # redundant as FP + 'Cd114': (True, 1, 0, None), + 'Cd115': (False, 1, 0, None), + 'In113': (True, 1, 0, None), + 'In115': (True, 1, 0, None), # redundant as FP + 'Sm152': (True, 1, 0, None), # redundant as FP + 'Sm153': (False, 1, 0, None), # redundant as FP + 'Eu151': (True, 1, 0, None), # redundant as FP + 'Eu152': (False, 1, 0, None), + 'Eu152_m1': (False, 1, 0, None), + 'Eu153': (True, 1, 0, None), # redundant as FP + 'Eu154': (False, 1, 0, None), # redundant as FP + 'Eu155': (False, 1, 0, None), # redundant as FP + 'Eu156': (False, 1, 0, None), # redundant as FP + 'Eu157': (False, 1, 0, None), # redundant as FP + 'Gd152': (True, 1, 0, None), + 'Gd154': (True, 1, 0, None), # redundant as FP + 'Gd155': (True, 1, 0, None), # redundant as FP + 'Gd156': (True, 1, 0, None), # redundant as FP + 'Gd157': (True, 1, 0, None), # redundant as FP + 'Gd158': (True, 1, 0, None), # redundant as FP + 'Gd159': (False, 1, 0, None), # redundant as FP + 'Gd160': (True, 1, 0, None), # redundant as FP + 'Gd161': (False, 1, 0, None), # redundant as FP + 'Tb159': (True, 1, 0, None), # redundant as FP + 'Tb160': (False, 1, 0, None), # redundant as FP + 'Tb161': (False, 1, 0, None), # redundant as FP + 'Dy160': (True, 1, 0, None), # redundant as FP + 'Dy161': (True, 1, 0, None), # redundant as FP + 'Dy162': (True, 1, 0, None), # redundant as FP + 'Dy163': (True, 1, 0, None), # redundant as FP + 'Dy164': (True, 1, 0, None), # redundant as FP + 'Dy165': (False, 1, 0, None), # redundant as FP + 'Ho165': (True, 1, 0, None), # redundant as FP + 'Er162': (True, 1, 0, None), + 'Er164': (True, 1, 0, None), + 'Er166': (True, 1, 0, None), + 'Er167': (True, 1, 0, None), + 'Er168': (True, 1, 0, None), + 'Er169': (False, 1, 0, None), + 'Er170': (True, 1, 0, None), + 'Er171': (False, 1, 0, None), + 'Tm169': (True, 1, 0, None), + 'Tm170': (False, 1, 0, None), + 'Tm171': (False, 1, 0, None), + 'Hf174': (True, 1, 0, None), + 'Hf176': (True, 1, 0, None), + 'Hf177': (True, 1, 0, None), + 'Hf178': (True, 1, 0, None), + 'Hf179': (True, 1, 0, None), + 'Hf180': (True, 1, 0, None), + 'Hf181': (False, 1, 0, None), + 'Ta181': (True, 1, 0, None), + 'Ta182': (False, 1, 0, None), + 'Th230': (False, 2, 0, None), + 'Th231': (False, 2, 0, None), + 'Th232': (False, 2, 0, None), + 'Th233': (False, 2, 0, None), + 'Th234': (False, 2, 0, None), + 'Pa231': (False, 2, 0, None), + 'Pa232': (False, 2, 0, None), + 'Pa233': (False, 2, 0, None), + 'Pa234': (False, 2, 0, None), + 'U232': (False, 2, 0, None), + 'U233': (False, 2, 0, None), + 'U234': (False, 2, 0, None), + 'U235': (False, 2, 0, None), + 'U236': (False, 2, 0, None), + 'U237': (False, 2, 0, None), + 'U238': (False, 2, 0, None), + 'U239': (False, 2, 0, None), + 'Np236': (False, 2, 0, None), + 'Np237': (False, 2, 0, None), + 'Np238': (False, 2, 0, None), + 'Np239': (False, 2, 0, None), + 'Np240': (False, 2, 0, None), + 'Np240_m1': (False, 2, 0, None), + 'Pu236': (False, 2, 0, None), + 'Pu237': (False, 2, 0, None), + 'Pu238': (False, 2, 0, None), + 'Pu239': (False, 2, 0, None), + 'Pu240': (False, 2, 0, None), + 'Pu241': (False, 2, 0, None), + 'Pu242': (False, 2, 0, None), + 'Pu243': (False, 2, 0, None), + 'Am241': (False, 2, 0, None), + 'Am242': (False, 2, 0, None), + 'Am242_m1': (False, 2, 0, None), + 'Am243': (False, 2, 0, None), + 'Am244': (False, 2, 0, None), + 'Am244_m1': (False, 2, 0, None), + 'Cm242': (False, 2, 0, None), + 'Cm243': (False, 2, 0, None), + 'Cm244': (False, 2, 0, None), + 'Cm245': (False, 2, 0, None), + 'Cm246': (False, 2, 0, None), + 'Br81': (True, 3, 2, None), + 'Br82': (False, 3, 2, None), + 'Kr82': (True, 3, 3, [('Br82_m1', 0.024, 1), ('Kr82', 1.000, 1)]), + 'Kr83': (True, 3, 2, None), + 'Kr84': (True, 3, 2, None), + 'Kr85': (False, 3, 2, None), + 'Kr86': (True, 3, 2, None), + 'Sr89': (False, 3, 2, None), + 'Sr90': (False, 3, 2, None), + 'Y89': (True, 3, 1, None), + 'Y90': (False, 3, 1, None), + 'Y91': (False, 3, 2, None), + 'Zr91': (True, 3, 1, None), + 'Zr93': (False, 3, 2, None), + 'Zr95': (False, 3, 2, None), + 'Zr96': (True, 3, 2, None), + 'Nb95': (False, 3, 3, [('Nb95',1.000, 1), ('Nb95_m1', 0.944, 1)]), + 'Mo95': (True, 3, 3, [('Nb95_m1',0.056, 1), ('Mo95', 1.000, 1)]), + 'Mo96': (True, 3, 3, [('Nb96',1.000, 1), ('Mo96', 1.000, 1)]), + 'Mo97': (True, 3, 2, None), + 'Mo98': (True, 3, 2, None), + 'Mo99': (False, 3, 2, None), + 'Mo100': (True, 3, 2, None), + 'Tc99': (False, 3, 1, None), + 'Tc99_m1': (False, 3, 1, None), + 'Tc100': (False, 3, 1, None), + 'Ru100': (True, 3, 1, None), + 'Ru101': (True, 3, 2, None), + 'Ru102': (True, 3, 2, None), + 'Ru103': (False, 3, 2, None), + 'Ru104': (True, 3, 2, None), + 'Ru105': (False, 3, 2, None), + 'Ru106': (False, 3, 2, None), + 'Rh102': (False, 3, 1, None), + 'Rh102_m1': (False, 3, 1, None), + 'Rh103': (True, 3, 1, None), + 'Rh103_m1': (False, 3, 1, None), + 'Rh104': (False, 3, 1, None), + 'Rh105': (False, 3, 1, None), + 'Rh105_m1': (False, 3, 1, None), + 'Rh106': (False, 3, 1, None), + 'Rh106_m1': (False, 3, 1, None), + 'Pd104': (True, 3, 1, None), + 'Pd105': (True, 3, 1, None), + 'Pd106': (True, 3, 1, None), + 'Pd107': (False, 3, 2, None), + 'Pd108': (True, 3, 2, None), + 'Pd109': (False, 3, 2, None), + 'Ag109': (True, 3, 1, None), + 'Ag109_m1': (False, 3, 1, None), + 'Ag110': (False, 3, 2, None), + 'Ag110_m1': (False, 3, 2, None), + 'Ag111': (False, 3, 2, None), + 'Cd110': (True, 3, 1, None), + 'Cd111': (True, 3, 3, [('Ag110', -1.000, 2), ('Cd110', 1.000, 2), ('Cd111', 1.000, 1)]), + 'Cd113': (True, 3, 2, None), + 'In115': (True, 3, 2, None), + 'Sb121': (True, 3, 2, None), + 'Sb123': (False, 3, 2, None), + 'Sb125': (False, 3, 2, None), + 'Sb127': (False, 3, 2, None), + 'Te127': (False, 3, -1, None), + 'Te127_m1': (False, 3, -1, None), + 'Te129_m1': (False, 3, 2, None), + 'Te132': (False, 3, 2, None), + 'I127': (True, 3, 1, None), + 'I128': (False, 3, 3, [('I128', 0.931, 2)]), + 'I129': (False, 3, 3, [('I129', 1.000, 2), ('I129', -1.000, 2)]), + 'I130': (False, 3, 2, None), + 'I131': (False, 3, 2, None), + 'I132': (False, 3, 1, None), + 'I135': (False, 3, 2, None), + 'Xe128': (True, 3, 1, None), + 'Xe130': (True, 3, 1, None), + 'Xe131': (True, 3, 1, None), + 'Xe132': (True, 3, 1, None), + 'Xe133': (False, 3, 2, None), + 'Xe134': (True, 3, 2, None), + 'Xe135': (False, 3, 1, None), + 'Xe135_m1': (False, 3, 1, None), + 'Xe136': (True, 3, 2, None), + 'Xe137': (False, 3, 2, None), + 'Cs133': (True, 3, 1, None), + 'Cs134': (False, 3, 1, None), + 'Cs135': (False, 3, 1, None), + 'Cs136': (False, 3, 1, None), + 'Cs137': (False, 3, 1, None), + 'Ba134': (True, 3, 1, None), + 'Ba137': (True, 3, 1, None), + 'Ba140': (False, 3, 2, None), + 'La139': (True, 3, 2, None), + 'La140': (False, 3, 1, None), + 'Ce140': (True, 3, 1, None), + 'Ce141': (False, 3, 2, None), + 'Ce142': (True, 3, 2, None), + 'Ce143': (False, 3, 2, None), + 'Ce144': (False, 3, 2, None), + 'Pr141': (True, 3, 1, None), + 'Pr142': (False, 3, 1, None), + 'Pr143': (False, 3, 1, None), + 'Pr144': (False, 3, 1, None), + 'Nd142': (True, 3, 1, None), + 'Nd143': (True, 3, 1, None), + 'Nd144': (False, 3, 1, None), + 'Nd145': (True, 3, 2, None), + 'Nd146': (True, 3, 2, None), + 'Nd147': (False, 3, 2, None), + 'Nd148': (True, 3, 2, None), + 'Nd149': (False, 3, 2, None), + 'Nd150': (True, 3, 2, None), + 'Nd151': (False, 3, 2, None), + 'Pm147': (False, 3, 1, None), + 'Pm148': (False, 3, -1, None), + 'Pm148_m1': (False, 3, -1, None), + 'Pm149': (False, 3, 1, None), + 'Pm150': (False, 3, 1, None), + 'Pm151': (False, 3, 1, None), + 'Sm147': (False, 3, 1, None), + 'Sm148': (False, 3, 1, None), + 'Sm149': (False, 3, 1, None), + 'Sm150': (True, 3, 1, None), + 'Sm151': (False, 3, 1, None), + 'Sm152': (True, 3, 2, None), + 'Sm153': (False, 3, 2, None), + 'Sm154': (True, 3, 2, None), + 'Sm155': (False, 3, 2, None), + 'Eu151': (True, 3, 1, None), + 'Eu153': (True, 3, 1, None), + 'Eu154': (False, 3, 1, None), + 'Eu155': (False, 3, 1, None), + 'Eu156': (False, 3, 2, None), + 'Eu157': (False, 3, 2, None), + 'Gd154': (True, 3, 1, None), + 'Gd155': (True, 3, 1, None), + 'Gd156': (True, 3, 1, None), + 'Gd157': (True, 3, 1, None), + 'Gd158': (True, 3, 2, None), + 'Gd159': (False, 3, 2, None), + 'Gd160': (True, 3, 2, None), + 'Gd161': (False, 3, 2, None), + 'Tb159': (True, 3, 1, None), + 'Tb160': (False, 3, 1, None), + 'Tb161': (False, 3, 1, None), + 'Dy160': (True, 3, 1, None), + 'Dy161': (True, 3, 1, None), + 'Dy162': (True, 3, 2, None), + 'Dy163': (True, 3, 2, None), + 'Dy164': (True, 3, 2, None), + 'Dy165': (False, 3, 2, None), + 'Ho165': (True, 3, 3, [('Dy165_m1', 0.022, 2), ('Ho165', 1.000, 1)]) +} diff --git a/scripts/openmc-make-depletion-chain-casl b/scripts/openmc-make-depletion-chain-casl new file mode 100755 index 0000000000..7fdeed072a --- /dev/null +++ b/scripts/openmc-make-depletion-chain-casl @@ -0,0 +1,240 @@ +#!/usr/bin/env python3 + +import glob +import os +from zipfile import ZipFile +from collections import OrderedDict, defaultdict +from io import StringIO +from itertools import chain + +try: + import lxml.etree as ET + _have_lxml = True +except ImportError: + import xml.etree.ElementTree as ET + _have_lxml = False + +import openmc.data +import openmc.deplete +from openmc._xml import clean_indentation +from openmc.deplete.chain import _REACTIONS +from openmc.deplete.nuclide import Nuclide, DecayTuple, ReactionTuple +from openmc._utils import download + +from casl_chain import CASL_CHAIN + +URLS = [ + 'http://www.nndc.bnl.gov/endf/b7.1/zips/ENDF-B-VII.1-neutrons.zip', + 'http://www.nndc.bnl.gov/endf/b7.1/zips/ENDF-B-VII.1-decay.zip', + 'http://www.nndc.bnl.gov/endf/b7.1/zips/ENDF-B-VII.1-nfy.zip' +] + +def main(): + if os.path.isdir('./decay') and os.path.isdir('./nfy') and os.path.isdir('./neutrons'): + endf_dir = '.' + elif 'OPENMC_ENDF_DATA' in os.environ: + endf_dir = os.environ['OPENMC_ENDF_DATA'] + else: + for url in URLS: + basename = download(url) + with ZipFile(basename, 'r') as zf: + print('Extracting {}...'.format(basename)) + zf.extractall() + endf_dir = '.' + + decay_files = glob.glob(os.path.join(endf_dir, 'decay', '*.endf')) + fpy_files = glob.glob(os.path.join(endf_dir, 'nfy', '*.endf')) + neutron_files = glob.glob(os.path.join(endf_dir, 'neutrons', '*.endf')) + + # Create a Chain + chain = openmc.deplete.Chain() + + print('Reading ENDF nuclear data from "{}"...'.format(os.path.abspath(endf_dir))) + + # Create dictionary mapping target to filename + print('Processing neutron sub-library files...') + reactions = {} + for f in neutron_files: + evaluation = openmc.data.endf.Evaluation(f) + name = evaluation.gnd_name + if name in CASL_CHAIN: + reactions[name] = {} + for mf, mt, nc, mod in evaluation.reaction_list: + if mf == 3: + file_obj = StringIO(evaluation.section[3, mt]) + openmc.data.endf.get_head_record(file_obj) + q_value = openmc.data.endf.get_cont_record(file_obj)[1] + reactions[name][mt] = q_value + + # Determine what decay and FPY nuclides are available + print('Processing decay sub-library files...') + decay_data = {} + for f in decay_files: + data = openmc.data.Decay(f) + name = data.nuclide['name'] + if name in CASL_CHAIN: + decay_data[name] = data + + for name in CASL_CHAIN: + if name not in decay_data: + print('WARNING: {} has no decay data!'.format(name)) + + print('Processing fission product yield sub-library files...') + fpy_data = {} + for f in fpy_files: + data = openmc.data.FissionProductYields(f) + name = data.nuclide['name'] + if name in CASL_CHAIN: + fpy_data[name] = data + + print('Creating depletion_chain...') + missing_daughter = [] + missing_rx_product = [] + missing_fpy = [] + + for idx, parent in enumerate(sorted(decay_data, key=openmc.data.zam)): + data = decay_data[parent] + + nuclide = Nuclide() + nuclide.name = parent + + chain.nuclides.append(nuclide) + chain.nuclide_dict[parent] = idx + + if not CASL_CHAIN[parent][0] and \ + not data.nuclide['stable'] and data.half_life.nominal_value != 0.0: + nuclide.half_life = data.half_life.nominal_value + nuclide.decay_energy = sum(E.nominal_value for E in + data.average_energies.values()) + sum_br = 0.0 + for i, mode in enumerate(data.modes): + type_ = ','.join(mode.modes) + if mode.daughter in decay_data: + target = mode.daughter + else: + print('missing {} {} {}'.format(parent, ','.join(mode.modes), mode.daughter)) + continue + + # Write branching ratio, taking care to ensure sum is unity + br = mode.branching_ratio.nominal_value + sum_br += br + if i == len(data.modes) - 1 and sum_br != 1.0: + br = 1.0 - sum(m.branching_ratio.nominal_value + for m in data.modes[:-1]) + + # Append decay mode + nuclide.decay_modes.append(DecayTuple(type_, target, br)) + + if parent in reactions: + reactions_available = set(reactions[parent].keys()) + for name, mts, changes in _REACTIONS: + if mts & reactions_available: + delta_A, delta_Z = changes + A = data.nuclide['mass_number'] + delta_A + Z = data.nuclide['atomic_number'] + delta_Z + daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A) + + if name not in chain.reactions: + chain.reactions.append(name) + + if daughter not in decay_data: + missing_rx_product.append((parent, name, daughter)) + daughter = 'Nothing' + + # Store Q value + for mt in sorted(mts): + if mt in reactions[parent]: + q_value = reactions[parent][mt] + break + else: + q_value = 0.0 + + nuclide.reactions.append(ReactionTuple( + name, daughter, q_value, 1.0)) + + if any(mt in reactions_available for mt in [18, 19, 20, 21, 38]): + if parent in fpy_data: + q_value = reactions[parent][18] + nuclide.reactions.append( + ReactionTuple('fission', 0, q_value, 1.0)) + + if 'fission' not in chain.reactions: + chain.reactions.append('fission') + else: + missing_fpy.append(parent) + + if parent in fpy_data: + fpy = fpy_data[parent] + + if fpy.energies is not None: + nuclide.yield_energies = fpy.energies + else: + nuclide.yield_energies = [0.0] + + for E, table_yd, table_yc in zip(nuclide.yield_energies, fpy.independent, fpy.cumulative): + yields = defaultdict(float) + for product in table_yd: + if product in decay_data: + # identifier + ifpy = CASL_CHAIN[product][2] + # 1 for independent + if ifpy == 1: + if product not in table_yd: + print('No independent fission yields found for {} in {}'.format(product, parent)) + else: + yields[product] += table_yd[product].nominal_value + # 2 for cumulative + elif ifpy == 2: + if product not in table_yc: + print('No cumulative fission yields found for {} in {}'.format(product, parent)) + else: + yields[product] += table_yc[product].nominal_value + # -1 for stable + unstable + elif ifpy == -1: + if product not in table_yd: + print('No independent fission yields found for {} in {}'.format(product, parent)) + else: + yields[product] += table_yc[product].nominal_value + product_meta = '{}_m1'.format(product) + if product_meta in table_yd: + yields[product] += table_yc[product_meta].nominal_value + # 3 for special treatment with weight fractions + elif ifpy == 3: + for tuple_i in CASL_CHAIN[product][3]: + name_i, weight_i, ifpy_i = tuple_i + if name_i not in table_yd: + print('No fission yields found for {} in {}'.format(name_i, parent)) + else: + if ifpy_i == 1: + yields[product] += weight_i * table_yd[name_i].nominal_value + elif ifpy_i == 2: + yields[product] += weight_i * table_yc[name_i].nominal_value + + nuclide.yield_data[E] = [] + for k in sorted(yields, key=openmc.data.zam): + nuclide.yield_data[E].append((k, yields[k])) + + # Display warnings + if missing_daughter: + print('The following decay modes have daughters with no decay data:') + for mode in missing_daughter: + print(' {}'.format(mode)) + print('') + + if missing_rx_product: + print('The following reaction products have no decay data:') + for vals in missing_rx_product: + print('{} {} -> {}'.format(*vals)) + print('') + + if missing_fpy: + print('The following fissionable nuclides have no fission product yields:') + for parent in missing_fpy: + print(' ' + parent) + print('') + + chain.export_to_xml('chain_casl.xml') + + +if __name__ == '__main__': + main() diff --git a/src/api.F90 b/src/api.F90 index 30310e5fce..f45ffb235d 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -81,6 +81,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/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 32fd9f2a14..67c418509e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2109,9 +2109,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 @@ -2293,7 +2290,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..f7ef87b56c 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1,5 +1,18 @@ #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/random_lcg.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 +20,518 @@ 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, int i_nuclide) +{ + // 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_); + i_nuclide_ = i_nuclide; + + // 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); + + // Read unresolved resonance probability tables if present + if (object_exists(group, "urr")) { + urr_present_ = true; + urr_data_.reserve(temps_to_read.size()); + + for (int i = 0; i < temps_to_read.size(); i++) { + // Get temperature as a string + std::string temp_str {std::to_string(temps_to_read[i]) + "K"}; + + // Read probability tables for i-th temperature + hid_t urr_group = open_group(group, ("urr/" + temp_str).c_str()); + urr_data_.emplace_back(urr_group); + close_group(urr_group); + + // Check for negative values + if (xt::any(urr_data_[i].prob_ < 0.) && mpi::master) { + warning("Negative value(s) found on probability table for nuclide " + + name_ + " at " + temp_str); + } + } + + // If the inelastic competition flag indicates that the inelastic cross + // section should be determined from a normal reaction cross section, we + // need to get the index of the reaction. + if (temps_to_read.size() > 0) { + if (urr_data_[0].inelastic_flag_ > 0) { + for (int i = 0; i < reactions_.size(); i++) { + if (reactions_[i]->mt_ == urr_data_[0].inelastic_flag_) { + urr_inelastic_ = i; + } + } + + // Abort if no corresponding inelastic reaction was found + if (urr_inelastic_ == C_NONE) { + fatal_error("Could no find inelastic reaction specified on " + "unresolved resonance probability table."); + } + } + } + } + + // 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() 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]; +} + +void Nuclide::calculate_urr_xs(int i_temp, double E) +{ + auto& micro = simulation::micro_xs[i_nuclide_]; + micro.use_ptable = true; + + // Create a shorthand for the URR data + const auto& urr = urr_data_[i_temp]; + + // Determine the energy table + int i_energy = 0; + while(E >= urr.energy_(i_energy + 1)) {++i_energy;}; + + // Sample the probability table using the cumulative distribution + + // Random nmbers for the xs calculation are sampled from a separate stream. + // This guarantees the randomness and, at the same time, makes sure we + // reuse random numbers for the same nuclide at different temperatures, + // therefore preserving correlation of temperature in probability tables. + prn_set_stream(STREAM_URR_PTABLE); + //TODO: to maintain the same random number stream as the Fortran code this + //replaces, the seed is set with i_nuclide_ + 1 instead of i_nuclide_ + double r = future_prn(static_cast(i_nuclide_ + 1)); + prn_set_stream(STREAM_TRACKING); + + int i_low = 0; + while (urr.prob_(i_energy, URR_CUM_PROB, i_low) <= r) {++i_low;}; + + int i_up = 0; + while (urr.prob_(i_energy + 1, URR_CUM_PROB, i_up) <= r) {++i_up;}; + + // Determine elastic, fission, and capture cross sections from the + // probability table + double elastic = 0.; + double fission = 0.; + double capture = 0.; + double f; + if (urr.interp_ == Interpolation::lin_lin) { + // Determine the interpolation factor on the table + f = (E - urr.energy_(i_energy)) / + (urr.energy_(i_energy + 1) - urr.energy_(i_energy)); + + elastic = (1. - f) * urr.prob_(i_energy, URR_ELASTIC, i_low) + + f * urr.prob_(i_energy + 1, URR_ELASTIC, i_up); + fission = (1. - f) * urr.prob_(i_energy, URR_FISSION, i_low) + + f * urr.prob_(i_energy + 1, URR_FISSION, i_up); + capture = (1. - f) * urr.prob_(i_energy, URR_N_GAMMA, i_low) + + f * urr.prob_(i_energy + 1, URR_N_GAMMA, i_up); + } else if (urr.interp_ == Interpolation::log_log) { + // Determine interpolation factor on the table + f = std::log(E / urr.energy_(i_energy)) / + std::log(urr.energy_(i_energy + 1) / urr.energy_(i_energy)); + + // Calculate the elastic cross section/factor + if ((urr.prob_(i_energy, URR_ELASTIC, i_low) > 0.) && + (urr.prob_(i_energy + 1, URR_ELASTIC, i_up) > 0.)) { + elastic = + std::exp((1. - f) * + std::log(urr.prob_(i_energy, URR_ELASTIC, i_low)) + + f * std::log(urr.prob_(i_energy + 1, URR_ELASTIC, i_up))); + } else { + elastic = 0.; + } + + // Calculate the fission cross section/factor + if ((urr.prob_(i_energy, URR_FISSION, i_low) > 0.) && + (urr.prob_(i_energy + 1, URR_FISSION, i_up) > 0.)) { + fission = + std::exp((1. - f) * + std::log(urr.prob_(i_energy, URR_FISSION, i_low)) + + f * std::log(urr.prob_(i_energy + 1, URR_FISSION, i_up))); + } else { + fission = 0.; + } + + // Calculate the capture cross section/factor + if ((urr.prob_(i_energy, URR_N_GAMMA, i_low) > 0.) && + (urr.prob_(i_energy + 1, URR_N_GAMMA, i_up) > 0.)) { + capture = + std::exp((1. - f) * + std::log(urr.prob_(i_energy, URR_N_GAMMA, i_low)) + + f * std::log(urr.prob_(i_energy + 1, URR_N_GAMMA, i_up))); + } else { + capture = 0.; + } + } + + // Determine the treatment of inelastic scattering + double inelastic = 0.; + if (urr.inelastic_flag_ != C_NONE) { + // get interpolation factor + f = micro.interp_factor; + + // Determine inelastic scattering cross section + Reaction* rx = reactions_[urr_inelastic_].get(); + int xs_index = micro.index_grid - rx->xs_[i_temp].threshold; + if (xs_index >= 0) { + inelastic = (1. - f) * rx->xs_[i_temp].value[xs_index] + + f * rx->xs_[i_temp].value[xs_index + 1]; + } + } + + // Multiply by smooth cross-section if needed + if (urr.multiply_smooth_) { + calculate_elastic_xs(); + elastic *= micro.elastic; + capture *= (micro.absorption - micro.fission); + fission *= micro.fission; + } + + // Check for negative values + if (elastic < 0.) {elastic = 0.;} + if (fission < 0.) {fission = 0.;} + if (capture < 0.) {capture = 0.;} + + // Set elastic, absorption, fission, and total x/s. Note that the total x/s + // is calculated as a sum of partials instead of the table-provided value + micro.elastic = elastic; + micro.absorption = capture + fission; + micro.fission = fission; + micro.total = elastic + inelastic + capture + fission; + + // Determine nu-fission cross-section + if (fissionable_) { + micro.nu_fission = nu(E, EmissionMode::total) * micro.fission; + } + +} + //============================================================================== // Fortran compatibility functions //============================================================================== @@ -24,4 +543,41 @@ 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, + data::nuclides.size())); + 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(); + } +} + +extern "C" void +nuclide_calculate_urr_xs(bool use_mp, int i_nuclide, int i_temp, double E) +{ + Nuclide* nuc = data::nuclides[i_nuclide - 1].get(); + if (settings::urr_ptables_on && (nuc->urr_present_ && !use_mp)) { + if ((E > nuc->urr_data_[i_temp - 1].energy_(0)) && + (E < nuc->urr_data_[i_temp - 1].energy_( + nuc->urr_data_[i_temp - 1].n_energy_ - 1))) { + nuc->calculate_urr_xs(i_temp - 1, E); + } + } +} + } // namespace openmc diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index dcd01da175..296a4a6a27 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -23,7 +23,6 @@ module nuclide_header use stl_vector, only: VectorInt, VectorReal use string use simulation_header, only: need_depletion_rx - use urr_header, only: UrrData implicit none @@ -71,31 +70,17 @@ 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 - ! Unresolved resonance data - logical :: urr_present = .false. - integer :: urr_inelastic - type(UrrData), allocatable :: urr_data(:) - ! Multipole data logical :: mp_present = .false. type(MultipoleArray), pointer :: multipole => null() ! 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 +90,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 +170,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 @@ -206,6 +192,14 @@ module nuclide_header character(kind=C_CHAR), intent(in) :: name(*) type(C_PTR) :: path end function + + subroutine nuclide_calculate_urr_xs(use_mp, i_nuclide, i_temp, E) bind(C) + import C_BOOL, C_INT, C_DOUBLE + logical(C_BOOL), value, intent(in) :: use_mp + integer(C_INT), value, intent(in) :: i_nuclide + integer(C_INT), value, intent(in) :: i_temp + real(C_DOUBLE), value, intent(in) :: E + end subroutine nuclide_calculate_urr_xs end interface contains @@ -231,78 +225,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 +236,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 +245,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 @@ -342,12 +255,11 @@ contains integer :: i integer :: i_closest integer :: n_temperature - integer(HID_T) :: urr_group, nu_group + integer(HID_T) :: nu_group integer(HID_T) :: energy_group, energy_dset 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 +273,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 +416,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,77 +433,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. - allocate(this % urr_data(n_temperature)) - - do i = 1, n_temperature - ! Get temperature as a string - temp_str = trim(to_str(temps_to_read % data(i))) // "K" - - ! Read probability tables for i-th temperature - urr_group = open_group(group_id, 'urr/' // trim(temp_str)) - call this % urr_data(i) % from_hdf5(urr_group) - call close_group(urr_group) - - ! Check for negative values - if (any(this % urr_data(i) % prob < ZERO) .and. master) then - call warning("Negative value(s) found on probability table & - &for nuclide " // this % name // " at " // trim(temp_str)) - end if - end do - - ! if the inelastic competition flag indicates that the inelastic cross - ! section should be determined from a normal reaction cross section, we - ! need to get the index of the reaction - if (n_temperature > 0) then - if (this % urr_data(1) % inelastic_flag > 0) then - do i = 1, size(this % reactions) - if (this % reactions(i) % MT == this % urr_data(1) % inelastic_flag) then - this % urr_inelastic = i - end if - end do - - ! Abort if no corresponding inelastic reaction was found - if (this % urr_inelastic == NONE) then - call fatal_error("Could not find inelastic reaction specified on & - &unresolved resonance probability table.") - end if - end if - end if - end if - ! Check for nu-total if (object_exists(group_id, 'total_nu')) then nu_group = open_group(group_id, 'total_nu') @@ -718,7 +574,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 +593,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 @@ -759,6 +613,7 @@ contains ! Calculate nu-fission cross section do t = 1, n_temperature if (this % fissionable) then + n_grid = size(this % grid(t) % energy) do i = 1, n_grid this % xs(t) % value(XS_NU_FISSION,i) = & this % nu(this % grid(t) % energy(i), EMISSION_TOTAL) * & @@ -882,7 +737,7 @@ contains real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache - logical :: use_mp ! true if XS can be calculated with windowed multipole + logical(C_BOOL) :: use_mp ! true if XS can be calculated with windowed multipole integer :: i_temp ! index for temperature integer :: i_grid ! index on nuclide energy grid integer :: i_low ! lower logarithmic mapping index @@ -937,7 +792,7 @@ contains ! 1. physics.F90 - scatter - For inelastic scatter. ! 2. physics.F90 - sample_fission - For partial fissions. ! 3. tally.F90 - score_general - For tallying on MTxxx reactions. - ! 4. nuclide_header.F90 - calculate_urr_xs - For unresolved purposes. + ! 4. nuclide.h - calculate_urr_xs - For unresolved purposes. ! It is worth noting that none of these occur in the resolved ! resonance range, so the value here does not matter. index_temp is ! set to -1 to force a segfault in case a developer messes up and tries @@ -1080,15 +935,10 @@ contains call calculate_sab_xs(this, i_sab, E, sqrtkT, sab_frac, micro_xs) end if + ! If the particle is in the unresolved resonance range and there are ! probability tables, we need to determine cross sections from the table - - if (urr_ptables_on .and. this % urr_present .and. .not. use_mp) then - if (E > this % urr_data(i_temp) % energy(1) .and. E < this % & - urr_data(i_temp) % energy(this % urr_data(i_temp) % n_energy)) then - call calculate_urr_xs(this, i_temp, E, micro_xs) - end if - end if + call nuclide_calculate_urr_xs(use_mp, this % i_nuclide, i_temp, E) micro_xs % last_E = E micro_xs % last_sqrtkT = sqrtkT @@ -1352,191 +1202,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 -!=============================================================================== - - subroutine calculate_urr_xs(this, i_temp, E, micro_xs) - class(Nuclide), intent(in) :: this ! Nuclide object - integer, intent(in) :: i_temp ! temperature index - real(8), intent(in) :: E ! energy - type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache - - integer :: i_energy ! index for energy - integer :: i_low ! band index at lower bounding energy - integer :: i_up ! band index at upper bounding energy - integer :: threshold ! threshold energy index - real(8) :: f ! interpolation factor - real(8) :: r ! pseudo-random number - real(8) :: elastic ! elastic cross section - real(8) :: capture ! (n,gamma) cross section - real(8) :: fission ! fission cross section - real(8) :: inelastic ! inelastic cross section - - micro_xs % use_ptable = .true. - - associate (urr => this % urr_data(i_temp)) - ! determine energy table - i_energy = 1 - do - if (E < urr % energy(i_energy + 1)) exit - i_energy = i_energy + 1 - end do - - ! determine interpolation factor on table - f = (E - urr % energy(i_energy)) / & - (urr % energy(i_energy + 1) - urr % energy(i_energy)) - - ! sample probability table using the cumulative distribution - - ! Random numbers for xs calculation are sampled from a separated stream. - ! This guarantees the randomness and, at the same time, makes sure we reuse - ! random number for the same nuclide at different temperatures, therefore - ! preserving correlation of temperature in probability tables. - call prn_set_stream(STREAM_URR_PTABLE) - r = future_prn(int(this % i_nuclide, 8)) - call prn_set_stream(STREAM_TRACKING) - - i_low = 1 - do - if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit - i_low = i_low + 1 - end do - i_up = 1 - do - if (urr % prob(i_energy + 1, URR_CUM_PROB, i_up) > r) exit - i_up = i_up + 1 - end do - - ! determine elastic, fission, and capture cross sections from probability - ! table - if (urr % interp == LINEAR_LINEAR) then - elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_low) + & - f * urr % prob(i_energy + 1, URR_ELASTIC, i_up) - fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_low) + & - f * urr % prob(i_energy + 1, URR_FISSION, i_up) - capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_low) + & - f * urr % prob(i_energy + 1, URR_N_GAMMA, i_up) - elseif (urr % interp == LOG_LOG) then - ! Get logarithmic interpolation factor - f = log(E / urr % energy(i_energy)) / & - log(urr % energy(i_energy + 1) / urr % energy(i_energy)) - - ! Calculate elastic cross section/factor - elastic = ZERO - if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then - elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, & - i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, & - i_up))) - end if - - ! Calculate fission cross section/factor - fission = ZERO - if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then - fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, & - i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, & - i_up))) - end if - - ! Calculate capture cross section/factor - capture = ZERO - if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then - capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, & - i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, & - i_up))) - end if - end if - - ! Determine treatment of inelastic scattering - inelastic = ZERO - if (urr % inelastic_flag > 0) then - ! Get index on energy grid and interpolation factor - i_energy = micro_xs % index_grid - f = micro_xs % interp_factor - - ! Determine inelastic scattering cross section - associate (rx => this % reactions(this % urr_inelastic)) - threshold = rx % xs_threshold(i_temp) - if (i_energy >= threshold) then - inelastic = (ONE - f) * rx % xs(i_temp, i_energy - threshold + 1) + & - f * rx % xs(i_temp, i_energy - threshold + 2) - end if - end associate - end if - - ! Multiply by smooth cross-section if needed - if (urr % multiply_smooth) then - call this % calculate_elastic_xs(micro_xs) - elastic = elastic * micro_xs % elastic - capture = capture * (micro_xs % absorption - micro_xs % fission) - fission = fission * micro_xs % fission - end if - - ! Check for negative values - if (elastic < ZERO) elastic = ZERO - if (fission < ZERO) fission = ZERO - if (capture < ZERO) capture = ZERO - - ! Set elastic, absorption, fission, and total cross sections. Note that the - ! total cross section is calculated as sum of partials rather than using the - ! table-provided value - micro_xs % elastic = elastic - micro_xs % absorption = capture + fission - micro_xs % fission = fission - micro_xs % total = elastic + inelastic + capture + fission - - ! Determine nu-fission cross section - if (this % fissionable) then - micro_xs % nu_fission = this % nu(E, EMISSION_TOTAL) * & - micro_xs % fission - end if - end associate - - end subroutine calculate_urr_xs - !=============================================================================== ! CHECK_DATA_VERSION checks for the right version of nuclear data within HDF5 ! files @@ -1572,6 +1237,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 @@ -1581,6 +1249,7 @@ contains call nuclides(i) % clear() end do deallocate(nuclides) + call nuclides_clear() end if n_nuclides = 0 @@ -1665,9 +1334,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) @@ -1705,4 +1371,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..27d0e7b81f --- /dev/null +++ b/src/physics.cpp @@ -0,0 +1,1183 @@ +#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(); + } + + 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 + // TODO: off-by-one + constexpr int neutron = static_cast(ParticleType::neutron) - 1; + 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 + // TODO: off-by-one + constexpr int neutron = static_cast(ParticleType::neutron) - 1; + 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..1d475aaa74 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" @@ -21,13 +22,13 @@ namespace openmc { void -collision_mg(Particle* p, const MaterialMacroXS* material_xs) +collision_mg(Particle* p) { // Add to the collision counter for the particle p->n_collision++; // Sample the reaction type - sample_reaction(p, material_xs); + sample_reaction(p); // Display information about collision if ((settings::verbosity >= 10) || (simulation::trace)) { @@ -38,7 +39,7 @@ collision_mg(Particle* p, const MaterialMacroXS* material_xs) } void -sample_reaction(Particle* p, const MaterialMacroXS* material_xs) +sample_reaction(Particle* p) { // 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 @@ -47,23 +48,20 @@ 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()); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { create_fission_sites(p, p->secondary_bank, &(p->n_secondary), - MAX_SECONDARY, material_xs); + MAX_SECONDARY); } } // If survival biasing is being used, the following subroutine adjusts the // weight of the particle. Otherwise, it checks to see if absorption occurs. - if (material_xs->absorption > 0.) { - absorption(p, material_xs); + if (simulation::material_xs.absorption > 0.) { + absorption(p); } else { p->absorb_wgt = 0.; } @@ -105,7 +103,7 @@ scatter(Particle* p) void create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, - int64_t bank_array_size, const MaterialMacroXS* material_xs) + int64_t bank_array_size) { // TODO: Heat generation from fission @@ -114,8 +112,8 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, 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 * material_xs->nu_fission / - material_xs->total; + double nu_t = p->wgt / simulation::keff * weight * + simulation::material_xs.nu_fission / simulation::material_xs.total; // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -203,11 +201,12 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, } void -absorption(Particle* p, const MaterialMacroXS* material_xs) +absorption(Particle* p) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->absorb_wgt = p->wgt * material_xs->absorption / material_xs->total; + p->absorb_wgt = p->wgt * + simulation::material_xs.absorption / simulation::material_xs.total; // Adjust weight of particle by the probability of absorption p->wgt -= p->absorb_wgt; @@ -215,13 +214,15 @@ absorption(Particle* p, const MaterialMacroXS* material_xs) // Score implicit absorpion estimate of keff #pragma omp atomic - global_tally_absorption += p->absorb_wgt * material_xs->nu_fission / - material_xs->absorption; + global_tally_absorption += p->absorb_wgt * + simulation::material_xs.nu_fission / + simulation::material_xs.absorption; } else { - if (material_xs->absorption > prn() * material_xs->total) { + if (simulation::material_xs.absorption > + prn() * simulation::material_xs.total) { #pragma omp atomic - global_tally_absorption += p->wgt * material_xs->nu_fission / - material_xs->absorption; + global_tally_absorption += p->wgt * simulation::material_xs.nu_fission / + simulation::material_xs.absorption; p->alive = false; p->event = EVENT_ABSORB; } 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 3ff6475eb4..3cfbd4286e 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -5,7 +5,9 @@ #include #include +#ifdef _OPENMP #include +#endif #include "openmc/capi.h" #include "openmc/constants.h" @@ -801,18 +803,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 4a4f632ebf..3591953cb8 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 @@ -22,11 +26,9 @@ namespace openmc { // data/functions from Fortran side extern "C" void accumulate_tallies(); -extern "C" void allocate_banks(); extern "C" void allocate_tally_results(); extern "C" void check_triggers(); 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(); @@ -91,6 +93,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) @@ -272,6 +275,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 @@ -438,13 +471,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 0d33d52d1c..d69217d8b7 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -34,16 +34,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 @@ -412,7 +410,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 @@ -588,8 +586,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/src/tracking.F90 b/src/tracking.F90 index e322f31f25..d345d7cc0b 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -33,10 +33,9 @@ module tracking implicit none interface - subroutine collision_mg(p, material_xs) bind(C) - import Particle, C_DOUBLE, MaterialMacroXS + subroutine collision_mg(p) bind(C) + import Particle, C_DOUBLE type(Particle), intent(inout) :: p - type(MaterialMacroXS), intent(in) :: material_xs end subroutine collision_mg end interface @@ -221,7 +220,7 @@ contains material_xs % nu_fission / material_xs % total end if - ! score surface current tallies -- this has to be done before the collision + ! Score surface current tallies -- this has to be done before the collision ! since the direction of the particle will change and we need to use the ! pre-collision direction to figure out what mesh surfaces were crossed @@ -234,7 +233,7 @@ contains if (run_CE) then call collision(p) else - call collision_mg(p, material_xs) + call collision_mg(p) end if ! Score collision estimator tallies -- this is done after a collision diff --git a/src/urr.cpp b/src/urr.cpp new file mode 100644 index 0000000000..28ad2505f0 --- /dev/null +++ b/src/urr.cpp @@ -0,0 +1,31 @@ +#include "openmc/urr.h" + +#include + +namespace openmc { + +UrrData::UrrData(hid_t group_id) +{ + // Read interpolation and other flags + int interp_temp; + read_attribute(group_id, "interpolation", interp_temp); + interp_ = static_cast(interp_temp); + + // read the metadata + read_attribute(group_id, "inelastic", inelastic_flag_); + read_attribute(group_id, "absorption", absorption_flag_); + int temp_multiply_smooth; + read_attribute(group_id, "multiply_smooth", temp_multiply_smooth); + multiply_smooth_ = (temp_multiply_smooth == 1); + + // read the energies at which tables exist + read_dataset(group_id, "energy", energy_); + + // Set n_energy_ + n_energy_ = energy_.shape()[0]; + + // Read URR tables + read_dataset(group_id, "table", prob_); +} + +} \ No newline at end of file diff --git a/src/urr_header.F90 b/src/urr_header.F90 deleted file mode 100644 index 0ec9e89deb..0000000000 --- a/src/urr_header.F90 +++ /dev/null @@ -1,72 +0,0 @@ -module urr_header - - use hdf5_interface, only: read_attribute, open_dataset, read_dataset, & - close_dataset, get_shape, HID_T, HSIZE_T - - implicit none - -!=============================================================================== -! URRDATA contains probability tables for the unresolved resonance range. -!=============================================================================== - - type UrrData - integer :: n_energy ! # of incident neutron energies - integer :: n_prob ! # of probabilities - integer :: interp ! inteprolation (2=lin-lin, 5=log-log) - integer :: inelastic_flag ! inelastic competition flag - integer :: absorption_flag ! other absorption flag - logical :: multiply_smooth ! multiply by smooth cross section? - real(8), allocatable :: energy(:) ! incident energies - real(8), allocatable :: prob(:,:,:) ! actual probabibility tables - contains - procedure :: from_hdf5 => urr_from_hdf5 - end type UrrData - -contains - - subroutine urr_from_hdf5(this, group_id) - class(UrrData), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer :: i, j, k - integer(HID_T) :: energy - integer(HID_T) :: table - integer(HSIZE_T) :: dims(1) - integer(HSIZE_T) :: dims3(3) - real(8), allocatable :: temp(:,:,:) - - ! Read interpolation and other flags - call read_attribute(this % interp, group_id, 'interpolation') - call read_attribute(this % inelastic_flag, group_id, 'inelastic') - call read_attribute(this % absorption_flag, group_id, 'absorption') - call read_attribute(i, group_id, 'multiply_smooth') - this % multiply_smooth = (i == 1) - - ! Read energies at which tables exist - energy = open_dataset(group_id, 'energy') - call get_shape(energy, dims) - this % n_energy = int(dims(1), 4) - allocate(this % energy(this % n_energy)) - call read_dataset(this % energy, energy) - call close_dataset(energy) - - ! Read URR tables - table = open_dataset(group_id, 'table') - call get_shape(table, dims3) - this % n_prob = int(dims3(1), 4) - allocate(temp(this % n_prob, 6, this % n_energy)) - call read_dataset(temp, table) - call close_dataset(table) - - ! Swap first and last indices - allocate(this % prob(this % n_energy, 6, this % n_prob)) - do i = 1, this % n_energy - do j = 1, 6 - do k = 1, this % n_prob - this % prob(i, j, k) = temp(k, j, i) - end do - end do - end do - end subroutine urr_from_hdf5 - -end module urr_header 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']