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

Bug fixes in geometry and secondary energy sampling.

Browse source
This commit is contained in:
Paul Romano 2011-03-24 21:46:20 +00:00
parent bf2b2df5b3
commit 948f2bab02
5 changed files with 153 additions and 67 deletions
Download patch file Download diff file Expand all files Collapse all files

9
ChangeLog
View file

@ -1,3 +1,12 @@
2011-03-24 Paul Romano <romano7@mit.edu>
* geometry.f90: Fixed significant bug in how quadratic surfaces
were handled. Changed x**2 to x*x in the interest of efficiency.
* global.f90: Added more interpolation types to prepare for
interpolation other than histogram or lin-lin.
* physics.f90: Fixed bug in E1 and EK energies for Law 4 (used NE
instead of NP).
2011-03-22 Paul Romano <romano7@mit.edu>
* benchmarks: Added PU-MET-FAST-002 and PU-MET-FAST-005 benchmark

128
src/geometry.f90
View file

@ -264,6 +264,7 @@ contains
real(8) :: a,b,c,k
real(8) :: quad
character(250) :: msg
logical :: on_surface
if (present(other_cell)) then
cell_p => cells(other_cell)
@ -301,7 +302,11 @@ contains
! check for coincident surface
index_surf = expression(i)
if (index_surf == current_surf) cycle
if (index_surf == current_surf) then
on_surface = .true.
else
on_surface = .false.
end if
! check for operators
index_surf = abs(index_surf)
@ -351,7 +356,7 @@ contains
end if
case (SURF_CYL_X)
a = ONE - u**2 ! v^2 + w^2
a = ONE - u*u ! v^2 + w^2
if (a == ZERO) then
d = INFINITY
else
@ -362,16 +367,32 @@ contains
y = y - y0
z = z - z0
k = y*v + z*w
c = y**2 + z**2 - r**2
quad = k**2 - a*c
c = y*y + z*z - r*r
quad = k*k - a*c
if (c < ZERO) then
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
elseif (on_surface) then
! particle is on the cylinder, thus one distance is
! positive/negative and the other is zero. The sign of
! k determines if we are facing in or out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cylinder, thus one distance
! must be negative and one must be positive. The
! positive distance will be the one with negative sign
! on sqrt(quad)
d = -(k - sqrt(quad))/a
d = (-k + sqrt(quad))/a
else
! particle is outside the cylinder, thus both
@ -379,14 +400,14 @@ contains
! positive, the smaller distance is the one with
! positive sign on sqrt(quad)
d = -(k + sqrt(quad))/a
d = (-k - sqrt(quad))/a
if (d < ZERO) d = INFINITY
end if
end if
case (SURF_CYL_Y)
a = ONE - v**2 ! u^2 + w^2
a = ONE - v*v ! u^2 + w^2
if (a == ZERO) then
d = INFINITY
else
@ -397,16 +418,32 @@ contains
x = x - x0
z = z - z0
k = x*u + z*w
c = x**2 + z**2 - r**2
quad = k**2 - a*c
c = x*x + z*z - r*r
quad = k*k - a*c
if (c < ZERO) then
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
elseif (on_surface) then
! particle is on the cylinder, thus one distance is
! positive/negative and the other is zero. The sign of
! k determines if we are facing in or out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cylinder, thus one distance
! must be negative and one must be positive. The
! positive distance will be the one with negative sign
! on sqrt(quad)
d = -(k - sqrt(quad))/a
d = (-k + sqrt(quad))/a
else
! particle is outside the cylinder, thus both
@ -414,14 +451,14 @@ contains
! positive, the smaller distance is the one with
! positive sign on sqrt(quad)
d = -(k + sqrt(quad))/a
d = (-k - sqrt(quad))/a
if (d < ZERO) d = INFINITY
end if
end if
case (SURF_CYL_Z)
a = ONE - w**2 ! u^2 + v^2
a = ONE - w*w ! u^2 + v^2
if (a == ZERO) then
d = INFINITY
else
@ -432,21 +469,32 @@ contains
x = x - x0
y = y - y0
k = x*u + y*v
c = x**2 + y**2 - r**2
quad = k**2 - a*c
c = x*x + y*y - r*r
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
elseif (on_surface) then
! particle is on the cylinder, thus one distance is
! positive/negative and the other is zero. The sign of
! k determines if we are facing in or out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cylinder, thus one distance
! must be negative and one must be positive. The
! positive distance will be the one with negative sign
! on sqrt(quad)
d = -(k - sqrt(quad))/a
d = (-k + sqrt(quad))/a
else
! particle is outside the cylinder, thus both
@ -454,8 +502,8 @@ contains
! positive, the smaller distance is the one with
! positive sign on sqrt(quad)
d = -(k + sqrt(quad))/a
if (d < ZERO) d = INFINITY
d = (-k - sqrt(quad))/a
if (d <= ZERO) d = INFINITY
end if
end if
@ -470,21 +518,32 @@ contains
y = y - y0
z = z - z0
k = x*u + y*v + z*w
c = x**2 + y**2 + z**2 - r**2
quad = k**2 - c
c = x*x + y*y + z*z - r*r
quad = k*k - c
if (quad < ZERO) then
! no intersection with sphere
d = INFINITY
elseif (on_surface) then
! particle is on the sphere, thus one distance is
! positive/negative and the other is zero. The sign of k
! determines if we are facing in or out
if (k >= ZERO) then
d = INFINITY
else
d = -k + sqrt(quad)
end if
elseif (c < ZERO) then
! particle is inside the sphere, thus one distance
! must be negative and one must be positive. The
! positive distance will be the one with negative sign
! on sqrt(quad)
d = -(k - sqrt(quad))
d = -k + sqrt(quad)
else
! particle is outside the sphere, thus both
@ -492,7 +551,7 @@ contains
! positive, the smaller distance is the one with
! positive sign on sqrt(quad)
d = -(k + sqrt(quad))
d = -k - sqrt(quad)
if (d < ZERO) d = INFINITY
end if
@ -562,26 +621,35 @@ contains
y0 = surf % coeffs(1)
z0 = surf % coeffs(2)
r = surf % coeffs(3)
func = (y-y0)**2 + (z-z0)**2 - r**2
y = y - y0
z = z - z0
func = y*y + z*z - r*r
case (SURF_CYL_Y)
x0 = surf % coeffs(1)
z0 = surf % coeffs(2)
r = surf % coeffs(3)
func = (x-x0)**2 + (z-z0)**2 - r**2
x = x - x0
z = z - z0
func = x*x + z*z - r*r
case (SURF_CYL_Z)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
r = surf % coeffs(3)
func = (x-x0)**2 + (y-y0)**2 - r**2
x = x - x0
y = y - y0
func = x*x + y*y - r*r
case (SURF_SPHERE)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
z0 = surf % coeffs(3)
r = surf % coeffs(4)
func = (x-x0)**2 + (y-y0)**2 + (z-z0)**2 - r**2
x = x - x0
y = y - y0
z = z - z0
func = x*x + y*y + z*z - r*r
case (SURF_BOX_X)
y0 = surf % coeffs(1)
@ -635,7 +703,7 @@ contains
return
case (SURF_GQ)
func = A*x**2 + B*y**2 + C*z**2 + D*x*y + E*y*z + F*x*z + G*x &
func = A*x*x + B*y*y + C*z*z + D*x*y + E*y*z + F*x*z + G*x &
& + H*y + I*z + J
end select

