crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-25 04:25:29 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

TTB implementation fixes

Browse source
This commit is contained in:
amandalund 2018-02-19 22:27:47 -06:00
parent 62985f6198
commit 150bedcb50
8 changed files with 136 additions and 55 deletions
Download patch file Download diff file Expand all files Collapse all files

1
openmc/data/photon.py
View file

@ -7,6 +7,7 @@ import os
import h5py
import numpy as np
import pandas as pd
from scipy.interpolate import CubicSpline
from openmc.mixin import EqualityMixin
import openmc.checkvalue as cv

9
openmc/filter.py
View file

@ -565,7 +565,7 @@ class ParticleFilter(WithIDFilter):
Attributes
----------
bins : Iterable of Integral
openmc.Materi IDs.
openmc.Material IDs.
id : int
Unique identifier for the filter
num_bins : Integral
@ -582,7 +582,12 @@ class ParticleFilter(WithIDFilter):
@bins.setter
def bins(self, bins):
bins = np.atleast_1d(bins)
cv.check_iterable_type('filter bins', bins, str)
cv.check_iterable_type('filter bins', bins, (Integral, str))
for edge in bins:
if isinstance(edge, Integral):
cv.check_value('filter bin', edge, _PARTICLE_IDS.values())
else:
cv.check_value('filter bin', edge, _PARTICLE_IDS.keys())
bins = np.atleast_1d([b if isinstance(b, Integral) else _PARTICLE_IDS[b]
for b in bins])
self._bins = bins

19
openmc/settings.py
View file

