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

Fixing bugs, partially there

Browse source
This commit is contained in:
Adam Nelson 2016-02-03 10:36:16 -05:00
parent e24a692c79
commit 63744db909
12 changed files with 113 additions and 125 deletions
Download patch file Download diff file Expand all files Collapse all files

8
src/ace.F90
View file

@ -1,6 +1,6 @@
module ace
use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing
use ace_header, only: Reaction
use constants
use distribution_univariate, only: Uniform, Equiprobable, Tabular
use endf, only: is_fission, is_disappearance
@ -11,13 +11,15 @@ module ace
use global
use list_header, only: ListInt
use material_header, only: Material
use nuclide_header
use output, only: write_message
use sab_header
use set_header, only: SetChar
use secondary_header, only: AngleEnergy
use secondary_correlated, only: CorrelatedAngleEnergy
use secondary_kalbach, only: KalbachMann
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
use string, only: to_str, to_lower
use simple_string, only: to_str, to_lower
implicit none
@ -1472,7 +1474,7 @@ contains
!===============================================================================
subroutine generate_nu_fission(nuc)
type(Nuclide), intent(inout) :: nuc
type(Nuclide_CE), intent(inout) :: nuc
integer :: i ! index on nuclide energy grid
real(8) :: E ! energy

6
src/distribution_multivariate.F90
View file

@ -1,9 +1,9 @@
module distribution_multivariate
use constants, only: ONE, TWO, PI
use constants, only: ONE, TWO, PI
use distribution_univariate, only: Distribution
use math, only: rotate_angle
use random_lcg, only: prn
use random_lcg, only: prn
use spectra, only: rotate_angle
implicit none

10
src/distribution_univariate.F90
View file

@ -1,11 +1,11 @@
module distribution_univariate
use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, &
use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, &
MAX_LINE_LEN, MAX_WORD_LEN
use error, only: fatal_error
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
use string, only: to_lower
use error, only: fatal_error
use random_lcg, only: prn
use spectra, only: maxwell_spectrum, watt_spectrum
use simple_string, only: to_lower
use xml_interface
implicit none

2
src/energy_distribution.F90
View file

@ -3,9 +3,9 @@ module energy_distribution
use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tab1
use interpolation, only: interpolate_tab1
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
use search, only: binary_search
use spectra, only: maxwell_spectrum, watt_spectrum
!===============================================================================
! ENERGYDISTRIBUTION (abstract) defines an energy distribution that is a

6
src/interpolation.F90
View file

@ -1,9 +1,9 @@
module interpolation
use constants
use endf_header, only: Tab1
use search, only: binary_search
use string, only: to_str
use endf_header, only: Tab1
use search, only: binary_search
use simple_string, only: to_str
implicit none

42
src/nuclide_header.F90
View file