9
src/global.f90
View file

@ -150,9 +150,12 @@ module global
& ANGLE_TABULAR = 3 ! Tabular angular distribution
! Interpolation flag
integer, parameter :: &
& HISTOGRAM = 1, & ! Histogram interpolation on angle
& LINEARLINEAR = 2 ! Linear-linear interpolation on angle
integer, parameter :: &
& HISTOGRAM = 1, & ! y is constant in x
& LINEAR_LINEAR = 2, & ! y is linear in x
& LINEAR_LOG = 3, & ! y is linear in ln(x)
& LOG_LINEAR = 4, & ! ln(y) is linear in x
& LOG_LOG = 5 ! ln(y) is linear in ln(x)
! Particle type
integer, parameter :: &

8
src/input_sample
View file

@ -5,9 +5,9 @@ cell 100 37 40 -1
cell 101 37 41 1
cell 2 0 41 2 -3
surface 1 cylz 0 0 3
surface 2 cylz 0 0 5
surface 3 cylz 0 0 7
surface 1 cylz 0 0 7
surface 2 cylz 0 0 9
surface 3 cylz 0 0 11
material 40 -4.5 &
92235.03c 1.0
@ -21,4 +21,4 @@ source box -4 -4 -4 4 4 4
xs_library endfb7
xs_data /opt/serpent/xsdata/endfb7
criticality 15 5 10000
verbosity 8
verbosity 7

66
src/physics.f90
View file