@ -45,6 +45,9 @@ class Settings(object):
weight assigned to particles that are not killed after Russian
roulette. Value of energy should be a float indicating energy in eV
below which particle type will be killed.
electron_treatment : {'led', 'ttb'}
Whether to deposit all energy from electrons locally ('led') or create
secondary bremsstrahlung photons ('ttb').
energy_mode : {'continuous-energy', 'multi-group'}
Set whether the calculation should be continuous-energy or multi-group.
entropy_mesh : openmc.Mesh
@ -185,6 +188,7 @@ class Settings(object):
self._confidence_intervals = None
self._cross_sections = None
self._electron_treatment = None
self._multipole_library = None
self._photon_transport = None
self._ptables = None
@ -282,6 +286,10 @@ class Settings(object):
def cross_sections(self):
return self._cross_sections
@property
def electron_treatment(self):
return self._electron_treatment
@property
def multipole_library(self):
return self._multipole_library
@ -547,6 +555,11 @@ class Settings(object):
cv.check_type('cross sections', cross_sections, string_types)
self._cross_sections = cross_sections
@electron_treatment.setter
def electron_treatment(self, electron_treatment):
cv.check_value('electron treatment', electron_treatment, ['led', 'ttb'])
self._electron_treatment = electron_treatment
@multipole_library.setter
def multipole_library(self, multipole_library):
warnings.warn('Settings.multipole_library has been deprecated and will '
@ -931,6 +944,11 @@ class Settings(object):
element = ET.SubElement(root, "cross_sections")
element.text = str(self._cross_sections)
def _create_electron_treatment_subelement(self, root):
if self._electron_treatment is not None:
element = ET.SubElement(root, "electron_treatment")
element.text = str(self._electron_treatment)
def _create_multipole_library_subelement(self, root):
if self._multipole_library is not None:
element = ET.SubElement(root, "multipole_library")
@ -1121,6 +1139,7 @@ class Settings(object):
self._create_confidence_intervals(root_element)
self._create_cross_sections_subelement(root_element)
self._create_multipole_library_subelement(root_element)
self._create_electron_treatment_subelement(root_element)
self._create_energy_mode_subelement(root_element)
self._create_max_order_subelement(root_element)
self._create_photon_transport_subelement(root_element)

15
src/input_xml.F90
View file

@ -4620,6 +4620,21 @@ contains
end if
end do
if (photon_transport .and. electron_treatment == ELECTRON_TTB) then
! Deallocate element bremsstrahlung DCS and stopping power data since
! only the material bremsstrahlung data is needed
do i = 1, size(elements)
if (allocated(elements(i) % stopping_power_collision)) &
deallocate(elements(i) % stopping_power_collision)
if (allocated(elements(i) % stopping_power_radiative)) &
deallocate(elements(i) % stopping_power_radiative)
if (allocated(elements(i) % dcs)) deallocate(elements(i) % dcs)
end do
! Take logarithm of electron energies since they are log-log interpolated
ttb_energy_electron = log(ttb_energy_electron)
end if
! Set up logarithmic grid for nuclides
do i = 1, size(nuclides)
call nuclides(i) % init_grid(energy_min_neutron, &

94
src/photon_header.F90
View file

@ -3,7 +3,7 @@ module photon_header
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use algorithm, only: binary_search
use constants, only: ZERO, HALF, SUBSHELLS
use constants
use dict_header, only: DictIntInt, DictCharInt
use endf_header, only: Tabulated1D
use hdf5_interface
@ -75,9 +75,9 @@ module photon_header
type Bremsstrahlung
integer :: i_material ! Index in materials array
real(8), allocatable :: yield(:) ! Photon number yield
real(8), allocatable :: dcs(:,:) ! Scaled bremsstrahlung DCS
real(8), allocatable :: cdf(:,:) ! Bremsstrahlung energy CDF
real(8), allocatable :: yield(:) ! Photon number yield
real(8), allocatable :: dcs(:,:) ! Bremsstrahlung scaled DCS
real(8), allocatable :: cdf(:,:) ! Bremsstrahlung energy CDF
contains
procedure :: init => bremsstrahlung_init
@ -366,10 +366,6 @@ contains
this % pair_production_total = -500.0_8
end where
if (electron_treatment == ELECTRON_TTB) then
ttb_energy_electron = log(ttb_energy_electron)
end if
end subroutine photon_from_hdf5
subroutine bremsstrahlung_init(this, i_material)
@ -383,16 +379,17 @@ contains
real(8) :: c
real(8) :: k, k_l, k_r, k_c
real(8) :: x_l, x_r, x_c
real(8) :: awr
real(8) :: density
real(8) :: Z_eq_sq
real(8) :: beta
real(8) :: atom_fraction
real(8) :: atom_sum
real(8) :: mass_sum
real(8), allocatable :: 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
@ -403,29 +400,50 @@ contains
! Allocate and initialize arrays
n_k = size(ttb_energy_photon)
n_e = size(ttb_energy_electron)
allocate(atom_fraction(mat % n_nuclides))
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
! Calculate the "equivalent" atomic number Zeq and the mass fraction of
! each element
! Calculate the "equivalent" atomic number Zeq, the atomic fraction and the
! mass fraction of each element, and the material density in atom/b-cm
Z_eq_sq = ZERO
do i = 1, mat % n_nuclides
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
awr = nuclides(mat % nuclide(i)) % awr
! Given atom percent
if (mat % atom_density(1) > ZERO) then
atom_fraction(i) = mat % atom_density(i)
mass_fraction(i) = mat % atom_density(i) * awr
! Given weight percent
else
atom_fraction(i) = -mat % atom_density(i) / awr
mass_fraction(i) = -mat % atom_density(i)
end if
Z_eq_sq = Z_eq_sq + atom_fraction(i) * nuclides(mat % nuclide(i)) % Z**2
end do
mass_fraction = mass_fraction / sum(mass_fraction)
atom_sum = sum(atom_fraction)
mass_sum = sum(mass_fraction)
! Given material density in g/cm^3
if (mat % density < ZERO) then
density = -mat % density * N_AVOGADRO / MASS_NEUTRON * (atom_sum / mass_sum)
! Given material density in atom/b-cm
else
density = mat % density
end if
Z_eq_sq = Z_eq_sq / atom_sum
atom_fraction = atom_fraction / atom_sum
mass_fraction = mass_fraction / mass_sum
! Calculate the molecular DCS and the molecular total stopping power using
! Bragg's additivity rule. Note: the collision stopping power cannot be
@ -440,13 +458,10 @@ contains
! Get pointer to current element
elm => elements(mat % element(i))
! 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 + atom_fraction * elm % Z**2 / Z_eq_sq * elm % dcs
this % dcs = this % dcs + atom_fraction(i) * elm % Z**2 / Z_eq_sq * elm % dcs
! Accumulate material total stopping power
stopping_power = stopping_power + mass_fraction(i) * elm % density * &
@ -455,11 +470,11 @@ contains
! Calculate inverse bremsstrahlung mean free path
do i = 1, n_e
e = exp(ttb_energy_electron(i))
e = ttb_energy_electron(i)
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_8)) / (e + MASS_ELECTRON/1.e6_8)
beta = sqrt(e*(e + TWO*MASS_ELECTRON)) / (e + MASS_ELECTRON)
! Integration lower bound
k_c = energy_cutoff(PHOTON) / e
@ -477,31 +492,36 @@ contains
! the DCS at the cutoff energy
x_c = (x_l * (k_r - k_c) + x_r * (k_c - k_l)) / (k_r - k_l)
! Integrate using the trapezoidal rule in log-log space
c = HALF * (log(k_r) - log(k_c)) * (x_c + x_r)
! Calculate the CDF
c = x_r - x_c + (log(k_r)-log(k_c))*(x_c - (k_c*(x_r-x_c)/(k_r-k_c)))
this % cdf(i_k,i) = c
do j = i_k, n_k - 1
c = c + HALF * (log(ttb_energy_photon(j+1)) - &
log(ttb_energy_photon(j))) * (this % dcs(j,i) + this % dcs(j+1,i))
k_l = ttb_energy_photon(j)
k_r = ttb_energy_photon(j+1)
x_l = this % dcs(j, i)
x_r = this % dcs(j+1, i)
c = c + x_r - x_l + (log(k_r)-log(k_l))*(x_l - (k_l*(x_r-x_l)/(k_r-k_l)))
this % cdf(j+1,i) = c
end do
! Calculate the inverse bremsstrahlung mean free path
mfp_inv(i) = c * mat % density * Z_eq_sq / beta**2 * 1.0e-3_8
mfp_inv(i) = c * density * Z_eq_sq / beta**2 * 1.0e-3_8
end do
! Calculate photon number yield
x = exp(ttb_energy_electron)
y = mfp_inv / stopping_power
call spline(x, y, z, n_e)
mfp_inv(:) = mfp_inv(:) / stopping_power(:)
call spline(ttb_energy_electron, mfp_inv, z, n_e)
do i = 1, n_e
this % yield(i) = spline_integrate(x, y, z, n_e, energy_cutoff(PHOTON), x(i))
this % yield(i) = spline_integrate(ttb_energy_electron, mfp_inv, z, &
n_e, energy_cutoff(PHOTON), ttb_energy_electron(i))
end do
! Use logarithm of number yield since it is log-log interpolated
this % yield = log(this % yield)
where (this % yield > ZERO)
this % yield = log(this % yield)
end where
deallocate(mass_fraction, stopping_power, mfp_inv, x, y, z)
deallocate(atom_fraction, mass_fraction, stopping_power, mfp_inv, z)
end subroutine bremsstrahlung_init

34
src/photon_physics.F90
View file

@ -7,6 +7,7 @@ module photon_physics
compton_profile_pz, ttb_energy_electron, &
ttb_energy_photon, ttb
use random_lcg, only: prn
use settings
contains
@ -374,25 +375,31 @@ contains
! THICK_TARGET_BREMSSTRAHLUNG
!===============================================================================
subroutine thick_target_bremsstrahlung(p)
subroutine thick_target_bremsstrahlung(p, E_lost)
type(Particle), intent(inout) :: p
real(8), intent(inout) :: E_lost
integer :: i, j
integer :: i_e, i_k
integer :: n
integer :: n_e, n_k
real(8) :: c_max
real(8) :: f
real(8) :: w
real(8) :: r
real(8) :: e, e_l, e_r
real(8) :: y, y_l, y_r
real(8) :: k, k_l, k_r
real(8) :: k, k_l, k_r, k_c
real(8) :: x, x_l, x_r
type(Bremsstrahlung), pointer :: mat
if (p % E < energy_cutoff(PHOTON)) return
! Get bremsstrahlung data for this material
mat => ttb(p % material)
k_c = energy_cutoff(PHOTON) / p % E
e = log(p % E)
n_e = size(ttb_energy_electron)
n_k = size(ttb_energy_photon)
@ -411,19 +418,20 @@ contains
y = exp((y_l * (e_r - e) + y_r * (e - e_l)) / (e_r - e_l))
! Sample number of secondary bremsstrahlung photons
n = floor(y + prn())
n = int(y + prn())
! 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
w = (e_r - e) / (e_r - e_l)
f = (e_r - e) / (e_r - e_l)
E_lost = ZERO
! 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() > w) then
if (prn() > f) then
i_e = i_e + 1
end if
@ -442,8 +450,12 @@ contains
k_r = ttb_energy_photon(i_k+1)
x_l = mat % dcs(i_k, i_e)
x_r = mat % dcs(i_k+1, i_e)
if (k_l < k_c) then
x_l = (x_l * (k_r - k_c) + x_r * (k_c - k_l)) / (k_r - k_l)
k_l = k_c
end if
! Sample the reduced photon energy k from the distribution k^-1 on the
! Sample the reduced photon energy k from the distribution 1/k on the
! interval (k(i), k(i+1))
k = k_l * (k_r / k_l)**r
@ -454,9 +466,13 @@ contains
if (prn() * max(x_l, x_r) < x) exit
end do
w = k * p % E
E_lost = E_lost + w
if (w < energy_cutoff(PHOTON)) cycle
! Create secondary photon
call p % create_secondary(p % coord(1) % uvw, k * p % E, PHOTON, &
run_ce=.true.)
call p % create_secondary(p % coord(1) % uvw, w, PHOTON, run_ce=.true.)
end do
end subroutine thick_target_bremsstrahlung

17
src/physics.F90
View file

@ -302,10 +302,12 @@ contains
subroutine sample_electron_reaction(p)
type(Particle), intent(inout) :: p
real(8) :: E_lost ! energy lost to bremsstrahlung photons
! TODO: create reaction types
if (electron_treatment == ELECTRON_TTB) then
call thick_target_bremsstrahlung(p)
call thick_target_bremsstrahlung(p, E_lost)
end if
p % E = ZERO
@ -324,16 +326,17 @@ contains
subroutine sample_positron_reaction(p)
type(Particle), intent(inout) :: p
real(8) :: mu ! scattering cosine
real(8) :: phi ! azimuthal angle
real(8) :: uvw(3) ! new direction
real(8) :: mu ! scattering cosine
real(8) :: phi ! azimuthal angle
real(8) :: uvw(3) ! new direction
real(8) :: E_lost ! energy lost to bremsstrahlung photons
! TODO: create reaction types
if (electron_treatment == ELECTRON_TTB) then
! TODO: implement thick-target bremsstrahlung model
call fatal_error("Thick-target bremsstrahlung treatment of positrons &
&is not yet implemented.")
! TODO: implement thick-target bremsstrahlung model for positrons
call thick_target_bremsstrahlung(p, E_lost)
end if
! Sample angle isotropically

2
src/tallies/tally_header.F90
View file

@ -339,6 +339,8 @@ contains
j = FILTER_AZIMUTHAL
type is (EnergyFunctionFilter)
j = FILTER_ENERGYFUNCTION
type is (ParticleFilter)
j = FILTER_PARTICLE
end select
this % find_filter(j) = i
end do