mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 21:25:36 -04:00
Remove Fortran surface reflection code
This commit is contained in:
Sterling Harper
parent
38aa76c9dc
commit
7af296902b
4 changed files with 26 additions and 234 deletions
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -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
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue