mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 21:55:41 -04:00
Implement C++ periodic BCs
This commit is contained in:
Sterling Harper
parent
163b6e3de5
commit
38aa76c9dc
4 changed files with 145 additions and 222 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -5,6 +5,7 @@
|
|||
#include "pugixml/pugixml.hpp"
|
||||
|
||||
// DEBUGGING
|
||||
#include <typeinfo>
|
||||
#include <iostream>
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -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<PeriodicSurface *>(surf1_gen);
|
||||
Surface *surf2_gen = surfaces_c[surf_ind2];
|
||||
PeriodicSurface *surf2 = dynamic_cast<PeriodicSurface *>(surf2_gen);
|
||||
|
||||
// Call the type-bound methods.
|
||||
return surf2->periodic_translate(surf1, xyz, uvw);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue