From dc2f2f8ae0927edb9696fef4763c9382728d396f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 3 Jan 2019 13:53:21 -0600 Subject: [PATCH 01/14] Make micro_photon_xs available on C++ side --- include/openmc/photon.h | 23 +++++++++++++++++++++-- src/nuclide.cpp | 3 +++ src/photon.cpp | 6 ++++-- src/photon_header.F90 | 31 +++++++++++++++++++++---------- 4 files changed, 49 insertions(+), 14 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 3f36decdb3..8ed1431d0c 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -22,16 +22,35 @@ public: int Z_; //! Atomic number }; +//============================================================================== +//! Cached microscopic photon cross sections for a particular element at the +//! current energy +//============================================================================== + +struct ElementMicroXS { + int index_grid; //!< index on element energy grid + double last_E {0.0}; //!< last evaluated energy in [eV] + double interp_factor; //!< interpolation factor on energy grid + double total; //!< microscopic total photon xs + double coherent; //!< microscopic coherent xs + double incoherent; //!< microscopic incoherent xs + double photoelectric; //!< microscopic photoelectric xs + double pair_production; //!< microscopic pair production xs +}; + //============================================================================== // Global variables //============================================================================== namespace data { - extern std::vector elements; - } // namespace data +namespace simulation { +extern ElementMicroXS* micro_photon_xs; +#pragma omp threadprivate(micro_photon_xs) +} // namespace simulation + } // namespace openmc #endif // OPENMC_PHOTON_H diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 030185f1b4..075e6f55b2 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -5,6 +5,7 @@ #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" +#include "openmc/photon.h" #include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/settings.h" @@ -935,12 +936,14 @@ extern "C" bool multipole_in_range(Nuclide* nuc, double E) extern "C" void nuclides_clear() { data::nuclides.clear(); } extern "C" NuclideMicroXS* micro_xs_ptr(); +extern "C" ElementMicroXS* micro_photon_xs_ptr(); void set_micro_xs() { #pragma omp parallel { simulation::micro_xs = micro_xs_ptr(); + simulation::micro_photon_xs = micro_photon_xs_ptr(); } } diff --git a/src/photon.cpp b/src/photon.cpp index b36f29a5c3..a13d287bf1 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -9,11 +9,13 @@ namespace openmc { //============================================================================== namespace data { - std::vector elements; - } // namespace data +namespace simulation { +ElementMicroXS* micro_photon_xs; +} // namespace simulation + //============================================================================== // PhotonInteraction implementation //============================================================================== diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 9da958a6f4..4be73831d0 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -1,5 +1,7 @@ module photon_header + use, intrinsic :: ISO_C_BINDING + use algorithm, only: binary_search use constants use dict_header, only: DictIntInt, DictCharInt @@ -92,18 +94,18 @@ module photon_header ! particular element at the current energy !=============================================================================== - type ElementMicroXS - integer :: index_grid ! index on element energy grid - real(8) :: last_E = ZERO ! last evaluated energy - real(8) :: interp_factor ! interpolation factor on energy grid - real(8) :: total ! microscropic total photon xs - real(8) :: coherent ! microscopic coherent xs - real(8) :: incoherent ! microscopic incoherent xs - real(8) :: photoelectric ! microscopic photoelectric xs - real(8) :: pair_production ! microscopic pair production xs + type, bind(C) :: ElementMicroXS + integer(C_INT) :: index_grid ! index on element energy grid + real(C_DOUBLE) :: last_E = ZERO ! last evaluated energy + real(C_DOUBLE) :: interp_factor ! interpolation factor on energy grid + real(C_DOUBLE) :: total ! microscropic total photon xs + real(C_DOUBLE) :: coherent ! microscopic coherent xs + real(C_DOUBLE) :: incoherent ! microscopic incoherent xs + real(C_DOUBLE) :: photoelectric ! microscopic photoelectric xs + real(C_DOUBLE) :: pair_production ! microscopic pair production xs end type ElementMicroXS - type(ElementMicroXS), allocatable :: micro_photon_xs(:) ! Cache for each element + type(ElementMicroXS), allocatable, target :: micro_photon_xs(:) ! Cache for each element !$omp threadprivate(micro_photon_xs) contains @@ -517,4 +519,13 @@ contains if (allocated(ttb)) deallocate(ttb) end subroutine free_memory_photon + function micro_photon_xs_ptr() result(ptr) bind(C) + type(C_PTR) :: ptr + if (size(micro_photon_xs) > 0) then + ptr = C_LOC(micro_photon_xs(1)) + else + ptr = C_NULL_PTR + end if + end function + end module photon_header From 6c3f35d0197fe0c826545a99c1043eb5bc56e0ca Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 4 Jan 2019 08:24:52 -0600 Subject: [PATCH 02/14] Read all photon data into PhotonInteraction class --- include/openmc/constants.h | 14 +-- include/openmc/photon.h | 63 ++++++++++ src/photon.cpp | 249 +++++++++++++++++++++++++++++++++++++ 3 files changed, 319 insertions(+), 7 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index d29abcdcfa..375ebc2767 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -4,8 +4,8 @@ #ifndef OPENMC_CONSTANTS_H #define OPENMC_CONSTANTS_H -#include #include +#include #include #include @@ -96,12 +96,12 @@ constexpr double N_AVOGADRO {0.6022140857}; // Avogadro's number in 10^24/ constexpr double K_BOLTZMANN {8.6173303e-5}; // Boltzmann constant in eV/K // Electron subshell labels -constexpr char SUBSHELLS[][4] { - "K ", "L1 ", "L2 ", "L3 ", "M1 ", "M2 ", "M3 ", "M4 ", "M5 ", - "N1 ", "N2 ", "N3 ", "N4 ", "N5 ", "N6 ", "N7 ", "O1 ", "O2 ", - "O3 ", "O4 ", "O5 ", "O6 ", "O7 ", "O8 ", "O9 ", "P1 ", "P2 ", - "P3 ", "P4 ", "P5 ", "P6 ", "P7 ", "P8 ", "P9 ", "P10", "P11", - "Q1 ", "Q2 ", "Q3 " +constexpr std::array SUBSHELLS = { + "K", "L1", "L2", "L3", "M1", "M2", "M3", "M4", "M5", + "N1", "N2", "N3", "N4", "N5", "N6", "N7", "O1", "O2", + "O3", "O4", "O5", "O6", "O7", "O8", "O9", "P1", "P2", + "P3", "P4", "P5", "P6", "P7", "P8", "P9", "P10", "P11", + "Q1", "Q2", "Q3" }; // Void material diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 8ed1431d0c..bc3f33e144 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -1,9 +1,13 @@ #ifndef OPENMC_PHOTON_H #define OPENMC_PHOTON_H +#include "openmc/endf.h" + #include +#include "xtensor/xtensor.hpp" #include +#include #include namespace openmc { @@ -12,6 +16,24 @@ namespace openmc { //! Photon interaction data for a single element //============================================================================== +class ElectronSubshell { +public: + // Constructors + ElectronSubshell() { }; + + int index_subshell; //!< index in SUBSHELLS + int threshold; + double n_electrons; + double binding_energy; + xt::xtensor cross_section; + + // Transition data + int n_transitions; + xt::xtensor transition_subshells; + xt::xtensor transition_energy; + xt::xtensor transition_probability; +}; + class PhotonInteraction { public: // Constructors @@ -20,6 +42,40 @@ public: // Data members std::string name_; //! Name of element, e.g. "Zr" int Z_; //! Atomic number + + // Microscopic cross sections + xt::xtensor energy_; + xt::xtensor coherent_; + xt::xtensor incoherent_; + xt::xtensor photoelectric_total_; + xt::xtensor pair_production_total_; + xt::xtensor pair_production_electron_; + xt::xtensor pair_production_nuclear_; + + // Form factors + Tabulated1D incoherent_form_factor_; + Tabulated1D coherent_int_form_factor_; + Tabulated1D coherent_anomalous_real_; + Tabulated1D coherent_anomalous_imag_; + + // Photoionization and atomic relaxation data + std::unordered_map shell_map_; // Given a shell designator, e.g. 3, this + // dictionary gives an index in shells_ + std::vector shells_; + + // Compton profile data + xt::xtensor profile_pdf_; + xt::xtensor profile_cdf_; + xt::xtensor binding_energy_; + xt::xtensor electron_pdf_; + + // Stopping power data + double I_; // mean excitation energy + xt::xtensor stopping_power_collision_; + xt::xtensor stopping_power_radiative_; + + // Bremsstrahlung scaled DCS + xt::xtensor dcs_; }; //============================================================================== @@ -43,7 +99,14 @@ struct ElementMicroXS { //============================================================================== namespace data { + +extern xt::xtensor compton_profile_pz; //! Compton profile momentum grid +extern xt::xtensor ttb_e_grid; //! energy T of incident electron in [eV] +extern xt::xtensor ttb_k_grid; //! reduced energy W/T of emitted photon + +//! Photon interaction data for each element extern std::vector elements; + } // namespace data namespace simulation { diff --git a/src/photon.cpp b/src/photon.cpp index a13d287bf1..edb0de4e47 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -1,6 +1,14 @@ #include "openmc/photon.h" +#include "openmc/constants.h" #include "openmc/hdf5_interface.h" +#include "openmc/particle.h" +#include "openmc/search.h" +#include "openmc/settings.h" + +#include "xtensor/xbuilder.hpp" +#include "xtensor/xoperation.hpp" +#include "xtensor/xview.hpp" namespace openmc { @@ -9,7 +17,13 @@ namespace openmc { //============================================================================== namespace data { + +xt::xtensor compton_profile_pz; +xt::xtensor ttb_e_grid; +xt::xtensor ttb_k_grid; + std::vector elements; + } // namespace data namespace simulation { @@ -25,7 +39,242 @@ PhotonInteraction::PhotonInteraction(hid_t group) // Get name of nuclide from group, removing leading '/' name_ = object_name(group).substr(1); + // Get atomic number read_attribute(group, "Z", Z_); + + // Determine number of energies and read energy grid + read_dataset(group, "energy", energy_); + + // Read coherent scattering + hid_t rgroup = open_group(group, "coherent"); + read_dataset(rgroup, "xs", coherent_); + + hid_t dset = open_dataset(rgroup, "integrated_scattering_factor"); + coherent_int_form_factor_ = Tabulated1D{dset}; + close_dataset(dset); + + if (object_exists(group, "anomalous_real")) { + dset = open_dataset(rgroup, "anomalous_real"); + coherent_anomalous_real_ = Tabulated1D{dset}; + close_dataset(dset); + } + + if (object_exists(group, "anomalous_imag")) { + dset = open_dataset(rgroup, "anomalous_imag"); + coherent_anomalous_imag_ = Tabulated1D{dset}; + close_dataset(dset); + } + close_group(rgroup); + + // Read incoherent scattering + rgroup = open_group(group, "incoherent"); + read_dataset(rgroup, "xs", incoherent_); + dset = open_dataset(rgroup, "scattering_factor"); + incoherent_form_factor_ = Tabulated1D{dset}; + close_dataset(dset); + close_group(rgroup); + + // Read pair production + rgroup = open_group(group, "pair_production_electron"); + read_dataset(rgroup, "xs", pair_production_electron_); + close_group(rgroup); + + // Read pair production + if (object_exists(group, "pair_production_nuclear")) { + rgroup = open_group(group, "pair_production_nuclear"); + read_dataset(rgroup, "xs", pair_production_nuclear_); + close_group(rgroup); + } else { + pair_production_nuclear_ = xt::zeros_like(energy_); + } + + // Read photoelectric + rgroup = open_group(group, "photoelectric"); + read_dataset(rgroup, "xs", photoelectric_total_); + close_group(rgroup); + + // Read subshell photoionization cross section and atomic relaxation data + rgroup = open_group(group, "subshells"); + std::vector designators; + read_attribute(rgroup, "designators", designators); + auto n_shell = designators.size(); + for (int i = 0; i < n_shell; ++i) { + const auto& designator {designators[i]}; + + // TODO: Move to ElectronSubshell constructor + + // Add empty shell + shells_.emplace_back(); + + // Create mapping from designator to index + int j = 0; + for (const auto& subshell : SUBSHELLS) { + if (designator == subshell) { + shell_map_[j] = i; + shells_[i].index_subshell = j; + break; + } + ++j; + } + + // Read binding energy and number of electrons + auto& shell {shells_.back()}; + hid_t tgroup = open_group(rgroup, designator.c_str()); + read_attribute(tgroup, "binding_energy", shell.binding_energy); + read_attribute(tgroup, "num_electrons", shell.n_electrons); + + // Read subshell cross section + dset = open_dataset(tgroup, "xs"); + read_attribute(dset, "threshold_idx", j); + // TODO: off-by-one + shell.threshold = j - 1; + close_dataset(dset); + read_dataset(tgroup, "xs", shell.cross_section); + + auto& xs = shell.cross_section; + xs = xt::where(xs > 0.0, xt::log(xs), -500.0); + + if (object_exists(tgroup, "transitions")) { + // Determine dimensions of transitions + dset = open_dataset(tgroup, "transitions"); + auto dims = object_shape(dset); + close_dataset(dset); + + int n_transition = dims[0]; + shell.n_transitions = n_transition; + if (n_transition > 0) { + xt::xtensor matrix; + read_dataset(tgroup, "transitions", matrix); + + shell.transition_subshells = xt::view(matrix, xt::all(), xt::range(0, 2)); + shell.transition_energy = xt::view(matrix, xt::all(), 2); + shell.transition_probability = xt::view(matrix, xt::all(), 3); + shell.transition_probability /= xt::sum(shell.transition_probability)(); + } + } + close_group(tgroup); + } + close_group(rgroup); + + // Determine number of electron shells + rgroup = open_group(group, "compton_profiles"); + + // Read electron shell PDF and binding energies + read_dataset(rgroup, "num_electrons", electron_pdf_); + electron_pdf_ /= xt::sum(electron_pdf_); + read_dataset(rgroup, "binding_energy", binding_energy_); + + // Read Compton profiles + read_dataset(rgroup, "J", profile_pdf_); + + // Get Compton profile momentum grid. By deafult, an xtensor has a size of 1. + // TODO: Update version of xtensor and change to 0 + if (data::compton_profile_pz.size() == 1) { + read_dataset(rgroup, "pz", data::compton_profile_pz); + } + close_group(rgroup); + + // Create Compton profile CDF + auto n_profile = data::compton_profile_pz.size(); + profile_cdf_ = xt::empty({n_shell, n_profile}); + for (int i = 0; i < n_shell; ++i) { + double c = 0.0; + profile_cdf_(i,0) = 0.0; + for (int j = 0; j < n_profile - 1; ++j) { + c += 0.5*(data::compton_profile_pz(j+1) - data::compton_profile_pz(j)) * + (profile_pdf_(i,j) + profile_pdf_(i,j+1)); + profile_cdf_(i,j+1) = c; + } + } + + // Calculate total pair production + pair_production_total_ = pair_production_nuclear_ + pair_production_electron_; + + if (settings::electron_treatment == ELECTRON_TTB) { + // Read bremsstrahlung scaled DCS + rgroup = open_group(group, "bremsstrahlung"); + read_dataset(rgroup, "dcs", dcs_); + auto n_e = dcs_.shape()[0]; + auto n_k = dcs_.shape()[1]; + + // Get energy grids used for bremsstrahlung DCS and for stopping powers + xt::xtensor electron_energy; + read_dataset(rgroup, "electron_energy", electron_energy); + if (data::ttb_k_grid.size() == 1) { + read_dataset(rgroup, "photon_energy", data::ttb_k_grid); + } + close_group(rgroup); + + // Read stopping power data + if (Z_ < 99) { + rgroup = open_group(group, "stopping_powers"); + read_dataset(rgroup, "s_collision", stopping_power_collision_); + read_dataset(rgroup, "s_radiative", stopping_power_radiative_); + read_attribute(rgroup, "I", I_); + close_group(rgroup); + } + + // Truncate the bremsstrahlung data at the cutoff energy + int photon = static_cast(ParticleType::photon) - 1; + const auto& E {electron_energy}; + double cutoff = settings::energy_cutoff[photon]; + if (cutoff > E(0)) { + size_t i_grid = lower_bound_index(E.cbegin(), E.cend(), + settings::energy_cutoff[photon]); + + // calculate interpolation factor + double f = (std::log(cutoff) - std::log(E(i_grid))) / + (std::log(E(i_grid+1)) - std::log(E(i_grid))); + + // Interpolate collision stopping power at the cutoff energy and truncate + auto& s_col {stopping_power_collision_}; + double y = std::exp(std::log(s_col(i_grid)) + f*(std::log(s_col(i_grid+1)) - + std::log(s_col(i_grid)))); + xt::xtensor frst {y}; + stopping_power_collision_ = xt::concatenate(xt::xtuple( + frst, xt::view(s_col, xt::range(i_grid+1, n_e)))); + + // Interpolate radiative stopping power at the cutoff energy and truncate + auto& s_rad {stopping_power_radiative_}; + y = std::exp(std::log(s_rad(i_grid)) + f*(std::log(s_rad(i_grid+1)) - + std::log(s_rad(i_grid)))); + frst(0) = y; + stopping_power_radiative_ = xt::concatenate(xt::xtuple( + frst, xt::view(s_rad, xt::range(i_grid+1, n_e)))); + + // Interpolate bremsstrahlung DCS at the cutoff energy and truncate + xt::xtensor dcs = xt::empty({n_e - i_grid, n_k}); + for (int i = 0; i < n_k; ++i) { + y = std::exp(std::log(dcs_(i_grid,i)) + + f*(std::log(dcs_(i_grid+1,i)) - std::log(dcs_(i_grid,i)))); + auto col_i = xt::view(dcs, xt::all(), i); + col_i(0) = y; + for (int j = i_grid + 1; j < n_e; ++j) { + col_i(j - i_grid) = dcs_(j, i); + } + } + dcs_ = dcs; + + frst(0) = cutoff; + electron_energy = xt::concatenate(xt::xtuple( + frst, xt::view(electron_energy, xt::range(i_grid+1, n_e)))); + } + + // Set incident particle energy grid + if (data::ttb_e_grid.size() == 1) { + data::ttb_e_grid = electron_energy; + } + } + + // Take logarithm of energies and cross sections since they are log-log + // interpolated + energy_ = xt::log(energy_); + coherent_ = xt::where(coherent_ > 0.0, xt::log(coherent_), -500.0); + incoherent_ = xt::where(incoherent_ > 0.0, xt::log(incoherent_), -500.0); + photoelectric_total_ = xt::where(photoelectric_total_ > 0.0, + xt::log(photoelectric_total_), -500.0); + pair_production_total_ = xt::where(pair_production_total_ > 0.0, + xt::log(pair_production_total_), -500.0); } //============================================================================== From 2e6ac034330c8b4d6e66d6cb911a69d64cfcc297 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 4 Jan 2019 16:27:45 -0600 Subject: [PATCH 03/14] Convert all photon physics except bremsstrahlung (not done yet) --- include/openmc/endf.h | 4 + include/openmc/math_functions.h | 2 +- include/openmc/photon.h | 24 ++ include/openmc/physics.h | 4 +- src/math_functions.cpp | 2 +- src/photon.cpp | 523 ++++++++++++++++++++++++++++++- src/photon_physics.F90 | 529 +------------------------------- src/physics.F90 | 188 +----------- src/physics.cpp | 313 ++++++++++--------- 9 files changed, 717 insertions(+), 872 deletions(-) diff --git a/include/openmc/endf.h b/include/openmc/endf.h index ae480f8cfc..7ec078e93f 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -78,6 +78,10 @@ public: //! \param[in] x independent variable //! \return Function evaluated at x double operator()(double x) const; + + // Accessors + const std::vector& x() const { return x_; } + const std::vector& y() const { return y_; } private: std::size_t n_regions_ {0}; //!< number of interpolation regions std::vector nbt_; //!< values separating interpolation regions diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index bc5ba93083..f80ed6d3b8 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -129,7 +129,7 @@ extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]); //! is passed //============================================================================== -extern "C" void rotate_angle_c(double uvw[3], double mu, double* phi); +extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi); Direction rotate_angle(Direction u, double mu, double* phi); diff --git a/include/openmc/photon.h b/include/openmc/photon.h index bc3f33e144..fce6a3626c 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -2,12 +2,14 @@ #define OPENMC_PHOTON_H #include "openmc/endf.h" +#include "openmc/particle.h" #include #include "xtensor/xtensor.hpp" #include #include +#include // for pair #include namespace openmc { @@ -39,6 +41,17 @@ public: // Constructors PhotonInteraction(hid_t group); + // Methods + void compton_scatter(double alpha, bool doppler, double* alpha_out, + double* mu, int* i_shell) const; + + double rayleigh_scatter(double alpha) const; + + void pair_production(double alpha, double* E_electron, double* E_positron, + double* mu_electron, double* mu_positron) const; + + void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const; + // Data members std::string name_; //! Name of element, e.g. "Zr" int Z_; //! Atomic number @@ -76,6 +89,9 @@ public: // Bremsstrahlung scaled DCS xt::xtensor dcs_; + +private: + void compton_doppler(double alpha, double mu, double* E_out, int* i_shell) const; }; //============================================================================== @@ -94,6 +110,14 @@ struct ElementMicroXS { double pair_production; //!< microscopic pair production xs }; +//============================================================================== +// Non-member functions +//============================================================================== + +std::pair klein_nishina(double alpha); + +void thick_target_bremsstrahlung(Particle& p, double* E_lost); + //============================================================================== // Global variables //============================================================================== diff --git a/include/openmc/physics.h b/include/openmc/physics.h index a51ed91b8b..8a3d7c4ef1 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -22,7 +22,7 @@ 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); +void sample_photon_reaction(Particle* p); //! Terminates the particle and either deposits all energy locally //! (electron_treatment = ELECTRON_LED) or creates secondary bremsstrahlung @@ -44,7 +44,7 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat); 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); +extern "C" int sample_element(Particle* p); Reaction* sample_fission(int i_nuclide, double E); diff --git a/src/math_functions.cpp b/src/math_functions.cpp index dcf8d38b16..041b6ff57a 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -617,7 +617,7 @@ void calc_zn_rad(int n, double rho, double zn_rad[]) { } -void rotate_angle_c(double uvw[3], double mu, double* phi) { +void rotate_angle_c(double uvw[3], double mu, const double* phi) { // Copy original directional cosines double u0 = uvw[0]; // original cosine in x direction double v0 = uvw[1]; // original cosine in y direction diff --git a/src/photon.cpp b/src/photon.cpp index edb0de4e47..56ff2042dd 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -3,6 +3,7 @@ #include "openmc/constants.h" #include "openmc/hdf5_interface.h" #include "openmc/particle.h" +#include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/settings.h" @@ -10,6 +11,10 @@ #include "xtensor/xoperation.hpp" #include "xtensor/xview.hpp" +#include +#include +#include // for tie + namespace openmc { //============================================================================== @@ -107,7 +112,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) shells_.emplace_back(); // Create mapping from designator to index - int j = 0; + int j = 1; for (const auto& subshell : SUBSHELLS) { if (designator == subshell) { shell_map_[j] = i; @@ -126,8 +131,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) // Read subshell cross section dset = open_dataset(tgroup, "xs"); read_attribute(dset, "threshold_idx", j); - // TODO: off-by-one - shell.threshold = j - 1; + shell.threshold = j; close_dataset(dset); read_dataset(tgroup, "xs", shell.cross_section); @@ -151,6 +155,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) shell.transition_probability = xt::view(matrix, xt::all(), 3); shell.transition_probability /= xt::sum(shell.transition_probability)(); } + } else { + shell.n_transitions = 0; } close_group(tgroup); } @@ -277,6 +283,517 @@ PhotonInteraction::PhotonInteraction(hid_t group) xt::log(pair_production_total_), -500.0); } +void PhotonInteraction::compton_scatter(double alpha, bool doppler, + double* alpha_out, double* mu, int* i_shell) const +{ + double form_factor_xmax = 0.0; + while (true) { + // Sample Klein-Nishina distribution for trial energy and angle + std::tie(*alpha_out, *mu) = klein_nishina(alpha); + + // Note that the parameter used here does not correspond exactly to the + // momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the + // parameter as defined by Hubbell, where the actual data comes from + double x = MASS_ELECTRON_EV/PLANCK_C*alpha*std::sqrt(0.5*(1.0 - *mu)); + + // Calculate S(x, Z) and S(x_max, Z) + double form_factor_x = incoherent_form_factor_(x); + if (form_factor_xmax == 0.0) { + form_factor_xmax = incoherent_form_factor_(MASS_ELECTRON_EV/PLANCK_C*alpha); + } + + // Perform rejection on form factor + if (prn() < form_factor_x / form_factor_xmax) { + if (doppler) { + double E_out; + this->compton_doppler(alpha, *mu, &E_out, i_shell); + *alpha_out = E_out/MASS_ELECTRON_EV; + } else { + *i_shell = -1; + } + break; + } + } +} + +void PhotonInteraction::compton_doppler(double alpha, double mu, + double* E_out, int* i_shell) const +{ + auto n = data::compton_profile_pz.size(); + + int shell; // index for shell + while (true) { + // Sample electron shell + double rn = prn(); + double c = 0.0; + for (shell = 0; shell < electron_pdf_.size(); ++shell) { + c += electron_pdf_(shell); + if (rn < c) break; + } + + // Determine binding energy of shell + double E_b = binding_energy_(shell); + + // Determine p_z,max + double E = alpha*MASS_ELECTRON_EV; + if (E < E_b) { + *E_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV; + break; + } + + double pz_max = -FINE_STRUCTURE*(E_b - (E - E_b)*alpha*(1.0 - mu)) / + std::sqrt(2.0*E*(E - E_b)*(1.0 - mu) + E_b*E_b); + if (pz_max < 0.0) { + *E_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV; + break; + } + + // Determine profile cdf value corresponding to p_z,max + double c_max; + if (pz_max > data::compton_profile_pz(n - 1)) { + c_max = profile_cdf_(shell, n - 1); + } else { + int i = lower_bound_index(data::compton_profile_pz.cbegin(), + data::compton_profile_pz.cend(), pz_max); + double pz_l = data::compton_profile_pz(i); + double pz_r = data::compton_profile_pz(i + 1); + double p_l = profile_pdf_(shell, i); + double p_r = profile_pdf_(shell, i + 1); + double c_l = profile_cdf_(shell, i); + if (pz_l == pz_r) { + c_max = c_l; + } else if (p_l == p_r) { + c_max = c_l + (pz_max - pz_l)*p_l; + } else { + double m = (p_l - p_r)/(pz_l - pz_r); + c_max = c_l + (std::pow((m*(pz_max - pz_l) + p_l), 2) - p_l*p_l)/(2.0*m); + } + } + + // Sample value on bounded cdf + c = prn()*c_max; + + // Determine pz corresponding to sampled cdf value + auto cdf_shell = xt::view(profile_cdf_, shell, xt::all()); + int i = lower_bound_index(cdf_shell.cbegin(), cdf_shell.cend(), c); + double pz_l = data::compton_profile_pz(i); + double pz_r = data::compton_profile_pz(i + 1); + double p_l = profile_pdf_(shell, i); + double p_r = profile_pdf_(shell, i + 1); + double c_l = profile_cdf_(shell, i); + double pz; + if (pz_l == pz_r) { + pz = pz_l; + } else if (p_l == p_r) { + pz = pz_l + (c - c_l)/p_l; + } else { + double m = (p_l - p_r)/(pz_l - pz_r); + pz = pz_l + (std::sqrt(p_l*p_l + 2.0*m*(c - c_l)) - p_l)/m; + } + + // Determine outgoing photon energy corresponding to electron momentum + // (solve Eq. 39 in LA-UR-04-0487 for E') + double momentum_sq = std::pow((pz/FINE_STRUCTURE), 2); + double f = 1.0 + alpha*(1.0 - mu); + double a = momentum_sq - f*f; + double b = 2.0*E*(f - momentum_sq*mu); + c = E*E*(momentum_sq - 1.0); + + double quad = b*b - 4.0*a*c; + if (quad < 0) { + *E_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV; + break; + } + quad = std::sqrt(quad); + double E_out1 = -(b + quad)/(2.0*a); + double E_out2 = -(b - quad)/(2.0*a); + + // Determine solution to quadratic equation that is positive + if (E_out1 > 0.0) { + if (E_out2 > 0.0) { + // If both are positive, pick one at random + *E_out = prn() < 0.5 ? E_out1 : E_out2; + } else { + *E_out = E_out1; + } + } else { + if (E_out2 > 0.0) { + *E_out = E_out2; + } else { + // No positive solution -- resample + continue; + } + } + if (*E_out < E - E_b) break; + } + + *i_shell = shell; +} + +double PhotonInteraction::rayleigh_scatter(double alpha) const +{ + double mu; + while (true) { + // Determine maximum value of x^2 + double x2_max = std::pow(MASS_ELECTRON_EV/PLANCK_C*alpha, 2); + + // Determine F(x^2_max, Z) + double F_max = coherent_int_form_factor_(x2_max); + + // Sample cumulative distribution + double F = prn()*F_max; + + // Determine x^2 corresponding to F + const auto& x {coherent_int_form_factor_.x()}; + const auto& y {coherent_int_form_factor_.y()}; + int i = lower_bound_index(y.cbegin(), y.cend(), F); + double r = (F - y[i]) / (y[i+1] - y[i]); + double x2 = x[i] + r*(x[i+1] - x[i]); + + // Calculate mu + mu = 1.0 - 2.0*x2/x2_max; + + if (prn() < 0.5*(1.0 + mu*mu)) break; + } + return mu; +} + +void PhotonInteraction::pair_production(double alpha, double* E_electron, + double* E_positron, double* mu_electron, double* mu_positron) const +{ + constexpr double r[] { + 122.81, 73.167, 69.228, 67.301, 64.696, 61.228, + 57.524, 54.033, 50.787, 47.851, 46.373, 45.401, + 44.503, 43.815, 43.074, 42.321, 41.586, 40.953, + 40.524, 40.256, 39.756, 39.144, 38.462, 37.778, + 37.174, 36.663, 35.986, 35.317, 34.688, 34.197, + 33.786, 33.422, 33.068, 32.740, 32.438, 32.143, + 31.884, 31.622, 31.438, 31.142, 30.950, 30.758, + 30.561, 30.285, 30.097, 29.832, 29.581, 29.411, + 29.247, 29.085, 28.930, 28.721, 28.580, 28.442, + 28.312, 28.139, 27.973, 27.819, 27.675, 27.496, + 27.285, 27.093, 26.911, 26.705, 26.516, 26.304, + 26.108, 25.929, 25.730, 25.577, 25.403, 25.245, + 25.100, 24.941, 24.790, 24.655, 24.506, 24.391, + 24.262, 24.145, 24.039, 23.922, 23.813, 23.712, + 23.621, 23.523, 23.430, 23.331, 23.238, 23.139, + 23.048, 22.967, 22.833, 22.694, 22.624, 22.545, + 22.446, 22.358, 22.264}; + + // The reduced screening radius r is the ratio of the screening radius to + // the Compton wavelength of the electron, where the screening radius is + // obtained under the assumption that the Coulomb field of the nucleus is + // exponentially screened by atomic electrons. This allows us to use a + // simplified atomic form factor and analytical approximations of the + // screening functions in the pair production DCS instead of computing the + // screening functions numerically. The reduced screening radii above for + // Z = 1-99 come from F. Salvat, J. M. Fernández-Varea, and J. Sempau, + // "PENELOPE-2011: A Code System for Monte Carlo Simulation of Electron and + // Photon Transport," OECD-NEA, Issy-les-Moulineaux, France (2011). + + // Compute the minimum and maximum values of the electron reduced energy, + // i.e. the fraction of the photon energy that is given to the electron + double e_min = 1.0/alpha; + double e_max = 1.0 - 1.0/alpha; + + // Compute the high-energy Coulomb correction + double a = Z_ / FINE_STRUCTURE; + double c = a*a*(1.0/(1.0 + a*a) + 0.202059 + a*a*(-0.03693 + a*a*(0.00835 + + a*a*(-0.00201 + a*a*(0.00049 + a*a*(-0.00012 + a*a*0.00003)))))); + + // The analytical approximation of the DCS underestimates the cross section + // at low energies. The correction factor f compensates for this. + double q = std::sqrt(2.0/alpha); + double f = q*(-0.1774 - 12.10*a + 11.18*a*a) + + q*q*(8.523 + 73.26*a - 44.41*a*a) + + q*q*q*(-13.52 - 121.1*a + 96.41*a*a) + + q*q*q*q*(8.946 + 62.05*a - 63.41*a*a); + + // Calculate phi_1(1/2) and phi_2(1/2). The unnormalized PDF for the reduced + // energy is given by p = 2*(1/2 - e)^2*phi_1(e) + phi_2(e), where phi_1 and + // phi_2 are non-negative and maximum at e = 1/2. + double b = 2.0*r[Z_]/alpha; + double t1 = 2.0*std::log(1.0 + b*b); + double t2 = b*std::atan(1.0/b); + double t3 = b*b*(4.0 - 4.0*t2 - 3.0*std::log(1.0 + 1.0/(b*b))); + double t4 = 4.0*std::log(r[Z_]) - 4.0*c + f; + double phi1_max = 7.0/3.0 - t1 - 6.0*t2 - t3 + t4; + double phi2_max = 11.0/6.0 - t1 - 3.0*t2 + 0.5*t3 + t4; + + // To aid sampling, the unnormalized PDF can be expressed as + // p = u_1*U_1(e)*pi_1(e) + u_2*U_2(e)*pi_2(e), where pi_1 and pi_2 are + // normalized PDFs on the interval (e_min, e_max) from which values of e can + // be sampled using the inverse transform method, and + // U_1 = phi_1(e)/phi_1(1/2) and U_2 = phi_2(e)/phi_2(1/2) are valid + // rejection functions. The reduced energy can now be sampled using a + // combination of the composition and rejection methods. + double u1 = 2.0/3.0*std::pow(0.5 - 1.0/alpha, 2)*phi1_max; + double u2 = phi2_max; + double e; + while (true) { + double rn = prn(); + + // Sample the index i in (1, 2) using the point probabilities + // p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2) + int i; + if (prn() < u1/(u1 + u2)) { + i = 1; + + // Sample e from pi_1 using the inverse transform method + e = rn >= 0.5 ? + 0.5 + (0.5 - 1.0/alpha)*std::pow(2.0*rn - 1.0, 1.0/3.0) : + 0.5 - (0.5 - 1.0/alpha)*std::pow(1.0 - 2.0*rn, 1.0/3.0); + } else { + i = 2; + + // Sample e from pi_2 using the inverse transform method + e = 1.0/alpha + (0.5 - 1.0/alpha)*2.0*rn; + } + + // Calculate phi_i(e) and deliver e if rn <= U_i(e) + b = r[Z_]/(2.0*alpha*e*(1.0 - e)); + t1 = 2.0*std::log(1.0 + b*b); + t2 = b*std::atan(1.0/b); + t3 = b*b*(4.0 - 4.0*t2 - 3.0*std::log(1.0 + 1.0/(b*b))); + if (i == 1) { + double phi1 = 7.0/3.0 - t1 - 6.0*t2 - t3 + t4; + if (prn() <= phi1/phi1_max) break; + } else { + double phi2 = 11.0/6.0 - t1 - 3.0*t2 + 0.5*t3 + t4; + if (prn() <= phi2/phi2_max) break; + } + } + + // Compute the kinetic energy of the electron and the positron + *E_electron = (alpha*e - 1.0)*MASS_ELECTRON_EV; + *E_positron = (alpha*(1.0 - e) - 1.0)*MASS_ELECTRON_EV; + + // Sample the scattering angle of the electron. The cosine of the polar + // angle of the direction relative to the incident photon is sampled from + // p(mu) = C/(1 - beta*mu)^2 using the inverse transform method. + double beta = std::sqrt(*E_electron*(*E_electron + 2.0*MASS_ELECTRON_EV)) + / (*E_electron + MASS_ELECTRON_EV) ; + double rn = 2.0*prn() - 1.0; + *mu_electron = (rn + beta)/(rn*beta + 1.0); + + // Sample the scattering angle of the positron + beta = std::sqrt(*E_positron*(*E_positron + 2.0*MASS_ELECTRON_EV)) + / (*E_positron + MASS_ELECTRON_EV); + rn = 2.0*prn() - 1.0; + *mu_positron = (rn + beta)/(rn*beta + 1.0); +} + +void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particle& p) const +{ + // If no transitions, assume fluorescent photon from captured free electron + if (shell.n_transitions == 0) { + 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); + double E = shell.binding_energy; + int photon = static_cast(ParticleType::photon); + p.create_secondary(uvw.data(), E, photon, true); + return; + } + + // Sample transition + double rn = prn(); + double c = 0.0; + int i_transition; + for (i_transition = 0; i_transition < shell.n_transitions; ++i_transition) { + c += shell.transition_probability(i_transition); + if (rn < c) break; + } + + // Get primary and secondary subshell designators + int primary = shell.transition_subshells(i_transition, 0); + int secondary = shell.transition_subshells(i_transition, 1); + + // 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); + + // Get the transition energy + double E = shell.transition_energy(i_transition); + + if (secondary != 0) { + // Non-radiative transition -- Auger/Coster-Kronig effect + + // Create auger electron + int electron = static_cast(ParticleType::electron); + p.create_secondary(uvw.data(), E, electron, true); + + // Fill hole left by emitted auger electron + int i_hole = shell_map_.at(secondary); + const auto& hole = shells_[i_hole]; + this->atomic_relaxation(hole, p); + } else { + // Radiative transition -- get X-ray energy + + // Create fluorescent photon + int photon = static_cast(ParticleType::photon); + p.create_secondary(uvw.data(), E, photon, true); + } + + // Fill hole created by electron transitioning to the photoelectron hole + int i_hole = shell_map_.at(primary); + const auto& hole = shells_[i_hole]; + this->atomic_relaxation(hole, p); +} + +//============================================================================== +// Non-member functions +//============================================================================== + +std::pair klein_nishina(double alpha) +{ + double alpha_out, mu; + double beta = 1.0 + 2.0*alpha; + if (alpha < 3.0) { + // Kahn's rejection method + double t = beta/(beta + 8.0); + double x; + while (true) { + if (prn() < t) { + // Left branch of flow chart + double r = 2.0*prn(); + x = 1.0 + alpha*r; + if (prn() < 4.0/x*(1.0 - 1.0/x)) { + mu = 1 - r; + break; + } + } else { + // Right branch of flow chart + x = beta/(1.0 + 2.0*alpha*prn()); + mu = 1.0 + (1.0 - x)/alpha; + if (prn() < 0.5*(mu*mu + 1.0/x)) break; + } + } + alpha_out = alpha/x; + + } else { + // Koblinger's direct method + double gamma = 1.0 - std::pow(beta, -2); + double s = prn()*(4.0/alpha + 0.5*gamma + + (1.0 - (1.0 + beta)/(alpha*alpha))*std::log(beta)); + if (s <= 2.0/alpha) { + // For first term, x = 1 + 2ar + // Therefore, a' = a/(1 + 2ar) + alpha_out = alpha/(1.0 + 2.0*alpha*prn()); + } else if (s <= 4.0/alpha) { + // For third term, x = beta/(1 + 2ar) + // Therefore, a' = a(1 + 2ar)/beta + alpha_out = alpha*(1.0 + 2.0*alpha*prn())/beta; + } else if (s <= 4.0/alpha + 0.5*gamma) { + // For fourth term, x = 1/sqrt(1 - gamma*r) + // Therefore, a' = a*sqrt(1 - gamma*r) + alpha_out = alpha*std::sqrt(1.0 - gamma*prn()); + } else { + // For third term, x = beta^r + // Therefore, a' = a/beta^r + alpha_out = alpha/std::pow(beta, prn()); + } + + // Calculate cosine of scattering angle based on basic relation + mu = 1.0 + 1.0/alpha - 1.0/alpha_out; + } + return {alpha_out, mu}; +} + +// void thick_target_bremsstrahlung(Particle& p, double* E_lost) +// { +// if (p.material == MATERIAL_VOID) return; + +// // TODO: off-by-one +// int photon = static_cast(ParticleType::photon) - 1; +// if (p.E < settings::energy_cutoff[photon]) return; + +// // // Get bremsstrahlung data for this material and particle type +// if (p.type == static_cast(ParticleType::positron)) { +// mat => ttb(p.material) % positron; +// } else { +// mat => ttb(p.material) % electron; +// } + +// double e = std::log(p.E); +// auto n_e = data::ttb_e_grid + +// // Find the lower bounding index of the incident electron energy +// int j = lower_bound_index(data::ttb_e_grid.cbegin(), data::ttb_e_grid.cend(), e); +// if (j == n_e - 1) --j; + +// // Get the interpolation bounds +// double e_l = data::ttb_e_grid(j); +// double e_r = data::ttb_e_grid(j+1); +// double y_l = mat % yield(j); +// double y_r = mat % yield(j+1); + +// // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF +// // interpolated in log energy, which can be interpreted as the probability +// // of index j+1 +// double f = (e - e_l)/(e_r - e_l); + +// // Get the photon number yield for the given energy using linear +// // interpolation on a log-log scale +// double y = std::exp(y_l + (y_r - y_l)*f); + +// // Sample number of secondary bremsstrahlung photons +// int n = y + prn(); + +// *E_lost = 0.0; +// if (n == 0) return; + +// // Sample index of the tabulated PDF in the energy grid, j or j+1 +// double c_max; +// if (prn() <= f || j == 0) { +// int i_e = j + 1; + +// // Interpolate the maximum value of the CDF at the incoming particle +// // energy on a log-log scale +// double p_l = mat % pdf(i_e-1, i_e); +// double p_r = mat % pdf(i_e, i_e); +// double c_l = mat % cdf(i_e-1, i_e); +// double a = std::log(p_r/p_l)/(e_r - e_l) + 1.0; +// c_max = c_l + std::exp(e_l)*p_l/a*(std::exp(a*(e - e_l)) - 1.0); +// } else { +// int i_e = j; + +// // Maximum value of the CDF +// c_max = mat % cdf(i_e, i_e) +// } + +// // Sample the energies of the emitted photons +// for (int i = 0; i < n; ++i) { +// // Generate a random number r and determine the index i for which +// // cdf(i) <= r*cdf,max <= cdf(i+1) +// double c = prn()*c_max; +// int i_w = lower_bound_index(mat % cdf(:i_e,i_e), c) + +// // Sample the photon energy +// double w_l = data::ttb_e_grid(i_w); +// double w_r = data::ttb_e_grid(i_w+1); +// double p_l = mat % pdf(i_w, i_e); +// double p_r = mat % pdf(i_w+1, i_e); +// double c_l = mat % cdf(i_w, i_e); +// double a = std::log(p_r/p_l)/(w_r - w_l) + 1.0; +// double w = std::exp(w_l)*std::pow(a*(c - c_l)/(std::exp(w_l)*p_l) + 1.0, 1.0/a); + +// if (w > settings::energy_cutoff[photon]) { +// // Create secondary photon +// int photon = static_cast(ParticleType::photon); +// p.create_secondary(p.coord[0].uvw, w, photon, true); +// *E_lost += w; +// } +// } +// } + //============================================================================== // Fortran compatibility //============================================================================== diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 1a9e4edeac..1575450747 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -3,539 +3,12 @@ module photon_physics use algorithm, only: binary_search use constants use particle_header - use photon_header, only: PhotonInteraction, BremsstrahlungData, & - compton_profile_pz, ttb_e_grid, ttb + use photon_header, only: BremsstrahlungData, ttb_e_grid, ttb use random_lcg, only: prn use settings contains -!=============================================================================== -! KLEIN_NISHINA -!=============================================================================== - - subroutine klein_nishina(alpha, alpha_out, mu) - real(8), intent(in) :: alpha - real(8), intent(out) :: alpha_out - real(8), intent(out) :: mu - - real(8) :: beta ! 1 + 2a - real(8) :: t ! (1 + 2a)/(9 + 2a) - real(8) :: r, s, x - real(8) :: gamma - - beta = ONE + TWO*alpha - if (alpha < THREE) then - ! Kahn's rejection method - t = beta/(beta + 8.0_8) - do - if (prn() < t) then - ! Left branch of flow chart - r = TWO*prn() - x = ONE + alpha*r - if (prn() < FOUR/x*(ONE - ONE/x)) then - mu = 1 - r - exit - end if - else - ! Right branch of flow chart - x = beta/(ONE + TWO*alpha*prn()) - mu = ONE + (ONE - x)/alpha - if (prn() < HALF*(mu**2 + ONE/x)) exit - end if - end do - alpha_out = alpha/x - - else - ! Koblinger's direct method - gamma = ONE - beta**(-2) - s = prn()*(FOUR/alpha + HALF*gamma + & - (ONE - (ONE + beta)/alpha**2)*log(beta)) - if (s <= 2./alpha) then - ! For first term, x = 1 + 2ar - ! Therefore, a' = a/(1 + 2ar) - alpha_out = alpha/(ONE + TWO*alpha*prn()) - elseif (s <= FOUR/alpha) then - ! For third term, x = beta/(1 + 2ar) - ! Therefore, a' = a(1 + 2ar)/beta - alpha_out = alpha*(ONE + TWO*alpha*prn())/beta - elseif (s <= FOUR/alpha + HALF*gamma) then - ! For fourth term, x = 1/sqrt(1 - gamma*r) - ! Therefore, a' = a*sqrt(1 - gamma*r) - alpha_out = alpha*sqrt(ONE - gamma*prn()) - else - ! For third term, x = beta^r - ! Therefore, a' = a/beta^r - alpha_out = alpha/beta**prn() - end if - - ! Calculate cosine of scattering angle based on basic relation - mu = ONE + ONE/alpha - ONE/alpha_out - end if - - end subroutine klein_nishina - -!=============================================================================== -! COMPTON_SCATTER -!=============================================================================== - - subroutine compton_scatter(el, alpha, alpha_out, mu, i_shell, use_doppler) - type(PhotonInteraction), intent(in) :: el - real(8), intent(in) :: alpha - real(8), intent(out) :: alpha_out - real(8), intent(out) :: mu - integer, intent(out) :: i_shell - logical, intent(in), optional :: use_doppler - - real(8) :: x - real(8) :: form_factor_xmax - real(8) :: form_factor_x - real(8) :: e_out - logical :: use_doppler_ - - if (present(use_doppler)) then - use_doppler_ = use_doppler - else - use_doppler_ = .false. - end if - - form_factor_xmax = ZERO - do - ! Sample Klein-Nishina distribution for trial energy and angle - call klein_nishina(alpha, alpha_out, mu) - - ! Note that the parameter used here does not correspond exactly to the - ! momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the - ! parameter as defined by Hubbell, where the actual data comes from - x = MASS_ELECTRON_EV/PLANCK_C*alpha*sqrt(HALF*(ONE - mu)) - - ! Calculate S(x, Z) and S(x_max, Z) - form_factor_x = el % incoherent_form_factor % evaluate(x) - if (form_factor_xmax == ZERO) then - form_factor_xmax = el % incoherent_form_factor % evaluate(& - MASS_ELECTRON_EV/PLANCK_C*alpha) - end if - - ! Perform rejection on form factor - if (prn() < form_factor_x / form_factor_xmax) then - if (use_doppler_) then - call compton_doppler(el, alpha, mu, e_out, i_shell) - alpha_out = e_out/MASS_ELECTRON_EV - else - i_shell = 0 - end if - exit - end if - end do - - end subroutine compton_scatter - -!=============================================================================== -! COMPTON_DOPPLER -!=============================================================================== - - subroutine compton_doppler(el, alpha, mu, e_out, i_shell) - type(PhotonInteraction), intent(in) :: el - real(8), intent(in) :: alpha - real(8), intent(in) :: mu - real(8), intent(out) :: e_out - integer, intent(out) :: i_shell - - integer :: i - integer :: n - real(8) :: rn, m - real(8) :: c, c_l, c_max - real(8) :: pz_l, pz_r, pz, pz_max - real(8) :: p_l, p_r - real(8) :: e, e_b - real(8) :: e_out1, e_out2 - real(8) :: a, b, quad - real(8) :: f - real(8) :: momentum_sq - - n = size(compton_profile_pz) - - do - ! Sample electron shell - rn = prn() - c = ZERO - do i_shell = 1, size(el % electron_pdf) - c = c + el % electron_pdf(i_shell) - if (rn < c) exit - end do - - ! Determine binding energy of shell - e_b = el % binding_energy(i_shell) - - ! Determine p_z,max - e = alpha*MASS_ELECTRON_EV - if (e < e_b) then - e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV - exit - end if - - pz_max = -FINE_STRUCTURE*(e_b - (e - e_b)*alpha*(ONE - mu)) / & - sqrt(TWO*e*(e - e_b)*(ONE - mu) + e_b**2) - if (pz_max < ZERO) then - e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV - exit - end if - - ! Determine profile cdf value corresponding to p_z,max - if (pz_max > compton_profile_pz(n)) then - c_max = el % profile_cdf(n, i_shell) - else - i = binary_search(compton_profile_pz, n, pz_max) - pz_l = compton_profile_pz(i) - pz_r = compton_profile_pz(i + 1) - p_l = el % profile_pdf(i, i_shell) - p_r = el % profile_pdf(i + 1, i_shell) - c_l = el % profile_cdf(i, i_shell) - if (pz_l == pz_r) then - c_max = c_l - elseif (p_l == p_r) then - c_max = c_l + (pz_max - pz_l)*p_l - else - m = (p_l - p_r)/(pz_l - pz_r) - c_max = c_l + ((m*(pz_max - pz_l) + p_l)**2 - p_l**2)/(TWO*m) - end if - end if - - ! Sample value on bounded cdf - c = prn()*c_max - - ! Determine pz corresponding to sampled cdf value - i = binary_search(el % profile_cdf(:, i_shell), n, c) - pz_l = compton_profile_pz(i) - pz_r = compton_profile_pz(i + 1) - p_l = el % profile_pdf(i, i_shell) - p_r = el % profile_pdf(i + 1, i_shell) - c_l = el % profile_cdf(i, i_shell) - if (pz_l == pz_r) then - pz = pz_l - elseif (p_l == p_r) then - pz = pz_l + (c - c_l)/p_l - else - m = (p_l - p_r)/(pz_l - pz_r) - pz = pz_l + (sqrt(p_l**2 + TWO*m*(c - c_l)) - p_l)/m - end if - - ! Determine outgoing photon energy corresponding to electron momentum - ! (solve Eq. 39 in LA-UR-04-0487 for E') - momentum_sq = (pz/FINE_STRUCTURE)**2 - f = ONE + alpha*(ONE - mu) - a = momentum_sq - f*f - b = TWO*e*(f - momentum_sq*mu) - c = e**2*(momentum_sq - ONE) - - quad = b**2 - FOUR*a*c - if (quad < 0) then - e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV - exit - end if - quad = sqrt(quad) - e_out1 = -(b + quad)/(TWO*a) - e_out2 = -(b - quad)/(TWO*a) - - ! Determine solution to quadratic equation that is positive - if (e_out1 > ZERO) then - if (e_out2 > ZERO) then - ! If both are positive, pick one at random - if (prn() < HALF) then - e_out = e_out1 - else - e_out = e_out2 - end if - else - e_out = e_out1 - end if - else - if (e_out2 > ZERO) then - e_out = e_out2 - else - ! No positive solution -- resample - cycle - end if - end if - if (e_out < e - e_b) exit - end do - - end subroutine compton_doppler - -!=============================================================================== -! RAYLEIGH_SCATTER -!=============================================================================== - - subroutine rayleigh_scatter(el, alpha, mu) - type(PhotonInteraction), intent(in) :: el - real(8), intent(in) :: alpha - real(8), intent(out) :: mu - - integer :: i - real(8) :: F - real(8) :: F_max - real(8) :: x2 - real(8) :: x2_max - real(8) :: r - - do - ! Determine maximum value of x^2 - x2_max = (MASS_ELECTRON_EV/PLANCK_C*alpha)**2 - - ! Determine F(x^2_max, Z) - F_max = el % coherent_int_form_factor % evaluate(x2_max) - - ! Sample cumulative distribution - F = prn()*F_max - - ! Determine x^2 corresponding to F - i = binary_search(el%coherent_int_form_factor%y, & - size(el%coherent_int_form_factor%y), F) - r = (F - el%coherent_int_form_factor%y(i)) / & - (el%coherent_int_form_factor%y(i+1) - el%coherent_int_form_factor%y(i)) - x2 = el%coherent_int_form_factor%x(i) + r*(el%coherent_int_form_factor%x(i+1) - & - el%coherent_int_form_factor%x(i)) - - ! Calculate mu - mu = ONE - TWO*x2/x2_max - - if (prn() < HALF*(ONE + mu**2)) exit - end do - - end subroutine rayleigh_scatter - -!=============================================================================== -! ATOMIC_RELAXATION -!=============================================================================== - - recursive subroutine atomic_relaxation(p, elm, i_shell) - type(Particle), intent(inout) :: p - type(PhotonInteraction), intent(in) :: elm - integer, intent(in) :: i_shell - - integer :: i_hole - integer :: i_transition - integer :: primary - integer :: secondary - real(8) :: c - real(8) :: rn - real(8) :: E - real(8) :: mu - real(8) :: phi - real(8) :: uvw(3) - - ! If no transitions, assume fluorescent photon from captured free electron - if (elm % shells(i_shell) % n_transitions == 0) then - 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) - E = elm % shells(i_shell) % binding_energy - call particle_create_secondary(p, uvw, E, PHOTON, run_ce=.true._C_BOOL) - return - end if - - ! Sample transition - rn = prn() - c = ZERO - do i_transition = 1, elm % shells(i_shell) % n_transitions - c = c + elm % shells(i_shell) % & - transition_probability(i_transition) - if (rn < c) exit - end do - - ! Get primary and secondary subshell designators - primary = elm % shells(i_shell) % transition_subshells(1, i_transition) - secondary = elm % shells(i_shell) % transition_subshells(2, i_transition) - - ! 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) - - ! Get the transition energy - E = elm % shells(i_shell) % transition_energy(i_transition) - - if (secondary /= 0) then - ! Non-radiative transition -- Auger/Coster-Kronig effect - - ! Create auger electron - call particle_create_secondary(p, uvw, E, ELECTRON, run_ce=.true._C_BOOL) - - ! Fill hole left by emitted auger electron - i_hole = elm % shell_dict % get(secondary) - call atomic_relaxation(p, elm, i_hole) - else - ! Radiative transition -- get X-ray energy - - ! Create fluorescent photon - call particle_create_secondary(p, uvw, E, PHOTON, run_ce=.true._C_BOOL) - - end if - - ! Fill hole created by electron transitioning to the photoelectron hole - i_hole = elm % shell_dict % get(primary) - call atomic_relaxation(p, elm, i_hole) - - end subroutine atomic_relaxation - -!=============================================================================== -! PAIR_PRODUCTION samples the kinetic energy and direction of the electron and -! positron created when a photon is absorbed near an atomic nucleus. The -! simulation procedure follows the semiempirical model outlined in F. Salvat, J. -! M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for Monte -! Carlo Simulation of Electron and Photon Transport," OECD-NEA, -! Issy-les-Moulineaux, France (2011). -!=============================================================================== - - subroutine pair_production(elm, alpha, E_electron, E_positron, mu_electron, & - mu_positron) - type(PhotonInteraction), intent(in) :: elm - real(8), intent(in) :: alpha - real(8), intent(out) :: E_electron - real(8), intent(out) :: E_positron - real(8), intent(out) :: mu_electron - real(8), intent(out) :: mu_positron - - integer :: i - real(8) :: f - real(8) :: c - real(8) :: a - real(8) :: b - real(8) :: q - real(8) :: rn - real(8) :: beta - real(8) :: e, e_min, e_max - real(8) :: t1, t2, t3, t4 - real(8) :: u1, u2 - real(8) :: phi1, phi2 - real(8) :: phi1_max, phi2_max - real(8), parameter :: r(99) = (/ & - 122.81_8, 73.167_8, 69.228_8, 67.301_8, 64.696_8, 61.228_8, & - 57.524_8, 54.033_8, 50.787_8, 47.851_8, 46.373_8, 45.401_8, & - 44.503_8, 43.815_8, 43.074_8, 42.321_8, 41.586_8, 40.953_8, & - 40.524_8, 40.256_8, 39.756_8, 39.144_8, 38.462_8, 37.778_8, & - 37.174_8, 36.663_8, 35.986_8, 35.317_8, 34.688_8, 34.197_8, & - 33.786_8, 33.422_8, 33.068_8, 32.740_8, 32.438_8, 32.143_8, & - 31.884_8, 31.622_8, 31.438_8, 31.142_8, 30.950_8, 30.758_8, & - 30.561_8, 30.285_8, 30.097_8, 29.832_8, 29.581_8, 29.411_8, & - 29.247_8, 29.085_8, 28.930_8, 28.721_8, 28.580_8, 28.442_8, & - 28.312_8, 28.139_8, 27.973_8, 27.819_8, 27.675_8, 27.496_8, & - 27.285_8, 27.093_8, 26.911_8, 26.705_8, 26.516_8, 26.304_8, & - 26.108_8, 25.929_8, 25.730_8, 25.577_8, 25.403_8, 25.245_8, & - 25.100_8, 24.941_8, 24.790_8, 24.655_8, 24.506_8, 24.391_8, & - 24.262_8, 24.145_8, 24.039_8, 23.922_8, 23.813_8, 23.712_8, & - 23.621_8, 23.523_8, 23.430_8, 23.331_8, 23.238_8, 23.139_8, & - 23.048_8, 22.967_8, 22.833_8, 22.694_8, 22.624_8, 22.545_8, & - 22.446_8, 22.358_8, 22.264_8 /) - - ! The reduced screening radius r is the ratio of the screening radius to - ! the Compton wavelength of the electron, where the screening radius is - ! obtained under the assumption that the Coulomb field of the nucleus is - ! exponentially screened by atomic electrons. This allows us to use a - ! simplified atomic form factor and analytical approximations of the - ! screening functions in the pair production DCS instead of computing the - ! screening functions numerically. The reduced screening radii above for - ! Z = 1-99 come from F. Salvat, J. M. Fernández-Varea, and J. Sempau, - ! "PENELOPE-2011: A Code System for Monte Carlo Simulation of Electron and - ! Photon Transport," OECD-NEA, Issy-les-Moulineaux, France (2011). - - ! Compute the minimum and maximum values of the electron reduced energy, - ! i.e. the fraction of the photon energy that is given to the electron - e_min = ONE/alpha - e_max = ONE - ONE/alpha - - ! Compute the high-energy Coulomb correction - a = elm % Z / FINE_STRUCTURE - c = a**2*(ONE/(ONE + a**2) + 0.202059_8 - 0.03693_8*a**2 + 0.00835_8*a**4 & - - 0.00201_8*a**6 + 0.00049_8*a**8 - 0.00012_8*a**10 + 0.00003_8*a**12) - - ! The analytical approximation of the DCS underestimates the cross section - ! at low energies. The correction factor f compensates for this. - q = sqrt(TWO/alpha) - f = q*(-0.1774_8 - 12.10_8*a + 11.18_8*a**2) & - + q**2*(8.523_8 + 73.26_8*a - 44.41_8*a**2) & - + q**3*(-13.52_8 - 121.1_8*a + 96.41_8*a**2) & - + q**4*(8.946_8 + 62.05_8*a - 63.41_8*a**2) - - ! Calculate phi_1(1/2) and phi_2(1/2). The unnormalized PDF for the reduced - ! energy is given by p = 2*(1/2 - e)^2*phi_1(e) + phi_2(e), where phi_1 and - ! phi_2 are non-negative and maximum at e = 1/2. - b = TWO*r(elm % Z)/alpha - t1 = TWO*log(ONE + b**2) - t2 = b*atan(ONE/b) - t3 = b**2*(FOUR - FOUR*t2 - THREE*log(ONE + ONE/b**2)) - t4 = FOUR*log(r(elm % Z)) - FOUR*c + f - phi1_max = 7.0_8/THREE - t1 - 6.0_8*t2 - t3 + t4 - phi2_max = 11.0_8/6.0_8 - t1 - THREE*t2 + HALF*t3 + t4 - - ! To aid sampling, the unnormalized PDF can be expressed as - ! p = u_1*U_1(e)*pi_1(e) + u_2*U_2(e)*pi_2(e), where pi_1 and pi_2 are - ! normalized PDFs on the interval (e_min, e_max) from which values of e can - ! be sampled using the inverse transform method, and - ! U_1 = phi_1(e)/phi_1(1/2) and U_2 = phi_2(e)/phi_2(1/2) are valid - ! rejection functions. The reduced energy can now be sampled using a - ! combination of the composition and rejection methods. - u1 = TWO/THREE*(HALF - ONE/alpha)**2*phi1_max - u2 = phi2_max - do - rn = prn() - - ! Sample the index i in (1, 2) using the point probabilities - ! p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2) - if (prn() < u1/(u1 + u2)) then - i = 1 - - ! Sample e from pi_1 using the inverse transform method - if (rn >= HALF) then - e = HALF + (HALF - ONE/alpha)*(TWO*rn - ONE)**(ONE/THREE) - else - e = HALF - (HALF - ONE/alpha)*(ONE - TWO*rn)**(ONE/THREE) - end if - else - i = 2 - - ! Sample e from pi_2 using the inverse transform method - e = ONE/alpha + (HALF - ONE/alpha)*TWO*rn - end if - - ! Calculate phi_i(e) and deliver e if rn <= U_i(e) - b = r(elm % Z)/(TWO*alpha*e*(ONE - e)) - t1 = TWO*log(ONE + b**2) - t2 = b*atan(ONE/b) - t3 = b**2*(FOUR - FOUR*t2 - THREE*log(ONE + ONE/b**2)) - if (i == 1) then - phi1 = 7.0_8/THREE - t1 - 6.0_8*t2 - t3 + t4 - if (prn() <= phi1/phi1_max) exit - else - phi2 = 11.0_8/6.0_8 - t1 - THREE*t2 + HALF*t3 + t4 - if (prn() <= phi2/phi2_max) exit - end if - end do - - ! Compute the kinetic energy of the electron and the positron - E_electron = (alpha*e - ONE)*MASS_ELECTRON_EV - E_positron = (alpha*(ONE - e) - ONE)*MASS_ELECTRON_EV - - ! Sample the scattering angle of the electron. The cosine of the polar - ! angle of the direction relative to the incident photon is sampled from - ! p(mu) = C/(1 - beta*mu)^2 using the inverse transform method. - beta = sqrt(E_electron*(E_electron + TWO*MASS_ELECTRON_EV)) & - / (E_electron + MASS_ELECTRON_EV) - rn = TWO*prn() - ONE - mu_electron = (rn + beta)/(rn*beta + ONE) - - ! Sample the scattering angle of the positron - beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON_EV)) & - / (E_positron + MASS_ELECTRON_EV) - rn = TWO*prn() - ONE - mu_positron = (rn + beta)/(rn*beta + ONE) - - end subroutine pair_production - !=============================================================================== ! THICK_TARGET_BREMSSTRAHLUNG !=============================================================================== diff --git a/src/physics.F90 b/src/physics.F90 index d77fb04dff..c6749a6a1f 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,21 +1,16 @@ module physics use constants - use error, only: fatal_error, warning, write_message + use error, only: fatal_error use material_header, only: Material, materials use math use message_passing use nuclide_header use particle_header use photon_header - use photon_physics, only: rayleigh_scatter, compton_scatter, & - atomic_relaxation, pair_production, & - thick_target_bremsstrahlung - use physics_common use random_lcg, only: prn use settings use simulation_header - use tally_header implicit none @@ -28,190 +23,13 @@ module physics contains -!=============================================================================== -! 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) bind(C) - type(Particle), intent(inout) :: p - - integer :: i_shell ! index in subshells - integer :: i_grid ! index on energy grid - integer :: i_element ! index in nuclides array - integer :: i_start ! threshold index - real(8) :: prob ! cumulative probability - real(8) :: cutoff ! sampled total cross section - real(8) :: f ! interpolation factor - real(8) :: xs ! photoionization cross section - real(8) :: r ! random number - real(8) :: prob_after - real(8) :: alpha ! photon energy divided by electron rest mass - real(8) :: alpha_out ! outgoing photon energy over electron rest mass - real(8) :: mu ! scattering cosine - real(8) :: mu_electron ! electron scattering cosine - real(8) :: mu_positron ! positron scattering cosine - real(8) :: phi ! azimuthal angle - real(8) :: uvw(3) ! new direction - real(8) :: rel_vel ! relative velocity of electron - real(8) :: e_b ! binding energy of electron - real(8) :: E_electron ! electron energy - real(8) :: E_positron ! positron energy - - ! 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)) then - p % E = ZERO - p % alive = .false. - return - end if - - ! 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 = ZERO - 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) then - call rayleigh_scatter(elm, alpha, mu) - p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) - p % event_MT = COHERENT - return - end if - - ! Incoherent (Compton) scattering - prob = prob + micro_photon_xs(i_element) % incoherent - if (prob > cutoff) then - call compton_scatter(elm, alpha, alpha_out, mu, i_shell, .true.) - - ! Determine binding energy of shell. The binding energy is zero if - ! doppler broadening is not used. - if (i_shell == 0) then - e_b = ZERO - else - e_b = elm % binding_energy(i_shell) - end if - - ! Create Compton electron - E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b - mu_electron = (alpha - alpha_out*mu) & - / sqrt(alpha**2 + alpha_out**2 - TWO*alpha*alpha_out*mu) - phi = TWO*PI*prn() - uvw = rotate_angle(p % coord(1) % uvw, mu_electron, phi) - call particle_create_secondary(p, uvw, E_electron, ELECTRON, .true._C_BOOL) - - ! 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) then - call atomic_relaxation(p, elm, i_shell) - end if - - 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 - end if - - ! Photoelectric effect - prob_after = prob + micro_photon_xs(i_element) % photoelectric - if (prob_after > cutoff) then - 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 = 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) then - 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 * (ONE - r) * r >= prn()) then - rel_vel = sqrt(E_electron * (E_electron + TWO * MASS_ELECTRON_EV))& - / (E_electron + MASS_ELECTRON_EV) - mu = (TWO * r + rel_vel - ONE) / & - (TWO * rel_vel * r - rel_vel + ONE) - exit SAMPLE_MU - end if - end do SAMPLE_MU - - phi = TWO*PI*prn() - uvw(1) = mu - uvw(2) = sqrt(ONE - mu*mu)*cos(phi) - uvw(3) = sqrt(ONE - mu*mu)*sin(phi) - - ! Create secondary electron - call particle_create_secondary(p, uvw, E_electron, ELECTRON, & - run_CE=.true._C_BOOL) - - ! Allow electrons to fill orbital and produce auger electrons - ! and fluorescent photons - call atomic_relaxation(p, elm, i_shell) - p % event_MT = 533 + elm % shells(i_shell) % index_subshell - p % alive = .false. - p % E = ZERO - - return - end if - end do - end if - prob = prob_after - - ! Pair production - prob = prob + micro_photon_xs(i_element) % pair_production - if (prob > cutoff) then - call pair_production(elm, alpha, E_electron, E_positron, mu_electron, & - mu_positron) - - ! Create secondary electron - uvw = rotate_angle(p % coord(1) % uvw, mu_electron) - call particle_create_secondary(p, uvw, E_electron, ELECTRON, .true._C_BOOL) - - ! Create secondary positron - uvw = rotate_angle(p % coord(1) % uvw, mu_positron) - call particle_create_secondary(p, uvw, E_positron, POSITRON, .true._C_BOOL) - - p % event_MT = PAIR_PROD - p % alive = .false. - p % E = ZERO - end if - - end associate - - end subroutine sample_photon_reaction - !=============================================================================== ! SAMPLE_ELEMENT !=============================================================================== - function sample_element(p) result(i_element) + function sample_element(p) result(i_element) bind(C) type(Particle), intent(in) :: p - integer :: i_element + integer(C_INT) :: i_element integer :: i real(8) :: prob diff --git a/src/physics.cpp b/src/physics.cpp index 8aec309e83..a4bd670660 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -214,154 +214,168 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ } } -// 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 + int photon = static_cast(ParticleType::photon) - 1; + if (p->E < settings::energy_cutoff[photon]) { + p->E = 0.0; + p->alive = false; + return; + } -// 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 + int i_element = sample_element(p); + p->event_nuclide = i_element; + // TODO: off-by-one + const auto& micro {simulation::micro_photon_xs[i_element - 1]}; + const auto& element {data::elements[i_element - 1]}; -// // Sample element within material -// i_element = sample_element(p) -// p->event_nuclide = i_element + // Calculate photon energy over electron rest mass equivalent + double alpha = p->E/MASS_ELECTRON_EV; -// // 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 + double prob = 0.0; + double cutoff = prn() * micro.total; -// // 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 + // Coherent (Rayleigh) scattering + prob += micro.coherent; + if (prob > cutoff) { + double mu = element.rayleigh_scatter(alpha); + rotate_angle_c(p->coord[0].uvw, mu, nullptr); + p->event_MT = COHERENT; + return; + } -// 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 += micro.incoherent; + if (prob > cutoff) { + double alpha_out, mu; + int i_shell; + element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell); -// // 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. + double e_b; + if (i_shell == -1) { + e_b = 0.0; + } else { + e_b = element.binding_energy_[i_shell]; + } -// // 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 + double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b; + double mu_electron = (alpha - alpha_out*mu) + / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu); + double phi = 2.0*PI*prn(); + double uvw[3]; + std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + rotate_angle_c(uvw, mu_electron, &phi); + int electron = static_cast(ParticleType::electron); + p->create_secondary(uvw, E_electron, electron, true); -// // 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) { + const auto& shell = element.shells_[i_shell]; + element.atomic_relaxation(shell, *p); + } -// // 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 += PI; + p->E = alpha_out*MASS_ELECTRON_EV; + rotate_angle_c(p->coord[0].uvw, mu, &phi); + p->event_MT = INCOHERENT; + return; + } -// 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 + double prob_after = prob + micro.photoelectric; + if (prob_after > cutoff) { + for (const auto& shell : element.shells_) { + // Get grid index and interpolation factor + // TODO: off-by-one + int i_grid = micro.index_grid - 1; + double f = micro.interp_factor; -// // 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 + int i_start = shell.threshold; + if (i_grid <= i_start) continue; -// // Check threshold of reaction -// i_start = elm % shells(i_shell) % threshold -// if (i_grid <= i_start) cycle + // Evaluation subshell photoionization cross section + double xs = std::exp(shell.cross_section(i_grid - i_start) + + f*(shell.cross_section(i_grid + 1 - i_start) - + shell.cross_section(i_grid - i_start))); -// // 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 += xs; + if (prob > cutoff) { + double E_electron = p->E - shell.binding_energy; -// 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 + double mu; + while (true) { + double r = prn(); + if (4.0*(1.0 - r)*r >= prn()) { + double 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); + break; + } + } -// // 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 + 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); -// 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 + int electron = static_cast(ParticleType::electron); + p->create_secondary(uvw.data(), E_electron, electron, true); -// // 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 + element.atomic_relaxation(shell, *p); + p->event_MT = 533 + shell.index_subshell; + p->alive = false; + p->E = 0.0; + return; + } + } + } + prob = prob_after; -// // 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 + // Pair production + prob += micro.pair_production; + if (prob > cutoff) { + double E_electron, E_positron; + double mu_electron, mu_positron; + element.pair_production(alpha, &E_electron, &E_positron, + &mu_electron, &mu_positron); -// return -// } -// end do -// } -// prob = prob_after + // Create secondary electron + double uvw[3]; + std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + rotate_angle_c(uvw, mu_electron, nullptr); + int electron = static_cast(ParticleType::electron); + p->create_secondary(uvw, E_electron, electron, true); -// // 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 positron + std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + rotate_angle_c(uvw, mu_positron, nullptr); + int positron = static_cast(ParticleType::positron); + p->create_secondary(uvw, E_positron, positron, true); -// // 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 -// } + p->event_MT = PAIR_PROD; + p->alive = false; + p->E = 0.0; + } +} void sample_electron_reaction(Particle* p) { @@ -465,36 +479,31 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat) } } -// TODO: Finish converting photon physics functions - -// void sample_element(Particle* p) +// int sample_element(Particle* p) // { -// associate (mat => materials(p->material)) -// // Sample cumulative distribution function -// cutoff = prn() * simulation::material_xs.total +// // Sample cumulative distribution function +// double cutoff = prn() * simulation::material_xs.total; -// i = 0 -// prob = 0.0 -// do while (prob < cutoff) -// i = i + 1 +// int i = 0; +// double prob = 0.0; +// while (prob < cutoff) { +// // Check to make sure that a nuclide was sampled +// if (i > mat % n_nuclides) { +// p->write_restart(); +// fatal_error("Did not sample any element during collision."); +// } -// // 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 +// int i_element = mat % element(i); +// double atom_density = mat % atom_density(i); -// // Find atom density -// i_element = mat % element(i) -// atom_density = mat % atom_density(i) +// // Determine microscopic cross section +// double sigma = atom_density * simulation::micro_photon_xs[i_element].total; -// // 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 +// // Increment probability to compare to cutoff +// prob += sigma; +// ++i; +// } // } Reaction* sample_fission(int i_nuclide, double E) From 39cf5a8ce7891fcd59e05ac006aee1d571eac0bc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 8 Jan 2019 08:36:44 -0600 Subject: [PATCH 04/14] Move calculate_xs for photon data to C++ --- include/openmc/photon.h | 13 +- src/material_header.F90 | 11 +- src/photon.cpp | 70 ++++++++- src/photon_header.F90 | 317 +--------------------------------------- src/physics.cpp | 3 +- 5 files changed, 89 insertions(+), 325 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index fce6a3626c..81675585ff 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -39,9 +39,11 @@ public: class PhotonInteraction { public: // Constructors - PhotonInteraction(hid_t group); + PhotonInteraction(hid_t group, int i_element); // Methods + void calculate_xs(double E) const; + void compton_scatter(double alpha, bool doppler, double* alpha_out, double* mu, int* i_shell) const; @@ -53,8 +55,9 @@ public: void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const; // Data members - std::string name_; //! Name of element, e.g. "Zr" - int Z_; //! Atomic number + std::string name_; //!< Name of element, e.g. "Zr" + int Z_; //!< Atomic number + int i_element_; //!< Index in global elements vector // Microscopic cross sections xt::xtensor energy_; @@ -72,8 +75,8 @@ public: Tabulated1D coherent_anomalous_imag_; // Photoionization and atomic relaxation data - std::unordered_map shell_map_; // Given a shell designator, e.g. 3, this - // dictionary gives an index in shells_ + std::unordered_map shell_map_; //!< Given a shell designator, e.g. 3, this + //!< dictionary gives an index in shells_ std::vector shells_; // Compton profile data diff --git a/src/material_header.F90 b/src/material_header.F90 index 0a2548df39..9341732419 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -487,6 +487,14 @@ contains integer :: i_element ! index into elements array real(8) :: atom_density ! atom density of a nuclide + interface + subroutine photon_calculate_xs(i_element, E) bind(C) + import C_INT, C_DOUBLE + integer(C_INT), value :: i_element + real(C_DOUBLE), value :: E + end subroutine + end interface + material_xs % coherent = ZERO material_xs % incoherent = ZERO material_xs % photoelectric = ZERO @@ -502,8 +510,7 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_photon_xs(i_element) % last_E) then - call elements(i_element) % calculate_xs(& - p % E, micro_photon_xs(i_element)) + call photon_calculate_xs(i_element, p % E) end if ! ======================================================================== diff --git a/src/photon.cpp b/src/photon.cpp index 56ff2042dd..7432a2dcd6 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -39,7 +39,8 @@ ElementMicroXS* micro_photon_xs; // PhotonInteraction implementation //============================================================================== -PhotonInteraction::PhotonInteraction(hid_t group) +PhotonInteraction::PhotonInteraction(hid_t group, int i_element) + : i_element_{i_element} { // Get name of nuclide from group, removing leading '/' name_ = object_name(group).substr(1); @@ -430,6 +431,66 @@ void PhotonInteraction::compton_doppler(double alpha, double mu, *i_shell = shell; } +void PhotonInteraction::calculate_xs(double E) const +{ + // Perform binary search on the element energy grid in order to determine + // which points to interpolate between + int n_grid = energy_.size(); + double log_E = std::log(E); + int i_grid; + if (log_E <= energy_[0]) { + i_grid = 0; + } else if (log_E > energy_(n_grid - 1)) { + i_grid = n_grid - 2; + } else { + // We use upper_bound_index here because sometimes photons are created with + // energies that exactly match a grid point + i_grid = upper_bound_index(energy_.cbegin(), energy_.cend(), log_E); + } + + // check for case where two energy points are the same + if (energy_(i_grid) == energy_(i_grid+1)) ++i_grid; + + // calculate interpolation factor + double f = (log_E - energy_(i_grid)) / (energy_(i_grid+1) - energy_(i_grid)); + + auto& xs {simulation::micro_photon_xs[i_element_]}; + xs.index_grid = i_grid; + xs.interp_factor = f; + + // Calculate microscopic coherent cross section + xs.coherent = std::exp(coherent_(i_grid) + + f*(coherent_(i_grid+1) - coherent_(i_grid))); + + // Calculate microscopic incoherent cross section + xs.incoherent = std::exp(incoherent_(i_grid) + + f*(incoherent_(i_grid+1) - incoherent_(i_grid))); + + // Calculate microscopic photoelectric cross section + xs.photoelectric = 0.0; + for (const auto& shell : shells_) { + // Check threshold of reaction + int i_start = shell.threshold; + if (i_grid < i_start) continue; + + // Evaluation subshell photoionization cross section + xs.photoelectric += + std::exp(shell.cross_section(i_grid-i_start) + + f*(shell.cross_section(i_grid+1-i_start) - + shell.cross_section(i_grid-i_start))); + } + + // Calculate microscopic pair production cross section + xs.pair_production = std::exp( + pair_production_total_(i_grid) + f*( + pair_production_total_(i_grid+1) - + pair_production_total_(i_grid))); + + // Calculate microscopic total cross section + xs.total = xs.coherent + xs.incoherent + xs.photoelectric + xs.pair_production; + xs.last_E = E; +} + double PhotonInteraction::rayleigh_scatter(double alpha) const { double mu; @@ -800,7 +861,12 @@ std::pair klein_nishina(double alpha) extern "C" void photon_from_hdf5_c(hid_t group) { - data::elements.emplace_back(group); + data::elements.emplace_back(group, data::elements.size()); +} + +extern "C" void photon_calculate_xs(int i_element, double E) +{ + data::elements[i_element - 1].calculate_xs(E); } } // namespace openmc diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 4be73831d0..5a887d259a 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -4,59 +4,20 @@ module photon_header use algorithm, only: binary_search use constants - use dict_header, only: DictIntInt, DictCharInt - use endf_header, only: Tabulated1D + use dict_header, only: DictCharInt use hdf5_interface - use nuclide_header, only: nuclides use settings real(8), allocatable :: compton_profile_pz(:) real(8), allocatable :: ttb_e_grid(:) ! energy T of incident electron real(8), allocatable :: ttb_k_grid(:) ! reduced energy W/T of emitted photon - type ElectronSubshell - integer :: index_subshell ! index in SUBSHELLS - integer :: threshold - real(8) :: n_electrons - real(8) :: binding_energy - real(8), allocatable :: cross_section(:) - - ! Transition data - integer :: n_transitions - integer, allocatable :: transition_subshells(:,:) - real(8), allocatable :: transition_energy(:) - real(8), allocatable :: transition_probability(:) - end type ElectronSubshell - type PhotonInteraction character(3) :: name ! atomic symbol, e.g. 'Zr' integer :: Z ! atomic number ! Microscopic cross sections real(8), allocatable :: energy(:) - real(8), allocatable :: coherent(:) - real(8), allocatable :: incoherent(:) - real(8), allocatable :: photoelectric_total(:) - real(8), allocatable :: pair_production_total(:) - real(8), allocatable :: pair_production_electron(:) - real(8), allocatable :: pair_production_nuclear(:) - - ! Form factors - type(Tabulated1D) :: incoherent_form_factor - type(Tabulated1D) :: coherent_int_form_factor - type(Tabulated1D) :: coherent_anomalous_real - type(Tabulated1D) :: coherent_anomalous_imag - - ! Photoionization and atomic relaxation data - type(DictIntInt) :: shell_dict ! Given a shell designator, e.g. 3, this - ! dictionary gives an index in shells(:) - type(ElectronSubshell), allocatable :: shells(:) - - ! Compton profile data - real(8), allocatable :: profile_pdf(:,:) - real(8), allocatable :: profile_cdf(:,:) - real(8), allocatable :: binding_energy(:) - real(8), allocatable :: electron_pdf(:) ! Stopping power data real(8) :: I ! mean excitation energy @@ -68,7 +29,6 @@ module photon_header contains procedure :: from_hdf5 => photon_from_hdf5 - procedure :: calculate_xs => photon_calculate_xs end type PhotonInteraction type BremsstrahlungData @@ -114,22 +74,16 @@ contains class(PhotonInteraction), intent(inout) :: this integer(HID_T), intent(in) :: group_id - integer :: i, j - integer(HID_T) :: rgroup, tgroup + integer :: i + integer(HID_T) :: rgroup integer(HID_T) :: dset_id integer(HSIZE_T) :: dims(1), dims2(2) integer :: n_energy - integer :: n_shell - integer :: n_profile - integer :: n_transition integer :: n_k integer :: n_e - character(3), allocatable :: designators(:) - real(8) :: c real(8) :: f real(8) :: y real(8), allocatable :: electron_energy(:) - real(8), allocatable :: matrix(:,:) real(8), allocatable :: dcs(:,:) interface @@ -159,171 +113,6 @@ contains call read_dataset(this % energy, dset_id) call close_dataset(dset_id) - ! Allocate arrays - allocate(this % coherent(n_energy)) - allocate(this % incoherent(n_energy)) - allocate(this % pair_production_total(n_energy)) - allocate(this % pair_production_nuclear(n_energy)) - allocate(this % pair_production_electron(n_energy)) - allocate(this % photoelectric_total(n_energy)) - - ! Read coherent scattering - rgroup = open_group(group_id, 'coherent') - call read_dataset(this % coherent, rgroup, 'xs') - - dset_id = open_dataset(rgroup, 'integrated_scattering_factor') - call this % coherent_int_form_factor % from_hdf5(dset_id) - call close_dataset(dset_id) - - if (object_exists(group_id, 'anomalous_real')) then - dset_id = open_dataset(rgroup, 'anomalous_real') - call this % coherent_anomalous_real % from_hdf5(dset_id) - call close_dataset(dset_id) - end if - - if (object_exists(group_id, 'anomalous_imag')) then - dset_id = open_dataset(rgroup, 'anomalous_imag') - call this % coherent_anomalous_imag % from_hdf5(dset_id) - call close_dataset(dset_id) - call close_group(rgroup) - end if - - ! Read incoherent scattering - rgroup = open_group(group_id, 'incoherent') - call read_dataset(this % incoherent, rgroup, 'xs') - dset_id = open_dataset(rgroup, 'scattering_factor') - call this % incoherent_form_factor % from_hdf5(dset_id) - call close_dataset(dset_id) - call close_group(rgroup) - - ! Read pair production - rgroup = open_group(group_id, 'pair_production_electron') - call read_dataset(this % pair_production_electron, rgroup, 'xs') - call close_group(rgroup) - - ! Read pair production - if (object_exists(group_id, 'pair_production_nuclear')) then - rgroup = open_group(group_id, 'pair_production_nuclear') - call read_dataset(this % pair_production_nuclear, rgroup, 'xs') - call close_group(rgroup) - else - this % pair_production_nuclear(:) = ZERO - end if - - ! Read photoelectric - rgroup = open_group(group_id, 'photoelectric') - call read_dataset(this % photoelectric_total, rgroup, 'xs') - call close_group(rgroup) - - ! Read subshell photoionization cross section and atomic relaxation data - rgroup = open_group(group_id, 'subshells') - call read_attribute(designators, rgroup, 'designators') - n_shell = size(designators) - allocate(this % shells(n_shell)) - do i = 1, n_shell - ! Create mapping from designator to index - do j = 1, size(SUBSHELLS) - if (designators(i) == SUBSHELLS(j)) then - call this % shell_dict % set(j, i) - this % shells(i) % index_subshell = j - exit - end if - end do - - ! Read binding energy and number of electrons - tgroup = open_group(rgroup, trim(designators(i))) - call read_attribute(this % shells(i) % binding_energy, tgroup, & - 'binding_energy') - call read_attribute(this % shells(i) % n_electrons, tgroup, & - 'num_electrons') - - ! Read subshell cross section - dset_id = open_dataset(tgroup, 'xs') - call read_attribute(j, dset_id, 'threshold_idx') - this % shells(i) % threshold = j - allocate(this % shells(i) % cross_section(n_energy - j)) - call read_dataset(this % shells(i) % cross_section, dset_id) - call close_dataset(dset_id) - where (this % shells(i) % cross_section > ZERO) - this % shells(i) % cross_section = log(this % shells(i) % cross_section) - elsewhere - this % shells(i) % cross_section = -500.0_8 - end where - - if (object_exists(tgroup, 'transitions')) then - dset_id = open_dataset(tgroup, 'transitions') - call get_shape(dset_id, dims2) - n_transition = int(dims2(2), 4) - this % shells(i) % n_transitions = n_transition - if (n_transition > 0) then - allocate(this % shells(i) % transition_subshells(2, n_transition)) - allocate(this % shells(i) % transition_energy(n_transition)) - allocate(this % shells(i) % transition_probability(n_transition)) - - allocate(matrix(dims2(1), dims2(2))) - call read_dataset(matrix, dset_id) - - this % shells(i) % transition_subshells(:,:) = int(matrix(1:2, :), 4) - this % shells(i) % transition_energy(:) = matrix(3, :) - this % shells(i) % transition_probability(:) = matrix(4, :) & - / sum(matrix(4, :)) - deallocate(matrix) - end if - call close_dataset(dset_id) - else - this % shells(i) % n_transitions = 0 - end if - call close_group(tgroup) - end do - call close_group(rgroup) - deallocate(designators) - - ! Determine number of electron shells - rgroup = open_group(group_id, 'compton_profiles') - - ! Determine number of shells - dset_id = open_dataset(rgroup, 'num_electrons') - call get_shape(dset_id, dims) - n_shell = int(dims(1), 4) - - ! Read electron shell PDF and binding energies - allocate(this % electron_pdf(n_shell), this % binding_energy(n_shell)) - call read_dataset(this % electron_pdf, dset_id) - call close_dataset(dset_id) - call read_dataset(this % binding_energy, rgroup, 'binding_energy') - this % electron_pdf(:) = this % electron_pdf / sum(this % electron_pdf) - - ! Read Compton profiles - dset_id = open_dataset(rgroup, 'J') - call get_shape(dset_id, dims2) - n_profile = int(dims2(1), 4) - allocate(this % profile_pdf(n_profile, n_shell)) - call read_dataset(this % profile_pdf, dset_id) - call close_dataset(dset_id) - - ! Get Compton profile momentum grid - if (.not. allocated(compton_profile_pz)) then - allocate(compton_profile_pz(n_profile)) - call read_dataset(compton_profile_pz, rgroup, 'pz') - end if - call close_group(rgroup) - - ! Create Compton profile CDF - allocate(this % profile_cdf(n_profile, n_shell)) - do i = 1, n_shell - c = ZERO - this % profile_cdf(1,i) = ZERO - do j = 1, n_profile - 1 - c = c + HALF*(compton_profile_pz(j+1) - compton_profile_pz(j)) * & - (this%profile_pdf(j,i) + this%profile_pdf(j+1,i)) - this % profile_cdf(j+1,i) = c - end do - end do - - ! Calculate total pair production - this % pair_production_total(:) = this % pair_production_nuclear + & - this % pair_production_electron - if (electron_treatment == ELECTRON_TTB) then ! Read bremsstrahlung scaled DCS rgroup = open_group(group_id, 'bremsstrahlung') @@ -401,108 +190,8 @@ contains ! interpolated this % energy = log(this % energy) - where (this % coherent > ZERO) - this % coherent = log(this % coherent) - elsewhere - this % coherent = -500.0_8 - end where - - where (this % incoherent > ZERO) - this % incoherent = log(this % incoherent) - elsewhere - this % incoherent = -500.0_8 - end where - - where (this % photoelectric_total > ZERO) - this % photoelectric_total = log(this % photoelectric_total) - elsewhere - this % photoelectric_total = -500.0_8 - end where - - where (this % pair_production_total > ZERO) - this % pair_production_total = log(this % pair_production_total) - elsewhere - this % pair_production_total = -500.0_8 - end where - end subroutine photon_from_hdf5 -!=============================================================================== -! CALCULATE_ELEMENT_XS determines microscopic photon cross sections for an -! element of a given index in the elements array at the energy of the given -! particle -!=============================================================================== - - subroutine photon_calculate_xs(this, E, xs) - class(PhotonInteraction), intent(in) :: this ! index into elements array - real(8), intent(in) :: E ! energy - type(ElementMicroXS), intent(inout) :: xs - - integer :: i_grid ! index on element energy grid - integer :: i_shell ! index in subshells - integer :: i_start ! threshold index - integer :: n_grid ! number of grid points - real(8) :: f ! interp factor on element energy grid - real(8) :: log_E ! logarithm of the energy - - ! Perform binary search on the element energy grid in order to determine - ! which points to interpolate between - n_grid = size(this % energy) - log_E = log(E) - if (log_E <= this % energy(1)) then - i_grid = 1 - elseif (log_E > this % energy(n_grid)) then - i_grid = n_grid - 1 - else - i_grid = binary_search(this % energy, n_grid, log_E) - end if - - ! check for case where two energy points are the same - if (this % energy(i_grid) == this % energy(i_grid+1)) i_grid = i_grid + 1 - - ! calculate interpolation factor - f = (log_E - this % energy(i_grid)) / & - (this % energy(i_grid+1) - this % energy(i_grid)) - - xs % index_grid = i_grid - xs % interp_factor = f - - ! Calculate microscopic coherent cross section - xs % coherent = exp(this % coherent(i_grid) + f * & - (this % coherent(i_grid+1) - this % coherent(i_grid))) - - ! Calculate microscopic incoherent cross section - xs % incoherent = exp(this % incoherent(i_grid) + & - f*(this % incoherent(i_grid+1) - this % incoherent(i_grid))) - - ! Calculate microscopic photoelectric cross section - xs % photoelectric = ZERO - do i_shell = 1, size(this % shells) - ! Check threshold of reaction - i_start = this % shells(i_shell) % threshold - if (i_grid <= i_start) cycle - - ! Evaluation subshell photoionization cross section - xs % photoelectric = xs % photoelectric + & - exp(this % shells(i_shell) % cross_section(i_grid-i_start) + & - f*(this % shells(i_shell) % cross_section(i_grid+1-i_start) - & - this % shells(i_shell) % cross_section(i_grid-i_start))) - end do - - ! Calculate microscopic pair production cross section - xs % pair_production = exp(& - this % pair_production_total(i_grid) + f*(& - this % pair_production_total(i_grid+1) - & - this % pair_production_total(i_grid))) - - ! Calculate microscopic total cross section - xs % total = xs % coherent + xs % incoherent + xs % photoelectric + & - xs % pair_production - - xs % last_E = E - - end subroutine photon_calculate_xs - !=============================================================================== ! FREE_MEMORY_PHOTON deallocates/resets global variables in this module !=============================================================================== diff --git a/src/physics.cpp b/src/physics.cpp index a4bd670660..dbc6252eb0 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -297,8 +297,7 @@ void sample_photon_reaction(Particle* p) if (prob_after > cutoff) { for (const auto& shell : element.shells_) { // Get grid index and interpolation factor - // TODO: off-by-one - int i_grid = micro.index_grid - 1; + int i_grid = micro.index_grid; double f = micro.interp_factor; // Check threshold of reaction From cf48405ee9b54644d9df1ca647b310b2635aeaf4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 9 Jan 2019 23:36:00 -0600 Subject: [PATCH 05/14] Fix particle restart runs --- src/particle_restart.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 78a558d346..9474f58896 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -10,6 +10,7 @@ module particle_restart use nuclide_header, only: micro_xs, n_nuclides use output, only: print_particle use particle_header + use photon_header, only: micro_photon_xs, n_elements use random_lcg, only: set_particle_seed use settings use simulation_header @@ -45,7 +46,8 @@ contains !$omp parallel allocate(micro_xs(n_nuclides)) - !$omp end parallel + allocate(micro_photon_xs(n_elements)) + !$omp end parallel call set_micro_xs() ! Initialize the particle to be tracked From 1992ebc810a3070fda88e23ab46b31eb2a7957cf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 9 Jan 2019 23:36:10 -0600 Subject: [PATCH 06/14] Update list of developers, publications --- docs/source/developers.rst | 5 +++++ docs/source/publications.rst | 7 ++++++- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/docs/source/developers.rst b/docs/source/developers.rst index 825fbd1c6d..7d32b5831c 100644 --- a/docs/source/developers.rst +++ b/docs/source/developers.rst @@ -18,6 +18,11 @@ code: * `Jingang Liang `_ * `Colin Josey `_ * `Amanda Lund `_ +* `Guillaume Giudicelli `_ +* `Isaac Meyer `_ +* `Patrick Shriwise `_ +* `Shikhar Kumar `_ +* `Andrew Davis `_ * `Benoit Forget `_ * `Kord Smith `_ * `Andrew Siegel `_ diff --git a/docs/source/publications.rst b/docs/source/publications.rst index b88bdc0f0c..49430d84b7 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -57,9 +57,14 @@ Benchmarking Coupling and Multi-physics -------------------------- +- Ze-Long Zhao, Yongwei Yang, and Shuang Hong, "`Application of FLUKA and OpenMC + in coupled physics calculation of target and subcritical reactor for ADS + `_," *Nucl. Sci. Tech.*, **30** + (2019). + - Jun Chen, Liangzhi Cao, Chuanqi Zhao, and Zhouyu Liu, "`Development of Subchannel Code SUBSC for high-fidelity multi-physics coupling application - `_", Energy Procedia, **127**, + `_", *Energy Procedia*, **127**, 264-274 (2017). - Tianliang Hu, Liangzhu Cao, Hongchun Wu, Xianan Du, and Mingtao He, "`Coupled From ea1f2e426dbb79abe309e7d21e80cf98d28ee642 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 11 Jan 2019 16:31:22 -0600 Subject: [PATCH 07/14] Start converting TTB data/physics --- include/openmc/bremsstrahlung.h | 44 +++++++ include/openmc/material.h | 8 ++ include/openmc/photon.h | 4 - src/bremsstrahlung.cpp | 109 +++++++++++++++++ src/material.cpp | 203 +++++++++++++++++++++++++++++++- src/material_header.F90 | 6 + src/photon.cpp | 87 -------------- 7 files changed, 368 insertions(+), 93 deletions(-) create mode 100644 include/openmc/bremsstrahlung.h create mode 100644 src/bremsstrahlung.cpp diff --git a/include/openmc/bremsstrahlung.h b/include/openmc/bremsstrahlung.h new file mode 100644 index 0000000000..ac74a2ab47 --- /dev/null +++ b/include/openmc/bremsstrahlung.h @@ -0,0 +1,44 @@ +#ifndef OPENMC_BREMSSTRAHLUNG_H +#define OPENMC_BREMSSTRAHLUNG_H + +#include "openmc/particle.h" + +#include "xtensor/xtensor.hpp" + +namespace openmc { + +//============================================================================== +// Bremsstrahlung classes +//============================================================================== + +class BremsstrahlungData { +public: + // Data + xt::xtensor pdf; //!< Bremsstrahlung energy PDF + xt::xtensor cdf; //!< Bremsstrahlung energy CDF + xt::xtensor yield; //!< Photon yield +}; + +class Bremsstrahlung { +public: + // Data + BremsstrahlungData electron; + BremsstrahlungData positron; +}; + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +extern xt::xtensor ttb_e_grid; //! energy T of incident electron in [eV] +extern xt::xtensor ttb_k_grid; //! reduced energy W/T of emitted photon + +} // namespace data + +void thick_target_bremsstrahlung(Particle& p, double* E_lost); + +} // namespace openmc + +#endif // OPENMC_BREMSSTRAHLUNG_H diff --git a/include/openmc/material.h b/include/openmc/material.h index ad9b9c9f5a..8df8fa3087 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -1,11 +1,14 @@ #ifndef OPENMC_MATERIAL_H #define OPENMC_MATERIAL_H +#include // for unique_ptr #include #include #include "pugixml.hpp" +#include "openmc/bremsstrahlung.h" +#include "openmc/particle.h" namespace openmc { @@ -38,9 +41,14 @@ public: //! A negative value indicates no default temperature was specified. double temperature_ {-1}; + std::unique_ptr ttb_; + Material() {}; explicit Material(pugi::xml_node material_node); + +private: + void init_bremsstrahlung(); }; //============================================================================== diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 81675585ff..592ae6b973 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -119,8 +119,6 @@ struct ElementMicroXS { std::pair klein_nishina(double alpha); -void thick_target_bremsstrahlung(Particle& p, double* E_lost); - //============================================================================== // Global variables //============================================================================== @@ -128,8 +126,6 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost); namespace data { extern xt::xtensor compton_profile_pz; //! Compton profile momentum grid -extern xt::xtensor ttb_e_grid; //! energy T of incident electron in [eV] -extern xt::xtensor ttb_k_grid; //! reduced energy W/T of emitted photon //! Photon interaction data for each element extern std::vector elements; diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp new file mode 100644 index 0000000000..1ccf1f3517 --- /dev/null +++ b/src/bremsstrahlung.cpp @@ -0,0 +1,109 @@ +#include "openmc/bremsstrahlung.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +xt::xtensor ttb_e_grid; +xt::xtensor ttb_k_grid; +std::vector ttb; + +} // namespace data + +//============================================================================== +// Non-member functions +//============================================================================== + + +// void thick_target_bremsstrahlung(Particle& p, double* E_lost) +// { +// if (p.material == MATERIAL_VOID) return; + +// // TODO: off-by-one +// int photon = static_cast(ParticleType::photon) - 1; +// if (p.E < settings::energy_cutoff[photon]) return; + +// // // Get bremsstrahlung data for this material and particle type +// if (p.type == static_cast(ParticleType::positron)) { +// mat => ttb(p.material) % positron; +// } else { +// mat => ttb(p.material) % electron; +// } + +// double e = std::log(p.E); +// auto n_e = data::ttb_e_grid + +// // Find the lower bounding index of the incident electron energy +// int j = lower_bound_index(data::ttb_e_grid.cbegin(), data::ttb_e_grid.cend(), e); +// if (j == n_e - 1) --j; + +// // Get the interpolation bounds +// double e_l = data::ttb_e_grid(j); +// double e_r = data::ttb_e_grid(j+1); +// double y_l = mat % yield(j); +// double y_r = mat % yield(j+1); + +// // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF +// // interpolated in log energy, which can be interpreted as the probability +// // of index j+1 +// double f = (e - e_l)/(e_r - e_l); + +// // Get the photon number yield for the given energy using linear +// // interpolation on a log-log scale +// double y = std::exp(y_l + (y_r - y_l)*f); + +// // Sample number of secondary bremsstrahlung photons +// int n = y + prn(); + +// *E_lost = 0.0; +// if (n == 0) return; + +// // Sample index of the tabulated PDF in the energy grid, j or j+1 +// double c_max; +// if (prn() <= f || j == 0) { +// int i_e = j + 1; + +// // Interpolate the maximum value of the CDF at the incoming particle +// // energy on a log-log scale +// double p_l = mat % pdf(i_e-1, i_e); +// double p_r = mat % pdf(i_e, i_e); +// double c_l = mat % cdf(i_e-1, i_e); +// double a = std::log(p_r/p_l)/(e_r - e_l) + 1.0; +// c_max = c_l + std::exp(e_l)*p_l/a*(std::exp(a*(e - e_l)) - 1.0); +// } else { +// int i_e = j; + +// // Maximum value of the CDF +// c_max = mat % cdf(i_e, i_e) +// } + +// // Sample the energies of the emitted photons +// for (int i = 0; i < n; ++i) { +// // Generate a random number r and determine the index i for which +// // cdf(i) <= r*cdf,max <= cdf(i+1) +// double c = prn()*c_max; +// int i_w = lower_bound_index(mat % cdf(:i_e,i_e), c) + +// // Sample the photon energy +// double w_l = data::ttb_e_grid(i_w); +// double w_r = data::ttb_e_grid(i_w+1); +// double p_l = mat % pdf(i_w, i_e); +// double p_r = mat % pdf(i_w+1, i_e); +// double c_l = mat % cdf(i_w, i_e); +// double a = std::log(p_r/p_l)/(w_r - w_l) + 1.0; +// double w = std::exp(w_l)*std::pow(a*(c - c_l)/(std::exp(w_l)*p_l) + 1.0, 1.0/a); + +// if (w > settings::energy_cutoff[photon]) { +// // Create secondary photon +// int photon = static_cast(ParticleType::photon); +// p.create_secondary(p.coord[0].uvw, w, photon, true); +// *E_lost += w; +// } +// } +// } + +} // namespace openmc diff --git a/src/material.cpp b/src/material.cpp index 7f8cfc50ba..9680949944 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1,11 +1,19 @@ #include "openmc/material.h" +#include #include #include -#include "openmc/error.h" -#include "openmc/xml_interface.h" +#include "xtensor/xbuilder.hpp" +#include "xtensor/xoperation.hpp" +#include "xtensor/xview.hpp" +#include "openmc/error.h" +#include "openmc/math_functions.h" +#include "openmc/nuclide.h" +#include "openmc/photon.h" +#include "openmc/search.h" +#include "openmc/xml_interface.h" namespace openmc { @@ -41,6 +49,197 @@ Material::Material(pugi::xml_node node) } } +int* material_element(int index); + +void Material::init_bremsstrahlung() +{ + // Create new object + ttb_ = std::make_unique(); + + // Get the size of the energy grids + auto n_k = data::ttb_k_grid.size(); + auto n_e = data::ttb_e_grid.size(); + + // Get pointers to nuclides, elements, densities + int32_t index; + openmc_get_material_index(id_, &index); + int* nuclide_; + double* atom_density_; + int n; + openmc_material_get_densities(index, &nuclide_, &atom_density_, &n); + int* element_ = material_element(index); + + for (int particle = 0; particle < 2; ++particle) { + // Loop over logic twice, once for electron, once for positron + BremsstrahlungData* ttb = (particle == 0) ? &ttb_->electron : &ttb_->positron; + bool positron = (particle == 1); + + // Allocate arrays for TTB data + ttb->pdf = xt::zeros({n_e, n_e}); + ttb->cdf = xt::zeros({n_e, n_e}); + ttb->yield = xt::empty({n_e}); + + // Allocate temporary arrays + std::array shape {n_e}; + xt::xtensor stopping_power_collision({n_e}, 0.0); + xt::xtensor stopping_power_radiative({n_e}, 0.0); + xt::xtensor dcs({n_e, n_k}, 0.0); + xt::xtensor f({n_e}, 0.0); + xt::xtensor z({n_e}, 0.0); + + double Z_eq_sq = 0.0; + double sum_density = 0.0; + + // Calculate the molecular DCS and the molecular total stopping power using + // Bragg's additivity rule. + // TODO: The collision stopping power cannot be accurately calculated using + // Bragg's additivity rule since the mean excitation energies and the + // density effect corrections cannot simply be summed together. Bragg's + // additivity rule fails especially when a higher-density compound is + // composed of elements that are in lower-density form at normal temperature + // and pressure (at which the NIST stopping powers are given). It will be + // used to approximate the collision stopping powers for now, but should be + // fixed in the future. + for (int i = 0; i < n; ++i) { + // Get pointer to current element + // TODO: off-by-one + const auto& elm = data::elements[element_[i] - 1]; + // TODO: off-by-one + double awr = data::nuclides[nuclide_[i] - 1]->awr_; + + // Get atomic density and mass density of nuclide given atom/weight percent + double atom_density = (atom_density_[0] > 0.0) ? + atom_density_[i] : -atom_density_[i] / awr; + double mass_density = atom_density * awr; + + // Calculate the "equivalent" atomic number Zeq of the material + Z_eq_sq += atom_density * elm.Z_ * elm.Z_; + sum_density += atom_density; + + // Accumulate material DCS + dcs += (atom_density * elm.Z_ * elm.Z_) * elm.dcs_; + + // Accumulate material collision stopping power + stopping_power_collision += (mass_density * MASS_NEUTRON / N_AVOGADRO) + * elm.stopping_power_collision_; + + // Accumulate material radiative stopping power + stopping_power_radiative += (mass_density * MASS_NEUTRON / N_AVOGADRO) + * elm.stopping_power_radiative_; + } + Z_eq_sq /= sum_density; + + // Calculate the positron DCS and radiative stopping power. These are + // obtained by multiplying the electron DCS and radiative stopping powers by + // a factor r, which is a numerical approximation of the ratio of the + // radiative stopping powers for positrons and electrons. Source: F. Salvat, + // J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for + // Monte Carlo Simulation of Electron and Photon Transport," OECD-NEA, + // Issy-les-Moulineaux, France (2011). + if (positron) { + for (int i = 0; i < n_e; ++i) { + double t = std::log(1.0 + 1.0e6*data::ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON_EV)); + double r = 1.0 - std::exp(-1.2359e-1*t + 6.1274e-2*std::pow(t, 2) + - 3.1516e-2*std::pow(t, 3) + 7.7446e-3*std::pow(t, 4) + - 1.0595e-3*std::pow(t, 5) + 7.0568e-5*std::pow(t, 6) + - 1.808e-6*std::pow(t, 7)); + stopping_power_radiative(i) *= r; + auto dcs_i = xt::view(dcs, i, xt::all()); + dcs_i *= r; + } + } + + // Total material stopping power + xt::xtensor stopping_power = stopping_power_collision + + stopping_power_radiative; + + // Loop over photon energies + for (int i = 0; i < n_e - 1; ++i) { + double w = data::ttb_e_grid(i); + + // Loop over incident particle energies + for (int j = 0; j < n_e; ++j) { + double e = data::ttb_e_grid(j); + + // Reduced photon energy + double k = w / e; + + // Find the lower bounding index of the reduced photon energy + int i_k = lower_bound_index(data::ttb_k_grid.cbegin(), data::ttb_k_grid.cend(), k); + + // Get the interpolation bounds + double k_l = data::ttb_k_grid(i_k); + double k_r = data::ttb_k_grid(i_k + 1); + double x_l = dcs(j, i_k); + double x_r = dcs(j, i_k+1); + + // Find the value of the DCS using linear interpolation in reduced + // photon energy k + double x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l); + + // Ratio of the velocity of the charged particle to the speed of light + double beta = std::sqrt(e*(e + 2.0*MASS_ELECTRON_EV)) / (e + MASS_ELECTRON_EV); + + // Compute the integrand of the PDF + f(j) = x / (beta*beta * stopping_power(j) * w); + } + + // Number of points to integrate + int n = n_e - i; + + // Integrate the PDF using cubic spline integration over the incident + // particle energy + if (n > 2) { + spline_c(n, &data::ttb_e_grid(i), &f(i), &z(i)); + + double c = 0.0; + for (int j = 0; j < n_e - 1; ++j) { + c += spline_integrate_c(n, &data::ttb_e_grid(i), &f(i), &z(i), + data::ttb_e_grid(j), data::ttb_e_grid(j+1)); + + ttb->pdf(j+1,i) = c; + } + + // Integrate the last two points using trapezoidal rule in log-log space + } else { + double e_l = std::log(data::ttb_e_grid(i)); + double e_r = std::log(data::ttb_e_grid(i+1)); + double x_l = std::log(f(i)); + double x_r = std::log(f(i+1)); + + ttb->pdf(i,i+1) = 0.5*(e_r - e_l)*(std::exp(e_l + x_l) + std::exp(e_r + x_r)); + } + } + + // Loop over incident particle energies + for (int j = 1; j < n_e; ++j) { + // Set last element of PDF to small non-zero value to enable log-log + // interpolation + ttb->pdf(j,j) = std::exp(-500.0); + + // Loop over photon energies + double c = 0.0; + for (int i = 0; i < j - 1; ++i) { + // Integrate the CDF from the PDF using the trapezoidal rule in log-log + // space + double w_l = std::log(data::ttb_e_grid(i)); + double w_r = std::log(data::ttb_e_grid(i+1)); + double x_l = std::log(ttb->pdf(j,i)); + double x_r = std::log(ttb->pdf(j,i+1)); + + c += 0.5*(w_r - w_l)*(std::exp(w_l + x_l) + std::exp(w_r + x_r)); + ttb->cdf(j,i+1) = c; + } + + // Set photon number yield + ttb->yield(j) = c; + } + + // Use logarithm of number yield since it is log-log interpolated + ttb->yield = xt::where(ttb->yield > 0.0, ttb->yield, -500.0); + } +} + //============================================================================== // Non-method functions //============================================================================== diff --git a/src/material_header.F90 b/src/material_header.F90 index 9341732419..ee3f1f8bf4 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -1065,4 +1065,10 @@ contains end associate end function + function material_element(i_material) result(ptr) bind(C) + integer(C_INT), value :: i_material + type(C_PTR) :: ptr + ptr = C_LOC(materials(i_material) % element(1)) + end function + end module material_header diff --git a/src/photon.cpp b/src/photon.cpp index 7432a2dcd6..b324f5ff3b 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -768,93 +768,6 @@ std::pair klein_nishina(double alpha) return {alpha_out, mu}; } -// void thick_target_bremsstrahlung(Particle& p, double* E_lost) -// { -// if (p.material == MATERIAL_VOID) return; - -// // TODO: off-by-one -// int photon = static_cast(ParticleType::photon) - 1; -// if (p.E < settings::energy_cutoff[photon]) return; - -// // // Get bremsstrahlung data for this material and particle type -// if (p.type == static_cast(ParticleType::positron)) { -// mat => ttb(p.material) % positron; -// } else { -// mat => ttb(p.material) % electron; -// } - -// double e = std::log(p.E); -// auto n_e = data::ttb_e_grid - -// // Find the lower bounding index of the incident electron energy -// int j = lower_bound_index(data::ttb_e_grid.cbegin(), data::ttb_e_grid.cend(), e); -// if (j == n_e - 1) --j; - -// // Get the interpolation bounds -// double e_l = data::ttb_e_grid(j); -// double e_r = data::ttb_e_grid(j+1); -// double y_l = mat % yield(j); -// double y_r = mat % yield(j+1); - -// // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF -// // interpolated in log energy, which can be interpreted as the probability -// // of index j+1 -// double f = (e - e_l)/(e_r - e_l); - -// // Get the photon number yield for the given energy using linear -// // interpolation on a log-log scale -// double y = std::exp(y_l + (y_r - y_l)*f); - -// // Sample number of secondary bremsstrahlung photons -// int n = y + prn(); - -// *E_lost = 0.0; -// if (n == 0) return; - -// // Sample index of the tabulated PDF in the energy grid, j or j+1 -// double c_max; -// if (prn() <= f || j == 0) { -// int i_e = j + 1; - -// // Interpolate the maximum value of the CDF at the incoming particle -// // energy on a log-log scale -// double p_l = mat % pdf(i_e-1, i_e); -// double p_r = mat % pdf(i_e, i_e); -// double c_l = mat % cdf(i_e-1, i_e); -// double a = std::log(p_r/p_l)/(e_r - e_l) + 1.0; -// c_max = c_l + std::exp(e_l)*p_l/a*(std::exp(a*(e - e_l)) - 1.0); -// } else { -// int i_e = j; - -// // Maximum value of the CDF -// c_max = mat % cdf(i_e, i_e) -// } - -// // Sample the energies of the emitted photons -// for (int i = 0; i < n; ++i) { -// // Generate a random number r and determine the index i for which -// // cdf(i) <= r*cdf,max <= cdf(i+1) -// double c = prn()*c_max; -// int i_w = lower_bound_index(mat % cdf(:i_e,i_e), c) - -// // Sample the photon energy -// double w_l = data::ttb_e_grid(i_w); -// double w_r = data::ttb_e_grid(i_w+1); -// double p_l = mat % pdf(i_w, i_e); -// double p_r = mat % pdf(i_w+1, i_e); -// double c_l = mat % cdf(i_w, i_e); -// double a = std::log(p_r/p_l)/(w_r - w_l) + 1.0; -// double w = std::exp(w_l)*std::pow(a*(c - c_l)/(std::exp(w_l)*p_l) + 1.0, 1.0/a); - -// if (w > settings::energy_cutoff[photon]) { -// // Create secondary photon -// int photon = static_cast(ParticleType::photon); -// p.create_secondary(p.coord[0].uvw, w, photon, true); -// *E_lost += w; -// } -// } -// } - //============================================================================== // Fortran compatibility //============================================================================== From 2eacd5f72ce8f9cd6ac68d5fb1717dfc7ece134a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 12 Jan 2019 17:04:36 -0600 Subject: [PATCH 08/14] Finish TTB implementation on C++ side, still a bug somewhere --- CMakeLists.txt | 1 + include/openmc/bremsstrahlung.h | 4 + include/openmc/material.h | 3 +- include/openmc/physics.h | 2 - src/bremsstrahlung.cpp | 162 ++++++++++++++++++-------------- src/input_xml.F90 | 13 ++- src/material.cpp | 32 ++++--- src/photon.cpp | 6 +- src/physics.cpp | 5 +- 9 files changed, 133 insertions(+), 95 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7e9e24161c..808cc8e2fb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -370,6 +370,7 @@ add_library(libopenmc SHARED src/tallies/trigger.F90 src/tallies/trigger_header.F90 src/bank.cpp + src/bremsstrahlung.cpp src/dagmc.cpp src/cell.cpp src/cmfd_solver.cpp diff --git a/include/openmc/bremsstrahlung.h b/include/openmc/bremsstrahlung.h index ac74a2ab47..bd9d012c1b 100644 --- a/include/openmc/bremsstrahlung.h +++ b/include/openmc/bremsstrahlung.h @@ -37,6 +37,10 @@ extern xt::xtensor ttb_k_grid; //! reduced energy W/T of emitted phot } // namespace data +//============================================================================== +// Global variables +//============================================================================== + void thick_target_bremsstrahlung(Particle& p, double* E_lost); } // namespace openmc diff --git a/include/openmc/material.h b/include/openmc/material.h index 8df8fa3087..d796b08ed6 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -47,7 +47,7 @@ public: explicit Material(pugi::xml_node material_node); -private: + //! Initialize bremsstrahlung data void init_bremsstrahlung(); }; @@ -55,6 +55,7 @@ private: // Fortran compatibility //============================================================================== +extern "C" int* material_element(int i_material); extern "C" bool material_isotropic(int i_material, int i_nuc_mat); } // namespace openmc diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 8a3d7c4ef1..921494db7a 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -82,8 +82,6 @@ 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/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index 1ccf1f3517..1ba1073a5d 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -1,5 +1,13 @@ #include "openmc/bremsstrahlung.h" +#include "openmc/constants.h" +#include "openmc/material.h" +#include "openmc/random_lcg.h" +#include "openmc/search.h" +#include "openmc/settings.h" + +#include "xtensor/xmath.hpp" + namespace openmc { //============================================================================== @@ -18,92 +26,104 @@ std::vector ttb; // Non-member functions //============================================================================== +void thick_target_bremsstrahlung(Particle& p, double* E_lost) +{ + if (p.material == MATERIAL_VOID) return; -// void thick_target_bremsstrahlung(Particle& p, double* E_lost) -// { -// if (p.material == MATERIAL_VOID) return; + // TODO: off-by-one + int photon = static_cast(ParticleType::photon) - 1; + if (p.E < settings::energy_cutoff[photon]) return; -// // TODO: off-by-one -// int photon = static_cast(ParticleType::photon) - 1; -// if (p.E < settings::energy_cutoff[photon]) return; + // Get bremsstrahlung data for this material and particle type + // TODO: off-by-one + BremsstrahlungData* mat; + if (p.type == static_cast(ParticleType::positron)) { + mat = &model::materials[p.material -1]->ttb_->positron; + } else { + mat = &model::materials[p.material -1]->ttb_->electron; + } -// // // Get bremsstrahlung data for this material and particle type -// if (p.type == static_cast(ParticleType::positron)) { -// mat => ttb(p.material) % positron; -// } else { -// mat => ttb(p.material) % electron; -// } + double e = std::log(p.E); + auto n_e = data::ttb_e_grid.size(); -// double e = std::log(p.E); -// auto n_e = data::ttb_e_grid + // Find the lower bounding index of the incident electron energy + int j = lower_bound_index(data::ttb_e_grid.cbegin(), + data::ttb_e_grid.cend(), e); + if (j == n_e - 1) --j; -// // Find the lower bounding index of the incident electron energy -// int j = lower_bound_index(data::ttb_e_grid.cbegin(), data::ttb_e_grid.cend(), e); -// if (j == n_e - 1) --j; + // Get the interpolation bounds + double e_l = data::ttb_e_grid(j); + double e_r = data::ttb_e_grid(j+1); + double y_l = mat->yield(j); + double y_r = mat->yield(j+1); -// // Get the interpolation bounds -// double e_l = data::ttb_e_grid(j); -// double e_r = data::ttb_e_grid(j+1); -// double y_l = mat % yield(j); -// double y_r = mat % yield(j+1); + // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF + // interpolated in log energy, which can be interpreted as the probability + // of index j+1 + double f = (e - e_l)/(e_r - e_l); -// // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF -// // interpolated in log energy, which can be interpreted as the probability -// // of index j+1 -// double f = (e - e_l)/(e_r - e_l); + // Get the photon number yield for the given energy using linear + // interpolation on a log-log scale + double y = std::exp(y_l + (y_r - y_l)*f); -// // Get the photon number yield for the given energy using linear -// // interpolation on a log-log scale -// double y = std::exp(y_l + (y_r - y_l)*f); + // Sample number of secondary bremsstrahlung photons + int n = y + prn(); -// // Sample number of secondary bremsstrahlung photons -// int n = y + prn(); + *E_lost = 0.0; + if (n == 0) return; -// *E_lost = 0.0; -// if (n == 0) return; + // Sample index of the tabulated PDF in the energy grid, j or j+1 + double c_max; + int i_e; + if (prn() <= f || j == 0) { + i_e = j + 1; -// // Sample index of the tabulated PDF in the energy grid, j or j+1 -// double c_max; -// if (prn() <= f || j == 0) { -// int i_e = j + 1; + // Interpolate the maximum value of the CDF at the incoming particle + // energy on a log-log scale + double p_l = mat->pdf(i_e, i_e - 1); + double p_r = mat->pdf(i_e, i_e); + double c_l = mat->cdf(i_e, i_e - 1); + double a = std::log(p_r/p_l)/(e_r - e_l) + 1.0; + c_max = c_l + std::exp(e_l)*p_l/a*(std::exp(a*(e - e_l)) - 1.0); + } else { + i_e = j; -// // Interpolate the maximum value of the CDF at the incoming particle -// // energy on a log-log scale -// double p_l = mat % pdf(i_e-1, i_e); -// double p_r = mat % pdf(i_e, i_e); -// double c_l = mat % cdf(i_e-1, i_e); -// double a = std::log(p_r/p_l)/(e_r - e_l) + 1.0; -// c_max = c_l + std::exp(e_l)*p_l/a*(std::exp(a*(e - e_l)) - 1.0); -// } else { -// int i_e = j; + // Maximum value of the CDF + c_max = mat->cdf(i_e, i_e); + } -// // Maximum value of the CDF -// c_max = mat % cdf(i_e, i_e) -// } + // Sample the energies of the emitted photons + for (int i = 0; i < n; ++i) { + // Generate a random number r and determine the index i for which + // cdf(i) <= r*cdf,max <= cdf(i+1) + double c = prn()*c_max; + int i_w = lower_bound_index(&mat->cdf(i_e, 0), &mat->cdf(i_e, 0) + i_e, c); -// // Sample the energies of the emitted photons -// for (int i = 0; i < n; ++i) { -// // Generate a random number r and determine the index i for which -// // cdf(i) <= r*cdf,max <= cdf(i+1) -// double c = prn()*c_max; -// int i_w = lower_bound_index(mat % cdf(:i_e,i_e), c) + // Sample the photon energy + double w_l = data::ttb_e_grid(i_w); + double w_r = data::ttb_e_grid(i_w + 1); + double p_l = mat->pdf(i_e, i_w); + double p_r = mat->pdf(i_e, i_w + 1); + double c_l = mat->cdf(i_e, i_w); + double a = std::log(p_r/p_l)/(w_r - w_l) + 1.0; + double w = std::exp(w_l)*std::pow(a*(c - c_l)/(std::exp(w_l)*p_l) + 1.0, 1.0/a); -// // Sample the photon energy -// double w_l = data::ttb_e_grid(i_w); -// double w_r = data::ttb_e_grid(i_w+1); -// double p_l = mat % pdf(i_w, i_e); -// double p_r = mat % pdf(i_w+1, i_e); -// double c_l = mat % cdf(i_w, i_e); -// double a = std::log(p_r/p_l)/(w_r - w_l) + 1.0; -// double w = std::exp(w_l)*std::pow(a*(c - c_l)/(std::exp(w_l)*p_l) + 1.0, 1.0/a); + if (w > settings::energy_cutoff[photon]) { + // Create secondary photon + int photon_ = static_cast(ParticleType::photon); + p.create_secondary(p.coord[0].uvw, w, photon_, true); + *E_lost += w; + } + } +} -// if (w > settings::energy_cutoff[photon]) { -// // Create secondary photon -// int photon = static_cast(ParticleType::photon); -// p.create_secondary(p.coord[0].uvw, w, photon, true); -// *E_lost += w; -// } -// } -// } +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" void set_log_ttb_e_grid() +{ + data::ttb_e_grid = xt::log(data::ttb_e_grid); +} } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index be93a164e0..4b6d774944 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2034,6 +2034,15 @@ contains type(SetChar) :: already_read type(SetChar) :: element_already_read + interface + subroutine material_init_bremsstrahlung(ptr) bind(C) + import C_PTR + type(C_PTR), value :: ptr + end subroutine + subroutine set_log_ttb_e_grid() bind(C) + end subroutine + end interface + allocate(nuclides(n_nuclides)) allocate(elements(n_elements)) if (photon_transport .and. electron_treatment == ELECTRON_TTB) then @@ -2127,8 +2136,7 @@ contains ! Generate material bremsstrahlung data for electrons and positrons if (photon_transport .and. electron_treatment == ELECTRON_TTB) then - call bremsstrahlung_init(ttb(i) % electron, i, ELECTRON) - call bremsstrahlung_init(ttb(i) % positron, i, POSITRON) + call material_init_bremsstrahlung(materials(i) % ptr) end if end do @@ -2155,6 +2163,7 @@ contains ! Take logarithm of energies since they are log-log interpolated ttb_e_grid = log(ttb_e_grid) + call set_log_ttb_e_grid() end if ! Set up logarithmic grid for nuclides diff --git a/src/material.cpp b/src/material.cpp index 9680949944..acbfc08f54 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -49,8 +49,6 @@ Material::Material(pugi::xml_node node) } } -int* material_element(int index); - void Material::init_bremsstrahlung() { // Create new object @@ -80,12 +78,9 @@ void Material::init_bremsstrahlung() ttb->yield = xt::empty({n_e}); // Allocate temporary arrays - std::array shape {n_e}; xt::xtensor stopping_power_collision({n_e}, 0.0); xt::xtensor stopping_power_radiative({n_e}, 0.0); xt::xtensor dcs({n_e, n_k}, 0.0); - xt::xtensor f({n_e}, 0.0); - xt::xtensor z({n_e}, 0.0); double Z_eq_sq = 0.0; double sum_density = 0.0; @@ -154,31 +149,35 @@ void Material::init_bremsstrahlung() stopping_power_radiative; // Loop over photon energies + xt::xtensor f({n_e}, 0.0); + xt::xtensor z({n_e}, 0.0); for (int i = 0; i < n_e - 1; ++i) { double w = data::ttb_e_grid(i); // Loop over incident particle energies - for (int j = 0; j < n_e; ++j) { + for (int j = i; j < n_e; ++j) { double e = data::ttb_e_grid(j); // Reduced photon energy double k = w / e; // Find the lower bounding index of the reduced photon energy - int i_k = lower_bound_index(data::ttb_k_grid.cbegin(), data::ttb_k_grid.cend(), k); + int i_k = lower_bound_index(data::ttb_k_grid.cbegin(), + data::ttb_k_grid.cend(), k); // Get the interpolation bounds double k_l = data::ttb_k_grid(i_k); double k_r = data::ttb_k_grid(i_k + 1); double x_l = dcs(j, i_k); - double x_r = dcs(j, i_k+1); + double x_r = dcs(j, i_k + 1); // Find the value of the DCS using linear interpolation in reduced // photon energy k - double x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l); + double x = x_l + (k - k_l)*(x_r - x_l)/(k_r - k_l); // Ratio of the velocity of the charged particle to the speed of light - double beta = std::sqrt(e*(e + 2.0*MASS_ELECTRON_EV)) / (e + MASS_ELECTRON_EV); + double beta = std::sqrt(e*(e + 2.0*MASS_ELECTRON_EV)) / + (e + MASS_ELECTRON_EV); // Compute the integrand of the PDF f(j) = x / (beta*beta * stopping_power(j) * w); @@ -193,7 +192,7 @@ void Material::init_bremsstrahlung() spline_c(n, &data::ttb_e_grid(i), &f(i), &z(i)); double c = 0.0; - for (int j = 0; j < n_e - 1; ++j) { + for (int j = i; j < n_e - 1; ++j) { c += spline_integrate_c(n, &data::ttb_e_grid(i), &f(i), &z(i), data::ttb_e_grid(j), data::ttb_e_grid(j+1)); @@ -207,7 +206,7 @@ void Material::init_bremsstrahlung() double x_l = std::log(f(i)); double x_r = std::log(f(i+1)); - ttb->pdf(i,i+1) = 0.5*(e_r - e_l)*(std::exp(e_l + x_l) + std::exp(e_r + x_r)); + ttb->pdf(i+1,i) = 0.5*(e_r - e_l)*(std::exp(e_l + x_l) + std::exp(e_r + x_r)); } } @@ -219,7 +218,7 @@ void Material::init_bremsstrahlung() // Loop over photon energies double c = 0.0; - for (int i = 0; i < j - 1; ++i) { + for (int i = 0; i < j; ++i) { // Integrate the CDF from the PDF using the trapezoidal rule in log-log // space double w_l = std::log(data::ttb_e_grid(i)); @@ -236,7 +235,7 @@ void Material::init_bremsstrahlung() } // Use logarithm of number yield since it is log-log interpolated - ttb->yield = xt::where(ttb->yield > 0.0, ttb->yield, -500.0); + ttb->yield = xt::where(ttb->yield > 0.0, xt::log(ttb->yield), -500.0); } } @@ -332,6 +331,11 @@ extern "C" { mat->fissionable = fissionable; } + void material_init_bremsstrahlung(Material* mat) + { + mat->init_bremsstrahlung(); + } + void extend_materials_c(int32_t n) { model::materials.reserve(model::materials.size() + n); diff --git a/src/photon.cpp b/src/photon.cpp index b324f5ff3b..3b489e2995 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -1,5 +1,6 @@ #include "openmc/photon.h" +#include "openmc/bremsstrahlung.h" #include "openmc/constants.h" #include "openmc/hdf5_interface.h" #include "openmc/particle.h" @@ -24,8 +25,6 @@ namespace openmc { namespace data { xt::xtensor compton_profile_pz; -xt::xtensor ttb_e_grid; -xt::xtensor ttb_k_grid; std::vector elements; @@ -250,7 +249,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) frst, xt::view(s_rad, xt::range(i_grid+1, n_e)))); // Interpolate bremsstrahlung DCS at the cutoff energy and truncate - xt::xtensor dcs = xt::empty({n_e - i_grid, n_k}); + xt::xtensor dcs({n_e - i_grid, n_k}); for (int i = 0; i < n_k; ++i) { y = std::exp(std::log(dcs_(i_grid,i)) + f*(std::log(dcs_(i_grid+1,i)) - std::log(dcs_(i_grid,i)))); @@ -268,6 +267,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) } // Set incident particle energy grid + // TODO: Change to zero when xtensor is updated if (data::ttb_e_grid.size() == 1) { data::ttb_e_grid = electron_energy; } diff --git a/src/physics.cpp b/src/physics.cpp index dbc6252eb0..564bc0f5cb 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -1,6 +1,7 @@ #include "openmc/physics.h" #include "openmc/bank.h" +#include "openmc/bremsstrahlung.h" #include "openmc/constants.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" @@ -382,7 +383,7 @@ void sample_electron_reaction(Particle* p) if (settings::electron_treatment == ELECTRON_TTB) { double E_lost; - thick_target_bremsstrahlung(p, &E_lost); + thick_target_bremsstrahlung(*p, &E_lost); } p->E = 0.0; @@ -395,7 +396,7 @@ void sample_positron_reaction(Particle* p) if (settings::electron_treatment == ELECTRON_TTB) { double E_lost; - thick_target_bremsstrahlung(p, &E_lost); + thick_target_bremsstrahlung(*p, &E_lost); } // Sample angle isotropically From 9f8ebf527b3d1f752e6c503165ee718089469816 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 15 Jan 2019 06:51:18 -0600 Subject: [PATCH 09/14] Remove bremsstrahlung code on Fortran side --- CMakeLists.txt | 1 - include/openmc/math_functions.h | 17 ++- include/openmc/settings.h | 2 +- src/constants.F90 | 26 +--- src/input_xml.F90 | 59 ++------- src/material.cpp | 31 ++++- src/material_header.F90 | 213 -------------------------------- src/math.F90 | 40 ------ src/math_functions.cpp | 14 +-- src/photon.cpp | 24 +++- src/photon_header.F90 | 173 +------------------------- src/photon_physics.F90 | 116 ----------------- src/settings.F90 | 1 - 13 files changed, 81 insertions(+), 636 deletions(-) delete mode 100644 src/photon_physics.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 808cc8e2fb..5bf510e084 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -331,7 +331,6 @@ add_library(libopenmc SHARED src/particle_header.F90 src/particle_restart.F90 src/photon_header.F90 - src/photon_physics.F90 src/physics_common.F90 src/physics.F90 src/pugixml/pugixml_f.F90 diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index f80ed6d3b8..87710fb6b7 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -167,7 +167,7 @@ extern "C" double watt_spectrum(double a, double b); //! //! Samples from a Normal distribution with a given mean and standard deviation //! The PDF is defined as s(x) = (1/2*sigma*sqrt(2) * e-((mu-x)/2*sigma)^2 -//! Its sampled according to +//! Its sampled according to //! http://www-pdg.lbl.gov/2009/reviews/rpp2009-rev-monte-carlo-techniques.pdf //! section 33.4.4 //! @@ -181,7 +181,7 @@ extern "C" double normal_variate(double mean, double std_dev); //============================================================================== //! Samples an energy from the Muir (Gaussian) energy-dependent distribution. //! -//! This is another form of the Gaussian distribution but with more easily +//! This is another form of the Gaussian distribution but with more easily //! modifiable parameters //! https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS //! @@ -193,8 +193,6 @@ extern "C" double normal_variate(double mean, double std_dev); extern "C" double muir_spectrum(double e0, double m_rat, double kt); - - //============================================================================== //! Doppler broadens the windowed multipole curvefit. //! @@ -224,7 +222,7 @@ extern "C" void broaden_wmp_polynomials(double E, double dopp, int n, double fac //! value of x. //============================================================================== -extern "C" void spline_c(int n, const double x[], const double y[], double z[]); +void spline(int n, const double x[], const double y[], double z[]); //============================================================================== //! Determine the cubic spline interpolated y-value for a given x-value. @@ -239,8 +237,8 @@ extern "C" void spline_c(int n, const double x[], const double y[], double z[]); //! \return Interpolated value //============================================================================== -extern "C" double spline_interpolate_c(int n, const double x[], const double y[], - const double z[], double xint); +double spline_interpolate(int n, const double x[], const double y[], + const double z[], double xint); //============================================================================== //! Evaluate the definite integral of the interpolating cubic spline between @@ -257,9 +255,8 @@ extern "C" double spline_interpolate_c(int n, const double x[], const double y[] //! \return Integral //============================================================================== -extern "C" double spline_integrate_c(int n, const double x[], const double y[], - const double z[], double xa, double xb); - +double spline_integrate(int n, const double x[], const double y[], + const double z[], double xa, double xb); //! Evaluate the Faddeeva function //! diff --git a/include/openmc/settings.h b/include/openmc/settings.h index bfd3d056f2..e95bcb107a 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -68,7 +68,7 @@ extern "C" int32_t n_inactive; //!< number of inactive batches extern "C" int32_t gen_per_batch; //!< number of generations per batch extern "C" int64_t n_particles; //!< number of particles per generation -extern "C" int electron_treatment; //!< how to treat secondary electrons +extern int electron_treatment; //!< how to treat secondary electrons extern "C" std::array energy_cutoff; //!< Energy cutoff in [eV] for each particle type extern "C" int legendre_to_tabular_points; //!< number of points to convert Legendres extern "C" int max_order; //!< Maximum Legendre order for multigroup data diff --git a/src/constants.F90 b/src/constants.F90 index efa733ba1f..93a40abdbd 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -68,13 +68,8 @@ module constants real(8), parameter :: & PI = 3.1415926535898_8, & ! pi - SQRT_PI = 1.7724538509055_8, & ! square root of pi MASS_NEUTRON = 1.00866491588_8, & ! mass of a neutron in amu MASS_NEUTRON_EV = 939.5654133e6_8, & ! mass of a neutron in eV/c^2 - MASS_PROTON = 1.007276466879_8, & ! mass of a proton in amu - MASS_ELECTRON_EV = 0.5109989461e6_8, & ! electron mass energy equivalent in eV/c^2 - FINE_STRUCTURE = 137.035999139_8, & ! inverse fine structure constant - PLANCK_C = 1.2398419739062977e4_8,& ! Planck's constant times c in eV-Angstroms AMU = 1.660539040e-27_8, & ! 1 amu in kg C_LIGHT = 2.99792458e8_8, & ! speed of light in m/s N_AVOGADRO = 0.6022140857_8, & ! Avogadro's number in 10^24/mol @@ -83,10 +78,7 @@ module constants ZERO = 0.0_8, & HALF = 0.5_8, & ONE = 1.0_8, & - TWO = 2.0_8, & - THREE = 3.0_8, & - FOUR = 4.0_8 - complex(8), parameter :: ONEI = (ZERO, ONE) + TWO = 2.0_8 ! Electron subshell labels character(3), parameter :: SUBSHELLS(39) = [ & @@ -185,12 +177,6 @@ module constants real(8), parameter :: & MACROSCOPIC_AWR = -TWO - ! Fission neutron emission (nu) type - integer, parameter :: & - NU_NONE = 0, & ! No nu values (non-fissionable) - NU_POLYNOMIAL = 1, & ! Nu values given by polynomial - NU_TABULAR = 2 ! Nu values given by tabular distribution - ! Secondary particle emission type integer, parameter :: & EMISSION_PROMPT = 1, & ! Prompt emission of secondary particle @@ -291,11 +277,6 @@ module constants FILTER_ZERNIKE_RADIAL = 21, & FILTER_PARTICLE = 22 - - ! Mesh types - integer, parameter :: & - MESH_REGULAR = 1 - ! Tally surface current directions integer, parameter :: & OUT_LEFT = 1, & ! x min @@ -382,11 +363,6 @@ module constants MODE_PARTICLE = 4, & ! Particle restart mode MODE_VOLUME = 5 ! Volume calculation mode - ! Electron treatments - integer, parameter :: & - ELECTRON_LED = 1, & ! Local Energy Deposition - ELECTRON_TTB = 2 ! Thick Target Bremsstrahlung - !============================================================================= ! DELAYED NEUTRON PRECURSOR CONSTANTS diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4b6d774944..0883776247 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2035,19 +2035,16 @@ contains type(SetChar) :: element_already_read interface - subroutine material_init_bremsstrahlung(ptr) bind(C) - import C_PTR - type(C_PTR), value :: ptr + subroutine photon_from_hdf5(group) bind(C) + import HID_T + integer(HID_T), value :: group end subroutine - subroutine set_log_ttb_e_grid() bind(C) + + subroutine read_ce_cross_sections_c() bind(C) end subroutine end interface allocate(nuclides(n_nuclides)) - allocate(elements(n_elements)) - if (photon_transport .and. electron_treatment == ELECTRON_TTB) then - allocate(ttb(n_materials)) - end if ! Read cross sections do i = 1, size(materials) @@ -2103,22 +2100,11 @@ contains ! Read element data from HDF5 group_id = open_group(file_id, element) - call elements(i_element) % from_hdf5(group_id) + + call photon_from_hdf5(group_id) call close_group(group_id) call file_close(file_id) - ! Determine if minimum/maximum energy for this element is - ! greater/less than the previous - if (size(elements(i_element) % energy) >= 1) then - energy_min(PHOTON) = max(energy_min(PHOTON), & - exp(elements(i_element) % energy(1))) - energy_max(PHOTON) = min(energy_max(PHOTON), & - exp(elements(i_element) % energy(size(elements(i_element) & - % energy)))) - call set_particle_energy_bounds(PHOTON, energy_min(PHOTON), & - energy_max(PHOTON)) - end if - ! Add element to set call element_already_read % add(element) end if @@ -2133,38 +2119,9 @@ contains call materials(i) % set_fissionable(.true.) end if end do - - ! Generate material bremsstrahlung data for electrons and positrons - if (photon_transport .and. electron_treatment == ELECTRON_TTB) then - call material_init_bremsstrahlung(materials(i) % ptr) - end if end do - if (photon_transport .and. electron_treatment == ELECTRON_TTB) then - ! Deallocate element bremsstrahlung DCS and stopping power data since - ! only the material bremsstrahlung data is needed - do i = 1, size(elements) - if (allocated(elements(i) % stopping_power_collision)) & - deallocate(elements(i) % stopping_power_collision) - if (allocated(elements(i) % stopping_power_radiative)) & - deallocate(elements(i) % stopping_power_radiative) - if (allocated(elements(i) % dcs)) deallocate(elements(i) % dcs) - if (allocated(ttb_k_grid)) deallocate(ttb_k_grid) - end do - - ! Determine if minimum/maximum energy for bremsstrahlung is greater/less - ! than the current minimum/maximum - if (size(ttb_e_grid) >= 1) then - energy_min(PHOTON) = max(energy_min(PHOTON), ttb_e_grid(1)) - energy_max(PHOTON) = min(energy_max(PHOTON), ttb_e_grid(size(ttb_e_grid))) - call set_particle_energy_bounds(PHOTON, energy_min(PHOTON), & - energy_max(PHOTON)) - end if - - ! Take logarithm of energies since they are log-log interpolated - ttb_e_grid = log(ttb_e_grid) - call set_log_ttb_e_grid() - end if + call read_ce_cross_sections_c() ! Set up logarithmic grid for nuclides do i = 1, size(nuclides) diff --git a/src/material.cpp b/src/material.cpp index acbfc08f54..53dc6779cc 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1,5 +1,6 @@ #include "openmc/material.h" +#include // for min, max #include #include #include @@ -13,6 +14,7 @@ #include "openmc/nuclide.h" #include "openmc/photon.h" #include "openmc/search.h" +#include "openmc/settings.h" #include "openmc/xml_interface.h" namespace openmc { @@ -189,11 +191,11 @@ void Material::init_bremsstrahlung() // Integrate the PDF using cubic spline integration over the incident // particle energy if (n > 2) { - spline_c(n, &data::ttb_e_grid(i), &f(i), &z(i)); + spline(n, &data::ttb_e_grid(i), &f(i), &z(i)); double c = 0.0; for (int j = i; j < n_e - 1; ++j) { - c += spline_integrate_c(n, &data::ttb_e_grid(i), &f(i), &z(i), + c += spline_integrate(n, &data::ttb_e_grid(i), &f(i), &z(i), data::ttb_e_grid(j), data::ttb_e_grid(j+1)); ttb->pdf(j+1,i) = c; @@ -266,6 +268,31 @@ read_materials(pugi::xml_node* node) } } +extern "C" void read_ce_cross_sections_c() +{ + for (auto& mat : model::materials) { + // Generate material bremsstrahlung data for electrons and positrons + if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) { + mat->init_bremsstrahlung(); + } + } + + if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) { + // Determine if minimum/maximum energy for bremsstrahlung is greater/less + // than the current minimum/maximum + if (data::ttb_e_grid.size() >= 1) { + // TODO: off-by-one + int photon = static_cast(ParticleType::photon) - 1; + int n_e = data::ttb_e_grid.size(); + data::energy_min[photon] = std::max(data::energy_min[photon], data::ttb_e_grid(1)); + data::energy_max[photon] = std::min(data::energy_max[photon], data::ttb_e_grid(n_e - 1)); + } + + // Take logarithm of energies since they are log-log interpolated + data::ttb_e_grid = xt::log(data::ttb_e_grid); + } +} + //============================================================================== // C API //============================================================================== diff --git a/src/material_header.F90 b/src/material_header.F90 index ee3f1f8bf4..0a14096123 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -5,7 +5,6 @@ module material_header use constants use dict_header, only: DictIntInt use error - use math, only: spline, spline_integrate use nuclide_header use particle_header, only: Particle use photon_header @@ -17,7 +16,6 @@ module material_header implicit none private - public :: bremsstrahlung_init public :: free_memory_material public :: openmc_extend_materials public :: openmc_get_material_index @@ -837,217 +835,6 @@ contains end function openmc_material_set_densities - subroutine bremsstrahlung_init(this, i_material, particle) - class(BremsstrahlungData), intent(inout) :: this - integer, intent(in) :: i_material - integer, intent(in) :: particle - - integer :: i, j - integer :: i_k - integer :: n, n_e, n_k - real(8) :: c - real(8) :: k, k_l, k_r - real(8) :: e, e_l, e_r - real(8) :: w, w_l, w_r - real(8) :: x, x_l, x_r - real(8) :: t - real(8) :: r - real(8) :: awr - real(8) :: beta - real(8) :: Z_eq_sq - real(8) :: atom_density - real(8) :: mass_density - real(8) :: sum_density - real(8), allocatable :: stopping_power_collision(:) - real(8), allocatable :: stopping_power_radiative(:) - real(8), allocatable :: stopping_power(:) - real(8), allocatable :: dcs(:,:) - real(8), allocatable :: f(:) - real(8), allocatable :: z(:) - logical :: positron_ - type(Material), pointer :: mat - type(PhotonInteraction), pointer :: elm - - ! Get pointer to this material - mat => materials(i_material) - - ! Determine whether we are generating electron or positron data - positron_ = (particle == POSITRON) - - ! Get the size of the energy grids - n_k = size(ttb_k_grid) - n_e = size(ttb_e_grid) - - ! Allocate arrays for TTB data - allocate(this % pdf(n_e, n_e), source=ZERO) - allocate(this % cdf(n_e, n_e), source=ZERO) - allocate(this % yield(n_e)) - - ! Allocate temporary arrays - allocate(stopping_power_collision(n_e), source=ZERO) - allocate(stopping_power_radiative(n_e), source=ZERO) - allocate(stopping_power(n_e)) - allocate(dcs(n_k, n_e), source=ZERO) - allocate(f(n_e)) - allocate(z(n_e)) - - Z_eq_sq = ZERO - sum_density = ZERO - - ! Calculate the molecular DCS and the molecular total stopping power using - ! Bragg's additivity rule. - ! TODO: The collision stopping power cannot be accurately calculated using - ! Bragg's additivity rule since the mean excitation energies and the - ! density effect corrections cannot simply be summed together. Bragg's - ! additivity rule fails especially when a higher-density compound is - ! composed of elements that are in lower-density form at normal temperature - ! and pressure (at which the NIST stopping powers are given). It will be - ! used to approximate the collision stopping powers for now, but should be - ! fixed in the future. - do i = 1, mat % n_nuclides - ! Get pointer to current element - elm => elements(mat % element(i)) - - awr = nuclides(mat % nuclide(i)) % awr - - ! Get atomic density and mass density of nuclide given atom percent - if (mat % atom_density(1) > ZERO) then - atom_density = mat % atom_density(i) - mass_density = mat % atom_density(i) * awr - ! Given weight percent - else - atom_density = -mat % atom_density(i) / awr - mass_density = -mat % atom_density(i) - end if - - ! Calculate the "equivalent" atomic number Zeq of the material - Z_eq_sq = Z_eq_sq + atom_density * elm % Z**2 - sum_density = sum_density + atom_density - - ! Accumulate material DCS - dcs = dcs + atom_density * elm % Z**2 * elm % dcs - - ! Accumulate material collision stopping power - stopping_power_collision = stopping_power_collision + mass_density & - * MASS_NEUTRON / N_AVOGADRO * elm % stopping_power_collision - - ! Accumulate material radiative stopping power - stopping_power_radiative = stopping_power_radiative + mass_density & - * MASS_NEUTRON / N_AVOGADRO * elm % stopping_power_radiative - end do - Z_eq_sq = Z_eq_sq / sum_density - - ! Calculate the positron DCS and radiative stopping power. These are - ! obtained by multiplying the electron DCS and radiative stopping powers by - ! a factor r, which is a numerical approximation of the ratio of the - ! radiative stopping powers for positrons and electrons. Source: F. Salvat, - ! J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for - ! Monte Carlo Simulation of Electron and Photon Transport," OECD-NEA, - ! Issy-les-Moulineaux, France (2011). - if (positron_) then - do i = 1, n_e - t = log(ONE + 1.0e6_8*ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON_EV)) - r = ONE - exp(-1.2359e-1_8*t + 6.1274e-2_8*t**2 - 3.1516e-2_8*t**3 + & - 7.7446e-3_8*t**4 - 1.0595e-3_8*t**5 + 7.0568e-5_8*t**6 - & - 1.808e-6_8*t**7) - stopping_power_radiative(i) = r*stopping_power_radiative(i) - dcs(:,i) = r*dcs(:,i) - end do - end if - - ! Total material stopping power - stopping_power = stopping_power_collision + stopping_power_radiative - - ! Loop over photon energies - do i = 1, n_e - 1 - w = ttb_e_grid(i) - - ! Loop over incident particle energies - do j = i, n_e - e = ttb_e_grid(j) - - ! Reduced photon energy - k = w / e - - ! Find the lower bounding index of the reduced photon energy - i_k = binary_search(ttb_k_grid, n_k, k) - - ! Get the interpolation bounds - k_l = ttb_k_grid(i_k) - k_r = ttb_k_grid(i_k+1) - x_l = dcs(i_k, j) - x_r = dcs(i_k+1, j) - - ! Find the value of the DCS using linear interpolation in reduced - ! photon energy k - x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l) - - ! Ratio of the velocity of the charged particle to the speed of light - beta = sqrt(e*(e + TWO*MASS_ELECTRON_EV)) / (e + MASS_ELECTRON_EV) - - ! Compute the integrand of the PDF - f(j) = x / (beta**2 * stopping_power(j) * w) - end do - - ! Number of points to integrate - n = n_e - i + 1 - - ! Integrate the PDF using cubic spline integration over the incident - ! particle energy - if (n > 2) then - call spline(n, ttb_e_grid(i:), f(i:), z(i:)) - - c = ZERO - do j = i, n_e - 1 - c = c + spline_integrate(n, ttb_e_grid(i:), f(i:), z(i:), & - ttb_e_grid(j), ttb_e_grid(j+1)) - this % pdf(i,j+1) = c - end do - - ! Integrate the last two points using trapezoidal rule in log-log space - else - e_l = log(ttb_e_grid(i)) - e_r = log(ttb_e_grid(i+1)) - x_l = log(f(i)) - x_r = log(f(i+1)) - - this % pdf(i,i+1) = HALF * (e_r - e_l) * (exp(e_l + x_l) + exp(e_r + x_r)) - end if - end do - - ! Loop over incident particle energies - do j = 2, n_e - ! Set last element of PDF to small non-zero value to enable log-log - ! interpolation - this % pdf(j,j) = exp(-500.0_8) - - ! Loop over photon energies - c = ZERO - do i = 1, j - 1 - ! Integrate the CDF from the PDF using the trapezoidal rule in log-log - ! space - w_l = log(ttb_e_grid(i)) - w_r = log(ttb_e_grid(i+1)) - x_l = log(this % pdf(i,j)) - x_r = log(this % pdf(i+1,j)) - - c = c + HALF * (w_r - w_l) * (exp(w_l + x_l) + exp(w_r + x_r)) - this % cdf(i+1,j) = c - end do - - ! Set photon number yield - this % yield(j) = c - end do - - ! Use logarithm of number yield since it is log-log interpolated - where (this % yield > ZERO) - this % yield = log(this % yield) - elsewhere - this % yield = -500.0_8 - end where - - end subroutine bremsstrahlung_init - !=============================================================================== ! Fortran compatibility !=============================================================================== diff --git a/src/math.F90 b/src/math.F90 index 5a1076d9d7..5607a95561 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -2,18 +2,12 @@ module math use, intrinsic :: ISO_C_BINDING - use constants - use random_lcg, only: prn - implicit none private public :: t_percentile public :: calc_pn public :: calc_rn public :: rotate_angle - public :: spline - public :: spline_interpolate - public :: spline_integrate interface @@ -49,40 +43,6 @@ module math real(C_DOUBLE), value, intent(in) :: mu real(C_DOUBLE), optional, intent(in) :: phi end subroutine rotate_angle_c_intfc - - subroutine spline(n, x, y, z) bind(C, name='spline_c') - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: n - real(C_DOUBLE), intent(in) :: x(n) - real(C_DOUBLE), intent(in) :: y(n) - real(C_DOUBLE), intent(in) :: z(n) - end subroutine spline - - function spline_interpolate(n, x, y, z, xint) & - bind(C, name='spline_interpolate_c') result(yint) - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: n - real(C_DOUBLE), intent(in) :: x(n) - real(C_DOUBLE), intent(in) :: y(n) - real(C_DOUBLE), intent(in) :: z(n) - real(C_DOUBLE), value, intent(in) :: xint - real(C_DOUBLE) :: yint - end function spline_interpolate - - function spline_integrate(n, x, y, z, xa, xb) & - bind(C, name='spline_integrate_c') result(s) - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: n - real(C_DOUBLE), intent(in) :: x(n) - real(C_DOUBLE), intent(in) :: y(n) - real(C_DOUBLE), intent(in) :: z(n) - real(C_DOUBLE), value, intent(in) :: xa - real(C_DOUBLE), value, intent(in) :: xb - real(C_DOUBLE) :: s - end function spline_integrate end interface contains diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 041b6ff57a..fe2c532307 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -681,7 +681,7 @@ double normal_variate(double mean, double standard_deviation) { while ( true ) { double v1 = 2 * prn() - 1.; double v2 = 2 * prn() - 1.; - + double r = std::pow(v1, 2) + std::pow(v2, 2); double r2 = std::pow(r, 2); if (r2 < 1) { @@ -697,7 +697,7 @@ double muir_spectrum(double e0, double m_rat, double kt) { // 8 in https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS double sigma = std::sqrt(2.*e0*kt/m_rat); return normal_variate(e0, sigma); -} +} double watt_spectrum(double a, double b) { @@ -753,7 +753,7 @@ void broaden_wmp_polynomials(double E, double dopp, int n, double factors[]) } -void spline_c(int n, const double x[], const double y[], double z[]) +void spline(int n, const double x[], const double y[], double z[]) { double c_new[n-1]; @@ -780,8 +780,8 @@ void spline_c(int n, const double x[], const double y[], double z[]) } -double spline_interpolate_c(int n, const double x[], const double y[], - const double z[], double xint) +double spline_interpolate(int n, const double x[], const double y[], + const double z[], double xint) { // Find the lower bounding index in x of xint int i = n - 1; @@ -801,8 +801,8 @@ double spline_interpolate_c(int n, const double x[], const double y[], } -double spline_integrate_c(int n, const double x[], const double y[], - const double z[], double xa, double xb) +double spline_integrate(int n, const double x[], const double y[], + const double z[], double xa, double xb) { // Find the lower bounding index in x of the lower limit of integration. int ia = n - 1; diff --git a/src/photon.cpp b/src/photon.cpp index 3b489e2995..e9c81af120 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -3,6 +3,7 @@ #include "openmc/bremsstrahlung.h" #include "openmc/constants.h" #include "openmc/hdf5_interface.h" +#include "openmc/nuclide.h" #include "openmc/particle.h" #include "openmc/random_lcg.h" #include "openmc/search.h" @@ -772,9 +773,22 @@ std::pair klein_nishina(double alpha) // Fortran compatibility //============================================================================== -extern "C" void photon_from_hdf5_c(hid_t group) +extern "C" void photon_from_hdf5(hid_t group) { data::elements.emplace_back(group, data::elements.size()); + + // Determine if minimum/maximum energy for this element is greater/less than + // the previous + const auto& element {data::elements.back()}; + if (element.energy_.size() >= 1) { + // TODO: off-by-one + int photon = static_cast(ParticleType::photon) - 1; + int n = element.energy_.size(); + data::energy_min[photon] = std::max(data::energy_min[photon], + std::exp(element.energy_(1))); + data::energy_max[photon] = std::min(data::energy_max[photon], + std::exp(element.energy_(n - 1))); + } } extern "C" void photon_calculate_xs(int i_element, double E) @@ -782,4 +796,12 @@ extern "C" void photon_calculate_xs(int i_element, double E) data::elements[i_element - 1].calculate_xs(E); } +extern "C" void free_memory_photon_c() +{ + data::elements.clear(); + data::compton_profile_pz.resize({0}); + data::ttb_e_grid.resize({0}); + data::ttb_k_grid.resize({0}); +} + } // namespace openmc diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 5a887d259a..d077b72a29 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -2,53 +2,13 @@ module photon_header use, intrinsic :: ISO_C_BINDING - use algorithm, only: binary_search use constants use dict_header, only: DictCharInt - use hdf5_interface - use settings - real(8), allocatable :: compton_profile_pz(:) - real(8), allocatable :: ttb_e_grid(:) ! energy T of incident electron - real(8), allocatable :: ttb_k_grid(:) ! reduced energy W/T of emitted photon - - type PhotonInteraction - character(3) :: name ! atomic symbol, e.g. 'Zr' - integer :: Z ! atomic number - - ! Microscopic cross sections - real(8), allocatable :: energy(:) - - ! Stopping power data - real(8) :: I ! mean excitation energy - real(8), allocatable :: stopping_power_collision(:) - real(8), allocatable :: stopping_power_radiative(:) - - ! Bremsstrahlung scaled DCS - real(8), allocatable :: dcs(:,:) - - contains - procedure :: from_hdf5 => photon_from_hdf5 - end type PhotonInteraction - - type BremsstrahlungData - real(8), allocatable :: pdf(:,:) ! Bremsstrahlung energy PDF - real(8), allocatable :: cdf(:,:) ! Bremsstrahlung energy CDF - real(8), allocatable :: yield(:) ! Photon number yield - end type BremsstrahlungData - - type Bremsstrahlung - type(BremsstrahlungData) :: electron - type(BremsstrahlungData) :: positron - end type Bremsstrahlung - - type(PhotonInteraction), allocatable, target :: elements(:) ! Photon cross sections integer :: n_elements ! Number of photon cross section tables type(DictCharInt) :: element_dict - type(Bremsstrahlung), allocatable, target :: ttb(:) ! Bremsstrahlung data - !=============================================================================== ! ELEMENTMICROXS contains cached microscopic photon cross sections for a ! particular element at the current energy @@ -70,142 +30,19 @@ module photon_header contains - subroutine photon_from_hdf5(this, group_id) - class(PhotonInteraction), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer :: i - integer(HID_T) :: rgroup - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(1), dims2(2) - integer :: n_energy - integer :: n_k - integer :: n_e - real(8) :: f - real(8) :: y - real(8), allocatable :: electron_energy(:) - 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) - - ! Get rid of leading '/' - this % name = trim(this % name(2:)) - - ! Get atomic number - call read_attribute(this % Z, group_id, 'Z') - - ! Determine number of energies and read energy grid - dset_id = open_dataset(group_id, 'energy') - call get_shape(dset_id, dims) - n_energy = int(dims(1), 4) - allocate(this % energy(dims(1))) - call read_dataset(this % energy, dset_id) - call close_dataset(dset_id) - - if (electron_treatment == ELECTRON_TTB) then - ! Read bremsstrahlung scaled DCS - rgroup = open_group(group_id, 'bremsstrahlung') - dset_id = open_dataset(rgroup, 'dcs') - call get_shape(dset_id, dims2) - n_k = int(dims2(1), 4) - n_e = int(dims2(2), 4) - allocate(this % dcs(n_k, n_e)) - call read_dataset(this % dcs, dset_id) - call close_dataset(dset_id) - - ! Get energy grids used for bremsstrahlung DCS and for stopping powers - allocate(electron_energy(n_e)) - call read_dataset(electron_energy, rgroup, 'electron_energy') - if (.not. allocated(ttb_k_grid)) then - allocate(ttb_k_grid(n_k)) - call read_dataset(ttb_k_grid, rgroup, 'photon_energy') - end if - call close_group(rgroup) - - ! Read stopping power data - if (this % Z < 99) then - rgroup = open_group(group_id, 'stopping_powers') - allocate(this % stopping_power_collision(n_e)) - allocate(this % stopping_power_radiative(n_e)) - call read_dataset(this % stopping_power_collision, rgroup, 's_collision') - call read_dataset(this % stopping_power_radiative, rgroup, 's_radiative') - call read_attribute(this % I, rgroup, 'I') - call close_group(rgroup) - end if - - ! Truncate the bremsstrahlung data at the cutoff energy - if (energy_cutoff(PHOTON) > electron_energy(1)) then - i_grid = binary_search(electron_energy, n_e, energy_cutoff(PHOTON)) - - ! calculate interpolation factor - f = (log(energy_cutoff(PHOTON)) - log(electron_energy(i_grid))) / & - (log(electron_energy(i_grid+1)) - log(electron_energy(i_grid))) - - ! Interpolate collision stopping power at the cutoff energy and - ! truncate - y = exp(log(this % stopping_power_collision(i_grid)) + & - f*(log(this % stopping_power_collision(i_grid+1)) - & - log(this % stopping_power_collision(i_grid)))) - this % stopping_power_collision = & - [y, this % stopping_power_collision(i_grid+1:n_e)] - - ! Interpolate radiative stopping power at the cutoff energy and - ! truncate - y = exp(log(this % stopping_power_radiative(i_grid)) + & - f*(log(this % stopping_power_radiative(i_grid+1)) - & - log(this % stopping_power_radiative(i_grid)))) - this % stopping_power_radiative = & - [y, this % stopping_power_radiative(i_grid+1:n_e)] - - ! Interpolate bremsstrahlung DCS at the cutoff energy and truncate - allocate(dcs(n_k, n_e-i_grid+1)) - do i = 1, n_k - y = exp(log(this % dcs(i,i_grid)) + & - f*(log(this % dcs(i,i_grid+1)) - log(this % dcs(i,i_grid)))) - dcs(i,:) = [y, this % dcs(i,i_grid+1:n_e)] - end do - call move_alloc(dcs, this % dcs) - - electron_energy = [energy_cutoff(PHOTON), electron_energy(i_grid+1:n_e)] - end if - - ! Set incident particle energy grid - if (.not. allocated(ttb_e_grid)) then - call move_alloc(electron_energy, ttb_e_grid) - end if - end if - - ! Take logarithm of energies and cross sections since they are log-log - ! interpolated - this % energy = log(this % energy) - - end subroutine photon_from_hdf5 - !=============================================================================== ! FREE_MEMORY_PHOTON deallocates/resets global variables in this module !=============================================================================== subroutine free_memory_photon() + interface + subroutine free_memory_photon_c() bind(C) + end subroutine + end interface + ! Deallocate photon cross section data - if (allocated(elements)) deallocate(elements) - if (allocated(compton_profile_pz)) deallocate(compton_profile_pz) n_elements = 0 call element_dict % clear() - - ! Clear TTB-related arrays - if (allocated(ttb_e_grid)) deallocate(ttb_e_grid) - if (allocated(ttb)) deallocate(ttb) end subroutine free_memory_photon function micro_photon_xs_ptr() result(ptr) bind(C) diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 deleted file mode 100644 index 1575450747..0000000000 --- a/src/photon_physics.F90 +++ /dev/null @@ -1,116 +0,0 @@ -module photon_physics - - use algorithm, only: binary_search - use constants - use particle_header - use photon_header, only: BremsstrahlungData, ttb_e_grid, ttb - use random_lcg, only: prn - use settings - -contains - -!=============================================================================== -! THICK_TARGET_BREMSSTRAHLUNG -!=============================================================================== - - subroutine thick_target_bremsstrahlung(p, E_lost) bind(C) - type(Particle), intent(inout) :: p - real(C_DOUBLE), intent(out) :: E_lost - - integer :: i, j - integer :: i_e, i_w - integer :: n - integer :: n_e - real(8) :: a - real(8) :: f - real(8) :: e, e_l, e_r - real(8) :: y, y_l, y_r - real(8) :: w, w_l, w_r - real(8) :: p_l, p_r - real(8) :: c, c_l, c_max - type(BremsstrahlungData), pointer :: mat - - if (p % material == MATERIAL_VOID) return - - if (p % E < energy_cutoff(PHOTON)) return - - ! Get bremsstrahlung data for this material and particle type - if (p % type == POSITRON) then - mat => ttb(p % material) % positron - else - mat => ttb(p % material) % electron - end if - - e = log(p % E) - n_e = size(ttb_e_grid) - - ! Find the lower bounding index of the incident electron energy - j = binary_search(ttb_e_grid, n_e, e) - if (j == n_e) j = j - 1 - - ! Get the interpolation bounds - e_l = ttb_e_grid(j) - e_r = ttb_e_grid(j+1) - y_l = mat % yield(j) - y_r = mat % yield(j+1) - - ! Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF - ! interpolated in log energy, which can be interpreted as the probability - ! of index j+1 - f = (e - e_l)/(e_r - e_l) - - ! Get the photon number yield for the given energy using linear - ! interpolation on a log-log scale - y = exp(y_l + (y_r - y_l)*f) - - ! Sample number of secondary bremsstrahlung photons - n = int(y + prn()) - - E_lost = ZERO - if (n == 0) return - - ! Sample index of the tabulated PDF in the energy grid, j or j+1 - if (prn() <= f .or. j == 1) then - i_e = j + 1 - - ! Interpolate the maximum value of the CDF at the incoming particle - ! energy on a log-log scale - p_l = mat % pdf(i_e-1, i_e) - p_r = mat % pdf(i_e, i_e) - c_l = mat % cdf(i_e-1, i_e) - a = log(p_r/p_l)/(e_r - e_l) + ONE - c_max = c_l + exp(e_l)*p_l/a*(exp(a*(e - e_l)) - ONE) - else - i_e = j - - ! Maximum value of the CDF - c_max = mat % cdf(i_e, i_e) - end if - - ! Sample the energies of the emitted photons - do i = 1, n - ! Generate a random number r and determine the index i for which - ! cdf(i) <= r*cdf,max <= cdf(i+1) - c = prn()*c_max - i_w = binary_search(mat % cdf(:i_e,i_e), i_e, c) - - ! Sample the photon energy - w_l = ttb_e_grid(i_w) - w_r = ttb_e_grid(i_w+1) - p_l = mat % pdf(i_w, i_e) - p_r = mat % pdf(i_w+1, i_e) - c_l = mat % cdf(i_w, i_e) - a = log(p_r/p_l)/(w_r - w_l) + ONE - w = exp(w_l)*(a*(c - c_l)/(exp(w_l)*p_l) + ONE)**(ONE/a) - - if (w > energy_cutoff(PHOTON)) then - ! Create secondary photon - call particle_create_secondary(p, p % coord(1) % uvw, w, PHOTON, & - run_ce=.true._C_BOOL) - E_lost = E_lost + w - end if - end do - - end subroutine thick_target_bremsstrahlung - -end module photon_physics diff --git a/src/settings.F90 b/src/settings.F90 index c26e9666be..f0b6a131fd 100644 --- a/src/settings.F90 +++ b/src/settings.F90 @@ -23,7 +23,6 @@ module settings integer(C_INT), bind(C) :: n_log_bins ! number of bins for logarithmic grid logical(C_BOOL), bind(C) :: photon_transport - integer(C_INT), bind(C) :: electron_treatment ! ============================================================================ ! MULTI-GROUP CROSS SECTION RELATED VARIABLES From 1138b5bb45531dea584f511eab992d3dbcc1cfe0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 15 Jan 2019 07:25:38 -0600 Subject: [PATCH 10/14] Convert sample_element to C++ and remove physics.F90 --- CMakeLists.txt | 1 - include/openmc/physics.h | 2 +- src/physics.F90 | 69 ---------------------------------------- src/physics.cpp | 61 ++++++++++++++++++++--------------- src/tracking.F90 | 7 ++-- 5 files changed, 42 insertions(+), 98 deletions(-) delete mode 100644 src/physics.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 5bf510e084..f2b54c578f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -332,7 +332,6 @@ add_library(libopenmc SHARED src/particle_restart.F90 src/photon_header.F90 src/physics_common.F90 - src/physics.F90 src/pugixml/pugixml_f.F90 src/random_lcg.F90 src/reaction_header.F90 diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 921494db7a..650f65feaf 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -44,7 +44,7 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat); void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_array, int64_t* bank_size, int64_t bank_capacity); -extern "C" int sample_element(Particle* p); +int sample_element(Particle* p); Reaction* sample_fission(int i_nuclide, double E); diff --git a/src/physics.F90 b/src/physics.F90 deleted file mode 100644 index c6749a6a1f..0000000000 --- a/src/physics.F90 +++ /dev/null @@ -1,69 +0,0 @@ -module physics - - use constants - use error, only: fatal_error - use material_header, only: Material, materials - use math - use message_passing - use nuclide_header - use particle_header - use photon_header - use random_lcg, only: prn - use settings - use simulation_header - - implicit none - - interface - subroutine collision(p) bind(C) - import Particle - type(Particle), intent(inout) :: p - end subroutine - end interface - -contains - -!=============================================================================== -! SAMPLE_ELEMENT -!=============================================================================== - - function sample_element(p) result(i_element) bind(C) - type(Particle), intent(in) :: p - integer(C_INT) :: i_element - - integer :: i - 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 - - associate (mat => materials(p % material)) - ! Sample cumulative distribution function - cutoff = prn() * material_xs % total - - i = 0 - prob = ZERO - do while (prob < cutoff) - i = i + 1 - - ! Check to make sure that a nuclide was sampled - if (i > mat % n_nuclides) then - call particle_write_restart(p) - call fatal_error("Did not sample any element during collision.") - end if - - ! 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 - - end function sample_element - -end module physics diff --git a/src/physics.cpp b/src/physics.cpp index 564bc0f5cb..817577cab9 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -229,10 +229,10 @@ void sample_photon_reaction(Particle* p) // Sample element within material int i_element = sample_element(p); - p->event_nuclide = i_element; // TODO: off-by-one - const auto& micro {simulation::micro_photon_xs[i_element - 1]}; - const auto& element {data::elements[i_element - 1]}; + p->event_nuclide = i_element + 1; + const auto& micro {simulation::micro_photon_xs[i_element]}; + const auto& element {data::elements[i_element]}; // Calculate photon energy over electron rest mass equivalent double alpha = p->E/MASS_ELECTRON_EV; @@ -479,32 +479,43 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat) } } -// int sample_element(Particle* p) -// { -// // Sample cumulative distribution function -// double cutoff = prn() * simulation::material_xs.total; +int sample_element(Particle* p) +{ + // Sample cumulative distribution function + double cutoff = prn() * simulation::material_xs.total; -// int i = 0; -// double prob = 0.0; -// while (prob < cutoff) { -// // Check to make sure that a nuclide was sampled -// if (i > mat % n_nuclides) { -// p->write_restart(); -// fatal_error("Did not sample any element during collision."); -// } + // Get pointers to elements, densities + int* nuclide; + double* density; + int n; + openmc_material_get_densities(p->material, &nuclide, &density, &n); + int* element = material_element(p->material); -// // Find atom density -// int i_element = mat % element(i); -// double atom_density = mat % atom_density(i); + int i = 0; + double prob = 0.0; + int i_element; + while (prob < cutoff) { + // Check to make sure that a nuclide was sampled + if (i >= n) { + p->write_restart(); + fatal_error("Did not sample any element during collision."); + } -// // Determine microscopic cross section -// double sigma = atom_density * simulation::micro_photon_xs[i_element].total; + // Find atom density + // TODO: off-by-one + i_element = element[i] - 1; + double atom_density = density[i]; -// // Increment probability to compare to cutoff -// prob += sigma; -// ++i; -// } -// } + // Determine microscopic cross section + double sigma = atom_density * simulation::micro_photon_xs[i_element].total; + + // Increment probability to compare to cutoff + prob += sigma; + ++i; + } + + return i_element; +} Reaction* sample_fission(int i_nuclide, double E) { diff --git a/src/tracking.F90 b/src/tracking.F90 index b9da99f127..bf8332852d 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -16,7 +16,6 @@ module tracking use mgxs_interface use nuclide_header use particle_header - use physics, only: collision use random_lcg, only: prn, prn_set_stream use settings use simulation_header @@ -33,11 +32,15 @@ module tracking implicit none interface + subroutine collision(p) bind(C) + import Particle + type(Particle), intent(inout) :: p + end subroutine + subroutine collision_mg(p) bind(C) import Particle, C_DOUBLE type(Particle), intent(inout) :: p end subroutine collision_mg - end interface contains From 6d4017f3a27f1d2cc3edc8d09386edebe9cf8060 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 15 Jan 2019 07:37:32 -0600 Subject: [PATCH 11/14] Remove physics_common.F90 --- CMakeLists.txt | 1 - src/physics_common.F90 | 20 -------------------- 2 files changed, 21 deletions(-) delete mode 100644 src/physics_common.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index f2b54c578f..e33979d2c8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -331,7 +331,6 @@ add_library(libopenmc SHARED src/particle_header.F90 src/particle_restart.F90 src/photon_header.F90 - src/physics_common.F90 src/pugixml/pugixml_f.F90 src/random_lcg.F90 src/reaction_header.F90 diff --git a/src/physics_common.F90 b/src/physics_common.F90 deleted file mode 100644 index 0d5b9fb324..0000000000 --- a/src/physics_common.F90 +++ /dev/null @@ -1,20 +0,0 @@ -module physics_common - - use, intrinsic :: ISO_C_BINDING - - use particle_header, only: Particle - - implicit none - -!=============================================================================== -! RUSSIAN_ROULETTE FROM C -!=============================================================================== - - interface - subroutine russian_roulette(p) bind(C) - import Particle - type(Particle), intent(inout) :: p - end subroutine russian_roulette - end interface - -end module physics_common From 0b6d25f75846fe34d3a05988d475372a7a194770 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 17 Jan 2019 15:07:04 -0600 Subject: [PATCH 12/14] Try using xvfb as a service (to fix Travis CI failures) --- .travis.yml | 3 ++- tools/ci/travis-before-script.sh | 4 ---- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/.travis.yml b/.travis.yml index 2c3a50559a..1b6ce78eb4 100644 --- a/.travis.yml +++ b/.travis.yml @@ -17,6 +17,8 @@ addons: - liblapack-dev config: retries: true +services: + - xvfb cache: directories: - $HOME/nndc_hdf5 @@ -31,7 +33,6 @@ env: - OPENMC_ENDF_DATA=$HOME/endf-b-vii.1 - LD_LIBRARY_PATH=$HOME/MOAB/lib:$HOME/DAGMC/lib - PATH=$PATH:$HOME/NJOY2016/build - - DISPLAY=:99.0 - COVERALLS_PARALLEL=true matrix: - OMP=n MPI=n PHDF5=n diff --git a/tools/ci/travis-before-script.sh b/tools/ci/travis-before-script.sh index 9af8212fa9..61f0bc0877 100755 --- a/tools/ci/travis-before-script.sh +++ b/tools/ci/travis-before-script.sh @@ -1,9 +1,5 @@ #!/bin/bash set -ex -# Allow tests that require GUI as described at: -# https://docs.travis-ci.com/user/gui-and-headless-browsers/#Using-xvfb-to-Run-Tests-That-Require-a-GUI -sh -e /etc/init.d/xvfb start - # Download NNDC HDF5 data, ENDF/B-VII.1 distribution, multipole library source tools/ci/download-xs.sh From 450359ac872d8c58e06a14a3bd0cd5ad1ec508b5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 18 Jan 2019 08:25:58 -0600 Subject: [PATCH 13/14] Make sure numpy gets upgraded on Travis --- tools/ci/travis-install.sh | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tools/ci/travis-install.sh b/tools/ci/travis-install.sh index 670f28f33b..10c168d129 100755 --- a/tools/ci/travis-install.sh +++ b/tools/ci/travis-install.sh @@ -9,11 +9,10 @@ if [[ $DAGMC = 'y' ]]; then ./tools/ci/travis-install-dagmc.sh fi -# Upgrade pip before doing anything else +# Upgrade pip, pytest, numpy before doing anything else pip install --upgrade pip - -# pytest installed by default -- make sure we get latest pip install --upgrade pytest +pip install --upgrade numpy # Install mpi4py for MPI configurations if [[ $MPI == 'y' ]]; then From 7650ea42d2bd253a20ae05c07eb64ce085836a20 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 22 Jan 2019 06:59:46 -0600 Subject: [PATCH 14/14] Fix threshold index check for photoelectric cross sections --- src/physics.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics.cpp b/src/physics.cpp index 817577cab9..90f0708f3f 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -303,7 +303,7 @@ void sample_photon_reaction(Particle* p) // Check threshold of reaction int i_start = shell.threshold; - if (i_grid <= i_start) continue; + if (i_grid < i_start) continue; // Evaluation subshell photoionization cross section double xs = std::exp(shell.cross_section(i_grid - i_start) +