@ -4,7 +4,7 @@ module physics
use geometry, only: find_cell, dist_to_boundary, cross_boundary
use types, only: Particle
use mcnp_random, only: rang
use output, only: error, message, print_particle
use output, only: error, warning, message, print_particle
use search, only: binary_search
use endf, only: reaction_name
@ -49,7 +49,7 @@ contains
if (verbosity >= 9) then
msg = "=== Particle " // trim(int_to_str(p % uid)) // " ==="
call message(msg, 10)
call message(msg, 9)
end if
if (verbosity >= 10) then
@ -90,6 +90,7 @@ contains
call cross_boundary(p)
else
! collision
p % surface = 0
call collision(p)
end if
@ -205,10 +206,10 @@ contains
& + f*(rxn%sigma(IE-rxn%IE+2)))
if (r1 < prob) exit
end do
if (verbosity >= 10) then
if (verbosity >= 10 .or. p % uid == 4593) then
msg = " " // trim(reaction_name(rxn%MT)) // " with nuclide " // &
& trim(table%name)
call message(msg, 10)
call message(msg, 8)
end if
! call appropriate subroutine
@ -227,6 +228,13 @@ contains
print *, 'warning: cant simulate reaction: ', rxn % MT
end select
! check for very low energy
if (p % E < 1.0e-100_8) then
p % alive = .false.
msg = "Killing neutron with extremely low energy"
call warning(msg)
end if
! find energy index, interpolation factor
call find_energy_index(p)
@ -764,7 +772,7 @@ contains
! Histogram interpolation
mu = mu0 + (xi - c_k)/p0
elseif (interp == LINEARLINEAR) then
elseif (interp == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about
! the formula given in the MCNP manual
p1 = rxn % adist % data(loc + NP + k+1)
@ -987,14 +995,16 @@ contains
l = i
end if
! interpolation on grid i for energy E1
loc = rxn%edist%data(2 + 2*NR + NE + i) + 2 ! start of EOUT(i)
E_i_1 = rxn%edist%data(loc + 1)
E_i_K = rxn%edist%data(loc + NE)
! interpolation for energy E1 and EK
loc = rxn%edist%data(2 + 2*NR + NE + i)
NP = rxn%edist%data(loc + 2)
E_i_1 = rxn%edist%data(loc + 2 + 1)
E_i_K = rxn%edist%data(loc + 2 + NP)
loc = rxn%edist%data(2 + 2*NR + NE + i + 1) + 2 ! start of EOUT(i+1)
E_i1_1 = rxn%edist%data(loc + 1)
E_i1_K = rxn%edist%data(loc + NE)
loc = rxn%edist%data(2 + 2*NR + NE + i + 1)
NP = rxn%edist%data(loc + 2)
E_i1_1 = rxn%edist%data(loc + 2 + 1)
E_i1_K = rxn%edist%data(loc + 2 + NP)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
@ -1036,7 +1046,7 @@ contains
! Histogram interpolation
E_out = E_l_k + (xi1 - c_k)/p_l_k
elseif (INTT == LINEARLINEAR) then
elseif (INTT == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about
! the formula given in the MCNP manual
E_l_k1 = rxn % edist % data(loc+k+1)
@ -1258,16 +1268,14 @@ contains
! determine endpoints on grid i
loc = rxn%edist%data(2+2*NR+NE + i) ! start of LDAT for i
NP = rxn%edist%data(loc + 2)
loc = loc + 2
E_i_1 = rxn%edist%data(loc + 1)
E_i_K = rxn%edist%data(loc + NP)
E_i_1 = rxn%edist%data(loc + 2 + 1)
E_i_K = rxn%edist%data(loc + 2 + NP)
! determine endpoints on grid i+1
loc = rxn%edist%data(2+2*NR+NE + i+1) ! start of LDAT for i+1
NP = rxn%edist%data(loc + 2)
loc = loc + 2
E_i1_1 = rxn%edist%data(loc + 1)
E_i1_K = rxn%edist%data(loc + NP)
NP = rxn%edist%data(loc + 2)
E_i1_1 = rxn%edist%data(loc + 2 + 1)
E_i1_K = rxn%edist%data(loc + 2 + NP)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
@ -1313,7 +1321,7 @@ contains
KM_R = rxn % edist % data(loc + 3*NP + k)
KM_A = rxn % edist % data(loc + 4*NP + k)
elseif (INTT == LINEARLINEAR) then
elseif (INTT == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about
! the formula given in the MCNP manual
E_l_k1 = rxn % edist % data(loc+k+1)
@ -1403,16 +1411,14 @@ contains
! determine endpoints on grid i
loc = rxn%edist%data(2+2*NR+NE + i) ! start of LDAT for i
NP = rxn%edist%data(loc + 2)
loc = loc + 2
E_i_1 = rxn%edist%data(loc + 1)
E_i_K = rxn%edist%data(loc + NP)
E_i_1 = rxn%edist%data(loc + 2 + 1)
E_i_K = rxn%edist%data(loc + 2 + NP)
! determine endpoints on grid i+1
loc = rxn%edist%data(2+2*NR+NE + i+1) ! start of LDAT for i+1
NP = rxn%edist%data(loc + 2)
loc = loc + 2
E_i1_1 = rxn%edist%data(loc + 1)
E_i1_K = rxn%edist%data(loc + NP)
NP = rxn%edist%data(loc + 2)
E_i1_1 = rxn%edist%data(loc + 2 + 1)
E_i1_K = rxn%edist%data(loc + 2 + NP)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
@ -1454,7 +1460,7 @@ contains
! Histogram interpolation
E_out = E_l_k + (xi1 - c_k)/p_l_k
elseif (INTT == LINEARLINEAR) then
elseif (INTT == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about
! the formula given in the MCNP manual
E_l_k1 = rxn % edist % data(loc+k+1)
@ -1509,7 +1515,7 @@ contains
! Histogram interpolation
mu_out = mu_k + (xi3 - c_k)/p_k
elseif (JJ == LINEARLINEAR) then
elseif (JJ == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about
! the formula given in the MCNP manual
p_k1 = rxn % edist % data(loc + NP + k+1)