From 7fc8c32295def1c83ee4f7b25e2b679ce7561b7d Mon Sep 17 00:00:00 2001 From: amandalund Date: Fri, 8 Jun 2018 09:56:59 -0500 Subject: [PATCH] Added TTB positron treatment --- src/input_xml.F90 | 5 +-- src/material_header.F90 | 73 +++++++++++++++++++++++++++++++---------- src/photon_header.F90 | 11 ++++--- src/photon_physics.F90 | 16 +++++---- 4 files changed, 76 insertions(+), 29 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4d29cdc811..6370fc9fb6 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4370,9 +4370,10 @@ contains end if end do - ! Generate material bremsstrahlung data + ! Generate material bremsstrahlung data for electrons and positrons if (photon_transport .and. electron_treatment == ELECTRON_TTB) then - call bremsstrahlung_init(ttb(i), i) + call bremsstrahlung_init(ttb(i) % electron, i, ELECTRON) + call bremsstrahlung_init(ttb(i) % positron, i, POSITRON) end if end do diff --git a/src/material_header.F90 b/src/material_header.F90 index 6cdff7397c..4afff94c8a 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -699,9 +699,10 @@ contains end function openmc_material_set_densities - subroutine bremsstrahlung_init(this, i_material) - class(Bremsstrahlung), intent(inout) :: this + 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 @@ -711,6 +712,8 @@ contains 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) :: density real(8) :: density_gpcc @@ -720,33 +723,44 @@ contains real(8) :: mass_sum real(8), allocatable :: atom_fraction(:) real(8), allocatable :: mass_fraction(:) + 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) - this % i_material = i_material - ! Allocate and initialize arrays + ! Determine whether we are generating electron or positron data + if (particle == POSITRON) then + positron_ = .true. + else + positron_ = .false. + end if + + ! Get the size of the energy grids n_k = size(ttb_k_grid) n_e = size(ttb_e_grid) - allocate(this % pdf(n_e, n_e)) - allocate(this % cdf(n_e, n_e)) + + ! 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(atom_fraction(mat % n_nuclides)) allocate(mass_fraction(mat % n_nuclides)) + 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)) + allocate(dcs(n_k, n_e), source=ZERO) allocate(f(n_e)) allocate(z(n_e)) - this % pdf(:,:) = ZERO - this % cdf(:,:) = ZERO - stopping_power(:) = ZERO - dcs(:,:) = ZERO ! Calculate the "equivalent" atomic number Zeq, the atomic fraction and the ! mass fraction of each element, and the material density in atom/b-cm and @@ -768,6 +782,7 @@ contains Z_eq_sq = Z_eq_sq + atom_fraction(i) * nuclides(mat % nuclide(i)) % Z**2 end do + atom_sum = sum(atom_fraction) mass_sum = sum(mass_fraction) @@ -800,16 +815,39 @@ contains ! Get pointer to current element elm => elements(mat % element(i)) - ! TODO: for molecular DCS, atom_fraction should actually be the number of - ! atoms in the molecule. ! Accumulate material DCS dcs = dcs + atom_fraction(i) * elm % Z**2 / Z_eq_sq * elm % dcs - ! Accumulate material total stopping power - stopping_power = stopping_power + mass_fraction(i) * density_gpcc * & - (elm % stopping_power_collision + elm % stopping_power_radiative) + ! Accumulate material collision stopping power + stopping_power_collision = stopping_power_collision + & + mass_fraction(i) * density_gpcc * elm % stopping_power_collision + + ! Accumulate material radiative stopping power + stopping_power_radiative = stopping_power_radiative + & + mass_fraction(i) * density_gpcc * elm % stopping_power_radiative end do + ! 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)) + 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) @@ -923,7 +961,8 @@ contains this % yield = log(this % yield) end where - deallocate(atom_fraction, mass_fraction, stopping_power, dcs, f, z) + deallocate(atom_fraction, mass_fraction, stopping_power_collision, & + stopping_power_radiative, stopping_power, dcs, f, z) end subroutine bremsstrahlung_init diff --git a/src/photon_header.F90 b/src/photon_header.F90 index bb48c7f40a..f3f2a3c8da 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -77,12 +77,15 @@ module photon_header procedure :: calculate_xs => photon_calculate_xs end type PhotonInteraction - type Bremsstrahlung - integer :: i_material ! Index in materials array - + 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 @@ -90,7 +93,7 @@ module photon_header type(DictCharInt) :: element_dict - type(Bremsstrahlung), allocatable, target :: ttb(:) ! Bremsstrahlung cross sections + type(Bremsstrahlung), allocatable, target :: ttb(:) ! Bremsstrahlung data !=============================================================================== ! ELEMENTMICROXS contains cached microscopic photon cross sections for a diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 1a2e12379c..29724bb636 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -3,7 +3,7 @@ module photon_physics use algorithm, only: binary_search use constants use particle_header, only: Particle - use photon_header, only: PhotonInteraction, Bremsstrahlung, & + use photon_header, only: PhotonInteraction, BremsstrahlungData, & compton_profile_pz, ttb_e_grid, ttb use random_lcg, only: prn use settings @@ -530,9 +530,9 @@ contains real(8) :: w, w_l, w_r real(8) :: p_l, p_r real(8) :: c, c_l, c_max - type(Bremsstrahlung), pointer :: mat + type(BremsstrahlungData), pointer :: mat - real(8) :: photon_energies(100) + !real(8) :: photon_energies(100) !p % E = 100.0e6_8 @@ -544,8 +544,12 @@ contains ! return !end if - ! Get bremsstrahlung data for this material - mat => ttb(p % material) + ! 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) @@ -620,7 +624,7 @@ contains w = exp(w_l)*(a*(c - c_l)/(exp(w_l)*p_l) + ONE)**(ONE/a) end if - photon_energies(i) = w + !photon_energies(i) = w if (w > energy_cutoff(PHOTON)) then ! Create secondary photon call p % create_secondary(p % coord(1) % uvw, w, PHOTON, run_ce=.true.)