Move calculate_xs for photon data to C++

This commit is contained in:
Paul Romano 2019-01-08 08:36:44 -06:00
parent 2e6ac03433
commit 39cf5a8ce7
5 changed files with 89 additions and 325 deletions

View file

@ -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<double, 1> energy_;
@ -72,8 +75,8 @@ public:
Tabulated1D coherent_anomalous_imag_;
// Photoionization and atomic relaxation data
std::unordered_map<int, int> shell_map_; // Given a shell designator, e.g. 3, this
// dictionary gives an index in shells_
std::unordered_map<int, int> shell_map_; //!< Given a shell designator, e.g. 3, this
//!< dictionary gives an index in shells_
std::vector<ElectronSubshell> shells_;
// Compton profile data

View file

@ -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
! ========================================================================

View file

@ -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<double, double> 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

View file

@ -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
!===============================================================================

View file

@ -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