From 7af296902bbb9ae7b0409eafb5060feacfb1d2ec Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 19 Jan 2018 01:27:04 -0500 Subject: [PATCH] Remove Fortran surface reflection code --- src/geometry.F90 | 15 ++- src/surface_header.C | 14 ++- src/surface_header.F90 | 226 +---------------------------------------- src/tracking.F90 | 5 +- 4 files changed, 26 insertions(+), 234 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index a0a19f1bf..c19a6c8dd 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -25,6 +25,15 @@ module geometry logical(C_BOOL) :: sense; end function surface_sense_c + pure subroutine surface_reflect_c(surf_ind, xyz, uvw) & + bind(C, name='surface_reflect') + use ISO_C_BINDING + implicit none + integer(C_INT), intent(in), value :: surf_ind; + real(C_DOUBLE), intent(in) :: xyz(3); + real(C_DOUBLE), intent(inout) :: uvw(3); + end subroutine surface_reflect_c + pure function surface_distance_c(surf_ind, xyz, uvw, coincident) & bind(C, name='surface_distance') result(d) use ISO_C_BINDING @@ -525,6 +534,7 @@ contains real(8) :: d_surf ! distance to surface real(8) :: x0,y0,z0 ! coefficients for surface real(8) :: xyz_cross(3) ! coordinates at projected surface crossing + real(8) :: surf_uvw(3) ! surface normal direction logical :: coincident ! is particle on surface? type(Cell), pointer :: c class(Surface), pointer :: surf @@ -782,9 +792,8 @@ contains ! traveling into if the surface is crossed if (.not. c % simple) then xyz_cross(:) = p % coord(j) % xyz + d_surf*p % coord(j) % uvw - surf => surfaces(abs(level_surf_cross)) % obj - if (dot_product(p % coord(j) % uvw, & - surf % normal(xyz_cross)) > ZERO) then + call surface_normal_c(abs(level_surf_cross)-1, xyz_cross, surf_uvw) + if (dot_product(p % coord(j) % uvw, surf_uvw) > ZERO) then surface_crossed = abs(level_surf_cross) else surface_crossed = -abs(level_surf_cross) diff --git a/src/surface_header.C b/src/surface_header.C index 23fe5adb3..45fb62399 100644 --- a/src/surface_header.C +++ b/src/surface_header.C @@ -233,8 +233,8 @@ class SurfaceXPlane : public PeriodicSurface public: SurfaceXPlane(pugi::xml_node surf_node); double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - const bool coincident) const; + double distance(const double xyz[3], const double uvw[3], bool coincident) + const; void normal(const double xyz[3], double uvw[3]) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; @@ -416,8 +416,8 @@ class SurfacePlane : public PeriodicSurface public: SurfacePlane(pugi::xml_node surf_node); double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - const bool coincident) const; + double distance(const double xyz[3], const double uvw[3], bool coincident) + const; void normal(const double xyz[3], double uvw[3]) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; @@ -1138,6 +1138,12 @@ surface_sense(int surf_ind, double xyz[3], double uvw[3]) return surfaces_c[surf_ind]->sense(xyz, uvw); } +extern "C" void +surface_reflect(int surf_ind, double xyz[3], double uvw[3]) +{ + surfaces_c[surf_ind]->reflect(xyz, uvw); +} + extern "C" double surface_distance(int surf_ind, double xyz[3], double uvw[3], bool coincident) { diff --git a/src/surface_header.F90 b/src/surface_header.F90 index d9a00170b..d63290387 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -2,7 +2,7 @@ module surface_header use, intrinsic :: ISO_C_BINDING - use constants, only: NONE, ONE, TWO, ZERO, HALF, INFINITY, FP_COINCIDENT + use constants, only: NONE use dict_header, only: DictIntInt implicit none @@ -20,20 +20,8 @@ module surface_header integer :: bc ! Boundary condition integer :: i_periodic = NONE ! Index of corresponding periodic surface character(len=104) :: name = "" ! User-defined name - contains - procedure :: reflect - procedure(surface_normal_), deferred :: normal end type Surface - abstract interface - pure function surface_normal_(this, xyz) result(uvw) - import Surface - class(Surface), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - end function surface_normal_ - end interface - !=============================================================================== ! SURFACECONTAINER allows us to store an array of different types of surfaces !=============================================================================== @@ -50,22 +38,16 @@ module surface_header type, extends(Surface) :: SurfaceXPlane ! x = x0 real(8) :: x0 - contains - procedure :: normal => x_plane_normal end type SurfaceXPlane type, extends(Surface) :: SurfaceYPlane ! y = y0 real(8) :: y0 - contains - procedure :: normal => y_plane_normal end type SurfaceYPlane type, extends(Surface) :: SurfaceZPlane ! z = z0 real(8) :: z0 - contains - procedure :: normal => z_plane_normal end type SurfaceZPlane type, extends(Surface) :: SurfacePlane @@ -74,8 +56,6 @@ module surface_header real(8) :: B real(8) :: C real(8) :: D - contains - procedure :: normal => plane_normal end type SurfacePlane type, extends(Surface) :: SurfaceXCylinder @@ -83,8 +63,6 @@ module surface_header real(8) :: y0 real(8) :: z0 real(8) :: r - contains - procedure :: normal => x_cylinder_normal end type SurfaceXCylinder type, extends(Surface) :: SurfaceYCylinder @@ -92,8 +70,6 @@ module surface_header real(8) :: x0 real(8) :: z0 real(8) :: r - contains - procedure :: normal => y_cylinder_normal end type SurfaceYCylinder type, extends(Surface) :: SurfaceZCylinder @@ -101,8 +77,6 @@ module surface_header real(8) :: x0 real(8) :: y0 real(8) :: r - contains - procedure :: normal => z_cylinder_normal end type SurfaceZCylinder type, extends(Surface) :: SurfaceSphere @@ -111,8 +85,6 @@ module surface_header real(8) :: y0 real(8) :: z0 real(8) :: r - contains - procedure :: normal => sphere_normal end type SurfaceSphere type, extends(Surface) :: SurfaceXCone @@ -121,8 +93,6 @@ module surface_header real(8) :: y0 real(8) :: z0 real(8) :: r2 - contains - procedure :: normal => x_cone_normal end type SurfaceXCone type, extends(Surface) :: SurfaceYCone @@ -131,8 +101,6 @@ module surface_header real(8) :: y0 real(8) :: z0 real(8) :: r2 - contains - procedure :: normal => y_cone_normal end type SurfaceYCone type, extends(Surface) :: SurfaceZCone @@ -141,15 +109,11 @@ module surface_header real(8) :: y0 real(8) :: z0 real(8) :: r2 - contains - procedure :: normal => z_cone_normal end type SurfaceZCone type, extends(Surface) :: SurfaceQuadric ! Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 real(8) :: A, B, C, D, E, F, G, H, J, K - contains - procedure :: normal => quadric_normal end type SurfaceQuadric integer(C_INT32_T), bind(C) :: n_surfaces ! # of surfaces @@ -161,194 +125,6 @@ module surface_header contains -!=============================================================================== -! REFLECT determines the direction a particle will travel if it is specularly -! reflected from the surface at a given position and direction. The position is -! needed because the reflection is performed using the surface normal, which -! depends on the position for second-order surfaces. -!=============================================================================== - - pure subroutine reflect(this, xyz, uvw) - class(Surface), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(inout) :: uvw(3) - - real(8) :: projection - real(8) :: magnitude - real(8) :: n(3) - - ! Construct normal vector - n(:) = this%normal(xyz) - - ! Determine projection of direction onto normal and squared magnitude of - ! normal - projection = n(1)*uvw(1) + n(2)*uvw(2) + n(3)*uvw(3) - magnitude = n(1)*n(1) + n(2)*n(2) + n(3)*n(3) - - ! Reflect direction according to normal - uvw(:) = uvw - TWO*projection/magnitude * n - end subroutine reflect - -!=============================================================================== -! SurfaceXPlane Implementation -!=============================================================================== - - pure function x_plane_normal(this, xyz) result(uvw) - class(SurfaceXPlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(:) = [ONE, ZERO, ZERO] - end function x_plane_normal - -!=============================================================================== -! SurfaceYPlane Implementation -!=============================================================================== - - pure function y_plane_normal(this, xyz) result(uvw) - class(SurfaceYPlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(:) = [ZERO, ONE, ZERO] - end function y_plane_normal - -!=============================================================================== -! SurfaceZPlane Implementation -!=============================================================================== - - pure function z_plane_normal(this, xyz) result(uvw) - class(SurfaceZPlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(:) = [ZERO, ZERO, ONE] - end function z_plane_normal - -!=============================================================================== -! SurfacePlane Implementation -!=============================================================================== - - pure function plane_normal(this, xyz) result(uvw) - class(SurfacePlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(:) = [this%A, this%B, this%C] - end function plane_normal - -!=============================================================================== -! SurfaceXCylinder Implementation -!=============================================================================== - - pure function x_cylinder_normal(this, xyz) result(uvw) - class(SurfaceXCylinder), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(1) = ZERO - uvw(2) = TWO*(xyz(2) - this%y0) - uvw(3) = TWO*(xyz(3) - this%z0) - end function x_cylinder_normal - -!=============================================================================== -! SurfaceYCylinder Implementation -!=============================================================================== - - pure function y_cylinder_normal(this, xyz) result(uvw) - class(SurfaceYCylinder), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(1) = TWO*(xyz(1) - this%x0) - uvw(2) = ZERO - uvw(3) = TWO*(xyz(3) - this%z0) - end function y_cylinder_normal - -!=============================================================================== -! SurfaceZCylinder Implementation -!=============================================================================== - - pure function z_cylinder_normal(this, xyz) result(uvw) - class(SurfaceZCylinder), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(1) = TWO*(xyz(1) - this%x0) - uvw(2) = TWO*(xyz(2) - this%y0) - uvw(3) = ZERO - end function z_cylinder_normal - -!=============================================================================== -! SurfaceSphere Implementation -!=============================================================================== - - pure function sphere_normal(this, xyz) result(uvw) - class(SurfaceSphere), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(:) = TWO*(xyz - [this%x0, this%y0, this%z0]) - end function sphere_normal - -!=============================================================================== -! SurfaceXCone Implementation -!=============================================================================== - - pure function x_cone_normal(this, xyz) result(uvw) - class(SurfaceXCone), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(1) = -TWO*this%r2*(xyz(1) - this%x0) - uvw(2) = TWO*(xyz(2) - this%y0) - uvw(3) = TWO*(xyz(3) - this%z0) - end function x_cone_normal - -!=============================================================================== -! SurfaceYCone Implementation -!=============================================================================== - - pure function y_cone_normal(this, xyz) result(uvw) - class(SurfaceYCone), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(1) = TWO*(xyz(1) - this%x0) - uvw(2) = -TWO*this%r2*(xyz(2) - this%y0) - uvw(3) = TWO*(xyz(3) - this%z0) - end function y_cone_normal - -!=============================================================================== -! SurfaceZCone Implementation -!=============================================================================== - - pure function z_cone_normal(this, xyz) result(uvw) - class(SurfaceZCone), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - uvw(1) = TWO*(xyz(1) - this%x0) - uvw(2) = TWO*(xyz(2) - this%y0) - uvw(3) = -TWO*this%r2*(xyz(3) - this%z0) - end function z_cone_normal - -!=============================================================================== -! SurfaceQuadric Implementation -!=============================================================================== - - pure function quadric_normal(this, xyz) result(uvw) - class(SurfaceQuadric), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: uvw(3) - - associate (x => xyz(1), y => xyz(2), z => xyz(3)) - uvw(1) = TWO*this%A*x + this%D*y + this%F*z + this%G - uvw(2) = TWO*this%B*y + this%D*x + this%E*z + this%H - uvw(3) = TWO*this%C*z + this%E*y + this%F*x + this%J - end associate - end function quadric_normal - !=============================================================================== ! FREE_MEMORY_SURFACES deallocates global arrays defined in this module !=============================================================================== diff --git a/src/tracking.F90 b/src/tracking.F90 index 4f0f110bc..d616212ec 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -5,7 +5,8 @@ module tracking use error, only: fatal_error, warning, write_message use geometry_header, only: cells use geometry, only: find_cell, distance_to_boundary, cross_lattice, & - check_cell_overlap, surface_periodic_c + check_cell_overlap, surface_reflect_c, & + surface_periodic_c use message_passing use mgxs_header use nuclide_header @@ -354,7 +355,7 @@ contains end if ! Reflect particle off surface - call surf%reflect(p%coord(1)%xyz, p%coord(1)%uvw) + call surface_reflect_c(i_surface-1, p%coord(1)%xyz, p%coord(1)%uvw) ! Make sure new particle direction is normalized u = p%coord(1)%uvw(1)