From 62985f61987b8c0e2d8cff6ebd2a318a590ce4f7 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 11 Jan 2018 21:43:21 -0600 Subject: [PATCH] Added cubic spline interpolation/integration; other fixes --- openmc/data/photon.py | 35 +++++----- src/math.F90 | 145 +++++++++++++++++++++++++++++++++++++++++ src/photon_header.F90 | 80 ++++++++++++++--------- src/photon_physics.F90 | 12 ++-- src/physics.F90 | 9 ++- 5 files changed, 225 insertions(+), 56 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 181f2f1fe6..d17aaccf7b 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -11,7 +11,7 @@ import pandas as pd from openmc.mixin import EqualityMixin import openmc.checkvalue as cv from . import HDF5_VERSION -from .data import ATOMIC_SYMBOL +from .data import ATOMIC_SYMBOL, EV_PER_MEV from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record from .function import Tabulated1D @@ -297,18 +297,18 @@ class IncidentPhoton(EqualityMixin): :math:`p_z` for each subshell). Note that subshell occupancies may not match the atomic relaxation data. stopping_powers : dict - Dictionary of stopping power data with keys 'energy', 'density' (mass - density in g/cm:sup:`3`), 'I' (mean excitation energy), 's_collision' - (collision stopping power in MeV cm:sup:`2`/g), 's_radiative' - (radiative stopping power in MeV cm:sup:`2`/g), and 'density_effect' - (density effect parameter). + Dictionary of stopping power data with keys 'energy' (in eV), 'density' + (mass density in g/cm:sup:`3`), 'I' (mean excitation energy), + 's_collision' (collision stopping power in eV cm:sup:`2`/g), + 's_radiative' (radiative stopping power in eV cm:sup:`2`/g), and + 'density_effect' (density effect parameter). bremsstrahlung : dict Dictionary of bremsstrahlung DCS data with keys 'electron_energy' - (incident electron kinetic energy values in MeV), 'photon_energy' + (incident electron kinetic energy values in eV), 'photon_energy' (ratio of the energy of the emitted photon to the incident electron - kinetic energy), and 'dcs' (cross sectin values). The cross sections - are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`, where - :math:`E_k` is the energy of the emitted photon. + kinetic energy), and 'dcs' (cross sectin values in mb). The cross + sections are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`, + where :math:`E_k` is the energy of the emitted photon. reactions : collections.OrderedDict Contains the cross sections for each photon reaction. The keys are MT values and the values are instances of :class:`PhotonReaction`. @@ -429,7 +429,8 @@ class IncidentPhoton(EqualityMixin): if not _STOPPING_POWERS: filename = os.path.join(os.path.dirname(__file__), 'stopping_powers.h5') with h5py.File(filename, 'r') as f: - _STOPPING_POWERS['energy'] = f['energy'].value + # Units are in MeV; convert to eV + _STOPPING_POWERS['energy'] = f['energy'].value*EV_PER_MEV for i in range(1, 99): group = f['{:03}'.format(i)] _STOPPING_POWERS[i] = {'density': group.attrs['density'], @@ -438,6 +439,10 @@ class IncidentPhoton(EqualityMixin): 's_radiative': group['s_radiative'].value, 'density_effect': group['density_effect'].value} + # Units are in MeV cm^2/g; convert to eV cm^2/g + _STOPPING_POWERS[i]['s_collision'] *= EV_PER_MEV + _STOPPING_POWERS[i]['s_radiative'] *= EV_PER_MEV + # Add stopping power data if Z < 99: data.stopping_powers['energy'] = _STOPPING_POWERS['energy'] @@ -448,8 +453,8 @@ class IncidentPhoton(EqualityMixin): filename = os.path.join(os.path.dirname(__file__), 'BREMX.DAT') brem = open(filename, 'r').read().split() - # Incident electron kinetic energy grid - _BREMSSTRAHLUNG['electron_energy'] = np.logspace(-3, 3, 200) + # Incident electron kinetic energy grid in eV + _BREMSSTRAHLUNG['electron_energy'] = np.logspace(3, 9, 200) log_energy = np.log(_BREMSSTRAHLUNG['electron_energy']) # Get number of tabulated electron and photon energy values @@ -460,8 +465,8 @@ class IncidentPhoton(EqualityMixin): p = 39 # Get log of incident electron kinetic energy values, used for cubic - # spline interpolation in log energy - logx = np.log(np.fromiter(brem[p:p+n], float, n)) + # spline interpolation in log energy. Units are in MeV, so convert to eV. + logx = np.log(np.fromiter(brem[p:p+n], float, n)*EV_PER_MEV) p += n # Get reduced photon energy values diff --git a/src/math.F90 b/src/math.F90 index 2fa0a6ce1e..5e4c8f81b1 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -2,6 +2,7 @@ module math use, intrinsic :: ISO_C_BINDING + use algorithm, only: binary_search use constants use random_lcg, only: prn @@ -827,4 +828,148 @@ contains end do end subroutine broaden_wmp_polynomials +!=============================================================================== +! SPLINE +!=============================================================================== + + subroutine spline(x, y, z, n) + + integer, intent(in) :: n + real(8), intent(in) :: x(n) + real(8), intent(in) :: y(n) + real(8), intent(out) :: z(n) + + integer :: i + real(8) :: a, b, c, d + real(8), allocatable :: c_new(:) + + allocate(c_new(n-1)) + + ! Set natural boundary conditions + c_new(1) = ZERO + z(1) = ZERO + z(n) = ZERO + + ! Solve using tridiagonal matrix algorithm; first do forward sweep + do i = 2, n - 1 + a = x(i) - x(i-1) + c = x(i+1) - x(i) + b = TWO * (a + c) + d = 6.0_8 * ((y(i+1) - y(i)) / c - (y(i) - y(i-1)) / a) + + c_new(i) = c / (b - a * c_new(i-1)) + z(i) = (d - a * z(i-1)) / (b - a * c_new(i-1)) + end do + + ! Back substitution + do i = n - 1, 1, -1 + z(i) = z(i) - c_new(i) * z(i+1) + end do + + deallocate(c_new) + + end subroutine spline + +!=============================================================================== +! SPLINE_INTERPOLATE +!=============================================================================== + + function spline_interpolate(x, y, z, n, xint) result(yint) + + integer, intent(in) :: n + real(8), intent(in) :: x(n) + real(8), intent(in) :: y(n) + real(8), intent(in) :: z(n) + real(8), intent(in) :: xint + real(8) :: yint + + integer :: i + real(8) :: h, r + real(8) :: b, c, d + + ! Find the lower bounding index in x of xint + if (xint < x(1)) then + i = 1 + else if (xint > x(n)) then + i = n - 1 + else + i = binary_search(x, n, xint) + end if + + h = x(i+1) - x(i) + r = xint - x(i) + + ! Compute the coefficients + b = (y(i+1) - y(i)) / h - (h / 6.0_8) * (z(i+1) + TWO * z(i)) + c = HALF * z(i) + d = (z(i+1) - z(i)) / (h * 6.0_8) + + yint = y(i) + b * r + c * r**2 + d * r**3 + + end function spline_interpolate + +!=============================================================================== +! SPLINE_INTEGRATE +!=============================================================================== + + function spline_integrate(x, y, z, n, xa, xb) result(s) + + integer, intent(in) :: n + real(8), intent(in) :: x(n) + real(8), intent(in) :: y(n) + real(8), intent(in) :: z(n) + real(8), intent(in) :: xa ! Lower limit of integration + real(8), intent(in) :: xb ! Upper limit of integration + real(8) :: s + + integer :: i + integer :: ia, ib + real(8) :: h, r + real(8) :: a, b, c, d + + ! Find the lower bounding index in x of the lower limit of integration. + if (xa < x(1)) then + ia = 1 + else if (xa > x(n)) then + ia = n - 1 + else + ia = binary_search(x, n, xa) + end if + + ! Find the lower bounding index in x of the upper limit of integration. + if (xb < x(1)) then + ib = 1 + else if (xb > x(n)) then + ib = n - 1 + else + ib = binary_search(x, n, xb) + end if + + ! Evaluate the integral + s = ZERO + do i = ia, ib + h = x(i+1) - x(i) + + ! Compute the coefficients + b = (y(i+1) - y(i)) / h - (h / 6.0_8) * (z(i+1) + TWO * z(i)) + c = HALF * z(i) + d = (z(i+1) - z(i)) / (h * 6.0_8) + + ! Subtract the integral from x(ia) to xa + if (i == ia) then + r = xa - x(ia) + s = s - (y(i) * r + b * r**2 * HALF + c * r**3 / THREE + d * r**4 / FOUR) + end if + + ! Integrate from x(ib) to xb in final interval + if (i == ib) then + h = xb - x(ib) + end if + + ! Accumulate the integral + s = s + y(i) * h + b * h**2 * HALF + c * h**3 / THREE + d * h**4 / FOUR + end do + + end function spline_integrate + end module math diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 70d5b3caf5..c00044d4c8 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -2,16 +2,19 @@ module photon_header use hdf5, only: HID_T, HSIZE_T, SIZE_T - use algorithm, only: binary_search - use constants, only: ZERO, HALF, SUBSHELLS - use dict_header, only: DictIntInt, DictCharInt - use endf_header, only: Tabulated1D + use algorithm, only: binary_search + use constants, only: ZERO, HALF, SUBSHELLS + use dict_header, only: DictIntInt, DictCharInt + use endf_header, only: Tabulated1D use hdf5_interface + use math, only: spline, spline_integrate + use material_header, only: Material, materials + use nuclide_header, only: nuclides + use settings real(8), allocatable :: compton_profile_pz(:) real(8), allocatable :: ttb_energy_electron(:) ! incident electron energy grid real(8), allocatable :: ttb_energy_photon(:) ! reduced photon energy grid - real(8) :: ttb_energy_cutoff type ElectronSubshell integer :: index_subshell ! index in SUBSHELLS @@ -72,7 +75,7 @@ module photon_header type Bremsstrahlung integer :: i_material ! Index in materials array - real(8), allocatable :: yield(:,:) ! Photon number yield + real(8), allocatable :: yield(:) ! Photon number yield real(8), allocatable :: dcs(:,:) ! Scaled bremsstrahlung DCS real(8), allocatable :: cdf(:,:) ! Bremsstrahlung energy CDF @@ -80,12 +83,12 @@ module photon_header procedure :: init => bremsstrahlung_init end type Bremsstrahlung - type(PhotonInteraction), allocatable :: elements(:) ! Photon cross sections + type(PhotonInteraction), allocatable, target :: elements(:) ! Photon cross sections integer :: n_elements ! Number of photon cross section tables type(DictCharInt) :: element_dict - type(Bremsstrahlung), allocatable :: ttb(:) ! Bremsstrahlung cross sections + type(Bremsstrahlung), allocatable, target :: ttb(:) ! Bremsstrahlung cross sections !=============================================================================== ! ELEMENTMICROXS contains cached microscopic photon cross sections for a @@ -302,9 +305,6 @@ contains this % pair_production_electron if (electron_treatment == ELECTRON_TTB) then - ! TODO: read in - ttb_energy_cutoff = 1.e-3 - ! Read bremsstrahlung scaled DCS rgroup = open_group(group_id, 'bremsstrahlung') dset_id = open_dataset(rgroup, 'dcs') @@ -385,9 +385,14 @@ contains real(8) :: x_l, x_r, x_c real(8) :: Z_eq_sq real(8) :: beta + real(8) :: atom_fraction real(8), allocatable :: mass_fraction(:) real(8), allocatable :: stopping_power(:) real(8), allocatable :: mfp_inv(:) + real(8), allocatable :: x(:) + real(8), allocatable :: y(:) + real(8), allocatable :: z(:) + type(DictIntInt) :: nuc_dict type(Material), pointer :: mat type(PhotonInteraction), pointer :: elm @@ -398,24 +403,29 @@ contains ! Allocate and initialize arrays n_k = size(ttb_energy_photon) n_e = size(ttb_energy_electron) - allocate(mass_fraction(size(mat % element))) + allocate(mass_fraction(mat % n_nuclides)) allocate(stopping_power(n_e)) allocate(mfp_inv(n_e)) allocate(this % yield(n_e)) allocate(this % dcs(n_k, n_e)) allocate(this % cdf(n_k, n_e)) + allocate(x(n_e)) + allocate(y(n_e)) + allocate(z(n_e)) stopping_power(:) = ZERO mfp_inv(:) = ZERO this % dcs(:,:) = ZERO this % cdf(:,:) = ZERO - ! TODO - ! Calculate the "equivalent" atomic number Zeq and get the mass fraction of + ! Calculate the "equivalent" atomic number Zeq and the mass fraction of ! each element - Z_eq_sq = 0 + Z_eq_sq = ZERO do i = 1, mat % n_nuclides - ! Zeq**2 = (atomic fraction of x)*Zx**2 + (atomic fraction of y)*Zy**2 + ... + atom_fraction = mat % atom_density(i) / mat % density + mass_fraction(i) = atom_fraction * nuclides(mat % nuclide(i)) % awr + Z_eq_sq = Z_eq_sq + atom_fraction * nuclides(mat % nuclide(i)) % Z**2 end do + mass_fraction = mass_fraction / sum(mass_fraction) ! Calculate the molecular DCS and the molecular total stopping power using ! Bragg's additivity rule. Note: the collision stopping power cannot be @@ -426,13 +436,17 @@ contains ! 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, size(mat % element) + do i = 1, mat % n_nuclides ! Get pointer to current element - elm => mat % element(i) + elm => elements(mat % element(i)) - ! TODO: n_atoms + ! Get atomic fraction + atom_fraction = mat % atom_density(i) / mat % density + + ! TODO: for molecular DCS, atom_fraction should actually be the number of + ! atoms in the molecule. ! Accumulate material DCS - this % dcs = this % dcs + n_atoms(i) * elm % Z**2 / Z_eq_sq * elm % dcs + this % dcs = this % dcs + atom_fraction * elm % Z**2 / Z_eq_sq * elm % dcs ! Accumulate material total stopping power stopping_power = stopping_power + mass_fraction(i) * elm % density * & @@ -442,13 +456,13 @@ contains ! Calculate inverse bremsstrahlung mean free path do i = 1, n_e e = exp(ttb_energy_electron(i)) - if (e <= ttb_energy_cutoff) cycle + if (e <= energy_cutoff(PHOTON)) cycle ! Ratio of the velocity of the charged particle to the speed of light - beta = sqrt(e*(e + TWO*MASS_ELECTRON/1.e6)) / (e + MASS_ELECTRON/1.e6) + beta = sqrt(e*(e + TWO*MASS_ELECTRON/1.e6_8)) / (e + MASS_ELECTRON/1.e6_8) ! Integration lower bound - k_c = ttb_energy_cutoff / e + k_c = energy_cutoff(PHOTON) / e ! Find the upper bounding index of the reduced photon cutoff energy i_k = binary_search(ttb_energy_photon, n_k, k_c) + 1 @@ -470,18 +484,24 @@ contains c = c + HALF * (log(ttb_energy_photon(j+1)) - & log(ttb_energy_photon(j))) * (this % dcs(j,i) + this % dcs(j+1,i)) this % cdf(j+1,i) = c + end do - ! TODO: density in atom/cm^3 - mfp_inv(i) = c * mat % density * Z_eq_sq / beta**2 * 1.0e-27_8 + ! Calculate the inverse bremsstrahlung mean free path + mfp_inv(i) = c * mat % density * Z_eq_sq / beta**2 * 1.0e-3_8 end do - ! TODO: ! Calculate photon number yield - ! cs = cubic_spline(exp(ttb_energy_electron), mfp_inv / stopping_power) - ! do i = 1, n_e - ! yield(i) = cs % integrate(ttb_energy_cutoff, exp(ttb_energy_electron(i))) + x = exp(ttb_energy_electron) + y = mfp_inv / stopping_power + call spline(x, y, z, n_e) + do i = 1, n_e + this % yield(i) = spline_integrate(x, y, z, n_e, energy_cutoff(PHOTON), x(i)) + end do + ! Use logarithm of number yield since it is log-log interpolated - ! yield = log(yield) + this % yield = log(this % yield) + + deallocate(mass_fraction, stopping_power, mfp_inv, x, y, z) end subroutine bremsstrahlung_init diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 4e3d815d15..bba86fcf69 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -5,7 +5,7 @@ module photon_physics use particle_header, only: Particle use photon_header, only: PhotonInteraction, Bremsstrahlung, & compton_profile_pz, ttb_energy_electron, & - ttb_energy_photon, ttb_energy_cutoff + ttb_energy_photon, ttb use random_lcg, only: prn contains @@ -382,7 +382,7 @@ contains integer :: n integer :: n_e, n_k real(8) :: c_max - real(8) :: pi + real(8) :: w real(8) :: r real(8) :: e, e_l, e_r real(8) :: y, y_l, y_r @@ -413,17 +413,17 @@ contains ! Sample number of secondary bremsstrahlung photons n = floor(y + prn()) - ! Calculate the interpolation weight pi_j of the bremsstrahlung energy PDF + ! Calculate the interpolation weight w_j of the bremsstrahlung energy PDF ! interpolated in log energy, which can be interpreted as the probability ! of index j - pi = (e_r - e) / (e_r - e_l) + w = (e_r - e) / (e_r - e_l) ! Sample the energies of the emitted photons do i = 1, n ! Sample index of the tabulated PDF in the energy grid, j or j+1 i_e = j - if (prn() > pi) + if (prn() > w) then i_e = i_e + 1 end if @@ -451,7 +451,7 @@ contains x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l) ! Determine whether to deliver k - if (prn()*max(x_l, x_r) < x) exit + if (prn() * max(x_l, x_r) < x) exit end do ! Create secondary photon diff --git a/src/physics.F90 b/src/physics.F90 index 57bb48968f..8364962ee0 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -15,7 +15,8 @@ module physics use particle_restart_write, only: write_particle_restart use photon_header use photon_physics, only: rayleigh_scatter, compton_scatter, & - atomic_relaxation + atomic_relaxation, & + thick_target_bremsstrahlung use physics_common use random_lcg, only: prn, advance_prn_seed, prn_set_stream use reaction_header, only: Reaction @@ -304,9 +305,7 @@ contains ! TODO: create reaction types if (electron_treatment == ELECTRON_TTB) then - ! TODO: implement thick-target bremsstrahlung model - call fatal_error("Thick-target bremsstrahlung treatment of electrons & - &is not yet implemented.") + call thick_target_bremsstrahlung(p) end if p % E = ZERO @@ -333,7 +332,7 @@ contains if (electron_treatment == ELECTRON_TTB) then ! TODO: implement thick-target bremsstrahlung model - call fatal_error("Thick-target bremsstrahlung treatment of electrons & + call fatal_error("Thick-target bremsstrahlung treatment of positrons & &is not yet implemented.") end if