diff --git a/src/geometry.F90 b/src/geometry.F90 index 0bacd7b5d1..10474848fa 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -206,8 +206,11 @@ contains end do j = p % n_coord - ! set size of list to search + ! Determine universe (if not set, use root universe i_universe = p % coord(j) % universe + if (i_universe == NONE) i_universe = root_universe + + ! set size of list to search if (present(search_cells)) then use_search_cells = .true. n = size(search_cells) diff --git a/src/initialize.F90 b/src/initialize.F90 index 269f7752ab..286ff97401 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -11,8 +11,7 @@ module initialize use constants use set_header, only: SetInt use error, only: fatal_error, warning, write_message - use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,& - root_universe + use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice use hdf5_interface, only: file_open, read_attribute, file_close, & hdf5_bank_t, hdf5_integer8_t use input_xml, only: read_input_xml diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cc0680ff12..973bb0e4dc 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4362,7 +4362,7 @@ contains ! Generate material bremsstrahlung data if (photon_transport .and. electron_treatment == ELECTRON_TTB) then - call ttb(i) % init(i) + call bremsstrahlung_init(ttb(i), i) end if end do diff --git a/src/material_header.F90 b/src/material_header.F90 index 43fc6679a3..70bf4bd21c 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -7,7 +7,7 @@ module material_header use error use nuclide_header use particle_header, only: Particle - use photon_header, only: micro_photon_xs, elements + use photon_header use sab_header use simulation_header, only: log_spacing use stl_vector, only: VectorReal, VectorInt @@ -16,6 +16,7 @@ module material_header implicit none private + public :: bremsstrahlung_init public :: free_memory_material public :: openmc_extend_materials public :: openmc_get_material_index @@ -324,7 +325,7 @@ contains ! If particle energy is greater than the highest energy for the ! S(a,b) table, then don't use the S(a,b) table - if (E > sab_tables(i_sab) % data(1) % threshold_inelastic) then + if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then i_sab = 0 end if @@ -383,7 +384,7 @@ contains subroutine calculate_photon_xs(this, p) class(Material), intent(in) :: this - type(Particle), intent(inout) :: p + type(Particle), intent(in) :: p integer :: i ! loop index over nuclides integer :: i_element ! index into elements array @@ -698,4 +699,166 @@ contains end function openmc_material_set_densities + subroutine bremsstrahlung_init(this, i_material) + class(Bremsstrahlung), intent(inout) :: this + integer, intent(in) :: i_material + + integer :: i, j + integer :: i_k + integer :: n_e, n_k + real(8) :: e + 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) :: density_gpcc + real(8) :: Z_eq_sq + real(8) :: beta + 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 :: z(:) + type(Material), pointer :: mat + type(PhotonInteraction), pointer :: elm + + ! Get pointer to this material + mat => materials(i_material) + this % i_material = i_material + + ! Allocate and initialize arrays + n_k = size(ttb_k_grid) + n_e = size(ttb_e_grid) + 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(z(n_e)) + stopping_power(:) = ZERO + mfp_inv(:) = ZERO + this % dcs(:,:) = ZERO + this % cdf(:,:) = ZERO + + ! Calculate the "equivalent" atomic number Zeq, the atomic fraction and the + ! mass fraction of each element, and the material density in atom/b-cm and + ! in g/cm^3 + Z_eq_sq = ZERO + do i = 1, mat % n_nuclides + 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 + 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 * (atom_sum / mass_sum) * N_AVOGADRO / MASS_NEUTRON + density_gpcc = -mat % density + + ! Given material density in atom/b-cm + else + density = mat % density + density_gpcc = mat % density * (mass_sum / atom_sum) * MASS_NEUTRON / & + N_AVOGADRO + 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 + ! accurately calculated using Bragg's additivity rule since the mean + ! excitation energies and the density effect corrections cannot simply be + ! summed together. Bragg's additivity rule fails especially when a + ! higher-density compound is composed of elements that are in lower-density + ! form at normal temperature and pressure (at which the NIST stopping + ! powers are given). It will be used to approximate the collision stopping + ! powers for now, but should be fixed in the future. + do i = 1, mat % n_nuclides + ! Get pointer to current element + elm => elements(mat % element(i)) + + ! 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(i) * elm % Z**2 / Z_eq_sq * elm % dcs + + ! Accumulate material total stopping power + stopping_power = stopping_power + mass_fraction(i) * density_gpcc * & + (elm % stopping_power_collision + elm % stopping_power_radiative) + end do + + ! Calculate inverse bremsstrahlung mean free path + do i = 1, n_e + e = ttb_e_grid(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)) / (e + MASS_ELECTRON) + + ! Integration lower bound + k_c = energy_cutoff(PHOTON) / e + + ! Find the upper bounding index of the reduced photon cutoff energy + i_k = binary_search(ttb_k_grid, n_k, k_c) + 1 + + ! Get the interpolation bounds + k_l = ttb_k_grid(i_k-1) + k_r = ttb_k_grid(i_k) + x_l = this % dcs(i_k-1, i) + x_r = this % dcs(i_k, i) + + ! Use linear interpolation in reduced photon energy k to find value of + ! the DCS at the cutoff energy + x_c = (x_l * (k_r - k_c) + x_r * (k_c - k_l)) / (k_r - k_l) + + ! Calculate the CDF using the trapezoidal rule in log-log space + c = HALF * (log(k_r) - log(k_c)) * (x_c + x_r) + this % cdf(i_k,i) = c + do j = i_k, n_k - 1 + c = c + HALF * (log(ttb_k_grid(j+1)) - log(ttb_k_grid(j))) * & + (this % dcs(j,i) + this % dcs(j+1,i)) + this % cdf(j+1,i) = c + end do + + ! Calculate the inverse bremsstrahlung mean free path + mfp_inv(i) = c * density * Z_eq_sq / beta**2 * 1.0e-3_8 + end do + + ! Calculate photon number yield + mfp_inv(:) = mfp_inv(:) / stopping_power(:) + call spline(ttb_e_grid, mfp_inv, z, n_e) + do i = 1, n_e + this % yield(i) = spline_integrate(ttb_e_grid, mfp_inv, z, n_e, & + energy_cutoff(PHOTON), ttb_e_grid(i)) + end do + + ! Use logarithm of number yield since it is log-log interpolated + where (this % yield > ZERO) + this % yield = log(this % yield) + end where + + deallocate(atom_fraction, mass_fraction, stopping_power, mfp_inv, z) + + end subroutine bremsstrahlung_init + end module material_header diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 0c0a3fc24f..08e69c1d83 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -5,7 +5,6 @@ module particle_header use bank_header, only: Bank, source_bank use constants use error, only: fatal_error, warning - use geometry_header, only: root_universe use hdf5_interface use settings use simulation_header @@ -177,14 +176,13 @@ contains end if n = this % n_secondary + 1 - this % secondary_bank(n) % particle = tyoe + 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 % secondary_bank(n) % E = E - this % secondary_bank(this % n_secondary) % E = this % 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 @@ -224,7 +222,7 @@ contains this % g = NONE ! Set up base level coordinates - this % coord(1) % universe = root_universe + this % coord(1) % universe = NONE this % n_coord = 1 this % last_n_coord = 1 diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 24f31f3dc9..27fe62f189 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -8,7 +8,6 @@ module photon_header use endf_header, only: Tabulated1D use hdf5_interface use math, only: spline, spline_integrate - use material_header, only: Material, materials use nuclide_header, only: nuclides use settings @@ -79,9 +78,6 @@ module photon_header 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 end type Bremsstrahlung type(PhotonInteraction), allocatable, target :: elements(:) ! Photon cross sections @@ -432,168 +428,6 @@ contains xs % last_E = E - end subroutine calculate_element_xs - - subroutine bremsstrahlung_init(this, i_material) - class(Bremsstrahlung), intent(inout) :: this - integer, intent(in) :: i_material - - integer :: i, j - integer :: i_k - integer :: n_e, n_k - real(8) :: e - 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) :: density_gpcc - real(8) :: Z_eq_sq - real(8) :: beta - 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 :: z(:) - type(Material), pointer :: mat - type(PhotonInteraction), pointer :: elm - - ! Get pointer to this material - mat => materials(i_material) - this % i_material = i_material - - ! Allocate and initialize arrays - n_k = size(ttb_k_grid) - n_e = size(ttb_e_grid) - 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(z(n_e)) - stopping_power(:) = ZERO - mfp_inv(:) = ZERO - this % dcs(:,:) = ZERO - this % cdf(:,:) = ZERO - - ! Calculate the "equivalent" atomic number Zeq, the atomic fraction and the - ! mass fraction of each element, and the material density in atom/b-cm and - ! in g/cm^3 - Z_eq_sq = ZERO - do i = 1, mat % n_nuclides - 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 - 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 * (atom_sum / mass_sum) * N_AVOGADRO / MASS_NEUTRON - density_gpcc = -mat % density - - ! Given material density in atom/b-cm - else - density = mat % density - density_gpcc = mat % density * (mass_sum / atom_sum) * MASS_NEUTRON / & - N_AVOGADRO - 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 - ! accurately calculated using Bragg's additivity rule since the mean - ! excitation energies and the density effect corrections cannot simply be - ! summed together. Bragg's additivity rule fails especially when a - ! higher-density compound is composed of elements that are in lower-density - ! form at normal temperature and pressure (at which the NIST stopping - ! powers are given). It will be used to approximate the collision stopping - ! powers for now, but should be fixed in the future. - do i = 1, mat % n_nuclides - ! Get pointer to current element - elm => elements(mat % element(i)) - - ! 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(i) * elm % Z**2 / Z_eq_sq * elm % dcs - - ! Accumulate material total stopping power - stopping_power = stopping_power + mass_fraction(i) * density_gpcc * & - (elm % stopping_power_collision + elm % stopping_power_radiative) - end do - - ! Calculate inverse bremsstrahlung mean free path - do i = 1, n_e - e = ttb_e_grid(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)) / (e + MASS_ELECTRON) - - ! Integration lower bound - k_c = energy_cutoff(PHOTON) / e - - ! Find the upper bounding index of the reduced photon cutoff energy - i_k = binary_search(ttb_k_grid, n_k, k_c) + 1 - - ! Get the interpolation bounds - k_l = ttb_k_grid(i_k-1) - k_r = ttb_k_grid(i_k) - x_l = this % dcs(i_k-1, i) - x_r = this % dcs(i_k, i) - - ! Use linear interpolation in reduced photon energy k to find value of - ! the DCS at the cutoff energy - x_c = (x_l * (k_r - k_c) + x_r * (k_c - k_l)) / (k_r - k_l) - - ! Calculate the CDF using the trapezoidal rule in log-log space - c = HALF * (log(k_r) - log(k_c)) * (x_c + x_r) - this % cdf(i_k,i) = c - do j = i_k, n_k - 1 - c = c + HALF * (log(ttb_k_grid(j+1)) - log(ttb_k_grid(j))) * & - (this % dcs(j,i) + this % dcs(j+1,i)) - this % cdf(j+1,i) = c - end do - - ! Calculate the inverse bremsstrahlung mean free path - mfp_inv(i) = c * density * Z_eq_sq / beta**2 * 1.0e-3_8 - end do - - ! Calculate photon number yield - mfp_inv(:) = mfp_inv(:) / stopping_power(:) - call spline(ttb_e_grid, mfp_inv, z, n_e) - do i = 1, n_e - this % yield(i) = spline_integrate(ttb_e_grid, mfp_inv, z, n_e, & - energy_cutoff(PHOTON), ttb_e_grid(i)) - end do - - ! Use logarithm of number yield since it is log-log interpolated - where (this % yield > ZERO) - this % yield = log(this % yield) - end where - - deallocate(atom_fraction, mass_fraction, stopping_power, mfp_inv, z) - - end subroutine bremsstrahlung_init + end subroutine photon_calculate_xs end module photon_header diff --git a/src/physics.F90 b/src/physics.F90 index 84cac725e6..0f27a5361c 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -439,7 +439,7 @@ contains ! Check to make sure that a nuclide was sampled if (i > mat % n_nuclides) then - call write_particle_restart(p) + call p % write_restart() call fatal_error("Did not sample any element during collision.") end if diff --git a/src/tracking.F90 b/src/tracking.F90 index 189e79abb0..b528f571a2 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -114,8 +114,7 @@ contains ! If the material is the same as the last material and the ! temperature hasn't changed, we don't need to lookup cross ! sections again. - call materials(p % material) % calculate_xs(p % E, p % sqrtkT, & - micro_xs, nuclides, material_xs) + call materials(p % material) % calculate_xs(p) end if else ! Get the MG data