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

Merged in latest changes to master

Browse source
This commit is contained in:
Sterling Harper 2013-06-11 13:51:08 -04:00
commit b3b80884df
20 changed files with 3490 additions and 2501 deletions
Download patch file Download diff file Expand all files Collapse all files

27
src/DEPENDENCIES
View file

@ -162,6 +162,7 @@ global.o: cmfd_header.o
global.o: constants.o
global.o: dict_header.o
global.o: geometry_header.o
global.o: hdf5_interface.o
global.o: material_header.o
global.o: mesh_header.o
global.o: particle_header.o
@ -171,17 +172,16 @@ global.o: source_header.o
global.o: tally_header.o
global.o: timer_header.o
hdf5_interface.o: ace_header.o
hdf5_interface.o: bank_header.o
hdf5_interface.o: constants.o
hdf5_interface.o: endf.o
hdf5_interface.o: error.o
hdf5_interface.o: geometry_header.o
hdf5_interface.o: global.o
hdf5_interface.o: material_header.o
hdf5_interface.o: output.o
hdf5_interface.o: string.o
hdf5_interface.o: tally_header.o
hdf5_summary.o: ace_header.o
hdf5_summary.o: constants.o
hdf5_summary.o: endf.o
hdf5_summary.o: geometry_header.o
hdf5_summary.o: global.o
hdf5_summary.o: hdf5_interface.o
hdf5_summary.o: material_header.o
hdf5_summary.o: output.o
hdf5_summary.o: string.o
hdf5_summary.o: tally_header.o
initialize.o: ace.o
initialize.o: bank_header.o
@ -262,6 +262,9 @@ output.o: plot_header.o
output.o: string.o
output.o: tally_header.o
output_interface.o: hdf5_interface.o
output_interface.o: mpiio_interface.o
particle_header.o: constants.o
particle_restart.o: bank_header.o
@ -337,8 +340,8 @@ source.o: string.o
state_point.o: error.o
state_point.o: global.o
state_point.o: math.o
state_point.o: output.o
state_point.o: output_interface.o
state_point.o: string.o
state_point.o: tally_header.o

3
src/OBJECTS
View file

@ -29,6 +29,7 @@ geometry.o \
geometry_header.o \
global.o \
hdf5_interface.o \
hdf5_summary.o \
initialize.o \
interpolation.o \
input_xml.o \
@ -38,7 +39,9 @@ material_header.o \
math.o \
mesh_header.o \
mesh.o \
mpiio_interface.o \
output.o \
output_interface.o \
particle_header.o \
particle_restart.o \
particle_restart_write.o \

18
src/constants.F90
View file

@ -11,8 +11,13 @@ module constants
integer, parameter :: VERSION_RELEASE = 1
! Revision numbers for binary files
integer, parameter :: REVISION_SOURCE = 1
integer, parameter :: REVISION_STATEPOINT = 8
integer, parameter :: REVISION_STATEPOINT = 9
integer, parameter :: REVISION_PARTICLE_RESTART = 1
! Binary file types
integer, parameter :: &
FILETYPE_STATEPOINT = -1, &
FILETYPE_PARTICLE_RESTART = -2
! ============================================================================
! ADJUSTABLE PARAMETERS
@ -359,11 +364,10 @@ module constants
integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file
integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file
integer, parameter :: UNIT_XS = 14 ! unit # for writing xs summary file
integer, parameter :: UNIT_SOURCE = 15 ! unit # for writing source file
integer, parameter :: UNIT_STATE = 16 ! unit # for writing state point
integer, parameter :: CMFD_BALANCE = 17 ! unit # for writing cmfd balance file
integer, parameter :: UNIT_PARTICLE = 18 ! unit # for writing particle restart
integer, parameter :: UNIT_TRACK = 19 ! unit # for writing particle tracks
integer, parameter :: CMFD_BALANCE = 15 ! unit # for writing cmfd balance file
integer, parameter :: UNIT_PARTICLE = 16 ! unit # for writing particle restart
integer, parameter :: UNIT_OUTPUT = 17 ! unit # for writing output
integer, parameter :: UNIT_TRACK = 18 ! unit # for writing particle tracks
!=============================================================================
! CMFD CONSTANTS

207
src/eigenvalue.F90
View file

