diff --git a/openmc/data/photon.py b/openmc/data/photon.py index d17aaccf7..be04ec911 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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 diff --git a/openmc/filter.py b/openmc/filter.py index 9d2f9c697..ef22090a8 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -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 diff --git a/openmc/settings.py b/openmc/settings.py index 760155015..563d2fefb 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -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) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a145d41bb..e4c68c3ea 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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, & diff --git a/src/photon_header.F90 b/src/photon_header.F90 index c00044d4c..3222dff00 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -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 diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index bba86fcf6..127ce0c26 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -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 diff --git a/src/physics.F90 b/src/physics.F90 index 8364962ee..49af53cb6 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index df9d935e9..b9e994750 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -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