From 38aa76c9dc9cdab0c10f35c2355d20f9842b65bc Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 19 Jan 2018 00:54:39 -0500 Subject: [PATCH] Implement C++ periodic BCs --- src/geometry.F90 | 11 +++ src/surface_header.C | 136 +++++++++++++++++++++++++++++++++-- src/surface_header.F90 | 157 +---------------------------------------- src/tracking.F90 | 63 +---------------- 4 files changed, 145 insertions(+), 222 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 873582dec9..a0a19f1bfc 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -44,6 +44,17 @@ module geometry real(C_DOUBLE), intent(in) :: xyz(3); real(C_DOUBLE), intent(out) :: uvw(3); end subroutine surface_normal_c + + function surface_periodic_c(surf_ind1, surf_ind2, xyz, uvw) & + bind(C, name="surface_periodic") result(rotational) + use ISO_C_BINDING + implicit none + integer(C_INT), intent(in), value :: surf_ind1; + integer(C_INT), intent(in), value :: surf_ind2; + real(C_DOUBLE), intent(inout) :: xyz(3); + real(C_DOUBLE), intent(inout) :: uvw(3); + logical(C_BOOL) :: rotational + end function surface_periodic_c end interface contains diff --git a/src/surface_header.C b/src/surface_header.C index cf200e5e46..23fe5adb35 100644 --- a/src/surface_header.C +++ b/src/surface_header.C @@ -5,6 +5,7 @@ #include "pugixml/pugixml.hpp" // DEBUGGING +#include #include //============================================================================== @@ -89,14 +90,14 @@ public: neighbor_neg[]; //!< List of cells on negative side int bc; //!< Boundary condition //TODO: switch that zero to a NONE constant. - int i_periodic = 0; //!< Index of corresponding periodic surface + //int i_periodic = 0; //!< Index of corresponding periodic surface char name[104]; //!< User-defined name //! Determine which side of a surface a point lies on. //! @param xyz[3] The 3D Cartesian coordinate of a point. //! @param uvw[3] A direction used to "break ties" and pick a sense when the //! point is very close to the surface. - //! @return True if the point is on the "positive" side of the surface and + //! @return true if the point is on the "positive" side of the surface and //! false otherwise. bool sense(const double xyz[3], const double uvw[3]) const; @@ -163,6 +164,29 @@ Surface::reflect(const double xyz[3], double uvw[3]) const uvw[2] -= 2.0 * projection / magnitude * norm[2]; } +//============================================================================== +//! A `Surface` that supports periodic boundary conditions. +//! +//! Translational periodicity is supported for the `XPlane`, `YPlane`, `ZPlane`, +//! and `Plane` types. Rotational periodicity is supported for +//! `XPlane`-`YPlane` pairs. +//============================================================================== + +class PeriodicSurface : public Surface +{ +public: + //! Translate a particle onto this surface from a periodic partner surface. + //! @param other A pointer to the partner surface in this periodic BC. + //! @param xyz[3] A point on the partner surface that will be translated onto + //! this surface. + //! @param uvw[3] A direction that will be rotated for systems with rotational + //! periodicity. + //! @return true if this surface and its partner make a rotationally-periodic + //! boundary condition. + virtual bool periodic_translate(PeriodicSurface *other, double xyz[3], + double uvw[3]) const = 0; +}; + //============================================================================== // Generic functions for x-, y-, and z-, planes. //============================================================================== @@ -203,7 +227,7 @@ axis_aligned_plane_normal(const double xyz[3], double uvw[3]) //! The plane is described by the equation \f$x - x_0 = 0\f$ //============================================================================== -class SurfaceXPlane : public Surface +class SurfaceXPlane : public PeriodicSurface { double x0; public: @@ -212,6 +236,8 @@ public: double distance(const double xyz[3], const double uvw[3], const 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; }; SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node) @@ -235,6 +261,32 @@ inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const axis_aligned_plane_normal<0, 1, 2>(xyz, uvw); } +bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], + double uvw[3]) const +{ + double other_norm[3]; + other->normal(xyz, other_norm); + if (other_norm[0] == 1 and other_norm[1] == 0 and other_norm[2] == 0) { + xyz[0] = x0; + return false; + } else { + // Assume the partner is an YPlane (the only supported partner). Use the + // evaluate function to find y0, then adjust xyz and uvw for rotational + // symmetry. + double xyz_test[3] {0, 0, 0}; + double y0 = -other->evaluate(xyz_test); + xyz[1] = xyz[0] - x0 + y0; + xyz[0] = x0; + xyz[2] = xyz[2]; + + double u = uvw[0]; + uvw[0] = -uvw[1]; + uvw[1] = u; + + return true; + } +}; + //============================================================================== // SurfaceYPlane //! A plane perpendicular to the y-axis. @@ -242,7 +294,7 @@ inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const //! The plane is described by the equation \f$y - y_0 = 0\f$ //============================================================================== -class SurfaceYPlane : public Surface +class SurfaceYPlane : public PeriodicSurface { double y0; public: @@ -251,6 +303,8 @@ public: 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; }; SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) @@ -274,6 +328,32 @@ inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const axis_aligned_plane_normal<1, 0, 2>(xyz, uvw); } +bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], + double uvw[3]) const +{ + double other_norm[3]; + other->normal(xyz, other_norm); + if (other_norm[0] == 0 and other_norm[1] == 1 and other_norm[2] == 0) { + // The periodic partner is also aligned along y. Just change the y coord. + xyz[1] = y0; + return false; + } else { + // Assume the partner is an XPlane (the only supported partner). Use the + // evaluate function to find x0, then adjust xyz and uvw for rotational + // symmetry. + double xyz_test[3] {0, 0, 0}; + double x0 = -other->evaluate(xyz_test); + xyz[0] = xyz[1] - y0 + x0; + xyz[1] = y0; + + double u = uvw[0]; + uvw[0] = uvw[1]; + uvw[1] = -u; + + return true; + } +}; + //============================================================================== // SurfaceZPlane //! A plane perpendicular to the z-axis. @@ -281,7 +361,7 @@ inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const //! The plane is described by the equation \f$z - z_0 = 0\f$ //============================================================================== -class SurfaceZPlane : public Surface +class SurfaceZPlane : public PeriodicSurface { double z0; public: @@ -290,6 +370,8 @@ public: 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; }; SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) @@ -313,6 +395,14 @@ inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const axis_aligned_plane_normal<2, 0, 1>(xyz, uvw); } +bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3], + double uvw[3]) const +{ + // Assume the other plane is aligned along z. Just change the z coord. + xyz[2] = z0; + return false; +}; + //============================================================================== // SurfacePlane //! A general plane. @@ -320,7 +410,7 @@ inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const //! The plane is described by the equation \f$A x + B y + C z - D = 0\f$ //============================================================================== -class SurfacePlane : public Surface +class SurfacePlane : public PeriodicSurface { double A, B, C, D; public: @@ -329,6 +419,8 @@ public: double distance(const double xyz[3], const double uvw[3], const 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; }; SurfacePlane::SurfacePlane(pugi::xml_node surf_node) @@ -365,6 +457,22 @@ SurfacePlane::normal(const double xyz[3], double uvw[3]) const uvw[2] = C; } +bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], + double uvw[3]) const +{ + // This function assumes the other plane shares this plane's normal direction. + + // Determine the distance to intersection. + double d = evaluate(xyz) / (A*A + B*B + C*C); + + // Move the particle that distance along the normal vector. + xyz[0] -= d * A; + xyz[1] -= d * B; + xyz[2] -= d * C; + + return false; +} + //============================================================================== // Generic functions for x-, y-, and z-, cylinders //============================================================================== @@ -1041,3 +1149,19 @@ surface_normal(int surf_ind, double xyz[3], double uvw[3]) { return surfaces_c[surf_ind]->normal(xyz, uvw); } + +extern "C" bool +surface_periodic(int surf_ind1, int surf_ind2, double xyz[3], double uvw[3]) +{ + // Hopefully this function has only been called on a pair of surfaces that + // support periodic BCs (checking should have been done when reading the + // geometry XML). Downcast the surfaces to the PeriodicSurface type so we + // can call the periodic_translate method. + Surface *surf1_gen = surfaces_c[surf_ind1]; + PeriodicSurface *surf1 = dynamic_cast(surf1_gen); + Surface *surf2_gen = surfaces_c[surf_ind2]; + PeriodicSurface *surf2 = dynamic_cast(surf2_gen); + + // Call the type-bound methods. + return surf2->periodic_translate(surf1, xyz, uvw); +} diff --git a/src/surface_header.F90 b/src/surface_header.F90 index b0da37c1e0..d9a00170ba 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -22,18 +22,10 @@ module surface_header character(len=104) :: name = "" ! User-defined name contains procedure :: reflect - procedure(surface_evaluate_), deferred :: evaluate procedure(surface_normal_), deferred :: normal end type Surface abstract interface - pure function surface_evaluate_(this, xyz) result(f) - import Surface - class(Surface), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - end function surface_evaluate_ - pure function surface_normal_(this, xyz) result(uvw) import Surface class(Surface), intent(in) :: this @@ -52,15 +44,13 @@ module surface_header !=============================================================================== ! All the derived types below are extensions of the abstract Surface type. They -! inherent the reflect() type-bound procedure and must implement evaluate(), -! and normal() +! inherent the reflect() type-bound procedure and must implement normal() !=============================================================================== type, extends(Surface) :: SurfaceXPlane ! x = x0 real(8) :: x0 contains - procedure :: evaluate => x_plane_evaluate procedure :: normal => x_plane_normal end type SurfaceXPlane @@ -68,7 +58,6 @@ module surface_header ! y = y0 real(8) :: y0 contains - procedure :: evaluate => y_plane_evaluate procedure :: normal => y_plane_normal end type SurfaceYPlane @@ -76,7 +65,6 @@ module surface_header ! z = z0 real(8) :: z0 contains - procedure :: evaluate => z_plane_evaluate procedure :: normal => z_plane_normal end type SurfaceZPlane @@ -87,7 +75,6 @@ module surface_header real(8) :: C real(8) :: D contains - procedure :: evaluate => plane_evaluate procedure :: normal => plane_normal end type SurfacePlane @@ -97,7 +84,6 @@ module surface_header real(8) :: z0 real(8) :: r contains - procedure :: evaluate => x_cylinder_evaluate procedure :: normal => x_cylinder_normal end type SurfaceXCylinder @@ -107,7 +93,6 @@ module surface_header real(8) :: z0 real(8) :: r contains - procedure :: evaluate => y_cylinder_evaluate procedure :: normal => y_cylinder_normal end type SurfaceYCylinder @@ -117,7 +102,6 @@ module surface_header real(8) :: y0 real(8) :: r contains - procedure :: evaluate => z_cylinder_evaluate procedure :: normal => z_cylinder_normal end type SurfaceZCylinder @@ -128,7 +112,6 @@ module surface_header real(8) :: z0 real(8) :: r contains - procedure :: evaluate => sphere_evaluate procedure :: normal => sphere_normal end type SurfaceSphere @@ -139,7 +122,6 @@ module surface_header real(8) :: z0 real(8) :: r2 contains - procedure :: evaluate => x_cone_evaluate procedure :: normal => x_cone_normal end type SurfaceXCone @@ -150,7 +132,6 @@ module surface_header real(8) :: z0 real(8) :: r2 contains - procedure :: evaluate => y_cone_evaluate procedure :: normal => y_cone_normal end type SurfaceYCone @@ -161,7 +142,6 @@ module surface_header real(8) :: z0 real(8) :: r2 contains - procedure :: evaluate => z_cone_evaluate procedure :: normal => z_cone_normal end type SurfaceZCone @@ -169,7 +149,6 @@ module surface_header ! 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 :: evaluate => quadric_evaluate procedure :: normal => quadric_normal end type SurfaceQuadric @@ -214,14 +193,6 @@ contains ! SurfaceXPlane Implementation !=============================================================================== - pure function x_plane_evaluate(this, xyz) result(f) - class(SurfaceXPlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - f = xyz(1) - this%x0 - end function x_plane_evaluate - pure function x_plane_normal(this, xyz) result(uvw) class(SurfaceXPlane), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -234,14 +205,6 @@ contains ! SurfaceYPlane Implementation !=============================================================================== - pure function y_plane_evaluate(this, xyz) result(f) - class(SurfaceYPlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - f = xyz(2) - this%y0 - end function y_plane_evaluate - pure function y_plane_normal(this, xyz) result(uvw) class(SurfaceYPlane), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -254,16 +217,6 @@ contains ! SurfaceZPlane Implementation !=============================================================================== - pure function z_plane_evaluate(this, xyz) result(f) - - class(SurfaceZPlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - f = xyz(3) - this%z0 - - end function z_plane_evaluate - pure function z_plane_normal(this, xyz) result(uvw) class(SurfaceZPlane), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -276,14 +229,6 @@ contains ! SurfacePlane Implementation !=============================================================================== - pure function plane_evaluate(this, xyz) result(f) - class(SurfacePlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - f = this%A*xyz(1) + this%B*xyz(2) + this%C*xyz(3) - this%D - end function plane_evaluate - pure function plane_normal(this, xyz) result(uvw) class(SurfacePlane), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -296,18 +241,6 @@ contains ! SurfaceXCylinder Implementation !=============================================================================== - pure function x_cylinder_evaluate(this, xyz) result(f) - class(SurfaceXCylinder), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - real(8) :: y, z - - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - f = y*y + z*z - this%r*this%r - end function x_cylinder_evaluate - pure function x_cylinder_normal(this, xyz) result(uvw) class(SurfaceXCylinder), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -322,18 +255,6 @@ contains ! SurfaceYCylinder Implementation !=============================================================================== - pure function y_cylinder_evaluate(this, xyz) result(f) - class(SurfaceYCylinder), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - real(8) :: x, z - - x = xyz(1) - this%x0 - z = xyz(3) - this%z0 - f = x*x + z*z - this%r*this%r - end function y_cylinder_evaluate - pure function y_cylinder_normal(this, xyz) result(uvw) class(SurfaceYCylinder), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -348,18 +269,6 @@ contains ! SurfaceZCylinder Implementation !=============================================================================== - pure function z_cylinder_evaluate(this, xyz) result(f) - class(SurfaceZCylinder), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - real(8) :: x, y - - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - f = x*x + y*y - this%r*this%r - end function z_cylinder_evaluate - pure function z_cylinder_normal(this, xyz) result(uvw) class(SurfaceZCylinder), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -374,19 +283,6 @@ contains ! SurfaceSphere Implementation !=============================================================================== - pure function sphere_evaluate(this, xyz) result(f) - class(SurfaceSphere), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - real(8) :: x, y, z - - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - f = x*x + y*y + z*z - this%r*this%r - end function sphere_evaluate - pure function sphere_normal(this, xyz) result(uvw) class(SurfaceSphere), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -399,19 +295,6 @@ contains ! SurfaceXCone Implementation !=============================================================================== - pure function x_cone_evaluate(this, xyz) result(f) - class(SurfaceXCone), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - real(8) :: x, y, z - - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - f = y*y + z*z - this%r2*x*x - end function x_cone_evaluate - pure function x_cone_normal(this, xyz) result(uvw) class(SurfaceXCone), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -426,19 +309,6 @@ contains ! SurfaceYCone Implementation !=============================================================================== - pure function y_cone_evaluate(this, xyz) result(f) - class(SurfaceYCone), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - real(8) :: x, y, z - - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - f = x*x + z*z - this%r2*y*y - end function y_cone_evaluate - pure function y_cone_normal(this, xyz) result(uvw) class(SurfaceYCone), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -453,19 +323,6 @@ contains ! SurfaceZCone Implementation !=============================================================================== - pure function z_cone_evaluate(this, xyz) result(f) - class(SurfaceZCone), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - real(8) :: x, y, z - - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - f = x*x + y*y - this%r2*z*z - end function z_cone_evaluate - pure function z_cone_normal(this, xyz) result(uvw) class(SurfaceZCone), intent(in) :: this real(8), intent(in) :: xyz(3) @@ -480,18 +337,6 @@ contains ! SurfaceQuadric Implementation !=============================================================================== - pure function quadric_evaluate(this, xyz) result(f) - class(SurfaceQuadric), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8) :: f - - associate (x => xyz(1), y => xyz(2), z => xyz(3)) - f = x*(this%A*x + this%D*y + this%G) + & - y*(this%B*y + this%E*z + this%H) + & - z*(this%C*z + this%F*x + this%J) + this%K - end associate - end function quadric_evaluate - pure function quadric_normal(this, xyz) result(uvw) class(SurfaceQuadric), intent(in) :: this real(8), intent(in) :: xyz(3) diff --git a/src/tracking.F90 b/src/tracking.F90 index 65769a7f18..4f0f110bc7 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -5,7 +5,7 @@ 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 + check_cell_overlap, surface_periodic_c use message_passing use mgxs_header use nuclide_header @@ -290,7 +290,6 @@ contains real(8) :: v ! y-component of direction real(8) :: w ! z-component of direction real(8) :: norm ! "norm" of surface normal - real(8) :: d ! distance between point and plane real(8) :: xyz(3) ! Saved global coordinate integer :: i_surface ! index in surfaces logical :: rotational ! if rotational periodic BC applied @@ -412,64 +411,8 @@ contains p % coord(1) % xyz = xyz end if - rotational = .false. - select type (surf) - type is (SurfaceXPlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfaceXPlane) - p % coord(1) % xyz(1) = opposite % x0 - type is (SurfaceYPlane) - rotational = .true. - - ! Rotate direction - u = p % coord(1) % uvw(1) - v = p % coord(1) % uvw(2) - p % coord(1) % uvw(1) = v - p % coord(1) % uvw(2) = -u - - ! Rotate position - p % coord(1) % xyz(1) = surf % x0 + p % coord(1) % xyz(2) - opposite % y0 - p % coord(1) % xyz(2) = opposite % y0 - end select - - type is (SurfaceYPlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfaceYPlane) - p % coord(1) % xyz(2) = opposite % y0 - type is (SurfaceXPlane) - rotational = .true. - - ! Rotate direction - u = p % coord(1) % uvw(1) - v = p % coord(1) % uvw(2) - p % coord(1) % uvw(1) = -v - p % coord(1) % uvw(2) = u - - ! Rotate position - p % coord(1) % xyz(2) = surf % y0 + p % coord(1) % xyz(1) - opposite % x0 - p % coord(1) % xyz(1) = opposite % x0 - end select - - type is (SurfaceZPlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfaceZPlane) - p % coord(1) % xyz(3) = opposite % z0 - end select - - type is (SurfacePlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfacePlane) - ! Get surface normal for opposite plane - xyz(:) = opposite % normal(p % coord(1) % xyz) - - ! Determine distance to plane - norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3) - d = opposite % evaluate(p % coord(1) % xyz) / norm - - ! Move particle along normal vector based on distance - p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz - end select - end select + rotational = surface_periodic_c(i_surface-1, surf % i_periodic-1, & + p % coord(1) % xyz, p % coord(1) % uvw) ! Reassign particle's surface if (rotational) then