Start moving surface->summary.h5 to C++

This commit is contained in:
Sterling Harper 2018-01-19 18:12:57 -05:00
parent 7af296902b
commit 57991b271d
23 changed files with 311 additions and 251 deletions

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,6 +435,7 @@ set(LIBOPENMC_FORTRAN_SRC
src/tallies/trigger_header.F90
)
set(LIBOPENMC_CXX_SRC
src/hdf5_interface.h
src/random_lcg.h
src/random_lcg.cpp
src/surface_header.C

View file

@ -171,7 +171,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

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

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)

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

View file

@ -537,7 +537,6 @@ contains
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)

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

47
src/hdf5_interface.h Normal file
View file

@ -0,0 +1,47 @@
#ifndef HDF5_INTERFACE_H
#define HDF5_INTERFACE_H
#include <array> // For std::array
#include <string.h> // For strlen
#include "hdf5.h"
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);
}
#endif //HDF5_INTERFACE_H

View file

@ -52,8 +52,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 +74,7 @@ contains
call time_total%start()
call time_initialize%start()
#ifdef MPI
#ifdef OPENMC_MPI
! Setup MPI
call initialize_mpi(comm)
#endif
@ -108,7 +108,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 +116,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 +124,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

View file

@ -54,7 +54,7 @@ module input_xml
implicit none
type(C_PTR) :: node_ptr
end subroutine read_surfaces
end interface
end interface
contains
@ -1081,77 +1081,23 @@ contains
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
xmin = min(xmin, coeffs(1))
xmax = max(xmax, coeffs(1))
if (xmin == coeffs(1)) i_xmin = i
if (xmax == coeffs(1)) 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
ymin = min(ymin, coeffs(1))
ymax = max(ymax, coeffs(1))
if (ymin == coeffs(1)) i_ymin = i
if (ymax == coeffs(1)) 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)
zmin = min(zmin, coeffs(1))
zmax = max(zmax, coeffs(1))
if (zmin == coeffs(1)) i_zmin = i
if (zmax == coeffs(1)) i_zmax = i
end select
! No longer need coefficients

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

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)

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

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

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)

View file

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

View file

@ -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
@ -849,7 +849,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 +897,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 +907,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 +933,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

View file

@ -24,6 +24,17 @@ module summary
public :: write_summary
interface
subroutine surface_to_hdf5_c(surf_ind, group) &
bind(C, name='surface_to_hdf5')
use ISO_C_BINDING
use hdf5
implicit none
integer(C_INT), intent(in), value :: surf_ind
integer(HID_T), intent(in), value :: group
end subroutine surface_to_hdf5_c
end interface
contains
!===============================================================================
@ -125,7 +136,6 @@ contains
integer(HID_T) :: surfaces_group, surface_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
@ -245,87 +255,11 @@ contains
surface_group = create_group(surfaces_group, "surface " // &
trim(to_str(s%id)))
call surface_to_hdf5_c(i-1, surface_group)
! 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)
end do SURFACE_LOOP

View file

