From 94dd50422ff3a4cc719ed776b312f2bf66beb694 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 28 Jun 2017 14:32:03 -0500 Subject: [PATCH] Initial implementation of photon physics --- CMakeLists.txt | 1 + openmc/source.py | 3 +- src/bank_header.F90 | 1 + src/cross_section.F90 | 155 ++++++++++++++++- src/global.F90 | 7 +- src/initialize.F90 | 18 +- src/nuclide_header.F90 | 8 +- src/particle_header.F90 | 13 +- src/photon_header.F90 | 16 ++ src/photon_physics.F90 | 365 ++++++++++++++++++++++++++++++++++++++++ src/physics.F90 | 163 +++++++++++++++++- src/simulation.F90 | 2 + src/source.F90 | 3 + 13 files changed, 732 insertions(+), 23 deletions(-) create mode 100644 src/photon_physics.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 367e1e5232..de01e2e543 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -342,6 +342,7 @@ set(LIBOPENMC_FORTRAN_SRC src/particle_restart.F90 src/particle_restart_write.F90 src/photon_header.F90 + src/photon_physics.F90 src/physics_common.F90 src/physics.F90 src/physics_mg.F90 diff --git a/openmc/source.py b/openmc/source.py index ca229638d3..4695f8ae03 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -128,7 +128,8 @@ class Source(object): """ element = ET.Element("source") element.set("strength", str(self.strength)) - element.set("particle", self.particle) + if self.particle != 'neutron': + element.set("particle", self.particle) if self.file is not None: element.set("file", self.file) if self.space is not None: diff --git a/src/bank_header.F90 b/src/bank_header.F90 index 97fb1f11fa..5444d6b146 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -16,6 +16,7 @@ module bank_header real(C_DOUBLE) :: uvw(3) ! diretional cosines real(C_DOUBLE) :: E ! energy / energy group if in MG mode. integer(C_INT) :: delayed_group ! delayed group + integer(C_INT) :: particle ! particle type (neutron, photon, etc.) end type Bank end module bank_header diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 28014f773a..918b06ad83 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -20,11 +20,26 @@ module cross_section contains !=============================================================================== -! CALCULATE_XS determines the macroscopic cross sections for the material the -! particle is currently traveling through. +! CALCULATE_XS determines the macroscopic neutron and/or photon cross sections +! for the material the particle is currently traveling through. !=============================================================================== subroutine calculate_xs(p) + type(Particle), intent(inout) :: p + + if (p % type == NEUTRON) then + call calculate_neutron_xs(p) + elseif (p % type == PHOTON) then + call calculate_photon_xs(p) + end if + end subroutine calculate_xs + +!=============================================================================== +! CALCULATE_NEUTRON_XS determines the macroscopic cross sections for the +! material the particle is currently traveling through. +!=============================================================================== + + subroutine calculate_neutron_xs(p) type(Particle), intent(inout) :: p @@ -125,7 +140,7 @@ contains end do end associate - end subroutine calculate_xs + end subroutine calculate_neutron_xs !=============================================================================== ! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a @@ -557,6 +572,140 @@ contains end subroutine calculate_urr_xs +!=============================================================================== +! CALCULATE_PHOTON_XS determines the macroscopic photon cross sections for the +! material the particle is currently traveling through. +!=============================================================================== + + subroutine calculate_photon_xs(p) + type(Particle), intent(inout) :: p + + integer :: i ! loop index over nuclides + integer :: i_element ! index into elements array + real(8) :: atom_density ! atom density of a nuclide + + ! Set all material macroscopic cross sections to zero + material_xs % total = ZERO + material_xs % coherent = ZERO + material_xs % incoherent = ZERO + material_xs % photoelectric = ZERO + material_xs % pair_production = ZERO + + ! Exit subroutine if material is void + if (p % material == MATERIAL_VOID) return + + associate (mat => materials(p % material)) + ! Add contribution from each nuclide in material + do i = 1, mat % n_nuclides + ! ======================================================================== + ! CALCULATE MICROSCOPIC CROSS SECTION + + ! Determine microscopic cross sections for this nuclide + i_element = mat % element(i) + + ! Calculate microscopic cross section for this nuclide + if (p % E /= micro_photon_xs(i_element) % last_E) then + call calculate_element_xs(i_element, p % E) + end if + + ! ======================================================================== + ! ADD TO MACROSCOPIC CROSS SECTION + + ! Copy atom density of nuclide in material + atom_density = mat % atom_density(i) + + ! Add contributions to material macroscopic total cross section + material_xs % total = material_xs % total + & + atom_density * micro_photon_xs(i_element) % total + + ! Add contributions to material macroscopic coherent cross section + material_xs % coherent = material_xs % coherent + & + atom_density * micro_photon_xs(i_element) % coherent + + ! Add contributions to material macroscopic incoherent cross section + material_xs % incoherent = material_xs % incoherent + & + atom_density * micro_photon_xs(i_element) % incoherent + + ! Add contributions to material macroscopic photoelectric cross section + material_xs % photoelectric = material_xs % photoelectric + & + atom_density * micro_photon_xs(i_element) % photoelectric + + ! Add contributions to material macroscopic pair production cross section + material_xs % pair_production = material_xs % pair_production + & + atom_density * micro_photon_xs(i_element) % pair_production + end do + end associate + + end subroutine calculate_photon_xs + +!=============================================================================== +! 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 calculate_element_xs(i_element, E) + integer, intent(in) :: i_element ! index into nuclides array + real(8), intent(in) :: E ! energy + + integer :: i_grid ! index on nuclide energy grid + integer :: n_grid ! number of grid points + real(8) :: f ! interp factor on nuclide energy grid + real(8) :: log_E ! logarithm of the energy + + associate (elm => elements(i_element)) + ! Perform binary search on the element energy grid in order to determine + ! which points to interpolate between + n_grid = size(elm % energy) + log_E = log(E) + if (log_E <= elm % energy(1)) then + i_grid = 1 + elseif (log_E > elm % energy(n_grid)) then + i_grid = n_grid - 1 + else + i_grid = binary_search(elm % energy, n_grid, log_E) + end if + + ! check for case where two energy points are the same + if (elm % energy(i_grid) == elm % energy(i_grid+1)) i_grid = i_grid + 1 + + ! calculate interpolation factor + f = (log_E - elm%energy(i_grid))/(elm%energy(i_grid+1) - elm%energy(i_grid)) + + micro_photon_xs(i_element) % index_grid = i_grid + micro_photon_xs(i_element) % interp_factor = f + + ! Calculate microscopic coherent cross section + micro_photon_xs(i_element) % coherent = exp(elm % coherent(i_grid) + f * & + (elm % coherent(i_grid+1) - elm % coherent(i_grid))) + + ! Calculate microscopic incoherent cross section + micro_photon_xs(i_element) % incoherent = exp(elm % incoherent(i_grid) + & + f*(elm % incoherent(i_grid+1) - elm % incoherent(i_grid))) + + ! Calculate microscopic photoelectric cross section + micro_photon_xs(i_element) % photoelectric = exp(elm % photoelectric_total(& + i_grid) + f*(elm % photoelectric_total(i_grid+1) - & + elm % photoelectric_total(i_grid))) + + ! Calculate microscopic pair production cross section + micro_photon_xs(i_element) % pair_production = exp(& + elm % pair_production_total(i_grid) + f*(& + elm % pair_production_total(i_grid+1) - & + elm % pair_production_total(i_grid))) + + ! Calculate microscopic total cross section + micro_photon_xs(i_element) % total = & + micro_photon_xs(i_element) % coherent + & + micro_photon_xs(i_element) % incoherent + & + micro_photon_xs(i_element) % photoelectric + & + micro_photon_xs(i_element) % pair_production + + micro_photon_xs(i_element) % last_E = E + end associate + + end subroutine calculate_element_xs + !=============================================================================== ! FIND_ENERGY_INDEX determines the index on the union energy grid at a certain ! energy diff --git a/src/global.F90 b/src/global.F90 index e42e04a8d7..6a37f2d748 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -15,7 +15,7 @@ module global use mesh_header, only: RegularMesh use mgxs_header, only: Mgxs, MgxsContainer use nuclide_header - use photon_header, only: PhotonInteraction + use photon_header, only: PhotonInteraction, ElementMicroXS use plot_header, only: ObjectPlot use sab_header, only: SAlphaBeta use set_header, only: SetInt @@ -79,6 +79,7 @@ module global ! Cross section caches type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide + type(ElementMicroXS), allocatable :: micro_photon_xs(:) ! Cache for each element type(MaterialMacroXS) :: material_xs ! Cache for current material ! Dictionaries to look up cross sections and listings @@ -452,8 +453,8 @@ module global real(8) :: res_scat_energy_max = 1000.0_8 character(10), allocatable :: res_scat_nuclides(:) -!$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, & -!$omp& trace, thread_id, current_work, filter_matches) +!$omp threadprivate(micro_xs, micro_photon_xs, material_xs, fission_bank, & +!$omp& n_bank, trace, thread_id, current_work, filter_matches) contains diff --git a/src/initialize.F90 b/src/initialize.F90 index 98f02863e2..8204307c99 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -163,13 +163,13 @@ contains integer, intent(in) :: intracomm ! MPI intracommunicator #endif - integer :: bank_blocks(5) ! Count for each datatype + integer :: bank_blocks(6) ! Count for each datatype #ifdef MPIF08 - type(MPI_Datatype) :: bank_types(5) + type(MPI_Datatype) :: bank_types(6) #else - integer :: bank_types(5) ! Datatypes + integer :: bank_types(6) ! Datatypes #endif - integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements + integer(MPI_ADDRESS_KIND) :: bank_disp(6) ! Displacements logical :: init_called type(Bank) :: b @@ -201,14 +201,16 @@ contains call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err) call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) call MPI_GET_ADDRESS(b % delayed_group, bank_disp(5), mpi_err) + call MPI_GET_ADDRESS(b % particle, bank_disp(6), mpi_err) ! Adjust displacements bank_disp = bank_disp - bank_disp(1) ! Define MPI_BANK for fission sites - bank_blocks = (/ 1, 3, 3, 1, 1 /) - bank_types = (/ MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_INTEGER /) - call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, & + bank_blocks = (/ 1, 3, 3, 1, 1, 1 /) + bank_types = (/ MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, & + MPI_INT, MPI_INT /) + call MPI_TYPE_CREATE_STRUCT(6, bank_blocks, bank_disp, & bank_types, MPI_BANK, mpi_err) call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) @@ -245,6 +247,8 @@ contains c_loc(tmpb(1)%E)), H5T_NATIVE_DOUBLE, hdf5_err) call h5tinsert_f(hdf5_bank_t, "delayed_group", h5offsetof(c_loc(tmpb(1)), & c_loc(tmpb(1)%delayed_group)), H5T_NATIVE_INTEGER, hdf5_err) + call h5tinsert_f(hdf5_bank_t, "particle", h5offsetof(c_loc(tmpb(1)), & + c_loc(tmpb(1)%particle)), H5T_NATIVE_INTEGER, hdf5_err) ! Determine type for integer(8) hdf5_integer8_t = h5kind_to_type(8, H5_INTEGER_KIND) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index c03ebee1a7..9cc2bc7e20 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -143,6 +143,12 @@ module nuclide_header real(8) :: absorption ! macroscopic absorption xs real(8) :: fission ! macroscopic fission xs real(8) :: nu_fission ! macroscopic production xs + + ! Photon cross sections + real(8) :: coherent ! macroscopic coherent xs + real(8) :: incoherent ! macroscopic incoherent xs + real(8) :: photoelectric ! macroscopic photoelectric xs + real(8) :: pair_production ! macroscopic pair production xs end type MaterialMacroXS !=============================================================================== @@ -155,7 +161,7 @@ module nuclide_header character(MAX_FILE_LEN) :: path end type Library - contains +contains !=============================================================================== ! NUCLIDE_CLEAR resets and deallocates data in Nuclide diff --git a/src/particle_header.F90 b/src/particle_header.F90 index f28e549761..f0942fc41c 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -199,6 +199,7 @@ contains call this % initialize() ! copy attributes from source bank site + this % type = src % particle this % wgt = src % wgt this % last_wgt = src % wgt this % coord(1) % xyz = src % xyz @@ -222,9 +223,10 @@ contains ! the secondary bank and increments the number of sites in the secondary bank. !=============================================================================== - subroutine create_secondary(this, uvw, type, run_CE) + subroutine create_secondary(this, uvw, E, type, run_CE) class(Particle), intent(inout) :: this real(8), intent(in) :: uvw(3) + real(8), intent(in) :: E integer, intent(in) :: type logical, intent(in) :: run_CE @@ -237,14 +239,15 @@ contains end if n = this % n_secondary + 1 - this % secondary_bank(n) % wgt = this % wgt + this % secondary_bank(n) % particle = type + this % secondary_bank(n) % wgt = this % wgt this % secondary_bank(n) % xyz(:) = this % coord(1) % xyz this % secondary_bank(n) % uvw(:) = uvw - this % n_secondary = n - this % secondary_bank(this % n_secondary) % E = this % E + this % secondary_bank(n) % E = E if (.not. run_CE) then - this % secondary_bank(this % n_secondary) % E = real(this % g, 8) + this % secondary_bank(n) % E = real(this % g, 8) end if + this % n_secondary = n end subroutine create_secondary diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 30f1b7cfed..dfedd5b0b9 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -56,6 +56,22 @@ module photon_header procedure :: from_hdf5 => photon_from_hdf5 end type PhotonInteraction +!=============================================================================== +! ELEMENTMICROXS contains cached microscopic photon cross sections for a +! 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 + end type ElementMicroXS + contains subroutine photon_from_hdf5(this, group_id) diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 new file mode 100644 index 0000000000..7b7244fe70 --- /dev/null +++ b/src/photon_physics.F90 @@ -0,0 +1,365 @@ +module photon_physics + + use algorithm, only: binary_search + use constants + use particle_header, only: Particle + use photon_header, only: PhotonInteraction, compton_profile_pz + use random_lcg, only: prn + +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, use_doppler) + type(PhotonInteraction), intent(in) :: el + real(8), intent(in) :: alpha + real(8), intent(out) :: alpha_out + real(8), intent(out) :: mu + 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/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/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) + alpha_out = e_out/MASS_ELECTRON + end if + exit + end if + end do + + end subroutine compton_scatter + +!=============================================================================== +! COMPTON_DOPPLER +!=============================================================================== + + subroutine compton_doppler(el, alpha, mu, e_out) + type(PhotonInteraction), intent(in) :: el + real(8), intent(in) :: alpha + real(8), intent(in) :: mu + real(8), intent(out) :: e_out + + integer :: i, i_shell + 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) - 1 + c = c + el % electron_pdf(i_shell + 1) + 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 + if (e < e_b) then + e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON + 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 + 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 + 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 + exit + end if + quad = sqrt(quad) + e_out1 = -(b + quad)/(TWO*a) + e_out2 = -(b - quad)/(TWO*a) + + if (e_out1 > ZERO) then + if (e_out2 > ZERO) then + 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) e_out = e_out2 + 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/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) + + ! Check for no transitions + if (elm % shells(i_shell) % n_transitions == 0) return + + ! Sample transition + rn = prn() + c = ZERO + do i_transition = 1, elm % shells(i_shell) % n_transitions - 1 + c = c + elm % shells(i_shell) % & + transition_probability(i_transition + 1) + 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) + + if (secondary == 0) then + ! Non-radiative trnasition -- Auger/Coster-Kronig effect + + ! TODO: Create electron + ! E_electron = transition_energy(i_transition) + + ! Fill secondary (higher) hole first + if (elm % shell_dict % has_key(secondary)) then + i_hole = elm % shell_dict % get_key(secondary) + call atomic_relaxation(p, elm, i_hole) + end if + else + ! Radiative transition -- get X-ray energy + E = elm % shells(i_shell) % transition_energy(i_transition) + + if (E > ZERO) then + ! Sample angle isotropically for X-ray + mu = TWO*prn() - ONE + phi = TWO*PI*prn() + uvw(1) = mu + uvw(2) = sqrt(ONE - mu*mu)*cos(phi) + uvw(3) = sqrt(ONE - mu*mu)*sin(phi) + + ! Create X-ray + call p % create_secondary(uvw, E, PHOTON, run_ce=.true.) + end if + end if + + ! Fill primary hole + if (elm % shell_dict % has_key(primary)) then + i_hole = elm % shell_dict % get_key(primary) + call atomic_relaxation(p, elm, i_hole) + end if + + end subroutine atomic_relaxation + +end module photon_physics diff --git a/src/physics.F90 b/src/physics.F90 index 9482772ef1..d6a4c0ac52 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -14,6 +14,8 @@ module physics use output, only: write_message use particle_header, only: Particle use particle_restart_write, only: write_particle_restart + use photon_physics, only: rayleigh_scatter, compton_scatter, & + atomic_relaxation use physics_common use random_lcg, only: prn, advance_prn_seed, prn_set_stream use reaction_header, only: Reaction @@ -41,8 +43,12 @@ contains ! Add to collision counter for particle p % n_collision = p % n_collision + 1 - ! Sample nuclide/reaction for the material the particle is in - call sample_reaction(p) + ! Sample reaction for the material the particle is in + if (p % type == NEUTRON) then + call sample_reaction(p) + elseif (p % type == PHOTON) then + call sample_photon_reaction(p) + end if ! Display information about collision if (verbosity >= 10 .or. trace) then @@ -137,6 +143,110 @@ contains end subroutine sample_reaction +!=============================================================================== +! SAMPLE_PHOTON_REACTION samples an element based on the macroscopic cross +! sections for each nuclide within a material and then samples a reaction for +! that element and calls the appropriate routine to process the physics. +!=============================================================================== + + subroutine sample_photon_reaction(p) + 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) :: 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) :: phi ! azimuthal angle + + ! 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 + + ! 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) + 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, .true.) + p % E = alpha_out*MASS_ELECTRON + p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) + 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 + ! TODO: Create electron + ! E_electron = p % E - elm % shells(i_shell) % binding_energy + + call atomic_relaxation(p, elm, i_shell) + p % alive = .false. + return + end if + end do + end if + prob = prob_after + end associate + + ! Pair production + prob = prob + micro_photon_xs(i_element) % pair_production + if (prob > cutoff) then + ! Sample angle isotropically + mu = TWO*prn() - ONE + phi = TWO*PI*prn() + p % coord(1) % uvw(1) = mu + p % coord(1) % uvw(2) = sqrt(ONE - mu*mu)*cos(phi) + p % coord(1) % uvw(3) = sqrt(ONE - mu*mu)*sin(phi) + + ! Set energy + p % E = MASS_ELECTRON + + ! Create photon in opposite direction + call p % create_secondary(-p % coord(1) % uvw, MASS_ELECTRON, & + PHOTON, .true.) + end if + + end subroutine sample_photon_reaction + !=============================================================================== ! SAMPLE_NUCLIDE !=============================================================================== @@ -199,6 +309,49 @@ contains end subroutine sample_nuclide +!=============================================================================== +! SAMPLE_ELEMENT +!=============================================================================== + + function sample_element(p) result(i_element) + type(Particle), intent(in) :: p + integer :: 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 write_particle_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 + !=============================================================================== ! SAMPLE_FISSION !=============================================================================== @@ -1140,6 +1293,9 @@ contains ! Bank source neutrons by copying particle data bank_array(i) % xyz = p % coord(1) % xyz + ! Set particle as neutron + bank_array(i) % particle = NEUTRON + ! Set weight of fission bank site bank_array(i) % wgt = ONE/weight @@ -1329,7 +1485,8 @@ contains if (mod(yield, ONE) == ZERO) then ! If yield is integral, create exactly that many secondary particles do i = 1, nint(yield) - 1 - call p % create_secondary(p % coord(1) % uvw, NEUTRON, run_CE=.true.) + call p % create_secondary(p % coord(1) % uvw, p % E, & + NEUTRON, run_CE=.true.) end do else ! Otherwise, change weight of particle based on yield diff --git a/src/simulation.F90 b/src/simulation.F90 index 87a17452d9..2c26bdc31d 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -378,6 +378,7 @@ contains !$omp parallel allocate(micro_xs(n_nuclides_total)) + allocate(micro_photon_xs(n_elements)) !$omp end parallel if (.not. restart_run) call initialize_source() @@ -403,6 +404,7 @@ contains !$omp parallel deallocate(micro_xs) + deallocate(micro_photon_xs) !$omp end parallel ! Increment total number of generations diff --git a/src/source.F90 b/src/source.F90 index 517715b563..cf7cd5b329 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -134,6 +134,9 @@ contains ! Set particle defaults call p % initialize() + ! Set particle type + site % particle = external_source(i) % particle + ! Sample spatial distribution site % xyz(:) = external_source(i) % space % sample()