From 9f8ebf527b3d1f752e6c503165ee718089469816 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 15 Jan 2019 06:51:18 -0600 Subject: [PATCH] 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