@ -1,8 +1,12 @@
#include <array> // For std::array
#include <cstring> // For strcmp
#include <limits> // For numeric_limits
#include <math.h> // For fabs
#include "pugixml/pugixml.hpp"
#include "hdf5.h"
#include "hdf5_interface.h"
// DEBUGGING
#include <typeinfo>
@ -15,6 +19,11 @@
const double FP_COINCIDENT = 1e-12;
const double INFTY = std::numeric_limits<double>::max();
const int BC_TRANSMIT{0};
const int BC_VACUUM{1};
const int BC_REFLECT{2};
const int BC_PERIODIC{3};
//==============================================================================
// Global array of surfaces
//==============================================================================
@ -127,6 +136,13 @@ public:
//! @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.
virtual void to_hdf5(hid_t group_id) const = 0;
protected:
void write_bc_to_hdf5(hid_t group_id) const;
};
bool
@ -164,6 +180,25 @@ Surface::reflect(const double xyz[3], double uvw[3]) const
uvw[2] -= 2.0 * projection / magnitude * norm[2];
}
void
Surface::write_bc_to_hdf5(hid_t group_id) const
{
switch(bc) {
case BC_TRANSMIT :
write_string(group_id, "boundary_type", "transmission");
break;
case BC_VACUUM :
write_string(group_id, "boundary_type", "vacuum");
break;
case BC_REFLECT :
write_string(group_id, "boundary_type", "reflective");
break;
case BC_PERIODIC :
write_string(group_id, "boundary_type", "periodic");
break;
}
}
//==============================================================================
//! A `Surface` that supports periodic boundary conditions.
//!
@ -236,6 +271,7 @@ public:
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(hid_t group_id) const;
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
const;
};
@ -261,6 +297,14 @@ 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(hid_t group_id) const
{
write_string(group_id, "type", "x-plane");
std::array<double, 1> coeffs{{x0}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
double uvw[3]) const
{
@ -285,7 +329,7 @@ bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
return true;
}
};
}
//==============================================================================
// SurfaceYPlane
@ -303,6 +347,7 @@ public:
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(hid_t group_id) const;
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
const;
};
@ -328,6 +373,14 @@ 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(hid_t group_id) const
{
write_string(group_id, "type", "y-plane");
std::array<double, 1> coeffs{{y0}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
double uvw[3]) const
{
@ -352,7 +405,7 @@ bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
return true;
}
};
}
//==============================================================================
// SurfaceZPlane
@ -370,6 +423,7 @@ public:
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(hid_t group_id) const;
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
const;
};
@ -395,13 +449,21 @@ 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(hid_t group_id) const
{
write_string(group_id, "type", "z-plane");
std::array<double, 1> coeffs{{z0}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
double uvw[3]) const
{
// Assume the other plane is aligned along z. Just change the z coord.
xyz[2] = z0;
return false;
};
}
//==============================================================================
// SurfacePlane
@ -419,6 +481,7 @@ public:
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(hid_t group_id) const;
bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
const;
};
@ -457,6 +520,14 @@ SurfacePlane::normal(const double xyz[3], double uvw[3]) const
uvw[2] = C;
}
void SurfacePlane::to_hdf5(hid_t group_id) const
{
write_string(group_id, "type", "plane");
std::array<double, 4> coeffs{{A, B, C, D}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3],
double uvw[3]) const
{
@ -564,6 +635,7 @@ public:
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(hid_t group_id) const;
};
SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node)
@ -588,6 +660,15 @@ 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(hid_t group_id) const
{
write_string(group_id, "type", "x-cylinder");
std::array<double, 3> coeffs{{y0, z0, r}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
//==============================================================================
// SurfaceYCylinder
//! A cylinder aligned along the y-axis.
@ -605,6 +686,7 @@ public:
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(hid_t group_id) const;
};
SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node)
@ -629,6 +711,14 @@ 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(hid_t group_id) const
{
write_string(group_id, "type", "y-cylinder");
std::array<double, 3> coeffs{{x0, z0, r}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
//==============================================================================
// SurfaceZCylinder
//! A cylinder aligned along the z-axis.
@ -646,6 +736,7 @@ public:
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(hid_t group_id) const;
};
SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node)
@ -670,6 +761,14 @@ 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(hid_t group_id) const
{
write_string(group_id, "type", "z-cylinder");
std::array<double, 3> coeffs{{x0, y0, r}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
//==============================================================================
// SurfaceSphere
//! A sphere.
@ -687,6 +786,7 @@ public:
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(hid_t group_id) const;
};
SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node)
@ -748,6 +848,14 @@ inline void SurfaceSphere::normal(const double xyz[3], double uvw[3]) const
uvw[2] = 2.0 * (xyz[2] - z0);
}
void SurfaceSphere::to_hdf5(hid_t group_id) const
{
write_string(group_id, "type", "sphere");
std::array<double, 4> coeffs{{x0, y0, z0, r}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
//==============================================================================
// Generic functions for x-, y-, and z-, cones
//==============================================================================
@ -848,6 +956,7 @@ public:
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(hid_t group_id) const;
};
SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node)
@ -872,6 +981,14 @@ 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(hid_t group_id) const
{
write_string(group_id, "type", "x-cone");
std::array<double, 4> coeffs{{x0, y0, z0, r_sq}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
//==============================================================================
// SurfaceYCone
//! A cone aligned along the y-axis.
@ -889,6 +1006,7 @@ public:
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(hid_t group_id) const;
};
SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node)
@ -913,6 +1031,14 @@ 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(hid_t group_id) const
{
write_string(group_id, "type", "y-cone");
std::array<double, 4> coeffs{{x0, y0, z0, r_sq}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
//==============================================================================
// SurfaceZCone
//! A cone aligned along the z-axis.
@ -930,6 +1056,7 @@ public:
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(hid_t group_id) const;
};
SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node)
@ -954,6 +1081,14 @@ 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(hid_t group_id) const
{
write_string(group_id, "type", "z-cone");
std::array<double, 4> coeffs{{x0, y0, z0, r_sq}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
//==============================================================================
// SurfaceQuadric
//! A general surface described by a quadratic equation.
@ -971,6 +1106,7 @@ public:
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(hid_t group_id) const;
};
SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node)
@ -1055,6 +1191,14 @@ SurfaceQuadric::normal(const double xyz[3], double uvw[3]) const
uvw[2] = 2.0*C*z + E*y + F*x + J;
}
void SurfaceQuadric::to_hdf5(hid_t group_id) const
{
write_string(group_id, "type", "quadric");
std::array<double, 10> coeffs{{A, B, C, D, E, F, G, H, J, K}};
write_double_1D(group_id, "coefficients", coeffs);
write_bc_to_hdf5(group_id);
}
//==============================================================================
extern "C" void
@ -1084,6 +1228,7 @@ read_surfaces(pugi::xml_node *node)
} else {
std::cout << "ERROR: Found a surface with no type attribute/node"
<< std::endl;
//TODO: call fatal_error
}
if (strcmp(surf_type, "x-plane") == 0) {
@ -1125,6 +1270,33 @@ read_surfaces(pugi::xml_node *node)
} else {
std::cout << "Call error or handle uppercase here!" << std::endl;
std::cout << surf_type << std::endl;
//TODO: call fatal_error
}
const pugi::char_t *surf_bc{""};
if (surf_node.attribute("boundary")) {
surf_bc = surf_node.attribute("boundary").value();
} else if (surf_node.child("boundary")) {
surf_bc = surf_node.attribute("boundary").value();
}
if (strcmp(surf_bc, "transmission") == 0
or strcmp(surf_bc, "transmit") == 0
or strcmp(surf_bc, "") == 0) {
surfaces_c[i_surf]->bc = BC_TRANSMIT;
} else if (strcmp(surf_bc, "vacuum") == 0) {
surfaces_c[i_surf]->bc = BC_VACUUM;
} else if (strcmp(surf_bc, "reflective") == 0
or strcmp(surf_bc, "reflect") == 0
or strcmp(surf_bc, "reflecting") == 0) {
surfaces_c[i_surf]->bc = BC_REFLECT;
} else if (strcmp(surf_bc, "periodic") == 0) {
surfaces_c[i_surf]->bc = BC_PERIODIC;
} else {
std::cout << "Unknown boundary condition" << std::endl;
//TODO: call fatal_error
}
}
}
@ -1156,6 +1328,12 @@ surface_normal(int surf_ind, double xyz[3], double uvw[3])
return surfaces_c[surf_ind]->normal(xyz, uvw);
}
extern "C" void
surface_to_hdf5(int surf_ind, hid_t group)
{
surfaces_c[surf_ind]->to_hdf5(group);
}
extern "C" bool
surface_periodic(int surf_ind1, int surf_ind2, double xyz[3], double uvw[3])
{

View file

@ -36,84 +36,39 @@ module surface_header
!===============================================================================
type, extends(Surface) :: SurfaceXPlane
! x = x0
real(8) :: x0
end type SurfaceXPlane
type, extends(Surface) :: SurfaceYPlane
! y = y0
real(8) :: y0
end type SurfaceYPlane
type, extends(Surface) :: SurfaceZPlane
! z = z0
real(8) :: z0
end type SurfaceZPlane
type, extends(Surface) :: SurfacePlane
! Ax + By + Cz = D
real(8) :: A
real(8) :: B
real(8) :: C
real(8) :: D
end type SurfacePlane
type, extends(Surface) :: SurfaceXCylinder
! (y - y0)^2 + (z - z0)^2 = R^2
real(8) :: y0
real(8) :: z0
real(8) :: r
end type SurfaceXCylinder
type, extends(Surface) :: SurfaceYCylinder
! (x - x0)^2 + (z - z0)^2 = R^2
real(8) :: x0
real(8) :: z0
real(8) :: r
end type SurfaceYCylinder
type, extends(Surface) :: SurfaceZCylinder
! (x - x0)^2 + (y - y0)^2 = R^2
real(8) :: x0
real(8) :: y0
real(8) :: r
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
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
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
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
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
end type SurfaceQuadric
integer(C_INT32_T), bind(C) :: n_surfaces ! # of surfaces

View file

@ -4241,7 +4241,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
@ -4289,7 +4289,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

View file

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

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