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 distance to surface for cones, again. Closes #143.

Browse source
This commit is contained in:
Paul Romano 2013-02-01 22:06:50 -05:00
parent dc01e73e15
commit b0fe887aeb
2 changed files with 64 additions and 55 deletions
Download patch file Download diff file Expand all files Collapse all files

118
src/geometry.F90
View file

@ -715,6 +715,10 @@ contains
! get pointer to surface
surf => surfaces(index_surf)
! TODO: Can probably combines a lot of the cases to reduce repetition
! since the algorithm is the same for (x-plane, y-plane, z-plane),
! (x-cylinder, y-cylinder, z-cylinder), etc.
select case (surf % type)
case (SURF_PX)
if (on_surface .or. u == ZERO) then
@ -964,39 +968,41 @@ contains
c = y*y + z*z - r*x*x
quad = k*k - a*c
if (quad < ZERO .or. a == ZERO) then
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
elseif (on_surface) then
! particle is on the cone, thus one distance is
! positive/negative and the other is zero. The sign of k
! determines if we are facing in or out
! particle is on the cone, thus one distance is positive/negative
! and the other is zero. The sign of k determines which distance is
! zero and which is not.
if (k >= ZERO) then
d = INFINITY
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cone, 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
else
! particle is outside the cone, thus both distances are either
! positive or negative. If positive, the smaller distance is the
! one with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d <= ZERO) d = INFINITY
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
case (SURF_CONE_Y)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
@ -1011,39 +1017,41 @@ contains
c = x*x + z*z - r*y*y
quad = k*k - a*c
if (quad < ZERO .or. a == ZERO) then
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
elseif (on_surface) then
! particle is on the cone, thus one distance is
! positive/negative and the other is zero. The sign of k
! determines if we are facing in or out
! particle is on the cone, thus one distance is positive/negative
! and the other is zero. The sign of k determines which distance is
! zero and which is not.
if (k >= ZERO) then
d = INFINITY
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cone, 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
else
! particle is outside the cone, thus both distances are either
! positive or negative. If positive, the smaller distance is the
! one with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d <= ZERO) d = INFINITY
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
case (SURF_CONE_Z)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
@ -1058,39 +1066,41 @@ contains
c = x*x + y*y - r*z*z
quad = k*k - a*c
if (quad < ZERO .or. a == ZERO) then
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
elseif (on_surface) then
! particle is on the cone, thus one distance is
! positive/negative and the other is zero. The sign of k
! determines if we are facing in or out
! particle is on the cone, thus one distance is positive/negative
! and the other is zero. The sign of k determines which distance is
! zero and which is not.
if (k >= ZERO) then
d = INFINITY
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cone, 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
else
! particle is outside the cone, thus both distances are either
! positive or negative. If positive, the smaller distance is the
! one with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d <= ZERO) d = INFINITY
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
case (SURF_GQ)
message = "Surface distance not yet implement for general quadratic."
call fatal_error()

1
tests/test_reflective_cone/geometry.xml
View file

@ -1,7 +1,6 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="z-cone" coeffs="0 0 0 5" boundary="reflective"/>
<cell id="1" material="1" surfaces="-1" />