crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 05:05:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #959 from smharper/c_surface

Browse source 
Transition surface implementation to C++
This commit is contained in:
Paul Romano 2018-02-03 12:30:53 -06:00 • committed by GitHub
commit 2d355266ee
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
33 changed files with 2181 additions and 1581 deletions
Download patch file Download diff file Expand all files Collapse all files

13
CMakeLists.txt
View file

@ -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)

5
openmc/capi/settings.py
View file

@ -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):

12
src/api.F90
View file

@ -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
!===============================================================================

8
src/cmfd_execute.F90
View file

@ -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

15
src/constants.F90
View file

@ -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 :: &

18
src/eigenvalue.F90
View file

@ -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)

12
src/error.F90
View file

@ -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.

34
src/error.h Normal file
View file

@ -0,0 +1,34 @@
#ifndef ERROR_H
#define ERROR_H
#include <cstring>
#include <string>
#include <sstream>
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

25
src/geometry.F90
View file

@ -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

4
src/hdf5_interface.F90
View file

@ -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

91
src/hdf5_interface.h Normal file
View file

@ -0,0 +1,91 @@
#ifndef HDF5_INTERFACE_H
#define HDF5_INTERFACE_H
#include <array>
#include <string.h>
#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<std::size_t array_len> void
write_double_1D(hid_t group_id, char const *name,
std::array<double, array_len> &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

15
src/initialize.F90
View file

@ -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

298
src/input_xml.F90
View file

@ -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 <surface>
call get_node_list(root, "surface", node_surf_list)
call read_surfaces(root % ptr)
! Get number of <surface> 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

8
src/main.F90
View file

@ -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

4
src/mesh.F90
View file

@ -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)

6
src/message_passing.F90
View file

@ -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

4
src/output.F90
View file

@ -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

13
src/random_lcg.F90
View file

@ -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

10
src/random_lcg.cpp
View file

@ -2,6 +2,9 @@
#include <cmath>
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

12
src/random_lcg.h
View file

@ -3,6 +3,9 @@
#include <cstdint>
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

8
src/simulation.F90
View file

@ -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)

2
src/source.F90
View file

@ -1,7 +1,7 @@
module source
use hdf5, only: HID_T
#ifdef MPI
#ifdef OPENMC_MPI
use message_passing
#endif

24
src/state_point.F90
View file

@ -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

93
src/summary.F90
View file

@ -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)

1235
src/surface.cpp Normal file
View file

File diff suppressed because it is too large Load diff

427
src/surface.h Normal file
View file

@ -0,0 +1,427 @@
#ifndef SURFACE_H
#define SURFACE_H
#include <map>
#include <limits> // For numeric_limits
#include <string>
#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<double>::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<int, int> 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

1207
src/surface_header.F90
View file

File diff suppressed because it is too large Load diff

4
src/tallies/tally.F90
View file

@ -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

2
src/tallies/tally_filter_surface.F90
View file

@ -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

2
src/tallies/trigger.F90
View file

@ -2,7 +2,7 @@ module trigger
use, intrinsic :: ISO_C_BINDING
#ifdef MPI
#ifdef OPENMC_MPI
use message_passing
#endif

97
src/tracking.F90
View file

@ -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

6
src/volume_calc.F90
View file

@ -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

48
src/xml_interface.h Normal file
View file

@ -0,0 +1,48 @@
#ifndef XML_INTERFACE_H
#define XML_INTERFACE_H
#include <algorithm> // for std::transform
#include <sstream>
#include <string>
#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