@ -17,18 +17,17 @@ module eigenvalue
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use source, only: get_source_particle
use state_point, only: write_state_point, replay_batch_history
use state_point, only: write_state_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
#ifdef HDF5
use hdf5_interface, only: hdf5_write_state_point
#endif
private
public :: run_eigenvalue
real(8) :: keff_generation ! Single-generation k on each processor
real(8) :: k_sum(2) = ZERO ! used to reduce sum and sum_sq
contains
!===============================================================================
@ -85,16 +84,7 @@ contains
! Accumulate time for transport
call time_transport % stop()
! Distribute fission bank across processors evenly
call time_bank % start()
call synchronize_bank()
call time_bank % stop()
! Calculate shannon entropy
if (entropy_on) call shannon_entropy()
! Write generation output
if (master .and. current_gen /= gen_per_batch) call print_generation()
call finalize_generation()
end do GENERATION_LOOP
@ -147,14 +137,43 @@ contains
subroutine initialize_generation()
! set overall generation number
overall_gen = gen_per_batch*(current_batch - 1) + current_gen
! Reset number of fission bank sites
n_bank = 0
! Count source sites if using uniform fission source weighting
if (ufs) call count_source_for_ufs()
! Store current value of tracklength k
keff_generation = global_tallies(K_TRACKLENGTH) % value
end subroutine initialize_generation
!===============================================================================
! FINALIZE_GENERATION
!===============================================================================
subroutine finalize_generation()
! Distribute fission bank across processors evenly
call time_bank % start()
call synchronize_bank()
call time_bank % stop()
! Calculate shannon entropy
if (entropy_on) call shannon_entropy()
! Collect results and statistics
call calculate_generation_keff()
call calculate_average_keff()
! Write generation output
if (master .and. current_gen /= gen_per_batch) call print_generation()
end subroutine finalize_generation
!===============================================================================
! FINALIZE_BATCH handles synchronization and accumulation of tallies,
! calculation of Shannon entropy, getting single-batch estimate of keff, and
@ -168,8 +187,11 @@ contains
call synchronize_tallies()
call time_tallies % stop()
! Collect results and statistics
call calculate_keff()
! Reset global tally results
if (.not. active_batches) then
call reset_result(global_tallies)
n_realizations = 0
end if
! Perform CMFD calculation if on
if (cmfd_on) call execute_cmfd()
@ -183,11 +205,7 @@ contains
if (master) call calculate_combined_keff()
! Create state point file
#ifdef HDF5
call hdf5_write_state_point()
#else
call write_state_point()
#endif
end if
if (master .and. current_batch == n_batches) then
@ -542,83 +560,59 @@ contains
end subroutine shannon_entropy
!===============================================================================
! CALCULATE_KEFF calculates the single batch estimate of keff as well as the
! mean and standard deviation of the mean for active batches
! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto
! the master processor and normalizes them. This should work whether or not the
! no-reduce method is being used.
!===============================================================================
subroutine calculate_keff()
subroutine calculate_generation_keff()
real(8) :: temp(2) ! used to reduce sum and sum_sq
real(8) :: alpha ! significance level for CI
real(8) :: t_value ! t-value for confidence intervals
message = "Calculate batch keff..."
call write_message(8)
! =========================================================================
! SINGLE-BATCH ESTIMATE OF K-EFFECTIVE
! Get keff for this generation by subtracting off the starting value
keff_generation = global_tallies(K_TRACKLENGTH) % value - keff_generation
#ifdef MPI
if (.not. reduce_tallies) then
! Reduce value of k_batch if running in parallel
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, k_batch(current_batch), 1, MPI_REAL8, &
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
else
! Receive buffer not significant at other processors
call MPI_REDUCE(k_batch(current_batch), temp, 1, MPI_REAL8, &
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
end if
end if
! Combine values across all processors
call MPI_ALLREDUCE(keff_generation, k_generation(overall_gen), 1, &
MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, mpi_err)
#else
k_generation(overall_gen) = keff_generation
#endif
! Normalize single batch estimate of k
if (master) then
k_batch(current_batch) = k_batch(current_batch) / &
(n_particles * gen_per_batch)
end if
! TODO: This should be normalized by total_weight, not by n_particles
k_generation(overall_gen) = k_generation(overall_gen) / n_particles
! =========================================================================
! ACCUMULATED K-EFFECTIVE AND ITS VARIANCE
end subroutine calculate_generation_keff
if (reduce_tallies) then
! In this case, global_tallies has already been reduced, so we don't
! need to perform any more reductions and just take the values from
! global_tallies directly
!===============================================================================
! CALCULATE_AVERAGE_KEFF calculates the mean and standard deviation of the mean
! of k-effective during active generations and broadcasts the mean to all
! processors
!===============================================================================
! Sample mean of keff
keff = global_tallies(K_TRACKLENGTH) % sum / n_realizations
subroutine calculate_average_keff()
if (n_realizations > 1) then
if (confidence_intervals) then
! Calculate t-value for confidence intervals
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1)
else
t_value = ONE
end if
integer :: n ! number of active generations
real(8) :: alpha ! significance level for CI
real(8) :: t_value ! t-value for confidence intervals
! Determine number of active generations
n = overall_gen - n_inactive*gen_per_batch
if (n <= 0) then
! For inactive generations, use current generation k as estimate for next
! generation
keff = k_generation(overall_gen)
! Standard deviation of the sample mean of k
keff_std = t_value * sqrt((global_tallies(K_TRACKLENGTH) % sum_sq / &
n_realizations - keff * keff) / (n_realizations - 1))
end if
else
! In this case, no reduce was ever done on global_tallies. Thus, we need
! to reduce the values in sum and sum^2 to get the sample mean and its
! standard deviation
! Sample mean of keff
k_sum(1) = k_sum(1) + k_generation(overall_gen)
k_sum(2) = k_sum(2) + k_generation(overall_gen)**2
#ifdef MPI
call MPI_REDUCE(global_tallies(K_TRACKLENGTH) % sum, temp, 2, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
#else
temp(1) = global_tallies(K_TRACKLENGTH) % sum
temp(2) = global_tallies(K_TRACKLENGTH) % sum_sq
#endif
! Determine mean
keff = k_sum(1) / n
! Sample mean of k
keff = temp(1) / n_realizations
if (n_realizations > 1) then
if (n > 1) then
if (confidence_intervals) then
! Calculate t-value for confidence intervals
alpha = ONE - CONFIDENCE_LEVEL
@ -628,23 +622,11 @@ contains
end if
! Standard deviation of the sample mean of k
keff_std = t_value * sqrt((temp(2)/n_realizations - keff*keff) / &
(n_realizations - 1))
keff_std = t_value * sqrt((k_sum(2)/n - keff**2) / (n - 1))
end if
end if
#ifdef MPI
! Broadcast new keff value to all processors
call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
#endif
! Reset global tally results
if (.not. active_batches) then
call reset_result(global_tallies)
n_realizations = 0
end if
end subroutine calculate_keff
end subroutine calculate_average_keff
!===============================================================================
! CALCULATE_COMBINED_KEFF calculates a minimum variance estimate of k-effective
@ -793,4 +775,37 @@ contains
end subroutine count_source_for_ufs
!===============================================================================
! REPLAY_BATCH_HISTORY displays keff and entropy for each generation within a
! batch using data read from a state point file
!===============================================================================
subroutine replay_batch_history
! Write message at beginning
if (current_batch == 1) then
message = "Replaying history from state point..."
call write_message(1)
end if
do current_gen = 1, gen_per_batch
overall_gen = overall_gen + 1
call calculate_average_keff()
! print out batch keff
if (current_gen < gen_per_batch) then
if (master) call print_generation()
else
if (master) call print_batch_keff()
end if
end do
! Write message at end
if (current_batch == restart_batch) then
message = "Resuming simulation..."
call write_message(1)
end if
end subroutine replay_batch_history
end module eigenvalue

10
src/finalize.F90
View file

@ -12,7 +12,7 @@ module finalize
#endif
#ifdef HDF5
use hdf5_interface, only: hdf5_finalize
use hdf5_interface, only: h5tclose_f, h5close_f, hdf5_err
#endif
implicit none
@ -56,8 +56,12 @@ contains
call free_memory()
#ifdef HDF5
! Close HDF5 interface and release memory
call hdf5_finalize()
! Release compound datatypes
call h5tclose_f(hdf5_tallyresult_t, hdf5_err)
call h5tclose_f(hdf5_bank_t, hdf5_err)
! Close FORTRAN interface.
call h5close_f(hdf5_err)
#endif
#ifdef MPI

8
src/fixed_source.F90
View file

@ -11,10 +11,6 @@ module fixed_source
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies
#ifdef HDF5
use hdf5_interface, only: hdf5_write_state_point
#endif
type(Bank), pointer :: source_site => null()
contains
@ -115,11 +111,7 @@ contains
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
#ifdef HDF5
call hdf5_write_state_point()
#else
call write_state_point()
#endif
end if
end subroutine finalize_batch

10
src/global.F90
View file

@ -21,7 +21,7 @@ module global
#endif
#ifdef HDF5
use hdf5
use hdf5_interface, only: HID_T
#endif
implicit none
@ -156,7 +156,8 @@ module global
integer :: n_active ! # of active batches
integer :: gen_per_batch = 1 ! # of generations per batch
integer :: current_batch = 0 ! current batch
integer :: current_gen = 0 ! current generation
integer :: current_gen = 0 ! current generation within a batch
integer :: overall_gen = 0 ! overall generation in the run
! External source
type(ExtSource), target :: external_source
@ -172,7 +173,7 @@ module global
integer(8) :: current_work ! index in source bank of current history simulated
! Temporary k-effective values
real(8), allocatable :: k_batch(:) ! batch estimates of k
real(8), allocatable :: k_generation(:) ! single-generation estimates of k
real(8) :: keff = ONE ! average k over active batches
real(8) :: keff_std ! standard deviation of average k
real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption
@ -182,7 +183,7 @@ module global
! Shannon entropy
logical :: entropy_on = .false.
real(8), allocatable :: entropy(:) ! shannon entropy at each batch
real(8), allocatable :: entropy(:) ! shannon entropy at each generation
real(8), allocatable :: entropy_p(:,:,:,:) ! % of source sites in each cell
type(StructuredMesh), pointer :: entropy_mesh
@ -244,7 +245,6 @@ module global
integer(HID_T) :: hdf5_tallyresult_t ! Compound type for TallyResult
integer(HID_T) :: hdf5_bank_t ! Compound type for Bank
integer(HID_T) :: hdf5_integer8_t ! type for integer(8)
integer :: hdf5_err ! error flag
#endif
! ============================================================================

1914
src/hdf5_interface.F90
View file

File diff suppressed because it is too large Load diff

715
src/hdf5_summary.F90 Normal file
View file

@ -0,0 +1,715 @@
module hdf5_summary
#ifdef HDF5
use ace_header, only: Reaction, UrrData, Nuclide
use constants
use endf, only: reaction_name
use geometry_header, only: Cell, Surface, Universe, Lattice
use global
use material_header, only: Material
use mesh_header, only: StructuredMesh
use output_interface
use output, only: time_stamp
use string, only: to_str
use tally_header, only: TallyObject
contains
!===============================================================================
! HDF5_WRITE_SUMMARY
!===============================================================================
subroutine hdf5_write_summary()
character(MAX_FILE_LEN) :: filename = "summary.h5"
! Create a new file using default properties.
call file_create(filename, "serial")
! Write header information
call hdf5_write_header()
! Write eigenvalue information
if (run_mode == MODE_EIGENVALUE) then
! Write number of particles
call write_data(n_particles, "n_particles")
! Use H5LT interface to write n_batches, n_inactive, and n_active
call write_data(n_batches, "n_batches")
call write_data(n_inactive, "n_inactive")
call write_data(n_active, "n_active")
call write_data(gen_per_batch, "gen_per_batch")
! Add description of each variable
call write_attribute_string("n_particles", &
"description", "Number of particles per generation")
call write_attribute_string("n_batches", &
"description", "Total number of batches")
call write_attribute_string("n_inactive", &
"description", "Number of inactive batches")
call write_attribute_string("n_active", &
"description", "Number of active batches")
call write_attribute_string("gen_per_batch", &
"description", "Number of generations per batch")
end if
call hdf5_write_geometry()
call hdf5_write_materials()
call hdf5_write_nuclides()
if (n_tallies > 0) then
call hdf5_write_tallies()
end if
! Terminate access to the file.
call file_close("serial")
end subroutine hdf5_write_summary
!===============================================================================
! HDF5_WRITE_HEADER
!===============================================================================
subroutine hdf5_write_header()
! Write version information
call write_data(VERSION_MAJOR, "version_major")
call write_data(VERSION_MINOR, "version_minor")
call write_data(VERSION_RELEASE, "version_release")
! Write current date and time
call write_data(time_stamp(), "date_and_time")
! Write MPI information
call write_data(n_procs, "n_procs")
call write_attribute_string("n_procs", "description", &
"Number of MPI processes")
end subroutine hdf5_write_header
!===============================================================================
! HDF5_WRITE_GEOMETRY
!===============================================================================
subroutine hdf5_write_geometry()
integer :: i, j, k, m
integer :: n_x, n_y, n_z
integer, allocatable :: lattice_universes(:,:,:)
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
type(Lattice), pointer :: lat => null()
! Use H5LT interface to write number of geometry objects
call write_data(n_cells, "n_cells", group="geometry")
call write_data(n_surfaces, "n_surfaces", group="geometry")
call write_data(n_universes, "n_universes", group="geometry")
call write_data(n_lattices, "n_lattices", group="geometry")
! ==========================================================================
! WRITE INFORMATION ON CELLS
! Create a cell group (nothing directly written in this group) then close
call hdf5_open_group("geometry/cells")
call hdf5_close_group()
! Write information on each cell
CELL_LOOP: do i = 1, n_cells
c => cells(i)
! Write universe for this cell
call write_data(universes(c % universe) % id, "universe", &
group="geometry/cells/cell " // trim(to_str(c % id)))
! Write information on what fills this cell
select case (c % type)
case (CELL_NORMAL)
call write_data("normal", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
if (c % material == MATERIAL_VOID) then
call write_data(-1, "material", &
group="geometry/cells/cell " // trim(to_str(c % id)))
else
call write_data(materials(c % material) % id, "material", &
group="geometry/cells/cell " // trim(to_str(c % id)))
end if
case (CELL_FILL)
call write_data("universe", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call write_data(universes(c % fill) % id, "material", &
group="geometry/cells/cell " // trim(to_str(c % id)))
case (CELL_LATTICE)
call write_data("lattice", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call write_data(lattices(c % fill) % id, "lattice", &
group="geometry/cells/cell " // trim(to_str(c % id)))
end select
! Write list of bounding surfaces
if (c % n_surfaces > 0) then
call write_data(c % surfaces, "surfaces", length= c % n_surfaces, &
group="geometry/cells/cell " // trim(to_str(c % id)))
end if
end do CELL_LOOP
! ==========================================================================
! WRITE INFORMATION ON SURFACES
! Create surfaces group (nothing directly written here) then close
call hdf5_open_group("geometry/surfaces")
call hdf5_close_group()
! Write information on each surface
SURFACE_LOOP: do i = 1, n_surfaces
s => surfaces(i)
! Write surface type
select case (s % type)
case (SURF_PX)
call write_data("X Plane", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_PY)
call write_data("Y Plane", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_PZ)
call write_data("Z Plane", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_PLANE)
call write_data("Plane", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CYL_X)
call write_data("X Cylinder", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CYL_Y)
call write_data("Y Cylinder", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CYL_Z)
call write_data("Z Cylinder", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_SPHERE)
call write_data("Sphere", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CONE_X)
call write_data("X Cone", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CONE_Y)
call write_data("Y Cone", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CONE_Z)
call write_data("Z Cone", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
end select
! Write coefficients for surface
call write_data(s % coeffs, "coefficients", length=size(s % coeffs), &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
! Write positive neighbors
if (allocated(s % neighbor_pos)) then
call write_data(s % neighbor_pos, "neighbors_positive", &
length=size(s % neighbor_pos), &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
end if
! Write negative neighbors
if (allocated(s % neighbor_neg)) then
call write_data(s % neighbor_neg, "neighbors_negative", &
length=size(s % neighbor_neg), &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
end if
! Write boundary condition
select case (s % bc)
case (BC_TRANSMIT)
call write_data("transmission", "boundary_condition", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (BC_VACUUM)
call write_data("vacuum", "boundary_condition", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (BC_REFLECT)
call write_data("reflective", "boundary_condition", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (BC_PERIODIC)
call write_data("periodic", "boundary_condition", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
end select
end do SURFACE_LOOP
! ==========================================================================
! WRITE INFORMATION ON UNIVERSES
! Create universes group (nothing directly written here) then close
call hdf5_open_group("geometry/universes")
call hdf5_close_group()
! Write information on each universe
UNIVERSE_LOOP: do i = 1, n_universes
u => universes(i)
! Write list of cells in this universe
if (u % n_cells > 0) then
call write_data(u % cells, "cells", length=u % n_cells, &
group="geometry/universes/universe " // trim(to_str(u % id)))
end if
end do UNIVERSE_LOOP
! ==========================================================================
! WRITE INFORMATION ON LATTICES
! Create lattices group (nothing directly written here) then close
call hdf5_open_group("geometry/lattices")
call hdf5_close_group()
! Write information on each lattice
LATTICE_LOOP: do i = 1, n_lattices
lat => lattices(i)
! Write lattice type
select case(lat % type)
case (LATTICE_RECT)
call write_data("rectangular", "type", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
case (LATTICE_HEX)
call write_data("hexagonal", "type", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end select
! Write lattice dimensions, lower left corner, and width of element
call write_data(lat % dimension, "dimension", &
length=lat % n_dimension, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
call write_data(lat % lower_left, "lower_left", &
length=lat % n_dimension, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
call write_data(lat % width, "width", &
length=lat % n_dimension, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
! Determine dimensions of lattice
n_x = lat % dimension(1)
n_y = lat % dimension(2)
if (lat % n_dimension == 3) then
n_z = lat % dimension(3)
else
n_z = 1
end if
! Write lattice universes
allocate(lattice_universes(n_x, n_y, n_z))
do j = 1, n_x
do k = 1, n_y
do m = 1, n_z
lattice_universes(j,k,m) = universes(lat % universes(j,k,m)) % id
end do
end do
end do
call write_data(lattice_universes, "universes", &
length=(/n_x, n_y, n_z/), &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
deallocate(lattice_universes)
end do LATTICE_LOOP
end subroutine hdf5_write_geometry
!===============================================================================
! HDF5_WRITE_MATERIALS
!===============================================================================
subroutine hdf5_write_materials()
integer :: i
integer :: j
integer, allocatable :: zaids(:)
type(Material), pointer :: m => null()
! Use H5LT interface to write number of materials
call write_data(n_materials, "n_materials", group="materials")
! Write information on each material
do i = 1, n_materials
m => materials(i)
! Write atom density with units
call write_data(m % density, "atom_density", &
group="materials/material " // trim(to_str(m % id)))
call write_attribute_string("atom_density", "units", "atom/b-cm", &
group="materials/material " // trim(to_str(m % id)))
! Copy ZAID for each nuclide to temporary array
allocate(zaids(m % n_nuclides))
do j = 1, m % n_nuclides
zaids(j) = nuclides(m % nuclide(j)) % zaid
end do
! Write temporary array to 'nuclides'
call write_data(zaids, "nuclides", length=m % n_nuclides, &
group="materials/material " // trim(to_str(m % id)))
! Deallocate temporary array
deallocate(zaids)
! Write atom densities
call write_data(m % atom_density, "nuclide_densities", &
length=m % n_nuclides, &
group="materials/material " // trim(to_str(m % id)))
! Write S(a,b) information if present
if (m % n_sab > 0) then
call write_data(m % i_sab_nuclides, "i_sab_nuclides", &
length=m % n_sab, &
group="materials/material " // trim(to_str(m % id)))
call write_data(m % i_sab_tables, "i_sab_tables", &
length=m % n_sab, &
group="materials/material " // trim(to_str(m % id)))
end if
end do
end subroutine hdf5_write_materials
!===============================================================================
! HDF5_WRITE_TALLIES
!===============================================================================
subroutine hdf5_write_tallies()
integer :: i, j
integer, allocatable :: temp_array(:) ! nuclide bin array
type(StructuredMesh), pointer :: m => null()
type(TallyObject), pointer :: t => null()
! Write total number of meshes
call write_data(n_meshes, "n_meshes", group="tallies")
! Write information for meshes
MESH_LOOP: do i = 1, n_meshes
m => meshes(i)
! Write type and number of dimensions
call write_data(m % type, "type", &
group="tallies/mesh " // trim(to_str(m % id)))
call write_data(m % n_dimension, "n_dimension", &
group="tallies/mesh " // trim(to_str(m % id)))
! Write mesh information
call write_data(m % dimension, "dimension", &
length=m % n_dimension, &
group="tallies/mesh " // trim(to_str(m % id)))
call write_data(m % lower_left, "lower_left", &
length=m % n_dimension, &
group="tallies/mesh " // trim(to_str(m % id)))
call write_data(m % upper_right, "upper_right", &
length=m % n_dimension, &
group="tallies/mesh " // trim(to_str(m % id)))
call write_data(m % width, "width", &
length=m % n_dimension, &
group="tallies/mesh " // trim(to_str(m % id)))
end do MESH_LOOP
! Write number of tallies
call write_data(n_tallies, "n_tallies", group="tallies")
TALLY_METADATA: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
! Write size of each tally
call write_data(t % total_score_bins, "total_score_bins", &
group="tallies/tally " // trim(to_str(t % id)))
call write_data(t % total_filter_bins, "total_filter_bins", &
group="tallies/tally " // trim(to_str(t % id)))
! Write number of filters
call write_data(t % n_filters, "n_filters", &
group="tallies/tally " // trim(to_str(t % id)))
FILTER_LOOP: do j = 1, t % n_filters
! Write type of filter
call write_data(t % filters(j) % type, "type", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
! Write number of bins for this filter
call write_data(t % filters(j) % n_bins, "n_bins", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
! Write filter bins
if (t % filters(j) % type == FILTER_ENERGYIN .or. &
t % filters(j) % type == FILTER_ENERGYOUT) then
call write_data(t % filters(j) % real_bins, "bins", &
length=size(t % filters(j) % real_bins), &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
else
call write_data(t % filters(j) % int_bins, "bins", &
length=size(t % filters(j) % int_bins), &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
end if
! Write name of type
select case (t % filters(j) % type)
case(FILTER_UNIVERSE)
call write_data("universe", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_MATERIAL)
call write_data("material", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_CELL)
call write_data("cell", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_CELLBORN)
call write_data("cellborn", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_SURFACE)
call write_data("surface", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_MESH)
call write_data("mesh", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_ENERGYIN)
call write_data("energy", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_ENERGYOUT)
call write_data("energyout", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
end select
end do FILTER_LOOP
! Write number of nuclide bins
call write_data(t % n_nuclide_bins, "n_nuclide_bins", &
group="tallies/tally " // trim(to_str(t % id)))
! Create temporary array for nuclide bins
allocate(temp_array(t % n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins
if (t % nuclide_bins(j) > 0) then
temp_array(j) = nuclides(t % nuclide_bins(j)) % zaid
else
temp_array(j) = t % nuclide_bins(j)
end if
end do NUCLIDE_LOOP
! Write and deallocate nuclide bins
call write_data(temp_array, "nuclide_bins", length=t % n_nuclide_bins, &
group="tallies/tally " // trim(to_str(t % id)))
deallocate(temp_array)
! Write number of score bins
call write_data(t % n_score_bins, "n_score_bins", &
group="tallies/tally " // trim(to_str(t % id)))
call write_data(t % score_bins, "score_bins", length=t % n_score_bins, &
group="tallies/tally " // trim(to_str(t % id)))
end do TALLY_METADATA
end subroutine hdf5_write_tallies
!===============================================================================
! HDF5_WRITE_NUCLIDES
!===============================================================================
subroutine hdf5_write_nuclides()
integer :: i, j
integer :: size_total
integer :: size_xs
integer :: size_angle
integer :: size_energy
type(Nuclide), pointer :: nuc => null()
type(Reaction), pointer :: rxn => null()
type(UrrData), pointer :: urr => null()
! Use H5LT interface to write number of nuclides
call write_data(n_nuclides_total, "n_nuclides", group="nuclides")
! Write information on each nuclide
NUCLIDE_LOOP: do i = 1, n_nuclides_total
nuc => nuclides(i)
! Determine size of cross-sections
size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8
size_total = size_xs
! Write some basic attributes
call write_data(nuc % zaid, "zaid", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % awr, "awr", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % kT, "kT", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % n_grid, "n_grid", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % n_reaction, "n_reactions", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % n_fission, "n_fission", &
group="nuclides/" // trim(nuc % name))
call write_data(size_xs, "size_xs", &
group="nuclides/" // trim(nuc % name))
! =======================================================================
! WRITE INFORMATION ON EACH REACTION
! Create overall group for reactions and close it
call hdf5_open_group("nuclides/" // trim(nuc % name) // "/reactions")
call hdf5_close_group()
RXN_LOOP: do j = 1, nuc % n_reaction
! Information on each reaction
rxn => nuc % reactions(j)
! Determine size of angle distribution
if (rxn % has_angle_dist) then
size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8
else
size_angle = 0
end if
! Determine size of energy distribution
if (rxn % has_energy_dist) then
size_energy = size(rxn % edist % data) * 8
else
size_energy = 0
end if
! Write information on reaction
call write_data(rxn % Q_value, "Q_value", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
call write_data(rxn % multiplicity, "multiplicity", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
call write_data(rxn % threshold, "threshold", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
call write_data(size_angle, "size_angle", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
call write_data(size_energy, "size_energy", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
! Accumulate data size
size_total = size_total + size_angle + size_energy
end do RXN_LOOP
! =======================================================================
! WRITE INFORMATION ON URR PROBABILITY TABLES
if (nuc % urr_present) then
urr => nuc % urr_data
call write_data(urr % n_energy, "urr_n_energy", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % n_prob, "urr_n_prob", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % interp, "urr_interp", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % inelastic_flag, "urr_inelastic", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % absorption_flag, "urr_absorption", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % energy(1), "urr_min_E", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % energy(urr % n_energy), "urr_max_E", &
group="nuclides/" // trim(nuc % name))
end if
! Write total memory used
call write_data(size_total, "size_total", &
group="nuclides/" // trim(nuc % name))
end do NUCLIDE_LOOP
end subroutine hdf5_write_nuclides
!===============================================================================
! HDF5_WRITE_TIMING
!===============================================================================
subroutine hdf5_write_timing()
integer(8) :: total_particles
real(8) :: speed
! Write timing data
call write_data(time_initialize % elapsed, "time_initialize", &
group="timing")
call write_data(time_read_xs % elapsed, "time_read_xs", &
group="timing")
call write_data(time_unionize % elapsed, "time_unionize", &
group="timing")
call write_data(time_transport % elapsed, "time_transport", &
group="timing")
call write_data(time_bank % elapsed, "time_bank", &
group="timing")
call write_data(time_bank_sample % elapsed, "time_bank_sample", &
group="timing")
call write_data(time_bank_sendrecv % elapsed, "time_bank_sendrecv", &
group="timing")
call write_data(time_tallies % elapsed, "time_tallies", &
group="timing")
call write_data(time_inactive % elapsed, "time_inactive", &
group="timing")
call write_data(time_active % elapsed, "time_active", &
group="timing")
call write_data(time_finalize % elapsed, "time_finalize", &
group="timing")
call write_data(time_total % elapsed, "time_total", &
group="timing")
! Add descriptions to timing data
call write_attribute_string("time_initialize", "description", &
"Total time elapsed for initialization (s)", group="timing")
call write_attribute_string("time_read_xs", "description", &
"Time reading cross-section libraries (s)", group="timing")
call write_attribute_string("time_unionize", "description", &
"Time unionizing energy grid (s)", group="timing")
call write_attribute_string("time_transport", "description", &
"Time in transport only (s)", group="timing")
call write_attribute_string("time_bank", "description", &
"Total time synchronizing fission bank (s)", group="timing")
call write_attribute_string("time_bank_sample", "description", &
"Time between generations sampling source sites (s)", group="timing")
call write_attribute_string("time_bank_sendrecv", "description", &
"Time between generations SEND/RECVing source sites (s)", &
group="timing")
call write_attribute_string("time_tallies", "description", &
"Time between batches accumulating tallies (s)", group="timing")
call write_attribute_string("time_inactive", "description", &
"Total time in inactive batches (s)", group="timing")
call write_attribute_string("time_active", "description", &
"Total time in active batches (s)", group="timing")
call write_attribute_string("time_finalize", "description", &
"Total time for finalization (s)", group="timing")
call write_attribute_string("time_total", "description", &
"Total time elapsed (s)", group="timing")
! Write calculation rate
total_particles = n_particles * n_batches * gen_per_batch
speed = real(total_particles) / (time_inactive % elapsed + &
time_active % elapsed)
call write_data(speed, "neutrons_per_second", group="timing")
end subroutine hdf5_write_timing
#endif
end module hdf5_summary

82
src/initialize.F90
View file

@ -13,11 +13,12 @@ module initialize
cells_in_univ_dict, read_plots_xml
use output, only: title, header, write_summary, print_version, &
print_usage, write_xs_summary, print_plot
use output_interface, only: file_open, file_close, read_data
use random_lcg, only: initialize_prng
use source, only: initialize_source
use state_point, only: load_state_point
use string, only: to_str, str_to_int, starts_with, ends_with
use tally_header, only: TallyObject
use tally_header, only: TallyObject, TallyResult
use tally_initialize, only: configure_tallies
#ifdef MPI
@ -25,8 +26,8 @@ module initialize
#endif
#ifdef HDF5
use hdf5_interface, only: hdf5_initialize, hdf5_write_summary, &
hdf5_load_state_point
use hdf5_interface
use hdf5_summary, only: hdf5_write_summary
#endif
implicit none
@ -121,11 +122,7 @@ contains
! If this is a restart run, load the state point data and binary source
! file
#ifdef HDF5
if (restart_run) call hdf5_load_state_point()
#else
if (restart_run) call load_state_point()
#endif
end if
if (master) then
@ -239,6 +236,52 @@ contains
end subroutine initialize_mpi
#endif
#ifdef HDF5
!===============================================================================
! HDF5_INITIALIZE
!===============================================================================
subroutine hdf5_initialize()
type(TallyResult), target :: tmp(2) ! temporary TallyResult
type(Bank), target :: tmpb(2) ! temporary Bank
integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals
integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates
! Initialize FORTRAN interface.
call h5open_f(hdf5_err)
! Create the compound datatype for TallyResult
call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmp(1)), &
c_loc(tmp(2))), hdf5_tallyresult_t, hdf5_err)
call h5tinsert_f(hdf5_tallyresult_t, "sum", h5offsetof(c_loc(tmp(1)), &
c_loc(tmp(1)%sum)), H5T_NATIVE_DOUBLE, hdf5_err)
call h5tinsert_f(hdf5_tallyresult_t, "sum_sq", h5offsetof(c_loc(tmp(1)), &
c_loc(tmp(1)%sum_sq)), H5T_NATIVE_DOUBLE, hdf5_err)
! Create compound type for xyz and uvw
call h5tarray_create_f(H5T_NATIVE_DOUBLE, 1, dims, coordinates_t, hdf5_err)
! Create the compound datatype for Bank
call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(2))), hdf5_bank_t, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "wgt", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%wgt)), H5T_NATIVE_DOUBLE, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "xyz", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%xyz)), coordinates_t, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "uvw", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%uvw)), coordinates_t, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "E", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%E)), H5T_NATIVE_DOUBLE, hdf5_err)
! Determine type for integer(8)
hdf5_integer8_t = h5kind_to_type(8, H5_INTEGER_KIND)
end subroutine hdf5_initialize
#endif
!===============================================================================
! READ_COMMAND_LINE reads all parameters from the command line
!===============================================================================
@ -248,6 +291,7 @@ contains
integer :: i ! loop index
integer :: argc ! number of command line arguments
integer :: last_flag ! index of last flag
integer :: filetype
character(MAX_FILE_LEN) :: pwd ! present working directory
character(MAX_WORD_LEN), allocatable :: argv(:) ! command line arguments
@ -284,10 +328,23 @@ contains
call fatal_error()
end if
case ('-r', '-restart', '--restart')
! Read path for state point
! Read path for state point/particle restart
i = i + 1
path_state_point = argv(i)
restart_run = .true.
! Check what type of file this is
call file_open(argv(i), 'parallel', 'r')
call read_data(filetype, 'filetype')
call file_close('parallel')
! Set path and flag for type of run
select case (filetype)
case (FILETYPE_STATEPOINT)
path_state_point = argv(i)
restart_run = .true.
case (FILETYPE_PARTICLE_RESTART)
path_particle_restart = argv(i)
particle_restart_run = .true.
end select
case ('-t', '-tallies', '--tallies')
run_mode = MODE_TALLIES
@ -306,11 +363,6 @@ contains
case ('-eps_tol', '-ksp_gmres_restart')
! Handle options that would be based to PETSC
i = i + 1
case ('-s','-particle','--particle')
! Read in path for particle restart
i = i + 1
path_particle_restart = argv(i)
particle_restart_run = .true.
case ('-tr', '-track', '--track')
write_track = .true.
case default

2
src/input_xml.F90
View file

@ -148,7 +148,7 @@ contains
gen_per_batch = eigenvalue_ % generations_per_batch
! Allocate array for batch keff and entropy
allocate(k_batch(n_batches))
allocate(k_generation(n_batches*gen_per_batch))
allocate(entropy(n_batches*gen_per_batch))
entropy = ZERO
end if

239
src/mpiio_interface.F90 Normal file
View file

@ -0,0 +1,239 @@
module mpiio_interface
#ifdef MPI
use mpi
implicit none
integer :: mpi_fh ! MPI file handle
integer :: mpiio_err ! MPI error code
contains
!===============================================================================
! MPI_CREATE_FILE creates a file using MPI file I/O
!===============================================================================
subroutine mpi_create_file(filename, fh)
character(*), intent(in) :: filename ! name of file to create
integer, intent(inout) :: fh ! file handle
! Create the file
call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, MPI_MODE_CREATE + &
MPI_MODE_WRONLY, MPI_INFO_NULL, fh, mpiio_err)
end subroutine mpi_create_file
!===============================================================================
! MPI_OPEN_FILE opens a file using MPI file I/O
!===============================================================================
subroutine mpi_open_file(filename, fh, mode)
character(*), intent(in) :: filename ! name of file to open
character(*), intent(in) :: mode ! open 'r' read, 'w' write
integer, intent(inout) :: fh ! file handle
integer :: open_mode
! Determine access mode
open_mode = MPI_MODE_RDONLY
if (mode == 'w') then
open_mode = MPI_MODE_WRONLY
end if
! Create the file
call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, &
open_mode, MPI_INFO_NULL, fh, mpiio_err)
end subroutine mpi_open_file
!===============================================================================
! MPI_CLOSE_FILE closes a file using MPI file I/O
!===============================================================================
subroutine mpi_close_file(fh)
integer, intent(inout) :: fh ! file handle
call MPI_FILE_CLOSE(fh, mpiio_err)
end subroutine mpi_close_file
!===============================================================================
! MPI_WRITE_INTEGER writes integer scalar data using MPI File I/O
!===============================================================================
subroutine mpi_write_integer(fh, buffer)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: buffer ! data to write
call MPI_FILE_WRITE(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_integer
!===============================================================================
! MPI_WRITE_INTEGER_1DARRAY writes integer 1-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_integer_1Darray(fh, buffer, length)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
integer, intent(in) :: buffer(:) ! data to write
call MPI_FILE_WRITE(fh, buffer, length, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_integer_1Darray
!===============================================================================
! MPI_WRITE_LONG writes long integer scalar data using MPI file I/O
!===============================================================================
subroutine mpi_write_long(fh, buffer)
integer, intent(in) :: fh ! file handle
integer(8), intent(in) :: buffer ! data to write
call MPI_FILE_WRITE(fh, buffer, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_long
!===============================================================================
! MPI_WRITE_DOUBLE writes double precision scalar data using MPI file I/O
!===============================================================================
subroutine mpi_write_double(fh, buffer)
integer, intent(in) :: fh ! file handle
real(8), intent(in) :: buffer ! data to write
call MPI_FILE_WRITE(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_double
!===============================================================================
! MPI_WRITE_DOUBLE_1DARRAY writes double precision 1-D array using MPI file I/O
!===============================================================================
subroutine mpi_write_double_1Darray(fh, buffer, length)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
real(8), intent(in) :: buffer(:) ! data to write
call MPI_FILE_WRITE(fh, buffer, length, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_double_1Darray
!===============================================================================
! MPI_WRITE_STRING writes string data using MPI file I/O
!===============================================================================
subroutine mpi_write_string(fh, buffer, length)
character(*), intent(in) :: buffer ! data to write
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of data
call MPI_FILE_WRITE(fh, buffer, length, MPI_CHARACTER, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_string
!===============================================================================
! MPI_READ_INTEGER reads integer scalar data using MPI file I/O
!===============================================================================
subroutine mpi_read_integer(fh, buffer)
integer, intent(in) :: fh ! file handle
integer, intent(inout) :: buffer ! read data to here
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_read_integer
!===============================================================================
! MPI_READ_INTEGER_1DARRAY reads integer 1-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_integer_1Darray(fh, buffer, length)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
integer, intent(inout) :: buffer(:) ! read data to here
call MPI_FILE_READ(fh, buffer, length, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_read_integer_1Darray
!===============================================================================
! MPI_READ_LONG reads long integer scalar data using MPI file I/O
!===============================================================================
subroutine mpi_read_long(fh, buffer)
integer, intent(in) :: fh ! file handle
integer(8), intent(inout) :: buffer ! read data to here
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_read_long
!===============================================================================
! MPI_READ_DOUBLE reads double precision scalar data using MPI file I/O
!===============================================================================
subroutine mpi_read_double(fh, buffer)
integer, intent(in) :: fh ! file handle
real(8), intent(inout) :: buffer ! read data to here
call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_read_double
!===============================================================================
! MPI_READ_DOUBLE_1DARRAY reads double precision 1-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_double_1Darray(fh, buffer, length)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
real(8), intent(inout) :: buffer(:) ! read data to here
call MPI_FILE_READ(fh, buffer, length, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_read_double_1Darray
!===============================================================================
! MPI_READ_STRING reads string data using MPI file I/O
!===============================================================================
subroutine mpi_read_string(fh, buffer, length)
character(*), intent(inout) :: buffer ! read data to here
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of string
call MPI_FILE_READ(fh, buffer, length, MPI_CHARACTER, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_read_string
#endif
end module mpiio_interface

26
src/output.F90
View file

@ -166,8 +166,8 @@ contains
write(OUTPUT_UNIT,*)
write(OUTPUT_UNIT,*) 'Options:'
write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode'
write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run'
write(OUTPUT_UNIT,*) ' -s, --particle Run a single particle history'
write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run from a state point'
write(OUTPUT_UNIT,*) ' or a particle restart file'
write(OUTPUT_UNIT,*) ' -t, --tallies Write tally results from state point'
write(OUTPUT_UNIT,*) ' -v, --version Show version information'
write(OUTPUT_UNIT,*) ' -?, --help Show this message'
@ -1215,24 +1215,24 @@ contains
if (entropy_on) then
if (cmfd_run) then
message = " Bat./Gen. k(batch) Entropy Average k CMFD k CMFD Ent"
message = " Bat./Gen. k Entropy Average k CMFD k CMFD Ent"
call write_message(1)
message = " ========= ======== ======== ==================== ======== ========"
call write_message(1)
else
message = " Bat./Gen. k(batch) Entropy Average k"
message = " Bat./Gen. k Entropy Average k"
call write_message(1)
message = " ========= ======== ======== ===================="
call write_message(1)
end if
else
if (cmfd_run) then
message = " Bat./Gen. k(batch) Average k CMFD k"
message = " Bat./Gen. k Average k CMFD k"
call write_message(1)
message = " ========= ======== ==================== ========"
call write_message(1)
else
message = " Bat./Gen. k(batch) Average k"
message = " Bat./Gen. k Average k"
call write_message(1)
message = " ========= ======== ===================="
call write_message(1)
@ -1251,11 +1251,17 @@ contains
! write out information about batch and generation
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
trim(to_str(current_batch)) // "/" // trim(to_str(current_gen))
write(UNIT=OUTPUT_UNIT, FMT='(11X)', ADVANCE='NO')
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') &
k_generation(overall_gen)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(current_gen + gen_per_batch*(current_batch - 1))
entropy(overall_gen)
if (overall_gen - n_inactive*gen_per_batch > 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
keff, keff_std
end if
! next line
write(UNIT=OUTPUT_UNIT, FMT=*)
@ -1273,14 +1279,14 @@ contains
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
trim(to_str(current_batch)) // "/" // trim(to_str(gen_per_batch))
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') &
k_batch(current_batch)
k_generation(overall_gen)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(current_batch*gen_per_batch)
! write out accumulated k-effective if after first active batch
if (current_batch > n_inactive + 1) then
if (overall_gen - n_inactive*gen_per_batch > 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
keff, keff_std
else

1127
src/output_interface.F90 Normal file
View file

File diff suppressed because it is too large Load diff

19
src/particle_restart.F90
View file

@ -7,7 +7,7 @@ module particle_restart
use geometry_header, only: BASE_UNIVERSE
use global
use particle_header, only: deallocate_coord
use output, only: write_message
use output, only: write_message, print_particle
use physics, only: transport
use random_lcg, only: set_particle_seed
use source, only: initialize_particle
@ -40,8 +40,6 @@ contains
subroutine read_hdf5_particle_restart()
integer(HSIZE_T) :: dims1(1)
! write meessage
message = "Loading particle restart file " // trim(path_particle_restart) &
// "..."
@ -55,8 +53,8 @@ contains
call hdf5_read_integer(hdf5_particle_file, 'current_batch', current_batch)
call hdf5_read_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch)
call hdf5_read_integer(hdf5_particle_file, 'current_gen', current_gen)
call hdf5_read_long(hdf5_particle_file, 'n_particles', n_particles)
call hdf5_read_long(hdf5_particle_file, 'id', p % id)
call hdf5_read_long(hdf5_particle_file, 'n_particles', n_particles, hdf5_integer8_t)
call hdf5_read_long(hdf5_particle_file, 'id', p % id, hdf5_integer8_t)
call hdf5_read_double(hdf5_particle_file, 'weight', p % wgt)
call hdf5_read_double(hdf5_particle_file, 'energy', p % E)
dims1 = (/3/)
@ -83,6 +81,9 @@ contains
subroutine read_binary_particle_restart()
integer :: filetype
integer :: revision
! write meessage
message = "Loading particle restart file " // trim(path_particle_restart) &
// "..."
@ -93,6 +94,8 @@ contains
ACCESS='stream')
! read data from file
read(UNIT_PARTICLE) filetype
read(UNIT_PARTICLE) revision
read(UNIT_PARTICLE) current_batch
read(UNIT_PARTICLE) gen_per_batch
read(UNIT_PARTICLE) current_gen
@ -158,11 +161,7 @@ contains
endif
! write output if particle made it
write(ou,*) 'Particle Successfully Transport:'
write(ou,*) 'WEIGHT:', p % wgt
write(ou,*) 'ENERGY:', p % E
write(ou,*) 'LOCATION:', p % coord % xyz
write(ou,*) 'ANGLE:', p % coord % uvw
call print_particle()
end subroutine run_particle_restart

11
src/particle_restart_write.F90
View file

@ -47,7 +47,6 @@ contains
subroutine write_particle_restart_hdf5()
character(MAX_FILE_LEN) :: filename
integer(HSIZE_T) :: dims1(1)
type(Bank), pointer :: src => null()
! set up file name
@ -61,11 +60,15 @@ contains
src => source_bank(current_work)
! write data to file
call hdf5_write_integer(hdf5_particle_file, 'filetype', &
FILETYPE_PARTICLE_RESTART)
call hdf5_write_integer(hdf5_particle_file, 'revision', &
REVISION_PARTICLE_RESTART)
call hdf5_write_integer(hdf5_particle_file, 'current_batch', current_batch)
call hdf5_write_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch)
call hdf5_write_integer(hdf5_particle_file, 'current_gen', current_gen)
call hdf5_write_long(hdf5_particle_file, 'n_particles', n_particles)
call hdf5_write_long(hdf5_particle_file, 'id', p % id)
call hdf5_write_long(hdf5_particle_file, 'n_particles', n_particles, hdf5_integer8_t)
call hdf5_write_long(hdf5_particle_file, 'id', p % id, hdf5_integer8_t)
call hdf5_write_double(hdf5_particle_file, 'weight', src % wgt)
call hdf5_write_double(hdf5_particle_file, 'energy', src % E)
dims1 = (/3/)
@ -102,6 +105,8 @@ contains
src => source_bank(current_work)
! write data to file
write(UNIT_PARTICLE) FILETYPE_PARTICLE_RESTART
write(UNIT_PARTICLE) REVISION_PARTICLE_RESTART
write(UNIT_PARTICLE) current_batch
write(UNIT_PARTICLE) gen_per_batch
write(UNIT_PARTICLE) current_gen

7
src/source.F90
View file

@ -9,7 +9,6 @@ module source
use particle_header, only: deallocate_coord
use physics, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn, set_particle_seed
use state_point, only: read_source_binary
use string, only: to_str
#ifdef MPI
@ -38,7 +37,8 @@ contains
! Read the source from a binary file instead of sampling from some
! assumed source distribution
call read_source_binary()
message = 'This feature is currently disabled and will be added back in.'
call fatal_error()
else
! Generation source sites from specified distribution in user input
@ -167,8 +167,7 @@ contains
p % id = bank_first + index_source - 1
! set random number seed
particle_seed = ((current_batch - 1)*gen_per_batch + &
current_gen - 1)*n_particles + p % id
particle_seed = (overall_gen - 1)*n_particles + p % id
call set_particle_seed(particle_seed)
! set particle trace

1552
src/state_point.F90
View file

File diff suppressed because it is too large Load diff

6
src/tally.F90
View file

@ -1757,12 +1757,6 @@ contains
end do
if (run_mode == MODE_EIGENVALUE) then
! Get the current batch estimate of k_analog for displaying to output
! --- this has to be performed after reduce_tally_values and before
! accumulate_result
k_batch(current_batch) = global_tallies(K_TRACKLENGTH) % value
if (active_batches) then
! Accumulate products of different estimators of k
k_col = global_tallies(K_COLLISION) % value / total_weight

8
src/utils/statepoint.py
View file

@ -147,8 +147,11 @@ class StatePoint(object):
self._read_metadata()
def _read_metadata(self):
# Read filetype
self.filetype = self._get_int(path='filetype')[0]
# Read statepoint revision
self.revision = self._get_int(path='revision_statepoint')[0]
self.revision = self._get_int(path='revision')[0]
# Read OpenMC version
if self._hdf5:
@ -179,7 +182,8 @@ class StatePoint(object):
if self.run_mode == 2:
self.n_inactive = self._get_int(path='n_inactive')[0]
self.gen_per_batch = self._get_int(path='gen_per_batch')[0]
self.k_batch = self._get_double(self.current_batch, path='k_batch')
self.k_batch = self._get_double(
self.current_batch*self.gen_per_batch, path='k_generation')
self.entropy = self._get_double(
self.current_batch*self.gen_per_batch, path='entropy')
if self.revision >= 8: