diff --git a/CMakeLists.txt b/CMakeLists.txt index 51bdcc4d9..1f5b2168f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -48,14 +48,14 @@ add_definitions(-DMAX_COORD=${maxcoord}) set(MPI_ENABLED FALSE) if($ENV{FC} MATCHES "(mpi[^/]*|ftn)$") message("-- Detected MPI wrapper: $ENV{FC}") - add_definitions(-DMPI) + add_definitions(-DOPENMC_MPI) set(MPI_ENABLED TRUE) endif() # Check for Fortran 2008 MPI interface if(MPI_ENABLED AND mpif08) message("-- Using Fortran 2008 MPI bindings") - add_definitions(-DMPIF08) + add_definitions(-DOPENMC_MPIF08) endif() #=============================================================================== @@ -435,8 +435,15 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/trigger_header.F90 ) set(LIBOPENMC_CXX_SRC + src/error.h + src/hdf5_interface.h + src/random_lcg.cpp src/random_lcg.h - src/random_lcg.cpp) + src/surface.cpp + src/surface.h + src/xml_interface.h + src/pugixml/pugixml.cpp + src/pugixml/pugixml.hpp) 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/openmc/capi/settings.py b/openmc/capi/settings.py index eede6cf47..1063d6463 100644 --- a/openmc/capi/settings.py +++ b/openmc/capi/settings.py @@ -11,8 +11,7 @@ _RUN_MODES = {1: 'fixed source', 5: 'volume'} _dll.openmc_set_seed.argtypes = [c_int64] -_dll.openmc_set_seed.restype = c_int -_dll.openmc_set_seed.errcheck = _error_handler +_dll.openmc_get_seed.restype = c_int64 class _Settings(object): @@ -43,7 +42,7 @@ class _Settings(object): @property def seed(self): - return c_int64.in_dll(_dll, 'seed').value + return _dll.openmc_get_seed() @seed.setter def seed(self, seed): diff --git a/src/api.F90 b/src/api.F90 index eb9a99695..f07e5e38a 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -18,7 +18,7 @@ module openmc_api use initialize, only: openmc_init use particle_header, only: Particle use plot, only: openmc_plot_geometry - use random_lcg, only: seed, openmc_set_seed + use random_lcg, only: openmc_get_seed, openmc_set_seed use settings use simulation_header use source_header, only: openmc_extend_sources, openmc_source_set_strength @@ -60,6 +60,7 @@ module openmc_api public :: openmc_get_filter_next_id public :: openmc_get_material_index public :: openmc_get_nuclide_index + public :: openmc_get_seed public :: openmc_get_tally_index public :: openmc_global_tallies public :: openmc_hard_reset @@ -79,6 +80,7 @@ module openmc_api public :: openmc_plot_geometry public :: openmc_reset public :: openmc_run + public :: openmc_set_seed public :: openmc_simulation_finalize public :: openmc_simulation_init public :: openmc_source_bank @@ -142,7 +144,7 @@ contains run_CE = .true. run_mode = NONE satisfy_triggers = .false. - seed = 1_8 + call openmc_set_seed(DEFAULT_SEED) source_latest = .false. source_separate = .false. source_write = .true. @@ -171,7 +173,7 @@ contains ! Close FORTRAN interface. call h5close_f(err) -#ifdef MPI +#ifdef OPENMC_MPI ! Free all MPI types call MPI_TYPE_FREE(MPI_BANK, err) #endif @@ -228,8 +230,6 @@ contains !=============================================================================== subroutine openmc_hard_reset() bind(C) - integer :: err - ! Reset all tallies and timers call openmc_reset() @@ -238,7 +238,7 @@ contains total_gen = 0 ! Reset the random number generator state - err = openmc_set_seed(1_8) + call openmc_set_seed(DEFAULT_SEED) end subroutine openmc_hard_reset !=============================================================================== diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 9cdb751a4..af8d57a80 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -105,7 +105,7 @@ contains real(8) :: hxyz(3) ! cell dimensions of current ijk cell real(8) :: vol ! volume of cell real(8),allocatable :: source(:,:,:,:) ! tmp source array for entropy -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif @@ -197,7 +197,7 @@ contains end if -#ifdef MPI +#ifdef OPENMC_MPI ! Broadcast full source to all procs call MPI_BCAST(cmfd % cmfd_src, n, MPI_REAL8, 0, mpi_intracomm, mpi_err) #endif @@ -234,7 +234,7 @@ contains real(8) :: norm ! normalization factor logical :: outside ! any source sites outside mesh logical :: in_mesh ! source site is inside mesh -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err #endif @@ -291,7 +291,7 @@ contains if (.not. cmfd_feedback) return ! Broadcast weight factors to all procs -#ifdef MPI +#ifdef OPENMC_MPI call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, & mpi_intracomm, mpi_err) #endif diff --git a/src/constants.F90 b/src/constants.F90 index 331a644f0..3334a4cc7 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -43,9 +43,9 @@ module constants real(8), parameter :: TINY_BIT = 1e-8_8 ! User for precision in geometry - real(8), parameter :: FP_PRECISION = 1e-14_8 - real(8), parameter :: FP_REL_PRECISION = 1e-5_8 - real(8), parameter :: FP_COINCIDENT = 1e-12_8 + real(C_DOUBLE), bind(C, name='FP_PRECISION') :: FP_PRECISION = 1e-14_8 + real(C_DOUBLE), bind(C, name='FP_REL_PRECISION') :: FP_REL_PRECISION = 1e-5_8 + real(C_DOUBLE), bind(C, name='FP_COINCIDENT') :: FP_COINCIDENT = 1e-12_8 ! Maximum number of collisions/crossings integer, parameter :: MAX_EVENTS = 1000000 @@ -95,11 +95,10 @@ module constants ! GEOMETRY-RELATED CONSTANTS ! Boundary conditions - integer, parameter :: & - BC_TRANSMIT = 0, & ! Transmission boundary condition (default) - BC_VACUUM = 1, & ! Vacuum boundary condition - BC_REFLECT = 2, & ! Reflecting boundary condition - BC_PERIODIC = 3 ! Periodic boundary condition + integer(C_INT), bind(C, name="BC_TRANSMIT") :: BC_TRANSMIT + integer(C_INT), bind(C, name="BC_VACUUM") :: BC_VACUUM + integer(C_INT), bind(C, name="BC_REFLECT") :: BC_REFLECT + integer(C_INT), bind(C, name="BC_PERIODIC") :: BC_PERIODIC ! Logical operators for cell definitions integer, parameter :: & diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 27ae4063e..7a11d1fb4 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -42,11 +42,11 @@ contains type(Bank), save, allocatable :: & & temp_sites(:) ! local array of extra sites on each node -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code integer(8) :: n ! number of sites to send/recv integer :: neighbor ! processor to send/recv data from -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 type(MPI_Request) :: request(20) #else integer :: request(20) ! communication request for send/recving sites @@ -66,7 +66,7 @@ contains ! fission bank its own sites starts in order to ensure reproducibility by ! skipping ahead to the proper seed. -#ifdef MPI +#ifdef OPENMC_MPI start = 0_8 call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & mpi_intracomm, mpi_err) @@ -148,7 +148,7 @@ contains ! neighboring processors, we have to perform an ALLGATHER to determine the ! indices for all processors -#ifdef MPI +#ifdef OPENMC_MPI ! First do an exclusive scan to get the starting indices for start = 0_8 call MPI_EXSCAN(index_temp, start, 1, MPI_INTEGER8, MPI_SUM, & @@ -191,7 +191,7 @@ contains call time_bank_sample % stop() call time_bank_sendrecv % start() -#ifdef MPI +#ifdef OPENMC_MPI ! ========================================================================== ! SEND BANK SITES TO NEIGHBORS @@ -343,14 +343,14 @@ contains subroutine calculate_generation_keff() real(8) :: keff_reduced -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif ! Get keff for this generation by subtracting off the starting value keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation -#ifdef MPI +#ifdef OPENMC_MPI ! Combine values across all processors call MPI_ALLREDUCE(keff_generation, keff_reduced, 1, MPI_REAL8, & MPI_SUM, mpi_intracomm, mpi_err) @@ -584,7 +584,7 @@ contains real(8) :: total ! total weight in source bank logical :: sites_outside ! were there sites outside the ufs mesh? -#ifdef MPI +#ifdef OPENMC_MPI integer :: n ! total number of ufs mesh cells integer :: mpi_err ! MPI error code #endif @@ -608,7 +608,7 @@ contains call fatal_error("Source sites outside of the UFS mesh!") end if -#ifdef MPI +#ifdef OPENMC_MPI ! Send source fraction to all processors n = product(m % dimension) call MPI_BCAST(source_frac, n, MPI_REAL8, 0, mpi_intracomm, mpi_err) diff --git a/src/error.F90 b/src/error.F90 index 27004d729..0a5af302d 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -128,7 +128,7 @@ contains integer :: line_wrap ! length of line integer :: length ! length of message integer :: indent ! length of indentation -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err #endif @@ -180,7 +180,7 @@ contains end if end do -#ifdef MPI +#ifdef OPENMC_MPI ! Abort MPI call MPI_ABORT(mpi_intracomm, code, mpi_err) #endif @@ -194,6 +194,14 @@ contains end subroutine fatal_error + subroutine fatal_error_from_c(message, message_len) bind(C) + integer(C_INT), intent(in), value :: message_len + character(kind=C_CHAR), intent(in) :: message(message_len) + character(message_len+1) :: message_out + write(message_out, *) message + call fatal_error(message_out) + end subroutine + !=============================================================================== ! WRITE_MESSAGE displays an informational message to the log file and the ! standard output stream. diff --git a/src/error.h b/src/error.h new file mode 100644 index 000000000..4c3373b3e --- /dev/null +++ b/src/error.h @@ -0,0 +1,34 @@ +#ifndef ERROR_H +#define ERROR_H + +#include +#include +#include + + +namespace openmc { + + +extern "C" void fatal_error_from_c(const char *message, int message_len); + + +void fatal_error(const char *message) +{ + fatal_error_from_c(message, strlen(message)); +} + + +void fatal_error(const std::string &message) +{ + fatal_error_from_c(message.c_str(), message.length()); +} + + +void fatal_error(const std::stringstream &message) +{ + std::string out {message.str()}; + fatal_error_from_c(out.c_str(), out.length()); +} + +} // namespace openmc +#endif // ERROR_H diff --git a/src/geometry.F90 b/src/geometry.F90 index ced62e681..0bacd7b5d 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -10,6 +10,8 @@ module geometry use stl_vector, only: VectorInt use string, only: to_str + use, intrinsic :: ISO_C_BINDING + implicit none contains @@ -63,7 +65,7 @@ contains in_cell = .false. exit else - actual_sense = surfaces(abs(token))%obj%sense(& + actual_sense = surfaces(abs(token)) % sense( & p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) if (actual_sense .neqv. (token > 0)) then in_cell = .false. @@ -112,7 +114,7 @@ contains elseif (-token == p%surface) then stack(i_stack) = .false. else - actual_sense = surfaces(abs(token))%obj%sense(& + actual_sense = surfaces(abs(token)) % sense( & p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) stack(i_stack) = (actual_sense .eqv. (token > 0)) end if @@ -480,9 +482,9 @@ 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 class(Lattice), pointer :: lat ! inialize distance to infinity (huge) @@ -517,8 +519,8 @@ contains if (index_surf >= OP_UNION) cycle ! Calculate distance to surface - surf => surfaces(index_surf) % obj - d = surf % distance(p % coord(j) % xyz, p % coord(j) % uvw, coincident) + d = surfaces(index_surf) % distance(p % coord(j) % xyz, & + p % coord(j) % uvw, logical(coincident, kind=C_BOOL)) ! Check if calculated distance is new minimum if (d < d_surf) then @@ -737,9 +739,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 surfaces(abs(level_surf_cross)) % normal(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) @@ -804,15 +805,15 @@ contains ! Copy positive neighbors to Surface instance n = neighbor_pos(i)%size() if (n > 0) then - allocate(surfaces(i)%obj%neighbor_pos(n)) - surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:n) + allocate(surfaces(i)%neighbor_pos(n)) + surfaces(i)%neighbor_pos(:) = neighbor_pos(i)%data(1:n) end if ! Copy negative neighbors to Surface instance n = neighbor_neg(i)%size() if (n > 0) then - allocate(surfaces(i)%obj%neighbor_neg(n)) - surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:n) + allocate(surfaces(i)%neighbor_neg(n)) + surfaces(i)%neighbor_neg(:) = neighbor_neg(i)%data(1:n) end if end do diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index e8c39b923..410149d9e 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -124,7 +124,7 @@ contains ! Setup file access property list with parallel I/O access call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) #ifdef PHDF5 -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, & MPI_INFO_NULL%MPI_VAL, hdf5_err) #else @@ -174,7 +174,7 @@ contains ! Setup file access property list with parallel I/O access call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) #ifdef PHDF5 -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, & MPI_INFO_NULL%MPI_VAL, hdf5_err) #else diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h new file mode 100644 index 000000000..7ec9ada9b --- /dev/null +++ b/src/hdf5_interface.h @@ -0,0 +1,91 @@ +#ifndef HDF5_INTERFACE_H +#define HDF5_INTERFACE_H + +#include +#include + +#include "hdf5.h" + +#include "error.h" + + +namespace openmc { + + +hid_t +create_group(hid_t parent_id, char const *name) +{ + hid_t out = H5Gcreate(parent_id, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + if (out < 0) { + std::string err_msg{"Failed to create HDF5 group \""}; + err_msg += name; + err_msg += "\""; + fatal_error(err_msg); + } + return out; +} + + +hid_t +create_group(hid_t parent_id, const std::string &name) +{ + return create_group(parent_id, name.c_str()); +} + + +void +close_group(hid_t group_id) +{ + herr_t err = H5Gclose(group_id); + if (err < 0) { + fatal_error("Failed to close HDF5 group"); + } +} + + +template void +write_double_1D(hid_t group_id, char const *name, + std::array &buffer) +{ + hsize_t dims[1]{array_len}; + hid_t dataspace = H5Screate_simple(1, dims, NULL); + + hid_t dataset = H5Dcreate(group_id, name, H5T_NATIVE_DOUBLE, dataspace, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, + &buffer[0]); + + H5Sclose(dataspace); + H5Dclose(dataset); +} + + +void +write_string(hid_t group_id, char const *name, char const *buffer) +{ + size_t buffer_len = strlen(buffer); + hid_t datatype = H5Tcopy(H5T_C_S1); + H5Tset_size(datatype, buffer_len); + + hid_t dataspace = H5Screate(H5S_SCALAR); + + hid_t dataset = H5Dcreate(group_id, name, datatype, dataspace, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + H5Dwrite(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer); + + H5Tclose(datatype); + H5Sclose(dataspace); + H5Dclose(dataset); +} + + +void +write_string(hid_t group_id, char const *name, const std::string &buffer) +{ + write_string(group_id, name, buffer.c_str()); +} + +} // namespace openmc +#endif //HDF5_INTERFACE_H diff --git a/src/initialize.F90 b/src/initialize.F90 index dd7894f96..7b699b9cb 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -44,7 +44,6 @@ contains subroutine openmc_init(intracomm) bind(C) integer, intent(in), optional :: intracomm ! MPI intracommunicator - integer :: err #ifdef _OPENMP character(MAX_WORD_LEN) :: envvar #endif @@ -52,8 +51,8 @@ contains ! Copy the communicator to a new variable. This is done to avoid changing ! the signature of this subroutine. If MPI is being used but no communicator ! was passed, assume MPI_COMM_WORLD. -#ifdef MPI -#ifdef MPIF08 +#ifdef OPENMC_MPI +#ifdef OPENMC_MPIF08 type(MPI_Comm), intent(in) :: comm ! MPI intracommunicator if (present(intracomm)) then comm % MPI_VAL = intracomm @@ -74,7 +73,7 @@ contains call time_total%start() call time_initialize%start() -#ifdef MPI +#ifdef OPENMC_MPI ! Setup MPI call initialize_mpi(comm) #endif @@ -95,7 +94,7 @@ contains ! Initialize random number generator -- if the user specifies a seed, it ! will be re-initialized later - err = openmc_set_seed(DEFAULT_SEED) + call openmc_set_seed(DEFAULT_SEED) ! Read XML input files call read_input_xml() @@ -108,7 +107,7 @@ contains end subroutine openmc_init -#ifdef MPI +#ifdef OPENMC_MPI !=============================================================================== ! INITIALIZE_MPI starts up the Message Passing Interface (MPI) and determines ! the number of processors the problem is being run with as well as the rank of @@ -116,7 +115,7 @@ contains !=============================================================================== subroutine initialize_mpi(intracomm) -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 type(MPI_Comm), intent(in) :: intracomm ! MPI intracommunicator #else integer, intent(in) :: intracomm ! MPI intracommunicator @@ -124,7 +123,7 @@ contains integer :: mpi_err ! MPI error code integer :: bank_blocks(5) ! Count for each datatype -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 type(MPI_Datatype) :: bank_types(5) #else integer :: bank_types(5) ! Datatypes diff --git a/src/input_xml.F90 b/src/input_xml.F90 index bd56422cb..0b7790a51 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 !=============================================================================== @@ -341,7 +349,7 @@ contains ! Copy random number seed if specified if (check_for_node(root, "seed")) then call get_node_value(root, "seed", seed) - err = openmc_set_seed(seed) + call openmc_set_seed(seed) end if ! Number of bins for logarithmic grid @@ -911,28 +919,20 @@ contains integer :: id integer :: univ_id integer :: n_cells_in_univ - integer :: coeffs_reqd - integer :: i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax - real(8) :: xmin, xmax, ymin, ymax, zmin, zmax integer, allocatable :: temp_int_array(:) real(8) :: phi, theta, psi - real(8), allocatable :: coeffs(:) logical :: file_exists logical :: boundary_exists character(MAX_LINE_LEN) :: filename - character(MAX_WORD_LEN) :: word character(MAX_WORD_LEN), allocatable :: sarray(:) character(:), allocatable :: region_spec type(Cell), pointer :: c - class(Surface), pointer :: s class(Lattice), pointer :: lat type(XMLDocument) :: doc type(XMLNode) :: root type(XMLNode) :: node_cell - type(XMLNode) :: node_surf type(XMLNode) :: node_lat type(XMLNode), allocatable :: node_cell_list(:) - type(XMLNode), allocatable :: node_surf_list(:) type(XMLNode), allocatable :: node_rlat_list(:) type(XMLNode), allocatable :: node_hlat_list(:) type(VectorInt) :: tokens @@ -964,218 +964,18 @@ contains ! applied to a surface boundary_exists = .false. - ! 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) - - ! Check for no surfaces - if (n_surfaces == 0) then - call fatal_error("No surfaces found in geometry.xml!") - end if - - xmin = INFINITY - xmax = -INFINITY - ymin = INFINITY - ymax = -INFINITY - zmin = INFINITY - zmax = -INFINITY - - ! Allocate cells array + ! Allocate surfaces array allocate(surfaces(n_surfaces)) do i = 1, n_surfaces - ! Get pointer to i-th surface node - node_surf = node_surf_list(i) + surfaces(i) % ptr = surface_pointer_c(i - 1); - ! Copy and interpret surface type - word = '' - if (check_for_node(node_surf, "type")) & - call get_node_value(node_surf, "type", word) - select case(to_lower(word)) - case ('x-plane') - coeffs_reqd = 1 - allocate(SurfaceXPlane :: surfaces(i)%obj) - case ('y-plane') - coeffs_reqd = 1 - allocate(SurfaceYPlane :: surfaces(i)%obj) - case ('z-plane') - coeffs_reqd = 1 - allocate(SurfaceZPlane :: surfaces(i)%obj) - case ('plane') - coeffs_reqd = 4 - allocate(SurfacePlane :: surfaces(i)%obj) - case ('x-cylinder') - coeffs_reqd = 3 - allocate(SurfaceXCylinder :: surfaces(i)%obj) - case ('y-cylinder') - coeffs_reqd = 3 - allocate(SurfaceYCylinder :: surfaces(i)%obj) - case ('z-cylinder') - coeffs_reqd = 3 - allocate(SurfaceZCylinder :: surfaces(i)%obj) - case ('sphere') - coeffs_reqd = 4 - allocate(SurfaceSphere :: surfaces(i)%obj) - case ('x-cone') - coeffs_reqd = 4 - allocate(SurfaceXCone :: surfaces(i)%obj) - case ('y-cone') - coeffs_reqd = 4 - allocate(SurfaceYCone :: surfaces(i)%obj) - case ('z-cone') - coeffs_reqd = 4 - allocate(SurfaceZCone :: surfaces(i)%obj) - case ('quadric') - coeffs_reqd = 10 - allocate(SurfaceQuadric :: surfaces(i)%obj) - case default - call fatal_error("Invalid surface type: " // trim(word)) - end select + if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true. - s => surfaces(i)%obj - - ! Copy data into cells - if (check_for_node(node_surf, "id")) then - call get_node_value(node_surf, "id", s%id) - else - call fatal_error("Must specify id of surface in geometry XML file.") - end if - - ! Check to make sure 'id' hasn't been used - if (surface_dict % has(s%id)) then - call fatal_error("Two or more surfaces use the same unique ID: " & - // to_str(s%id)) - end if - - ! Copy surface name - if (check_for_node(node_surf, "name")) then - call get_node_value(node_surf, "name", s%name) - end if - - ! Check to make sure that the proper number of coefficients - ! have been specified for the given type of surface. Then copy - ! surface coordinates. - - n = node_word_count(node_surf, "coeffs") - if (n < coeffs_reqd) then - call fatal_error("Not enough coefficients specified for surface: " & - // trim(to_str(s%id))) - elseif (n > coeffs_reqd) then - call fatal_error("Too many coefficients specified for surface: " & - // trim(to_str(s%id))) - end if - - allocate(coeffs(n)) - call get_node_array(node_surf, "coeffs", coeffs) - - select type(s) - type is (SurfaceXPlane) - s%x0 = coeffs(1) - - ! Determine outer surfaces - xmin = min(xmin, s % x0) - xmax = max(xmax, s % x0) - if (xmin == s % x0) i_xmin = i - if (xmax == s % x0) i_xmax = i - type is (SurfaceYPlane) - s%y0 = coeffs(1) - - ! Determine outer surfaces - ymin = min(ymin, s % y0) - ymax = max(ymax, s % y0) - if (ymin == s % y0) i_ymin = i - if (ymax == s % y0) i_ymax = i - type is (SurfaceZPlane) - s%z0 = coeffs(1) - - ! Determine outer surfaces - zmin = min(zmin, s % z0) - zmax = max(zmax, s % z0) - if (zmin == s % z0) i_zmin = i - if (zmax == s % z0) i_zmax = i - type is (SurfacePlane) - s%A = coeffs(1) - s%B = coeffs(2) - s%C = coeffs(3) - s%D = coeffs(4) - type is (SurfaceXCylinder) - s%y0 = coeffs(1) - s%z0 = coeffs(2) - s%r = coeffs(3) - type is (SurfaceYCylinder) - s%x0 = coeffs(1) - s%z0 = coeffs(2) - s%r = coeffs(3) - type is (SurfaceZCylinder) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%r = coeffs(3) - type is (SurfaceSphere) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%z0 = coeffs(3) - s%r = coeffs(4) - type is (SurfaceXCone) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%z0 = coeffs(3) - s%r2 = coeffs(4) - type is (SurfaceYCone) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%z0 = coeffs(3) - s%r2 = coeffs(4) - type is (SurfaceZCone) - s%x0 = coeffs(1) - s%y0 = coeffs(2) - s%z0 = coeffs(3) - s%r2 = coeffs(4) - type is (SurfaceQuadric) - s%A = coeffs(1) - s%B = coeffs(2) - s%C = coeffs(3) - s%D = coeffs(4) - s%E = coeffs(5) - s%F = coeffs(6) - s%G = coeffs(7) - s%H = coeffs(8) - s%J = coeffs(9) - s%K = coeffs(10) - end select - - ! No longer need coefficients - deallocate(coeffs) - - ! Boundary conditions - word = '' - if (check_for_node(node_surf, "boundary")) & - call get_node_value(node_surf, "boundary", word) - select case (to_lower(word)) - case ('transmission', 'transmit', '') - s%bc = BC_TRANSMIT - case ('vacuum') - s%bc = BC_VACUUM - boundary_exists = .true. - case ('reflective', 'reflect', 'reflecting') - s%bc = BC_REFLECT - boundary_exists = .true. - case ('periodic') - s%bc = BC_PERIODIC - boundary_exists = .true. - - ! Check for specification of periodic surface - if (check_for_node(node_surf, "periodic_surface_id")) then - call get_node_value(node_surf, "periodic_surface_id", & - s % i_periodic) - end if - case default - call fatal_error("Unknown boundary condition '" // trim(word) // & - &"' specified on surface " // trim(to_str(s%id))) - end select ! Add surface to dictionary - call surface_dict % set(s%id, i) + call surface_dict % set(surfaces(i) % id(), i) end do ! Check to make sure a boundary condition was applied to at least one @@ -1186,76 +986,6 @@ contains end if end if - ! Determine opposite side for periodic boundaries - do i = 1, size(surfaces) - if (surfaces(i) % obj % bc == BC_PERIODIC) then - select type (surf => surfaces(i) % obj) - type is (SurfaceXPlane) - if (surf % i_periodic == NONE) then - if (i == i_xmin) then - surf % i_periodic = i_xmax - elseif (i == i_xmax) then - surf % i_periodic = i_xmin - else - call fatal_error("Periodic boundary condition applied to & - &interior surface.") - end if - else - surf % i_periodic = surface_dict % get(surf % i_periodic) - end if - - type is (SurfaceYPlane) - if (surf % i_periodic == NONE) then - if (i == i_ymin) then - surf % i_periodic = i_ymax - elseif (i == i_ymax) then - surf % i_periodic = i_ymin - else - call fatal_error("Periodic boundary condition applied to & - &interior surface.") - end if - else - surf % i_periodic = surface_dict % get(surf % i_periodic) - end if - - type is (SurfaceZPlane) - if (surf % i_periodic == NONE) then - if (i == i_zmin) then - surf % i_periodic = i_zmax - elseif (i == i_zmax) then - surf % i_periodic = i_zmin - else - call fatal_error("Periodic boundary condition applied to & - &interior surface.") - end if - else - surf % i_periodic = surface_dict % get(surf % i_periodic) - end if - - type is (SurfacePlane) - if (surf % i_periodic == NONE) then - call fatal_error("No matching periodic surface specified for & - &periodic boundary condition on surface " // & - trim(to_str(surf % id)) // ".") - else - surf % i_periodic = surface_dict % get(surf % i_periodic) - end if - - class default - call fatal_error("Periodic boundary condition applied to & - &non-planar surface.") - end select - - ! Make sure opposite surface is also periodic - associate (surf => surfaces(i) % obj) - if (surfaces(surf % i_periodic) % obj % bc /= BC_PERIODIC) then - call fatal_error("Could not find matching surface for periodic & - &boundary on surface " // trim(to_str(surf % id)) // ".") - end if - end associate - end if - end do - ! ========================================================================== ! READ CELLS FROM GEOMETRY.XML diff --git a/src/main.F90 b/src/main.F90 index d8e52643c..120bdd1e0 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -9,13 +9,13 @@ program main implicit none -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif ! Initialize run -- when run with MPI, pass communicator -#ifdef MPI -#ifdef MPIF08 +#ifdef OPENMC_MPI +#ifdef OPENMC_MPIF08 call openmc_init(MPI_COMM_WORLD % MPI_VAL) #else call openmc_init(MPI_COMM_WORLD) @@ -39,7 +39,7 @@ program main ! finalize run call openmc_finalize() -#ifdef MPI +#ifdef OPENMC_MPI ! If MPI is in use and enabled, terminate it call MPI_FINALIZE(mpi_err) #endif diff --git a/src/mesh.F90 b/src/mesh.F90 index e997890a4..790dc7d88 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -34,7 +34,7 @@ contains integer :: n ! number of energy groups / size integer :: mesh_bin ! mesh bin integer :: e_bin ! energy bin -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif logical :: outside ! was any site outside mesh? @@ -86,7 +86,7 @@ contains cnt_(e_bin, mesh_bin) = cnt_(e_bin, mesh_bin) + bank_array(i) % wgt end do FISSION_SITES -#ifdef MPI +#ifdef OPENMC_MPI ! collect values from all processors n = size(cnt_) call MPI_REDUCE(cnt_, cnt, n, MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err) diff --git a/src/message_passing.F90 b/src/message_passing.F90 index 0f2b94d1a..7391a632c 100644 --- a/src/message_passing.F90 +++ b/src/message_passing.F90 @@ -1,7 +1,7 @@ module message_passing -#ifdef MPI -#ifdef MPIF08 +#ifdef OPENMC_MPI +#ifdef OPENMC_MPIF08 use mpi_f08 #else use mpi @@ -16,7 +16,7 @@ module message_passing integer :: rank = 0 ! rank of process logical :: master = .true. ! master process? logical :: mpi_enabled = .false. ! is MPI in use and initialized? -#ifdef MPIF08 +#ifdef OPENMC_MPIF08 type(MPI_Datatype) :: MPI_BANK ! MPI datatype for fission bank type(MPI_Comm) :: mpi_intracomm ! MPI intra-communicator #else diff --git a/src/output.F90 b/src/output.F90 index 9a144721d..316fe9594 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -88,7 +88,7 @@ contains ! Write the date and time write(UNIT=OUTPUT_UNIT, FMT='(9X,"Date/Time | ",A)') time_stamp() -#ifdef MPI +#ifdef OPENMC_MPI ! Write number of processors write(UNIT=OUTPUT_UNIT, FMT='(5X,"MPI Processes | ",A)') & trim(to_str(n_procs)) @@ -259,7 +259,7 @@ contains ! Print surface if (p % surface /= NONE) then - write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%obj%id, p % surface)) + write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%id(), p % surface)) end if ! Display weight, energy, grid index, and interpolation factor diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 839bd729a..d64717cf4 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -4,8 +4,6 @@ module random_lcg implicit none - integer(C_INT64_T), bind(C) :: seed - interface function prn() result(pseudo_rn) bind(C) use ISO_C_BINDING @@ -38,11 +36,16 @@ module random_lcg integer(C_INT), value :: n end subroutine prn_set_stream - function openmc_set_seed(new_seed) result(err) bind(C) + function openmc_get_seed() result(seed) bind(C) + use ISO_C_BINDING + implicit none + integer(C_INT64_T) :: seed + end function openmc_get_seed + + subroutine openmc_set_seed(new_seed) bind(C) use ISO_C_BINDING implicit none integer(C_INT64_T), value :: new_seed - integer(C_INT) :: err - end function openmc_set_seed + end subroutine openmc_set_seed end interface end module random_lcg diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index 047bda8e5..2a57316cb 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -2,6 +2,9 @@ #include +namespace openmc { + + // Constants extern "C" const int N_STREAMS {5}; extern "C" const int STREAM_TRACKING {0}; @@ -130,7 +133,9 @@ prn_set_stream(int i) // API FUNCTIONS //============================================================================== -extern "C" int +extern "C" int64_t openmc_get_seed() {return seed;} + +extern "C" void openmc_set_seed(int64_t new_seed) { seed = new_seed; @@ -141,5 +146,6 @@ openmc_set_seed(int64_t new_seed) } prn_set_stream(STREAM_TRACKING); } - return 0; } + +} // namespace openmc diff --git a/src/random_lcg.h b/src/random_lcg.h index 04191254f..29b3fca47 100644 --- a/src/random_lcg.h +++ b/src/random_lcg.h @@ -3,6 +3,9 @@ #include + +namespace openmc { + //============================================================================== // Module constants. //============================================================================== @@ -74,12 +77,19 @@ extern "C" void prn_set_stream(int n); // API FUNCTIONS //============================================================================== +//============================================================================== +//! Get OpenMC's master seed. +//============================================================================== + +extern "C" int64_t openmc_get_seed(); + //============================================================================== //! Set OpenMC's master seed. //! @param new_seed The master seed. All other seeds will be derived from this //! one. //============================================================================== -extern "C" int openmc_set_seed(int64_t new_seed); +extern "C" void openmc_set_seed(int64_t new_seed); +} // namespace openmc #endif // RANDOM_LCG_H diff --git a/src/simulation.F90 b/src/simulation.F90 index ef29e7832..24e9b29fc 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -320,7 +320,7 @@ contains subroutine finalize_batch() -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code #endif @@ -342,7 +342,7 @@ contains ! Check_triggers if (master) call check_triggers() -#ifdef MPI +#ifdef OPENMC_MPI call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, & mpi_intracomm, mpi_err) #endif @@ -469,7 +469,7 @@ contains subroutine openmc_simulation_finalize() bind(C) integer :: i ! loop index -#ifdef MPI +#ifdef OPENMC_MPI integer :: n ! size of arrays integer :: mpi_err ! MPI error code integer(8) :: temp @@ -494,7 +494,7 @@ contains ! Increment total number of generations total_gen = total_gen + current_batch*gen_per_batch -#ifdef MPI +#ifdef OPENMC_MPI ! Broadcast tally results so that each process has access to results if (allocated(tallies)) then do i = 1, size(tallies) diff --git a/src/source.F90 b/src/source.F90 index 557120e97..5beecd887 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -1,7 +1,7 @@ module source use hdf5, only: HID_T -#ifdef MPI +#ifdef OPENMC_MPI use message_passing #endif diff --git a/src/state_point.F90 b/src/state_point.F90 index f7a892898..980a7333c 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -26,7 +26,7 @@ module state_point use mgxs_header, only: nuclides_MG use nuclide_header, only: nuclides use output, only: time_stamp - use random_lcg, only: seed + use random_lcg, only: openmc_get_seed, openmc_set_seed use settings use simulation_header use string, only: to_str, count_digits, zero_padded, to_f_string @@ -97,7 +97,7 @@ contains call write_attribute(file_id, "path", path_input) ! Write out random number seed - call write_dataset(file_id, "seed", seed) + call write_dataset(file_id, "seed", openmc_get_seed()) ! Write run information if (run_CE) then @@ -523,7 +523,7 @@ contains integer(HID_T) :: tallies_group, tally_group real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results real(8), target :: global_temp(3,N_GLOBAL_TALLIES) -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code real(8) :: dummy ! temporary receive buffer for non-root reduces #endif @@ -543,7 +543,7 @@ contains end if -#ifdef MPI +#ifdef OPENMC_MPI ! Reduce global tallies n_bins = size(global_tallies) call MPI_REDUCE(global_tallies, global_temp, n_bins, MPI_REAL8, MPI_SUM, & @@ -586,7 +586,7 @@ contains ! The MPI_IN_PLACE specifier allows the master to copy values into ! a receive buffer without having a temporary variable -#ifdef MPI +#ifdef OPENMC_MPI call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & MPI_SUM, 0, mpi_intracomm, mpi_err) #endif @@ -610,7 +610,7 @@ contains deallocate(dummy_tally % results) else ! Receive buffer not significant at other processors -#ifdef MPI +#ifdef OPENMC_MPI call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, & 0, mpi_intracomm, mpi_err) #endif @@ -642,6 +642,7 @@ contains integer :: n integer :: int_array(3) integer, allocatable :: array(:) + integer(C_INT64_T) :: seed integer(HID_T) :: file_id integer(HID_T) :: cmfd_group integer(HID_T) :: tallies_group @@ -672,6 +673,7 @@ contains ! Read and overwrite random number seed call read_dataset(seed, file_id, "seed") + call openmc_set_seed(seed) ! It is not impossible for a state point to be generated from a CE run but ! to be loaded in to an MG run (or vice versa), check to prevent that. @@ -849,7 +851,7 @@ contains integer(HID_T) :: plist ! property list #else integer :: i -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code type(Bank), allocatable, target :: temp_source(:) #endif @@ -897,7 +899,7 @@ contains dspace, dset, hdf5_err) ! Save source bank sites since the souce_bank array is overwritten below -#ifdef MPI +#ifdef OPENMC_MPI allocate(temp_source(work)) temp_source(:) = source_bank(:) #endif @@ -907,7 +909,7 @@ contains dims(1) = work_index(i+1) - work_index(i) call h5screate_simple_f(1, dims, memspace, hdf5_err) -#ifdef MPI +#ifdef OPENMC_MPI ! Receive source sites from other processes if (i > 0) then call MPI_RECV(source_bank, int(dims(1)), MPI_BANK, i, i, & @@ -933,12 +935,12 @@ contains call h5dclose_f(dset, hdf5_err) ! Restore state of source bank -#ifdef MPI +#ifdef OPENMC_MPI source_bank(:) = temp_source(:) deallocate(temp_source) #endif else -#ifdef MPI +#ifdef OPENMC_MPI call MPI_SEND(source_bank, int(work), MPI_BANK, 0, rank, & mpi_intracomm, mpi_err) #endif diff --git a/src/summary.F90 b/src/summary.F90 index 0f45ef1bb..3aeb42178 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -122,13 +122,11 @@ contains real(8), allocatable :: cell_temperatures(:) integer(HID_T) :: geom_group integer(HID_T) :: cells_group, cell_group - integer(HID_T) :: surfaces_group, surface_group + integer(HID_T) :: surfaces_group integer(HID_T) :: universes_group, univ_group integer(HID_T) :: lattices_group, lattice_group - real(8), allocatable :: coeffs(:) character(:), allocatable :: region_spec type(Cell), pointer :: c - class(Surface), pointer :: s class(Lattice), pointer :: lat ! Use H5LT interface to write number of geometry objects @@ -223,7 +221,7 @@ contains region_spec = trim(region_spec) // " |" case default region_spec = trim(region_spec) // " " // to_str(& - sign(surfaces(abs(k))%obj%id, k)) + sign(surfaces(abs(k))%id(), k)) end select end do call write_dataset(cell_group, "region", adjustl(region_spec)) @@ -241,92 +239,7 @@ contains ! Write information on each surface SURFACE_LOOP: do i = 1, n_surfaces - s => surfaces(i)%obj - surface_group = create_group(surfaces_group, "surface " // & - trim(to_str(s%id))) - - ! Write name for this surface - call write_dataset(surface_group, "name", s%name) - - ! Write surface type - select type (s) - type is (SurfaceXPlane) - call write_dataset(surface_group, "type", "x-plane") - allocate(coeffs(1)) - coeffs(1) = s%x0 - - type is (SurfaceYPlane) - call write_dataset(surface_group, "type", "y-plane") - allocate(coeffs(1)) - coeffs(1) = s%y0 - - type is (SurfaceZPlane) - call write_dataset(surface_group, "type", "z-plane") - allocate(coeffs(1)) - coeffs(1) = s%z0 - - type is (SurfacePlane) - call write_dataset(surface_group, "type", "plane") - allocate(coeffs(4)) - coeffs(:) = [s%A, s%B, s%C, s%D] - - type is (SurfaceXCylinder) - call write_dataset(surface_group, "type", "x-cylinder") - allocate(coeffs(3)) - coeffs(:) = [s%y0, s%z0, s%r] - - type is (SurfaceYCylinder) - call write_dataset(surface_group, "type", "y-cylinder") - allocate(coeffs(3)) - coeffs(:) = [s%x0, s%z0, s%r] - - type is (SurfaceZCylinder) - call write_dataset(surface_group, "type", "z-cylinder") - allocate(coeffs(3)) - coeffs(:) = [s%x0, s%y0, s%r] - - type is (SurfaceSphere) - call write_dataset(surface_group, "type", "sphere") - allocate(coeffs(4)) - coeffs(:) = [s%x0, s%y0, s%z0, s%r] - - type is (SurfaceXCone) - call write_dataset(surface_group, "type", "x-cone") - allocate(coeffs(4)) - coeffs(:) = [s%x0, s%y0, s%z0, s%r2] - - type is (SurfaceYCone) - call write_dataset(surface_group, "type", "y-cone") - allocate(coeffs(4)) - coeffs(:) = [s%x0, s%y0, s%z0, s%r2] - - type is (SurfaceZCone) - call write_dataset(surface_group, "type", "z-cone") - allocate(coeffs(4)) - coeffs(:) = [s%x0, s%y0, s%z0, s%r2] - - type is (SurfaceQuadric) - call write_dataset(surface_group, "type", "quadric") - allocate(coeffs(10)) - coeffs(:) = [s%A, s%B, s%C, s%D, s%E, s%F, s%G, s%H, s%J, s%K] - - end select - call write_dataset(surface_group, "coefficients", coeffs) - deallocate(coeffs) - - ! Write boundary type - select case (s%bc) - case (BC_TRANSMIT) - call write_dataset(surface_group, "boundary_type", "transmission") - case (BC_VACUUM) - call write_dataset(surface_group, "boundary_type", "vacuum") - case (BC_REFLECT) - call write_dataset(surface_group, "boundary_type", "reflective") - case (BC_PERIODIC) - call write_dataset(surface_group, "boundary_type", "periodic") - end select - - call close_group(surface_group) + call surfaces(i) % to_hdf5(surfaces_group) end do SURFACE_LOOP call close_group(surfaces_group) diff --git a/src/surface.cpp b/src/surface.cpp new file mode 100644 index 000000000..9abfd23d2 --- /dev/null +++ b/src/surface.cpp @@ -0,0 +1,1235 @@ +#include "surface.h" + +#include +#include +#include + +#include "error.h" +#include "hdf5_interface.h" +#include "xml_interface.h" + + +namespace openmc { + +//============================================================================== +// Helper functions for reading the "coeffs" node of an XML surface element +//============================================================================== + +int word_count(const std::string &text) +{ + bool in_word = false; + int count {0}; + for (auto c = text.begin(); c != text.end(); c++) { + if (std::isspace(*c)) { + if (in_word) { + in_word = false; + count++; + } + } else { + in_word = true; + } + } + if (in_word) count++; + return count; +} + +void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1) +{ + // Check the given number of coefficients. + std::string coeffs = get_node_value(surf_node, "coeffs"); + int n_words = word_count(coeffs); + if (n_words != 1) { + std::stringstream err_msg; + err_msg << "Surface " << surf_id << " expects 1 coeff but was given " + << n_words; + fatal_error(err_msg); + } + + // Parse the coefficients. + int stat = sscanf(coeffs.c_str(), "%lf", &c1); + if (stat != 1) { + fatal_error("Something went wrong reading surface coeffs"); + } +} + +void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1, double &c2, + double &c3) +{ + // Check the given number of coefficients. + std::string coeffs = get_node_value(surf_node, "coeffs"); + int n_words = word_count(coeffs); + if (n_words != 3) { + std::stringstream err_msg; + err_msg << "Surface " << surf_id << " expects 3 coeffs but was given " + << n_words; + fatal_error(err_msg); + } + + // Parse the coefficients. + int stat = sscanf(coeffs.c_str(), "%lf %lf %lf", &c1, &c2, &c3); + if (stat != 3) { + fatal_error("Something went wrong reading surface coeffs"); + } +} + +void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1, double &c2, + double &c3, double &c4) +{ + // Check the given number of coefficients. + std::string coeffs = get_node_value(surf_node, "coeffs"); + int n_words = word_count(coeffs); + if (n_words != 4) { + std::stringstream err_msg; + err_msg << "Surface " << surf_id << " expects 4 coeffs but was given " + << n_words; + fatal_error(err_msg); + } + + // Parse the coefficients. + int stat = sscanf(coeffs.c_str(), "%lf %lf %lf %lf", &c1, &c2, &c3, &c4); + if (stat != 4) { + fatal_error("Something went wrong reading surface coeffs"); + } +} + +void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1, double &c2, + double &c3, double &c4, double &c5, double &c6, double &c7, + double &c8, double &c9, double &c10) +{ + // Check the given number of coefficients. + std::string coeffs = get_node_value(surf_node, "coeffs"); + int n_words = word_count(coeffs); + if (n_words != 10) { + std::stringstream err_msg; + err_msg << "Surface " << surf_id << " expects 10 coeffs but was given " + << n_words; + fatal_error(err_msg); + } + + // Parse the coefficients. + int stat = sscanf(coeffs.c_str(), "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", + &c1, &c2, &c3, &c4, &c5, &c6, &c7, &c8, &c9, &c10); + if (stat != 10) { + fatal_error("Something went wrong reading surface coeffs"); + } +} + +//============================================================================== +// Surface implementation +//============================================================================== + +Surface::Surface(pugi::xml_node surf_node) +{ + if (check_for_node(surf_node, "id")) { + id = stoi(get_node_value(surf_node, "id")); + } else { + fatal_error("Must specify id of surface in geometry XML file."); + } + + if (check_for_node(surf_node, "name")) { + name = get_node_value(surf_node, "name"); + } + + if (check_for_node(surf_node, "boundary")) { + std::string surf_bc = get_node_value(surf_node, "boundary"); + + if (surf_bc == "transmission" || surf_bc == "transmit" ||surf_bc.empty()) { + bc = BC_TRANSMIT; + + } else if (surf_bc == "vacuum") { + bc = BC_VACUUM; + + } else if (surf_bc == "reflective" || surf_bc == "reflect" + || surf_bc == "reflecting") { + bc = BC_REFLECT; + } else if (surf_bc == "periodic") { + bc = BC_PERIODIC; + } else { + std::stringstream err_msg; + err_msg << "Unknown boundary condition \"" << surf_bc + << "\" specified on surface " << id; + fatal_error(err_msg); + } + + } else { + bc = BC_TRANSMIT; + } + +} + +bool +Surface::sense(const double xyz[3], const double uvw[3]) const +{ + // 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. + if (std::abs(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); + return uvw[0] * norm[0] + uvw[1] * norm[1] + uvw[2] * norm[2] > 0.0; + } + return f > 0.0; +} + +void +Surface::reflect(const double xyz[3], double uvw[3]) const +{ + // 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]; +} + +void +Surface::to_hdf5(hid_t group_id) const +{ + std::string group_name {"surface "}; + group_name += std::to_string(id); + + hid_t surf_group = create_group(group_id, group_name); + + switch(bc) { + case BC_TRANSMIT : + write_string(surf_group, "boundary_type", "transmission"); + break; + case BC_VACUUM : + write_string(surf_group, "boundary_type", "vacuum"); + break; + case BC_REFLECT : + write_string(surf_group, "boundary_type", "reflective"); + break; + case BC_PERIODIC : + write_string(surf_group, "boundary_type", "periodic"); + break; + } + + if (!name.empty()) { + write_string(surf_group, "name", name); + } + + to_hdf5_inner(surf_group); + + close_group(surf_group); +} + +//============================================================================== +// PeriodicSurface implementation +//============================================================================== + +PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node) + : Surface {surf_node} +{ + if (check_for_node(surf_node, "periodic_surface_id")) { + i_periodic = stoi(get_node_value(surf_node, "periodic_surface_id")); + } +} + +//============================================================================== +// Generic functions for x-, y-, and z-, planes. +//============================================================================== + +// The template parameter indicates the axis normal to the plane. +template double +axis_aligned_plane_evaluate(const double xyz[3], double offset) +{ + return xyz[i] - offset; +} + +// The template parameter indicates the axis normal to the plane. +template double +axis_aligned_plane_distance(const double xyz[3], const double uvw[3], + bool coincident, double offset) +{ + const double f = offset - xyz[i]; + if (coincident or std::abs(f) < FP_COINCIDENT or uvw[i] == 0.0) return INFTY; + const double d = f / uvw[i]; + if (d < 0.0) return INFTY; + return d; +} + +// The first template parameter indicates the axis normal to the plane. The +// other two parameters indicate the other two axes. +template void +axis_aligned_plane_normal(const double xyz[3], double uvw[3]) +{ + uvw[i1] = 1.0; + uvw[i2] = 0.0; + uvw[i3] = 0.0; +} + +//============================================================================== +// SurfaceXPlane implementation +//============================================================================== + +SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node) + : PeriodicSurface(surf_node) +{ + read_coeffs(surf_node, id, x0); +} + +inline double SurfaceXPlane::evaluate(const double xyz[3]) const +{ + return axis_aligned_plane_evaluate<0>(xyz, x0); +} + +inline double SurfaceXPlane::distance(const double xyz[3], const double uvw[3], + bool coincident) const +{ + return axis_aligned_plane_distance<0>(xyz, uvw, coincident, x0); +} + +inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const +{ + axis_aligned_plane_normal<0, 1, 2>(xyz, uvw); +} + +void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "x-plane"); + std::array coeffs {{x0}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +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; + + double u = uvw[0]; + uvw[0] = -uvw[1]; + uvw[1] = u; + + return true; + } +} + +BoundingBox +SurfaceXPlane::bounding_box() const +{ + BoundingBox out {x0, x0, -INFTY, INFTY, -INFTY, INFTY}; + return out; +} + +//============================================================================== +// SurfaceYPlane implementation +//============================================================================== + +SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) + : PeriodicSurface(surf_node) +{ + read_coeffs(surf_node, id, y0); +} + +inline double SurfaceYPlane::evaluate(const double xyz[3]) const +{ + return axis_aligned_plane_evaluate<1>(xyz, y0); +} + +inline double SurfaceYPlane::distance(const double xyz[3], const double uvw[3], + bool coincident) const +{ + return axis_aligned_plane_distance<1>(xyz, uvw, coincident, y0); +} + +inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const +{ + axis_aligned_plane_normal<1, 0, 2>(xyz, uvw); +} + +void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "y-plane"); + std::array coeffs {{y0}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +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; + } +} + +BoundingBox +SurfaceYPlane::bounding_box() const +{ + BoundingBox out {-INFTY, INFTY, y0, y0, -INFTY, INFTY}; + return out; +} + +//============================================================================== +// SurfaceZPlane implementation +//============================================================================== + +SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) + : PeriodicSurface(surf_node) +{ + read_coeffs(surf_node, id, z0); +} + +inline double SurfaceZPlane::evaluate(const double xyz[3]) const +{ + return axis_aligned_plane_evaluate<2>(xyz, z0); +} + +inline double SurfaceZPlane::distance(const double xyz[3], const double uvw[3], + bool coincident) const +{ + return axis_aligned_plane_distance<2>(xyz, uvw, coincident, z0); +} + +inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const +{ + axis_aligned_plane_normal<2, 0, 1>(xyz, uvw); +} + +void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "z-plane"); + std::array coeffs {{z0}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +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; +} + +BoundingBox +SurfaceZPlane::bounding_box() const +{ + BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, z0, z0}; + return out; +} + +//============================================================================== +// SurfacePlane implementation +//============================================================================== + +SurfacePlane::SurfacePlane(pugi::xml_node surf_node) + : PeriodicSurface(surf_node) +{ + read_coeffs(surf_node, id, A, B, C, D); +} + +double +SurfacePlane::evaluate(const double xyz[3]) const +{ + return A*xyz[0] + B*xyz[1] + C*xyz[2] - D; +} + +double +SurfacePlane::distance(const double xyz[3], const double uvw[3], + bool coincident) const +{ + const double f = A*xyz[0] + B*xyz[1] + C*xyz[2] - D; + const double projection = A*uvw[0] + B*uvw[1] + C*uvw[2]; + if (coincident or std::abs(f) < FP_COINCIDENT or projection == 0.0) { + return INFTY; + } else { + const double d = -f / projection; + if (d < 0.0) return INFTY; + return d; + } +} + +void +SurfacePlane::normal(const double xyz[3], double uvw[3]) const +{ + uvw[0] = A; + uvw[1] = B; + uvw[2] = C; +} + +void SurfacePlane::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "plane"); + std::array coeffs {{A, B, C, D}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +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; +} + +BoundingBox +SurfacePlane::bounding_box() const +{ + BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + return out; +} + +//============================================================================== +// Generic functions for x-, y-, and z-, cylinders +//============================================================================== + +// The template parameters indicate the axes perpendicular to the axis of the +// cylinder. offset1 and offset2 should correspond with i1 and i2, +// respectively. +template double +axis_aligned_cylinder_evaluate(const double xyz[3], double offset1, + double offset2, double radius) +{ + const double xyz1 = xyz[i1] - offset1; + const double xyz2 = xyz[i2] - offset2; + return xyz1*xyz1 + xyz2*xyz2 - radius*radius; +} + +// The first template parameter indicates which axis the cylinder is aligned to. +// The other two parameters indicate the other two axes. offset1 and offset2 +// should correspond with i2 and i3, respectively. +template double +axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3], + bool coincident, double offset1, double offset2, double radius) +{ + const double a = 1.0 - uvw[i1]*uvw[i1]; // u^2 + v^2 + if (a == 0.0) return INFTY; + + const double xyz2 = xyz[i2] - offset1; + const double xyz3 = xyz[i3] - offset2; + const double k = xyz2 * uvw[i2] + xyz3 * uvw[i3]; + const double c = xyz2*xyz2 + xyz3*xyz3 - radius*radius; + const double quad = k*k - a*c; + + if (quad < 0.0) { + // No intersection with cylinder. + return INFTY; + + } else if (coincident or std::abs(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). + const double d = (-k - sqrt(quad)) / a; + if (d < 0.0) return INFTY; + return d; + } +} + +// The first template parameter indicates which axis the cylinder is aligned to. +// The other two parameters indicate the other two axes. offset1 and offset2 +// should correspond with i2 and i3, respectively. +template void +axis_aligned_cylinder_normal(const double xyz[3], double uvw[3], double offset1, + double offset2) +{ + uvw[i2] = 2.0 * (xyz[i2] - offset1); + uvw[i3] = 2.0 * (xyz[i3] - offset2); + uvw[i1] = 0.0; +} + +//============================================================================== +// SurfaceXCylinder implementation +//============================================================================== + +SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node) + : Surface(surf_node) +{ + read_coeffs(surf_node, id, y0, z0, r); +} + +inline double SurfaceXCylinder::evaluate(const double xyz[3]) const +{ + return axis_aligned_cylinder_evaluate<1, 2>(xyz, y0, z0, r); +} + +inline double SurfaceXCylinder::distance(const double xyz[3], + const double uvw[3], bool coincident) const +{ + return axis_aligned_cylinder_distance<0, 1, 2>(xyz, uvw, coincident, y0, z0, + r); +} + +inline void SurfaceXCylinder::normal(const double xyz[3], double uvw[3]) const +{ + axis_aligned_cylinder_normal<0, 1, 2>(xyz, uvw, y0, z0); +} + + +void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "x-cylinder"); + std::array coeffs {{y0, z0, r}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +//============================================================================== +// SurfaceYCylinder implementation +//============================================================================== + +SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node) + : Surface(surf_node) +{ + read_coeffs(surf_node, id, x0, z0, r); +} + +inline double SurfaceYCylinder::evaluate(const double xyz[3]) const +{ + return axis_aligned_cylinder_evaluate<0, 2>(xyz, x0, z0, r); +} + +inline double SurfaceYCylinder::distance(const double xyz[3], + const double uvw[3], bool coincident) const +{ + return axis_aligned_cylinder_distance<1, 0, 2>(xyz, uvw, coincident, x0, z0, + r); +} + +inline void SurfaceYCylinder::normal(const double xyz[3], double uvw[3]) const +{ + axis_aligned_cylinder_normal<1, 0, 2>(xyz, uvw, x0, z0); +} + +void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "y-cylinder"); + std::array coeffs {{x0, z0, r}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +//============================================================================== +// SurfaceZCylinder implementation +//============================================================================== + +SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) + : Surface(surf_node) +{ + read_coeffs(surf_node, id, x0, y0, r); +} + +inline double SurfaceZCylinder::evaluate(const double xyz[3]) const +{ + return axis_aligned_cylinder_evaluate<0, 1>(xyz, x0, y0, r); +} + +inline double SurfaceZCylinder::distance(const double xyz[3], + const double uvw[3], bool coincident) const +{ + return axis_aligned_cylinder_distance<2, 0, 1>(xyz, uvw, coincident, x0, y0, + r); +} + +inline void SurfaceZCylinder::normal(const double xyz[3], double uvw[3]) const +{ + axis_aligned_cylinder_normal<2, 0, 1>(xyz, uvw, x0, y0); +} + +void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "z-cylinder"); + std::array coeffs {{x0, y0, r}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +//============================================================================== +// SurfaceSphere implementation +//============================================================================== + +SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node) + : Surface(surf_node) +{ + read_coeffs(surf_node, id, x0, y0, z0, r); +} + +double SurfaceSphere::evaluate(const double xyz[3]) const +{ + const double x = xyz[0] - x0; + const double y = xyz[1] - y0; + const double z = xyz[2] - z0; + return x*x + y*y + z*z - r*r; +} + +double SurfaceSphere::distance(const double xyz[3], const double uvw[3], + bool coincident) const +{ + const double x = xyz[0] - x0; + const double y = xyz[1] - y0; + const double z = xyz[2] - z0; + const double k = x*uvw[0] + y*uvw[1] + z*uvw[2]; + const double c = x*x + y*y + z*z - r*r; + const double quad = k*k - c; + + if (quad < 0.0) { + // No intersection with sphere. + return INFTY; + + } else if (coincident or std::abs(c) < FP_COINCIDENT) { + // Particle is on the sphere, 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); + } + + } else if (c < 0.0) { + // Particle is inside the sphere, 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); + + } else { + // Particle is outside the sphere, thus both distances are either positive + // or negative. If positive, the smaller distance is the one with positive + // sign on sqrt(quad). + const double d = -k - sqrt(quad); + if (d < 0.0) return INFTY; + return d; + } +} + +inline void SurfaceSphere::normal(const double xyz[3], double uvw[3]) const +{ + uvw[0] = 2.0 * (xyz[0] - x0); + uvw[1] = 2.0 * (xyz[1] - y0); + uvw[2] = 2.0 * (xyz[2] - z0); +} + +void SurfaceSphere::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "sphere"); + std::array coeffs {{x0, y0, z0, r}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +//============================================================================== +// Generic functions for x-, y-, and z-, cones +//============================================================================== + +// The first template parameter indicates which axis the cone is aligned to. +// The other two parameters indicate the other two axes. offset1, offset2, +// and offset3 should correspond with i1, i2, and i3, respectively. +template double +axis_aligned_cone_evaluate(const double xyz[3], double offset1, + double offset2, double offset3, double radius_sq) +{ + const double xyz1 = xyz[i1] - offset1; + const double xyz2 = xyz[i2] - offset2; + const double xyz3 = xyz[i3] - offset3; + return xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1; +} + +// The first template parameter indicates which axis the cone is aligned to. +// The other two parameters indicate the other two axes. offset1, offset2, +// and offset3 should correspond with i1, i2, and i3, respectively. +template double +axis_aligned_cone_distance(const double xyz[3], const double uvw[3], + bool coincident, double offset1, double offset2, double offset3, + double radius_sq) +{ + const double xyz1 = xyz[i1] - offset1; + const double xyz2 = xyz[i2] - offset2; + const double xyz3 = xyz[i3] - offset3; + const double a = uvw[i2]*uvw[i2] + uvw[i3]*uvw[i3] + - radius_sq*uvw[i1]*uvw[i1]; + const double k = xyz2*uvw[i2] + xyz3*uvw[i3] - radius_sq*xyz1*uvw[i1]; + const double c = xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1; + double quad = k*k - a*c; + + double d; + + if (quad < 0.0) { + // No intersection with cone. + return INFTY; + + } else if (coincident or std::abs(c) < FP_COINCIDENT) { + // Particle is on the cone, 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) { + d = (-k - sqrt(quad)) / a; + } else { + d = (-k + sqrt(quad)) / a; + } + + } else { + // Calculate both solutions to the quadratic. + quad = sqrt(quad); + d = (-k - quad) / a; + const double b = (-k + quad) / a; + + // Determine the smallest positive solution. + if (d < 0.0) { + if (b > 0.0) d = b; + } else { + if (b > 0.0) { + if (b < d) d = b; + } + } + } + + // If the distance was negative, set boundary distance to infinity. + if (d <= 0.0) return INFTY; + return d; +} + +// The first template parameter indicates which axis the cone is aligned to. +// The other two parameters indicate the other two axes. offset1, offset2, +// and offset3 should correspond with i1, i2, and i3, respectively. +template void +axis_aligned_cone_normal(const double xyz[3], double uvw[3], double offset1, + double offset2, double offset3, double radius_sq) +{ + uvw[i1] = -2.0 * radius_sq * (xyz[i1] - offset1); + uvw[i2] = 2.0 * (xyz[i2] - offset2); + uvw[i3] = 2.0 * (xyz[i3] - offset3); +} + +//============================================================================== +// SurfaceXCone implementation +//============================================================================== + +SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) + : Surface(surf_node) +{ + read_coeffs(surf_node, id, x0, y0, z0, r_sq); +} + +inline double SurfaceXCone::evaluate(const double xyz[3]) const +{ + return axis_aligned_cone_evaluate<0, 1, 2>(xyz, x0, y0, z0, r_sq); +} + +inline double SurfaceXCone::distance(const double xyz[3], + const double uvw[3], bool coincident) const +{ + return axis_aligned_cone_distance<0, 1, 2>(xyz, uvw, coincident, x0, y0, z0, + r_sq); +} + +inline void SurfaceXCone::normal(const double xyz[3], double uvw[3]) const +{ + axis_aligned_cone_normal<0, 1, 2>(xyz, uvw, x0, y0, z0, r_sq); +} + +void SurfaceXCone::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "x-cone"); + std::array coeffs {{x0, y0, z0, r_sq}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +//============================================================================== +// SurfaceYCone implementation +//============================================================================== + +SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) + : Surface(surf_node) +{ + read_coeffs(surf_node, id, x0, y0, z0, r_sq); +} + +inline double SurfaceYCone::evaluate(const double xyz[3]) const +{ + return axis_aligned_cone_evaluate<1, 0, 2>(xyz, y0, x0, z0, r_sq); +} + +inline double SurfaceYCone::distance(const double xyz[3], + const double uvw[3], bool coincident) const +{ + return axis_aligned_cone_distance<1, 0, 2>(xyz, uvw, coincident, y0, x0, z0, + r_sq); +} + +inline void SurfaceYCone::normal(const double xyz[3], double uvw[3]) const +{ + axis_aligned_cone_normal<1, 0, 2>(xyz, uvw, y0, x0, z0, r_sq); +} + +void SurfaceYCone::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "y-cone"); + std::array coeffs {{x0, y0, z0, r_sq}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +//============================================================================== +// SurfaceZCone implementation +//============================================================================== + +SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) + : Surface(surf_node) +{ + read_coeffs(surf_node, id, x0, y0, z0, r_sq); +} + +inline double SurfaceZCone::evaluate(const double xyz[3]) const +{ + return axis_aligned_cone_evaluate<2, 0, 1>(xyz, z0, x0, y0, r_sq); +} + +inline double SurfaceZCone::distance(const double xyz[3], + const double uvw[3], bool coincident) const +{ + return axis_aligned_cone_distance<2, 0, 1>(xyz, uvw, coincident, z0, x0, y0, + r_sq); +} + +inline void SurfaceZCone::normal(const double xyz[3], double uvw[3]) const +{ + axis_aligned_cone_normal<2, 0, 1>(xyz, uvw, z0, x0, y0, r_sq); +} + +void SurfaceZCone::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "z-cone"); + std::array coeffs {{x0, y0, z0, r_sq}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +//============================================================================== +// SurfaceQuadric implementation +//============================================================================== + +SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node) + : Surface(surf_node) +{ + read_coeffs(surf_node, id, A, B, C, D, E, F, G, H, J, K); +} + +double +SurfaceQuadric::evaluate(const double xyz[3]) const +{ + const double &x = xyz[0]; + const double &y = xyz[1]; + const double &z = xyz[2]; + return x*(A*x + D*y + G) + + y*(B*y + E*z + H) + + z*(C*z + F*x + J) + K; +} + +double +SurfaceQuadric::distance(const double xyz[3], + const double uvw[3], bool coincident) const +{ + const double &x = xyz[0]; + const double &y = xyz[1]; + const double &z = xyz[2]; + const double &u = uvw[0]; + const double &v = uvw[1]; + const double &w = uvw[2]; + + const double a = A*u*u + B*v*v + C*w*w + D*u*v + E*v*w + F*u*w; + const double k = (A*u*x + B*v*y + C*w*z + 0.5*(D*(u*y + v*x) + + E*(v*z + w*y) + F*(w*x + u*z) + G*u + H*v + J*w)); + const double c = A*x*x + B*y*y + C*z*z + D*x*y + E*y*z + F*x*z + G*x + H*y + + J*z + K; + double quad = k*k - a*c; + + double d; + + if (quad < 0.0) { + // No intersection with surface. + return INFTY; + + } else if (coincident or std::abs(c) < FP_COINCIDENT) { + // Particle is on the surface, thus one distance is positive/negative and + // the other is zero. The sign of k determines which distance is zero and + // which is not. + if (k >= 0.0) { + d = (-k - sqrt(quad)) / a; + } else { + d = (-k + sqrt(quad)) / a; + } + + } else { + // Calculate both solutions to the quadratic. + quad = sqrt(quad); + d = (-k - quad) / a; + double b = (-k + quad) / a; + + // Determine the smallest positive solution. + if (d < 0.0) { + if (b > 0.0) d = b; + } else { + if (b > 0.0) { + if (b < d) d = b; + } + } + } + + // If the distance was negative, set boundary distance to infinity. + if (d <= 0.0) return INFTY; + return d; +} + +void +SurfaceQuadric::normal(const double xyz[3], double uvw[3]) const +{ + const double &x = xyz[0]; + const double &y = xyz[1]; + const double &z = xyz[2]; + uvw[0] = 2.0*A*x + D*y + F*z + G; + uvw[1] = 2.0*B*y + D*x + E*z + H; + uvw[2] = 2.0*C*z + E*y + F*x + J; +} + +void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "quadric"); + std::array coeffs {{A, B, C, D, E, F, G, H, J, K}}; + write_double_1D(group_id, "coefficients", coeffs); +} + +//============================================================================== + +extern "C" void +read_surfaces(pugi::xml_node *node) +{ + // Count the number of surfaces. + for (pugi::xml_node surf_node: node->children("surface")) {n_surfaces++;} + if (n_surfaces == 0) { + fatal_error("No surfaces found in geometry.xml!"); + } + + // 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++) { + std::string surf_type = get_node_value(surf_node, "type"); + + if (surf_type == "x-plane") { + surfaces_c[i_surf] = new SurfaceXPlane(surf_node); + + } else if (surf_type == "y-plane") { + surfaces_c[i_surf] = new SurfaceYPlane(surf_node); + + } else if (surf_type == "z-plane") { + surfaces_c[i_surf] = new SurfaceZPlane(surf_node); + + } else if (surf_type == "plane") { + surfaces_c[i_surf] = new SurfacePlane(surf_node); + + } else if (surf_type == "x-cylinder") { + surfaces_c[i_surf] = new SurfaceXCylinder(surf_node); + + } else if (surf_type == "y-cylinder") { + surfaces_c[i_surf] = new SurfaceYCylinder(surf_node); + + } else if (surf_type == "z-cylinder") { + surfaces_c[i_surf] = new SurfaceZCylinder(surf_node); + + } else if (surf_type == "sphere") { + surfaces_c[i_surf] = new SurfaceSphere(surf_node); + + } else if (surf_type == "x-cone") { + surfaces_c[i_surf] = new SurfaceXCone(surf_node); + + } else if (surf_type == "y-cone") { + surfaces_c[i_surf] = new SurfaceYCone(surf_node); + + } else if (surf_type == "z-cone") { + surfaces_c[i_surf] = new SurfaceZCone(surf_node); + + } else if (surf_type == "quadric") { + surfaces_c[i_surf] = new SurfaceQuadric(surf_node); + + } else { + std::stringstream err_msg; + err_msg << "Invalid surface type, \"" << surf_type << "\""; + fatal_error(err_msg); + } + } + } + + // Fill the surface dictionary. + for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { + int id = surfaces_c[i_surf]->id; + auto in_dict = surface_dict.find(id); + if (in_dict == surface_dict.end()) { + surface_dict[id] = i_surf; + } else { + std::stringstream err_msg; + err_msg << "Two or more surfaces use the same unique ID: " << id; + fatal_error(err_msg); + } + } + + // Find the global bounding box (of periodic BC surfaces). + double xmin {INFTY}, xmax {-INFTY}, ymin {INFTY}, ymax {-INFTY}, + zmin {INFTY}, zmax {-INFTY}; + int i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax; + for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { + if (surfaces_c[i_surf]->bc == BC_PERIODIC) { + // Downcast to the PeriodicSurface type. + Surface *surf_base = surfaces_c[i_surf]; + PeriodicSurface *surf = dynamic_cast(surf_base); + + // Make sure this surface inherits from PeriodicSurface. + if (!surf) { + std::stringstream err_msg; + err_msg << "Periodic boundary condition not supported for surface " + << surf_base->id + << ". Periodic BCs are only supported for planar surfaces."; + fatal_error(err_msg); + } + + // See if this surface makes part of the global bounding box. + BoundingBox bb = surf->bounding_box(); + if (bb.xmin > -INFTY and bb.xmin < xmin) { + xmin = bb.xmin; + i_xmin = i_surf; + } + if (bb.xmax < INFTY and bb.xmax > xmax) { + xmax = bb.xmax; + i_xmax = i_surf; + } + if (bb.ymin > -INFTY and bb.ymin < ymin) { + ymin = bb.ymin; + i_ymin = i_surf; + } + if (bb.ymax < INFTY and bb.ymax > ymax) { + ymax = bb.ymax; + i_ymax = i_surf; + } + if (bb.zmin > -INFTY and bb.zmin < zmin) { + zmin = bb.zmin; + i_zmin = i_surf; + } + if (bb.zmax < INFTY and bb.zmax > zmax) { + zmax = bb.zmax; + i_zmax = i_surf; + } + } + } + + // Set i_periodic for periodic BC surfaces. + for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { + if (surfaces_c[i_surf]->bc == BC_PERIODIC) { + // Downcast to the PeriodicSurface type. + Surface *surf_base = surfaces_c[i_surf]; + PeriodicSurface *surf = dynamic_cast(surf_base); + + // Also try downcasting to the SurfacePlane type (which must be handled + // differently). + SurfacePlane *surf_p = dynamic_cast(surf); + + if (!surf_p) { + // This is not a SurfacePlane. + if (surf->i_periodic == C_NONE) { + // The user did not specify the matching periodic surface. See if we + // can find the partnered surface from the bounding box information. + if (i_surf == i_xmin) { + surf->i_periodic = i_xmax; + } else if (i_surf == i_xmax) { + surf->i_periodic = i_xmin; + } else if (i_surf == i_ymin) { + surf->i_periodic = i_ymax; + } else if (i_surf == i_ymax) { + surf->i_periodic = i_ymin; + } else if (i_surf == i_zmin) { + surf->i_periodic = i_zmax; + } else if (i_surf == i_zmax) { + surf->i_periodic = i_zmin; + } else { + fatal_error("Periodic boundary condition applied to interior " + "surface"); + } + } else { + // Convert the surface id to an index. + surf->i_periodic = surface_dict[surf->i_periodic]; + } + } else { + // This is a SurfacePlane. We won't try to find it's partner if the + // user didn't specify one. + if (surf->i_periodic == C_NONE) { + std::stringstream err_msg; + err_msg << "No matching periodic surface specified for periodic " + "boundary condition on surface " << surf->id; + fatal_error(err_msg); + } else { + // Convert the surface id to an index. + surf->i_periodic = surface_dict[surf->i_periodic]; + } + } + + // Make sure the opposite surface is also periodic. + if (surfaces_c[surf->i_periodic]->bc != BC_PERIODIC) { + std::stringstream err_msg; + err_msg << "Could not find matching surface for periodic boundary " + "condition on surface " << surf->id; + fatal_error(err_msg); + } + } + } +} + +} // namespace openmc diff --git a/src/surface.h b/src/surface.h new file mode 100644 index 000000000..627a8fbae --- /dev/null +++ b/src/surface.h @@ -0,0 +1,427 @@ +#ifndef SURFACE_H +#define SURFACE_H + +#include +#include // For numeric_limits +#include + +#include "hdf5.h" +#include "pugixml/pugixml.hpp" + + +namespace openmc { + +//============================================================================== +// Module constants +//============================================================================== + +extern "C" const int BC_TRANSMIT {0}; +extern "C" const int BC_VACUUM {1}; +extern "C" const int BC_REFLECT {2}; +extern "C" const int BC_PERIODIC {3}; + +//============================================================================== +// Constants that should eventually be moved out of this file +//============================================================================== + +extern "C" double FP_COINCIDENT; +constexpr double INFTY{std::numeric_limits::max()}; +constexpr int C_NONE {-1}; + +//============================================================================== +// Global variables +//============================================================================== + +// Braces force n_surfaces to be defined here, not just declared. +extern "C" {int32_t n_surfaces {0};} + +class Surface; +Surface **surfaces_c; + +std::map surface_dict; + +//============================================================================== +//! Coordinates for an axis-aligned cube that bounds a geometric object. +//============================================================================== + +struct BoundingBox +{ + double xmin; + double xmax; + double ymin; + double ymax; + double zmin; + double zmax; +}; + +//============================================================================== +//! A geometry primitive used to define regions of 3D space. +//============================================================================== + +class Surface +{ +public: + 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 + std::string name{""}; //!< User-defined name + + explicit Surface(pugi::xml_node surf_node); + + virtual ~Surface() {} + + //! 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 + //! false otherwise. + bool sense(const double xyz[3], const double uvw[3]) const; + + //! Determine the direction of a ray reflected from the surface. + //! @param xyz[3] The point at which the ray is incident. + //! @param uvw[3] A direction. This is both an input and an output parameter. + //! It specifies the icident direction on input and the reflected direction + //! on output. + void reflect(const double xyz[3], double uvw[3]) const; + + //! Evaluate the equation describing the surface. + //! + //! Surfaces can be described by some function f(x, y, z) = 0. This member + //! function evaluates that mathematical function. + //! @param xyz[3] A 3D Cartesian coordinate. + virtual double evaluate(const double xyz[3]) const = 0; + + //! Compute the distance between a point and the surface along a ray. + //! @param xyz[3] A 3D Cartesian coordinate. + //! @param uvw[3] The direction of the ray. + //! @param coincident A hint to the code that the given point should lie + //! exactly on the surface. + virtual double distance(const double xyz[3], const double uvw[3], + bool coincident) const = 0; + + //! Compute the local outward normal direction of the surface. + //! @param xyz[3] A 3D Cartesian coordinate. + //! @param uvw[3] This output argument provides the normal. + virtual void normal(const double xyz[3], double uvw[3]) const = 0; + + //! Write all information needed to reconstruct the surface to an HDF5 group. + //! @param group_id An HDF5 group id. + void to_hdf5(hid_t group_id) const; + +protected: + virtual void to_hdf5_inner(hid_t group_id) const = 0; +}; + +//============================================================================== +//! 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: + int i_periodic{C_NONE}; //!< Index of corresponding periodic surface + + explicit PeriodicSurface(pugi::xml_node surf_node); + + //! 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; + + //! Get the bounding box for this surface. + virtual BoundingBox bounding_box() const = 0; +}; + +//============================================================================== +//! A plane perpendicular to the x-axis. +// +//! The plane is described by the equation \f$x - x_0 = 0\f$ +//============================================================================== + +class SurfaceXPlane : public PeriodicSurface +{ + double x0; +public: + explicit SurfaceXPlane(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], bool coincident) + const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; + bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + const; + BoundingBox bounding_box() const; +}; + +//============================================================================== +//! A plane perpendicular to the y-axis. +// +//! The plane is described by the equation \f$y - y_0 = 0\f$ +//============================================================================== + +class SurfaceYPlane : public PeriodicSurface +{ + double y0; +public: + explicit SurfaceYPlane(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; + bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + const; + BoundingBox bounding_box() const; +}; + +//============================================================================== +//! A plane perpendicular to the z-axis. +// +//! The plane is described by the equation \f$z - z_0 = 0\f$ +//============================================================================== + +class SurfaceZPlane : public PeriodicSurface +{ + double z0; +public: + explicit SurfaceZPlane(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; + bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + const; + BoundingBox bounding_box() const; +}; + +//============================================================================== +//! A general plane. +// +//! The plane is described by the equation \f$A x + B y + C z - D = 0\f$ +//============================================================================== + +class SurfacePlane : public PeriodicSurface +{ + double A, B, C, D; +public: + explicit SurfacePlane(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], bool coincident) + const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; + bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + const; + BoundingBox bounding_box() const; +}; + +//============================================================================== +//! A cylinder aligned along the x-axis. +// +//! The cylinder is described by the equation +//! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$ +//============================================================================== + +class SurfaceXCylinder : public Surface +{ + double y0, z0, r; +public: + explicit SurfaceXCylinder(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; +}; + +//============================================================================== +//! A cylinder aligned along the y-axis. +// +//! The cylinder is described by the equation +//! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 = 0\f$ +//============================================================================== + +class SurfaceYCylinder : public Surface +{ + double x0, z0, r; +public: + explicit SurfaceYCylinder(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; +}; + +//============================================================================== +//! A cylinder aligned along the z-axis. +// +//! The cylinder is described by the equation +//! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 = 0\f$ +//============================================================================== + +class SurfaceZCylinder : public Surface +{ + double x0, y0, r; +public: + explicit SurfaceZCylinder(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; +}; + +//============================================================================== +//! A sphere. +// +//! The cylinder is described by the equation +//! \f$(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$ +//============================================================================== + +class SurfaceSphere : public Surface +{ + double x0, y0, z0, r; +public: + explicit SurfaceSphere(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; +}; + +//============================================================================== +//! A cone aligned along the x-axis. +// +//! The cylinder is described by the equation +//! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 (x - x_0)^2 = 0\f$ +//============================================================================== + +class SurfaceXCone : public Surface +{ + double x0, y0, z0, r_sq; +public: + explicit SurfaceXCone(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; +}; + +//============================================================================== +//! A cone aligned along the y-axis. +// +//! The cylinder is described by the equation +//! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 (y - y_0)^2 = 0\f$ +//============================================================================== + +class SurfaceYCone : public Surface +{ + double x0, y0, z0, r_sq; +public: + explicit SurfaceYCone(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; +}; + +//============================================================================== +//! A cone aligned along the z-axis. +// +//! The cylinder is described by the equation +//! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 (z - z_0)^2 = 0\f$ +//============================================================================== + +class SurfaceZCone : public Surface +{ + double x0, y0, z0, r_sq; +public: + explicit SurfaceZCone(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; +}; + +//============================================================================== +//! A general surface described by a quadratic equation. +// +//! \f$A x^2 + B y^2 + C z^2 + D x y + E y z + F x z + G x + H y + J z + K = 0\f$ +//============================================================================== + +class SurfaceQuadric : public Surface +{ + // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 + double A, B, C, D, E, F, G, H, J, K; +public: + explicit SurfaceQuadric(pugi::xml_node surf_node); + double evaluate(const double xyz[3]) const; + double distance(const double xyz[3], const double uvw[3], + bool coincident) const; + void normal(const double xyz[3], double uvw[3]) const; + void to_hdf5_inner(hid_t group_id) const; +}; + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];} + +extern "C" int surface_id(Surface *surf) {return surf->id;} + +extern "C" int surface_bc(Surface *surf) {return surf->bc;} + +extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3]) +{return surf->sense(xyz, uvw);} + +extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3]) +{surf->reflect(xyz, uvw);} + +extern "C" double +surface_distance(Surface *surf, double xyz[3], double uvw[3], bool coincident) +{return surf->distance(xyz, uvw, coincident);} + +extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3]) +{return surf->normal(xyz, uvw);} + +extern "C" void surface_to_hdf5(Surface *surf, hid_t group) +{surf->to_hdf5(group);} + +extern "C" int surface_i_periodic(PeriodicSurface *surf) +{return surf->i_periodic;} + +extern "C" bool +surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3], + double uvw[3]) +{return surf->periodic_translate(other, xyz, uvw);} + +extern "C" void free_memory_surfaces_c() +{ + for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];} + delete surfaces_c; + surfaces_c = nullptr; + n_surfaces = 0; + surface_dict.clear(); +} + +} // namespace openmc +#endif // SURFACE_H diff --git a/src/surface_header.F90 b/src/surface_header.F90 index ed4ef86e0..9ed89a14e 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -1,1073 +1,212 @@ module surface_header use, intrinsic :: ISO_C_BINDING + use hdf5 - use constants, only: NONE, ONE, TWO, ZERO, HALF, INFINITY, FP_COINCIDENT use dict_header, only: DictIntInt implicit none + interface + pure function surface_pointer_c(surf_ind) & + bind(C, name='surface_pointer') result(ptr) + use ISO_C_BINDING + implicit none + integer(C_INT), intent(in), value :: surf_ind + type(C_PTR) :: ptr + end function surface_pointer_c + + pure function surface_id_c(surf_ptr) bind(C, name='surface_id') result(id) + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: surf_ptr + integer(C_INT) :: id + end function surface_id_c + + pure function surface_bc_c(surf_ptr) bind(C, name='surface_bc') result(bc) + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: surf_ptr + integer(C_INT) :: bc + end function surface_bc_c + + pure function surface_sense_c(surf_ptr, xyz, uvw) & + bind(C, name='surface_sense') result(sense) + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: surf_ptr + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + logical(C_BOOL) :: sense + end function surface_sense_c + + pure subroutine surface_reflect_c(surf_ptr, xyz, uvw) & + bind(C, name='surface_reflect') + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: surf_ptr + 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_ptr, xyz, uvw, coincident) & + bind(C, name='surface_distance') result(d) + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: surf_ptr + real(C_DOUBLE), intent(in) :: xyz(3); + real(C_DOUBLE), intent(in) :: uvw(3); + logical(C_BOOL), intent(in), value :: coincident; + real(C_DOUBLE) :: d; + end function surface_distance_c + + pure subroutine surface_normal_c(surf_ptr, xyz, uvw) & + bind(C, name='surface_normal') + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: surf_ptr + real(C_DOUBLE), intent(in) :: xyz(3); + real(C_DOUBLE), intent(out) :: uvw(3); + end subroutine surface_normal_c + + subroutine surface_to_hdf5_c(surf_ptr, group) & + bind(C, name='surface_to_hdf5') + use ISO_C_BINDING + use hdf5 + implicit none + type(C_PTR), intent(in), value :: surf_ptr + integer(HID_T), intent(in), value :: group + end subroutine surface_to_hdf5_c + + pure function surface_i_periodic_c(surf_ptr) & + bind(C, name="surface_i_periodic") result(i_periodic) + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: surf_ptr + integer(C_INT) :: i_periodic + end function surface_i_periodic_c + + function surface_periodic_c(surf_ptr, other_ptr, xyz, uvw) & + bind(C, name="surface_periodic") result(rotational) + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: surf_ptr + type(C_PTR), intent(in), value :: other_ptr + real(C_DOUBLE), intent(inout) :: xyz(3); + real(C_DOUBLE), intent(inout) :: uvw(3); + logical(C_BOOL) :: rotational + end function surface_periodic_c + + subroutine free_memory_surfaces_c() bind(C) + end subroutine free_memory_surfaces_c + end interface + !=============================================================================== ! SURFACE type defines a first- or second-order surface that can be used to ! construct closed volumes (cells) !=============================================================================== - type, abstract :: Surface - integer :: id ! Unique ID + type :: Surface integer, allocatable :: & neighbor_pos(:), & ! List of cells on positive side neighbor_neg(:) ! List of cells on negative side - integer :: bc ! Boundary condition - integer :: i_periodic = NONE ! Index of corresponding periodic surface - character(len=104) :: name = "" ! User-defined name + type(C_PTR) :: ptr + contains - procedure :: sense - procedure :: reflect - procedure(surface_evaluate_), deferred :: evaluate - procedure(surface_distance_), deferred :: distance - procedure(surface_normal_), deferred :: normal + + procedure :: id => surface_id + procedure :: bc => surface_bc + procedure :: sense => surface_sense + procedure :: reflect => surface_reflect + procedure :: distance => surface_distance + procedure :: normal => surface_normal + procedure :: to_hdf5 => surface_to_hdf5 + procedure :: i_periodic => surface_i_periodic + procedure :: periodic_translate => surface_periodic + 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_distance_(this, xyz, uvw, coincident) result(d) - import Surface - class(Surface), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - end function surface_distance_ - - 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 -!=============================================================================== - - type :: SurfaceContainer - class(Surface), allocatable :: obj - end type SurfaceContainer - -!=============================================================================== -! All the derived types below are extensions of the abstract Surface type. They -! inherent the reflect() and sense() type-bound procedures and must implement -! evaluate(), distance(), and normal() -!=============================================================================== - - type, extends(Surface) :: SurfaceXPlane - ! x = x0 - real(8) :: x0 - contains - procedure :: evaluate => x_plane_evaluate - procedure :: distance => x_plane_distance - procedure :: normal => x_plane_normal - end type SurfaceXPlane - - type, extends(Surface) :: SurfaceYPlane - ! y = y0 - real(8) :: y0 - contains - procedure :: evaluate => y_plane_evaluate - procedure :: distance => y_plane_distance - procedure :: normal => y_plane_normal - end type SurfaceYPlane - - type, extends(Surface) :: SurfaceZPlane - ! z = z0 - real(8) :: z0 - contains - procedure :: evaluate => z_plane_evaluate - procedure :: distance => z_plane_distance - procedure :: normal => z_plane_normal - end type SurfaceZPlane - - type, extends(Surface) :: SurfacePlane - ! Ax + By + Cz = D - real(8) :: A - real(8) :: B - real(8) :: C - real(8) :: D - contains - procedure :: evaluate => plane_evaluate - procedure :: distance => plane_distance - procedure :: normal => plane_normal - end type SurfacePlane - - type, extends(Surface) :: SurfaceXCylinder - ! (y - y0)^2 + (z - z0)^2 = R^2 - real(8) :: y0 - real(8) :: z0 - real(8) :: r - contains - procedure :: evaluate => x_cylinder_evaluate - procedure :: distance => x_cylinder_distance - procedure :: normal => x_cylinder_normal - end type SurfaceXCylinder - - type, extends(Surface) :: SurfaceYCylinder - ! (x - x0)^2 + (z - z0)^2 = R^2 - real(8) :: x0 - real(8) :: z0 - real(8) :: r - contains - procedure :: evaluate => y_cylinder_evaluate - procedure :: distance => y_cylinder_distance - procedure :: normal => y_cylinder_normal - end type SurfaceYCylinder - - type, extends(Surface) :: SurfaceZCylinder - ! (x - x0)^2 + (y - y0)^2 = R^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: r - contains - procedure :: evaluate => z_cylinder_evaluate - procedure :: distance => z_cylinder_distance - procedure :: normal => z_cylinder_normal - end type SurfaceZCylinder - - type, extends(Surface) :: SurfaceSphere - ! (x - x0)^2 + (y - y0)^2 + (z - z0)^2 = R^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: z0 - real(8) :: r - contains - procedure :: evaluate => sphere_evaluate - procedure :: distance => sphere_distance - procedure :: normal => sphere_normal - end type SurfaceSphere - - type, extends(Surface) :: SurfaceXCone - ! (y - y0)^2 + (z - z0)^2 = R^2*(x - x0)^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: z0 - real(8) :: r2 - contains - procedure :: evaluate => x_cone_evaluate - procedure :: distance => x_cone_distance - procedure :: normal => x_cone_normal - end type SurfaceXCone - - type, extends(Surface) :: SurfaceYCone - ! (x - x0)^2 + (z - z0)^2 = R^2*(y - y0)^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: z0 - real(8) :: r2 - contains - procedure :: evaluate => y_cone_evaluate - procedure :: distance => y_cone_distance - procedure :: normal => y_cone_normal - end type SurfaceYCone - - type, extends(Surface) :: SurfaceZCone - ! (x - x0)^2 + (y - y0)^2 = R^2*(z - z0)^2 - real(8) :: x0 - real(8) :: y0 - real(8) :: z0 - real(8) :: r2 - contains - procedure :: evaluate => z_cone_evaluate - procedure :: distance => z_cone_distance - 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 :: evaluate => quadric_evaluate - procedure :: distance => quadric_distance - procedure :: normal => quadric_normal - end type SurfaceQuadric - integer(C_INT32_T), bind(C) :: n_surfaces ! # of surfaces - type(SurfaceContainer), allocatable, target :: surfaces(:) + type(Surface), allocatable, target :: surfaces(:) ! Dictionary that maps user IDs to indices in 'surfaces' type(DictIntInt) :: surface_dict contains -!=============================================================================== -! SENSE determines whether a point is on the 'positive' or 'negative' side of a -! surface. This routine is crucial for determining what cell a particular point -! is in. The positive side is indicated by a returned value of .true. and the -! negative side is indicated by a returned value of .false. -!=============================================================================== - - pure function sense(this, xyz, uvw) result(s) - class(Surface), intent(in) :: this ! surface - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical :: s ! sense of particle - - real(8) :: f ! surface function evaluated at point - - ! Evaluate the surface equation at the particle's coordinates to determine - ! which side the particle is on - f = this%evaluate(xyz) - - ! Check which side of surface the point is on - if (abs(f) < FP_COINCIDENT) then - ! 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. - s = (dot_product(uvw, this%normal(xyz)) > ZERO) - else - s = (f > ZERO) - end if - end function sense - -!=============================================================================== -! 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) + pure function surface_id(this) result(id) 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceXPlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: f - - f = this%x0 - xyz(1) - if (coincident .or. abs(f) < FP_COINCIDENT .or. uvw(1) == ZERO) then - d = INFINITY - else - d = f/uvw(1) - if (d < ZERO) d = INFINITY - end if - end function x_plane_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceYPlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: f - - f = this%y0 - xyz(2) - if (coincident .or. abs(f) < FP_COINCIDENT .or. uvw(2) == ZERO) then - d = INFINITY - else - d = f/uvw(2) - if (d < ZERO) d = INFINITY - end if - end function y_plane_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceZPlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: f - - f = this%z0 - xyz(3) - if (coincident .or. abs(f) < FP_COINCIDENT .or. uvw(3) == ZERO) then - d = INFINITY - else - d = f/uvw(3) - if (d < ZERO) d = INFINITY - end if - end function z_plane_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfacePlane), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: f - real(8) :: projection - - f = this%A*xyz(1) + this%B*xyz(2) + this%C*xyz(3) - this%D - projection = this%A*uvw(1) + this%B*uvw(2) + this%C*uvw(3) - if (coincident .or. abs(f) < FP_COINCIDENT .or. projection == ZERO) then - d = INFINITY - else - d = -f/projection - if (d < ZERO) d = INFINITY - end if - end function plane_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceXCylinder), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: y, z, k, a, c, quad - - a = ONE - uvw(1)*uvw(1) ! v^2 + w^2 - if (a == ZERO) then - d = INFINITY - else - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - k = y*uvw(2) + z*uvw(3) - c = y*y + z*z - this%r*this%r - quad = k*k - a*c - - if (quad < ZERO) then - ! no intersection with cylinder - - d = INFINITY - - elseif (coincident .or. abs(c) < FP_COINCIDENT) then - ! 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 >= ZERO) then - d = INFINITY - else - d = (-k + sqrt(quad))/a - end if - - elseif (c < ZERO) then - ! 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) - - d = (-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) - - d = (-k - sqrt(quad))/a - if (d < ZERO) d = INFINITY - - end if - end if - end function x_cylinder_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceYCylinder), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: x, z, k, a, c, quad - - a = ONE - uvw(2)*uvw(2) ! u^2 + w^2 - if (a == ZERO) then - d = INFINITY - else - x = xyz(1) - this%x0 - z = xyz(3) - this%z0 - k = x*uvw(1) + z*uvw(3) - c = x*x + z*z - this%r*this%r - quad = k*k - a*c - - if (quad < ZERO) then - ! no intersection with cylinder - - d = INFINITY - - elseif (coincident .or. abs(c) < FP_COINCIDENT) then - ! 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 >= ZERO) then - d = INFINITY - else - d = (-k + sqrt(quad))/a - end if - - elseif (c < ZERO) then - ! 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) - - d = (-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) - - d = (-k - sqrt(quad))/a - if (d < ZERO) d = INFINITY - - end if - end if - end function y_cylinder_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceZCylinder), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: x, y, k, a, c, quad - - a = ONE - uvw(3)*uvw(3) ! u^2 + v^2 - if (a == ZERO) then - d = INFINITY - else - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - k = x*uvw(1) + y*uvw(2) - c = x*x + y*y - this%r*this%r - quad = k*k - a*c - - if (quad < ZERO) then - ! no intersection with cylinder - - d = INFINITY - - elseif (coincident .or. abs(c) < FP_COINCIDENT) then - ! 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 >= ZERO) then - d = INFINITY - else - d = (-k + sqrt(quad))/a - end if - - elseif (c < ZERO) then - ! 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) - - d = (-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) - - d = (-k - sqrt(quad))/a - if (d < ZERO) d = INFINITY - - end if - end if - end function z_cylinder_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceSphere), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: x, y, z, k, c, quad - - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - k = x*uvw(1) + y*uvw(2) + z*uvw(3) - c = x*x + y*y + z*z - this%r*this%r - quad = k*k - c - - if (quad < ZERO) then - ! no intersection with sphere - - d = INFINITY - - elseif (coincident .or. abs(c) < FP_COINCIDENT) then - ! particle is on the sphere, 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 >= ZERO) then - d = INFINITY - else - d = -k + sqrt(quad) - end if - - elseif (c < ZERO) then - ! particle is inside the sphere, thus one distance must be negative and - ! one must be positive. The positive distance will be the one with - ! negative sign on sqrt(quad) - - d = -k + sqrt(quad) - - else - ! particle is outside the sphere, thus both distances are either positive - ! or negative. If positive, the smaller distance is the one with positive - ! sign on sqrt(quad) - - d = -k - sqrt(quad) - if (d < ZERO) d = INFINITY - - end if - end function sphere_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceXCone), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: x, y, z, k, a, b, c, quad - - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - a = uvw(2)*uvw(2) + uvw(3)*uvw(3) - this%r2*uvw(1)*uvw(1) - k = y*uvw(2) + z*uvw(3) - this%r2*x*uvw(1) - c = y*y + z*z - this%r2*x*x - quad = k*k - a*c - - if (quad < ZERO) then - ! no intersection with cone - - d = INFINITY - - elseif (coincident .or. abs(c) < FP_COINCIDENT) then - ! particle is on the cone, thus one distance is positive/negative and the - ! other is zero. The sign of k determines which distance is zero and which - ! is not. - - if (k >= ZERO) then - d = (-k - sqrt(quad))/a - else - d = (-k + sqrt(quad))/a - end if - - else - ! calculate both solutions to the quadratic - quad = sqrt(quad) - d = (-k - quad)/a - b = (-k + quad)/a - - ! determine the smallest positive solution - if (d < ZERO) then - if (b > ZERO) then - d = b - end if - else - if (b > ZERO) d = min(d, b) - end if - end if - - ! If the distance was negative, set boundary distance to infinity - if (d <= ZERO) d = INFINITY - end function x_cone_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceYCone), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: x, y, z, k, a, b, c, quad - - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - a = uvw(1)*uvw(1) + uvw(3)*uvw(3) - this%r2*uvw(2)*uvw(2) - k = x*uvw(1) + z*uvw(3) - this%r2*y*uvw(2) - c = x*x + z*z - this%r2*y*y - quad = k*k - a*c - - if (quad < ZERO) then - ! no intersection with cone - - d = INFINITY - - elseif (coincident .or. abs(c) < FP_COINCIDENT) then - ! particle is on the cone, thus one distance is positive/negative and the - ! other is zero. The sign of k determines which distance is zero and which - ! is not. - - if (k >= ZERO) then - d = (-k - sqrt(quad))/a - else - d = (-k + sqrt(quad))/a - end if - - else - ! calculate both solutions to the quadratic - quad = sqrt(quad) - d = (-k - quad)/a - b = (-k + quad)/a - - ! determine the smallest positive solution - if (d < ZERO) then - if (b > ZERO) then - d = b - end if - else - if (b > ZERO) d = min(d, b) - end if - end if - - ! If the distance was negative, set boundary distance to infinity - if (d <= ZERO) d = INFINITY - end function y_cone_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceZCone), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: x, y, z, k, a, b, c, quad - - x = xyz(1) - this%x0 - y = xyz(2) - this%y0 - z = xyz(3) - this%z0 - a = uvw(1)*uvw(1) + uvw(2)*uvw(2) - this%r2*uvw(3)*uvw(3) - k = x*uvw(1) + y*uvw(2) - this%r2*z*uvw(3) - c = x*x + y*y - this%r2*z*z - quad = k*k - a*c - - if (quad < ZERO) then - ! no intersection with cone - - d = INFINITY - - elseif (coincident .or. abs(c) < FP_COINCIDENT) then - ! particle is on the cone, thus one distance is positive/negative and the - ! other is zero. The sign of k determines which distance is zero and which - ! is not. - - if (k >= ZERO) then - d = (-k - sqrt(quad))/a - else - d = (-k + sqrt(quad))/a - end if - - else - ! calculate both solutions to the quadratic - quad = sqrt(quad) - d = (-k - quad)/a - b = (-k + quad)/a - - ! determine the smallest positive solution - if (d < ZERO) then - if (b > ZERO) then - d = b - end if - else - if (b > ZERO) d = min(d, b) - end if - end if - - ! If the distance was negative, set boundary distance to infinity - if (d <= ZERO) d = INFINITY - end function z_cone_distance - - 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_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_distance(this, xyz, uvw, coincident) result(d) - class(SurfaceQuadric), intent(in) :: this - real(8), intent(in) :: xyz(3) - real(8), intent(in) :: uvw(3) - logical, intent(in) :: coincident - real(8) :: d - - real(8) :: a, k, c - real(8) :: quad, b - - associate (x => xyz(1), y => xyz(2), z => xyz(3), & - u => uvw(1), v => uvw(2), w => uvw(3)) - - a = this%A*u*u + this%B*v*v + this%C*w*w + this%D*u*v + this%E*v*w + & - this%F*u*w - k = (this%A*u*x + this%B*v*y + this%C*w*z + HALF*(this%D*(u*y + v*x) + & - this%E*(v*z + w*y) + this%F*(w*x + u*z) + this%G*u + this%H*v + & - this%J*w)) - c = this%A*x*x + this%B*y*y + this%C*z*z + this%D*x*y + this%E*y*z + & - this%F*x*z + this%G*x + this%H*y + this%J*z + this%K - quad = k*k - a*c - - if (quad < ZERO) then - ! no intersection with surface - - d = INFINITY - - elseif (coincident .or. abs(c) < FP_COINCIDENT) then - ! particle is on the surface, thus one distance is positive/negative and - ! the other is zero. The sign of k determines which distance is zero and - ! which is not. - - if (k >= ZERO) then - d = (-k - sqrt(quad))/a - else - d = (-k + sqrt(quad))/a - end if - - else - ! calculate both solutions to the quadratic - quad = sqrt(quad) - d = (-k - quad)/a - b = (-k + quad)/a - - ! determine the smallest positive solution - if (d < ZERO) then - if (b > ZERO) then - d = b - end if - else - if (b > ZERO) d = min(d, b) - end if - end if - - ! If the distance was negative, set boundary distance to infinity - if (d <= ZERO) d = INFINITY - end associate - end function quadric_distance - - 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 + integer(C_INT) :: id + id = surface_id_c(this % ptr) + end function surface_id + + pure function surface_bc(this) result(bc) + class(Surface), intent(in) :: this + integer(C_INT) :: bc + bc = surface_bc_c(this % ptr) + end function surface_bc + + pure function surface_sense(this, xyz, uvw) result(sense) + class(Surface), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + logical(C_BOOL) :: sense + sense = surface_sense_c(this % ptr, xyz, uvw) + end function surface_sense + + pure subroutine surface_reflect(this, xyz, uvw) + class(Surface), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3); + real(C_DOUBLE), intent(inout) :: uvw(3); + call surface_reflect_c(this % ptr, xyz, uvw) + end subroutine surface_reflect + + pure function surface_distance(this, xyz, uvw, coincident) result(d) + class(Surface), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3); + real(C_DOUBLE), intent(in) :: uvw(3); + logical(C_BOOL), intent(in) :: coincident; + real(C_DOUBLE) :: d; + d = surface_distance_c(this % ptr, xyz, uvw, coincident) + end function surface_distance + + pure subroutine surface_normal(this, xyz, uvw) + class(Surface), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3); + real(C_DOUBLE), intent(out) :: uvw(3); + call surface_normal_c(this % ptr, xyz, uvw) + end subroutine surface_normal + + subroutine surface_to_hdf5(this, group) + class(Surface), intent(in) :: this + integer(HID_T), intent(in) :: group + call surface_to_hdf5_c(this % ptr, group) + end subroutine surface_to_hdf5 + + pure function surface_i_periodic(this) result(i_periodic) + class(Surface), intent(in) :: this + integer(C_INT) :: i_periodic + i_periodic = surface_i_periodic_c(this % ptr) + end function surface_i_periodic + + function surface_periodic(this, other, xyz, uvw) result(rotational) + class(Surface), intent(in) :: this + class(Surface), intent(in) :: other + real(C_DOUBLE), intent(inout) :: xyz(3); + real(C_DOUBLE), intent(inout) :: uvw(3); + logical(C_BOOL) :: rotational + rotational = surface_periodic_c(this % ptr, other % ptr, xyz, uvw) + end function surface_periodic !=============================================================================== ! FREE_MEMORY_SURFACES deallocates global arrays defined in this module !=============================================================================== subroutine free_memory_surfaces() - n_surfaces = 0 if (allocated(surfaces)) deallocate(surfaces) call surface_dict % clear() + call free_memory_surfaces_c() end subroutine free_memory_surfaces end module surface_header diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 027e335e9..da0d8929c 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -4258,7 +4258,7 @@ contains real(C_DOUBLE) :: k_tra ! Copy of batch tracklength estimate of keff real(C_DOUBLE) :: val -#ifdef MPI +#ifdef OPENMC_MPI ! Combine tally results onto master process if (reduce_tallies) call reduce_tally_results() #endif @@ -4306,7 +4306,7 @@ contains ! REDUCE_TALLY_RESULTS collects all the results from tallies onto one processor !=============================================================================== -#ifdef MPI +#ifdef OPENMC_MPI subroutine reduce_tally_results() integer :: i diff --git a/src/tallies/tally_filter_surface.F90 b/src/tallies/tally_filter_surface.F90 index daecd8891..df3bf3603 100644 --- a/src/tallies/tally_filter_surface.F90 +++ b/src/tallies/tally_filter_surface.F90 @@ -105,7 +105,7 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Surface " // to_str(surfaces(this % surfaces(bin)) % obj % id) + label = "Surface " // to_str(surfaces(this % surfaces(bin)) % id()) end function text_label_surface end module tally_filter_surface diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 9e9b06d20..ee2259ed4 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -2,7 +2,7 @@ module trigger use, intrinsic :: ISO_C_BINDING -#ifdef MPI +#ifdef OPENMC_MPI use message_passing #endif diff --git a/src/tracking.F90 b/src/tracking.F90 index f65151701..127da0222 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -4,7 +4,7 @@ module tracking use cross_section, only: calculate_xs use error, only: warning, write_message use geometry_header, only: cells - use geometry, only: find_cell, distance_to_boundary, cross_lattice, & + use geometry, only: find_cell, distance_to_boundary, cross_lattice,& check_cell_overlap use message_passing use mgxs_header @@ -292,20 +292,20 @@ 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 logical :: found ! particle found in universe? class(Surface), pointer :: surf + class(Surface), pointer :: surf2 ! periodic partner surface i_surface = abs(p % surface) - surf => surfaces(i_surface)%obj + surf => surfaces(i_surface) if (verbosity >= 10 .or. trace) then - call write_message(" Crossing surface " // trim(to_str(surf % id))) + call write_message(" Crossing surface " // trim(to_str(surf % id()))) end if - if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then + if (surf % bc() == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then ! ======================================================================= ! PARTICLE LEAKS OUT OF PROBLEM @@ -330,11 +330,11 @@ contains ! Display message if (verbosity >= 10 .or. trace) then call write_message(" Leaked out of surface " & - &// trim(to_str(surf % id))) + &// trim(to_str(surf % id()))) end if return - elseif (surf % bc == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then + elseif (surf % bc() == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then ! ======================================================================= ! PARTICLE REFLECTS FROM SURFACE @@ -357,7 +357,7 @@ contains end if ! Reflect particle off surface - call surf%reflect(p%coord(1)%xyz, p%coord(1)%uvw) + call surf % reflect(p%coord(1)%xyz, p%coord(1)%uvw) ! Make sure new particle direction is normalized u = p%coord(1)%uvw(1) @@ -378,7 +378,7 @@ contains call find_cell(p, found) if (.not. found) then call p % mark_as_lost("Couldn't find particle after reflecting& - & from surface " // trim(to_str(surf % id)) // ".") + & from surface " // trim(to_str(surf % id())) // ".") return end if @@ -388,10 +388,10 @@ contains ! Diagnostic message if (verbosity >= 10 .or. trace) then call write_message(" Reflected from surface " & - &// trim(to_str(surf%id))) + &// trim(to_str(surf%id()))) end if return - elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then + elseif (surf % bc() == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then ! ======================================================================= ! PERIODIC BOUNDARY @@ -406,7 +406,6 @@ contains ! Score surface currents since reflection causes the direction of the ! particle to change -- artificially move the particle slightly back in ! case the surface crossing is coincident with a mesh boundary - if (active_current_tallies % size() > 0) then xyz = p % coord(1) % xyz p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw @@ -414,70 +413,19 @@ 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. + ! Get a pointer to the partner periodic surface. Offset the index to + ! correct for C vs. Fortran indexing. + surf2 => surfaces(surf % i_periodic() + 1) - ! 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 + ! Adjust the particle's location and direction. + rotational = surf2 % periodic_translate(surf, p % coord(1) % xyz, & + p % coord(1) % uvw) ! Reassign particle's surface if (rotational) then - p % surface = surf % i_periodic + p % surface = surf % i_periodic() + 1 else - p % surface = sign(surf % i_periodic, p % surface) + p % surface = sign(surf % i_periodic() + 1, p % surface) end if ! Figure out what cell particle is in now @@ -485,7 +433,8 @@ contains call find_cell(p, found) if (.not. found) then call p % mark_as_lost("Couldn't find particle after hitting & - &periodic boundary on surface " // trim(to_str(surf % id)) // ".") + &periodic boundary on surface " // trim(to_str(surf % id())) & + // ".") return end if @@ -495,7 +444,7 @@ contains ! Diagnostic message if (verbosity >= 10 .or. trace) then call write_message(" Hit periodic boundary on surface " & - // trim(to_str(surf%id))) + // trim(to_str(surf%id()))) end if return end if @@ -542,7 +491,7 @@ contains if (.not. found) then call p % mark_as_lost("After particle " // trim(to_str(p % id)) & - // " crossed surface " // trim(to_str(surf % id)) & + // " crossed surface " // trim(to_str(surf % id())) & // " it could not be located in any cell and it did not leak.") return end if diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index 5985a236b..07dc5b418 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -136,7 +136,7 @@ contains integer :: total_hits ! total hits for a single domain (summed over materials) integer :: min_samples ! minimum number of samples per process integer :: remainder ! leftover samples from uneven divide -#ifdef MPI +#ifdef OPENMC_MPI integer :: mpi_err ! MPI error code integer :: m ! index over materials integer :: n ! number of materials @@ -279,7 +279,7 @@ contains total_hits = 0 if (master) then -#ifdef MPI +#ifdef OPENMC_MPI do j = 1, n_procs - 1 call MPI_RECV(n, 1, MPI_INTEGER, j, 0, mpi_intracomm, & MPI_STATUS_IGNORE, mpi_err) @@ -341,7 +341,7 @@ contains end do else -#ifdef MPI +#ifdef OPENMC_MPI n = master_indices(i_domain) % size() allocate(data(2*n)) do k = 0, n - 1 diff --git a/src/xml_interface.h b/src/xml_interface.h new file mode 100644 index 000000000..778497562 --- /dev/null +++ b/src/xml_interface.h @@ -0,0 +1,48 @@ +#ifndef XML_INTERFACE_H +#define XML_INTERFACE_H + +#include // for std::transform +#include +#include + +#include "pugixml/pugixml.hpp" + + +namespace openmc { + +bool +check_for_node(pugi::xml_node node, const char *name) +{ + return node.attribute(name) || node.child(name); +} + + +std::string +get_node_value(pugi::xml_node node, const char *name) +{ + // Search for either an attribute or child tag and get the data as a char*. + const pugi::char_t *value_char; + if (node.attribute(name)) { + value_char = node.attribute(name).value(); + } else if (node.child(name)) { + value_char = node.child_value(name); + } else { + std::stringstream err_msg; + err_msg << "Node \"" << name << "\" is not a member of the \"" + << node.name() << "\" XML node"; + fatal_error(err_msg); + } + + // Convert to lowercase string. + std::string value(value_char); + std::transform(value.begin(), value.end(), value.begin(), ::tolower); + + // Remove whitespace. + value.erase(0, value.find_first_not_of(" \t\r\n")); + value.erase(value.find_last_not_of(" \t\r\n") + 1); + + return value; +} + +} // namespace openmc +#endif // XML_INTERFACE_H