From 948f2bab02b9631af61de93d94ce1c2cbc4494b1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 24 Mar 2011 21:46:20 +0000 Subject: [PATCH] Bug fixes in geometry and secondary energy sampling. --- ChangeLog | 9 ++++ src/geometry.f90 | 128 ++++++++++++++++++++++++++++++++++++----------- src/global.f90 | 9 ++-- src/input_sample | 8 +-- src/physics.f90 | 66 +++++++++++++----------- 5 files changed, 153 insertions(+), 67 deletions(-) diff --git a/ChangeLog b/ChangeLog index bc4be65a0..6c2b1d731 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,12 @@ +2011-03-24 Paul Romano + + * 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 * benchmarks: Added PU-MET-FAST-002 and PU-MET-FAST-005 benchmark diff --git a/src/geometry.f90 b/src/geometry.f90 index 5873b4b55..761464e53 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -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 diff --git a/src/global.f90 b/src/global.f90 index a1c3dec0b..f68feb91c 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -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 :: & diff --git a/src/input_sample b/src/input_sample index 53c7b920f..cc4ade038 100644 --- a/src/input_sample +++ b/src/input_sample @@ -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 diff --git a/src/physics.f90 b/src/physics.f90 index eaf413e72..3c4c810ff 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -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)