mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 21:55:41 -04:00
Fixed sampling of secondary angle for elastic scattering. Added 2 benchmarks.
This commit is contained in:
Paul Romano
parent
eee77ec2f3
commit
bf2b2df5b3
7 changed files with 150 additions and 51 deletions
12
ChangeLog
12
ChangeLog
|
|
@ -1,3 +1,15 @@
|
|||
2011-03-22 Paul Romano <romano7@mit.edu>
|
||||
|
||||
* benchmarks: Added PU-MET-FAST-002 and PU-MET-FAST-005 benchmark
|
||||
problems.
|
||||
* mpi_routines.f90: Added check for no fission bank sites.
|
||||
* output.f90: Changed error subroutine to work when run in
|
||||
parallel.
|
||||
* physics.f90: Fixed sampling of seconary angle for elastic
|
||||
scattering. Changed sampling of reaction. Changed n_gamma to
|
||||
n_absorption and now includes (n,p), etc. Reverted sample_angle
|
||||
routine to arguments u,v,w,mu.
|
||||
|
||||
2011-03-22 Paul Romano <romano7@mit.edu>
|
||||
|
||||
* cross_section.f90: Calculation of material total cross section
|
||||
|
|
|
|||
|
|
@ -10,7 +10,7 @@ cell 1 0 1 -1
|
|||
|
||||
surface 1 sph 0. 0. 0. 6.3849
|
||||
|
||||
source_box -4 -4 -4 4 4 4
|
||||
source_box -1 -1 -1 1 1 1
|
||||
|
||||
material 1 0.04029014 &
|
||||
94239.03c 0.037047 &
|
||||
|
|
@ -20,5 +20,5 @@ material 1 0.04029014 &
|
|||
|
||||
xs_library endfb7
|
||||
xs_data /opt/serpent/xsdata/endfb7
|
||||
criticality 150 25 1000
|
||||
criticality 3000 20 10000
|
||||
verbosity 7
|
||||
|
|
|
|||
25
benchmarks/pu-met-fast-002
Normal file
25
benchmarks/pu-met-fast-002
Normal file
|
|
@ -0,0 +1,25 @@
|
|||
# ====================================================================
|
||||
# Description: Bare Pu-240 Jezebel Benchmark
|
||||
# Case: PU-MET-FAST-002
|
||||
# Written By: Paul Romano
|
||||
# Date: 3/22/2011
|
||||
# ====================================================================
|
||||
|
||||
# cell no. universe material surfaces
|
||||
cell 1 0 1 -1
|
||||
|
||||
surface 1 sph 0. 0. 0. 6.6595
|
||||
|
||||
source_box -1 -1 -1 1 1 1
|
||||
|
||||
material 1 0.04055292 &
|
||||
94239.03c 0.029934 &
|
||||
94240.03c 0.0078754 &
|
||||
94241.03c 0.0012146 &
|
||||
94242.03c 0.00015672 &
|
||||
31000.03c 0.0013722
|
||||
|
||||
xs_library endfb7
|
||||
xs_data /opt/serpent/xsdata/endfb7
|
||||
criticality 3000 20 10000
|
||||
verbosity 7
|
||||
32
benchmarks/pu-met-fast-005
Normal file
32
benchmarks/pu-met-fast-005
Normal file
|
|
@ -0,0 +1,32 @@
|
|||
# ====================================================================
|
||||
# Description: Unmoderated Plutonium Metal Button Array Benchmark
|
||||
# Case: PU-MET-FAST-005
|
||||
# Written By: Paul Romano
|
||||
# Date: 3/22/2011
|
||||
# ====================================================================
|
||||
|
||||
# cell no. universe material surfaces
|
||||
cell 1 0 1 -1
|
||||
cell 2 0 2 1 -2
|
||||
|
||||
surface 1 sph 0. 0. 0. 5.0419
|
||||
surface 2 sph 0. 0. 0. 9.7409
|
||||
|
||||
source_box -1 -1 -1 1 1 1
|
||||
|
||||
material 1 0.04070346 &
|
||||
94239.03c 0.037291 &
|
||||
94240.03c 0.0019277 &
|
||||
94241.03c 0.00012196 &
|
||||
31000.03c 0.0013628
|
||||
|
||||
material 2 0.06605308 &
|
||||
74000.03c 0.051468 &
|
||||
28000.03c 0.0097124 &
|
||||
29000.03c 0.0040774 &
|
||||
40000.03c 0.00079528
|
||||
|
||||
xs_library endfb7
|
||||
xs_data /opt/serpent/xsdata/endfb7
|
||||
criticality 3000 20 10000
|
||||
verbosity 7
|
||||
|
|
@ -119,6 +119,12 @@ contains
|
|||
call MPI_BCAST(total, 1, MPI_INTEGER, n_procs - 1, &
|
||||
& MPI_COMM_WORLD, ierr)
|
||||
|
||||
! Check if there are no fission sites
|
||||
if (total == 0) then
|
||||
msg = "No fission sites banked!"
|
||||
call error(msg)
|
||||
end if
|
||||
|
||||
allocate(temp_sites(2*work))
|
||||
count = 0 ! Index for source_bank
|
||||
index = 0 ! Index for uid -- must account for all nodes
|
||||
|
|
|
|||
|
|
@ -159,20 +159,21 @@ module output
|
|||
integer :: i
|
||||
|
||||
! Only allow master to print to screen
|
||||
if (.not. master) return
|
||||
if (master) then
|
||||
write(eu, fmt='(1X,A7)', advance='no') 'ERROR: '
|
||||
|
||||
write(eu, fmt='(1X,A7)', advance='no') 'ERROR: '
|
||||
|
||||
n_lines = (len_trim(msg)-1)/72 + 1
|
||||
do i = 1, n_lines
|
||||
if ( i == 1 ) then
|
||||
write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i)
|
||||
else
|
||||
write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i)
|
||||
end if
|
||||
end do
|
||||
write(eu,*)
|
||||
n_lines = (len_trim(msg)-1)/72 + 1
|
||||
do i = 1, n_lines
|
||||
if ( i == 1 ) then
|
||||
write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i)
|
||||
else
|
||||
write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i)
|
||||
end if
|
||||
end do
|
||||
write(eu,*)
|
||||
end if
|
||||
|
||||
! All processors abort
|
||||
call free_memory()
|
||||
|
||||
end subroutine error
|
||||
|
|
|
|||
|
|
@ -22,7 +22,8 @@ contains
|
|||
type(Particle), pointer :: p
|
||||
|
||||
character(250) :: msg
|
||||
integer :: i, surf
|
||||
integer :: i
|
||||
integer :: surf
|
||||
logical :: alive
|
||||
logical :: found_cell
|
||||
|
||||
|
|
@ -184,11 +185,11 @@ contains
|
|||
end do
|
||||
|
||||
! Get table, total xs, interpolation factor
|
||||
density_i = cMaterial%atom_percent(i)*density
|
||||
table => xs_continuous(cMaterial%table(i))
|
||||
Sigma = Sigma_t(i)
|
||||
Sigma = Sigma_t(i) / density_i
|
||||
IE = table%grid_index(p % IE)
|
||||
f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
|
||||
density = cMaterial%atom_percent(i)*density
|
||||
|
||||
! free memory
|
||||
deallocate(Sigma_t)
|
||||
|
|
@ -198,9 +199,9 @@ contains
|
|||
prob = ZERO
|
||||
do i = 1, table%n_reaction
|
||||
rxn => table%reactions(i)
|
||||
if (rxn%MT >= 200) cycle
|
||||
if (rxn%MT >= 200 .or. rxn%MT == 4) cycle
|
||||
if (IE < rxn%IE) cycle
|
||||
prob = prob + density * ((ONE-f)*rxn%sigma(IE-rxn%IE+1) &
|
||||
prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) &
|
||||
& + f*(rxn%sigma(IE-rxn%IE+2)))
|
||||
if (r1 < prob) exit
|
||||
end do
|
||||
|
|
@ -213,17 +214,17 @@ contains
|
|||
! call appropriate subroutine
|
||||
select case (rxn%MT)
|
||||
case (2)
|
||||
call elastic_scatter(p, table%awr)
|
||||
call elastic_scatter(p, table, rxn)
|
||||
case (11, 16, 17, 24, 25, 30, 37, 41, 42)
|
||||
call n_xn(p, table, rxn)
|
||||
case (22, 23, 28, 29, 32:36, 43, 44, 51:91)
|
||||
call inelastic_scatter(p, table, rxn)
|
||||
case (18:21, 38)
|
||||
call n_fission(p, table, rxn)
|
||||
case (102)
|
||||
call n_gamma(p)
|
||||
case (102:117)
|
||||
call n_absorption(p)
|
||||
case default
|
||||
call elastic_scatter(p, table%awr)
|
||||
print *, 'warning: cant simulate reaction: ', rxn % MT
|
||||
end select
|
||||
|
||||
! find energy index, interpolation factor
|
||||
|
|
@ -235,11 +236,13 @@ contains
|
|||
! ELASTIC_SCATTER
|
||||
!=====================================================================
|
||||
|
||||
subroutine elastic_scatter(p, awr)
|
||||
subroutine elastic_scatter(p, table, rxn)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
real(8), intent(in) :: awr
|
||||
type(Particle), pointer :: p
|
||||
type(AceContinuous), pointer :: table
|
||||
type(AceReaction), pointer :: rxn
|
||||
|
||||
real(8) :: awr
|
||||
real(8) :: phi ! azimuthal angle
|
||||
real(8) :: mu ! cosine of polar angle
|
||||
real(8) :: vx, vy, vz
|
||||
|
|
@ -250,12 +253,12 @@ contains
|
|||
integer :: IE
|
||||
|
||||
vel = sqrt(p % E)
|
||||
awr = table % awr
|
||||
|
||||
vx = vel*p % uvw(1)
|
||||
vy = vel*p % uvw(2)
|
||||
vz = vel*p % uvw(3)
|
||||
|
||||
|
||||
vcx = vx/(awr + ONE)
|
||||
vcy = vy/(awr + ONE)
|
||||
vcz = vz/(awr + ONE)
|
||||
|
|
@ -267,13 +270,16 @@ contains
|
|||
|
||||
vel = sqrt(vx*vx + vy*vy + vz*vz)
|
||||
|
||||
! Select isotropic direcion -- this is only valid for s-wave
|
||||
! scattering
|
||||
phi = TWO*PI*rang()
|
||||
mu = TWO*rang() - ONE
|
||||
u = mu
|
||||
v = sqrt(ONE - mu**2) * cos(phi)
|
||||
w = sqrt(ONE - mu**2) * sin(phi)
|
||||
! Sample scattering angle
|
||||
mu = sample_angle(rxn, p % E)
|
||||
|
||||
! Determine direction cosines
|
||||
u = vx/vel
|
||||
v = vy/vel
|
||||
w = vz/vel
|
||||
|
||||
! Change direction cosines according to mu
|
||||
call rotate_angle(u, v, w, mu)
|
||||
|
||||
vx = u*vel
|
||||
vy = v*vel
|
||||
|
|
@ -523,6 +529,7 @@ contains
|
|||
real(8) :: mu ! cosine of scattering angle
|
||||
real(8) :: E ! outgoing energy in laboratory
|
||||
real(8) :: E_cm ! outgoing energy in center-of-mass
|
||||
real(8) :: u,v,w ! direction cosines
|
||||
character(250) :: msg ! error message
|
||||
|
||||
! copy energy of neutron
|
||||
|
|
@ -557,8 +564,14 @@ contains
|
|||
mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E)
|
||||
end if
|
||||
|
||||
! copy directional cosines
|
||||
u = p % uvw(1)
|
||||
v = p % uvw(2)
|
||||
w = p % uvw(3)
|
||||
|
||||
! change direction of particle
|
||||
call rotate_angle(p, mu)
|
||||
call rotate_angle(u, v, w, mu)
|
||||
p % uvw = (/ u, v, w /)
|
||||
|
||||
! change energy of particle
|
||||
p % E = E
|
||||
|
|
@ -583,6 +596,7 @@ contains
|
|||
real(8) :: mu ! cosine of scattering angle
|
||||
real(8) :: E ! outgoing energy in laboratory
|
||||
real(8) :: E_cm ! outgoing energy in center-of-mass
|
||||
real(8) :: u,v,w ! direction cosines
|
||||
character(250) :: msg ! error message
|
||||
|
||||
! copy energy of neutron
|
||||
|
|
@ -617,8 +631,14 @@ contains
|
|||
mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E)
|
||||
end if
|
||||
|
||||
! copy directional cosines
|
||||
u = p % uvw(1)
|
||||
v = p % uvw(2)
|
||||
w = p % uvw(3)
|
||||
|
||||
! change direction of particle
|
||||
call rotate_angle(p, mu)
|
||||
call rotate_angle(u, v, w, mu)
|
||||
p % uvw = (/ u, v, w /)
|
||||
|
||||
! change energy of particle
|
||||
p % E = E
|
||||
|
|
@ -632,10 +652,10 @@ contains
|
|||
end subroutine n_xn
|
||||
|
||||
!=====================================================================
|
||||
! N_GAMMA
|
||||
! N_ABSORPTION
|
||||
!=====================================================================
|
||||
|
||||
subroutine n_gamma(p)
|
||||
subroutine n_absorption(p)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
|
||||
|
|
@ -649,7 +669,7 @@ contains
|
|||
call message(msg, 10)
|
||||
end if
|
||||
|
||||
end subroutine n_gamma
|
||||
end subroutine n_absorption
|
||||
|
||||
!=====================================================================
|
||||
! SAMPLE_ANGLE samples the cosine of the angle between incident and
|
||||
|
|
@ -775,19 +795,22 @@ contains
|
|||
! done in MCNP and SERPENT.
|
||||
!=====================================================================
|
||||
|
||||
subroutine rotate_angle(p, mu)
|
||||
subroutine rotate_angle(u, v, w, mu)
|
||||
|
||||
type(Particle), pointer :: p
|
||||
real(8), intent(in) :: mu ! cosine of angle in lab
|
||||
! type(Particle), pointer :: p
|
||||
real(8), intent(inout) :: u
|
||||
real(8), intent(inout) :: v
|
||||
real(8), intent(inout) :: w
|
||||
real(8), intent(in) :: mu ! cosine of angle in lab
|
||||
|
||||
real(8) :: phi, sinphi, cosphi
|
||||
real(8) :: a,b
|
||||
real(8) :: u0, v0, w0
|
||||
|
||||
! Copy original directional cosines
|
||||
u0 = p % uvw(1)
|
||||
v0 = p % uvw(2)
|
||||
w0 = p % uvw(3)
|
||||
u0 = u
|
||||
v0 = v
|
||||
w0 = w
|
||||
|
||||
! Sample azimuthal angle in [0,2pi)
|
||||
phi = TWO * PI * rang()
|
||||
|
|
@ -801,14 +824,14 @@ contains
|
|||
! 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
|
||||
p % uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
|
||||
p % uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
|
||||
p % uvw(3) = mu*w0 - a*b*cosphi
|
||||
u = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
|
||||
v = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
|
||||
w = mu*w0 - a*b*cosphi
|
||||
else
|
||||
b = sqrt(ONE - v0*v0)
|
||||
p % uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
|
||||
p % uvw(2) = mu*v0 - a*b*cosphi
|
||||
p % uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
|
||||
u = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
|
||||
v = mu*v0 - a*b*cosphi
|
||||
w = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
|
||||
end if
|
||||
|
||||
end subroutine rotate_angle
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue