From 91ccc4da0051a3d1f75286bbc4c2305cd720867f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 8 Aug 2017 18:10:53 -0400 Subject: [PATCH] Implement z planes and cylinders in C++ --- CMakeLists.txt | 5 +- src/geometry.F90 | 47 +++++++- src/input_xml.F90 | 9 ++ src/surface_header.C | 280 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 338 insertions(+), 3 deletions(-) create mode 100644 src/surface_header.C diff --git a/CMakeLists.txt b/CMakeLists.txt index c5150745b..200d836fe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -436,7 +436,10 @@ set(LIBOPENMC_FORTRAN_SRC ) set(LIBOPENMC_CXX_SRC src/random_lcg.h - src/random_lcg.cpp) + src/random_lcg.cpp + src/surface_header.C + src/pugixml/pugixml.hpp + src/pugixml/pugixml.cpp) add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC} ${LIBOPENMC_CXX_SRC}) set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc) add_executable(${program} src/main.F90) diff --git a/src/geometry.F90 b/src/geometry.F90 index ced62e681..bc1eca2ae 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -10,8 +10,33 @@ module geometry use stl_vector, only: VectorInt use string, only: to_str + use, intrinsic :: ISO_C_BINDING + implicit none + interface + function surface_sense_c(surf_ind, xyz, uvw) & + bind(C, name='surface_sense') result(sense) + use ISO_C_BINDING + implicit none + integer(C_INT), value :: surf_ind; + real(C_DOUBLE) :: xyz(3); + real(C_DOUBLE) :: uvw(3); + logical(C_BOOL) :: sense; + end function surface_sense_c + + function surface_distance_c(surf_ind, xyz, uvw, coincident) & + bind(C, name='surface_distance') result(d) + use ISO_C_BINDING + implicit none + integer(C_INT), value :: surf_ind; + real(C_DOUBLE) :: xyz(3); + real(C_DOUBLE) :: uvw(3); + logical(C_BOOL), value :: coincident; + real(C_DOUBLE) :: d; + end function surface_distance_c + end interface + contains !=============================================================================== @@ -28,7 +53,8 @@ contains ! expression using a stack, similar to how a RPN calculator would work. !=============================================================================== - pure function cell_contains(c, p) result(in_cell) + !pure function cell_contains(c, p) result(in_cell) + function cell_contains(c, p) result(in_cell) type(Cell), intent(in) :: c type(Particle), intent(in) :: p logical :: in_cell @@ -40,7 +66,8 @@ contains end if end function cell_contains - pure function simple_cell_contains(c, p) result(in_cell) + !pure function simple_cell_contains(c, p) result(in_cell) + function simple_cell_contains(c, p) result(in_cell) type(Cell), intent(in) :: c type(Particle), intent(in) :: p logical :: in_cell @@ -48,6 +75,7 @@ contains integer :: i integer :: token logical :: actual_sense ! sense of particle wrt surface + logical :: alt_sense in_cell = .true. do i = 1, size(c%rpn) @@ -65,6 +93,15 @@ contains else actual_sense = surfaces(abs(token))%obj%sense(& p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) + alt_sense = surface_sense_c(abs(token)-1, & + p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) + if (actual_sense .neqv. alt_sense) then + write(*, *) surfaces(abs(token)) % obj % id + write(*, *) p % coord (p % n_coord) % xyz + write(*, *) p % coord (p % n_coord) % uvw + write(*, *) actual_sense, alt_sense + call fatal_error("") + end if if (actual_sense .neqv. (token > 0)) then in_cell = .false. exit @@ -484,6 +521,7 @@ contains type(Cell), pointer :: c class(Surface), pointer :: surf class(Lattice), pointer :: lat + real(8) :: d_alt ! inialize distance to infinity (huge) dist = INFINITY @@ -519,6 +557,11 @@ contains ! Calculate distance to surface surf => surfaces(index_surf) % obj d = surf % distance(p % coord(j) % xyz, p % coord(j) % uvw, coincident) + d_alt = surface_distance_c(index_surf-1, p % coord(j) % xyz, & + p % coord(j) % uvw, logical(coincident, kind=C_BOOL)) + if (d /= d_alt) then + write(*, *) d, d_alt + end if ! Check if calculated distance is new minimum if (d < d_surf) then diff --git a/src/input_xml.F90 b/src/input_xml.F90 index bd56422cb..4994ff76b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -48,6 +48,14 @@ module input_xml implicit none save + interface + subroutine read_surfaces(node_ptr) bind(C, name='read_surfaces') + use ISO_C_BINDING + implicit none + type(C_PTR) :: node_ptr + end subroutine read_surfaces + end interface + contains !=============================================================================== @@ -966,6 +974,7 @@ contains ! get pointer to list of xml call get_node_list(root, "surface", node_surf_list) + call read_surfaces(root % ptr) ! Get number of tags n_surfaces = size(node_surf_list) diff --git a/src/surface_header.C b/src/surface_header.C new file mode 100644 index 000000000..08287da88 --- /dev/null +++ b/src/surface_header.C @@ -0,0 +1,280 @@ +#include // For strcmp +#include +#include // For numeric_limits +#include // For fabs +#include "pugixml/pugixml.hpp" + +//============================================================================== +// Constants +//============================================================================== + +const double FP_COINCIDENT = 1e-12; +const double INFTY = std::numeric_limits::max(); + +//============================================================================== +// Global array of surfaces +//============================================================================== + +class Surface; +Surface **surfaces_c; + +//============================================================================== +// SURFACE type defines a first- or second-order surface that can be used to +// construct closed volumes (cells) +//============================================================================== + +class Surface +{ + int id; // Unique ID + int neighbor_pos[], // List of cells on positive side + neighbor_neg[]; // List of cells on negative side + int bc; // Boundary condition + //TODO: swith that zero to a NONE constant. + int i_periodic = 0; // Index of corresponding periodic surface + char name[104]; // User-defined name + +public: + bool sense(double xyz[3], double uvw[3]); + void reflect(double xyz[3], double uvw[3]); + virtual double evaluate(double xyz[3]) = 0; + virtual double distance(double xyz[3], double uvw[3], bool coincident) = 0; + virtual void normal(double xyz[3], double uvw[3]) = 0; +}; + +bool +Surface::sense(double xyz[3], double uvw[3]) +{ + // Evaluate the surface equation at the particle's coordinates to determine + // which side the particle is on. + const double f = evaluate(xyz); + + // Check which side of surface the point is on. + bool s; + if (fabs(f) < FP_COINCIDENT) + { + // Particle may be coincident with this surface. To determine the sense, we + // look at the direction of the particle relative to the surface normal (by + // default in the positive direction) via their dot product. + double norm[3]; + normal(xyz, norm); + s = (uvw[0] * norm[0] + uvw[1] * norm[1] + uvw[2] * norm[2] > 0.0); + } + else + { + s = (f > 0.0); + } + return s; +} + +void +Surface::reflect(double xyz[3], double uvw[3]) +{ + + // Determine projection of direction onto normal and squared magnitude of + // normal. + double norm[3]; + normal(xyz, norm); + const double projection = norm[0]*uvw[0] + norm[1]*uvw[1] + norm[2]*uvw[2]; + const double magnitude = norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2]; + + // Reflect direction according to normal. + uvw[0] -= 2.0 * projection / magnitude * norm[0]; + uvw[1] -= 2.0 * projection / magnitude * norm[1]; + uvw[2] -= 2.0 * projection / magnitude * norm[2]; +} + +//============================================================================== + +class SurfaceZPlane : public Surface +{ + double z0; +public: + SurfaceZPlane(pugi::xml_node surf_node); + double evaluate(double xyz[3]); + double distance(double xyz[3], double uvw[3], bool coincident); + void normal(double xyz[3], double uvw[3]); +}; + +SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) +{ + const char *coeffs = surf_node.attribute("coeffs").value(); + int stat = sscanf(coeffs, "%lf", &z0); + if (stat != 1) + { + std::cout << "Something went wrong reading surface coeffs!" << std::endl; + } +} + +double +SurfaceZPlane::evaluate(double xyz[3]) +{ + double f = xyz[2] - z0; + return f; +} + +double +SurfaceZPlane::distance(double xyz[3], double uvw[3], bool coincident) +{ + double f = z0 - xyz[2]; + double d; + if (coincident or fabs(f) < FP_COINCIDENT or uvw[2] == 0.0) + { + d = INFTY; + } + else + { + d = f / uvw[2]; + if (d < 0.0) d = INFTY; + } + return d; +} + +void +SurfaceZPlane::normal(double xyz[3], double uvw[3]) +{ + uvw[0] = 0.0; + uvw[1] = 0.0; + uvw[2] = 1.0; +} + +//============================================================================== + +class SurfaceZCylinder : public Surface +{ + double x0, y0, r; +public: + SurfaceZCylinder(pugi::xml_node surf_node); + double evaluate(double xyz[3]); + double distance(double xyz[3], double uvw[3], bool coincident); + void normal(double xyz[3], double uvw[3]); +}; + +SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) +{ + const char *coeffs = surf_node.attribute("coeffs").value(); + int stat = sscanf(coeffs, "%lf %lf %lf", &x0, &y0, &r); + if (stat != 3) + { + std::cout << "Something went wrong reading surface coeffs!" << std::endl; + } +} + +double +SurfaceZCylinder::evaluate(double xyz[3]) +{ + const double x = xyz[0] - x0; + const double y = xyz[1] - y0; + double f = x*x + y*y - r*r; + return f; +} + +double +SurfaceZCylinder::distance(double xyz[3], double uvw[3], bool coincident) +{ + double a = 1.0 - uvw[2]*uvw[2]; // u^2 + v^2 + + if (a == 0.0) return INFTY; + + double x = xyz[0] - x0; + double y = xyz[1] - y0; + double k = x*uvw[0] + y*uvw[1]; + double c = x*x + y*y - r*r; + double quad = k*k - a*c; + + if (quad < 0.0) + { + // No intersection with cylinder. + return INFTY; + } + else if (coincident or fabs(c) < FP_COINCIDENT) + { + // Particle is on the cylinder, thus one distance is positive/negative + // and the other is zero. The sign of k determines if we are facing in or + // out. + if (k >= 0.0) return INFTY; + else return (-k + sqrt(quad)) / a; + } + else if (c < 0.0) + { + // Particle is inside the cylinder, thus one distance must be negative + // and one must be positive. The positive distance will be the one with + // negative sign on sqrt(quad) + return (-k + sqrt(quad)) / a; + } + else + { + // Particle is outside the cylinder, thus both distances are either + // positive or negative. If positive, the smaller distance is the one + // with positive sign on sqrt(quad) + + double d = (-k - sqrt(quad)) / a; + if (d < 0.0) return INFTY; + return d; + } +} + +void +SurfaceZCylinder::normal(double xyz[3], double uvw[3]) +{ + uvw[0] = 2.0 * (xyz[0] - x0); + uvw[1] = 2.0 * (xyz[1] - y0); + uvw[2] = 0.0; +} + +//============================================================================== + +extern "C" void +read_surfaces(pugi::xml_node *node) +{ + // Count the number of surfaces. + int n_surfaces = 0; + for (pugi::xml_node surf_node = node->child("surface"); surf_node; + surf_node = surf_node.next_sibling("surface")) + { + n_surfaces += 1; + } + + // Allocate the array of Surface pointers. + surfaces_c = new Surface* [n_surfaces]; + + // Loop over XML surface elements and populate the array. + { + pugi::xml_node surf_node; + int i_surf; + for (surf_node = node->child("surface"), i_surf = 0; surf_node; + surf_node = surf_node.next_sibling("surface"), i_surf++) + { + if (surf_node.attribute("type")) + { + const pugi::char_t *surf_type = surf_node.attribute("type").value(); + if (strcmp(surf_type, "z-cylinder") == 0) + { + surfaces_c[i_surf] = new SurfaceZCylinder(surf_node); + } + else if (strcmp(surf_type, "z-plane") == 0) + { + surfaces_c[i_surf] = new SurfaceZPlane(surf_node); + } + else + { + std::cout << "Call error or handle uppercase here!" << std::endl; + std::cout << surf_type << std::endl; + } + } + } + } +} + +//============================================================================== + +extern "C" bool +surface_sense(int surf_ind, double xyz[3], double uvw[3]) +{ + return surfaces_c[surf_ind]->sense(xyz, uvw); +} + +extern "C" double +surface_distance(int surf_ind, double xyz[3], double uvw[3], bool coincident) +{ + return surfaces_c[surf_ind]->distance(xyz, uvw, coincident); +}