diff --git a/src/Makefile b/src/Makefile index 1649c82a7..41d0ab1e0 100644 --- a/src/Makefile +++ b/src/Makefile @@ -35,7 +35,7 @@ global.o: types.o string.o: global.o output.o output.o: global.o source.o: global.o mcnp_random.o -physics.o: types.o global.o +physics.o: types.o global.o mcnp_random.o geometry.o output.o geometry.o: types.o global.o output.o string.o fileio.o: types.o global.o string.o output.o main.o: global.o fileio.o output.o geometry.o mcnp_random.o \ diff --git a/src/geometry.f90 b/src/geometry.f90 index 5f41246e2..34dee79cb 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -188,12 +188,12 @@ contains integer, intent(out) :: surf integer, optional, intent(in) :: other_cell - type(Cell), pointer :: cl => null() - type(Surface), pointer :: surf => null() + type(Cell), pointer :: cell_p => null() + type(Surface), pointer :: surf_p => null() integer :: i integer :: n_boundaries integer, allocatable :: expression(:) - integer :: surf_num + integer :: index_surf real(8) :: x,y,z,u,v,w real(8) :: d real(8) :: x0,y0,z0,r @@ -203,9 +203,9 @@ contains character(250) :: msg if ( present(other_cell) ) then - cl => cells(other_cell) + cell_p => cells(other_cell) else - cl => cells(neut%cell) + cell_p => cells(neut%cell) end if x = neut%xyz(1) @@ -216,18 +216,20 @@ contains z = neut%uvw(3) dist = INFINITY - n_boundaries = size(cl%boundary_list) + n_boundaries = size(cell_p%boundary_list) + allocate( expression(n_boundaries) ) + expression = cell_p%boundary_list do i = 1, n_boundaries - surf_num = abs(expression(i)) - if ( surf_num >= OP_DIFFERENCE ) cycle + index_surf = abs(expression(i)) + if ( index_surf >= OP_DIFFERENCE ) cycle - surf => surfaces(surf_num) - select case ( surf%type ) + surf_p => surfaces(index_surf) + select case ( surf_p%type ) case ( SURF_PX ) if ( u == 0.0 ) then d = INFINITY else - x0 = surf%coeffs(1) + x0 = surf_p%coeffs(1) d = (x0 - x)/u if ( d < 0 ) d = INFINITY end if @@ -236,7 +238,7 @@ contains if ( v == 0.0 ) then d = INFINITY else - y0 = surf%coeffs(1) + y0 = surf_p%coeffs(1) d = (y0 - y)/v if ( d < 0 ) d = INFINITY end if @@ -245,16 +247,16 @@ contains if ( w == 0.0 ) then d = INFINITY else - z0 = surf%coeffs(1) + z0 = surf_p%coeffs(1) d = (z0 - z)/w if ( d < 0.0 ) d = INFINITY end if case ( SURF_PLANE ) - A = surf%coeffs(1) - B = surf%coeffs(2) - C = surf%coeffs(3) - D = surf%coeffs(4) + A = surf_p%coeffs(1) + B = surf_p%coeffs(2) + C = surf_p%coeffs(3) + D = surf_p%coeffs(4) tmp = A*u + B*v + C*w if ( tmp == 0.0 ) then @@ -269,9 +271,9 @@ contains if ( a == 0.0 ) then d = INFINITY else - y0 = surf%coeffs(1) - z0 = surf%coeffs(2) - r = surf%coeffs(3) + y0 = surf_p%coeffs(1) + z0 = surf_p%coeffs(2) + r = surf_p%coeffs(3) y = y - y0 z = z - z0 @@ -304,9 +306,9 @@ contains if ( a == 0.0 ) then d = INFINITY else - x0 = surf%coeffs(1) - z0 = surf%coeffs(2) - r = surf%coeffs(3) + x0 = surf_p%coeffs(1) + z0 = surf_p%coeffs(2) + r = surf_p%coeffs(3) x = x - x0 z = z - z0 @@ -339,9 +341,9 @@ contains if ( a == 0.0 ) then d = INFINITY else - x0 = surf%coeffs(1) - y0 = surf%coeffs(2) - r = surf%coeffs(3) + x0 = surf_p%coeffs(1) + y0 = surf_p%coeffs(2) + r = surf_p%coeffs(3) x = x - x0 y = y - y0 @@ -375,10 +377,10 @@ contains end if case ( SURF_SPHERE ) - x0 = surf%coeffs(1) - y0 = surf%coeffs(2) - z0 = surf%coeffs(3) - r = surf%coeffs(4) + x0 = surf_p%coeffs(1) + y0 = surf_p%coeffs(2) + z0 = surf_p%coeffs(3) + r = surf_p%coeffs(4) x = x - x0 y = y - y0 @@ -425,6 +427,9 @@ contains end do + ! deallocate expression + deallocate( expression ) + end subroutine dist_to_boundary !===================================================================== diff --git a/src/global.f90 b/src/global.f90 index b17ebfa8d..5d3a720ac 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -30,10 +30,10 @@ module global & INFINITY = huge(0.0_8) ! positive infinity ! Boundary conditions - integer, parameter :: & - & TRANSMIT = 0, & ! Transmission boundary condition (default) - & VACUUM = 1, & ! Vacuum boundary condition - & REFLECT = 2 ! Reflecting boundary condition + integer, parameter :: & + & BC_TRANSMIT = 0, & ! Transmission boundary condition (default) + & BC_VACUUM = 1, & ! Vacuum boundary condition + & BC_REFLECT = 2 ! Reflecting boundary condition ! Logical operators for cell definitions integer, parameter :: & diff --git a/src/input_sample b/src/input_sample index 333ddb426..fc2720427 100644 --- a/src/input_sample +++ b/src/input_sample @@ -1,9 +1,9 @@ # test input file cell 1 1 -2 -cell 2 1 2 -22 +cell 2 1 2 -1 -surface 1 sph 0 0 0 1 +surface 1 sph 0 0 0 5 surface 2 sph 0 0.0 0.0 3 material 1 U235 1.00 diff --git a/src/physics.f90 b/src/physics.f90 index 6130786df..c9b044aef 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -1,8 +1,10 @@ module physics use global - use geometry, only: find_cell, dist_to_boundary - use types, only: Neutron + use geometry, only: find_cell, dist_to_boundary, cross_boundary + use types, only: Neutron + use mcnp_random, only: rang + use output, only: error implicit none @@ -33,12 +35,11 @@ contains ! Determine distance neutron moves call dist_to_boundary( neut, d_to_boundary, surf ) - print *, d_to_boundary, surf d_to_collision = -log(rang()) / 1.0 distance = min( d_to_boundary, d_to_collision ) ! Advance neutron - neut%xyz = p%xyz + distance*p%uvw + neut%xyz = neut%xyz + distance*neut%uvw ! Add pathlength tallies @@ -48,6 +49,8 @@ contains call cross_boundary( neut ) else ! collision + msg = "Collision not implemented yet!" + call error( msg ) end if end do