@ -94,7 +94,7 @@ module nuclide_header
integer :: n_precursor ! # of delayed neutron precursors
real(8), allocatable :: nu_d_data(:)
real(8), allocatable :: nu_d_precursor_data(:)
type(DistEnergy), pointer :: nu_d_edist(:) => null()
type(AngleEnergyContainer), allocatable :: nu_d_edist(:)
! Unresolved resonance data
logical :: urr_present
@ -308,16 +308,7 @@ module nuclide_header
integer :: i ! Loop counter
if (associated(this % nu_d_edist)) then
do i = 1, size(this % nu_d_edist)
call this % nu_d_edist(i) % clear()
end do
deallocate(this % nu_d_edist)
end if
if (associated(this % urr_data)) then
deallocate(this % urr_data)
end if
if (associated(this % urr_data)) deallocate(this % urr_data)
if (allocated(this % reactions)) then
do i = 1, size(this % reactions)
@ -389,7 +380,6 @@ module nuclide_header
subroutine nuclide_ce_print(this, unit)
class(Nuclide_CE), intent(in) :: this
type(Nuclide), intent(in) :: nuc
integer, intent(in), optional :: unit
integer :: i ! loop index over nuclides
@ -410,36 +400,36 @@ module nuclide_header
size_xs = 0
! Basic nuclide information
write(unit_,*) 'Nuclide ' // trim(nuc % name)
write(unit_,*) ' zaid = ' // trim(to_str(nuc % zaid))
write(unit_,*) ' awr = ' // trim(to_str(nuc % awr))
write(unit_,*) ' kT = ' // trim(to_str(nuc % kT))
write(unit_,*) ' # of grid points = ' // trim(to_str(nuc % n_grid))
write(unit_,*) ' Fissionable = ', nuc % fissionable
write(unit_,*) ' # of fission reactions = ' // trim(to_str(nuc % n_fission))
write(unit_,*) ' # of reactions = ' // trim(to_str(nuc % n_reaction))
write(unit_,*) 'Nuclide ' // trim(this % name)
write(unit_,*) ' zaid = ' // trim(to_str(this % zaid))
write(unit_,*) ' awr = ' // trim(to_str(this % awr))
write(unit_,*) ' kT = ' // trim(to_str(this % kT))
write(unit_,*) ' # of grid points = ' // trim(to_str(this % n_grid))
write(unit_,*) ' Fissionable = ', this % fissionable
write(unit_,*) ' # of fission reactions = ' // trim(to_str(this % n_fission))
write(unit_,*) ' # of reactions = ' // trim(to_str(this % n_reaction))
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, nuc % n_reaction
associate (rxn => nuc % reactions(i))
do i = 1, this % n_reaction
associate (rxn => this % reactions(i))
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
rxn % threshold
! Accumulate data size
size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8
size_xs = size_xs + (this % n_grid - rxn%threshold + 1) * 8
end associate
end do
! Add memory required for summary reactions (total, absorption, fission,
! nu-fission)
size_xs = 8 * nuc % n_grid * 4
size_xs = 8 * this % n_grid * 4
! Write information about URR probability tables
size_urr = 0
if (nuc % urr_present) then
urr => nuc % urr_data
if (this % urr_present) then
urr => this % urr_data
write(unit_,*) ' Unresolved resonance probability table:'
write(unit_,*) ' # of energies = ' // trim(to_str(urr % n_energy))
write(unit_,*) ' # of probabilities = ' // trim(to_str(urr % n_prob))

2
src/output.F90
View file

@ -2,7 +2,7 @@ module output
use, intrinsic :: ISO_FORTRAN_ENV
use ace_header, only: Nuclide, Reaction, UrrData
use ace_header, only: Reaction, UrrData
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning

4
src/physics.F90
View file

@ -1,6 +1,6 @@
module physics
use ace_header, only: Reaction, DistEnergy
use ace_header, only: Reaction
use constants
use cross_section, only: elastic_xs_0K
use endf, only: reaction_name
@ -1343,7 +1343,7 @@ contains
p % wgt = yield * p % wgt
else
do i = 1, rxn % multiplicity - 1
call p % create_secondary(p % coord(1) % uvw, NEUTRON)
call p % create_secondary(p % coord(1) % uvw, NEUTRON, run_CE=.True.)
end do
end if

53
src/physics_common.F90
View file

