From 39cf5a8ce7891fcd59e05ac006aee1d571eac0bc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 8 Jan 2019 08:36:44 -0600 Subject: [PATCH] Move calculate_xs for photon data to C++ --- include/openmc/photon.h | 13 +- src/material_header.F90 | 11 +- src/photon.cpp | 70 ++++++++- src/photon_header.F90 | 317 +--------------------------------------- src/physics.cpp | 3 +- 5 files changed, 89 insertions(+), 325 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index fce6a3626..81675585f 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -39,9 +39,11 @@ public: class PhotonInteraction { public: // Constructors - PhotonInteraction(hid_t group); + PhotonInteraction(hid_t group, int i_element); // Methods + void calculate_xs(double E) const; + void compton_scatter(double alpha, bool doppler, double* alpha_out, double* mu, int* i_shell) const; @@ -53,8 +55,9 @@ public: void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const; // Data members - std::string name_; //! Name of element, e.g. "Zr" - int Z_; //! Atomic number + std::string name_; //!< Name of element, e.g. "Zr" + int Z_; //!< Atomic number + int i_element_; //!< Index in global elements vector // Microscopic cross sections xt::xtensor energy_; @@ -72,8 +75,8 @@ public: Tabulated1D coherent_anomalous_imag_; // Photoionization and atomic relaxation data - std::unordered_map shell_map_; // Given a shell designator, e.g. 3, this - // dictionary gives an index in shells_ + std::unordered_map shell_map_; //!< Given a shell designator, e.g. 3, this + //!< dictionary gives an index in shells_ std::vector shells_; // Compton profile data diff --git a/src/material_header.F90 b/src/material_header.F90 index 0a2548df3..934173241 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -487,6 +487,14 @@ contains integer :: i_element ! index into elements array real(8) :: atom_density ! atom density of a nuclide + interface + subroutine photon_calculate_xs(i_element, E) bind(C) + import C_INT, C_DOUBLE + integer(C_INT), value :: i_element + real(C_DOUBLE), value :: E + end subroutine + end interface + material_xs % coherent = ZERO material_xs % incoherent = ZERO material_xs % photoelectric = ZERO @@ -502,8 +510,7 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_photon_xs(i_element) % last_E) then - call elements(i_element) % calculate_xs(& - p % E, micro_photon_xs(i_element)) + call photon_calculate_xs(i_element, p % E) end if ! ======================================================================== diff --git a/src/photon.cpp b/src/photon.cpp index 56ff2042d..7432a2dcd 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -39,7 +39,8 @@ ElementMicroXS* micro_photon_xs; // PhotonInteraction implementation //============================================================================== -PhotonInteraction::PhotonInteraction(hid_t group) +PhotonInteraction::PhotonInteraction(hid_t group, int i_element) + : i_element_{i_element} { // Get name of nuclide from group, removing leading '/' name_ = object_name(group).substr(1); @@ -430,6 +431,66 @@ void PhotonInteraction::compton_doppler(double alpha, double mu, *i_shell = shell; } +void PhotonInteraction::calculate_xs(double E) const +{ + // Perform binary search on the element energy grid in order to determine + // which points to interpolate between + int n_grid = energy_.size(); + double log_E = std::log(E); + int i_grid; + if (log_E <= energy_[0]) { + i_grid = 0; + } else if (log_E > energy_(n_grid - 1)) { + i_grid = n_grid - 2; + } else { + // We use upper_bound_index here because sometimes photons are created with + // energies that exactly match a grid point + i_grid = upper_bound_index(energy_.cbegin(), energy_.cend(), log_E); + } + + // check for case where two energy points are the same + if (energy_(i_grid) == energy_(i_grid+1)) ++i_grid; + + // calculate interpolation factor + double f = (log_E - energy_(i_grid)) / (energy_(i_grid+1) - energy_(i_grid)); + + auto& xs {simulation::micro_photon_xs[i_element_]}; + xs.index_grid = i_grid; + xs.interp_factor = f; + + // Calculate microscopic coherent cross section + xs.coherent = std::exp(coherent_(i_grid) + + f*(coherent_(i_grid+1) - coherent_(i_grid))); + + // Calculate microscopic incoherent cross section + xs.incoherent = std::exp(incoherent_(i_grid) + + f*(incoherent_(i_grid+1) - incoherent_(i_grid))); + + // Calculate microscopic photoelectric cross section + xs.photoelectric = 0.0; + for (const auto& shell : shells_) { + // Check threshold of reaction + int i_start = shell.threshold; + if (i_grid < i_start) continue; + + // Evaluation subshell photoionization cross section + xs.photoelectric += + std::exp(shell.cross_section(i_grid-i_start) + + f*(shell.cross_section(i_grid+1-i_start) - + shell.cross_section(i_grid-i_start))); + } + + // Calculate microscopic pair production cross section + xs.pair_production = std::exp( + pair_production_total_(i_grid) + f*( + pair_production_total_(i_grid+1) - + pair_production_total_(i_grid))); + + // Calculate microscopic total cross section + xs.total = xs.coherent + xs.incoherent + xs.photoelectric + xs.pair_production; + xs.last_E = E; +} + double PhotonInteraction::rayleigh_scatter(double alpha) const { double mu; @@ -800,7 +861,12 @@ std::pair klein_nishina(double alpha) extern "C" void photon_from_hdf5_c(hid_t group) { - data::elements.emplace_back(group); + data::elements.emplace_back(group, data::elements.size()); +} + +extern "C" void photon_calculate_xs(int i_element, double E) +{ + data::elements[i_element - 1].calculate_xs(E); } } // namespace openmc diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 4be73831d..5a887d259 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -4,59 +4,20 @@ module photon_header use algorithm, only: binary_search use constants - use dict_header, only: DictIntInt, DictCharInt - use endf_header, only: Tabulated1D + use dict_header, only: DictCharInt use hdf5_interface - use nuclide_header, only: nuclides use settings real(8), allocatable :: compton_profile_pz(:) real(8), allocatable :: ttb_e_grid(:) ! energy T of incident electron real(8), allocatable :: ttb_k_grid(:) ! reduced energy W/T of emitted photon - type ElectronSubshell - integer :: index_subshell ! index in SUBSHELLS - integer :: threshold - real(8) :: n_electrons - real(8) :: binding_energy - real(8), allocatable :: cross_section(:) - - ! Transition data - integer :: n_transitions - integer, allocatable :: transition_subshells(:,:) - real(8), allocatable :: transition_energy(:) - real(8), allocatable :: transition_probability(:) - end type ElectronSubshell - type PhotonInteraction character(3) :: name ! atomic symbol, e.g. 'Zr' integer :: Z ! atomic number ! Microscopic cross sections real(8), allocatable :: energy(:) - real(8), allocatable :: coherent(:) - real(8), allocatable :: incoherent(:) - real(8), allocatable :: photoelectric_total(:) - real(8), allocatable :: pair_production_total(:) - real(8), allocatable :: pair_production_electron(:) - real(8), allocatable :: pair_production_nuclear(:) - - ! Form factors - type(Tabulated1D) :: incoherent_form_factor - type(Tabulated1D) :: coherent_int_form_factor - type(Tabulated1D) :: coherent_anomalous_real - type(Tabulated1D) :: coherent_anomalous_imag - - ! Photoionization and atomic relaxation data - type(DictIntInt) :: shell_dict ! Given a shell designator, e.g. 3, this - ! dictionary gives an index in shells(:) - type(ElectronSubshell), allocatable :: shells(:) - - ! Compton profile data - real(8), allocatable :: profile_pdf(:,:) - real(8), allocatable :: profile_cdf(:,:) - real(8), allocatable :: binding_energy(:) - real(8), allocatable :: electron_pdf(:) ! Stopping power data real(8) :: I ! mean excitation energy @@ -68,7 +29,6 @@ module photon_header contains procedure :: from_hdf5 => photon_from_hdf5 - procedure :: calculate_xs => photon_calculate_xs end type PhotonInteraction type BremsstrahlungData @@ -114,22 +74,16 @@ contains class(PhotonInteraction), intent(inout) :: this integer(HID_T), intent(in) :: group_id - integer :: i, j - integer(HID_T) :: rgroup, tgroup + integer :: i + integer(HID_T) :: rgroup integer(HID_T) :: dset_id integer(HSIZE_T) :: dims(1), dims2(2) integer :: n_energy - integer :: n_shell - integer :: n_profile - integer :: n_transition integer :: n_k integer :: n_e - character(3), allocatable :: designators(:) - real(8) :: c real(8) :: f real(8) :: y real(8), allocatable :: electron_energy(:) - real(8), allocatable :: matrix(:,:) real(8), allocatable :: dcs(:,:) interface @@ -159,171 +113,6 @@ contains call read_dataset(this % energy, dset_id) call close_dataset(dset_id) - ! Allocate arrays - allocate(this % coherent(n_energy)) - allocate(this % incoherent(n_energy)) - allocate(this % pair_production_total(n_energy)) - allocate(this % pair_production_nuclear(n_energy)) - allocate(this % pair_production_electron(n_energy)) - allocate(this % photoelectric_total(n_energy)) - - ! Read coherent scattering - rgroup = open_group(group_id, 'coherent') - call read_dataset(this % coherent, rgroup, 'xs') - - dset_id = open_dataset(rgroup, 'integrated_scattering_factor') - call this % coherent_int_form_factor % from_hdf5(dset_id) - call close_dataset(dset_id) - - if (object_exists(group_id, 'anomalous_real')) then - dset_id = open_dataset(rgroup, 'anomalous_real') - call this % coherent_anomalous_real % from_hdf5(dset_id) - call close_dataset(dset_id) - end if - - if (object_exists(group_id, 'anomalous_imag')) then - dset_id = open_dataset(rgroup, 'anomalous_imag') - call this % coherent_anomalous_imag % from_hdf5(dset_id) - call close_dataset(dset_id) - call close_group(rgroup) - end if - - ! Read incoherent scattering - rgroup = open_group(group_id, 'incoherent') - call read_dataset(this % incoherent, rgroup, 'xs') - dset_id = open_dataset(rgroup, 'scattering_factor') - call this % incoherent_form_factor % from_hdf5(dset_id) - call close_dataset(dset_id) - call close_group(rgroup) - - ! Read pair production - rgroup = open_group(group_id, 'pair_production_electron') - call read_dataset(this % pair_production_electron, rgroup, 'xs') - call close_group(rgroup) - - ! Read pair production - if (object_exists(group_id, 'pair_production_nuclear')) then - rgroup = open_group(group_id, 'pair_production_nuclear') - call read_dataset(this % pair_production_nuclear, rgroup, 'xs') - call close_group(rgroup) - else - this % pair_production_nuclear(:) = ZERO - end if - - ! Read photoelectric - rgroup = open_group(group_id, 'photoelectric') - call read_dataset(this % photoelectric_total, rgroup, 'xs') - call close_group(rgroup) - - ! Read subshell photoionization cross section and atomic relaxation data - rgroup = open_group(group_id, 'subshells') - call read_attribute(designators, rgroup, 'designators') - n_shell = size(designators) - allocate(this % shells(n_shell)) - do i = 1, n_shell - ! Create mapping from designator to index - do j = 1, size(SUBSHELLS) - if (designators(i) == SUBSHELLS(j)) then - call this % shell_dict % set(j, i) - this % shells(i) % index_subshell = j - exit - end if - end do - - ! Read binding energy and number of electrons - tgroup = open_group(rgroup, trim(designators(i))) - call read_attribute(this % shells(i) % binding_energy, tgroup, & - 'binding_energy') - call read_attribute(this % shells(i) % n_electrons, tgroup, & - 'num_electrons') - - ! Read subshell cross section - dset_id = open_dataset(tgroup, 'xs') - call read_attribute(j, dset_id, 'threshold_idx') - this % shells(i) % threshold = j - allocate(this % shells(i) % cross_section(n_energy - j)) - call read_dataset(this % shells(i) % cross_section, dset_id) - call close_dataset(dset_id) - where (this % shells(i) % cross_section > ZERO) - this % shells(i) % cross_section = log(this % shells(i) % cross_section) - elsewhere - this % shells(i) % cross_section = -500.0_8 - end where - - if (object_exists(tgroup, 'transitions')) then - dset_id = open_dataset(tgroup, 'transitions') - call get_shape(dset_id, dims2) - n_transition = int(dims2(2), 4) - this % shells(i) % n_transitions = n_transition - if (n_transition > 0) then - allocate(this % shells(i) % transition_subshells(2, n_transition)) - allocate(this % shells(i) % transition_energy(n_transition)) - allocate(this % shells(i) % transition_probability(n_transition)) - - allocate(matrix(dims2(1), dims2(2))) - call read_dataset(matrix, dset_id) - - this % shells(i) % transition_subshells(:,:) = int(matrix(1:2, :), 4) - this % shells(i) % transition_energy(:) = matrix(3, :) - this % shells(i) % transition_probability(:) = matrix(4, :) & - / sum(matrix(4, :)) - deallocate(matrix) - end if - call close_dataset(dset_id) - else - this % shells(i) % n_transitions = 0 - end if - call close_group(tgroup) - end do - call close_group(rgroup) - deallocate(designators) - - ! Determine number of electron shells - rgroup = open_group(group_id, 'compton_profiles') - - ! Determine number of shells - dset_id = open_dataset(rgroup, 'num_electrons') - call get_shape(dset_id, dims) - n_shell = int(dims(1), 4) - - ! Read electron shell PDF and binding energies - allocate(this % electron_pdf(n_shell), this % binding_energy(n_shell)) - call read_dataset(this % electron_pdf, dset_id) - call close_dataset(dset_id) - call read_dataset(this % binding_energy, rgroup, 'binding_energy') - this % electron_pdf(:) = this % electron_pdf / sum(this % electron_pdf) - - ! Read Compton profiles - dset_id = open_dataset(rgroup, 'J') - call get_shape(dset_id, dims2) - n_profile = int(dims2(1), 4) - allocate(this % profile_pdf(n_profile, n_shell)) - call read_dataset(this % profile_pdf, dset_id) - call close_dataset(dset_id) - - ! Get Compton profile momentum grid - if (.not. allocated(compton_profile_pz)) then - allocate(compton_profile_pz(n_profile)) - call read_dataset(compton_profile_pz, rgroup, 'pz') - end if - call close_group(rgroup) - - ! Create Compton profile CDF - allocate(this % profile_cdf(n_profile, n_shell)) - do i = 1, n_shell - c = ZERO - this % profile_cdf(1,i) = ZERO - do j = 1, n_profile - 1 - c = c + HALF*(compton_profile_pz(j+1) - compton_profile_pz(j)) * & - (this%profile_pdf(j,i) + this%profile_pdf(j+1,i)) - this % profile_cdf(j+1,i) = c - end do - end do - - ! Calculate total pair production - this % pair_production_total(:) = this % pair_production_nuclear + & - this % pair_production_electron - if (electron_treatment == ELECTRON_TTB) then ! Read bremsstrahlung scaled DCS rgroup = open_group(group_id, 'bremsstrahlung') @@ -401,108 +190,8 @@ contains ! interpolated this % energy = log(this % energy) - where (this % coherent > ZERO) - this % coherent = log(this % coherent) - elsewhere - this % coherent = -500.0_8 - end where - - where (this % incoherent > ZERO) - this % incoherent = log(this % incoherent) - elsewhere - this % incoherent = -500.0_8 - end where - - where (this % photoelectric_total > ZERO) - this % photoelectric_total = log(this % photoelectric_total) - elsewhere - this % photoelectric_total = -500.0_8 - end where - - where (this % pair_production_total > ZERO) - this % pair_production_total = log(this % pair_production_total) - elsewhere - this % pair_production_total = -500.0_8 - end where - end subroutine photon_from_hdf5 -!=============================================================================== -! CALCULATE_ELEMENT_XS determines microscopic photon cross sections for an -! element of a given index in the elements array at the energy of the given -! particle -!=============================================================================== - - subroutine photon_calculate_xs(this, E, xs) - class(PhotonInteraction), intent(in) :: this ! index into elements array - real(8), intent(in) :: E ! energy - type(ElementMicroXS), intent(inout) :: xs - - integer :: i_grid ! index on element energy grid - integer :: i_shell ! index in subshells - integer :: i_start ! threshold index - integer :: n_grid ! number of grid points - real(8) :: f ! interp factor on element energy grid - real(8) :: log_E ! logarithm of the energy - - ! Perform binary search on the element energy grid in order to determine - ! which points to interpolate between - n_grid = size(this % energy) - log_E = log(E) - if (log_E <= this % energy(1)) then - i_grid = 1 - elseif (log_E > this % energy(n_grid)) then - i_grid = n_grid - 1 - else - i_grid = binary_search(this % energy, n_grid, log_E) - end if - - ! check for case where two energy points are the same - if (this % energy(i_grid) == this % energy(i_grid+1)) i_grid = i_grid + 1 - - ! calculate interpolation factor - f = (log_E - this % energy(i_grid)) / & - (this % energy(i_grid+1) - this % energy(i_grid)) - - xs % index_grid = i_grid - xs % interp_factor = f - - ! Calculate microscopic coherent cross section - xs % coherent = exp(this % coherent(i_grid) + f * & - (this % coherent(i_grid+1) - this % coherent(i_grid))) - - ! Calculate microscopic incoherent cross section - xs % incoherent = exp(this % incoherent(i_grid) + & - f*(this % incoherent(i_grid+1) - this % incoherent(i_grid))) - - ! Calculate microscopic photoelectric cross section - xs % photoelectric = ZERO - do i_shell = 1, size(this % shells) - ! Check threshold of reaction - i_start = this % shells(i_shell) % threshold - if (i_grid <= i_start) cycle - - ! Evaluation subshell photoionization cross section - xs % photoelectric = xs % photoelectric + & - exp(this % shells(i_shell) % cross_section(i_grid-i_start) + & - f*(this % shells(i_shell) % cross_section(i_grid+1-i_start) - & - this % shells(i_shell) % cross_section(i_grid-i_start))) - end do - - ! Calculate microscopic pair production cross section - xs % pair_production = exp(& - this % pair_production_total(i_grid) + f*(& - this % pair_production_total(i_grid+1) - & - this % pair_production_total(i_grid))) - - ! Calculate microscopic total cross section - xs % total = xs % coherent + xs % incoherent + xs % photoelectric + & - xs % pair_production - - xs % last_E = E - - end subroutine photon_calculate_xs - !=============================================================================== ! FREE_MEMORY_PHOTON deallocates/resets global variables in this module !=============================================================================== diff --git a/src/physics.cpp b/src/physics.cpp index a4bd67066..dbc6252eb 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -297,8 +297,7 @@ void sample_photon_reaction(Particle* p) if (prob_after > cutoff) { for (const auto& shell : element.shells_) { // Get grid index and interpolation factor - // TODO: off-by-one - int i_grid = micro.index_grid - 1; + int i_grid = micro.index_grid; double f = micro.interp_factor; // Check threshold of reaction