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

Fix compiler errors and circular dependencies

Browse source
This commit is contained in:
Paul Romano 2018-03-11 13:40:01 -05:00
parent 460452d4f7
commit 712a904636
8 changed files with 178 additions and 182 deletions
Download patch file Download diff file Expand all files Collapse all files

5
src/geometry.F90
View file

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

3
src/initialize.F90
View file

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

2
src/input_xml.F90
View file

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

169
src/material_header.F90
View file

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

8
src/particle_header.F90
View file

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

168
src/photon_header.F90
View file

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

2
src/physics.F90
View file

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

3
src/tracking.F90
View file

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