@ -30,57 +30,4 @@ contains
end subroutine russian_roulette
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done
! with direct sampling rather than rejection as is done in MCNP and SERPENT.
!===============================================================================
function rotate_angle(uvw0, mu, phi) result(uvw)
real(8), intent(in) :: uvw0(3) ! directional cosine
real(8), intent(in) :: mu ! cosine of angle in lab or CM
real(8), optional :: phi ! azimuthal angle
real(8) :: uvw(3) ! rotated directional cosine
real(8) :: phi_ ! azimuthal angle
real(8) :: sinphi ! sine of azimuthal angle
real(8) :: cosphi ! cosine of azimuthal angle
real(8) :: a ! sqrt(1 - mu^2)
real(8) :: b ! sqrt(1 - w^2)
real(8) :: u0 ! original cosine in x direction
real(8) :: v0 ! original cosine in y direction
real(8) :: w0 ! original cosine in z direction
! Copy original directional cosines
u0 = uvw0(1)
v0 = uvw0(2)
w0 = uvw0(3)
! Sample azimuthal angle in [0,2pi) if none provided
if (present(phi)) then
phi_ = phi
else
phi_ = TWO * PI * prn()
end if
! Precompute factors to save flops
sinphi = sin(phi_)
cosphi = cos(phi_)
a = sqrt(max(ZERO, ONE - mu*mu))
b = sqrt(max(ZERO, ONE - w0*w0))
! Need to treat special case where sqrt(1 - w**2) is close to zero by
! expanding about the v component rather than the w component
if (b > 1e-10) then
uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
uvw(3) = mu*w0 - a*b*cosphi
else
b = sqrt(ONE - v0*v0)
uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
uvw(2) = mu*v0 - a*b*cosphi
uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
end if
end function rotate_angle
end module physics_common

6
src/physics_mg.F90
View file

@ -16,6 +16,7 @@ module physics_mg
use random_lcg, only: prn
use scattdata_header
use simple_string, only: to_str
use spectra, only: rotate_angle
implicit none
@ -71,7 +72,6 @@ contains
! absorption (including fission)
if (mat % fissionable) then
call sample_fission(i_nuclide, i_reaction)
if (run_mode == MODE_EIGENVALUE) then
call create_fission_sites(p, fission_bank, n_bank)
elseif (run_mode == MODE_FIXEDSOURCE) then
@ -234,10 +234,6 @@ contains
! Bank source neutrons
if (nu == 0 .or. size_bank == size(bank_array)) return
! Initialize counter of delayed neutrons encountered for each delayed group
! to zero.
nu_d(:) = 0
p % fission = .true. ! Fission neutrons will be banked
do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4)
! Bank source neutrons by copying particle data

43
src/sab_header.F90
View file

@ -66,19 +66,18 @@ module sab_header
contains
!===============================================================================
! PRINT_SAB_TABLE displays information about a S(a,b) table containing data
! describing thermal scattering from bound materials such as hydrogen in water.
!===============================================================================
subroutine print_sab_table(this, unit)
type(SAlphaBeta), intent(in) :: sab
class(SAlphaBeta), intent(in) :: this
integer, intent(in), optional :: unit
integer :: size_sab ! memory used by S(a,b) table
integer :: unit_ ! unit to write to
integer :: i ! Loop counter for parsing through sab % zaid
integer :: i ! Loop counter for parsing through this % zaid
integer :: char_count ! Counter for the number of characters on a line
! set default unit for writing information
@ -89,11 +88,11 @@ module sab_header
end if
! Basic S(a,b) table information
write(unit_,*) 'S(a,b) Table ' // trim(sab % name)
write(unit_,*) 'S(a,b) Table ' // trim(this % name)
write(unit_,'(A)',advance="no") ' zaids = '
! Initialize the counter based on the above string
char_count = 11
do i = 1, sab % n_zaid
do i = 1, this % n_zaid
! Deal with a line thats too long
if (char_count >= 73) then ! 73 = 80 - (5 ZAID chars + 1 space + 1 comma)
! End the line
@ -103,44 +102,44 @@ module sab_header
! reset the counter to 11
char_count = 11
end if
if (i < sab % n_zaid) then
if (i < this % n_zaid) then
! Include a comma
write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i))) // ", "
char_count = char_count + len(trim(to_str(sab % zaid(i)))) + 2
write(unit_,'(A)',advance="no") trim(to_str(this % zaid(i))) // ", "
char_count = char_count + len(trim(to_str(this % zaid(i)))) + 2
else
! Don't include a comma, since we are all done
write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i)))
write(unit_,'(A)',advance="no") trim(to_str(this % zaid(i)))
end if
end do
write(unit_,*) "" ! Move to next line
write(unit_,*) ' awr = ' // trim(to_str(sab % awr))
write(unit_,*) ' kT = ' // trim(to_str(sab % kT))
write(unit_,*) ' awr = ' // trim(to_str(this % awr))
write(unit_,*) ' kT = ' // trim(to_str(this % kT))
! Inelastic data
write(unit_,*) ' # of Incoming Energies (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_e_in))
trim(to_str(this % n_inelastic_e_in))
write(unit_,*) ' # of Outgoing Energies (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_e_out))
trim(to_str(this % n_inelastic_e_out))
write(unit_,*) ' # of Outgoing Angles (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_mu))
trim(to_str(this % n_inelastic_mu))
write(unit_,*) ' Threshold for Inelastic = ' // &
trim(to_str(sab % threshold_inelastic))
trim(to_str(this % threshold_inelastic))
! Elastic data
if (sab % n_elastic_e_in > 0) then
if (this % n_elastic_e_in > 0) then
write(unit_,*) ' # of Incoming Energies (Elastic) = ' // &
trim(to_str(sab % n_elastic_e_in))
trim(to_str(this % n_elastic_e_in))
write(unit_,*) ' # of Outgoing Angles (Elastic) = ' // &
trim(to_str(sab % n_elastic_mu))
trim(to_str(this % n_elastic_mu))
write(unit_,*) ' Threshold for Elastic = ' // &
trim(to_str(sab % threshold_elastic))
trim(to_str(this % threshold_elastic))
end if
! Determine memory used by S(a,b) table and write out
size_sab = 8 * (sab % n_inelastic_e_in * (2 + sab % n_inelastic_e_out * &
(1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * &
(2 + sab % n_elastic_mu))
size_sab = 8 * (this % n_inelastic_e_in * (2 + this % n_inelastic_e_out * &
(1 + this % n_inelastic_mu)) + this % n_elastic_e_in * &
(2 + this % n_elastic_mu))
write(unit_,*) ' Memory Used = ' // trim(to_str(size_sab)) // ' bytes'
! Blank line at end

56
src/spectra.F90
View file

@ -1,6 +1,6 @@
module spectra
use constants, only: ONE, TWO, PI
use constants, only: ZERO, ONE, TWO, PI
use random_lcg, only: prn
implicit none
@ -55,4 +55,58 @@ contains
end function watt_spectrum
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done
! with direct sampling rather than rejection as is done in MCNP and SERPENT.
!===============================================================================
function rotate_angle(uvw0, mu, phi) result(uvw)
real(8), intent(in) :: uvw0(3) ! directional cosine
real(8), intent(in) :: mu ! cosine of angle in lab or CM
real(8), optional :: phi ! azimuthal angle
real(8) :: uvw(3) ! rotated directional cosine
real(8) :: phi_ ! azimuthal angle
real(8) :: sinphi ! sine of azimuthal angle
real(8) :: cosphi ! cosine of azimuthal angle
real(8) :: a ! sqrt(1 - mu^2)
real(8) :: b ! sqrt(1 - w^2)
real(8) :: u0 ! original cosine in x direction
real(8) :: v0 ! original cosine in y direction
real(8) :: w0 ! original cosine in z direction
! Copy original directional cosines
u0 = uvw0(1)
v0 = uvw0(2)
w0 = uvw0(3)
! Sample azimuthal angle in [0,2pi) if none provided
if (present(phi)) then
phi_ = phi
else
phi_ = TWO * PI * prn()
end if
! Precompute factors to save flops
sinphi = sin(phi_)
cosphi = cos(phi_)
a = sqrt(max(ZERO, ONE - mu*mu))
b = sqrt(max(ZERO, ONE - w0*w0))
! Need to treat special case where sqrt(1 - w**2) is close to zero by
! expanding about the v component rather than the w component
if (b > 1e-10) then
uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
uvw(3) = mu*w0 - a*b*cosphi
else
b = sqrt(ONE - v0*v0)
uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
uvw(2) = mu*v0 - a*b*cosphi
uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
end if
end function rotate_angle
end module spectra