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

Merge remote-tracking branch 'upstream/develop' into develop

Browse source
This commit is contained in:
Will Boyd 2017-02-27 11:14:32 -05:00
commit 647bf77a57
185 changed files with 5132 additions and 5834 deletions
Download patch file Download diff file Expand all files Collapse all files

99
CMakeLists.txt
View file

@ -266,15 +266,104 @@ add_library(faddeeva STATIC src/Faddeeva.c)
#===============================================================================
set(program "openmc")
file(GLOB source src/*.F90 src/xml/openmc_fox.F90)
add_executable(${program} ${source})
set(LIBOPENMC_FORTRAN_SRC
src/algorithm.F90
src/angle_distribution.F90
src/angleenergy_header.F90
src/bank_header.F90
src/cmfd_data.F90
src/cmfd_execute.F90
src/cmfd_header.F90
src/cmfd_input.F90
src/cmfd_loss_operator.F90
src/cmfd_prod_operator.F90
src/cmfd_solver.F90
src/constants.F90
src/cross_section.F90
src/dict_header.F90
src/distribution_multivariate.F90
src/distribution_univariate.F90
src/doppler.F90
src/eigenvalue.F90
src/endf.F90
src/endf_header.F90
src/energy_distribution.F90
src/energy_grid.F90
src/error.F90
src/finalize.F90
src/geometry.F90
src/geometry_header.F90
src/global.F90
src/hdf5_interface.F90
src/initialize.F90
src/input_xml.F90
src/list_header.F90
src/material_header.F90
src/math.F90
src/matrix_header.F90
src/mesh.F90
src/mesh_header.F90
src/message_passing.F90
src/mgxs_data.F90
src/mgxs_header.F90
src/multipole.F90
src/multipole_header.F90
src/nuclide_header.F90
src/output.F90
src/particle_header.F90
src/particle_restart.F90
src/particle_restart_write.F90
src/physics_common.F90
src/physics.F90
src/physics_mg.F90
src/plot.F90
src/plot_header.F90
src/product_header.F90
src/progress_header.F90
src/random_lcg.F90
src/reaction_header.F90
src/relaxng
src/sab_header.F90
src/scattdata_header.F90
src/secondary_correlated.F90
src/secondary_kalbach.F90
src/secondary_nbody.F90
src/secondary_uncorrelated.F90
src/set_header.F90
src/simulation.F90
src/source.F90
src/source_header.F90
src/state_point.F90
src/stl_vector.F90
src/string.F90
src/summary.F90
src/surface_header.F90
src/tally.F90
src/tally_filter.F90
src/tally_filter_header.F90
src/tally_header.F90
src/tally_initialize.F90
src/timer_header.F90
src/tracking.F90
src/track_output.F90
src/trigger.F90
src/trigger_header.F90
src/urr_header.F90
src/vector_header.F90
src/volume_calc.F90
src/volume_header.F90
src/xml_interface.F90
src/xml/openmc_fox.F90)
add_library(libopenmc STATIC ${LIBOPENMC_FORTRAN_SRC})
set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc)
add_executable(${program} src/main.F90)
# target_include_directories was added in CMake 2.8.11 and is the recommended
# way to set include directories. For lesser versions, we revert to set_property
if(CMAKE_VERSION VERSION_LESS 2.8.11)
include_directories(${HDF5_INCLUDE_DIRS})
else()
target_include_directories(${program} PUBLIC ${HDF5_INCLUDE_DIRS})
target_include_directories(libopenmc PUBLIC ${HDF5_INCLUDE_DIRS})
endif()
# target_compile_options was added in CMake 2.8.12 and is the recommended way to
@ -288,6 +377,7 @@ if (CMAKE_VERSION VERSION_LESS 2.8.12)
set_property(TARGET faddeeva PROPERTY COMPILE_FLAGS "${cflags}")
else()
target_compile_options(${program} PUBLIC ${f90flags})
target_compile_options(libopenmc PUBLIC ${f90flags})
target_compile_options(faddeeva PRIVATE ${cflags})
endif()
@ -298,7 +388,8 @@ endforeach()
# target_link_libraries treats any arguments starting with - but not -l as
# linker flags. Thus, we can pass both linker flags and libraries together.
target_link_libraries(${program} ${ldflags} ${HDF5_LIBRARIES} fox_dom faddeeva)
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} fox_dom faddeeva)
target_link_libraries(${program} ${ldflags} libopenmc)
#===============================================================================
# Install executable, scripts, manpage, license

2
docs/source/conf.py
View file

@ -248,7 +248,7 @@ napoleon_use_ivar = True
intersphinx_mapping = {
'python': ('https://docs.python.org/3', None),
'numpy': ('http://docs.scipy.org/doc/numpy/', None),
'numpy': ('https://docs.scipy.org/doc/numpy/', None),
'pandas': ('http://pandas.pydata.org/pandas-docs/stable/', None),
'matplotlib': ('http://matplotlib.org/', None)
}

74
docs/source/io_formats/particle_restart.rst
View file

@ -4,55 +4,31 @@
Particle Restart File Format
============================
The current revision of the particle restart file format is 1.
The current version of the particle restart file format is 2.0.
**/filetype** (*char[]*)
**/**
String indicating the type of file.
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
- **version** (*int[2]*) -- Major and minor version of the particle
restart file format.
- **openmc_version** (*int[3]*) -- Major, minor, and release
version number for OpenMC.
- **git_sha1** (*char[40]*) -- Git commit SHA-1 hash.
**/revision** (*int*)
Revision of the particle restart file format. Any time a change is made in
the format, this integer is incremented.
**/current_batch** (*int*)
The number of batches already simulated.
**/gen_per_batch** (*int*)
Number of generations per batch.
**/current_gen** (*int*)
The number of generations already simulated.
**/n_particles** (*int8_t*)
Number of particles used per generation.
**/run_mode** (*int*)
Run mode used. A value of 1 indicates a fixed-source run and a value of 2
indicates an eigenvalue run.
**/id** (*int8_t*)
Unique identifier of the particle.
**/weight** (*double*)
Weight of the particle.
**/energy** (*double*)
Energy of the particle in eV for continuous-energy mode, or the energy
group of the particle for multi-group mode.
**/xyz** (*double[3]*)
Position of the particle.
**/uvw** (*double[3]*)
Direction of the particle.
:Datasets: - **current_batch** (*int*) -- The number of batches already
simulated.
- **generations_per_batch** (*int*) -- Number of generations per
batch.
- **current_generation** (*int*) -- The number of generations already
simulated.
- **n_particles** (*int8_t*) -- Number of particles used per
generation.
- **run_mode** (*char[]*) -- Run mode used, either 'fixed source',
'eigenvalue', or 'particle restart'.
- **id** (*int8_t*) -- Unique identifier of the particle.
- **weight** (*double*) -- Weight of the particle.
- **energy** (*double*) -- Energy of the particle in eV for
continuous-energy mode, or the energy group of the particle for
multi-group mode.
- **xyz** (*double[3]*) -- Position of the particle.
- **uvw** (*double[3]*) -- Direction of the particle.

519
docs/source/io_formats/statepoint.rst
View file

@ -4,354 +4,171 @@
State Point File Format
=======================
The current revision of the statepoint file format is 15.
**/filetype** (*char[]*)
String indicating the type of file.
**/revision** (*int*)
Revision of the state point file format. Any time a change is made in the
format, this integer is incremented.
**/version_major** (*int*)
Major version number for OpenMC
**/version_minor** (*int*)
Minor version number for OpenMC
**/version_release** (*int*)
Release version number for OpenMC
**/date_and_time** (*char[]*)
Date and time the state point was written.
**/path** (*char[]*)
Absolute path to directory containing input files.
**/seed** (*int8_t*)
Pseudo-random number generator seed.
**/run_CE** (*int*)
Flag to denote continuous-energy or multi-group mode. A value of 1
indicates a continuous-energy run while a value of 0 indicates a
multi-group run.
**/run_mode** (*char[]*)
Run mode used. A value of 1 indicates a fixed-source run and a value of 2
indicates an eigenvalue run.
**/n_particles** (*int8_t*)
Number of particles used per generation.
**/n_batches** (*int*)
Number of batches to simulate.
**/current_batch** (*int*)
The number of batches already simulated.
if run_mode == 'k-eigenvalue':
**/n_inactive** (*int*)
Number of inactive batches.
**/gen_per_batch** (*int*)
Number of generations per batch.
**/k_generation** (*double[]*)
k-effective for each generation simulated.
**/entropy** (*double[]*)
Shannon entropy for each generation simulated
**/k_col_abs** (*double*)
Sum of product of collision/absorption estimates of k-effective
**/k_col_tra** (*double*)
Sum of product of collision/track-length estimates of k-effective
**/k_abs_tra** (*double*)
Sum of product of absorption/track-length estimates of k-effective
**/k_combined** (*double[2]*)
Mean and standard deviation of a combined estimate of k-effective
**/cmfd_on** (*int*)
Flag indicating whether CMFD is on (1) or off (0).
if (cmfd_on)
**/cmfd/indices** (*int[4]*)
Indices for cmfd mesh (i,j,k,g)
**/cmfd/k_cmfd** (*double[]*)
CMFD eigenvalues
**/cmfd/cmfd_src** (*double[][][][]*)
CMFD fission source
**/cmfd/cmfd_entropy** (*double[]*)
CMFD estimate of Shannon entropy
**/cmfd/cmfd_balance** (*double[]*)
RMS of the residual neutron balance equation on CMFD mesh
**/cmfd/cmfd_dominance** (*double[]*)
CMFD estimate of dominance ratio
**/cmfd/cmfd_srccmp** (*double[]*)
RMS comparison of difference between OpenMC and CMFD fission source
**/tallies/n_meshes** (*int*)
Number of meshes in tallies.xml file
**/tally/meshes/ids** (*int[]*)
Internal unique ID of each mesh.
**/tally/meshes/keys** (*int[]*)
User-identified unique ID of each mesh.
**/tallies/meshes/mesh <uid>/type** (*char[]*)
Type of mesh.
**/tallies/meshes/mesh <uid>/dimension** (*int*)
Number of mesh cells in each dimension.
**/tallies/meshes/mesh <uid>/lower_left** (*double[]*)
Coordinates of lower-left corner of mesh.
**/tallies/meshes/mesh <uid>/upper_right** (*double[]*)
Coordinates of upper-right corner of mesh.
**/tallies/meshes/mesh <uid>/width** (*double[]*)
Width of each mesh cell in each dimension.
**/tallies/derivatives/derivative <id>/independent variable** (*char[]*)
Independent variable of tally derivative
**/tallies/derivatives/derivative <id>/material** (*int*)
ID of the perturbed material
**/tallies/derivatives/derivative <id>/nuclide** (*char[]*)
Alias of the perturbed nuclide
**/tallies/n_tallies** (*int*)
Number of user-defined tallies.
**/tallies/ids** (*int[]*)
Internal unique ID of each tally.
**/tallies/keys** (*int[]*)
User-identified unique ID of each tally.
**/tallies/tally <uid>/estimator** (*char[]*)
Type of tally estimator, either 'analog', 'tracklength', or 'collision'.
**/tallies/tally <uid>/n_realizations** (*int*)
Number of realizations.
**/tallies/tally <uid>/n_filters** (*int*)
Number of filters used.
**/tallies/tally <uid>/filter <j>/type** (*char[]*)
Type of the j-th filter. Can be 'universe', 'material', 'cell', 'cellborn',
'surface', 'mesh', 'energy', 'energyout', 'distribcell', 'mu', 'polar',
'azimuthal', 'delayedgroup', or 'energyfunction'.
**/tallies/tally <uid>/filter <j>/n_bins** (*int*)
Number of bins for the j-th filter. Not present for 'energyfunction'
filters.
**/tallies/tally <uid>/filter <j>/bins** (*int[]* or *double[]*)
Value for each filter bin of this type. Not present for 'energyfunction'
filters.
**/tallies/tally <uid>/filter <j>/energy** (*double[]*)
Energy grid points for energyfunction interpolation. Only used for
'energyfunction' filters.
**/tallies/tally <uid>/filter <j>/y** (*double[]*)
Interpolant values for energyfunction interpolation. Only used for
'energyfunction' filters.
**/tallies/tally <uid>/nuclides** (*char[][]*)
Array of nuclides to tally. Note that if no nuclide is specified in the user
input, a single 'total' nuclide appears here.
**/tallies/tally <uid>/derivative** (*int*)
ID of the derivative applied to the tally.
**/tallies/tally <uid>/n_score_bins** (*int*)
Number of scoring bins for a single nuclide. In general, this can be greater
than the number of user-specified scores since each score might have
multiple scoring bins, e.g., scatter-PN.
**/tallies/tally <uid>/score_bins** (*char[][]*)
Values of specified scores.
**/tallies/tally <uid>/n_user_scores** (*int*)
Number of scores without accounting for those added by expansions,
e.g. scatter-PN.
**/tallies/tally <uid>/moment_orders** (*char[][]*)
Tallying moment orders for Legendre and spherical harmonic tally expansions
(*e.g.*, 'P2', 'Y1,2', etc.).
**/tallies/tally <uid>/results** (*double[][][2]*)
Accumulated sum and sum-of-squares for each bin of the i-th tally. The first
dimension represents combinations of filter bins, the second dimensions
represents scoring bins, and the third dimension has two entries for the sum
and the sum-of-squares.
**/source_present** (*int*)
Flag indicated if source bank is present in the file
**/n_realizations** (*int*)
Number of realizations for global tallies.
**/n_global_tallies** (*int*)
Number of global tally scores.
**/global_tallies** (Compound type)
Accumulated sum and sum-of-squares for each global tally. The compound type
has fields named ``sum`` and ``sum_sq``.
**/tallies_present** (*int*)
Flag indicated if tallies are present in the file.
if (run_mode == 'k-eigenvalue' and source_present > 0)
**/source_bank** (Compound type)
Source bank information for each particle. The compound type has fields
``wgt``, ``xyz``, ``uvw``, ``E``, ``g``, and ``delayed_group``, which
represent the weight, position, direction, energy, energy group, and
delayed_group of the source particle, respectively.
**/runtime/total initialization** (*double*)
Time (in seconds on the master process) spent reading inputs, allocating
arrays, etc.
**/runtime/reading cross sections** (*double*)
Time (in seconds on the master process) spent loading cross section
libraries (this is a subset of initialization).
**/runtime/simulation** (*double*)
Time (in seconds on the master process) spent between initialization and
finalization.
**/runtime/transport** (*double*)
Time (in seconds on the master process) spent transporting particles.
**/runtime/inactive batches** (*double*)
Time (in seconds on the master process) spent in the inactive batches
(including non-transport activities like communcating sites).
**/runtime/active batches** (*double*)
Time (in seconds on the master process) spent in the active batches
(including non-transport activities like communicating sites).
**/runtime/synchronizing fission bank** (*double*)
Time (in seconds on the master process) spent sampling source particles
from fission sites and communicating them to other processes for load
balancing.
**/runtime/sampling source sites** (*double*)
Time (in seconds on the master process) spent sampling source particles
from fission sites.
**/runtime/SEND-RECV source sites** (*double*)
Time (in seconds on the master process) spent communicating source sites
between processes for load balancing.
**/runtime/accumulating tallies** (*double*)
Time (in seconds on the master process) spent communicating tally results
and evaluating their statistics.
**/runtime/CMFD** (*double*)
Time (in seconds on the master process) spent evaluating CMFD.
**/runtime/CMFD building matrices** (*double*)
Time (in seconds on the master process) spent buliding CMFD matrices.
**/runtime/CMFD solving matrices** (*double*)
Time (in seconds on the master process) spent solving CMFD matrices.
**/runtime/total** (*double*)
Total time spent (in seconds on the master process) in the program.
The current version of the statepoint file format is 16.0.
**/**
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
- **version** (*int[2]*) -- Major and minor version of the
statepoint file format.
- **openmc_version** (*int[3]*) -- Major, minor, and release
version number for OpenMC.
- **git_sha1** (*char[40]*) -- Git commit SHA-1 hash.
- **date_and_time** (*char[]*) -- Date and time the summary was
written.
- **path** (*char[]*) -- Path to directory containing input files.
- **cmfd_on** (*int*) -- Flag indicating whether CMFD is on (1) or
off (0).
- **tallies_present** (*int*) -- Flag indicating whether tallies
are present (1) or not (0).
- **source_present** (*int*) -- Flag indicating whether the source
bank is present (1) or not (0).
:Datasets: - **seed** (*int8_t*) -- Pseudo-random number generator seed.
- **energy_mode** (*char[]*) -- Energy mode of the run, either
'continuous-energy' or 'multi-group'.
- **run_mode** (*char[]*) -- Run mode used, either 'eigenvalue' or
'fixed source'.
- **n_particles** (*int8_t*) -- Number of particles used per generation.
- **n_batches** (*int*) -- Number of batches to simulate.
- **current_batch** (*int*) -- The number of batches already simulated.
- **n_inactive** (*int*) -- Number of inactive batches. Only present
when `run_mode` is 'eigenvalue'.
- **generations_per_batch** (*int*) -- Number of generations per
batch. Only present when `run_mode` is 'eigenvalue'.
- **k_generation** (*double[]*) -- k-effective for each generation
simulated.
- **entropy** (*double[]*) -- Shannon entropy for each generation
simulated.
- **k_col_abs** (*double*) -- Sum of product of collision/absorption
estimates of k-effective.
- **k_col_tra** (*double*) -- Sum of product of
collision/track-length estimates of k-effective.
- **k_abs_tra** (*double*) -- Sum of product of
absorption/track-length estimates of k-effective.
- **k_combined** (*double[2]*) -- Mean and standard deviation of a
combined estimate of k-effective.
- **n_realizations** (*int*) -- Number of realizations for global
tallies.
- **global_tallies** (*double[][2]*) -- Accumulated sum and
sum-of-squares for each global tally.
- **source_bank** (Compound type) -- Source bank information for each
particle. The compound type has fields ``wgt``, ``xyz``, ``uvw``,
``E``, ``g``, and ``delayed_group``, which represent the weight,
position, direction, energy, energy group, and delayed_group of the
source particle, respectively. Only present when `run_mode` is
'eigenvalue'.
**/cmfd/**
:Datasets: - **indices** (*int[4]*) -- Indices for cmfd mesh (i,j,k,g)
- **k_cmfd** (*double[]*) -- CMFD eigenvalues
- **cmfd_src** (*double[][][][]*) -- CMFD fission source
- **cmfd_entropy** (*double[]*) -- CMFD estimate of Shannon entropy
- **cmfd_balance** (*double[]*) -- RMS of the residual neutron
balance equation on CMFD mesh
- **cmfd_dominance** (*double[]*) -- CMFD estimate of dominance ratio
- **cmfd_srccmp** (*double[]*) -- RMS comparison of difference
between OpenMC and CMFD fission source
**/tallies/**
:Attributes: - **n_tallies** (*int*) -- Number of user-defined tallies.
- **ids** (*int[]*) -- User-defined unique ID of each tally.
**/tallies/meshes/**
:Attributes: - **n_meshes** (*int*) -- Number of meshes in the problem.
- **ids** (*int[]*) -- User-defined unique ID of each mesh.
**/tallies/meshes/mesh <uid>/**
:Datasets: - **type** (*char[]*) -- Type of mesh.
- **dimension** (*int*) -- Number of mesh cells in each dimension.
- **lower_left** (*double[]*) -- Coordinates of lower-left corner of
mesh.
- **upper_right** (*double[]*) -- Coordinates of upper-right corner
of mesh.
- **width** (*double[]*) -- Width of each mesh cell in each
dimension.
**/tallies/derivatives/derivative <id>/**
:Datasets: - **independent variable** (*char[]*) -- Independent variable of
tally derivative.
- **material** (*int*) -- ID of the perturbed material.
- **nuclide** (*char[]*) -- Alias of the perturbed nuclide.
- **estimator** (*char[]*) -- Type of tally estimator, either
'analog', 'tracklength', or 'collision'.
**/tallies/tally <uid>/**
:Datasets: - **n_realizations** (*int*) -- Number of realizations.
- **n_filters** (*int*) -- Number of filters used.
- **nuclides** (*char[][]*) -- Array of nuclides to tally. Note that
if no nuclide is specified in the user input, a single 'total'
nuclide appears here.
- **derivative** (*int*) -- ID of the derivative applied to the
tally.
- **n_score_bins** (*int*) -- Number of scoring bins for a single
nuclide. In general, this can be greater than the number of
user-specified scores since each score might have multiple scoring
bins, e.g., scatter-PN.
- **score_bins** (*char[][]*) -- Values of specified scores.
- **n_user_scores** (*int*) -- Number of scores without accounting
for those added by expansions, e.g. scatter-PN.
- **moment_orders** (*char[][]*) -- Tallying moment orders for
Legendre and spherical harmonic tally expansions (e.g., 'P2',
'Y1,2', etc.).
- **results** (*double[][][2]*) -- Accumulated sum and sum-of-squares
for each bin of the i-th tally. The first dimension represents
combinations of filter bins, the second dimensions represents
scoring bins, and the third dimension has two entries for the sum
and the sum-of-squares.
**/tallies/tally <uid>/filter <j>/**
:Datasets: - **type** (*char[]*) -- Type of the j-th filter. Can be 'universe',
'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
'energyout', 'distribcell', 'mu', 'polar', 'azimuthal',
'delayedgroup', or 'energyfunction'.
- **n_bins** (*int*) -- Number of bins for the j-th filter. Not
present for 'energyfunction' filters.
- **bins** (*int[]* or *double[]*) -- Value for each filter bin of
this type. Not present for 'energyfunction' filters.
- **energy** (*double[]*) -- Energy grid points for energyfunction
interpolation. Only used for 'energyfunction' filters.
- **y** (*double[]*) -- Interpolant values for energyfunction
interpolation. Only used for 'energyfunction' filters.
**/runtime/**
All values are given in seconds and are measured on the master process.
:Datasets: - **total initialization** (*double*) -- Time spent reading inputs,
allocating arrays, etc.
- **reading cross sections** (*double*) -- Time spent loading cross
section libraries (this is a subset of initialization).
- **simulation** (*double*) -- Time spent between initialization and
finalization.
- **transport** (*double*) -- Time spent transporting particles.
- **inactive batches** (*double*) -- Time spent in the inactive
batches (including non-transport activities like communcating
sites).
- **active batches** (*double*) -- Time spent in the active batches
(including non-transport activities like communicating sites).
- **synchronizing fission bank** (*double*) -- Time spent sampling
source particles from fission sites and communicating them to other
processes for load balancing.
- **sampling source sites** (*double*) -- Time spent sampling source
particles from fission sites.
- **SEND-RECV source sites** (*double*) -- Time spent communicating
source sites between processes for load balancing.
- **accumulating tallies** (*double*) -- Time spent communicating
tally results and evaluating their statistics.
- **CMFD** (*double*) -- Time spent evaluating CMFD.
- **CMFD building matrices** (*double*) -- Time spent buliding CMFD
matrices.
- **CMFD solving matrices** (*double*) -- Time spent solving CMFD
matrices.
- **total** (*double*) -- Total time spent in the program.

380
docs/source/io_formats/summary.rst
View file

@ -4,255 +4,131 @@
Summary File Format
===================
The current revision of the summary file format is 4.
**/filetype** (*char[]*)
String indicating the type of file.
**/revision** (*int*)
Revision of the summary file format. Any time a change is made in the
format, this integer is incremented.
**/version_major** (*int*)
Major version number for OpenMC
**/version_minor** (*int*)
Minor version number for OpenMC
**/version_release** (*int*)
Release version number for OpenMC
**/date_and_time** (*char[]*)
Date and time the summary was written.
**/n_procs** (*int*)
Number of MPI processes used.
**/n_particles** (*int8_t*)
Number of particles used per generation.
**/n_batches** (*int*)
Number of batches to simulate.
**/n_inactive** (*int*)
Number of inactive batches. Only present if /run_mode is set to
'k-eigenvalue'.
**/n_active** (*int*)
Number of active batches. Only present if /run_mode is set to
'k-eigenvalue'.
**/gen_per_batch** (*int*)
Number of generations per batch. Only present if /run_mode is set to
'k-eigenvalue'.
**/geometry/n_cells** (*int*)
Number of cells in the problem.
**/geometry/n_surfaces** (*int*)
Number of surfaces in the problem.
**/geometry/n_universes** (*int*)
Number of unique universes in the problem.
**/geometry/n_lattices** (*int*)
Number of lattices in the problem.
**/geometry/cells/cell <uid>/index** (*int*)
Index in cells array used internally in OpenMC.
**/geometry/cells/cell <uid>/name** (*char[]*)
Name of the cell.
**/geometry/cells/cell <uid>/universe** (*int*)
Universe assigned to the cell. If none is specified, the default
universe (0) is assigned.
**/geometry/cells/cell <uid>/fill_type** (*char[]*)
Type of fill for the cell. Can be 'normal', 'universe', or 'lattice'.
**/geometry/cells/cell <uid>/material** (*int* or *int[]*)
Unique ID of the material(s) assigned to the cell. This dataset is present
only if fill_type is set to 'normal'. The value '-1' signifies void
material. The data is an array if the cell uses distributed materials,
otherwise it is a scalar.
**/geometry/cells/cell <uid>/temperature** (*double[]*)
Temperature of the cell in Kelvin.
**/geometry/cells/cell <uid>/offset** (*int[]*)
Offsets used for distribcell tally filter. This dataset is present only if
fill_type is set to 'universe'.
**/geometry/cells/cell <uid>/translation** (*double[3]*)
Translation applied to the fill universe. This dataset is present only if
fill_type is set to 'universe'.
**/geometry/cells/cell <uid>/rotation** (*double[3]*)
Angles in degrees about the x-, y-, and z-axes for which the fill universe
should be rotated. This dataset is present only if fill_type is set to
'universe'.
**/geometry/cells/cell <uid>/lattice** (*int*)
Unique ID of the lattice which fills the cell. Only present if fill_type is
set to 'lattice'.
**/geometry/cells/cell <uid>/region** (*char[]*)
Region specification for the cell.
**/geometry/cells/cell <uid>/distribcell_index** (*int*)
Index of this cell in distribcell arrays. Only present if this cell is
listed in a distribcell filter or if it uses distributed materials.
**/geometry/cells/cell <uid>/paths** (*char[][]*)
The paths traversed through the CSG tree to reach each distribcell
instance. This consists of the integer IDs for each universe, cell and
lattice delimited by '->'. Each lattice cell is specified by its (x,y) or
(x,y,z) indices. Only present if this cell is listed in a distribcell filter
or if it uses distributed materials.
**/geometry/surfaces/surface <uid>/index** (*int*)
Index in surfaces array used internally in OpenMC.
**/geometry/surfaces/surface <uid>/name** (*char[]*)
Name of the surface.
**/geometry/surfaces/surface <uid>/type** (*char[]*)
Type of the surface. Can be 'x-plane', 'y-plane', 'z-plane', 'plane',
'x-cylinder', 'y-cylinder', 'sphere', 'x-cone', 'y-cone', 'z-cone', or
'quadric'.
**/geometry/surfaces/surface <uid>/coefficients** (*double[]*)
Array of coefficients that define the surface. See :ref:`surface_element`
for what coefficients are defined for each surface type.
**/geometry/surfaces/surface <uid>/boundary_condition** (*char[]*)
Boundary condition applied to the surface. Can be 'transmission', 'vacuum',
'reflective', or 'periodic'.
**/geometry/universes/universe <uid>/index** (*int*)
Index in the universes array used internally in OpenMC.
**/geometry/universes/universe <uid>/cells** (*int[]*)
Array of unique IDs of cells that appear in the universe.
**/geometry/lattices/lattice <uid>/index** (*int*)
Index in the lattices array used internally in OpenMC.
**/geometry/lattices/lattice <uid>/name** (*char[]*)
Name of the lattice.
**/geometry/lattices/lattice <uid>/type** (*char[]*)
Type of the lattice, either 'rectangular' or 'hexagonal'.
**/geometry/lattices/lattice <uid>/pitch** (*double[]*)
Pitch of the lattice.
**/geometry/lattices/lattice <uid>/outer** (*int*)
Outer universe assigned to lattice cells outside the defined range.
**/geometry/lattices/lattice <uid>/offsets** (*int[]*)
Offsets used for distribcell tally filter.
**/geometry/lattices/lattice <uid>/universes** (*int[]*)
Three-dimensional array of universes assigned to each cell of the lattice.
**/geometry/lattices/lattice <uid>/dimension** (*int[]*)
The number of lattice cells in each direction. This dataset is present only
when the 'type' dataset is set to 'rectangular'.
**/geometry/lattices/lattice <uid>/lower_left** (*double[]*)
The coordinates of the lower-left corner of the lattice. This dataset is
present only when the 'type' dataset is set to 'rectangular'.
**/geometry/lattices/lattice <uid>/n_rings** (*int*)
Number of radial ring positions in the xy-plane. This dataset is present
only when the 'type' dataset is set to 'hexagonal'.
**/geometry/lattices/lattice <uid>/n_axial** (*int*)
Number of lattice positions along the z-axis. This dataset is present only
when the 'type' dataset is set to 'hexagonal'.
**/geometry/lattices/lattice <uid>/center** (*double[]*)
Coordinates of the center of the lattice. This dataset is present only when
the 'type' dataset is set to 'hexagonal'.
**/n_materials** (*int*)
Number of materials in the problem.
**/materials/material <uid>/index** (*int*)
Index in materials array used internally in OpenMC.
**/materials/material <uid>/name** (*char[]*)
Name of the material.
**/materials/material <uid>/atom_density** (*double[]*)
Total atom density of the material in atom/b-cm.
**/materials/material <uid>/nuclides** (*char[][]*)
Array of nuclides present in the material, e.g., 'U-235.71c'.
**/materials/material <uid>/nuclide_densities** (*double[]*)
Atom density of each nuclide.
**/materials/material <uid>/sab_names** (*char[][]*)
Names of S(:math:`\alpha`,:math:`\beta`) tables assigned to the material.
**/tallies/tally <uid>/name** (*char[]*)
Name of the tally.
The current version of the summary file format is 5.0.
**/**
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
- **version** (*int[2]*) -- Major and minor version of the summary
file format.
- **openmc_version** (*int[3]*) -- Major, minor, and release
version number for OpenMC.
- **git_sha1** (*char[40]*) -- Git commit SHA-1 hash.
- **date_and_time** (*char[]*) -- Date and time the summary was
written.
**/geometry/**
:Attributes: - **n_cells** (*int*) -- Number of cells in the problem.
- **n_surfaces** (*int*) -- Number of surfaces in the problem.
- **n_universes** (*int*) -- Number of unique universes in the
problem.
- **n_lattices** (*int*) -- Number of lattices in the problem.
**/geometry/cells/cell <uid>/**
:Datasets: - **name** (*char[]*) -- User-defined name of the cell.
- **universe** (*int*) -- Universe assigned to the cell. If none is
specified, the default universe (0) is assigned.
- **fill_type** (*char[]*) -- Type of fill for the cell. Can be
'material', 'universe', or 'lattice'.
- **material** (*int* or *int[]*) -- Unique ID of the material(s)
assigned to the cell. This dataset is present only if fill_type is
set to 'normal'. The value '-1' signifies void material. The data
is an array if the cell uses distributed materials, otherwise it is
a scalar.
- **temperature** (*double[]*) -- Temperature of the cell in Kelvin.
- **offset** (*int[]*) -- Offsets used for distribcell tally
filter. This dataset is present only if fill_type is set to
'universe'.
- **translation** (*double[3]*) -- Translation applied to the fill
universe. This dataset is present only if fill_type is set to
'universe'.
- **rotation** (*double[3]*) -- Angles in degrees about the x-, y-,
and z-axes for which the fill universe should be rotated. This
dataset is present only if fill_type is set to 'universe'.
- **lattice** (*int*) -- Unique ID of the lattice which fills the
cell. Only present if fill_type is set to 'lattice'.
- **region** (*char[]*) -- Region specification for the cell.
- **distribcell_index** (*int*) -- Index of this cell in distribcell
arrays. Only present if this cell is listed in a distribcell filter
or if it uses distributed materials.
- **paths** (*char[][]*) -- The paths traversed through the CSG tree
to reach each distribcell instance. This consists of the integer
IDs for each universe, cell and lattice delimited by '->'. Each
lattice cell is specified by its (x,y) or (x,y,z) indices. Only
present if this cell is listed in a distribcell filter or if it
uses distributed materials.
**/geometry/surfaces/surface <uid>/**
:Datasets: - **name** (*char[]*) -- Name of the surface.
- **type** (*char[]*) -- Type of the surface. Can be 'x-plane',
'y-plane', 'z-plane', 'plane', 'x-cylinder', 'y-cylinder',
'z-cylinder', 'sphere', 'x-cone', 'y-cone', 'z-cone', or 'quadric'.
- **coefficients** (*double[]*) -- Array of coefficients that define
the surface. See :ref:`surface_element` for what coefficients are
defined for each surface type.
- **boundary_condition** (*char[]*) -- Boundary condition applied to
the surface. Can be 'transmission', 'vacuum', 'reflective', or
'periodic'.
**/geometry/universes/universe <uid>/**
:Datasets:
- **cells** (*int[]*) -- Array of unique IDs of cells that appear in
the universe.
**/geometry/lattices/lattice <uid>/**
:Datasets: - **name** (*char[]*) -- Name of the lattice.
- **type** (*char[]*) -- Type of the lattice, either 'rectangular' or
'hexagonal'.
- **pitch** (*double[]*) -- Pitch of the lattice in centimeters.
- **outer** (*int*) -- Outer universe assigned to lattice cells
outside the defined range.
- **offsets** (*int[]*) -- Offsets used for distribcell tally filter.
- **universes** (*int[][][]*) -- Three-dimensional array of universes
assigned to each cell of the lattice.
- **dimension** (*int[]*) -- The number of lattice cells in each
direction. This dataset is present only when the 'type' dataset is
set to 'rectangular'.
- **lower_left** (*double[]*) -- The coordinates of the lower-left
corner of the lattice. This dataset is present only when the 'type'
dataset is set to 'rectangular'.
- **n_rings** (*int*) -- Number of radial ring positions in the
xy-plane. This dataset is present only when the 'type' dataset is
set to 'hexagonal'.
- **n_axial** (*int*) -- Number of lattice positions along the
z-axis. This dataset is present only when the 'type' dataset is set
to 'hexagonal'.
- **center** (*double[]*) -- Coordinates of the center of the
lattice. This dataset is present only when the 'type' dataset is
set to 'hexagonal'.
**/materials/**
:Attributes: - **n_materials** (*int*) -- Number of materials in the problem.
**/materials/material <uid>/**
:Datasets: - **name** (*char[]*) -- Name of the material.
- **atom_density** (*double[]*) -- Total atom density of the material
in atom/b-cm.
- **nuclides** (*char[][]*) -- Array of nuclides present in the
material, e.g., 'U235'.
- **nuclide_densities** (*double[]*) -- Atom density of each nuclide.
- **sab_names** (*char[][]*) -- Names of
S(:math:`\alpha,\beta`) tables assigned to the material.
**/nuclides/**
:Attributes: - **n_nuclides** (*int*) -- Number of nuclides in the problem.
:Datasets: - **names** (*char[][]*) -- Names of nuclides.
- **awrs** (*float[]*) -- Atomic weight ratio of each nuclide.
**/tallies/tally <uid>/**
:Datasets: - **name** (*char[]*) -- Name of the tally.

33
docs/source/io_formats/track.rst
View file

@ -4,27 +4,18 @@
Track File Format
=================
The current revision of the particle track file format is 1.
The current revision of the particle track file format is 2.0.
**/filetype** (*char[]*)
**/**
String indicating the type of file.
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
- **version** (*int[2]*) -- Major and minor version of the track
file format.
- **n_particles** (*int*) -- Number of particles for which tracks
are recorded.
- **n_coords** (*int[]*) -- Number of coordinates for each
particle.
**/revision** (*int*)
Revision of the track file format. Any time a change is made in the format,
this integer is incremented.
**/n_particles** (*int*)
Number of particles for which tracks are recorded.
**/n_coords** (*int[]*)
Number of coordinates for each particle.
*do i = 1, n_particles*
**/coordinates_i** (*double[][3]*)
(x,y,z) coordinates for the *i*-th particle.
:Datasets:
- **coordinates_<i>** (*double[][3]*) -- (x,y,z) coordinates for the
*i*-th particle.

18
docs/source/io_formats/volume.rst
View file

@ -4,19 +4,31 @@
Volume File Format
==================
The current version of the volume file format is 1.0.
**/**
:Attributes: - **samples** (*int*) -- Number of samples
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
- **version** (*int[2]*) -- Major and minor version of the summary
file format.
- **openmc_version** (*int[3]*) -- Major, minor, and release
version number for OpenMC.
- **git_sha1** (*char[40]*) -- Git commit SHA-1 hash.
- **date_and_time** (*char[]*) -- Date and time the summary was
written.
- **domain_type** (*char[]*) -- The type of domain for which
volumes are calculated, either 'cell', 'material', or 'universe'.
- **samples** (*int*) -- Number of samples
- **lower_left** (*double[3]*) -- Lower-left coordinates of
bounding box
- **upper_right** (*double[3]*) -- Upper-right coordinates of
bounding box
**/cell_<id>/**
**/domain_<id>/**
:Datasets: - **volume** (*double[2]*) -- Calculated volume and its uncertainty
in cubic centimeters
- **nuclides** (*char[][]*) -- Names of nuclides identified in the
cell
domain
- **atoms** (*double[][2]*) -- Total number of atoms of each nuclide
and its uncertainty

34
docs/source/io_formats/voxel.rst
View file

@ -4,22 +4,24 @@
Voxel Plot File Format
======================
**/filetype** (*char[]*)
The current version of the voxel file format is 1.0.
String indicating the type of file.
**/**
**/num_voxels** (*int[3]*)
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
- **version** (*int[2]*) -- Major and minor version of the voxel
file format.
- **openmc_version** (*int[3]*) -- Major, minor, and release
version number for OpenMC.
- **git_sha1** (*char[40]*) -- Git commit SHA-1 hash.
- **date_and_time** (*char[]*) -- Date and time the summary was
written.
- **num_voxels** (*int[3]*) -- Number of voxels in the x-, y-, and
z- directions.
- **voxel_width** (*double[3]*) -- Width of a voxel in centimeters.
- **lower_left** (*double[3]*) -- Cartesian coordinates of the
lower-left corner of the plot.
Number of voxels in the x-, y-, and z- directions.
**/voxel_width** (*double[3]*)
Width of a voxel in centimeters.
**/lower_left** (*double[3]*)
Cartesian coordinates of the lower-left corner of the plot.
**/data** (*int[][][]*)
Data for each voxel that represents a material or cell ID.
:Datasets:
- **data** (*int[][][]*) -- Data for each voxel that represents a
material or cell ID.

227
docs/source/pythonapi/examples/mdgxs-part-i.ipynb
View file

File diff suppressed because one or more lines are too long

56
docs/source/pythonapi/examples/mgxs-part-i.ipynb
View file

@ -369,21 +369,21 @@
"\n",
"* `TotalXS`\n",
"* `TransportXS`\n",
"* `NuTransportXS`\n",
"* `AbsorptionXS`\n",
"* `CaptureXS`\n",
"* `FissionXS`\n",
"* `NuFissionXS`\n",
"* `KappaFissionXS`\n",
"* `ScatterXS`\n",
"* `NuScatterXS`\n",
"* `ScatterMatrixXS`\n",
"* `NuScatterMatrixXS`\n",
"* `Chi`\n",
"* `ChiPrompt`\n",
"* `InverseVelocity`\n",
"* `PromptNuFissionXS`\n",
"\n",
"Of course, we are aware that the fission cross section (`FissionXS`) can sometimes be paired with the fission neutron multiplication to become $\\nu\\sigma_f$. This can be accomodated in to the `FissionXS` class by setting the `nu` parameter to `True` as shown below.\n",
"\n",
"Additionally, scattering reactions (like (n,2n)) can also be defined to take in to account the neutron multiplication to become $\\nu\\sigma_s$. This can be accomodated in the the transport (`TransportXS`), scattering (`ScatterXS`), and scattering-matrix (`ScatterMatrixXS`) cross sections types by setting the `nu` parameter to `True` as shown below.\n",
"\n",
"These classes provide us with an interface to generate the tally inputs as well as perform post-processing of OpenMC's tally data to compute the respective multi-group cross sections. In this case, let's create the multi-group total, absorption and scattering cross sections with our 2-group structure."
]
},
@ -398,7 +398,11 @@
"# Instantiate a few different sections\n",
"total = mgxs.TotalXS(domain=cell, groups=groups)\n",
"absorption = mgxs.AbsorptionXS(domain=cell, groups=groups)\n",
"scattering = mgxs.ScatterXS(domain=cell, groups=groups)"
"scattering = mgxs.ScatterXS(domain=cell, groups=groups)\n",
"\n",
"# Note that if we wanted to incorporate neutron multiplication in the\n",
"# scattering cross section we would write the previous line as:\n",
"# scattering = mgxs.ScatterXS(domain=cell, groups=groups, nu=True)"
]
},
{
@ -520,8 +524,8 @@
" Copyright | 2011-2017 Massachusetts Institute of Technology\n",
" License | http://openmc.readthedocs.io/en/latest/license.html\n",
" Version | 0.8.0\n",
" Git SHA1 | 54b65c8bda6af5788bd762b8cf9855d1a8008238\n",
" Date/Time | 2017-02-12 13:36:24\n",
" Git SHA1 | 60a1f157dae88b62e1865a5fe3efd7ef0773a068\n",
" Date/Time | 2017-02-25 14:26:54\n",
" OpenMP Threads | 8\n",
"\n",
" ===========================================================================\n",
@ -592,7 +596,7 @@
" 42/1 1.13779 1.16177 +/- 0.00531\n",
" 43/1 1.15066 1.16143 +/- 0.00516\n",
" 44/1 1.12174 1.16026 +/- 0.00514\n",
" 45/1 1.17479 1.16068 +/- 0.00501\n",
" 45/1 1.17478 1.16068 +/- 0.00501\n",
" 46/1 1.14146 1.16014 +/- 0.00489\n",
" 47/1 1.20464 1.16135 +/- 0.00491\n",
" 48/1 1.15119 1.16108 +/- 0.00479\n",
@ -607,20 +611,20 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 4.1327E-01 seconds\n",
" Reading cross sections = 3.2638E-01 seconds\n",
" Total time in simulation = 2.2324E+00 seconds\n",
" Time in transport only = 2.1226E+00 seconds\n",
" Time in inactive batches = 3.0650E-01 seconds\n",
" Time in active batches = 1.9259E+00 seconds\n",
" Time synchronizing fission bank = 2.7640E-03 seconds\n",
" Sampling source sites = 2.0198E-03 seconds\n",
" SEND/RECV source sites = 7.0929E-04 seconds\n",
" Time accumulating tallies = 4.5355E-05 seconds\n",
" Total time for finalization = 4.1885E-04 seconds\n",
" Total time elapsed = 2.6534E+00 seconds\n",
" Calculation Rate (inactive) = 81567.1 neutrons/second\n",
" Calculation Rate (active) = 51923.2 neutrons/second\n",
" Total time for initialization = 3.5070E-01 seconds\n",
" Reading cross sections = 2.4151E-01 seconds\n",
" Total time in simulation = 2.3276E+00 seconds\n",
" Time in transport only = 2.2350E+00 seconds\n",
" Time in inactive batches = 2.5677E-01 seconds\n",
" Time in active batches = 2.0708E+00 seconds\n",
" Time synchronizing fission bank = 2.7683E-03 seconds\n",
" Sampling source sites = 2.0233E-03 seconds\n",
" SEND/RECV source sites = 7.1007E-04 seconds\n",
" Time accumulating tallies = 5.0753E-05 seconds\n",
" Total time for finalization = 3.8695E-04 seconds\n",
" Total time elapsed = 2.6857E+00 seconds\n",
" Calculation Rate (inactive) = 97364.6 neutrons/second\n",
" Calculation Rate (active) = 48290.8 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
@ -901,7 +905,7 @@
" <td>6.250000e-01</td>\n",
" <td>total</td>\n",
" <td>(((total / flux) - (absorption / flux)) - (sca...</td>\n",
" <td>-2.664535e-15</td>\n",
" <td>-5.551115e-15</td>\n",
" <td>0.011292</td>\n",
" </tr>\n",
" <tr>\n",
@ -911,7 +915,7 @@
" <td>2.000000e+07</td>\n",
" <td>total</td>\n",
" <td>(((total / flux) - (absorption / flux)) - (sca...</td>\n",
" <td>-3.330669e-16</td>\n",
" <td>-1.110223e-16</td>\n",
" <td>0.002570</td>\n",
" </tr>\n",
" </tbody>\n",
@ -924,8 +928,8 @@
"1 1 6.25e-01 2.00e+07 total \n",
"\n",
" score mean std. dev. \n",
"0 (((total / flux) - (absorption / flux)) - (sca... -2.66e-15 1.13e-02 \n",
"1 (((total / flux) - (absorption / flux)) - (sca... -3.33e-16 2.57e-03 "
"0 (((total / flux) - (absorption / flux)) - (sca... -5.55e-15 1.13e-02 \n",
"1 (((total / flux) - (absorption / flux)) - (sca... -1.11e-16 2.57e-03 "
]
},
"execution_count": 22,

1532
docs/source/pythonapi/examples/mgxs-part-ii.ipynb
View file

File diff suppressed because one or more lines are too long

463
docs/source/pythonapi/examples/mgxs-part-iii.ipynb
View file

File diff suppressed because one or more lines are too long

344
docs/source/pythonapi/examples/mgxs-part-iv.ipynb
View file

File diff suppressed because one or more lines are too long

8
docs/source/pythonapi/index.rst
View file

@ -115,15 +115,17 @@ Many of the above classes are derived from several abstract classes:
openmc.Region
openmc.Lattice
One function is also available to create a hexagonal region defined by the
intersection of six surface half-spaces.
Two helper function are also available to create rectangular and hexagonal
prisms defined by the intersection of four and six surface half-spaces,
respectively.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.make_hexagon_region
openmc.get_hexagonal_prism
openmc.get_rectangular_prism
Constructing Tallies
--------------------

160
docs/source/usersguide/input.rst
View file

@ -94,6 +94,18 @@ Settings Specification -- settings.xml
All simulation parameters and miscellaneous options are specified in the
settings.xml file.
``<batches>`` Element
---------------------
The ``<batches>`` element indicates the total number of batches to execute,
where each batch corresponds to a tally realization. In a fixed source
calculation, each batch consists of a number of source particles. In an
eigenvalue calculation, each batch consists of one or many fission source
iterations (generations), where each generation itself consists of a number of
source neutrons.
*Default*: None
``<confidence_intervals>`` Element
----------------------------------
@ -132,67 +144,6 @@ you care. This element has the following attributes/sub-elements:
*Default*: 0.0
.. _eigenvalue:
``<eigenvalue>`` Element
------------------------
The ``<eigenvalue>`` element indicates that a :math:`k`-eigenvalue calculation
should be performed. It has the following attributes/sub-elements:
:batches:
The total number of batches, where each batch corresponds to multiple
fission source iterations. Batching is done to eliminate correlation between
realizations of random variables.
*Default*: None
:generations_per_batch:
The number of total fission source iterations per batch.
*Default*: 1
:inactive:
The number of inactive batches. In general, the starting cycles in a
criticality calculation can not be used to contribute to tallies since the
fission source distribution and eigenvalue are generally not converged
immediately.
*Default*: None
:particles:
The number of neutrons to simulate per fission source iteration.
*Default*: None
:keff_trigger:
This tag specifies a precision trigger on the combined :math:`k_{eff}`. The
trigger is a convergence criterion on the uncertainty of the estimated
eigenvalue. It has the following attributes/sub-elements:
:type:
The type of precision trigger. Accepted options are "variance", "std_dev",
and "rel_err".
:variance:
Variance of the batch mean :math:`\sigma^2`
:std_dev:
Standard deviation of the batch mean :math:`\sigma`
:rel_err:
Relative error of the batch mean :math:`\frac{\sigma}{\mu}`
*Default*: None
:threshold:
The precision trigger's convergence criterion for the
combined :math:`k_{eff}`.
*Default*: None
.. note:: See section on the :ref:`trigger` for more information.
``<energy_grid>`` Element
-------------------------
@ -247,23 +198,58 @@ problem. It has the following attributes/sub-elements:
*Default*: None
``<fixed_source>`` Element
``<generations_per_batch>`` Element
-----------------------------------
The ``<generations_per_batch>`` element indicates the number of total fission
source iterations per batch for an eigenvalue calculation. This element is
ignored for all run modes other than "eigenvalue".
*Default*: 1
``<inactive>`` Element
----------------------
The ``<inactive>`` element indicates the number of inactive batches used in a
k-eigenvalue calculation. In general, the starting fission source iterations in
an eigenvalue calculation can not be used to contribute to tallies since the
fission source distribution and eigenvalue are generally not converged
immediately. This element is ignored for all run modes other than "eigenvalue".
*Default*: 0
``<keff_trigger>`` Element
--------------------------
The ``<fixed_source>`` element indicates that a fixed source calculation should
be performed. It has the following attributes/sub-elements:
The ``<keff_trigger>`` element (ignored for all run modes other than
"eigenvalue".) specifies a precision trigger on the combined
:math:`k_{eff}`. The trigger is a convergence criterion on the uncertainty of
the estimated eigenvalue. It has the following attributes/sub-elements:
:batches:
The total number of batches. For fixed source calculations, each batch
represents a realization of random variables for tallies.
:type:
The type of precision trigger. Accepted options are "variance", "std_dev",
and "rel_err".
:variance:
Variance of the batch mean :math:`\sigma^2`
:std_dev:
Standard deviation of the batch mean :math:`\sigma`
:rel_err:
Relative error of the batch mean :math:`\frac{\sigma}{\mu}`
*Default*: None
:particles:
The number of particles to simulate per batch.
:threshold:
The precision trigger's convergence criterion for the
combined :math:`k_{eff}`.
*Default*: None
.. note:: See section on the :ref:`trigger` for more information.
``<log_grid_bins>`` Element
---------------------------
@ -336,6 +322,15 @@ will abort.
*Default*: Current working directory
``<particles>`` Element
-----------------------
This element indicates the number of neutrons to simulate per fission source
iteration when a k-eigenvalue calculation is performed or the number of neutrons
per batch for a fixed source simulation.
*Default*: None
``<ptables>`` Element
---------------------
@ -408,7 +403,16 @@ The ``<run_cmfd>`` element indicates whether or not CMFD acceleration should be
turned on or off. This element has no attributes or sub-elements and can be set
to either "false" or "true".
*Defualt*: false
*Default*: false
``<run_mode>`` Element
----------------------
The ``<run_mode>`` element indicates which run mode should be used when OpenMC
is executed. This element has no attributes or sub-elements and can be set to
"eigenvalue", "fixed source", "plot", "volume", or "particle restart".
*Default*: None
``<seed>`` Element
------------------
@ -774,13 +778,13 @@ number, and particle number, respectively.
-------------------------
OpenMC includes tally precision triggers which allow the user to define
uncertainty thresholds on :math:`k_{eff}` in the ``<eigenvalue>`` subelement of
``settings.xml``, and/or tallies in ``tallies.xml``. When using triggers,
uncertainty thresholds on :math:`k_{eff}` in the ``<keff_trigger>`` subelement
of ``settings.xml``, and/or tallies in ``tallies.xml``. When using triggers,
OpenMC will run until it completes as many batches as defined by ``<batches>``.
At this point, the uncertainties on all tallied values are computed and
compared with their corresponding trigger thresholds. If any triggers have not
been met, OpenMC will continue until either all trigger thresholds have been
satisfied or ``<max_batches>`` has been reached.
At this point, the uncertainties on all tallied values are computed and compared
with their corresponding trigger thresholds. If any triggers have not been met,
OpenMC will continue until either all trigger thresholds have been satisfied or
``<max_batches>`` has been reached.
The ``<trigger>`` element provides an active "toggle switch" for tally
precision trigger(s), the maximum number of batches and the batch interval. It
@ -793,8 +797,8 @@ has the following attributes/sub-elements:
:max_batches:
This describes the maximum number of batches allowed when using trigger(s).
.. note:: When max_batches is set, the number of ``batches`` shown in
``<eigenvalue>`` element represents minimum number of batches to
.. note:: When max_batches is set, the number of ``batches`` shown in the
``<batches>`` element represents minimum number of batches to
simulate when using the trigger(s).
:batch_interval:

10
examples/xml/basic/settings.xml
View file

@ -1,12 +1,10 @@
<?xml version="1.0"?>
<settings>
<!-- Parameters for k-eigenvalue calculation -->
<eigenvalue>
<batches>15</batches>
<inactive>5</inactive>
<particles>10000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>15</batches>
<inactive>5</inactive>
<particles>10000</particles>
<!-- Starting source -->
<source>

9
examples/xml/boxes/settings.xml
View file

@ -2,11 +2,10 @@
<settings>
<!-- Parameters for k-eigenvalue calculation -->
<eigenvalue>
<batches>15</batches>
<inactive>5</inactive>
<particles>10000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>15</batches>
<inactive>5</inactive>
<particles>10000</particles>
<!-- Starting source -->
<source>

9
examples/xml/lattice/nested/settings.xml
View file

@ -2,11 +2,10 @@
<settings>
<!-- Parameters for k-eigenvalue calculation -->
<eigenvalue>
<batches>20</batches>
<inactive>10</inactive>
<particles>10000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>20</batches>
<inactive>10</inactive>
<particles>10000</particles>
<!-- Starting source -->
<source>

9
examples/xml/lattice/simple/settings.xml
View file

@ -2,11 +2,10 @@
<settings>
<!-- Parameters for k-eigenvalue calculation -->
<eigenvalue>
<batches>20</batches>
<inactive>10</inactive>
<particles>10000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>20</batches>
<inactive>10</inactive>
<particles>10000</particles>
<!-- Starting source -->
<source>

11
examples/xml/pincell/settings.xml
View file

@ -2,11 +2,10 @@
<settings>
<!-- Define how many particles to run and for how many batches -->
<eigenvalue>
<batches>100</batches>
<inactive>10</inactive>
<particles>1000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>100</batches>
<inactive>10</inactive>
<particles>1000</particles>
<!-- The starting source is a uniform distribution over the entire pin
cell. Note that since this is effectively a 2D model, the z coordinates
@ -29,4 +28,4 @@
<dimension>10 10 1</dimension>
</entropy>
</settings>
</settings>

18
examples/xml/pincell_multigroup/settings.xml
View file

@ -1,14 +1,10 @@
<?xml version='1.0' encoding='utf-8'?>
<settings>
<eigenvalue>
<particles>1000</particles>
<batches>100</batches>
<inactive>10</inactive>
</eigenvalue>
<source strength="1.0">
<space type="box">
<parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
<run_mode>eigenvalue</run_mode>
<source strength="1.0">
<space type="box">
<parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
</settings>

9
examples/xml/reflective/settings.xml
View file

@ -2,11 +2,10 @@
<settings>
<!-- Parameters for k-eigenvalue calculation -->
<eigenvalue>
<batches>500</batches>
<inactive>10</inactive>
<particles>10000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>500</batches>
<inactive>10</inactive>
<particles>10000</particles>
<!-- Starting source -->
<source>

25
openmc/cell.py
View file

@ -433,7 +433,7 @@ class Cell(object):
Returns
-------
nuclides : dict
nuclides : collections.OrderedDict
Dictionary whose keys are nuclide names and values are 2-tuples of
(nuclide, density)
@ -467,7 +467,7 @@ class Cell(object):
Returns
-------
cells : dict
cells : collections.orderedDict
Dictionary whose keys are cell IDs and values are :class:`Cell`
instances
@ -485,20 +485,23 @@ class Cell(object):
Returns
-------
materials : dict
materials : collections.OrderedDict
Dictionary whose keys are material IDs and values are
:class:`Material` instances
"""
materials = OrderedDict()
if self.fill_type == 'material':
materials[self.fill.id] = self.fill
# Append all Cells in each Cell in the Universe to the dictionary
cells = self.get_all_cells()
for cell in cells.values():
materials.update(cell.get_all_materials())
elif self.fill_type == 'distribmat':
for m in self.fill:
if m is not None:
materials[m.id] = m
else:
# Append all Cells in each Cell in the Universe to the dictionary
cells = self.get_all_cells()
for cell in cells.values():
materials.update(cell.get_all_materials())
return materials
@ -508,7 +511,7 @@ class Cell(object):
Returns
-------
universes : dict
universes : collections.OrderedDict
Dictionary whose keys are universe IDs and values are
:class:`Universe` instances
@ -563,7 +566,7 @@ class Cell(object):
if isinstance(node, Halfspace):
path = "./surface[@id='{}']".format(node.surface.id)
if xml_element.find(path) is None:
xml_element.append(node.surface.create_xml_subelement())
xml_element.append(node.surface.to_xml_element())
elif isinstance(node, Complement):
create_surface_elements(node.node, element)
else:

34
openmc/checkvalue.py
View file

@ -239,6 +239,40 @@ def check_greater_than(name, value, minimum, equality=False):
raise ValueError(msg)
def check_filetype_version(obj, expected_type, expected_version):
"""Check filetype and version of an HDF5 file.
Parameters
----------
obj : h5py.File
HDF5 file to check
expected_type
Expected file type, e.g. 'statepoint'
expected_version
Expected major version number.
"""
try:
this_filetype = obj.attrs['filetype'].decode()
this_version = obj.attrs['version']
# Check filetype
if this_filetype != expected_type:
raise IOError('{} is not a {} file.'.format(
obj.filename, expected_type))
# Check version
if this_version[0] != expected_version:
raise IOError('{} file has a version of {} which is not '
'consistent with the version expected by OpenMC, {}'
.format(this_filetype, '.'.join(this_version),
expected_version))
except AttributeError:
raise IOError('Could not read {} file. This most likely means the {} '
'file was produced by a different version of OpenMC than '
'the one you are using.'.format(expected_type))
class CheckedList(list):
"""A list for which each element is type-checked as it's added

46
openmc/executor.py
View file

@ -5,21 +5,33 @@ import sys
from six import string_types
_summary_indicator = "TIMING STATISTICS"
def _run(command, output, cwd):
# Launch a subprocess
p = subprocess.Popen(command, shell=True, cwd=cwd, stdout=subprocess.PIPE,
stderr=subprocess.STDOUT, universal_newlines=True)
storage_flag = False
# Capture and re-print OpenMC output in real-time
while True:
# If OpenMC is finished, break loop
line = p.stdout.readline()
if not line and p.poll() != None:
if not line and p.poll() is not None:
break
# If user requested output, print to screen
if output:
if output == 'full':
# If user requested output, print to screen
print(line, end='')
elif output == 'summary' and _summary_indicator in line:
# If they requested a summary, look for the start of the summary
storage_flag = True
if storage_flag:
# If a summary is requested, and we have reached the summary,
# then print it
print(line, end='')
# Return the returncode (integer, zero if no problems encountered)
@ -44,8 +56,8 @@ def plot_geometry(output=True, openmc_exec='openmc', cwd='.'):
def run(particles=None, threads=None, geometry_debug=False,
restart_file=None, tracks=False, mpi_procs=1, output=True,
openmc_exec='openmc', mpi_exec='mpiexec', cwd='.'):
restart_file=None, tracks=False, output='full', cwd='.',
openmc_exec='openmc', mpi_args=None):
"""Run an OpenMC simulation.
Parameters
@ -56,23 +68,25 @@ def run(particles=None, threads=None, geometry_debug=False,
Number of OpenMP threads. If OpenMC is compiled with OpenMP threading
enabled, the default is implementation-dependent but is usually equal to
the number of hardware threads available (or a value set by the
OMP_NUM_THREADS environment variable).
:envvar:`OMP_NUM_THREADS` environment variable).
geometry_debug : bool, optional
Turn on geometry debugging during simulation. Defaults to False.
restart_file : str, optional
Path to restart file to use
tracks : bool, optional
Write tracks for all particles. Defaults to False.
mpi_procs : int, optional
Number of MPI processes.
output : bool, optional
Capture OpenMC output from standard out. Defaults to True.
output : {"full", "summary", "none", False}, optional
Degree of OpenMC output captured from standard out. "full" prints all
output; "summary" prints only the results summary, and "none" or False
does not show the output. Defaults to "full".
cwd : str, optional
Path to working directory to run in. Defaults to the current working
directory.
openmc_exec : str, optional
Path to OpenMC executable. Defaults to 'openmc'.
mpi_exec : str, optional
MPI execute command. Defaults to 'mpiexec'.
cwd : str, optional
Path to working directory to run in. Defaults to the current working directory.
mpi_args : list of str, optional
MPI execute command and any additional MPI arguments to pass,
e.g. ['mpiexec', '-n', '8'].
"""
@ -94,8 +108,8 @@ def run(particles=None, threads=None, geometry_debug=False,
if tracks:
post_args += '-t'
if isinstance(mpi_procs, Integral) and mpi_procs > 1:
pre_args += '{} -n {} '.format(mpi_exec, mpi_procs)
if mpi_args is not None:
pre_args = ' '.join(mpi_args) + ' '
command = pre_args + openmc_exec + ' ' + post_args

108
openmc/geometry.py
View file

@ -151,102 +151,83 @@ class Geometry(object):
return offset
def get_all_cells(self):
"""Return all cells defined
"""Return all cells in the geometry.
Returns
-------
list of openmc.Cell
Cells in the geometry
collections.OrderedDict
Dictionary mapping cell IDs to :class:`openmc.Cell` instances
"""
all_cells = self.root_universe.get_all_cells()
cells = list(set(all_cells.values()))
cells.sort(key=lambda x: x.id)
return cells
return self.root_universe.get_all_cells()
def get_all_universes(self):
"""Return all universes defined
"""Return all universes in the geometry.
Returns
-------
list of openmc.Universe
Universes in the geometry
collections.OrderedDict
Dictionary mapping universe IDs to :class:`openmc.Universe`
instances
"""
all_universes = self._root_universe.get_all_universes()
universes = list(set(all_universes.values()))
universes.sort(key=lambda x: x.id)
universes = OrderedDict()
universes[self.root_universe.id] = self.root_universe
universes.update(self.root_universe.get_all_universes())
return universes
def get_all_materials(self):
"""Return all materials assigned to a cell
"""Return all materials within the geometry.
Returns
-------
list of openmc.Material
Materials in the geometry
collections.OrderedDict
Dictionary mapping material IDs to :class:`openmc.Material`
instances
"""
material_cells = self.get_all_material_cells()
materials = []
for cell in material_cells:
if cell.fill_type == 'distribmat':
for m in cell.fill:
if m is not None and m not in materials:
materials.append(m)
elif cell.fill_type == 'material':
if cell.fill not in materials:
materials.append(cell.fill)
materials.sort(key=lambda x: x.id)
return materials
return self.root_universe.get_all_materials()
def get_all_material_cells(self):
"""Return all cells filled by a material
Returns
-------
list of openmc.Cell
Cells filled by Materials in the geometry
collections.OrderedDict
Dictionary mapping cell IDs to :class:`openmc.Cell` instances that
are filled with materials or distributed materials.
"""
material_cells = OrderedDict()
all_cells = self.get_all_cells()
material_cells = []
for cell in all_cells:
for cell in self.get_all_cells().values():
if cell.fill_type in ('material', 'distribmat'):
if cell not in material_cells:
material_cells.append(cell)
material_cells[cell.id] = cell
material_cells.sort(key=lambda x: x.id)
return material_cells
def get_all_material_universes(self):
"""Return all universes composed of at least one non-fill cell
"""Return all universes having at least one material-filled cell.
This method can be used to find universes that have at least one cell
that is filled with a material or is void.
Returns
-------
list of openmc.Universe
Universes with non-fill cells
collections.OrderedDict
Dictionary mapping universe IDs to :class:`openmc.Universe`
instances with at least one material-filled cell
"""
material_universes = OrderedDict()
all_universes = self.get_all_universes()
material_universes = []
for universe in all_universes:
cells = universe.cells
for cell in cells:
for universe in self.get_all_universes():
for cell in universe.cells:
if cell.fill_type in ('material', 'distribmat', 'void'):
if universe not in material_universes:
material_universes.append(universe)
material_universes[universe.id] = universe
material_universes.sort(key=lambda x: x.id)
return material_universes
def get_all_lattices(self):
@ -254,20 +235,17 @@ class Geometry(object):
Returns
-------
list of openmc.Lattice
Lattices in the geometry
collections.OrderedDict
Dictionary mapping lattice IDs to :class:`openmc.Lattice` instances
"""
lattices = OrderedDict()
cells = self.get_all_cells()
lattices = []
for cell in cells:
for cell in self.get_all_cells():
if cell.fill_type == 'lattice':
if cell.fill not in lattices:
lattices.append(cell.fill)
lattices[cell.fill.id] = cell.fill
lattices.sort(key=lambda x: x.id)
return lattices
def get_materials_by_name(self, name, case_sensitive=False, matching=False):
@ -293,7 +271,7 @@ class Geometry(object):
if not case_sensitive:
name = name.lower()
all_materials = self.get_all_materials()
all_materials = self.get_all_materials().values()
materials = set()
for material in all_materials:
@ -333,7 +311,7 @@ class Geometry(object):
if not case_sensitive:
name = name.lower()
all_cells = self.get_all_cells()
all_cells = self.get_all_cells().values()
cells = set()
for cell in all_cells:
@ -373,7 +351,7 @@ class Geometry(object):
if not case_sensitive:
name = name.lower()
all_cells = self.get_all_cells()
all_cells = self.get_all_cells().values()
cells = set()
for cell in all_cells:
@ -413,7 +391,7 @@ class Geometry(object):
if not case_sensitive:
name = name.lower()
all_universes = self.get_all_universes()
all_universes = self.get_all_universes().values()
universes = set()
for universe in all_universes:
@ -453,7 +431,7 @@ class Geometry(object):
if not case_sensitive:
name = name.lower()
all_lattices = self.get_all_lattices()
all_lattices = self.get_all_lattices().values()
lattices = set()
for lattice in all_lattices:

151
openmc/lattice.py
View file

@ -110,6 +110,157 @@ class Lattice(object):
cv.check_type('outer universe', outer, openmc.Universe)
self._outer = outer
@staticmethod
def from_hdf5(group, universes):
"""Create lattice from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
universes : dict
Dictionary mapping universe IDs to instances of
:class:`openmc.Universe`.
Returns
-------
openmc.Lattice
Instance of lattice subclass
"""
lattice_id = int(group.name.split('/')[-1].lstrip('lattice '))
name = group['name'].value.decode()
lattice_type = group['type'].value.decode()
if 'offsets' in group:
offsets = group['offsets'][...]
else:
offsets = None
if lattice_type == 'rectangular':
dimension = group['dimension'][...]
lower_left = group['lower_left'][...]
pitch = group['pitch'][...]
outer = group['outer'].value
universe_ids = group['universes'][...]
# Create the Lattice
lattice = openmc.RectLattice(lattice_id, name)
lattice.lower_left = lower_left
lattice.pitch = pitch
# If the Universe specified outer the Lattice is not void
if outer >= 0:
lattice.outer = universes[outer]
# Build array of Universe pointers for the Lattice
uarray = np.empty(universe_ids.shape, dtype=openmc.Universe)
for z in range(universe_ids.shape[0]):
for y in range(universe_ids.shape[1]):
for x in range(universe_ids.shape[2]):
uarray[z, y, x] = universes[universe_ids[z, y, x]]
# Use 2D NumPy array to store lattice universes for 2D lattices
if len(dimension) == 2:
uarray = np.squeeze(uarray)
uarray = np.atleast_2d(uarray)
# Set the universes for the lattice
lattice.universes = uarray
# Set the distribcell offsets for the lattice
if offsets is not None:
lattice.offsets = offsets
elif lattice_type == 'hexagonal':
n_rings = group['n_rings'].value
n_axial = group['n_axial'].value
center = group['center'][...]
pitch = group['pitch'][...]
outer = group['outer'].value
universe_ids = group['universes'][...]
# Create the Lattice
lattice = openmc.HexLattice(lattice_id, name)
lattice.center = center
lattice.pitch = pitch
# If the Universe specified outer the Lattice is not void
if outer >= 0:
lattice.outer = universes[outer]
# Build array of Universe pointers for the Lattice. Note that
# we need to convert between the HDF5's square array of
# (x, alpha, z) to the Python API's format of a ragged nested
# list of (z, ring, theta).
uarray = []
for z in range(n_axial):
# Add a list for this axial level.
uarray.append([])
x = n_rings - 1
a = 2*n_rings - 2
for r in range(n_rings - 1, 0, -1):
# Add a list for this ring.
uarray[-1].append([])
# Climb down the top-right.
for i in range(r):
uarray[-1][-1].append(universe_ids[z, a, x])
x += 1
a -= 1
# Climb down the right.
for i in range(r):
uarray[-1][-1].append(universe_ids[z, a, x])
a -= 1
# Climb down the bottom-right.
for i in range(r):
uarray[-1][-1].append(universe_ids[z, a, x])
x -= 1
# Climb up the bottom-left.
for i in range(r):
uarray[-1][-1].append(universe_ids[z, a, x])
x -= 1
a += 1
# Climb up the left.
for i in range(r):
uarray[-1][-1].append(universe_ids[z, a, x])
a += 1
# Climb up the top-left.
for i in range(r):
uarray[-1][-1].append(universe_ids[z, a, x])
x += 1
# Move down to the next ring.
a -= 1
# Convert the ids into Universe objects.
uarray[-1][-1] = [universes[u_id]
for u_id in uarray[-1][-1]]
# Handle the degenerate center ring separately.
u_id = universe_ids[z, a, x]
uarray[-1].append([universes[u_id]])
# Add the universes to the lattice.
if len(pitch) == 2:
# Lattice is 3D
lattice.universes = uarray
else:
# Lattice is 2D; extract the only axial level
lattice.universes = uarray[0]
if offsets is not None:
lattice.offsets = offsets
return lattice
def get_unique_universes(self):
"""Determine all unique universes in the lattice

165
openmc/material.py
View file

@ -56,6 +56,9 @@ class Material(object):
Units used for `density`. Can be one of 'g/cm3', 'g/cc', 'kg/cm3',
'atom/b-cm', 'atom/cm3', 'sum', or 'macro'. The 'macro' unit only
applies in the case of a multi-group calculation.
depletable : bool
Indicate whether the material is depletable. This attribute can be used
by downstream depletion applications.
elements : list of tuple
List in which each item is a 4-tuple consisting of an
:class:`openmc.Element` instance, the percent density, the percent
@ -78,6 +81,7 @@ class Material(object):
self.temperature = temperature
self._density = None
self._density_units = ''
self._depletable = False
# A list of tuples (nuclide, percent, percent type)
self._nuclides = []
@ -127,37 +131,36 @@ class Material(object):
def __repr__(self):
string = 'Material\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
string += '{0: <16}{1}{2}\n'.format('\tTemperature', '=\t',
self._temperature)
string += '{: <16}=\t{}\n'.format('\tID', self._id)
string += '{: <16}=\t{}\n'.format('\tName', self._name)
string += '{: <16}=\t{}\n'.format('\tTemperature', self._temperature)
string += '{0: <16}{1}{2}'.format('\tDensity', '=\t', self._density)
string += ' [{0}]\n'.format(self._density_units)
string += '{: <16}=\t{}'.format('\tDensity', self._density)
string += ' [{}]\n'.format(self._density_units)
string += '{0: <16}\n'.format('\tS(a,b) Tables')
string += '{: <16}\n'.format('\tS(a,b) Tables')
for sab in self._sab:
string += '{0: <16}{1}{2}\n'.format('\tS(a,b)', '=\t', sab)
string += '{: <16}=\t{}\n'.format('\tS(a,b)', sab)
string += '{0: <16}\n'.format('\tNuclides')
string += '{: <16}\n'.format('\tNuclides')
for nuclide, percent, percent_type in self._nuclides:
string += '{0: <16}'.format('\t{0.name}'.format(nuclide))
string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type)
string += '=\t{: <12} [{}]\n'.format(percent, percent_type)
if self._macroscopic is not None:
string += '{0: <16}\n'.format('\tMacroscopic Data')
string += '{0: <16}'.format('\t{0}'.format(self._macroscopic))
string += '{: <16}\n'.format('\tMacroscopic Data')
string += '{: <16}'.format('\t{}'.format(self._macroscopic))
string += '{0: <16}\n'.format('\tElements')
string += '{: <16}\n'.format('\tElements')
for element, percent, percent_type, enr in self._elements:
string += '{0: <16}'.format('\t{0.name}'.format(element))
if enr is None:
string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type)
string += '=\t{: <12} [{}]\n'.format(percent, percent_type)
else:
string += '=\t{0: <12} [{1}] @ {2} w/o enrichment\n'\
string += '=\t{: <12} [{}] @ {} w/o enrichment\n'\
.format(percent, percent_type, enr)
return string
@ -182,6 +185,10 @@ class Material(object):
def density_units(self):
return self._density_units
@property
def depletable(self):
return self._depletable
@property
def elements(self):
return self._elements
@ -234,7 +241,7 @@ class Material(object):
@name.setter
def name(self, name):
if name is not None:
cv.check_type('name for Material ID="{0}"'.format(self._id),
cv.check_type('name for Material ID="{}"'.format(self._id),
name, string_types)
self._name = name
else:
@ -242,10 +249,59 @@ class Material(object):
@temperature.setter
def temperature(self, temperature):
cv.check_type('Temperature for Material ID="{0}"'.format(self._id),
cv.check_type('Temperature for Material ID="{}"'.format(self._id),
temperature, (Real, type(None)))
self._temperature = temperature
@depletable.setter
def depletable(self, depletable):
cv.check_type('Depletable flag for Material ID="{}"'.format(self.id),
depletable, bool)
self._depletable = depletable
@classmethod
def from_hdf5(cls, group):
"""Create material from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
Returns
-------
openmc.Material
Material instance
"""
mat_id = int(group.name.split('/')[-1].lstrip('material '))
name = group['name'].value.decode()
density = group['atom_density'].value
nuc_densities = group['nuclide_densities'][...]
nuclides = group['nuclides'].value
# Create the Material
material = cls(mat_id, name)
material.depletable = bool(group.attrs['depletable'])
# Read the names of the S(a,b) tables for this Material and add them
if 'sab_names' in group:
sab_tables = group['sab_names'].value
for sab_table in sab_tables:
name = sab_table.decode()
material.add_s_alpha_beta(name)
# Set the Material's density to atom/b-cm as used by OpenMC
material.set_density(density=density, units='atom/b-cm')
# Add all nuclides to the Material
for fullname, density in zip(nuclides, nuc_densities):
name = fullname.decode().strip()
material.add_nuclide(name, percent=density, percent_type='ao')
return material
def set_density(self, units, density=None):
"""Set the density of the material
@ -264,17 +320,17 @@ class Material(object):
if units is 'sum':
if density is not None:
msg = 'Density "{0}" for Material ID="{1}" is ignored ' \
msg = 'Density "{}" for Material ID="{}" is ignored ' \
'because the unit is "sum"'.format(density, self.id)
warnings.warn(msg)
else:
if density is None:
msg = 'Unable to set the density for Material ID="{0}" ' \
msg = 'Unable to set the density for Material ID="{}" ' \
'because a density value must be given when not using ' \
'"sum" unit'.format(self.id)
raise ValueError(msg)
cv.check_type('the density for Material ID="{0}"'.format(self.id),
cv.check_type('the density for Material ID="{}"'.format(self.id),
density, Real)
self._density = density
@ -285,8 +341,8 @@ class Material(object):
'version of openmc')
if not isinstance(filename, string_types) and filename is not None:
msg = 'Unable to add OTF material file to Material ID="{0}" with a ' \
'non-string name "{1}"'.format(self._id, filename)
msg = 'Unable to add OTF material file to Material ID="{}" with a ' \
'non-string name "{}"'.format(self._id, filename)
raise ValueError(msg)
self._distrib_otf_file = filename
@ -314,23 +370,23 @@ class Material(object):
"""
if self._macroscopic is not None:
msg = 'Unable to add a Nuclide to Material ID="{0}" as a ' \
msg = 'Unable to add a Nuclide to Material ID="{}" as a ' \
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if not isinstance(nuclide, string_types + (openmc.Nuclide,)):
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'non-Nuclide value "{1}"'.format(self._id, nuclide)
msg = 'Unable to add a Nuclide to Material ID="{}" with a ' \
'non-Nuclide value "{}"'.format(self._id, nuclide)
raise ValueError(msg)
elif not isinstance(percent, Real):
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'non-floating point value "{1}"'.format(self._id, percent)
msg = 'Unable to add a Nuclide to Material ID="{}" with a ' \
'non-floating point value "{}"'.format(self._id, percent)
raise ValueError(msg)
elif percent_type not in ['ao', 'wo', 'at/g-cm']:
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'percent type "{1}"'.format(self._id, percent_type)
msg = 'Unable to add a Nuclide to Material ID="{}" with a ' \
'percent type "{}"'.format(self._id, percent_type)
raise ValueError(msg)
if isinstance(nuclide, openmc.Nuclide):
@ -353,7 +409,7 @@ class Material(object):
"""
if not isinstance(nuclide, openmc.Nuclide):
msg = 'Unable to remove a Nuclide "{0}" in Material ID="{1}" ' \
msg = 'Unable to remove a Nuclide "{}" in Material ID="{}" ' \
'since it is not a Nuclide'.format(self._id, nuclide)
raise ValueError(msg)
@ -377,15 +433,15 @@ class Material(object):
# Ensure no nuclides, elements, or sab are added since these would be
# incompatible with macroscopics
if self._nuclides or self._elements or self._sab:
msg = 'Unable to add a Macroscopic data set to Material ID="{0}" ' \
'with a macroscopic value "{1}" as an incompatible data ' \
msg = 'Unable to add a Macroscopic data set to Material ID="{}" ' \
'with a macroscopic value "{}" as an incompatible data ' \
'member (i.e., nuclide, element, or S(a,b) table) ' \
'has already been added'.format(self._id, macroscopic)
raise ValueError(msg)
if not isinstance(macroscopic, string_types + (openmc.Macroscopic,)):
msg = 'Unable to add a Macroscopic to Material ID="{0}" with a ' \
'non-Macroscopic value "{1}"'.format(self._id, macroscopic)
msg = 'Unable to add a Macroscopic to Material ID="{}" with a ' \
'non-Macroscopic value "{}"'.format(self._id, macroscopic)
raise ValueError(msg)
if isinstance(macroscopic, openmc.Macroscopic):
@ -398,7 +454,7 @@ class Material(object):
if self._macroscopic is None:
self._macroscopic = macroscopic
else:
msg = 'Unable to add a Macroscopic to Material ID="{0}". ' \
msg = 'Unable to add a Macroscopic to Material ID="{}". ' \
'Only one Macroscopic allowed per ' \
'Material.'.format(self._id)
raise ValueError(msg)
@ -422,7 +478,7 @@ class Material(object):
"""
if not isinstance(macroscopic, openmc.Macroscopic):
msg = 'Unable to remove a Macroscopic "{0}" in Material ID="{1}" ' \
msg = 'Unable to remove a Macroscopic "{}" in Material ID="{}" ' \
'since it is not a Macroscopic'.format(self._id, macroscopic)
raise ValueError(msg)
@ -450,23 +506,23 @@ class Material(object):
"""
if self._macroscopic is not None:
msg = 'Unable to add an Element to Material ID="{0}" as a ' \
msg = 'Unable to add an Element to Material ID="{}" as a ' \
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if not isinstance(element, string_types + (openmc.Element,)):
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'non-Element value "{1}"'.format(self._id, element)
msg = 'Unable to add an Element to Material ID="{}" with a ' \
'non-Element value "{}"'.format(self._id, element)
raise ValueError(msg)
if not isinstance(percent, Real):
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'non-floating point value "{1}"'.format(self._id, percent)
msg = 'Unable to add an Element to Material ID="{}" with a ' \
'non-floating point value "{}"'.format(self._id, percent)
raise ValueError(msg)
if percent_type not in ['ao', 'wo']:
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'percent type "{1}"'.format(self._id, percent_type)
msg = 'Unable to add an Element to Material ID="{}" with a ' \
'percent type "{}"'.format(self._id, percent_type)
raise ValueError(msg)
# Copy this Element to separate it from same Element in other Materials
@ -477,14 +533,14 @@ class Material(object):
if enrichment is not None:
if not isinstance(enrichment, Real):
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'non-floating point enrichment value "{1}"'\
msg = 'Unable to add an Element to Material ID="{}" with a ' \
'non-floating point enrichment value "{}"'\
.format(self._id, enrichment)
raise ValueError(msg)
elif element.name != 'U':
msg = 'Unable to use enrichment for element {0} which is not ' \
'uranium for Material ID="{1}"'.format(element.name,
msg = 'Unable to use enrichment for element {} which is not ' \
'uranium for Material ID="{}"'.format(element.name,
self._id)
raise ValueError(msg)
@ -493,8 +549,8 @@ class Material(object):
cv.check_greater_than('enrichment', enrichment, 0., equality=True)
if enrichment > 5.0:
msg = 'A uranium enrichment of {0} was given for Material ID='\
'"{1}". OpenMC assumes the U234/U235 mass ratio is '\
msg = 'A uranium enrichment of {} was given for Material ID='\
'"{}". OpenMC assumes the U234/U235 mass ratio is '\
'constant at 0.008, which is only valid at low ' \
'enrichments. Consider setting the isotopic ' \
'composition manually for enrichments over 5%.'.\
@ -514,7 +570,7 @@ class Material(object):
"""
if not isinstance(element, openmc.Element):
msg = 'Unable to remove "{0}" in Material ID="{1}" ' \
msg = 'Unable to remove "{}" in Material ID="{}" ' \
'since it is not an Element'.format(self.id, element)
raise ValueError(msg)
@ -534,13 +590,13 @@ class Material(object):
"""
if self._macroscopic is not None:
msg = 'Unable to add an S(a,b) table to Material ID="{0}" as a ' \
msg = 'Unable to add an S(a,b) table to Material ID="{}" as a ' \
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if not isinstance(name, string_types):
msg = 'Unable to add an S(a,b) table to Material ID="{0}" with a ' \
'non-string table name "{1}"'.format(self._id, name)
msg = 'Unable to add an S(a,b) table to Material ID="{}" with a ' \
'non-string table name "{}"'.format(self._id, name)
raise ValueError(msg)
new_name = openmc.data.get_thermal_name(name)
@ -758,6 +814,9 @@ class Material(object):
if len(self._name) > 0:
element.set("name", str(self._name))
if self._depletable:
element.set("depletable", "true")
# Create temperature XML subelement
if self.temperature is not None:
subelement = ET.SubElement(element, "temperature")

66
openmc/mesh.py
View file

@ -160,23 +160,50 @@ class Mesh(EqualityMixin):
string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._width)
return string
@classmethod
def from_hdf5(cls, group):
"""Create mesh from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
Returns
-------
openmc.Mesh
Mesh instance
"""
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
# Read and assign mesh properties
mesh = cls(mesh_id)
mesh.type = group['type'].value.decode()
mesh.dimension = group['dimension'].value
mesh.lower_left = group['lower_left'].value
mesh.upper_right = group['upper_right'].value
mesh.width = group['width'].value
return mesh
def cell_generator(self):
"""Generator function to traverse through every [i,j,k] index
of the mesh.
"""Generator function to traverse through every [i,j,k] index of the
mesh
For example the following code:
.. code-block:: python
for mesh_index in mymesh.cell_generator():
print mesh_index
print(mesh_index)
will produce the following output for a 3-D 2x2x2 mesh in mymesh::
[1, 1, 1]
[1, 1, 2]
[2, 1, 1]
[1, 2, 1]
[1, 2, 2]
[2, 2, 1]
...
@ -186,13 +213,13 @@ class Mesh(EqualityMixin):
for x in range(self.dimension[0]):
yield [x + 1, 1, 1]
elif len(self.dimension) == 2:
for x in range(self.dimension[0]):
for y in range(self.dimension[1]):
for y in range(self.dimension[1]):
for x in range(self.dimension[0]):
yield [x + 1, y + 1, 1]
else:
for x in range(self.dimension[0]):
for z in range(self.dimension[2]):
for y in range(self.dimension[1]):
for z in range(self.dimension[2]):
for x in range(self.dimension[0]):
yield [x + 1, y + 1, z + 1]
def to_xml_element(self):
@ -294,25 +321,17 @@ class Mesh(EqualityMixin):
# Build the universes which will be used for each of the [i,j,k]
# locations within the mesh.
# We will also have to build cells to assign to these universes
universes = np.ndarray(self.dimension[::-1], dtype=np.object)
# We will concurrently build cells to assign to these universes
cells = []
universes = []
for [i, j, k] in self.cell_generator():
if len(self.dimension) == 1:
universes[i - 1] = openmc.Universe()
cells.append(openmc.Cell())
universes[i - 1].add_cells([cells[-1]])
elif len(self.dimension) == 2:
universes[j - 1, i - 1] = openmc.Universe()
cells.append(openmc.Cell())
universes[j - 1, i - 1].add_cells([cells[-1]])
else:
universes[k - 1, j - 1, i - 1] = openmc.Universe()
cells.append(openmc.Cell())
universes[k - 1, j - 1, i - 1].add_cells([cells[-1]])
cells.append(openmc.Cell())
universes.append(openmc.Universe())
universes[-1].add_cell(cells[-1])
lattice = openmc.RectLattice()
lattice.lower_left = self.lower_left
lattice.universes = np.reshape(universes, self.dimension)
if self.width is not None:
lattice.pitch = self.width
@ -332,7 +351,6 @@ class Mesh(EqualityMixin):
dz = ((self.upper_right[2] - self.lower_left[2]) /
self.dimension[2])
lattice.pitch = [dx, dy, dz]
lattice.universes = universes
# Fill Cell with the Lattice
root_cell.fill = lattice

69
openmc/mgxs/library.py
View file

@ -3,7 +3,7 @@ import os
import copy
import pickle
from numbers import Integral
from collections import OrderedDict
from collections import OrderedDict, Iterable
from warnings import warn
from six import string_types
@ -31,7 +31,7 @@ class Library(object):
Parameters
----------
openmc_geometry : openmc.Geometry
geometry : openmc.Geometry
A geometry which has been initialized with a root universe
by_nuclide : bool
If true, computes cross sections for each nuclide in each domain
@ -43,12 +43,8 @@ class Library(object):
Attributes
----------
openmc_geometry : openmc.Geometry
geometry : openmc.Geometry
An geometry which has been initialized with a root universe
opencg_geometry : opencg.Geometry
An OpenCG geometry object equivalent to the OpenMC geometry
encapsulated by the summary file. Use of this attribute requires
installation of the OpenCG Python module.
by_nuclide : bool
If true, computes cross sections for each nuclide in each domain
mgxs_types : Iterable of str
@ -100,12 +96,11 @@ class Library(object):
"""
def __init__(self, openmc_geometry, by_nuclide=False,
def __init__(self, geometry, by_nuclide=False,
mgxs_types=None, name=''):
self._name = ''
self._openmc_geometry = None
self._opencg_geometry = None
self._geometry = None
self._by_nuclide = None
self._mgxs_types = []
self._domain_type = None
@ -126,7 +121,7 @@ class Library(object):
self._estimator = None
self.name = name
self.openmc_geometry = openmc_geometry
self.geometry = geometry
self.by_nuclide = by_nuclide
if mgxs_types is not None:
@ -139,8 +134,7 @@ class Library(object):
if existing is None:
clone = type(self).__new__(type(self))
clone._name = self.name
clone._openmc_geometry = self.openmc_geometry
clone._opencg_geometry = None
clone._geometry = self.geometry
clone._by_nuclide = self.by_nuclide
clone._mgxs_types = self.mgxs_types
clone._domain_type = self.domain_type
@ -175,15 +169,8 @@ class Library(object):
return existing
@property
def openmc_geometry(self):
return self._openmc_geometry
@property
def opencg_geometry(self):
if self._opencg_geometry is None:
from openmc.opencg_compatible import get_opencg_geometry
self._opencg_geometry = get_opencg_geometry(self._openmc_geometry)
return self._opencg_geometry
def geometry(self):
return self._geometry
@property
def name(self):
@ -205,11 +192,11 @@ class Library(object):
def domains(self):
if self._domains == 'all':
if self.domain_type == 'material':
return self.openmc_geometry.get_all_materials()
return list(self.geometry.get_all_materials().values())
elif self.domain_type in ['cell', 'distribcell']:
return self.openmc_geometry.get_all_material_cells()
return list(self.geometry.get_all_material_cells().values())
elif self.domain_type == 'universe':
return self.openmc_geometry.get_all_universes()
return list(self.geometry.get_all_universes().values())
elif self.domain_type == 'mesh':
raise ValueError('Unable to get domains for Mesh domain type')
else:
@ -277,11 +264,10 @@ class Library(object):
def sparse(self):
return self._sparse
@openmc_geometry.setter
def openmc_geometry(self, openmc_geometry):
cv.check_type('openmc_geometry', openmc_geometry, openmc.Geometry)
self._openmc_geometry = openmc_geometry
self._opencg_geometry = None
@geometry.setter
def geometry(self, geometry):
cv.check_type('geometry', geometry, openmc.Geometry)
self._geometry = geometry
@name.setter
def name(self, name):
@ -329,16 +315,16 @@ class Library(object):
# User specified a list of material, cell or universe domains
else:
if self.domain_type == 'material':
cv.check_iterable_type('domain', domains, openmc.Material)
all_domains = self.openmc_geometry.get_all_materials()
cv.check_type('domain', domains, Iterable, openmc.Material)
all_domains = self.geometry.get_all_materials().values()
elif self.domain_type in ['cell', 'distribcell']:
cv.check_iterable_type('domain', domains, openmc.Cell)
all_domains = self.openmc_geometry.get_all_material_cells()
cv.check_type('domain', domains, Iterable, openmc.Cell)
all_domains = self.geometry.get_all_material_cells().values()
elif self.domain_type == 'universe':
cv.check_iterable_type('domain', domains, openmc.Universe)
all_domains = self.openmc_geometry.get_all_universes()
cv.check_type('domain', domains, Iterable, openmc.Universe)
all_domains = self.geometry.get_all_universes().values()
elif self.domain_type == 'mesh':
cv.check_iterable_type('domain', domains, openmc.Mesh)
cv.check_type('domain', domains, Iterable, openmc.Mesh)
# The mesh and geometry are independent, so set all_domains
# to the input domains
@ -353,7 +339,7 @@ class Library(object):
raise ValueError('Domain "{}" could not be found in the '
'geometry.'.format(domain))
self._domains = domains
self._domains = list(domains)
@energy_groups.setter
def energy_groups(self, energy_groups):
@ -597,7 +583,7 @@ class Library(object):
raise ValueError(msg)
self._sp_filename = statepoint._f.filename
self._openmc_geometry = statepoint.summary.openmc_geometry
self._geometry = statepoint.summary.geometry
self._nuclides = statepoint.summary.nuclides
if statepoint.run_mode == 'k-eigenvalue':
@ -1331,6 +1317,7 @@ class Library(object):
root_cell, cells = \
self.domains[0].build_cells(bc)
root.add_cell(root_cell)
geometry = openmc.Geometry()
geometry.root_universe = root
materials = openmc.Materials()
@ -1352,11 +1339,11 @@ class Library(object):
else:
# Create a copy of the Geometry for these Macroscopics
geometry = copy.deepcopy(self.openmc_geometry)
geometry = copy.deepcopy(self.geometry)
materials = openmc.Materials()
# Get all Cells from the Geometry for differentiation
all_cells = geometry.get_all_material_cells()
all_cells = geometry.get_all_material_cells().values()
# Create the xsdata object and add it to the mgxs_file
for i, domain in enumerate(self.domains):

1167
openmc/mgxs/mgxs.py
View file

File diff suppressed because it is too large Load diff

454
openmc/mgxs_library.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
import copy
from numbers import Real, Integral
import os
@ -14,10 +14,17 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \
# Supported incoming particle MGXS angular treatment representations
_REPRESENTATIONS = ['isotropic', 'angle']
# Supported scattering angular distribution representations
_SCATTER_TYPES = ['tabular', 'legendre', 'histogram']
# List of MGXS indexing schemes
_XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[DG][G]",
"[DG][G']", "[DG][G][G']"]
# Number of mu points for conversion between scattering formats
_NMU = 257
class XSdata(object):
"""A multi-group cross section data set providing all the
@ -185,6 +192,52 @@ class XSdata(object):
self._inverse_velocity = len(temperatures) * [None]
self._xs_shapes = None
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, copy it
if existing is None:
clone = type(self).__new__(type(self))
clone._name = self.name
clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
clone._num_delayed_groups = self.num_delayed_groups
clone._temperatures = copy.deepcopy(self.temperatures, memo)
clone._representation = self.representation
clone._atomic_weight_ratio = self._atomic_weight_ratio
clone._fissionable = self._fissionable
clone._scatter_format = self._scatter_format
clone._order = self._order
clone._num_polar = self._num_polar
clone._num_azimuthal = self._num_azimuthal
clone._total = copy.deepcopy(self._total, memo)
clone._absorption = copy.deepcopy(self._absorption, memo)
clone._scatter_matrix = copy.deepcopy(self._scatter_matrix, memo)
clone._multiplicity_matrix = \
copy.deepcopy(self._multiplicity_matrix, memo)
clone._fission = copy.deepcopy(self._fission, memo)
clone._nu_fission = copy.deepcopy(self._nu_fission, memo)
clone._prompt_nu_fission = \
copy.deepcopy(self._prompt_nu_fission, memo)
clone._delayed_nu_fission = \
copy.deepcopy(self._delayed_nu_fission, memo)
clone._kappa_fission = copy.deepcopy(self._kappa_fission, memo)
clone._chi = copy.deepcopy(self._chi, memo)
clone._chi_prompt = copy.deepcopy(self._chi_prompt, memo)
clone._chi_delayed = copy.deepcopy(self._chi_delayed, memo)
clone._beta = copy.deepcopy(self._beta, memo)
clone._decay_rate = copy.deepcopy(self._decay_rate, memo)
clone._inverse_velocity = \
copy.deepcopy(self._inverse_velocity, memo)
clone._xs_shapes = copy.deepcopy(self._xs_shapes, memo)
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
@property
def name(self):
return self._name
@ -298,8 +351,8 @@ class XSdata(object):
self._xs_shapes["[DG]"] = (self.num_delayed_groups,)
self._xs_shapes["[DG][G]"] = (self.num_delayed_groups,
self.energy_groups.num_groups)
self._xs_shapes["[DG'][G']"] = (self.num_delayed_groups,
self.energy_groups.num_groups)
self._xs_shapes["[DG][G']"] = (self.num_delayed_groups,
self.energy_groups.num_groups)
self._xs_shapes["[DG][G][G']"] = (self.num_delayed_groups,
self.energy_groups.num_groups,
self.energy_groups.num_groups)
@ -318,15 +371,14 @@ class XSdata(object):
@name.setter
def name(self, name):
check_type('name for XSdata', name, string_types)
self._name = name
@energy_groups.setter
def energy_groups(self, energy_groups):
# Check validity of energy_groups
check_type('energy_groups', energy_groups, openmc.mgxs.EnergyGroups)
if energy_groups.group_edges is None:
msg = 'Unable to assign an EnergyGroups object ' \
'with uninitialized group edges'
@ -337,7 +389,6 @@ class XSdata(object):
@num_delayed_groups.setter
def num_delayed_groups(self, num_delayed_groups):
# Check validity of num_delayed_groups
check_type('num_delayed_groups', num_delayed_groups, Integral)
check_less_than('num_delayed_groups', num_delayed_groups,
openmc.mgxs.MAX_DELAYED_GROUPS, equality=True)
@ -348,14 +399,12 @@ class XSdata(object):
@representation.setter
def representation(self, representation):
# Check it is of valid type.
check_value('representation', representation, _REPRESENTATIONS)
self._representation = representation
@atomic_weight_ratio.setter
def atomic_weight_ratio(self, atomic_weight_ratio):
# Check validity of type and that the atomic_weight_ratio value is > 0
check_type('atomic_weight_ratio', atomic_weight_ratio, Real)
check_greater_than('atomic_weight_ratio', atomic_weight_ratio, 0.0)
self._atomic_weight_ratio = atomic_weight_ratio
@ -369,14 +418,12 @@ class XSdata(object):
@scatter_format.setter
def scatter_format(self, scatter_format):
# check to see it is of a valid type and value
check_value('scatter_format', scatter_format, _SCATTER_TYPES)
self._scatter_format = scatter_format
@order.setter
def order(self, order):
# Check type and value
check_type('order', order, Integral)
check_greater_than('order', order, 0, equality=True)
self._order = order
@ -384,7 +431,6 @@ class XSdata(object):
@num_polar.setter
def num_polar(self, num_polar):
# Make sure we have positive ints
check_type('num_polar', num_polar, Integral)
check_greater_than('num_polar', num_polar, 0)
self._num_polar = num_polar
@ -1028,12 +1074,12 @@ class XSdata(object):
def set_nu_fission_mgxs(self, nu_fission, temperature=294.,
nuclide='total', xs_type='macro', subdomain=None):
"""This method allows for an openmc.mgxs.NuFissionXS
"""This method allows for an openmc.mgxs.FissionXS
to be used to set the nu-fission cross section for this XSdata object.
Parameters
----------
nu_fission: openmc.mgxs.NuFissionXS
nu_fission: openmc.mgxs.FissionXS
MGXS Object containing the nu-fission cross section
for the domain of interest.
temperature : float
@ -1056,8 +1102,10 @@ class XSdata(object):
"""
check_type('nu_fission', nu_fission, (openmc.mgxs.NuFissionXS,
check_type('nu_fission', nu_fission, (openmc.mgxs.FissionXS,
openmc.mgxs.NuFissionMatrixXS))
if isinstance(nu_fission, openmc.mgxs.FissionXS):
check_value('nu', nu_fission.nu, [True])
check_value('energy_groups', nu_fission.energy_groups,
[self.energy_groups])
check_value('domain_type', nu_fission.domain_type,
@ -1078,13 +1126,13 @@ class XSdata(object):
subdomain=None):
"""Sets the prompt-nu-fission cross section.
This method allows for an openmc.mgxs.PromptNuFissionXS or
openmc.mgxs.PromptNuFissionMatrixXS to be used to set the
prompt-nu-fission cross section for this XSdata object.
This method allows for an openmc.mgxs.FissionXS or
openmc.mgxs.NuFissionMatrixXS to be used to set the prompt-nu-fission
cross section for this XSdata object.
Parameters
----------
prompt_nu_fission: openmc.mgxs.PromptNuFissionXS or openmc.mgxs.PromptNuFissionMatrixXS
prompt_nu_fission: openmc.mgxs.FissionXS or openmc.mgxs.NuFissionMatrixXS
MGXS Object containing the prompt-nu-fission cross section
for the domain of interest.
temperature : float
@ -1108,8 +1156,8 @@ class XSdata(object):
"""
check_type('prompt_nu_fission', prompt_nu_fission,
(openmc.mgxs.PromptNuFissionXS,
openmc.mgxs.PromptNuFissionMatrixXS))
(openmc.mgxs.FissionXS, openmc.mgxs.NuFissionMatrixXS))
check_value('prompt', prompt_nu_fission.prompt, [True])
check_value('energy_groups', prompt_nu_fission.energy_groups,
[self.energy_groups])
check_value('domain_type', prompt_nu_fission.domain_type,
@ -1261,12 +1309,12 @@ class XSdata(object):
def set_chi_prompt_mgxs(self, chi_prompt, temperature=294.,
nuclide='total', xs_type='macro', subdomain=None):
"""This method allows for an openmc.mgxs.ChiPrompt
to be used to set chi-prompt for this XSdata object.
"""This method allows for an openmc.mgxs.Chi to be used to set
chi-prompt for this XSdata object.
Parameters
----------
chi_prompt: openmc.mgxs.ChiPrompt
chi_prompt: openmc.mgxs.Chi
MGXS Object containing chi-prompt for the domain of interest.
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
@ -1288,7 +1336,8 @@ class XSdata(object):
"""
check_type('chi_prompt', chi_prompt, openmc.mgxs.ChiPrompt)
check_type('chi_prompt', chi_prompt, openmc.mgxs.Chi)
check_value('prompt', chi_prompt.prompt, [True])
check_value('energy_groups', chi_prompt.energy_groups,
[self.energy_groups])
check_value('domain_type', chi_prompt.domain_type,
@ -1552,7 +1601,7 @@ class XSdata(object):
"""
check_type('nuscatter', nuscatter, (openmc.mgxs.NuScatterMatrixXS,
check_type('nuscatter', nuscatter, (openmc.mgxs.ScatterMatrixXS,
openmc.mgxs.MultiplicityMatrixXS))
check_value('energy_groups', nuscatter.energy_groups,
[self.energy_groups])
@ -1622,7 +1671,8 @@ class XSdata(object):
"""
check_type('inverse_velocity', inverse_velocity, openmc.mgxs.InverseVelocity)
check_type('inverse_velocity', inverse_velocity,
openmc.mgxs.InverseVelocity)
check_value('energy_groups', inverse_velocity.energy_groups,
[self.energy_groups])
check_value('domain_type', inverse_velocity.domain_type,
@ -1634,6 +1684,269 @@ class XSdata(object):
self._inverse_velocity[i] = inverse_velocity.get_xs(
nuclides=nuclide, xs_type=xs_type, subdomains=subdomain)
def convert_representation(self, target_representation, num_polar=None,
num_azimuthal=None):
"""Produce a new XSdata object with the same data, but converted to the
new representation (isotropic or angle-dependent).
This method cannot be used to change the number of polar or
azimuthal bins of an XSdata object that already uses an angular
representation. Finally, this method simply uses an arithmetic mean to
convert from an angular to isotropic representation; no flux-weighting
is applied and therefore reaction rates will not be preserved.
Parameters
----------
target_representation : {'isotropic', 'angle'}
Representation of the MGXS (isotropic or angle-dependent flux
weighting).
num_polar : int, optional
Number of equal width angular bins that the polar angular
domain is subdivided into. This is required when
:param:`target_representation` is "angle".
num_azimuthal : int, optional
Number of equal width angular bins that the azimuthal angular
domain is subdivided into. This is required when
:param:`target_representation` is "angle".
Returns
-------
openmc.XSdata
Multi-group cross section data with the same data as self, but
represented as specified in :param:`target_representation`.
"""
check_value('target_representation', target_representation,
_REPRESENTATIONS)
if target_representation == 'angle':
check_type('num_polar', num_polar, Integral)
check_type('num_azimuthal', num_azimuthal, Integral)
check_greater_than('num_polar', num_polar, 0)
check_greater_than('num_azimuthal', num_azimuthal, 0)
xsdata = copy.deepcopy(self)
# First handle the case where the current and requested
# representations are the same
if target_representation == self.representation:
# Check to make sure the num_polar and num_azimuthal values match
if target_representation == 'angle':
if num_polar != self.num_polar or num_azimuthal != self.num_azimuthal:
raise ValueError("Cannot translate between `angle`"
" representations with different angle"
" bin structures")
# Nothing to do as the same structure was requested
return xsdata
xsdata.representation = target_representation
# We have different actions depending on the representation conversion
if target_representation == 'isotropic':
# This is not needed for the correct functionality, but these
# values are changed back to None for clarity
xsdata._num_polar = None
xsdata._num_azimuthal = None
elif target_representation == 'angle':
xsdata.num_polar = num_polar
xsdata.num_azimuthal = num_azimuthal
# Reset xs_shapes so it is recalculated the next time it is needed
xsdata._xs_shapes = None
for i, temp in enumerate(xsdata.temperatures):
for xs in ['total', 'absorption', 'fission', 'nu_fission',
'scatter_matrix', 'multiplicity_matrix',
'prompt_nu_fission', 'delayed_nu_fission',
'kappa_fission', 'chi', 'chi_prompt', 'chi_delayed',
'beta', 'decay_rate', 'inverse_velocity']:
# Get the original data
orig_data = getattr(self, '_' + xs)[i]
if orig_data is not None:
if target_representation == 'isotropic':
# Since we are going from angle to isotropic, the
# current data is just the average over the angle bins
new_data = orig_data.mean(axis=(0, 1))
elif target_representation == 'angle':
# Since we are going from isotropic to angle, the
# current data is just copied for every angle bin
new_shape = (num_polar, num_azimuthal) + \
orig_data.shape
new_data = np.resize(orig_data, new_shape)
setter = getattr(xsdata, 'set_' + xs)
setter(new_data, temp)
return xsdata
def convert_scatter_format(self, target_format, target_order=None):
"""Produce a new MGXSLibrary object with the same data, but converted
to the new scatter format and order
Parameters
----------
target_format : {'tabular', 'legendre', 'histogram'}
Representation of the scattering angle distribution
target_order : int
Either the Legendre target_order, number of bins, or number of
points used to describe the angular distribution associated with
each group-to-group transfer probability
Returns
-------
openmc.XSdata
Multi-group cross section data with the same data as in self, but
represented as specified in :param:`target_format`.
"""
from scipy.interpolate import interp1d
from scipy.integrate import simps
from scipy.special import eval_legendre
check_value('target_format', target_format, _SCATTER_TYPES)
check_type('target_order', target_order, Integral)
if target_format == 'legendre':
check_greater_than('target_order', target_order, 0, equality=True)
else:
check_greater_than('target_order', target_order, 0)
xsdata = copy.deepcopy(self)
xsdata.scatter_format = target_format
xsdata.order = target_order
# Reset and re-generate XSdata.xs_shapes with the new scattering format
xsdata._xs_shapes = None
for i, temp in enumerate(xsdata.temperatures):
orig_data = self._scatter_matrix[i]
new_shape = orig_data.shape[:-1] + (xsdata.num_orders,)
new_data = np.zeros(new_shape)
if self.scatter_format == 'legendre':
if target_format == 'legendre':
# Then we are changing orders and only need to change
# dimensionality of the mu data and pad/truncate as needed
order = min(xsdata.num_orders, self.num_orders)
new_data[..., :order] = orig_data[..., :order]
elif target_format == 'tabular':
mu = np.linspace(-1, 1, xsdata.num_orders)
# Evaluate the legendre on the mu grid
for imu in range(len(mu)):
new_data[..., imu] = \
np.sum((l + 0.5) * eval_legendre(l, mu[imu]) *
orig_data[..., l]
for l in range(self.num_orders))
elif target_format == 'histogram':
# This code uses the vectorized integration capabilities
# instead of having an isotropic and angle representation
# path.
# Set the histogram mu grid
mu = np.linspace(-1, 1, xsdata.num_orders + 1)
# For every bin perform simpson integration of a finely
# sampled orig_data
for h_bin in range(xsdata.num_orders):
mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU)
table_fine = np.zeros(new_data.shape[:-1] + (_NMU,))
for imu in range(len(mu_fine)):
table_fine[..., imu] = \
np.sum((l + 0.5) *
eval_legendre(l, mu_fine[imu]) *
orig_data[..., l]
for l in range(self.num_orders))
new_data[..., h_bin] = simps(table_fine, mu_fine)
elif self.scatter_format == 'tabular':
# Calculate the mu points of the current data
mu_self = np.linspace(-1, 1, self.num_orders)
if target_format == 'legendre':
# Find the Legendre coefficients via integration. To best
# use the vectorized integration capabilities of scipy,
# this is done with fixed sample integration routines.
mu_fine = np.linspace(-1, 1, _NMU)
y = [interp1d(mu_self, orig_data)(mu_fine) *
eval_legendre(l, mu_fine)
for l in range(xsdata.num_orders)]
for l in range(xsdata.num_orders):
new_data[..., l] = simps(y[l], mu_fine)
elif target_format == 'tabular':
# Simply use an interpolating function to get the new data
mu = np.linspace(-1, 1, xsdata.num_orders)
new_data[..., :] = interp1d(mu_self, orig_data)(mu)
elif target_format == 'histogram':
# Use an interpolating function to do the bin-wise
# integrals
mu = np.linspace(-1, 1, xsdata.num_orders + 1)
# Like the tabular -> legendre path above, this code will
# be written to utilize the vectorized integration
# capabilities instead of having an isotropic and
# angle representation path.
interp = interp1d(mu_self, orig_data)
for h_bin in range(xsdata.num_orders):
mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU)
new_data[..., h_bin] = simps(interp(mu_fine), mu_fine)
elif self.scatter_format == 'histogram':
# The histogram format does not have enough information to
# convert to the other forms without inducing some amount of
# error. We will make the assumption that the center of the bin
# has the value of the bin. The mu=-1 and 1 points will be
# extrapolated from the shape.
mu_midpoint = np.linspace(-1, 1, self.num_orders,
endpoint=False)
mu_midpoint += (mu_midpoint[1] - mu_midpoint[0]) * 0.5
interp = interp1d(mu_midpoint, orig_data,
fill_value='extrapolate')
# Now get the distribution normalization factor to take from
# an integral quantity to a point-wise quantity
norm = float(self.num_orders) / 2.0
# We now have a tabular distribution in tab_data on mu_self.
# We now proceed just like the tabular branch above.
if target_format == 'legendre':
# find the legendre coefficients via integration. To best
# use the vectorized integration capabilities of scipy,
# this will be done with fixed sample integration routines.
mu_fine = np.linspace(-1, 1, _NMU)
y = [interp(mu_fine) * norm * eval_legendre(l, mu_fine)
for l in range(xsdata.num_orders)]
for l in range(xsdata.num_orders):
new_data[..., l] = simps(y[l], mu_fine)
elif target_format == 'tabular':
# Simply use an interpolating function to get the new data
mu = np.linspace(-1, 1, xsdata.num_orders)
new_data[..., :] = interp(mu) * norm
elif target_format == 'histogram':
# Use an interpolating function to do the bin-wise
# integrals
mu = np.linspace(-1, 1, xsdata.num_orders + 1)
# Like the tabular -> legendre path above, this code will
# be written to utilize the vectorized integration
# capabilities instead of having an isotropic and
# angle representation path.
for h_bin in range(xsdata.num_orders):
mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU)
new_data[..., h_bin] = \
norm * simps(interp(mu_fine), mu_fine)
# Remove small values resulting from numerical precision issues
new_data[..., np.abs(new_data) < 1.E-10] = 0.
xsdata.set_scatter_matrix(new_data, temp)
return xsdata
def to_hdf5(self, file):
"""Write XSdata to an HDF5 file
@ -1757,7 +2070,7 @@ class XSdata(object):
elif self.representation == 'angle':
matrix = \
self._scatter_matrix[i][p, a, g_in, :, 0]
elif self.scatter_format == 'histogram':
else:
if self.representation == 'isotropic':
matrix = \
np.sum(self._scatter_matrix[i][g_in, :, :],
@ -1995,6 +2308,24 @@ class MGXSLibrary(object):
self.num_delayed_groups = num_delayed_groups
self._xsdatas = []
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, copy it
if existing is None:
clone = type(self).__new__(type(self))
clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
clone._num_delayed_groups = self.num_delayed_groups
clone._xsdatas = copy.deepcopy(self.xsdatas, memo)
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
@property
def energy_groups(self):
return self._energy_groups
@ -2099,6 +2430,75 @@ class MGXSLibrary(object):
result = xsdata
return result
def convert_representation(self, target_representation, num_polar=None,
num_azimuthal=None):
"""Produce a new XSdata object with the same data, but converted to the
new representation (isotropic or angle-dependent).
This method cannot be used to change the number of polar or
azimuthal bins of an XSdata object that already uses an angular
representation. Finally, this method simply uses an arithmetic mean to
convert from an angular to isotropic representation; no flux-weighting
is applied and therefore the reaction rates will not be preserved.
Parameters
----------
target_representation : {'isotropic', 'angle'}
Representation of the MGXS (isotropic or angle-dependent flux
weighting).
num_polar : int, optional
Number of equal width angular bins that the polar angular
domain is subdivided into. This is required when
:param:`target_representation` is "angle".
num_azimuthal : int, optional
Number of equal width angular bins that the azimuthal angular
domain is subdivided into. This is required when
:param:`target_representation` is "angle".
Returns
-------
openmc.MGXSLibrary
Multi-group Library with the same data as self, but represented as
specified in :param:`target_representation`.
"""
library = copy.deepcopy(self)
for i, xsdata in enumerate(self.xsdatas):
library.xsdatas[i] = \
xsdata.convert_representation(target_representation,
num_polar, num_azimuthal)
return library
def convert_scatter_format(self, target_format, target_order):
"""Produce a new MGXSLibrary object with the same data, but converted
to the new scatter format and order
Parameters
----------
target_format : {'tabular', 'legendre', 'histogram'}
Representation of the scattering angle distribution
target_order : int
Either the Legendre target_order, number of bins, or number of
points used to describe the angular distribution associated with
each group-to-group transfer probability
Returns
-------
openmc.MGXSLibrary
Multi-group Library with the same data as self, but with the
scatter format represented as specified in :param:`target_format`
and :param:`target_order`.
"""
library = copy.deepcopy(self)
for i, xsdata in enumerate(self.xsdatas):
library.xsdatas[i] = \
xsdata.convert_scatter_format(target_format, target_order)
return library
def export_to_hdf5(self, filename='mgxs.h5'):
"""Create an hdf5 file that can be used for a simulation.

28
openmc/particle_restart.py
View file

@ -1,5 +1,9 @@
import h5py
import openmc.checkvalue as cv
_VERSION_PARTICLE_RESTART = 2
class Particle(object):
"""Information used to restart a specific particle that caused a simulation to
fail.
@ -13,9 +17,9 @@ class Particle(object):
----------
current_batch : int
The batch containing the particle
gen_per_batch : int
generations_per_batch : int
Number of generations per batch
current_gen : int
current_generation : int
The generation containing the particle
n_particles : int
Number of particles per generation
@ -37,31 +41,25 @@ class Particle(object):
def __init__(self, filename):
self._f = h5py.File(filename, 'r')
# Ensure filetype and revision are correct
if 'filetype' not in self._f or self._f[
'filetype'].value.decode() != 'particle restart':
raise IOError('{} is not a particle restart file.'.format(filename))
if self._f['revision'].value != 1:
raise IOError('Particle restart file has a file revision of {} '
'which is not consistent with the revision this '
'version of OpenMC expects ({}).'.format(
self._f['revision'].value, 1))
# Ensure filetype and version are correct
cv.check_filetype_version(self._f, 'particle restart',
_VERSION_PARTICLE_RESTART)
@property
def current_batch(self):
return self._f['current_batch'].value
@property
def current_gen(self):
return self._f['current_gen'].value
def current_generation(self):
return self._f['current_generation'].value
@property
def energy(self):
return self._f['energy'].value
@property
def gen_per_batch(self):
return self._f['gen_per_batch'].value
def generations_per_batch(self):
return self._f['generations_per_batch'].value
@property
def id(self):

4
openmc/plots.py
View file

@ -351,11 +351,11 @@ class Plot(object):
# Generate random colors for each feature
self.col_spec = {}
for domain in domains:
for domain_id in domains:
r = np.random.randint(0, 256)
g = np.random.randint(0, 256)
b = np.random.randint(0, 256)
self.col_spec[domain] = (r, g, b)
self.col_spec[domain_id] = (r, g, b)
def highlight_domains(self, geometry, domains, seed=1,
alpha=0.5, background='gray'):

24
openmc/plotter.py
View file

@ -676,7 +676,7 @@ def calculate_mgxs(this, types, orders=None, temperature=294.,
for line in range(len(types)):
for g in range(library.energy_groups.num_groups):
data[g * 2: g * 2 + 2] = mgxs[line, g]
data[line, g * 2: g * 2 + 2] = mgxs[line, g]
return energy_grid[::-1], data
@ -775,28 +775,28 @@ def _calculate_mgxs_nuc_macro(this, types, library, orders=None,
data[i, :] = temp_data[orders[i]]
else:
data[i, :] = np.sum(temp_data[:])
elif shape in (xsdata.xs_shapes["[G'][DG]"],
xsdata.xs_shapes["[G][DG]"]):
elif shape in (xsdata.xs_shapes["[DG][G']"],
xsdata.xs_shapes["[DG][G]"]):
# Then we have an array vs groups with values for each
# delayed group. The user-provided value of orders tells us
# which delayed group we want. If none are provided, then
# we sum all the delayed groups together.
if orders[i]:
if orders[i] < len(shape[1]):
data[i, :] = temp_data[:, orders[i]]
if orders[i] < len(shape[0]):
data[i, :] = temp_data[orders[i], :]
else:
data[i, :] = np.sum(temp_data[:, :], axis=1)
elif shape == xsdata.xs_shapes["[G][G'][DG]"]:
data[i, :] = np.sum(temp_data[:, :], axis=0)
elif shape == xsdata.xs_shapes["[DG][G][G']"]:
# Then we have a delayed group matrix. We will first
# remove the outgoing group dependency
temp_data = np.sum(temp_data, axis=1)
temp_data = np.sum(temp_data, axis=-1)
# And then proceed in exactly the same manner as the
# "[G'][DG]" of "[G][DG]" shapes in the previous block.
# "[DG][G']" or "[DG][G]" shapes in the previous block.
if orders[i]:
if orders[i] < len(shape[1]):
data[i, :] = temp_data[:, orders[i]]
if orders[i] < len(shape[0]):
data[i, :] = temp_data[orders[i], :]
else:
data[i, :] = np.sum(temp_data[:, :], axis=1)
data[i, :] = np.sum(temp_data[:, :], axis=0)
elif shape == xsdata.xs_shapes["[G][G'][Order]"]:
# This is a scattering matrix with angular data
# First remove the outgoing group dependence

32
openmc/settings.py
View file

@ -11,6 +11,9 @@ from openmc.clean_xml import clean_xml_indentation
import openmc.checkvalue as cv
from openmc import Nuclide, VolumeCalculation, Source, Mesh
_RUN_MODES = ['eigenvalue', 'fixed source', 'plot', 'volume',
'particle restart']
class Settings(object):
"""Settings used for an OpenMC simulation.
@ -81,7 +84,7 @@ class Settings(object):
The elastic scattering model to use for resonant isotopes
run_cmfd : bool
Indicate if coarse mesh finite difference acceleration is to be used
run_mode : {'eigenvalue' or 'fixed source'}
run_mode : {'eigenvalue', 'fixed source', 'plot', 'volume', 'particle restart'}
The type of calculation to perform (default is 'eigenvalue')
seed : int
Seed for the linear congruential pseudorandom number generator
@ -388,10 +391,7 @@ class Settings(object):
@run_mode.setter
def run_mode(self, run_mode):
if run_mode not in ['eigenvalue', 'fixed source']:
msg = 'Unable to set run mode to "{0}". Only "eigenvalue" ' \
'and "fixed source" are supported."'.format(run_mode)
raise ValueError(msg)
cv.check_value('run mode', run_mode, _RUN_MODES)
self._run_mode = run_mode
@batches.setter
@ -794,18 +794,8 @@ class Settings(object):
self._create_fission_neutrons = create_fission_neutrons
def _create_run_mode_subelement(self, root):
if self.run_mode == 'eigenvalue':
elem = ET.SubElement(root, "eigenvalue")
self._create_particles_subelement(elem)
self._create_batches_subelement(elem)
self._create_inactive_subelement(elem)
self._create_generations_per_batch_subelement(elem)
self._create_keff_trigger_subelement(elem)
else:
elem = ET.SubElement(root, "fixed_source")
self._create_particles_subelement(elem)
self._create_batches_subelement(elem)
elem = ET.SubElement(root, "run_mode")
elem.text = self._run_mode
def _create_batches_subelement(self, run_mode_element):
if self._batches is not None:
@ -814,8 +804,7 @@ class Settings(object):
def _create_generations_per_batch_subelement(self, run_mode_element):
if self._generations_per_batch is not None:
element = ET.SubElement(run_mode_element,
"generations_per_batch")
element = ET.SubElement(run_mode_element, "generations_per_batch")
element.text = str(self._generations_per_batch)
def _create_inactive_subelement(self, run_mode_element):
@ -1081,6 +1070,11 @@ class Settings(object):
root_element = ET.Element("settings")
self._create_run_mode_subelement(root_element)
self._create_particles_subelement(root_element)
self._create_batches_subelement(root_element)
self._create_inactive_subelement(root_element)
self._create_generations_per_batch_subelement(root_element)
self._create_keff_trigger_subelement(root_element)
self._create_source_subelement(root_element)
self._create_output_subelement(root_element)
self._create_statepoint_subelement(root_element)

222
openmc/statepoint.py
View file

@ -10,6 +10,8 @@ import h5py
import openmc
import openmc.checkvalue as cv
_VERSION_STATEPOINT = 16
class StatePoint(object):
"""State information on a simulation at a certain point in time (at the end
@ -50,7 +52,7 @@ class StatePoint(object):
Date and time when simulation began
entropy : numpy.ndarray
Shannon entropy of fission source at each batch
gen_per_batch : Integral
generations_per_batch : int
Number of fission generations per batch
global_tallies : numpy.ndarray of compound datatype
Global tallies for k-effective estimates and leakage. The compound
@ -78,7 +80,7 @@ class StatePoint(object):
path : str
Working directory for simulation
run_mode : str
Simulation run mode, e.g. 'k-eigenvalue'
Simulation run mode, e.g. 'eigenvalue'
runtime : dict
Dictionary whose keys are strings describing various runtime metrics
and whose values are time values in seconds.
@ -109,22 +111,12 @@ class StatePoint(object):
def __init__(self, filename, autolink=True):
self._f = h5py.File(filename, 'r')
self._meshes = {}
self._tallies = {}
self._derivs = {}
# Ensure filetype and revision are correct
try:
if 'filetype' not in self._f or self._f[
'filetype'].value.decode() != 'statepoint':
raise IOError('{} is not a statepoint file.'.format(filename))
except AttributeError:
raise IOError('Could not read statepoint file. This most likely '
'means the statepoint file was produced by a '
'different version of OpenMC than the one you are '
'using.')
if self._f['revision'].value != 15:
raise IOError('Statepoint file has a file revision of {} '
'which is not consistent with the revision this '
'version of OpenMC expects ({}).'.format(
self._f['revision'].value, 15))
# Check filetype and version
cv.check_filetype_version(self._f, 'statepoint', _VERSION_STATEPOINT)
# Set flags for what data has been read
self._meshes_read = False
@ -147,12 +139,9 @@ class StatePoint(object):
vol = openmc.VolumeCalculation.from_hdf5(path_i)
self.add_volume_information(vol)
def close(self):
self._f.close()
@property
def cmfd_on(self):
return self._f['cmfd_on'].value > 0
return self._f.attrs['cmfd_on'] > 0
@property
def cmfd_balance(self):
@ -188,19 +177,19 @@ class StatePoint(object):
@property
def date_and_time(self):
return self._f['date_and_time'].value.decode()
return self._f.attrs['date_and_time'].decode()
@property
def entropy(self):
if self.run_mode == 'k-eigenvalue':
if self.run_mode == 'eigenvalue':
return self._f['entropy'].value
else:
return None
@property
def gen_per_batch(self):
if self.run_mode == 'k-eigenvalue':
return self._f['gen_per_batch'].value
def generations_per_batch(self):
if self.run_mode == 'eigenvalue':
return self._f['generations_per_batch'].value
else:
return None
@ -234,35 +223,35 @@ class StatePoint(object):
@property
def k_generation(self):
if self.run_mode == 'k-eigenvalue':
if self.run_mode == 'eigenvalue':
return self._f['k_generation'].value
else:
return None
@property
def k_combined(self):
if self.run_mode == 'k-eigenvalue':
if self.run_mode == 'eigenvalue':
return self._f['k_combined'].value
else:
return None
@property
def k_col_abs(self):
if self.run_mode == 'k-eigenvalue':
if self.run_mode == 'eigenvalue':
return self._f['k_col_abs'].value
else:
return None
@property
def k_col_tra(self):
if self.run_mode == 'k-eigenvalue':
if self.run_mode == 'eigenvalue':
return self._f['k_col_tra'].value
else:
return None
@property
def k_abs_tra(self):
if self.run_mode == 'k-eigenvalue':
if self.run_mode == 'eigenvalue':
return self._f['k_abs_tra'].value
else:
return None
@ -270,46 +259,12 @@ class StatePoint(object):
@property
def meshes(self):
if not self._meshes_read:
# Initialize dictionaries for the Meshes
# Keys - Mesh IDs
# Values - Mesh objects
self._meshes = {}
# Read the number of Meshes
n_meshes = self._f['tallies/meshes/n_meshes'].value
# Read a list of the IDs for each Mesh
if n_meshes > 0:
# User-defined Mesh IDs
mesh_keys = self._f['tallies/meshes/keys'].value
else:
mesh_keys = []
# Build dictionary of Meshes
base = 'tallies/meshes/mesh '
mesh_group = self._f['tallies/meshes']
# Iterate over all Meshes
for mesh_key in mesh_keys:
# Read the mesh type
mesh_type = self._f['{0}{1}/type'.format(base, mesh_key)].value.decode()
# Read the mesh dimensions, lower-left coordinates,
# upper-right coordinates, and width of each mesh cell
dimension = self._f['{0}{1}/dimension'.format(base, mesh_key)].value
lower_left = self._f['{0}{1}/lower_left'.format(base, mesh_key)].value
upper_right = self._f['{0}{1}/upper_right'.format(base, mesh_key)].value
width = self._f['{0}{1}/width'.format(base, mesh_key)].value
# Create the Mesh and assign properties to it
mesh = openmc.Mesh(mesh_key)
mesh.dimension = dimension
mesh.width = width
mesh.lower_left = lower_left
mesh.upper_right = upper_right
mesh.type = mesh_type
# Add mesh to the global dictionary of all Meshes
self._meshes[mesh_key] = mesh
for group in mesh_group.values():
mesh = openmc.Mesh.from_hdf5(group)
self._meshes[mesh.id] = mesh
self._meshes_read = True
@ -321,7 +276,7 @@ class StatePoint(object):
@property
def n_inactive(self):
if self.run_mode == 'k-eigenvalue':
if self.run_mode == 'eigenvalue':
return self._f['n_inactive'].value
else:
return None
@ -336,7 +291,7 @@ class StatePoint(object):
@property
def path(self):
return self._f['path'].value.decode()
return self._f.attrs['path'].decode()
@property
def run_mode(self):
@ -357,7 +312,7 @@ class StatePoint(object):
@property
def source_present(self):
return self._f['source_present'].value > 0
return self._f.attrs['source_present'] > 0
@property
def sparse(self):
@ -365,68 +320,55 @@ class StatePoint(object):
@property
def tallies(self):
if not self._tallies_read:
# Initialize dictionary for tallies
self._tallies = {}
if self.tallies_present and not self._tallies_read:
# Read the number of tallies
n_tallies = self._f['tallies/n_tallies'].value
tallies_group = self._f['tallies']
n_tallies = tallies_group.attrs['n_tallies']
# Read a list of the IDs for each Tally
if n_tallies > 0:
# OpenMC Tally IDs (redefined internally from user definitions)
tally_keys = self._f['tallies/keys'].value
# Tally user-defined IDs
tally_ids = tallies_group.attrs['ids']
else:
tally_keys = []
tally_ids = []
base = 'tallies/tally '
# Iterate over all tallies
for tally_id in tally_ids:
group = tallies_group['tally {}'.format(tally_id)]
# Iterate over all Tallies
for tally_key in tally_keys:
# Read the Tally size specifications
n_realizations = \
self._f['{0}{1}/n_realizations'.format(base, tally_key)].value
# Read the number of realizations
n_realizations = group['n_realizations'].value
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_id=tally_key)
tally = openmc.Tally(tally_id)
tally._sp_filename = self._f.filename
tally.estimator = self._f['{0}{1}/estimator'.format(
base, tally_key)].value.decode()
tally.name = group['name'].value.decode()
tally.estimator = group['estimator'].value.decode()
tally.num_realizations = n_realizations
# Read derivative information.
if 'derivative' in self._f['{0}{1}'.format(base, tally_key)]:
deriv_id = self._f['{0}{1}/derivative'.format(
base, tally_key)].value
if 'derivative' in group:
deriv_id = group['derivative'].value
tally.derivative = self.tally_derivatives[deriv_id]
# Read the number of Filters
n_filters = \
self._f['{0}{1}/n_filters'.format(base, tally_key)].value
subbase = '{0}{1}/filter '.format(base, tally_key)
# Read all filters
for j in range(1, n_filters+1):
subsubbase = '{0}{1}'.format(subbase, j)
new_filter = openmc.Filter.from_hdf5(self._f[subsubbase],
n_filters = group['n_filters'].value
for j in range(1, n_filters + 1):
filter_group = group['filter {}'.format(j)]
new_filter = openmc.Filter.from_hdf5(filter_group,
meshes=self.meshes)
tally.filters.append(new_filter)
# Read Nuclide bins
nuclide_names = \
self._f['{0}{1}/nuclides'.format(base, tally_key)].value
# Read nuclide bins
nuclide_names = group['nuclides'].value
# Add all Nuclides to the Tally
# Add all nuclides to the Tally
for name in nuclide_names:
nuclide = openmc.Nuclide(name.decode().strip())
tally.nuclides.append(nuclide)
scores = self._f['{0}{1}/score_bins'.format(
base, tally_key)].value
n_score_bins = self._f['{0}{1}/n_score_bins'
.format(base, tally_key)].value
scores = group['score_bins'].value
n_score_bins = group['n_score_bins'].value
# Compute and set the filter strides
for i in range(n_filters):
@ -437,8 +379,7 @@ class StatePoint(object):
tally_filter.stride *= tally.filters[j].num_bins
# Read scattering moment order strings (e.g., P3, Y1,2, etc.)
moments = self._f['{0}{1}/moment_orders'.format(
base, tally_key)].value
moments = group['moment_orders'].value
# Add the scores to the Tally
for j, score in enumerate(scores):
@ -452,7 +393,7 @@ class StatePoint(object):
# Add Tally to the global dictionary of all Tallies
tally.sparse = self.sparse
self._tallies[tally_key] = tally
self._tallies[tally_id] = tally
self._tallies_read = True
@ -460,16 +401,11 @@ class StatePoint(object):
@property
def tallies_present(self):
return self._f['tallies/tallies_present'].value
return self._f.attrs['tallies_present'] > 0
@property
def tally_derivatives(self):
if not self._derivs_read:
# Initialize dictionaries for the Meshes
# Keys - Derivative IDs
# Values - TallyDerivative objects
self._derivs = {}
# Populate the dictionary if any derivatives are present.
if 'derivatives' in self._f['tallies']:
# Read the derivative ids.
@ -478,21 +414,17 @@ class StatePoint(object):
# Create each derivative object and add it to the dictionary.
for d_id in deriv_ids:
base = 'tallies/derivatives/derivative {:d}'.format(d_id)
group = self._f['tallies/derivatives/derivative {}'
.format(d_id)]
deriv = openmc.TallyDerivative(derivative_id=d_id)
deriv.variable = \
self._f[base + '/independent variable'].value.decode()
deriv.variable = group['independent variable'].value.decode()
if deriv.variable == 'density':
deriv.material = self._f[base + '/material'].value
deriv.material = group['material'].value
elif deriv.variable == 'nuclide_density':
deriv.material = self._f[base + '/material'].value
deriv.nuclide = \
self._f[base + '/nuclide'].value.decode()
deriv.material = group['material'].value
deriv.nuclide = group['nuclide'].value.decode()
elif deriv.variable == 'temperature':
deriv.material = self._f[base + '/material'].value
else:
raise RuntimeError('Unrecognized tally differential '
'variable')
deriv.material = group['material'].value
self._derivs[d_id] = deriv
self._derivs_read = True
@ -501,9 +433,7 @@ class StatePoint(object):
@property
def version(self):
return (self._f['version_major'].value,
self._f['version_minor'].value,
self._f['version_release'].value)
return tuple(self._f.attrs['version'])
@property
def summary(self):
@ -705,38 +635,20 @@ class StatePoint(object):
RuntimeWarning)
return
if not isinstance(summary, openmc.summary.Summary):
if not isinstance(summary, openmc.Summary):
msg = 'Unable to link statepoint with "{0}" which ' \
'is not a Summary object'.format(summary)
raise ValueError(msg)
cells = summary.geometry.get_all_cells()
for tally_id, tally in self.tallies.items():
tally.name = summary.tally_names[tally_id]
tally.with_summary = True
for tally_filter in tally.filters:
if isinstance(tally_filter, (openmc.CellFilter,
openmc.DistribcellFilter)):
distribcell_ids = []
for bin in tally_filter.bins:
distribcell_ids.append(summary.cells[bin].id)
tally_filter.bins = distribcell_ids
if isinstance(tally_filter, (openmc.DistribcellFilter)):
cell_id = tally_filter.bins[0]
cell = summary.get_cell_by_id(cell_id)
cell = cells[cell_id]
tally_filter.distribcell_paths = cell.distribcell_paths
if isinstance(tally_filter, openmc.UniverseFilter):
universe_ids = []
for bin in tally_filter.bins:
universe_ids.append(summary.universes[bin].id)
tally_filter.bins = universe_ids
if isinstance(tally_filter, openmc.MaterialFilter):
material_ids = []
for bin in tally_filter.bins:
material_ids.append(summary.materials[bin].id)
tally_filter.bins = material_ids
self._summary = summary

644
openmc/summary.py
View file

@ -5,21 +5,28 @@ import numpy as np
import h5py
import openmc
import openmc.checkvalue as cv
from openmc.region import Region
_VERSION_SUMMARY = 5
class Summary(object):
"""Information summarizing the geometry, materials, and tallies used in a
simulation.
"""Summary of geometry, materials, and tallies used in a simulation.
Attributes
----------
openmc_geometry : openmc.Geometry
An OpenMC geometry object reconstructed from the summary file
opencg_geometry : opencg.Geometry
An OpenCG geometry object equivalent to the OpenMC geometry
encapsulated by the summary file. Use of this attribute requires
installation of the OpenCG Python module.
date_and_time : str
Date and time when simulation began
geometry : openmc.Geometry
The geometry reconstructed from the summary file
materials : openmc.Materials
The materials reconstructed from the summary file
nuclides : dict
Dictionary whose keys are nuclide names and values are atomic weight
ratios.
version: tuple of int
Version of OpenMC
"""
@ -31,229 +38,85 @@ class Summary(object):
raise ValueError(msg)
self._f = h5py.File(filename, 'r')
self._openmc_geometry = None
self._opencg_geometry = None
cv.check_filetype_version(self._f, 'summary', _VERSION_SUMMARY)
self._geometry = openmc.Geometry()
self._materials = openmc.Materials()
self._nuclides = {}
self._read_metadata()
self._read_nuclides()
self._read_geometry()
self._read_tallies()
self._f.close()
@property
def openmc_geometry(self):
return self._openmc_geometry
def date_and_time(self):
return self._f.attrs['date_and_time']
@property
def opencg_geometry(self):
if self._opencg_geometry is None:
from openmc.opencg_compatible import get_opencg_geometry
self._opencg_geometry = get_opencg_geometry(self.openmc_geometry)
return self._opencg_geometry
def geometry(self):
return self._geometry
def _read_metadata(self):
# Read OpenMC version
self.version = [self._f['version_major'].value,
self._f['version_minor'].value,
self._f['version_release'].value]
# Read date and time
self.date_and_time = self._f['date_and_time'][...]
@property
def materials(self):
return self._materials
# Read if continuous-energy or multi-group
self.run_CE = (self._f['run_CE'].value == 1)
@property
def nuclides(self):
return self._nuclides
self.n_batches = self._f['n_batches'].value
self.n_particles = self._f['n_particles'].value
if 'n_inactive' in self._f:
self.n_active = self._f['n_active'].value
self.n_inactive = self._f['n_inactive'].value
self.gen_per_batch = self._f['gen_per_batch'].value
self.n_procs = self._f['n_procs'].value
@property
def version(self):
return tuple(self._f.attrs['openmc_version'])
def _read_nuclides(self):
self.nuclides = {}
n_nuclides = self._f['nuclides/n_nuclides_total'].value
names = self._f['nuclides/names'].value
awrs = self._f['nuclides/awrs'].value
for n in range(n_nuclides):
name = names[n].decode()
name = name[:name.find('.')]
self.nuclides[name] = awrs[n]
for name, awr in zip(names, awrs):
self._nuclides[name.decode()] = awr
def _read_geometry(self):
# Read in and initialize the Materials and Geometry
self._read_materials()
self._read_surfaces()
self._read_cells()
self._read_universes()
self._read_lattices()
self._finalize_geometry()
surfaces = self._read_surfaces()
cells, cell_fills = self._read_cells(surfaces)
universes = self._read_universes(cells)
lattices = self._read_lattices(universes)
self._finalize_geometry(cells, cell_fills, universes, lattices)
def _read_materials(self):
self.n_materials = self._f['n_materials'].value
for group in self._f['materials'].values():
material = openmc.Material.from_hdf5(group)
# Initialize dictionary for each Material
# Keys - Material keys
# Values - Material objects
self.materials = {}
for key in self._f['materials'].keys():
if key == 'n_materials':
continue
material_id = int(key.lstrip('material '))
index = self._f['materials'][key]['index'].value
name = self._f['materials'][key]['name'].value.decode()
density = self._f['materials'][key]['atom_density'].value
nuc_densities = self._f['materials'][key]['nuclide_densities'][...]
nuclides = self._f['materials'][key]['nuclides'].value
# Create the Material
material = openmc.Material(material_id=material_id, name=name)
# Read the names of the S(a,b) tables for this Material and add them
if 'sab_names' in self._f['materials'][key]:
sab_tables = self._f['materials'][key]['sab_names'].value
for sab_table in sab_tables:
name = sab_table.decode()
material.add_s_alpha_beta(name)
# Set the Material's density to atom/b-cm as used by OpenMC
material.set_density(density=density, units='atom/b-cm')
# Add all nuclides to the Material
for fullname, density in zip(nuclides, nuc_densities):
name = fullname.decode().strip()
if 'nat' in name:
material.add_element(openmc.Element(name=name),
percent=density, percent_type='ao')
else:
material.add_nuclide(openmc.Nuclide(name=name),
percent=density, percent_type='ao')
# Add the Material to the global dictionary of all Materials
self.materials[index] = material
# Add the material to the Materials collection
self.materials.append(material)
def _read_surfaces(self):
self.n_surfaces = self._f['geometry/n_surfaces'].value
surfaces = {}
for group in self._f['geometry/surfaces'].values():
surface = openmc.Surface.from_hdf5(group)
surfaces[surface.id] = surface
# Initialize dictionary for each Surface
# Keys - Surface keys
# Values - Surfacee objects
self.surfaces = {}
return surfaces
for key in self._f['geometry/surfaces'].keys():
if key == 'n_surfaces':
continue
surface_id = int(key.lstrip('surface '))
index = self._f['geometry/surfaces'][key]['index'].value
name = self._f['geometry/surfaces'][key]['name'].value.decode()
surf_type = self._f['geometry/surfaces'][key]['type'].value.decode()
bc = self._f['geometry/surfaces'][key]['boundary_condition'].value.decode()
coeffs = self._f['geometry/surfaces'][key]['coefficients'][...]
# Create the Surface based on its type
if surf_type == 'x-plane':
x0 = coeffs[0]
surface = openmc.XPlane(surface_id, bc, x0, name)
elif surf_type == 'y-plane':
y0 = coeffs[0]
surface = openmc.YPlane(surface_id, bc, y0, name)
elif surf_type == 'z-plane':
z0 = coeffs[0]
surface = openmc.ZPlane(surface_id, bc, z0, name)
elif surf_type == 'plane':
A = coeffs[0]
B = coeffs[1]
C = coeffs[2]
D = coeffs[3]
surface = openmc.Plane(surface_id, bc, A, B, C, D, name)
elif surf_type == 'x-cylinder':
y0 = coeffs[0]
z0 = coeffs[1]
R = coeffs[2]
surface = openmc.XCylinder(surface_id, bc, y0, z0, R, name)
elif surf_type == 'y-cylinder':
x0 = coeffs[0]
z0 = coeffs[1]
R = coeffs[2]
surface = openmc.YCylinder(surface_id, bc, x0, z0, R, name)
elif surf_type == 'z-cylinder':
x0 = coeffs[0]
y0 = coeffs[1]
R = coeffs[2]
surface = openmc.ZCylinder(surface_id, bc, x0, y0, R, name)
elif surf_type == 'sphere':
x0 = coeffs[0]
y0 = coeffs[1]
z0 = coeffs[2]
R = coeffs[3]
surface = openmc.Sphere(surface_id, bc, x0, y0, z0, R, name)
elif surf_type in ['x-cone', 'y-cone', 'z-cone']:
x0 = coeffs[0]
y0 = coeffs[1]
z0 = coeffs[2]
R2 = coeffs[3]
if surf_type == 'x-cone':
surface = openmc.XCone(surface_id, bc, x0, y0, z0, R2, name)
if surf_type == 'y-cone':
surface = openmc.YCone(surface_id, bc, x0, y0, z0, R2, name)
if surf_type == 'z-cone':
surface = openmc.ZCone(surface_id, bc, x0, y0, z0, R2, name)
elif surf_type == 'quadric':
a, b, c, d, e, f, g, h, j, k = coeffs
surface = openmc.Quadric(surface_id, bc, a, b, c, d, e, f,
g, h, j, k, name)
# Add Surface to global dictionary of all Surfaces
self.surfaces[index] = surface
def _read_cells(self):
self.n_cells = self._f['geometry/n_cells'].value
# Initialize dictionary for each Cell
# Keys - Cell keys
# Values - Cell objects
self.cells = {}
def _read_cells(self, surfaces):
cells = {}
# Initialize dictionary for each Cell's fill
# (e.g., Material, Universe or Lattice ID)
# This dictionary is used later to link the fills with
# the corresponding objects
# Keys - Cell keys
# Values - Filling Material, Universe or Lattice ID
self._cell_fills = {}
for key in self._f['geometry/cells'].keys():
if key == 'n_cells':
continue
cell_fills = {}
for key, group in self._f['geometry/cells'].items():
cell_id = int(key.lstrip('cell '))
index = self._f['geometry/cells'][key]['index'].value
name = self._f['geometry/cells'][key]['name'].value.decode()
fill_type = self._f['geometry/cells'][key]['fill_type'].value.decode()
name = group['name'].value.decode()
fill_type = group['fill_type'].value.decode()
if fill_type == 'normal':
fill = self._f['geometry/cells'][key]['material'].value
if fill_type == 'material':
fill = group['material'].value
elif fill_type == 'universe':
fill = self._f['geometry/cells'][key]['fill'].value
fill = group['fill'].value
else:
fill = self._f['geometry/cells'][key]['lattice'].value
fill = group['lattice'].value
if 'region' in self._f['geometry/cells'][key].keys():
region = self._f['geometry/cells'][key]['region'].value.decode()
if 'region' in group.keys():
region = group['region'].value.decode()
else:
region = []
@ -261,283 +124,79 @@ class Summary(object):
cell = openmc.Cell(cell_id=cell_id, name=name)
if fill_type == 'universe':
if 'offset' in self._f['geometry/cells'][key]:
offset = self._f['geometry/cells'][key]['offset'][...]
if 'offset' in group:
offset = group['offset'][...]
cell.offsets = offset
if 'translation' in self._f['geometry/cells'][key]:
translation = \
self._f['geometry/cells'][key]['translation'][...]
if 'translation' in group:
translation = group['translation'][...]
translation = np.asarray(translation, dtype=np.float64)
cell.translation = translation
if 'rotation' in self._f['geometry/cells'][key]:
rotation = \
self._f['geometry/cells'][key]['rotation'][...]
if 'rotation' in group:
rotation = group['rotation'][...]
rotation = np.asarray(rotation, dtype=np.int)
cell._rotation = rotation
elif fill_type == 'normal':
cell.temperature = \
self._f['geometry/cells'][key]['temperature'][...]
elif fill_type == 'material':
cell.temperature = group['temperature'][...]
# Store Cell fill information for after Universe/Lattice creation
self._cell_fills[index] = (fill_type, fill)
cell_fills[cell.id] = (fill_type, fill)
# Generate Region object given infix expression
if region:
cell.region = Region.from_expression(
region, {s.id: s for s in self.surfaces.values()})
cell.region = Region.from_expression(region, surfaces)
# Get the distribcell data
if 'distribcell_index' in self._f['geometry/cells'][key]:
ind = self._f['geometry/cells'][key]['distribcell_index'].value
if 'distribcell_index' in group:
ind = group['distribcell_index'].value
cell.distribcell_index = ind
paths = self._f['geometry/cells'][key]['paths'][...]
paths = group['paths'][...]
paths = [str(path.decode()) for path in paths]
cell.distribcell_paths = paths
# Add the Cell to the global dictionary of all Cells
self.cells[index] = cell
cells[cell.id] = cell
def _read_universes(self):
self.n_universes = self._f['geometry/n_universes'].value
return cells, cell_fills
# Initialize dictionary for each Universe
# Keys - Universe keys
# Values - Universe objects
self.universes = {}
def _read_universes(self, cells):
universes = {}
for group in self._f['geometry/universes'].values():
universe = openmc.Universe.from_hdf5(group, cells)
universes[universe.id] = universe
return universes
for key in self._f['geometry/universes'].keys():
if key == 'n_universes':
continue
def _read_lattices(self, universes):
lattices = {}
for group in self._f['geometry/lattices'].values():
lattice = openmc.Lattice.from_hdf5(group, universes)
lattices[lattice.id] = lattice
return lattices
universe_id = int(key.lstrip('universe '))
index = self._f['geometry/universes'][key]['index'].value
cells = self._f['geometry/universes'][key]['cells'][...]
# Create this Universe
universe = openmc.Universe(universe_id=universe_id)
# Add each Cell to the Universe
for cell_id in cells:
cell = self.cells[cell_id]
universe.add_cell(cell)
# Add the Universe to the global list of Universes
self.universes[index] = universe
def _read_lattices(self):
self.n_lattices = self._f['geometry/n_lattices'].value
# Initialize lattices for each Lattice
# Keys - Lattice keys
# Values - Lattice objects
self.lattices = {}
for key in self._f['geometry/lattices'].keys():
if key == 'n_lattices':
continue
lattice_id = int(key.lstrip('lattice '))
index = self._f['geometry/lattices'][key]['index'].value
name = self._f['geometry/lattices'][key]['name'].value.decode()
lattice_type = self._f['geometry/lattices'][key]['type'].value.decode()
if 'offsets' in self._f['geometry/lattices'][key]:
offsets = self._f['geometry/lattices'][key]['offsets'][...]
else:
offsets = None
if lattice_type == 'rectangular':
dimension = self._f['geometry/lattices'][key]['dimension'][...]
lower_left = \
self._f['geometry/lattices'][key]['lower_left'][...]
pitch = self._f['geometry/lattices'][key]['pitch'][...]
outer = self._f['geometry/lattices'][key]['outer'].value
universe_ids = \
self._f['geometry/lattices'][key]['universes'][...]
# Create the Lattice
lattice = openmc.RectLattice(lattice_id=lattice_id, name=name)
lattice.lower_left = lower_left
lattice.pitch = pitch
# If the Universe specified outer the Lattice is not void (-22)
if outer != -22:
lattice.outer = self.universes[outer]
# Build array of Universe pointers for the Lattice
universes = \
np.empty(tuple(universe_ids.shape), dtype=openmc.Universe)
for z in range(universe_ids.shape[0]):
for y in range(universe_ids.shape[1]):
for x in range(universe_ids.shape[2]):
universes[z, y, x] = \
self.get_universe_by_id(universe_ids[z, y, x])
# Use 2D NumPy array to store lattice universes for 2D lattices
if len(dimension) == 2:
universes = np.squeeze(universes)
universes = np.atleast_2d(universes)
# Set the universes for the lattice
lattice.universes = universes
# Set the distribcell offsets for the lattice
if offsets is not None:
lattice.offsets = offsets
# Add the Lattice to the global dictionary of all Lattices
self.lattices[index] = lattice
if lattice_type == 'hexagonal':
n_rings = self._f['geometry/lattices'][key]['n_rings'].value
n_axial = self._f['geometry/lattices'][key]['n_axial'].value
center = self._f['geometry/lattices'][key]['center'][...]
pitch = self._f['geometry/lattices'][key]['pitch'][...]
outer = self._f['geometry/lattices'][key]['outer'].value
universe_ids = self._f[
'geometry/lattices'][key]['universes'][...]
# Create the Lattice
lattice = openmc.HexLattice(lattice_id=lattice_id, name=name)
lattice.center = center
lattice.pitch = pitch
# If the Universe specified outer the Lattice is not void (-22)
if outer != -22:
lattice.outer = self.universes[outer]
# Build array of Universe pointers for the Lattice. Note that
# we need to convert between the HDF5's square array of
# (x, alpha, z) to the Python API's format of a ragged nested
# list of (z, ring, theta).
universes = []
for z in range(n_axial):
# Add a list for this axial level.
universes.append([])
x = n_rings - 1
a = 2*n_rings - 2
for r in range(n_rings - 1, 0, -1):
# Add a list for this ring.
universes[-1].append([])
# Climb down the top-right.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
x += 1
a -= 1
# Climb down the right.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
a -= 1
# Climb down the bottom-right.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
x -= 1
# Climb up the bottom-left.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
x -= 1
a += 1
# Climb up the left.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
a += 1
# Climb up the top-left.
for i in range(r):
universes[-1][-1].append(universe_ids[z, a, x])
x += 1
# Move down to the next ring.
a -= 1
# Convert the ids into Universe objects.
universes[-1][-1] = [self.get_universe_by_id(u_id)
for u_id in universes[-1][-1]]
# Handle the degenerate center ring separately.
u_id = universe_ids[z, a, x]
universes[-1].append([self.get_universe_by_id(u_id)])
# Add the universes to the lattice.
if len(pitch) == 2:
# Lattice is 3D
lattice.universes = universes
else:
# Lattice is 2D; extract the only axial level
lattice.universes = universes[0]
if offsets is not None:
lattice.offsets = offsets
# Add the Lattice to the global dictionary of all Lattices
self.lattices[index] = lattice
def _finalize_geometry(self):
# Initialize Geometry object
self._openmc_geometry = openmc.Geometry()
def _finalize_geometry(self, cells, cell_fills, universes, lattices):
materials = {m.id: m for m in self.materials}
# Iterate over all Cells and add fill Materials, Universes and Lattices
for cell_key in self._cell_fills.keys():
# Determine fill type ('normal', 'universe', or 'lattice') and ID
fill_type = self._cell_fills[cell_key][0]
fill_id = self._cell_fills[cell_key][1]
for cell_id, (fill_type, fill_id) in cell_fills.items():
# Retrieve the object corresponding to the fill type and ID
if fill_type == 'normal':
if fill_type == 'material':
if isinstance(fill_id, Iterable):
fill = [self.get_material_by_id(mat) if mat > 0 else None
fill = [materials[mat] if mat > 0 else None
for mat in fill_id]
else:
if fill_id > 0:
fill = self.get_material_by_id(fill_id)
else:
fill = None
fill = materials[fill_id] if fill_id > 0 else None
elif fill_type == 'universe':
fill = self.get_universe_by_id(fill_id)
fill = universes[fill_id]
else:
fill = self.get_lattice_by_id(fill_id)
fill = lattices[fill_id]
# Set the fill for the Cell
self.cells[cell_key].fill = fill
cells[cell_id].fill = fill
# Set the root universe for the Geometry
root_universe = self.get_universe_by_id(0)
self.openmc_geometry.root_universe = root_universe
def _read_tallies(self):
# Initialize a dictionary for the tally names
# Keys - Tally IDs
# Values - Tally names
self.tally_names = {}
# Read the number of tallies
if 'tallies' not in self._f:
return
# OpenMC Tally keys
all_keys = self._f['tallies/'].keys()
tally_keys = [key for key in all_keys if 'tally' in key]
base = 'tallies/tally '
# Iterate over all Tallies
for tally_key in tally_keys:
tally_id = int(tally_key.strip('tally '))
subbase = '{0}{1}'.format(base, tally_id)
# Read Tally name metadata
tally_name = self._f['{0}/name'.format(subbase)].value.decode()
self.tally_names[tally_id] = tally_name
self.geometry.root_universe = universes[0]
def add_volume_information(self, volume_calc):
"""Add volume information to the geometry within the summary file
@ -548,109 +207,4 @@ class Summary(object):
Results from a stochastic volume calculation
"""
self.openmc_geometry.add_volume_information(volume_calc)
def get_material_by_id(self, material_id):
"""Return a Material object given the material id
Parameters
----------
id : int
Unique identifier for the material
Returns
-------
material : openmc.Material
Material with given id
"""
for material in self.materials.values():
if material.id == material_id:
return material
return None
def get_surface_by_id(self, surface_id):
"""Return a Surface object given the surface id
Parameters
----------
id : int
Unique identifier for the surface
Returns
-------
surface : openmc.Surface
Surface with given id
"""
for surface in self.surfaces.values():
if surface.id == surface_id:
return surface
return None
def get_cell_by_id(self, cell_id):
"""Return a Cell object given the cell id
Parameters
----------
id : int
Unique identifier for the cell
Returns
-------
cell : openmc.Cell
Cell with given id
"""
for cell in self.cells.values():
if cell.id == cell_id:
return cell
return None
def get_universe_by_id(self, universe_id):
"""Return a Universe object given the universe id
Parameters
----------
id : int
Unique identifier for the universe
Returns
-------
universe : openmc.Universe
Universe with given id
"""
for universe in self.universes.values():
if universe.id == universe_id:
return universe
return None
def get_lattice_by_id(self, lattice_id):
"""Return a Lattice object given the lattice id
Parameters
----------
id : int
Unique identifier for the lattice
Returns
-------
lattice : openmc.Lattice
Lattice with given id
"""
for lattice in self.lattices.values():
if lattice.id == lattice_id:
return lattice
return None
self.geometry.add_volume_information(volume_calc)

96
openmc/surface.py
View file

@ -14,7 +14,7 @@ from openmc.region import Region, Intersection, Union
# A static variable for auto-generated Surface IDs
AUTO_SURFACE_ID = 10000
_BC_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic']
_BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic']
def reset_auto_surface_id():
@ -141,7 +141,7 @@ class Surface(object):
@boundary_type.setter
def boundary_type(self, boundary_type):
check_type('boundary type', boundary_type, string_types)
check_value('boundary type', boundary_type, _BC_TYPES)
check_value('boundary type', boundary_type, _BOUNDARY_TYPES)
self._boundary_type = boundary_type
def bounding_box(self, side):
@ -171,7 +171,15 @@ class Surface(object):
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
def create_xml_subelement(self):
def to_xml_element(self):
"""Return XML representation of the surface
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing source data
"""
element = ET.Element("surface")
element.set("id", str(self._id))
@ -186,6 +194,76 @@ class Surface(object):
return element
@staticmethod
def from_hdf5(group):
"""Create surface from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
Returns
-------
openmc.Surface
Instance of surface subclass
"""
surface_id = int(group.name.split('/')[-1].lstrip('surface '))
name = group['name'].value.decode()
surf_type = group['type'].value.decode()
bc = group['boundary_type'].value.decode()
coeffs = group['coefficients'][...]
# Create the Surface based on its type
if surf_type == 'x-plane':
x0 = coeffs[0]
surface = XPlane(surface_id, bc, x0, name)
elif surf_type == 'y-plane':
y0 = coeffs[0]
surface = YPlane(surface_id, bc, y0, name)
elif surf_type == 'z-plane':
z0 = coeffs[0]
surface = ZPlane(surface_id, bc, z0, name)
elif surf_type == 'plane':
A, B, C, D = coeffs
surface = Plane(surface_id, bc, A, B, C, D, name)
elif surf_type == 'x-cylinder':
y0, z0, R = coeffs
surface = XCylinder(surface_id, bc, y0, z0, R, name)
elif surf_type == 'y-cylinder':
x0, z0, R = coeffs
surface = YCylinder(surface_id, bc, x0, z0, R, name)
elif surf_type == 'z-cylinder':
x0, y0, R = coeffs
surface = ZCylinder(surface_id, bc, x0, y0, R, name)
elif surf_type == 'sphere':
x0, y0, z0, R = coeffs
surface = Sphere(surface_id, bc, x0, y0, z0, R, name)
elif surf_type in ['x-cone', 'y-cone', 'z-cone']:
x0, y0, z0, R2 = coeffs
if surf_type == 'x-cone':
surface = XCone(surface_id, bc, x0, y0, z0, R2, name)
elif surf_type == 'y-cone':
surface = YCone(surface_id, bc, x0, y0, z0, R2, name)
elif surf_type == 'z-cone':
surface = ZCone(surface_id, bc, x0, y0, z0, R2, name)
elif surf_type == 'quadric':
a, b, c, d, e, f, g, h, j, k = coeffs
surface = Quadric(surface_id, bc, a, b, c, d, e, f, g,
h, j, k, name)
return surface
class Plane(Surface):
"""An arbitrary plane of the form :math:`Ax + By + Cz = D`.
@ -314,8 +392,16 @@ class Plane(Surface):
x, y, z = point
return self.a*x + self.b*y + self.c*z - self.d
def create_xml_subelement(self):
element = super(Plane, self).create_xml_subelement()
def to_xml_element(self):
"""Return XML representation of the surface
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing source data
"""
element = super(Plane, self).to_xml_element()
# Add periodic surface pair information
if self.boundary_type == 'periodic':

2
openmc/tallies.py
View file

@ -195,7 +195,7 @@ class Tally(object):
if isinstance(nuclide, openmc.Nuclide):
string += nuclide.name + ' '
else:
string += nuclide + ' '
string += str(nuclide) + ' '
string += '\n'

40
openmc/universe.py
View file

@ -118,6 +118,35 @@ class Universe(object):
else:
self._name = ''
@classmethod
def from_hdf5(cls, group, cells):
"""Create universe from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
cells : dict
Dictionary mapping cell IDs to instances of :class:`openmc.Cell`.
Returns
-------
openmc.Universe
Universe instance
"""
universe_id = int(group.name.split('/')[-1].lstrip('universe '))
cell_ids = group['cells'].value
# Create this Universe
universe = cls(universe_id)
# Add each Cell to the Universe
for cell_id in cell_ids:
universe.add_cell(cells[cell_id])
return universe
def find(self, point):
"""Find cells/universes/lattices which contain a given point
@ -411,7 +440,7 @@ class Universe(object):
Returns
-------
materials : Collections.OrderedDict
materials : collections.OrderedDict
Dictionary whose keys are material IDs and values are
:class:`Material` instances
@ -436,14 +465,9 @@ class Universe(object):
:class:`Universe` instances
"""
# Get all Cells in this Universe
cells = self.get_all_cells()
# Append all Universes within each Cell to the dictionary
universes = OrderedDict()
# Append all Universes containing each Cell to the dictionary
for cell in cells.values():
for cell in self.get_all_cells().values():
universes.update(cell.get_all_universes())
return universes

4
openmc/volume.py
View file

@ -10,6 +10,8 @@ import h5py
import openmc
import openmc.checkvalue as cv
_VERSION_VOLUME = 1
class VolumeCalculation(object):
"""Stochastic volume calculation specifications and results.
@ -189,6 +191,8 @@ class VolumeCalculation(object):
"""
with h5py.File(filename, 'r') as f:
cv.check_filetype_version(f, "volume", _VERSION_VOLUME)
domain_type = f.attrs['domain_type'].decode()
samples = f.attrs['samples']
lower_left = f.attrs['lower_left']

8
scripts/openmc-voxel-to-silovtk
View file

@ -27,9 +27,9 @@ def parse_options():
def main(filename, o):
# Read data from voxel file
fh = h5py.File(filename, 'r')
dimension = fh['num_voxels'].value
width = fh['voxel_width'].value
lower_left = fh['lower_left'].value
dimension = fh.attrs['num_voxels']
width = fh.attrs['voxel_width']
lower_left = fh.attrs['lower_left']
voxel_data = fh['data'].value
nx, ny, nz = dimension
@ -40,7 +40,7 @@ def main(filename, o):
import vtk
except:
print('The vtk python bindings do not appear to be installed '
'properly.\nOn Ubuntu: sudo apt-get install python-vtk\n'
'properly.\nOn Ubuntu: sudo apt install python-vtk\n'
'See: http://www.vtk.org/')
return

22
src/cmfd_execute.F90
View file

@ -22,6 +22,7 @@ contains
use cmfd_data, only: set_up_cmfd
use cmfd_solver, only: cmfd_solver_execute
use error, only: warning, fatal_error
use message_passing, only: master
! CMFD single processor on master
if (master) then
@ -90,13 +91,9 @@ contains
subroutine calc_fission_source()
use constants, only: CMFD_NOACCEL, ZERO, TWO
use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch
use string, only: to_str
#ifdef MPI
use global, only: mpi_err
use global, only: cmfd, cmfd_coremap, entropy_on, current_batch
use message_passing
#endif
use string, only: to_str
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
@ -202,7 +199,7 @@ contains
#ifdef MPI
! Broadcast full source to all procs
call MPI_BCAST(cmfd % cmfd_src, n, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
call MPI_BCAST(cmfd % cmfd_src, n, MPI_REAL8, 0, mpi_intracomm, mpi_err)
#endif
end subroutine calc_fission_source
@ -216,16 +213,11 @@ contains
use algorithm, only: binary_search
use constants, only: ZERO, ONE
use error, only: warning, fatal_error
use global, only: meshes, source_bank, work, n_user_meshes, cmfd, &
master
use global, only: meshes, source_bank, work, n_user_meshes, cmfd
use mesh_header, only: RegularMesh
use mesh, only: count_bank_sites, get_mesh_indices
use string, only: to_str
#ifdef MPI
use global, only: mpi_err
use message_passing
#endif
use string, only: to_str
logical, intent(in) :: new_weights ! calcualte new weights
@ -289,7 +281,7 @@ contains
! Broadcast weight factors to all procs
#ifdef MPI
call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, &
MPI_COMM_WORLD, mpi_err)
mpi_intracomm, mpi_err)
#endif
end if

1
src/cmfd_input.F90
View file

@ -15,6 +15,7 @@ contains
subroutine configure_cmfd()
use cmfd_header, only: allocate_cmfd
use message_passing, only: master
integer :: color ! color group of processor

6
src/cmfd_solver.F90
View file

@ -301,10 +301,12 @@ contains
subroutine convergence(iter, innerits, iconv)
use constants, only: ONE, ZERO
use global, only: cmfd_power_monitor, master
use, intrinsic :: ISO_FORTRAN_ENV
use constants, only: ONE, ZERO
use global, only: cmfd_power_monitor
use message_passing, only: master
integer, intent(in) :: iter ! outer iteration number
integer, intent(in) :: innerits ! inner iteration nubmer
logical, intent(out) :: iconv ! convergence logical

32
src/constants.F90
View file

@ -9,17 +9,20 @@ module constants
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 8
integer, parameter :: VERSION_RELEASE = 0
integer, parameter :: &
VERSION(3) = [VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE]
! HDF5 data format
integer, parameter :: HDF5_VERSION_MAJOR = 1
integer, parameter :: HDF5_VERSION_MINOR = 0
integer, parameter :: HDF5_VERSION(2) = [1, 0]
! Revision numbers for binary files
integer, parameter :: REVISION_STATEPOINT = 15
integer, parameter :: REVISION_PARTICLE_RESTART = 1
integer, parameter :: REVISION_TRACK = 1
integer, parameter :: REVISION_SUMMARY = 4
character(10), parameter :: MULTIPOLE_VERSION = "v0.2"
! Version numbers for binary files
integer, parameter :: VERSION_STATEPOINT(2) = [16, 0]
integer, parameter :: VERSION_PARTICLE_RESTART(2) = [2, 0]
integer, parameter :: VERSION_TRACK(2) = [2, 0]
integer, parameter :: VERSION_SUMMARY(2) = [5, 0]
integer, parameter :: VERSION_VOLUME(2) = [1, 0]
integer, parameter :: VERSION_VOXEL(2) = [1, 0]
character(10), parameter :: VERSION_MULTIPOLE = "v0.2"
! ============================================================================
! ADJUSTABLE PARAMETERS
@ -104,11 +107,11 @@ module constants
OP_INTERSECTION = huge(0) - 3, & ! Intersection operator
OP_UNION = huge(0) - 4 ! Union operator (^)
! Cell types
! Cell fill types
integer, parameter :: &
CELL_NORMAL = 1, & ! Cell with a specified material
CELL_FILL = 2, & ! Cell filled by a separate universe
CELL_LATTICE = 3 ! Cell filled with a lattice
FILL_MATERIAL = 1, & ! Cell with a specified material
FILL_UNIVERSE = 2, & ! Cell filled by a separate universe
FILL_LATTICE = 3 ! Cell filled with a lattice
! Void material
integer, parameter :: MATERIAL_VOID = -1
@ -142,7 +145,7 @@ module constants
SURF_CONE_Z = 11 ! Cone parallel to z-axis
! Flag to say that the outside of a lattice is not defined
integer, parameter :: NO_OUTER_UNIVERSE = -22
integer, parameter :: NO_OUTER_UNIVERSE = -1
! Maximum number of lost particles
integer, parameter :: MAX_LOST_PARTICLES = 10
@ -424,7 +427,8 @@ module constants
MODE_FIXEDSOURCE = 1, & ! Fixed source mode
MODE_EIGENVALUE = 2, & ! K eigenvalue mode
MODE_PLOTTING = 3, & ! Plotting mode
MODE_PARTICLE = 4 ! Particle restart mode
MODE_PARTICLE = 4, & ! Particle restart mode
MODE_VOLUME = 5 ! Volume calculation mode
!=============================================================================
! CMFD CONSTANTS

187
src/eigenvalue.F90
View file

@ -1,8 +1,5 @@
module eigenvalue
#ifdef MPI
use message_passing
#endif
use algorithm, only: binary_search
use constants, only: ZERO
@ -11,6 +8,7 @@ module eigenvalue
use math, only: t_percentile
use mesh, only: count_bank_sites
use mesh_header, only: RegularMesh
use message_passing
use random_lcg, only: prn, set_particle_seed, advance_prn_seed
use string, only: to_str
@ -66,7 +64,7 @@ contains
#ifdef MPI
start = 0_8
call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, &
MPI_COMM_WORLD, mpi_err)
mpi_intracomm, mpi_err)
! While we would expect the value of start on rank 0 to be 0, the MPI
! standard says that the receive buffer on rank 0 is undefined and not
@ -76,7 +74,7 @@ contains
finish = start + n_bank
total = finish
call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, &
MPI_COMM_WORLD, mpi_err)
mpi_intracomm, mpi_err)
#else
start = 0_8
@ -150,13 +148,13 @@ contains
! 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, &
MPI_COMM_WORLD, mpi_err)
mpi_intracomm, mpi_err)
finish = start + index_temp
! Allocate space for bank_position if this hasn't been done yet
if (.not. allocated(bank_position)) allocate(bank_position(n_procs))
call MPI_ALLGATHER(start, 1, MPI_INTEGER8, bank_position, 1, &
MPI_INTEGER8, MPI_COMM_WORLD, mpi_err)
MPI_INTEGER8, mpi_intracomm, mpi_err)
#else
start = 0_8
finish = index_temp
@ -210,7 +208,7 @@ contains
if (neighbor /= rank) then
n_request = n_request + 1
call MPI_ISEND(temp_sites(index_local), int(n), MPI_BANK, neighbor, &
rank, MPI_COMM_WORLD, request(n_request), mpi_err)
rank, mpi_intracomm, request(n_request), mpi_err)
end if
! Increment all indices
@ -254,7 +252,7 @@ contains
n_request = n_request + 1
call MPI_IRECV(source_bank(index_local), int(n), MPI_BANK, &
neighbor, neighbor, MPI_COMM_WORLD, request(n_request), mpi_err)
neighbor, neighbor, mpi_intracomm, request(n_request), mpi_err)
else
! If the source sites are on this procesor, we can simply copy them
@ -379,7 +377,7 @@ contains
#ifdef MPI
! Combine values across all processors
call MPI_ALLREDUCE(keff_generation, k_generation(overall_gen), 1, &
MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, mpi_err)
MPI_REAL8, MPI_SUM, mpi_intracomm, mpi_err)
#else
k_generation(overall_gen) = keff_generation
#endif
@ -448,7 +446,7 @@ contains
integer :: l ! loop index
integer :: i, j, k ! indices referring to collision, absorption, or track
integer :: n ! number of realizations
real(8) :: n ! number of realizations
real(8) :: kv(3) ! vector of k-effective estimates
real(8) :: cov(3,3) ! sample covariance matrix
real(8) :: f ! weighting factor
@ -460,73 +458,132 @@ contains
if (n_realizations <= 3) return
! Initialize variables
n = n_realizations
g = ZERO
S = ZERO
k_combined = ZERO
n = real(n_realizations, 8)
! Copy estimates of k-effective and its variance (not variance of the mean)
kv(1) = global_tallies(RESULT_SUM, K_COLLISION) / n
kv(2) = global_tallies(RESULT_SUM, K_ABSORPTION) / n
kv(3) = global_tallies(RESULT_SUM, K_TRACKLENGTH) / n
cov(1,1) = (global_tallies(RESULT_SUM_SQ, K_COLLISION) - &
n * kv(1) * kv(1)) / (n - 1)
cov(2,2) = (global_tallies(RESULT_SUM_SQ, K_ABSORPTION) - &
n * kv(2) * kv(2)) / (n - 1)
cov(3,3) = (global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) - &
n * kv(3) * kv(3)) / (n - 1)
cov(1, 1) = (global_tallies(RESULT_SUM_SQ, K_COLLISION) - &
n * kv(1) * kv(1)) / (n - ONE)
cov(2, 2) = (global_tallies(RESULT_SUM_SQ, K_ABSORPTION) - &
n * kv(2) * kv(2)) / (n - ONE)
cov(3, 3) = (global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) - &
n * kv(3) * kv(3)) / (n - ONE)
! Calculate covariances based on sums with Bessel's correction
cov(1,2) = (k_col_abs - n * kv(1) * kv(2))/(n - 1)
cov(1,3) = (k_col_tra - n * kv(1) * kv(3))/(n - 1)
cov(2,3) = (k_abs_tra - n * kv(2) * kv(3))/(n - 1)
cov(2,1) = cov(1,2)
cov(3,1) = cov(1,3)
cov(3,2) = cov(2,3)
cov(1, 2) = (k_col_abs - n * kv(1) * kv(2)) / (n - ONE)
cov(1, 3) = (k_col_tra - n * kv(1) * kv(3)) / (n - ONE)
cov(2, 3) = (k_abs_tra - n * kv(2) * kv(3)) / (n - ONE)
cov(2, 1) = cov(1, 2)
cov(3, 1) = cov(1, 3)
cov(3, 2) = cov(2, 3)
do l = 1, 3
! Permutations of estimates
if (l == 1) then
! i = collision, j = absorption, k = tracklength
i = 1
j = 2
k = 3
elseif (l == 2) then
! i = absortion, j = tracklength, k = collision
i = 2
j = 3
k = 1
elseif (l == 3) then
! i = tracklength, j = collision, k = absorption
i = 3
j = 1
k = 2
end if
! Check to see if two estimators are the same; this is guaranteed to happen
! in MG-mode with survival biasing when the collision and absorption
! estimators are the same, but can theoretically happen at anytime.
! If it does, the standard estimators will produce floating-point
! exceptions and an expression specifically derived for the combination of
! two estimators (vice three) should be used instead.
! Calculate weighting
f = cov(j,j)*(cov(k,k) - cov(i,k)) - cov(k,k)*cov(i,j) + &
cov(j,k)*(cov(i,j) + cov(i,k) - cov(j,k))
! First we will identify if there are any matching estimators
if ((abs(kv(1) - kv(2)) / kv(1) < FP_REL_PRECISION) .and. &
(abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) then
! 1 and 2 match, so only use 1 and 3 in our comparisons
i = 1
j = 3
! Add to S sums for variance of combined estimate
S(1) = S(1) + f * cov(1,l)
S(2) = S(2) + (cov(j,j) + cov(k,k) - TWO*cov(j,k))*kv(l)*kv(l)
S(3) = S(3) + (cov(k,k) + cov(i,j) - cov(j,k) - cov(i,k))*kv(l)*kv(j)
else if ((abs(kv(1) - kv(3)) / kv(1) < FP_REL_PRECISION) .and. &
(abs(cov(1, 1) - cov(3, 3)) / cov(1, 1) < FP_REL_PRECISION)) then
! 1 and 3 match, so only use 1 and 2 in our comparisons
i = 1
j = 2
! Add to sum for combined k-effective
k_combined(1) = k_combined(1) + f * kv(l)
g = g + f
end do
else if ((abs(kv(2) - kv(3)) / kv(2) < FP_REL_PRECISION) .and. &
(abs(cov(2, 2) - cov(3, 3)) / cov(2, 2) < FP_REL_PRECISION)) then
! 2 and 3 match, so only use 1 and 2 in our comparisons
i = 1
j = 2
! Complete calculations of S sums
S = (n - 1)*S
S(1) = (n - 1)**2 * S(1)
else
! No two estimators match, so set i to 0 and this will be the indicator
! to use all three estimators.
i = 0
end if
! Calculate combined estimate of k-effective
k_combined(1) = k_combined(1) / g
if (i == 0) then
! Use three estimators as derived in the paper by Urbatsch
! Calculate standard deviation of combined estimate
g = (n - 1)**2 * g
k_combined(2) = sqrt(S(1)/(g*n*(n-3)) * (ONE + n*((S(2) - TWO*S(3))/g)))
! Initialize variables
g = ZERO
S = ZERO
k_combined = ZERO
do l = 1, 3
! Permutations of estimates
if (l == 1) then
! i = collision, j = absorption, k = tracklength
i = 1
j = 2
k = 3
elseif (l == 2) then
! i = absortion, j = tracklength, k = collision
i = 2
j = 3
k = 1
elseif (l == 3) then
! i = tracklength, j = collision, k = absorption
i = 3
j = 1
k = 2
end if
! Calculate weighting
f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) + &
cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k))
! Add to S sums for variance of combined estimate
S(1) = S(1) + f * cov(1, l)
S(2) = S(2) + (cov(j, j) + cov(k, k) - TWO * cov(j, k)) * kv(l) * kv(l)
S(3) = S(3) + (cov(k, k) + cov(i, j) - cov(j, k) - &
cov(i, k)) * kv(l) * kv(j)
! Add to sum for combined k-effective
k_combined(1) = k_combined(1) + f * kv(l)
g = g + f
end do
! Complete calculations of S sums
S = (n - ONE) * S
S(1) = (n - ONE)**2 * S(1)
! Calculate combined estimate of k-effective
k_combined(1) = k_combined(1) / g
! Calculate standard deviation of combined estimate
g = (n - ONE)**2 * g
k_combined(2) = sqrt(S(1) / &
(g * n * (n - THREE)) * (ONE + n * ((S(2) - TWO * S(3)) / g)))
else
! Use only two estimators
! These equations are derived analogously to that done in the paper by
! Urbatsch, but are simpler than for the three estimators case since the
! block matrices of the three estimator equations reduces to scalars here
! Store the commonly used term
f = kv(i) - kv(j)
g = cov(i, i) + cov(j, j) - TWO * cov(i, j)
! Calculate combined estimate of k-effective
k_combined(1) = kv(i) - (cov(i, i) - cov(i, j)) / g * f
! Calculate standard deviation of combined estimate
k_combined(2) = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) * &
(g + n * f * f) / (n * (n - TWO) * g * g)
k_combined(2) = sqrt(k_combined(2))
end if
end subroutine calculate_combined_keff
@ -565,7 +622,7 @@ contains
#ifdef MPI
! Send source fraction to all processors
n = product(ufs_mesh % dimension)
call MPI_BCAST(source_frac, n, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
call MPI_BCAST(source_frac, n, MPI_REAL8, 0, mpi_intracomm, mpi_err)
#endif
! Normalize to total weight to get fraction of source in each cell

4
src/error.F90
View file

@ -3,9 +3,7 @@ module error
use, intrinsic :: ISO_FORTRAN_ENV
use constants
#ifdef MPI
use message_passing
#endif
implicit none
@ -139,7 +137,7 @@ contains
#ifdef MPI
! Abort MPI
call MPI_ABORT(MPI_COMM_WORLD, code, mpi_err)
call MPI_ABORT(mpi_intracomm, code, mpi_err)
#endif
! Abort program

60
src/finalize.F90
View file

@ -1,17 +1,11 @@
module finalize
use global
use output, only: print_runtime, print_results, &
print_overlap_check, write_tallies
use tally, only: tally_statistics
#ifdef MPI
use message_passing
#endif
use hdf5_interface, only: hdf5_bank_t
use hdf5, only: h5tclose_f, h5close_f
use global
use hdf5_interface, only: hdf5_bank_t
use message_passing
implicit none
contains
@ -21,34 +15,10 @@ contains
! statistics and writing out tallies
!===============================================================================
subroutine finalize_run()
subroutine openmc_finalize()
integer :: hdf5_err
! Start finalization timer
call time_finalize%start()
if (run_mode /= MODE_PLOTTING .and. run_mode /= MODE_PARTICLE) then
! Calculate statistics for tallies and write to tallies.out
if (master) then
if (n_realizations > 1) call tally_statistics()
end if
if (output_tallies) then
if (master) call write_tallies()
end if
if (check_overlaps) call reduce_overlap_count()
end if
! Stop timers and show timing statistics
call time_finalize%stop()
call time_total%stop()
if (master .and. (run_mode /= MODE_PLOTTING .and. &
run_mode /= MODE_PARTICLE)) then
call print_runtime()
call print_results()
if (check_overlaps) call print_overlap_check()
end if
! Deallocate arrays
call free_memory()
@ -66,24 +36,6 @@ contains
call MPI_FINALIZE(mpi_err)
#endif
end subroutine finalize_run
!===============================================================================
! REDUCE_OVERLAP_COUNT accumulates cell overlap check counts to master
!===============================================================================
subroutine reduce_overlap_count()
#ifdef MPI
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, overlap_check_cnt, n_cells, &
MPI_INTEGER8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
else
call MPI_REDUCE(overlap_check_cnt, overlap_check_cnt, n_cells, &
MPI_INTEGER8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
end if
#endif
end subroutine reduce_overlap_count
end subroutine openmc_finalize
end module finalize

32
src/geometry.F90
View file

@ -244,7 +244,7 @@ contains
call write_message(" Entering cell " // trim(to_str(c % id)))
end if
CELL_TYPE: if (c % type == CELL_NORMAL) then
CELL_TYPE: if (c % type == FILL_MATERIAL) then
! ======================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
@ -260,10 +260,10 @@ contains
distribcell_index = c % distribcell_index
offset = 0
do k = 1, p % n_coord
if (cells(p % coord(k) % cell) % type == CELL_FILL) then
if (cells(p % coord(k) % cell) % type == FILL_UNIVERSE) then
offset = offset + cells(p % coord(k) % cell) % &
offset(distribcell_index)
elseif (cells(p % coord(k) % cell) % type == CELL_LATTICE) then
elseif (cells(p % coord(k) % cell) % type == FILL_LATTICE) then
if (lattices(p % coord(k + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(k + 1) % lattice_x, &
@ -293,7 +293,7 @@ contains
p % sqrtkT = c % sqrtkT(1)
end if
elseif (c % type == CELL_FILL) then CELL_TYPE
elseif (c % type == FILL_UNIVERSE) then CELL_TYPE
! ======================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
@ -322,7 +322,7 @@ contains
j = p % n_coord
if (.not. found) exit
elseif (c % type == CELL_LATTICE) then CELL_TYPE
elseif (c % type == FILL_LATTICE) then CELL_TYPE
! ======================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
@ -1116,11 +1116,11 @@ contains
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type == CELL_NORMAL) then
if (c % type == FILL_MATERIAL) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type == CELL_FILL) then
elseif (c % type == FILL_UNIVERSE) then
! Set offset for the cell on this level
c % offset(map) = offset
@ -1134,7 +1134,7 @@ contains
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
elseif (c % type == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
@ -1232,11 +1232,11 @@ contains
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type == CELL_NORMAL) then
if (c % type == FILL_MATERIAL) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type == CELL_FILL) then
elseif (c % type == FILL_UNIVERSE) then
next_univ => universes(c % fill)
@ -1251,7 +1251,7 @@ contains
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
elseif (c % type == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
@ -1346,11 +1346,11 @@ contains
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type == CELL_NORMAL) then
if (c % type == FILL_MATERIAL) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type == CELL_FILL) then
elseif (c % type == FILL_UNIVERSE) then
next_univ => universes(c % fill)
@ -1359,7 +1359,7 @@ contains
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
elseif (c % type == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
@ -1428,14 +1428,14 @@ contains
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type == CELL_FILL) then
if (c % type == FILL_UNIVERSE) then
next_univ => universes(c % fill)
levels_below = max(levels_below, maximum_levels(next_univ))
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
elseif (c % type == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj

15
src/global.F90
View file

@ -267,21 +267,6 @@ module global
! ============================================================================
! PARALLEL PROCESSING VARIABLES
! The defaults set here for the number of processors, rank, and master and
! mpi_enabled flag are for when MPI is not being used at all, i.e. a serial
! run. In this case, these variables are still used at times.
integer :: n_procs = 1 ! number of processes
integer :: rank = 0 ! rank of process
logical :: master = .true. ! master process?
logical :: mpi_enabled = .false. ! is MPI in use and initialized?
integer :: mpi_err ! MPI error code
#ifdef MPIF08
type(MPI_Datatype) :: MPI_BANK
#else
integer :: MPI_BANK ! MPI datatype for fission bank
#endif
#ifdef _OPENMP
integer :: n_threads = NONE ! number of OpenMP threads
integer :: thread_id ! ID of a given thread

89
src/hdf5_interface.F90
View file

@ -17,7 +17,7 @@ module hdf5_interface
use error, only: fatal_error
#ifdef PHDF5
use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL
use message_passing, only: mpi_intracomm, MPI_INFO_NULL
#endif
implicit none
@ -75,6 +75,8 @@ module hdf5_interface
module procedure write_attribute_double
module procedure write_attribute_double_1D
module procedure write_attribute_integer
module procedure write_attribute_integer_1D
module procedure write_attribute_string
end interface write_attribute
public :: write_dataset
@ -93,7 +95,6 @@ module hdf5_interface
public :: close_dataset
public :: get_shape
public :: get_ndims
public :: write_attribute_string
public :: get_groups
public :: get_datasets
public :: get_name
@ -124,10 +125,10 @@ contains
call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err)
#ifdef PHDF5
#ifdef MPIF08
call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD%MPI_VAL, &
call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, &
MPI_INFO_NULL%MPI_VAL, hdf5_err)
#else
call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD, MPI_INFO_NULL, hdf5_err)
call h5pset_fapl_mpio_f(plist, mpi_intracomm, MPI_INFO_NULL, hdf5_err)
#endif
#endif
@ -174,10 +175,10 @@ contains
call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err)
#ifdef PHDF5
#ifdef MPIF08
call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD%MPI_VAL, &
call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, &
MPI_INFO_NULL%MPI_VAL, hdf5_err)
#else
call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD, MPI_INFO_NULL, hdf5_err)
call h5pset_fapl_mpio_f(plist, mpi_intracomm, MPI_INFO_NULL, hdf5_err)
#endif
#endif
@ -2046,15 +2047,46 @@ contains
! WRITE_ATTRIBUTE_STRING
!===============================================================================
subroutine write_attribute_string(group_id, var, attr_type, attr_str)
integer(HID_T), intent(in) :: group_id
character(*), intent(in) :: var ! variable name for attr
character(*), intent(in) :: attr_type ! attr identifier type
character(*), intent(in) :: attr_str ! string for attr id type
subroutine write_attribute_string(obj_id, name, buffer)
integer(HID_T), intent(in) :: obj_id ! object to write attribute to
character(*), intent(in) :: name ! name of attribute
character(*), intent(in), target :: buffer ! string to write
integer :: hdf5_err
integer :: hdf5_err
integer(HID_T) :: dspace_id
integer(HID_T) :: attr_id
integer(HID_T) :: filetype
integer(SIZE_T) :: i
integer(SIZE_T) :: n
character(kind=C_CHAR), allocatable, target :: temp_buffer(:)
type(c_ptr) :: f_ptr
call h5ltset_attribute_string_f(group_id, var, attr_type, attr_str, hdf5_err)
! Create datatype for HDF5 file based on C char
n = len_trim(buffer)
if (n > 0) then
call h5tcopy_f(H5T_C_S1, filetype, hdf5_err)
call h5tset_size_f(filetype, n, hdf5_err)
! Create memory space and attribute
call h5screate_f(H5S_SCALAR_F, dspace_id, hdf5_err)
call h5acreate_f(obj_id, trim(name), filetype, dspace_id, &
attr_id, hdf5_err)
! Copy string to temporary buffer
allocate(temp_buffer(n))
do i = 1, n
temp_buffer(i) = buffer(i:i)
end do
! Write attribute
f_ptr = c_loc(buffer(1:1))
call h5awrite_f(attr_id, filetype, f_ptr, hdf5_err)
! Close attribute
call h5aclose_f(attr_id, hdf5_err)
call h5sclose_f(dspace_id, hdf5_err)
call h5tclose_f(filetype, hdf5_err)
end if
end subroutine write_attribute_string
subroutine read_attribute_double(buffer, obj_id, name)
@ -2270,6 +2302,37 @@ contains
call h5aread_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err)
end subroutine read_attribute_integer_1D_explicit
subroutine write_attribute_integer_1D(obj_id, name, buffer)
integer(HID_T), intent(in) :: obj_id
character(*), intent(in) :: name
integer, target, intent(in) :: buffer(:)
integer(HSIZE_T) :: dims(1)
dims(:) = shape(buffer)
call write_attribute_integer_1D_explicit(obj_id, dims, name, buffer)
end subroutine write_attribute_integer_1D
subroutine write_attribute_integer_1D_explicit(obj_id, dims, name, buffer)
integer(HID_T), intent(in) :: obj_id
integer(HSIZE_T), intent(in) :: dims(1)
character(*), intent(in) :: name
integer, target, intent(in) :: buffer(dims(1))
integer :: hdf5_err
integer(HID_T) :: dspace_id
integer(HID_T) :: attr_id
type(C_PTR) :: f_ptr
call h5screate_simple_f(1, dims, dspace_id, hdf5_err)
call h5acreate_f(obj_id, trim(name), H5T_NATIVE_INTEGER, dspace_id, &
attr_id, hdf5_err)
f_ptr = c_loc(buffer)
call h5awrite_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err)
call h5aclose_f(attr_id, hdf5_err)
call h5sclose_f(dspace_id, hdf5_err)
end subroutine write_attribute_integer_1D_explicit
subroutine read_attribute_integer_2D(buffer, obj_id, name)
integer, target, allocatable, intent(inout) :: buffer(:,:)
integer(HID_T), intent(in) :: obj_id

71
src/initialize.F90
View file

@ -1,22 +1,30 @@
module initialize
use, intrinsic :: ISO_C_BINDING, only: c_loc
use hdf5
#ifdef _OPENMP
use omp_lib
#endif
use bank_header, only: Bank
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
use set_header, only: SetInt
use energy_grid, only: logarithmic_grid, grid_method
use error, only: fatal_error, warning
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
maximum_levels
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
&BASE_UNIVERSE
use global
use hdf5_interface, only: file_open, read_dataset, file_close, hdf5_bank_t,&
hdf5_integer8_t
use hdf5_interface, only: file_open, read_attribute, file_close, &
hdf5_bank_t, hdf5_integer8_t
use input_xml, only: read_input_xml, cells_in_univ_dict, read_plots_xml
use material_header, only: Material
use message_passing
use mgxs_data, only: read_mgxs, create_macro_xs
use output, only: title, header, print_version, write_message, &
use output, only: title, header, print_version, write_message, &
print_usage, print_plot
use random_lcg, only: initialize_prng
use state_point, only: load_state_point
@ -27,30 +35,23 @@ module initialize
use tally_filter
use tally, only: init_tally_routines
#ifdef MPI
use message_passing
#endif
#ifdef _OPENMP
use omp_lib
#endif
use hdf5
use, intrinsic :: ISO_C_BINDING, only: c_loc
implicit none
contains
!===============================================================================
! INITIALIZE_RUN takes care of all initialization tasks, i.e. reading
! OPENMC_INIT takes care of all initialization tasks, i.e. reading
! from command line, reading xml input files, initializing random
! number seeds, reading cross sections, initializing starting source,
! setting up timers, etc.
!===============================================================================
subroutine initialize_run()
subroutine openmc_init(intracomm)
#ifdef MPIF08
type(MPI_Comm), intent(in) :: intracomm ! MPI intracommunicator
#else
integer, intent(in), optional :: intracomm ! MPI intracommunicator
#endif
! Start total and initialization timer
call time_total%start()
@ -58,7 +59,7 @@ contains
#ifdef MPI
! Setup MPI
call initialize_mpi()
call initialize_mpi(intracomm)
#endif
! Initialize HDF5 interface
@ -155,7 +156,7 @@ contains
! Stop initialization timer
call time_initialize%stop()
end subroutine initialize_run
end subroutine openmc_init
#ifdef MPI
!===============================================================================
@ -164,7 +165,12 @@ contains
! each processor.
!===============================================================================
subroutine initialize_mpi()
subroutine initialize_mpi(intracomm)
#ifdef MPIF08
type(MPI_Comm), intent(in) :: intracomm ! MPI intracommunicator
#else
integer, intent(in) :: intracomm ! MPI intracommunicator
#endif
integer :: bank_blocks(5) ! Count for each datatype
#ifdef MPIF08
@ -182,8 +188,9 @@ contains
call MPI_INIT(mpi_err)
! Determine number of processors and rank of each processor
call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, mpi_err)
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_err)
mpi_intracomm = intracomm
call MPI_COMM_SIZE(mpi_intracomm, n_procs, mpi_err)
call MPI_COMM_RANK(mpi_intracomm, rank, mpi_err)
! Determine master
if (rank == 0) then
@ -300,11 +307,11 @@ contains
! Check what type of file this is
file_id = file_open(argv(i), 'r', parallel=.true.)
call read_dataset(filetype, file_id, 'filetype')
call read_attribute(filetype, file_id, 'filetype')
call file_close(file_id)
! Set path and flag for type of run
select case (filetype)
select case (trim(filetype))
case ('statepoint')
path_state_point = argv(i)
restart_run = .true.
@ -312,7 +319,7 @@ contains
path_particle_restart = argv(i)
particle_restart_run = .true.
case default
call fatal_error("Unrecognized file after restart flag.")
call fatal_error("Unrecognized file after restart flag: " // filetype // ".")
end select
! If its a restart run check for additional source file
@ -326,7 +333,7 @@ contains
! Check file type is a source file
file_id = file_open(argv(i), 'r', parallel=.true.)
call read_dataset(filetype, file_id, 'filetype')
call read_attribute(filetype, file_id, 'filetype')
call file_close(file_id)
if (filetype /= 'source') then
call fatal_error("Second file after restart flag must be a &
@ -553,11 +560,11 @@ contains
if (c % material(1) == NONE) then
id = c % fill
if (universe_dict % has_key(id)) then
c % type = CELL_FILL
c % type = FILL_UNIVERSE
c % fill = universe_dict % get_key(id)
elseif (lattice_dict % has_key(id)) then
lid = lattice_dict % get_key(id)
c % type = CELL_LATTICE
c % type = FILL_LATTICE
c % fill = lid
else
call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "&
@ -568,9 +575,9 @@ contains
do j = 1, size(c % material)
id = c % material(j)
if (id == MATERIAL_VOID) then
c % type = CELL_NORMAL
c % type = FILL_MATERIAL
else if (material_dict % has_key(id)) then
c % type = CELL_NORMAL
c % type = FILL_MATERIAL
c % material(j) = material_dict % get_key(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
@ -927,7 +934,7 @@ contains
! Allocate offset table for fill cells
do i = 1, n_cells
if (cells(i) % type /= CELL_NORMAL) then
if (cells(i) % type /= FILL_MATERIAL) then
allocate(cells(i) % offset(n_maps))
end if
end do

254
src/input_xml.F90
View file

@ -15,6 +15,7 @@ module input_xml
use hdf5_interface
use list_header, only: ListChar, ListInt, ListReal
use mesh_header, only: RegularMesh
use message_passing
use mgxs_data, only: create_macro_xs, read_mgxs
use multipole, only: multipole_read
use output, only: write_message
@ -78,7 +79,6 @@ contains
integer :: temp_int
integer :: temp_int_array3(3)
integer, allocatable :: temp_int_array(:)
integer(8) :: temp_long
real(8), allocatable :: temp_real(:)
integer :: n_tracks
logical :: file_exists
@ -99,7 +99,6 @@ contains
type(Node), pointer :: node_res_scat => null()
type(Node), pointer :: node_scatterer => null()
type(Node), pointer :: node_trigger => null()
type(Node), pointer :: node_keff_trigger => null()
type(Node), pointer :: node_vol => null()
type(Node), pointer :: node_tab_leg => null()
type(NodeList), pointer :: node_scat_list => null()
@ -223,119 +222,62 @@ contains
end if
end if
! Make sure that either eigenvalue or fixed source was specified
if (.not. check_for_node(doc, "eigenvalue") .and. &
.not. check_for_node(doc, "fixed_source")) then
call fatal_error("<eigenvalue> or <fixed_source> not specified.")
end if
! Check run mode if it hasn't been set from the command line
if (run_mode == NONE) then
if (check_for_node(doc, "run_mode")) then
call get_node_value(doc, "run_mode", temp_str)
select case (to_lower(temp_str))
case ("eigenvalue")
run_mode = MODE_EIGENVALUE
case ("fixed source")
run_mode = MODE_FIXEDSOURCE
case ("plot")
run_mode = MODE_PLOTTING
case ("particle restart")
run_mode = MODE_PARTICLE
case ("volume")
run_mode = MODE_VOLUME
end select
! Eigenvalue information
if (check_for_node(doc, "eigenvalue")) then
! Set run mode
if (run_mode == NONE) run_mode = MODE_EIGENVALUE
! Assume XML specifics <particles>, <batches>, etc. directly
node_mode => doc
else
call warning("<run_mode> should be specified.")
! Get pointer to eigenvalue XML block
call get_node_ptr(doc, "eigenvalue", node_mode)
! Check number of particles
if (.not. check_for_node(node_mode, "particles")) then
call fatal_error("Need to specify number of particles per generation.")
end if
! Get number of particles
call get_node_value(node_mode, "particles", temp_long)
! If the number of particles was specified as a command-line argument, we
! don't set it here
if (n_particles == 0) n_particles = temp_long
! Get number of basic batches
call get_node_value(node_mode, "batches", n_batches)
if (.not. trigger_on) then
n_max_batches = n_batches
end if
! Get number of inactive batches
call get_node_value(node_mode, "inactive", n_inactive)
n_active = n_batches - n_inactive
if (check_for_node(node_mode, "generations_per_batch")) then
call get_node_value(node_mode, "generations_per_batch", gen_per_batch)
end if
! Allocate array for batch keff and entropy
allocate(k_generation(n_max_batches*gen_per_batch))
allocate(entropy(n_max_batches*gen_per_batch))
entropy = ZERO
! Get the trigger information for keff
if (check_for_node(node_mode, "keff_trigger")) then
call get_node_ptr(node_mode, "keff_trigger", node_keff_trigger)
if (check_for_node(node_keff_trigger, "type")) then
call get_node_value(node_keff_trigger, "type", temp_str)
temp_str = trim(to_lower(temp_str))
select case (temp_str)
case ('std_dev')
keff_trigger % trigger_type = STANDARD_DEVIATION
case ('variance')
keff_trigger % trigger_type = VARIANCE
case ('rel_err')
keff_trigger % trigger_type = RELATIVE_ERROR
case default
call fatal_error("Unrecognized keff trigger type " // temp_str)
end select
else
call fatal_error("Specify keff trigger type in settings XML")
! Make sure that either eigenvalue or fixed source was specified
if (.not. check_for_node(doc, "eigenvalue") .and. &
.not. check_for_node(doc, "fixed_source")) then
call fatal_error("<eigenvalue> or <fixed_source> not specified.")
end if
if (check_for_node(node_keff_trigger, "threshold")) then
call get_node_value(node_keff_trigger, "threshold", &
keff_trigger % threshold)
else
call fatal_error("Specify keff trigger threshold in settings XML")
if (check_for_node(doc, "eigenvalue")) then
! Set run mode
if (run_mode == NONE) run_mode = MODE_EIGENVALUE
! Get pointer to eigenvalue XML block
call get_node_ptr(doc, "eigenvalue", node_mode)
elseif (check_for_node(doc, "fixed_source")) then
! Set run mode
if (run_mode == NONE) run_mode = MODE_FIXEDSOURCE
! Get pointer to fixed_source XML block
call get_node_ptr(doc, "fixed_source", node_mode)
end if
end if
end if
! Fixed source calculation information
if (check_for_node(doc, "fixed_source")) then
! Set run mode
if (run_mode == NONE) run_mode = MODE_FIXEDSOURCE
if (run_mode == MODE_EIGENVALUE .or. run_mode == MODE_FIXEDSOURCE) then
! Read run parameters
call get_run_parameters(node_mode)
! Get pointer to fixed_source XML block
call get_node_ptr(doc, "fixed_source", node_mode)
! Check number of particles
if (.not. check_for_node(node_mode, "particles")) then
call fatal_error("Need to specify number of particles per batch.")
! Check number of active batches, inactive batches, and particles
if (n_active <= 0) then
call fatal_error("Number of active batches must be greater than zero.")
elseif (n_inactive < 0) then
call fatal_error("Number of inactive batches must be non-negative.")
elseif (n_particles <= 0) then
call fatal_error("Number of particles must be greater than zero.")
end if
! Get number of particles
call get_node_value(node_mode, "particles", temp_long)
! If the number of particles was specified as a command-line argument, we
! don't set it here
if (n_particles == 0) n_particles = temp_long
! Copy batch information
call get_node_value(node_mode, "batches", n_batches)
if (.not. trigger_on) then
n_max_batches = n_batches
end if
n_active = n_batches
n_inactive = 0
gen_per_batch = 1
end if
! Check number of active batches, inactive batches, and particles
if (n_active <= 0) then
call fatal_error("Number of active batches must be greater than zero.")
elseif (n_inactive < 0) then
call fatal_error("Number of inactive batches must be non-negative.")
elseif (n_particles <= 0) then
call fatal_error("Number of particles must be greater than zero.")
end if
! Copy random number seed if specified
@ -379,7 +321,10 @@ contains
! Get point to list of <source> elements and make sure there is at least one
call get_node_list(doc, "source", node_source_list)
n = get_list_size(node_source_list)
if (n == 0) call fatal_error("No source specified in settings XML file.")
if (run_mode == MODE_EIGENVALUE .or. run_mode == MODE_FIXEDSOURCE) then
if (n == 0) call fatal_error("No source specified in settings XML file.")
end if
! Allocate array for sources
allocate(external_source(n))
@ -1092,6 +1037,86 @@ contains
end subroutine read_settings_xml
!===============================================================================
! GET_RUN_PARAMETERS
!===============================================================================
subroutine get_run_parameters(node_base)
type(Node), pointer :: node_base
integer(8) :: temp_long
character(MAX_LINE_LEN) :: temp_str
type(Node), pointer :: node_keff_trigger => null()
! Check number of particles
if (.not. check_for_node(node_base, "particles")) then
call fatal_error("Need to specify number of particles.")
end if
! Get number of particles
call get_node_value(node_base, "particles", temp_long)
! If the number of particles was specified as a command-line argument, we
! don't set it here
if (n_particles == 0) n_particles = temp_long
! Get number of basic batches
call get_node_value(node_base, "batches", n_batches)
if (.not. trigger_on) then
n_max_batches = n_batches
end if
n_inactive = 0
gen_per_batch = 1
! Get number of inactive batches
if (run_mode == MODE_EIGENVALUE) then
call get_node_value(node_base, "inactive", n_inactive)
if (check_for_node(node_base, "generations_per_batch")) then
call get_node_value(node_base, "generations_per_batch", gen_per_batch)
end if
! Allocate array for batch keff and entropy
allocate(k_generation(n_max_batches*gen_per_batch))
allocate(entropy(n_max_batches*gen_per_batch))
entropy = ZERO
! Get the trigger information for keff
if (check_for_node(node_base, "keff_trigger")) then
call get_node_ptr(node_base, "keff_trigger", node_keff_trigger)
if (check_for_node(node_keff_trigger, "type")) then
call get_node_value(node_keff_trigger, "type", temp_str)
temp_str = trim(to_lower(temp_str))
select case (temp_str)
case ('std_dev')
keff_trigger % trigger_type = STANDARD_DEVIATION
case ('variance')
keff_trigger % trigger_type = VARIANCE
case ('rel_err')
keff_trigger % trigger_type = RELATIVE_ERROR
case default
call fatal_error("Unrecognized keff trigger type " // temp_str)
end select
else
call fatal_error("Specify keff trigger type in settings XML")
end if
if (check_for_node(node_keff_trigger, "threshold")) then
call get_node_value(node_keff_trigger, "threshold", &
keff_trigger % threshold)
else
call fatal_error("Specify keff trigger threshold in settings XML")
end if
end if
end if
! Determine number of active batches
n_active = n_batches - n_inactive
end subroutine get_run_parameters
!===============================================================================
! READ_GEOMETRY_XML reads data from a geometry.xml file and parses it, checking
! for errors and placing properly-formatted data in the right data structures
@ -2226,6 +2251,13 @@ contains
call fatal_error("Must specify id of material in materials XML file")
end if
! Check if material is depletable
if (check_for_node(node_mat, "depletable")) then
call get_node_value(node_mat, "depletable", temp_str)
if (to_lower(temp_str) == "true" .or. temp_str == "1") &
mat % depletable = .true.
end if
! Check to make sure 'id' hasn't been used
if (material_dict % has_key(mat % id)) then
call fatal_error("Two or more materials use the same unique ID: " &
@ -4229,14 +4261,14 @@ contains
end select
! Set output file path
filename = trim(to_str(pl % id)) // "_plot"
filename = "plot_" // trim(to_str(pl % id))
if (check_for_node(node_plot, "filename")) &
call get_node_value(node_plot, "filename", filename)
select case (pl % type)
case (PLOT_TYPE_SLICE)
pl % path_plot = trim(path_input) // trim(filename) // ".ppm"
case (PLOT_TYPE_VOXEL)
pl % path_plot = trim(path_input) // trim(filename) // ".voxel"
pl % path_plot = trim(path_input) // trim(filename) // ".h5"
end select
! Copy plot pixel size
@ -5380,16 +5412,16 @@ contains
if (attribute_exists(file_id, 'version')) then
call read_attribute(version, file_id, 'version')
if (version(1) /= HDF5_VERSION_MAJOR) then
if (version(1) /= HDF5_VERSION(1)) then
call fatal_error("HDF5 data format uses version " // trim(to_str(&
version(1))) // "." // trim(to_str(version(2))) // " whereas &
&your installation of OpenMC expects version " // trim(to_str(&
HDF5_VERSION_MAJOR)) // ".x data.")
HDF5_VERSION(1))) // ".x data.")
end if
else
call fatal_error("HDF5 data does not indicate a version. Your &
&installation of OpenMC expects version " // trim(to_str(&
HDF5_VERSION_MAJOR)) // ".x data.")
HDF5_VERSION(1))) // ".x data.")
end if
end subroutine check_data_version

18
src/main.F90
View file

@ -1,17 +1,23 @@
program main
use constants
use finalize, only: finalize_run
use finalize, only: openmc_finalize
use global
use initialize, only: initialize_run
use initialize, only: openmc_init
use message_passing
use particle_restart, only: run_particle_restart
use plot, only: run_plot
use simulation, only: run_simulation
use volume_calc, only: run_volume_calculations
implicit none
! set up problem
call initialize_run()
! Initialize run -- when run with MPI, pass communicator
#ifdef MPI
call openmc_init(MPI_COMM_WORLD)
#else
call openmc_init()
#endif
! start problem based on mode
select case (run_mode)
@ -21,9 +27,11 @@ program main
call run_plot()
case (MODE_PARTICLE)
if (master) call run_particle_restart()
case (MODE_VOLUME)
call run_volume_calculations()
end select
! finalize run
call finalize_run()
call openmc_finalize()
end program main

3
src/material_header.F90
View file

@ -31,8 +31,9 @@ module material_header
character(20), allocatable :: names(:) ! isotope names
character(20), allocatable :: sab_names(:) ! name of S(a,b) table
! Does this material contain fissionable nuclides?
! Does this material contain fissionable nuclides? Is it depletable?
logical :: fissionable = .false.
logical :: depletable = .false.
! enforce isotropic scattering in lab
logical, allocatable :: p0(:)

11
src/mesh.F90
View file

@ -1,13 +1,10 @@
module mesh
#ifdef MPI
use message_passing
#endif
use algorithm, only: binary_search
use constants
use global
use mesh_header
use message_passing
implicit none
@ -207,17 +204,17 @@ contains
! collect values from all processors
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, cnt, n, MPI_REAL8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
mpi_intracomm, mpi_err)
else
! Receive buffer not significant at other processors
call MPI_REDUCE(cnt, dummy, n, MPI_REAL8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
mpi_intracomm, mpi_err)
end if
! Check if there were sites outside the mesh for any processor
if (present(sites_outside)) then
call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, &
MPI_COMM_WORLD, mpi_err)
mpi_intracomm, mpi_err)
end if
#else

17
src/message_passing.F90
View file

@ -8,4 +8,21 @@ module message_passing
#endif
#endif
! The defaults set here for the number of processors, rank, and master and
! mpi_enabled flag are for when MPI is not being used at all, i.e. a serial
! run. In this case, these variables are still used at times.
integer :: n_procs = 1 ! number of processes
integer :: rank = 0 ! rank of process
logical :: master = .true. ! master process?
logical :: mpi_enabled = .false. ! is MPI in use and initialized?
integer :: mpi_err ! MPI error code
#ifdef MPIF08
type(MPI_Datatype) :: MPI_BANK ! MPI datatype for fission bank
type(MPI_Comm) :: mpi_intracomm ! MPI intra-communicator
#else
integer :: MPI_BANK ! MPI datatype for fission bank
integer :: mpi_intracomm ! MPI intra-communicator
#endif
end module message_passing

8
src/mgxs_header.F90
View file

@ -2340,7 +2340,7 @@ module mgxs_header
! Now need to compare this material maximum scattering order with
! the problem wide max scatt order and use whichever is lower
order = min(mat_max_order, max_order)
order = min(mat_max_order, max_order + 1)
! Ok, got our order, store the dimensionality
order_dim = order
@ -3467,9 +3467,7 @@ module mgxs_header
type(MaterialMacroXS), intent(inout) :: xs ! Resultant Mgxs Data
xs % total = this % xs(this % index_temp) % total(gin)
xs % elastic = this % xs(this % index_temp) % scatter % scattxs(gin)
xs % absorption = this % xs(this % index_temp) % absorption(gin)
xs % fission = this % xs(this % index_temp) % fission(gin)
xs % nu_fission = &
this % xs(this % index_temp) % prompt_nu_fission(gin) + &
sum(this % xs(this % index_temp) % delayed_nu_fission(:, gin))
@ -3487,12 +3485,8 @@ module mgxs_header
call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol)
xs % total = this % xs(this % index_temp) % &
total(gin, iazi, ipol)
xs % elastic = this % xs(this % index_temp) % &
scatter(iazi, ipol) % obj % scattxs(gin)
xs % absorption = this % xs(this % index_temp) % &
absorption(gin, iazi, ipol)
xs % fission = this % xs(this % index_temp) % &
fission(gin, iazi, ipol)
xs % nu_fission = this % xs(this % index_temp) % &
prompt_nu_fission(gin, iazi, ipol) + &
sum(this % xs(this % index_temp) % &

9
src/multipole.F90
View file

@ -1,8 +1,9 @@
module multipole
use hdf5
use constants
use global
use hdf5
use hdf5_interface
use multipole_header, only: MultipoleArray, FIT_T, FIT_A, FIT_F, &
MP_FISS, FORM_MLBW, FORM_RM
@ -39,9 +40,9 @@ contains
! Check the file version number.
call read_dataset(version, file_id, "version")
if (version /= MULTIPOLE_VERSION) call fatal_error("The current multipole&
& format version is " // trim(MULTIPOLE_VERSION) // " but the file "&
// trim(filename) // " uses version " // trim(version))
if (version /= VERSION_MULTIPOLE) call fatal_error("The current multipole&
& format version is " // trim(VERSION_MULTIPOLE) // " but the file "&
// trim(filename) // " uses version " // trim(version) // ".")
! Load in all the array size scalars
call read_dataset(multipole % length, group_id, "length")

6
src/output.F90
View file

@ -11,6 +11,7 @@ module output
use math, only: t_percentile
use mesh_header, only: RegularMesh
use mesh, only: mesh_indices_to_bin, bin_to_mesh_indices
use message_passing, only: master, n_procs
use nuclide_header
use particle_header, only: LocalCoord, Particle
use plot_header
@ -172,10 +173,13 @@ contains
if (master) then
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I1,".",I1)') &
"OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,A)') "Git SHA1: ", GIT_SHA1
#endif
write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2015 &
&Massachusetts Institute of Technology"
write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at &
&<http://mit-crpg.github.io/openmc/license.html>"
&<http://openmc.readthedocs.io/en/latest/license.html>"
end if
end subroutine print_version

9
src/particle_restart.F90
View file

@ -69,7 +69,6 @@ contains
type(Particle), intent(inout) :: p
integer, intent(inout) :: previous_run_mode
integer :: int_scalar
integer(HID_T) :: file_id
character(MAX_WORD_LEN) :: tempstr
@ -81,15 +80,13 @@ contains
file_id = file_open(path_particle_restart, 'r')
! Read data from file
call read_dataset(tempstr, file_id, 'filetype')
call read_dataset(int_scalar, file_id, 'revision')
call read_dataset(current_batch, file_id, 'current_batch')
call read_dataset(gen_per_batch, file_id, 'gen_per_batch')
call read_dataset(current_gen, file_id, 'current_gen')
call read_dataset(gen_per_batch, file_id, 'generations_per_batch')
call read_dataset(current_gen, file_id, 'current_generation')
call read_dataset(n_particles, file_id, 'n_particles')
call read_dataset(tempstr, file_id, 'run_mode')
select case (tempstr)
case ('k-eigenvalue')
case ('eigenvalue')
previous_run_mode = MODE_EIGENVALUE
case ('fixed source')
previous_run_mode = MODE_FIXEDSOURCE

50
src/particle_restart_write.F90
View file

@ -23,7 +23,6 @@ contains
integer(HID_T) :: file_id
character(MAX_FILE_LEN) :: filename
type(Bank), pointer :: src
! Dont write another restart file if in particle restart mode
if (run_mode == MODE_PARTICLE) return
@ -36,29 +35,34 @@ contains
! Create file
file_id = file_create(filename)
! Get information about source particle
src => source_bank(current_work)
associate (src => source_bank(current_work))
! Write filetype and version info
call write_attribute(file_id, 'filetype', 'particle restart')
call write_attribute(file_id, 'version', VERSION_PARTICLE_RESTART)
call write_attribute(file_id, "openmc_version", VERSION)
#ifdef GIT_SHA1
call write_attribute(file_id, "git_sha1", GIT_SHA1)
#endif
! Write data to file
call write_dataset(file_id, 'filetype', 'particle restart')
call write_dataset(file_id, 'revision', REVISION_PARTICLE_RESTART)
call write_dataset(file_id, 'current_batch', current_batch)
call write_dataset(file_id, 'gen_per_batch', gen_per_batch)
call write_dataset(file_id, 'current_gen', current_gen)
call write_dataset(file_id, 'n_particles', n_particles)
select case(run_mode)
case (MODE_FIXEDSOURCE)
call write_dataset(file_id, 'run_mode', 'fixed source')
case (MODE_EIGENVALUE)
call write_dataset(file_id, 'run_mode', 'k-eigenvalue')
case (MODE_PARTICLE)
call write_dataset(file_id, 'run_mode', 'particle restart')
end select
call write_dataset(file_id, 'id', p%id)
call write_dataset(file_id, 'weight', src%wgt)
call write_dataset(file_id, 'energy', src%E)
call write_dataset(file_id, 'xyz', src%xyz)
call write_dataset(file_id, 'uvw', src%uvw)
! Write data to file
call write_dataset(file_id, 'current_batch', current_batch)
call write_dataset(file_id, 'generations_per_batch', gen_per_batch)
call write_dataset(file_id, 'current_generation', current_gen)
call write_dataset(file_id, 'n_particles', n_particles)
select case(run_mode)
case (MODE_FIXEDSOURCE)
call write_dataset(file_id, 'run_mode', 'fixed source')
case (MODE_EIGENVALUE)
call write_dataset(file_id, 'run_mode', 'eigenvalue')
case (MODE_PARTICLE)
call write_dataset(file_id, 'run_mode', 'particle restart')
end select
call write_dataset(file_id, 'id', p%id)
call write_dataset(file_id, 'weight', src%wgt)
call write_dataset(file_id, 'energy', src%E)
call write_dataset(file_id, 'xyz', src%xyz)
call write_dataset(file_id, 'uvw', src%uvw)
end associate
! Close file
call file_close(file_id)

1
src/physics.F90
View file

@ -9,6 +9,7 @@ module physics
use material_header, only: Material
use math
use mesh, only: get_mesh_indices
use message_passing
use nuclide_header
use output, only: write_message
use particle_header, only: Particle

1
src/physics_mg.F90
View file

@ -9,6 +9,7 @@ module physics_mg
use math, only: rotate_angle
use mgxs_header, only: Mgxs, MgxsContainer
use mesh, only: get_mesh_indices
use message_passing
use output, only: write_message
use particle_header, only: Particle
use particle_restart_write, only: write_particle_restart

270
src/plot.F90
View file

@ -8,17 +8,22 @@ module plot
use hdf5_interface
use mesh, only: get_mesh_indices
use mesh_header, only: RegularMesh
use output, only: write_message
use output, only: write_message, time_stamp
use particle_header, only: LocalCoord, Particle
use plot_header
use ppmlib, only: Image, init_image, allocate_image, &
deallocate_image, set_pixel
use progress_header, only: ProgressBar
use string, only: to_str
use hdf5
implicit none
private
public :: run_plot
integer, parameter :: RED = 1
integer, parameter :: GREEN = 2
integer, parameter :: BLUE = 3
contains
@ -29,22 +34,21 @@ contains
subroutine run_plot()
integer :: i ! loop index for plots
type(ObjectPlot), pointer :: pl => null()
do i = 1, n_plots
pl => plots(i)
associate (pl => plots(i))
! Display output message
call write_message("Processing plot " // trim(to_str(pl % id)) &
// ": " // trim(pl % path_plot) // " ...", 5)
! Display output message
call write_message("Processing plot " // trim(to_str(pl % id)) &
&// ": " // trim(pl % path_plot) // " ...", 5)
if (pl % type == PLOT_TYPE_SLICE) then
! create 2d image
call create_ppm(pl)
else if (pl % type == PLOT_TYPE_VOXEL) then
! create dump for 3D silomesh utility script
call create_3d_dump(pl)
end if
if (pl % type == PLOT_TYPE_SLICE) then
! create 2d image
call create_ppm(pl)
else if (pl % type == PLOT_TYPE_VOXEL) then
! create dump for 3D silomesh utility script
call create_voxel(pl)
end if
end associate
end do
end subroutine run_plot
@ -55,15 +59,13 @@ contains
!===============================================================================
subroutine position_rgb(p, pl, rgb, id)
type(Particle), intent(inout) :: p
type(ObjectPlot), pointer, intent(in) :: pl
integer, intent(out) :: rgb(3)
integer, intent(out) :: id
type(Particle), intent(inout) :: p
type(ObjectPlot), intent(in) :: pl
integer, intent(out) :: rgb(3)
integer, intent(out) :: id
integer :: j
logical :: found_cell
type(Cell), pointer :: c
p % n_coord = 1
@ -81,19 +83,20 @@ contains
else
if (pl % color_by == PLOT_COLOR_MATS) then
! Assign color based on material
c => cells(p % coord(j) % cell)
if (c % type == CELL_FILL) then
! If we stopped on a middle universe level, treat as if not found
rgb = pl % not_found % rgb
id = -1
else if (p % material == MATERIAL_VOID) then
! By default, color void cells white
rgb = 255
id = -1
else
rgb = pl % colors(p % material) % rgb
id = materials(p % material) % id
end if
associate (c => cells(p % coord(j) % cell))
if (c % type == FILL_UNIVERSE) then
! If we stopped on a middle universe level, treat as if not found
rgb = pl % not_found % rgb
id = -1
else if (p % material == MATERIAL_VOID) then
! By default, color void cells white
rgb = 255
id = -1
else
rgb = pl % colors(p % material) % rgb
id = materials(p % material) % id
end if
end associate
else if (pl % color_by == PLOT_COLOR_CELLS) then
! Assign color based on cell
rgb = pl % colors(p % coord(j) % cell) % rgb
@ -112,27 +115,29 @@ contains
!===============================================================================
subroutine create_ppm(pl)
type(ObjectPlot), pointer :: pl
type(ObjectPlot), intent(in) :: pl
integer :: in_i
integer :: out_i
integer :: x, y ! pixel location
integer :: rgb(3) ! colors (red, green, blue) from 0-255
integer :: id
integer :: height, width
real(8) :: in_pixel
real(8) :: out_pixel
real(8) :: xyz(3)
type(Image) :: img
integer, allocatable :: data(:,:,:)
type(Particle) :: p
type(ProgressBar) :: progress
! Initialize and allocate space for image
call init_image(img)
call allocate_image(img, pl % pixels(1), pl % pixels(2))
width = pl % pixels(1)
height = pl % pixels(2)
in_pixel = pl % width(1)/dble(pl % pixels(1))
out_pixel = pl % width(2)/dble(pl % pixels(2))
in_pixel = pl % width(1)/dble(width)
out_pixel = pl % width(2)/dble(height)
! Allocate and initialize results array
allocate(data(3, width, height))
data(:,:,:) = 0
if (pl % basis == PLOT_BASIS_XY) then
in_i = 1
@ -160,50 +165,42 @@ contains
p % coord(1) % uvw = [ HALF, HALF, HALF ]
p % coord(1) % universe = BASE_UNIVERSE
do y = 1, img % height
call progress % set_value(dble(y)/dble(img % height)*100)
do x = 1, img % width
!$omp parallel do firstprivate(p) private(x, rgb, id) reduction(+ : data)
do y = 1, height
! Set y coordinate
p % coord(1) % xyz(out_i) = xyz(out_i) - out_pixel*(y - 1)
do x = 1, width
! Set x coordinate
p % coord(1) % xyz(in_i) = xyz(in_i) + in_pixel*(x - 1)
! get pixel color
call position_rgb(p, pl, rgb, id)
! Create a pixel at (x,y) with color (r,g,b)
call set_pixel(img, x-1, y-1, rgb(1), rgb(2), rgb(3))
! Advance pixel in first direction
p % coord(1) % xyz(in_i) = p % coord(1) % xyz(in_i) + in_pixel
data(:, x, y) = rgb
end do
! Advance pixel in second direction
p % coord(1) % xyz(in_i) = xyz(in_i)
p % coord(1) % xyz(out_i) = p % coord(1) % xyz(out_i) - out_pixel
end do
!$omp end parallel do
! Draw tally mesh boundaries on the image if requested
if (associated(pl % meshlines_mesh)) call draw_mesh_lines(pl, img)
if (associated(pl % meshlines_mesh)) call draw_mesh_lines(pl, data)
! Write out the ppm to a file
call output_ppm(pl,img)
! Free up space
call deallocate_image(img)
! Clear particle
call p % clear()
call output_ppm(pl, data)
end subroutine create_ppm
!===============================================================================
! DRAW_MESH_LINES draws mesh line boundaries on an image
!===============================================================================
subroutine draw_mesh_lines(pl, img)
type(ObjectPlot), pointer, intent(in) :: pl
type(Image), intent(inout) :: img
subroutine draw_mesh_lines(pl, data)
type(ObjectPlot), intent(in) :: pl
integer, intent(inout) :: data(:,:,:)
logical :: in_mesh
integer :: out_, in_ ! pixel location
integer :: r, g, b ! RGB color for meshlines pixels
integer :: rgb(3) ! RGB color for meshlines pixels
integer :: outrange(2), inrange(2) ! range of pixel locations
integer :: i, j ! loop indices
integer :: plus
@ -216,13 +213,8 @@ contains
real(8) :: xyz_ur_plot(3) ! upper right xyz of plot image
real(8) :: xyz_ll(3) ! lower left xyz
real(8) :: xyz_ur(3) ! upper right xyz
type(RegularMesh), pointer :: m
m => pl % meshlines_mesh
r = pl % meshlines_color % rgb(1)
g = pl % meshlines_color % rgb(2)
b = pl % meshlines_color % rgb(3)
rgb(:) = pl % meshlines_color % rgb
select case (pl % basis)
case(PLOT_BASIS_XY)
@ -246,55 +238,57 @@ contains
width = xyz_ur_plot - xyz_ll_plot
call get_mesh_indices(m, xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh)
call get_mesh_indices(m, xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh)
associate (m => pl % meshlines_mesh)
call get_mesh_indices(m, xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh)
call get_mesh_indices(m, xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh)
! sweep through all meshbins on this plane and draw borders
do i = ijk_ll(outer), ijk_ur(outer)
do j = ijk_ll(inner), ijk_ur(inner)
! check if we're in the mesh for this ijk
if (i > 0 .and. i <= m % dimension(outer) .and. &
j > 0 .and. j <= m % dimension(inner)) then
! sweep through all meshbins on this plane and draw borders
do i = ijk_ll(outer), ijk_ur(outer)
do j = ijk_ll(inner), ijk_ur(inner)
! check if we're in the mesh for this ijk
if (i > 0 .and. i <= m % dimension(outer) .and. &
j > 0 .and. j <= m % dimension(inner)) then
! get xyz's of lower left and upper right of this mesh cell
xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1)
xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1)
xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i
xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j
! get xyz's of lower left and upper right of this mesh cell
xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1)
xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1)
xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i
xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j
! map the xyz ranges to pixel ranges
! map the xyz ranges to pixel ranges
frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer)
outrange(1) = int(frac * real(img % width, 8))
frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer)
outrange(2) = int(frac * real(img % width, 8))
frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer)
outrange(1) = int(frac * real(pl % pixels(1), 8))
frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer)
outrange(2) = int(frac * real(pl % pixels(1), 8))
frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner)
inrange(1) = int((ONE - frac) * real(img % height, 8))
frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner)
inrange(2) = int((ONE - frac) * real(img % height, 8))
frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner)
inrange(1) = int((ONE - frac) * real(pl % pixels(2), 8))
frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner)
inrange(2) = int((ONE - frac) * real(pl % pixels(2), 8))
! draw lines
do out_ = outrange(1), outrange(2)
do plus = 0, pl % meshlines_width
call set_pixel(img, out_, inrange(1) + plus, r, g, b)
call set_pixel(img, out_, inrange(2) + plus, r, g, b)
call set_pixel(img, out_, inrange(1) - plus, r, g, b)
call set_pixel(img, out_, inrange(2) - plus, r, g, b)
! draw lines
do out_ = outrange(1), outrange(2)
do plus = 0, pl % meshlines_width
data(:, out_ + 1, inrange(1) + plus + 1) = rgb
data(:, out_ + 1, inrange(2) + plus + 1) = rgb
data(:, out_ + 1, inrange(1) - plus + 1) = rgb
data(:, out_ + 1, inrange(2) - plus + 1) = rgb
end do
end do
end do
do in_ = inrange(1), inrange(2)
do plus = 0, pl % meshlines_width
call set_pixel(img, outrange(1) + plus, in_, r, g, b)
call set_pixel(img, outrange(2) + plus, in_, r, g, b)
call set_pixel(img, outrange(1) - plus, in_, r, g, b)
call set_pixel(img, outrange(2) - plus, in_, r, g, b)
do in_ = inrange(1), inrange(2)
do plus = 0, pl % meshlines_width
data(:, outrange(1) + plus + 1, in_ + 1) = rgb
data(:, outrange(2) + plus + 1, in_ + 1) = rgb
data(:, outrange(1) - plus + 1, in_ + 1) = rgb
data(:, outrange(2) - plus + 1, in_ + 1) = rgb
end do
end do
end do
end if
end if
end do
end do
end do
end associate
end subroutine draw_mesh_lines
@ -302,13 +296,12 @@ contains
! OUTPUT_PPM writes out a previously generated image to a PPM file
!===============================================================================
subroutine output_ppm(pl, img)
subroutine output_ppm(pl, data)
type(ObjectPlot), intent(in) :: pl
integer, intent(in) :: data(:,:,:)
type(ObjectPlot), pointer :: pl
type(Image), intent(in) :: img
integer :: i ! loop index for height
integer :: j ! loop index for width
integer :: y ! loop index for height
integer :: x ! loop index for width
integer :: unit_plot
! Open PPM file for writing
@ -316,24 +309,23 @@ contains
! Write header
write(unit_plot, '(A2)') 'P6'
write(unit_plot, '(I0,'' '',I0)') img%width, img%height
write(unit_plot, '(I0,'' '',I0)') pl % pixels(1), pl % pixels(2)
write(unit_plot, '(A)') '255'
! Write color for each pixel
do j = 1, img % height
do i = 1, img % width
write(unit_plot, '(3A1)', advance='no') achar(img%red(i,j)), &
achar(img%green(i,j)), achar(img%blue(i,j))
do y = 1, pl % pixels(2)
do x = 1, pl % pixels(1)
write(unit_plot, '(3A1)', advance='no') achar(data(RED, x, y)), &
achar(data(GREEN, x, y)), achar(data(BLUE, x, y))
end do
end do
! Close plot file
close(UNIT=unit_plot)
end subroutine output_ppm
!===============================================================================
! CREATE_3D_DUMP outputs a binary file that can be input into silomesh for 3D
! CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D
! geometry visualization. It works the same way as create_ppm by dragging a
! particle across the geometry for the specified number of voxels. The first
! 3 int(4)'s in the binary are the number of x, y, and z voxels. The next 3
@ -344,9 +336,8 @@ contains
! id. For 1 million voxels this produces a file of approximately 15MB.
!===============================================================================
subroutine create_3d_dump(pl)
type(ObjectPlot), pointer :: pl
subroutine create_voxel(pl)
type(ObjectPlot), intent(in) :: pl
integer :: x, y, z ! voxel location indices
integer :: rgb(3) ! colors (red, green, blue) from 0-255
@ -355,7 +346,7 @@ contains
integer, target :: data(pl%pixels(3),pl%pixels(2))
integer(HID_T) :: file_id
integer(HID_T) :: dspace
integeR(HID_T) :: memspace
integer(HID_T) :: memspace
integer(HID_T) :: dset
integer(HSIZE_T) :: dims(3)
integer(HSIZE_T) :: dims_slab(3)
@ -381,11 +372,20 @@ contains
! Open binary plot file for writing
file_id = file_create(pl%path_plot)
! write plot header info
call write_dataset(file_id, "filetype", 'voxel')
call write_dataset(file_id, "num_voxels", pl%pixels)
call write_dataset(file_id, "voxel_width", vox)
call write_dataset(file_id, "lower_left", ll)
! write header info
call write_attribute(file_id, "filetype", 'voxel')
call write_attribute(file_id, "version", VERSION_VOXEL)
call write_attribute(file_id, "openmc_version", VERSION)
#ifdef GIT_SHA1
call write_attribute(file_id, "git_sha1", GIT_SHA1)
#endif
! Write current date and time
call write_attribute(file_id, "date_and_time", time_stamp())
call write_attribute(file_id, "num_voxels", pl%pixels)
call write_attribute(file_id, "voxel_width", vox)
call write_attribute(file_id, "lower_left", ll)
! Create dataset for voxel data -- note that the dimensions are reversed
! since we want the order in the file to be z, y, x
@ -443,6 +443,6 @@ contains
call h5sclose_f(memspace, hdf5_err)
call file_close(file_id)
end subroutine create_3d_dump
end subroutine create_voxel
end module plot

127
src/ppmlib.F90
View file

@ -1,127 +0,0 @@
module ppmlib
implicit none
!===============================================================================
! Image holds RGB information for output PPM image
!===============================================================================
type Image
integer, dimension(:,:), pointer :: red, green, blue
integer :: width, height
end type Image
contains
!===============================================================================
! INIT_IMAGE initializes the Image derived type
!===============================================================================
subroutine init_image(img)
type(Image), intent(out) :: img
nullify(img % red)
nullify(img % green)
nullify(img % blue)
img % width = 0
img % height = 0
end subroutine init_image
!===============================================================================
! ALLOCATE_IMAGE sets the width and height of an image and allocates color
! arrays
!===============================================================================
subroutine allocate_image(img, w, h)
type(Image), intent(inout) :: img
integer, intent(in) :: w ! width of image
integer, intent(in) :: h ! height of image
! allocate red, green, and blue array
allocate(img % red(w, h))
allocate(img % green(w, h))
allocate(img % blue(w, h))
! set width and height
img % width = w
img % height = h
end subroutine allocate_image
!===============================================================================
! DEALLOCATE_IMAGE
!===============================================================================
subroutine deallocate_image(img)
type(Image) :: img
if (associated(img % red)) deallocate(img % red)
if (associated(img % green)) deallocate(img % green)
if (associated(img % blue)) deallocate(img % blue)
end subroutine deallocate_image
!===============================================================================
! INSIDE_IMAGE determines whether a point (x,y) is inside the image
!===============================================================================
function inside_image(img, x, y) result(inside)
type(Image), intent(in) :: img
integer, intent(in) :: x, y
logical :: inside
inside = .false.
if ((x < img % width) .and. (y < img % height) .and. &
(x >= 0) .and. (y >= 0)) inside = .true.
end function inside_image
!===============================================================================
! VALID_IMAGE checks whether the image has a width and height and if its color
! arrays are allocated
!===============================================================================
function valid_image(img) result(valid)
type(Image), intent(in) :: img
logical :: valid
valid = .false.
if (img % width == 0) return
if (img % height == 0) return
if (.not. associated(img % red) .or. &
.not. associated(img % green) .or. &
.not. associated(img % blue)) return
valid = .true.
end function valid_image
!===============================================================================
! SET_PIXEL sets the colors for a given pixel
!===============================================================================
subroutine set_pixel(img, x, y, r, g, b)
type(Image), intent(inout) :: img
integer, intent(in) :: x, y ! coordinates
integer, intent(in) :: r, g, b ! red, green, and blue
if (inside_image(img, x, y) .and. valid_image(img)) then
img % red(x+1,y+1) = mod(abs(r), 256)
img % green(x+1, y+1) = mod(abs(g), 256)
img % blue(x+1, y+1) = mod(abs(b), 256)
end if
end subroutine set_pixel
end module ppmlib

4
src/ppmlib.LICENSE
View file

@ -1,4 +0,0 @@
The source code in ppmlib.F90 is adapted from code found on Rosetta Code
<http://rosettacode.org/wiki/Bitmap/Write_a_PPM_file#Fortran> which was authored
by Mauro Panigada. The authors of OpenMC have obtained permission from Mauro
Panigada to use and distribution this source code as part of OpenMC.

4
src/relaxng/geometry.rnc
View file

@ -6,8 +6,8 @@ element geometry {
(element universe { xsd:int } | attribute universe { xsd:int })? &
(
(element fill { xsd:int } | attribute fill { xsd:int }) |
(element material { ( xsd:int | "void" )+ } |
attribute material { ( xsd:int | "void" )+ })
(element material { list { ( xsd:int | "void" )+ } } |
attribute material { list { ( xsd:int | "void" )+ } })
) &
(element temperature { list { xsd:double+ } } |
attribute temperature { list { xsd:double+ } } )? &

28
src/relaxng/geometry.rng
View file

@ -47,20 +47,24 @@
</choice>
<choice>
<element name="material">
<oneOrMore>
<choice>
<data type="int"/>
<value>void</value>
</choice>
</oneOrMore>
<list>
<oneOrMore>
<choice>
<data type="int"/>
<value>void</value>
</choice>
</oneOrMore>
</list>
</element>
<attribute name="material">
<oneOrMore>
<choice>
<data type="int"/>
<value>void</value>
</choice>
</oneOrMore>
<list>
<oneOrMore>
<choice>
<data type="int"/>
<value>void</value>
</choice>
</oneOrMore>
</list>
</attribute>
</choice>
</choice>

43
src/relaxng/settings.rnc
View file

@ -1,29 +1,7 @@
element settings {
element confidence_intervals { xsd:boolean }? &
element batches { xsd:positiveInteger }? &
(
element eigenvalue {
(element batches { xsd:positiveInteger } |
attribute batches { xsd:positiveInteger }) &
(element inactive { xsd:nonNegativeInteger } |
attribute inactive { xsd:nonNegativeInteger }) &
(element particles { xsd:positiveInteger } |
attribute particles { xsd:positiveInteger }) &
(element generations_per_batch { xsd:positiveInteger } |
attribute generations_per_batch { xsd:positiveInteger })? &
(element keff_trigger {
(element type { xsd:string } | attribute type { xsd:string }) &
(element threshold { xsd:double} | attribute threshold { xsd:double })
}
)?
} |
element fixed_source {
(element batches { xsd:positiveInteger } |
attribute batches { xsd:positiveInteger }) &
(element particles { xsd:positiveInteger } |
attribute particles { xsd:positiveInteger })
}
) &
element confidence_intervals { xsd:boolean }? &
element cutoff {
(element weight { xsd:double } | attribute weight { xsd:double })? &
@ -43,12 +21,19 @@ element settings {
attribute upper_right { list { xsd:double+ } })
}? &
element generations_per_batch { xsd:positiveInteger }? &
element inactive { xsd:nonNegativeInteger }? &
element keff_trigger {
(element type { xsd:string } | attribute type { xsd:string }) &
(element threshold { xsd:double} | attribute threshold { xsd:double })
}? &
element log_grid_bins { xsd:positiveInteger }? &
element max_order { xsd:nonNegativeInteger }? &
element natural_elements { xsd:string { maxLength = "20" } }? &
element no_reduce { xsd:boolean }? &
element output {
@ -60,10 +45,14 @@ element settings {
element output_path { xsd:string { maxLength = "255" } }? &
element particles { xsd:positiveInteger }? &
element ptables { xsd:boolean }? &
element run_cmfd { xsd:boolean }? &
element run_mode { xsd:string }? &
element seed { xsd:positiveInteger }? &
element source {
@ -179,5 +168,5 @@ element settings {
(element E_max { xsd:double } |
attribute E_max { xsd:double })?
}*
}? &
}?
}

150
src/relaxng/settings.rng
View file

@ -1,107 +1,16 @@
<?xml version="1.0" encoding="UTF-8"?>
<element name="settings" xmlns="http://relaxng.org/ns/structure/1.0" datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
<interleave>
<optional>
<element name="batches">
<data type="positiveInteger"/>
</element>
</optional>
<optional>
<element name="confidence_intervals">
<data type="boolean"/>
</element>
</optional>
<choice>
<element name="eigenvalue">
<interleave>
<choice>
<element name="batches">
<data type="positiveInteger"/>
</element>
<attribute name="batches">
<data type="positiveInteger"/>
</attribute>
</choice>
<choice>
<element name="inactive">
<data type="nonNegativeInteger"/>
</element>
<attribute name="inactive">
<data type="nonNegativeInteger"/>
</attribute>
</choice>
<choice>
<element name="particles">
<data type="positiveInteger"/>
</element>
<attribute name="particles">
<data type="positiveInteger"/>
</attribute>
</choice>
<optional>
<choice>
<element name="generations_per_batch">
<data type="positiveInteger"/>
</element>
<attribute name="generations_per_batch">
<data type="positiveInteger"/>
</attribute>
</choice>
</optional>
<optional>
<element name="keff_trigger">
<interleave>
<choice>
<element name="type">
<data type="string"/>
</element>
<attribute name="type">
<data type="string"/>
</attribute>
</choice>
<choice>
<element name="threshold">
<data type="double"/>
</element>
<attribute name="threshold">
<data type="double"/>
</attribute>
</choice>
</interleave>
</element>
</optional>
</interleave>
</element>
<element name="fixed_source">
<interleave>
<choice>
<element name="batches">
<data type="positiveInteger"/>
</element>
<attribute name="batches">
<data type="positiveInteger"/>
</attribute>
</choice>
<choice>
<element name="particles">
<data type="positiveInteger"/>
</element>
<attribute name="particles">
<data type="positiveInteger"/>
</attribute>
</choice>
</interleave>
</element>
</choice>
<optional>
<element name="cross_sections">
<data type="string">
<param name="maxLength">255</param>
</data>
</element>
</optional>
<optional>
<element name="multipole_library">
<data type="string">
<param name="maxLength">255</param>
</data>
</element>
</optional>
<optional>
<element name="cutoff">
<interleave>
@ -208,6 +117,38 @@
</interleave>
</element>
</optional>
<optional>
<element name="generations_per_batch">
<data type="positiveInteger"/>
</element>
</optional>
<optional>
<element name="inactive">
<data type="nonNegativeInteger"/>
</element>
</optional>
<optional>
<element name="keff_trigger">
<interleave>
<choice>
<element name="type">
<data type="string"/>
</element>
<attribute name="type">
<data type="string"/>
</attribute>
</choice>
<choice>
<element name="threshold">
<data type="double"/>
</element>
<attribute name="threshold">
<data type="double"/>
</attribute>
</choice>
</interleave>
</element>
</optional>
<optional>
<element name="log_grid_bins">
<data type="positiveInteger"/>
@ -218,13 +159,6 @@
<data type="nonNegativeInteger"/>
</element>
</optional>
<optional>
<element name="natural_elements">
<data type="string">
<param name="maxLength">20</param>
</data>
</element>
</optional>
<optional>
<element name="no_reduce">
<data type="boolean"/>
@ -273,6 +207,11 @@
</data>
</element>
</optional>
<optional>
<element name="particles">
<data type="positiveInteger"/>
</element>
</optional>
<optional>
<element name="ptables">
<data type="boolean"/>
@ -283,6 +222,11 @@
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="run_mode">
<data type="string"/>
</element>
</optional>
<optional>
<element name="seed">
<data type="positiveInteger"/>

6
src/scattdata_header.F90
View file

@ -522,7 +522,7 @@ contains
gout = this % gmin(gin)
prob = this % energy(gin) % data(gout)
do while (prob < xi)
do while ((prob < xi) .and. (gout < this % gmax(gin)))
gout = gout + 1
prob = prob + this % energy(gin) % data(gout)
end do
@ -568,7 +568,7 @@ contains
gout = this % gmin(gin)
prob = this % energy(gin) % data(gout)
do while (prob < xi)
do while ((prob < xi) .and. (gout < this % gmax(gin)))
gout = gout + 1
prob = prob + this % energy(gin) % data(gout)
end do
@ -605,7 +605,7 @@ contains
gout = this % gmin(gin)
prob = this % energy(gin) % data(gout)
do while (prob < xi)
do while ((prob < xi) .and. (gout < this % gmax(gin)))
gout = gout + 1
prob = prob + this % energy(gin) % data(gout)
end do

66
src/simulation.F90
View file

@ -1,9 +1,5 @@
module simulation
#ifdef MPI
use message_passing
#endif
use cmfd_execute, only: cmfd_init_batch, execute_cmfd
use constants, only: ZERO
use eigenvalue, only: count_source_for_ufs, calculate_average_keff, &
@ -13,14 +9,17 @@ module simulation
use eigenvalue, only: join_bank_from_threads
#endif
use global
use message_passing
use output, only: write_message, header, print_columns, &
print_batch_keff, print_generation
print_batch_keff, print_generation, print_runtime, &
print_results, print_overlap_check, write_tallies
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use source, only: initialize_source, sample_external_source
use state_point, only: write_state_point, write_source_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies
use tally, only: synchronize_tallies, setup_active_usertallies, &
tally_statistics
use trigger, only: check_triggers
use tracking, only: transport
use volume_calc, only: run_volume_calculations
@ -42,9 +41,6 @@ contains
type(Particle) :: p
integer(8) :: i_work
! Volume calculations
if (size(volume_calcs) > 0) call run_volume_calculations()
if (.not. restart_run) call initialize_source()
! Display header
@ -116,6 +112,8 @@ contains
if (master) call header("SIMULATION FINISHED", level=1)
call finalize_simulation()
! Clear particle
call p % clear()
@ -322,7 +320,7 @@ contains
if (master) call check_triggers()
#ifdef MPI
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
MPI_COMM_WORLD, mpi_err)
mpi_intracomm, mpi_err)
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
@ -382,4 +380,52 @@ contains
end subroutine replay_batch_history
!===============================================================================
! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays
! execution time and results
!===============================================================================
subroutine finalize_simulation
! Start finalization timer
call time_finalize%start()
! Calculate statistics for tallies and write to tallies.out
if (master) then
if (n_realizations > 1) call tally_statistics()
end if
if (output_tallies) then
if (master) call write_tallies()
end if
if (check_overlaps) call reduce_overlap_count()
! Stop timers and show timing statistics
call time_finalize%stop()
call time_total%stop()
if (master) then
call print_runtime()
call print_results()
if (check_overlaps) call print_overlap_check()
end if
end subroutine finalize_simulation
!===============================================================================
! REDUCE_OVERLAP_COUNT accumulates cell overlap check counts to master
!===============================================================================
subroutine reduce_overlap_count()
#ifdef MPI
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, overlap_check_cnt, n_cells, &
MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err)
else
call MPI_REDUCE(overlap_check_cnt, overlap_check_cnt, n_cells, &
MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err)
end if
#endif
end subroutine reduce_overlap_count
end module simulation

1
src/source.F90
View file

@ -15,6 +15,7 @@ module source
use geometry_header, only: BASE_UNIVERSE
use global
use hdf5_interface, only: file_create, file_open, file_close, read_dataset
use message_passing, only: rank
use output, only: write_message
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_set_stream

320
src/state_point.F90
View file

@ -11,26 +11,21 @@ module state_point
! intervals, using the <state_point ... /> tag.
!===============================================================================
use, intrinsic :: ISO_C_BINDING, only: c_loc, c_ptr
use hdf5
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning
use global
use hdf5_interface
use mesh_header, only: RegularMesh
use message_passing
use output, only: write_message, time_stamp
use random_lcg, only: seed
use string, only: to_str, count_digits, zero_padded
use tally_header, only: TallyObject
use mesh_header, only: RegularMesh
use dict_header, only: ElemKeyValueII, ElemKeyValueCI
use random_lcg, only: seed
#ifdef MPI
use message_passing
#endif
use hdf5
use, intrinsic :: ISO_C_BINDING, only: c_loc, c_ptr
implicit none
@ -47,17 +42,13 @@ contains
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
integer, allocatable :: id_array(:)
integer, allocatable :: key_array(:)
integer(HID_T) :: file_id
integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, &
mesh_group, filter_group, derivs_group, deriv_group, &
runtime_group
character(MAX_WORD_LEN), allocatable :: str_array(:)
character(MAX_FILE_LEN) :: filename
type(RegularMesh), pointer :: meshp
type(TallyObject), pointer :: tally
type(ElemKeyValueII), pointer :: current
type(ElemKeyValueII), pointer :: next
! Set filename for state point
filename = trim(path_output) // 'statepoint.' // &
@ -72,36 +63,37 @@ contains
file_id = file_create(filename)
! Write file type
call write_dataset(file_id, "filetype", 'statepoint')
call write_attribute(file_id, "filetype", "statepoint")
! Write revision number for state point file
call write_dataset(file_id, "revision", REVISION_STATEPOINT)
call write_attribute(file_id, "version", VERSION_STATEPOINT)
! Write OpenMC version
call write_dataset(file_id, "version_major", VERSION_MAJOR)
call write_dataset(file_id, "version_minor", VERSION_MINOR)
call write_dataset(file_id, "version_release", VERSION_RELEASE)
call write_attribute(file_id, "openmc_version", VERSION)
#ifdef GIT_SHA1
call write_attribute(file_id, "git_sha1", GIT_SHA1)
#endif
! Write current date and time
call write_dataset(file_id, "date_and_time", time_stamp())
call write_attribute(file_id, "date_and_time", time_stamp())
! Write path to input
call write_dataset(file_id, "path", path_input)
call write_attribute(file_id, "path", path_input)
! Write out random number seed
call write_dataset(file_id, "seed", seed)
! Write run information
if (run_CE) then
call write_dataset(file_id, "run_CE", 1)
call write_dataset(file_id, "energy_mode", "continuous-energy")
else
call write_dataset(file_id, "run_CE", 0)
call write_dataset(file_id, "energy_mode", "multi-group")
end if
select case(run_mode)
case (MODE_FIXEDSOURCE)
call write_dataset(file_id, "run_mode", "fixed source")
case (MODE_EIGENVALUE)
call write_dataset(file_id, "run_mode", "k-eigenvalue")
call write_dataset(file_id, "run_mode", "eigenvalue")
end select
call write_dataset(file_id, "n_particles", n_particles)
call write_dataset(file_id, "n_batches", n_batches)
@ -111,15 +103,15 @@ contains
! Indicate whether source bank is stored in statepoint
if (source_separate) then
call write_dataset(file_id, "source_present", 0)
call write_attribute(file_id, "source_present", 0)
else
call write_dataset(file_id, "source_present", 1)
call write_attribute(file_id, "source_present", 1)
end if
! Write out information for eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call write_dataset(file_id, "n_inactive", n_inactive)
call write_dataset(file_id, "gen_per_batch", gen_per_batch)
call write_dataset(file_id, "generations_per_batch", gen_per_batch)
call write_dataset(file_id, "k_generation", k_generation)
call write_dataset(file_id, "entropy", entropy)
call write_dataset(file_id, "k_col_abs", k_col_abs)
@ -129,7 +121,7 @@ contains
! Write out CMFD info
if (cmfd_on) then
call write_dataset(file_id, "cmfd_on", 1)
call write_attribute(file_id, "cmfd_on", 1)
cmfd_group = create_group(file_id, "cmfd")
call write_dataset(cmfd_group, "indices", cmfd % indices)
@ -141,7 +133,7 @@ contains
call write_dataset(cmfd_group, "cmfd_srccmp", cmfd % src_cmp)
call close_group(cmfd_group)
else
call write_dataset(file_id, "cmfd_on", 0)
call write_attribute(file_id, "cmfd_on", 0)
end if
end if
@ -149,52 +141,35 @@ contains
! Write number of meshes
meshes_group = create_group(tallies_group, "meshes")
call write_dataset(meshes_group, "n_meshes", n_meshes)
call write_attribute(meshes_group, "n_meshes", n_meshes)
if (n_meshes > 0) then
! Print list of mesh IDs
current => mesh_dict % keys()
! Write IDs of meshes
allocate(id_array(n_meshes))
allocate(key_array(n_meshes))
i = 1
do while (associated(current))
key_array(i) = current % key
id_array(i) = current % value
! Move to next mesh
next => current % next
deallocate(current)
current => next
i = i + 1
do i = 1, n_meshes
id_array(i) = meshes(i) % id
end do
call write_dataset(meshes_group, "ids", id_array)
call write_dataset(meshes_group, "keys", key_array)
deallocate(key_array)
call write_attribute(meshes_group, "ids", id_array)
deallocate(id_array)
! Write information for meshes
MESH_LOOP: do i = 1, n_meshes
meshp => meshes(id_array(i))
mesh_group = create_group(meshes_group, "mesh " &
// trim(to_str(meshp % id)))
associate (m => meshes(i))
mesh_group = create_group(meshes_group, "mesh " &
// trim(to_str(m % id)))
select case (meshp % type)
case (MESH_REGULAR)
call write_dataset(mesh_group, "type", "regular")
end select
call write_dataset(mesh_group, "dimension", meshp % dimension)
call write_dataset(mesh_group, "lower_left", meshp % lower_left)
call write_dataset(mesh_group, "upper_right", meshp % upper_right)
call write_dataset(mesh_group, "width", meshp % width)
select case (m % type)
case (MESH_REGULAR)
call write_dataset(mesh_group, "type", "regular")
end select
call write_dataset(mesh_group, "dimension", m % dimension)
call write_dataset(mesh_group, "lower_left", m % lower_left)
call write_dataset(mesh_group, "upper_right", m % upper_right)
call write_dataset(mesh_group, "width", m % width)
call close_group(mesh_group)
call close_group(mesh_group)
end associate
end do MESH_LOOP
deallocate(id_array)
end if
call close_group(meshes_group)
@ -232,25 +207,16 @@ contains
end if
! Write number of tallies
call write_dataset(tallies_group, "n_tallies", n_tallies)
call write_attribute(tallies_group, "n_tallies", n_tallies)
if (n_tallies > 0) then
! Print list of tally IDs
! Write array of tally IDs
allocate(id_array(n_tallies))
allocate(key_array(n_tallies))
! Write all tally information except results
do i = 1, n_tallies
tally => tallies(i)
key_array(i) = tally % id
id_array(i) = i
id_array(i) = tallies(i) % id
end do
call write_dataset(tallies_group, "ids", id_array)
call write_dataset(tallies_group, "keys", key_array)
deallocate(key_array)
call write_attribute(tallies_group, "ids", id_array)
deallocate(id_array)
! Write all tally information except results
TALLY_METADATA: do i = 1, n_tallies
@ -260,6 +226,9 @@ contains
tally_group = create_group(tallies_group, "tally " // &
trim(to_str(tally % id)))
! Write the name for this tally
call write_dataset(tally_group, "name", tally % name)
select case(tally % estimator)
case (ESTIMATOR_ANALOG)
call write_dataset(tally_group, "estimator", "analog")
@ -376,14 +345,14 @@ contains
call write_dataset(file_id, "n_realizations", n_realizations)
! Write global tallies
call write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES)
call write_dataset(file_id, "global_tallies", global_tallies)
! Write tallies
tallies_group = open_group(file_id, "tallies")
if (tallies_on) then
! Indicate that tallies are on
call write_dataset(tallies_group, "tallies_present", 1)
call write_attribute(file_id, "tallies_present", 1)
tallies_group = open_group(file_id, "tallies")
! Write all tally results
TALLY_RESULTS: do i = 1, n_tallies
@ -397,12 +366,12 @@ contains
call close_group(tally_group)
end do TALLY_RESULTS
call close_group(tallies_group)
else
! Indicate tallies are off
call write_dataset(tallies_group, "tallies_present", 0)
call write_attribute(file_id, "tallies_present", 0)
end if
call close_group(tallies_group)
! Write out the runtime metrics.
runtime_group = create_group(file_id, "runtime")
@ -525,10 +494,6 @@ contains
#ifdef MPI
real(8) :: dummy ! temporary receive buffer for non-root reduces
#endif
integer, allocatable :: id_array(:)
type(ElemKeyValueII), pointer :: current
type(ElemKeyValueII), pointer :: next
type(TallyObject), pointer :: tally
type(TallyObject) :: dummy_tally
! ==========================================================================
@ -553,7 +518,7 @@ contains
! receive buffer without having a temporary variable
#ifdef MPI
call MPI_REDUCE(MPI_IN_PLACE, global_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
0, mpi_intracomm, mpi_err)
#endif
! Transfer values to value on master
@ -567,99 +532,79 @@ contains
! Receive buffer not significant at other processors
#ifdef MPI
call MPI_REDUCE(global_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
0, mpi_intracomm, mpi_err)
#endif
end if
if (tallies_on) then
! Indicate that tallies are on
if (master) then
call write_dataset(tallies_group, "tallies_present", 1)
! Build list of tally IDs
current => tally_dict%keys()
allocate(id_array(n_tallies))
i = 1
do while (associated(current))
id_array(i) = current%value
! Move to next tally
next => current%next
deallocate(current)
current => next
i = i + 1
end do
call write_attribute(file_id, "tallies_present", 1)
end if
! Write all tally results
TALLY_RESULTS: do i = 1, n_tallies
associate (t => tallies(i))
! Determine size of tally results array
m = size(t % results, 2)
n = size(t % results, 3)
n_bins = m*n*2
tally => tallies(i)
! Allocate array for storing sums and sums of squares, but
! contiguously in memory for each
allocate(tally_temp(2,m,n))
tally_temp(1,:,:) = t % results(RESULT_SUM,:,:)
tally_temp(2,:,:) = t % results(RESULT_SUM_SQ,:,:)
! Determine size of tally results array
m = size(tally%results, 2)
n = size(tally%results, 3)
n_bins = m*n*2
if (master) then
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(t % id)))
! Allocate array for storing sums and sums of squares, but
! contiguously in memory for each
allocate(tally_temp(2,m,n))
tally_temp(1,:,:) = tally%results(RESULT_SUM,:,:)
tally_temp(2,:,:) = tally%results(RESULT_SUM_SQ,:,:)
if (master) then
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(tally%id)))
! The MPI_IN_PLACE specifier allows the master to copy values into
! a receive buffer without having a temporary variable
! The MPI_IN_PLACE specifier allows the master to copy values into
! a receive buffer without having a temporary variable
#ifdef MPI
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, &
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, &
MPI_SUM, 0, mpi_intracomm, mpi_err)
#endif
! At the end of the simulation, store the results back in the
! regular TallyResults array
if (current_batch == n_max_batches .or. satisfy_triggers) then
tally%results(RESULT_SUM,:,:) = tally_temp(1,:,:)
tally%results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:)
! At the end of the simulation, store the results back in the
! regular TallyResults array
if (current_batch == n_max_batches .or. satisfy_triggers) then
t % results(RESULT_SUM,:,:) = tally_temp(1,:,:)
t % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:)
end if
! Put in temporary tally result
allocate(dummy_tally % results(3,m,n))
dummy_tally % results(RESULT_SUM,:,:) = tally_temp(1,:,:)
dummy_tally % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:)
! Write reduced tally results to file
call dummy_tally % write_results_hdf5(tally_group)
! Deallocate temporary tally result
deallocate(dummy_tally % results)
else
! Receive buffer not significant at other processors
#ifdef MPI
call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, &
0, mpi_intracomm, mpi_err)
#endif
end if
! Put in temporary tally result
allocate(dummy_tally % results(3,m,n))
dummy_tally % results(RESULT_SUM,:,:) = tally_temp(1,:,:)
dummy_tally % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:)
! Deallocate temporary copy of tally results
deallocate(tally_temp)
! Write reduced tally results to file
call dummy_tally % write_results_hdf5(tally_group)
! Deallocate temporary tally result
deallocate(dummy_tally % results)
else
! Receive buffer not significant at other processors
#ifdef MPI
call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
#endif
end if
! Deallocate temporary copy of tally results
deallocate(tally_temp)
if (master) call close_group(tally_group)
if (master) call close_group(tally_group)
end associate
end do TALLY_RESULTS
deallocate(id_array)
if (master) call close_group(tallies_group)
else
if (master) then
! Indicate that tallies are off
call write_dataset(tallies_group, "tallies_present", 0)
end if
! Indicate that tallies are off
if (master) call write_dataset(file_id, "tallies_present", 0)
end if
if (master) call close_group(tallies_group)
end subroutine write_tally_results_nr
@ -671,15 +616,14 @@ contains
integer :: i
integer :: int_array(3)
integer, allocatable :: array(:)
integer(HID_T) :: file_id
integer(HID_T) :: cmfd_group
integer(HID_T) :: tallies_group
integer(HID_T) :: tally_group
real(8) :: real_array(3)
logical :: source_present
integer :: sp_run_CE
character(MAX_WORD_LEN) :: word
type(TallyObject), pointer :: tally
! Write message
call write_message("Loading state point " // trim(path_state_point) &
@ -689,15 +633,15 @@ contains
file_id = file_open(path_state_point, 'r', parallel=.true.)
! Read filetype
call read_dataset(word, file_id, "filetype")
call read_attribute(word, file_id, "filetype")
if (word /= 'statepoint') then
call fatal_error("OpenMC tried to restart from a non-statepoint file.")
end if
! Read revision number for state point file and make sure it matches with
! current version
call read_dataset(int_array(1), file_id, "revision")
if (int_array(1) /= REVISION_STATEPOINT) then
call read_attribute(array, file_id, "version")
if (any(array /= VERSION_STATEPOINT)) then
call fatal_error("State point version does not match current version &
&in OpenMC.")
end if
@ -707,11 +651,11 @@ contains
! It is not impossible for a state point to be generated from a CE run but
! to be loaded in to an MG run (or vice versa), check to prevent that.
call read_dataset(sp_run_CE, file_id, "run_CE")
if (sp_run_CE == 0 .and. run_CE) then
call read_dataset(word, file_id, "energy_mode")
if (word == "multi-group" .and. run_CE) then
call fatal_error("State point file is from multi-group run but &
& current run is continous-energy!")
else if (sp_run_CE == 1 .and. .not. run_CE) then
else if (word == "continuous-energy" .and. .not. run_CE) then
call fatal_error("State point file is from continuous-energy run but &
& current run is multi-group!")
end if
@ -721,7 +665,7 @@ contains
select case(word)
case ('fixed source')
run_mode = MODE_FIXEDSOURCE
case ('k-eigenvalue')
case ('eigenvalue')
run_mode = MODE_EIGENVALUE
end select
call read_dataset(n_particles, file_id, "n_particles")
@ -734,7 +678,7 @@ contains
call read_dataset(restart_batch, file_id, "current_batch")
! Check for source in statepoint if needed
call read_dataset(int_array(1), file_id, "source_present")
call read_attribute(int_array(1), file_id, "source_present")
if (int_array(1) == 1) then
source_present = .true.
else
@ -749,7 +693,7 @@ contains
! Read information specific to eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call read_dataset(int_array(1), file_id, "n_inactive")
call read_dataset(gen_per_batch, file_id, "gen_per_batch")
call read_dataset(gen_per_batch, file_id, "generations_per_batch")
call read_dataset(k_generation(1:restart_batch*gen_per_batch), &
file_id, "k_generation")
call read_dataset(entropy(1:restart_batch*gen_per_batch), &
@ -763,7 +707,7 @@ contains
n_inactive = max(n_inactive, int_array(1))
! Read in to see if CMFD was on
call read_dataset(int_array(1), file_id, "cmfd_on")
call read_attribute(int_array(1), file_id, "cmfd_on")
! Read in CMFD info
if (int_array(1) == 1) then
@ -804,30 +748,28 @@ contains
call read_dataset(global_tallies, file_id, "global_tallies")
! Check if tally results are present
tallies_group = open_group(file_id, "tallies")
call read_dataset(int_array(1), tallies_group, "tallies_present", &
indep=.true.)
call read_attribute(int_array(1), file_id, "tallies_present")
! Read in sum and sum squared
if (int_array(1) == 1) then
tallies_group = open_group(file_id, "tallies")
TALLY_RESULTS: do i = 1, n_tallies
! Set pointer to tally
tally => tallies(i)
! Read sum, sum_sq, and N for each bin
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(tally % id)))
call tally % read_results_hdf5(tally_group)
call read_dataset(tally % n_realizations, tally_group, &
"n_realizations")
call close_group(tally_group)
associate (t => tallies(i))
! Read sum, sum_sq, and N for each bin
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(t % id)))
call t % read_results_hdf5(tally_group)
call read_dataset(t % n_realizations, tally_group, &
"n_realizations")
call close_group(tally_group)
end associate
end do TALLY_RESULTS
call close_group(tallies_group)
end if
call close_group(tallies_group)
end if
! Read source if in eigenvalue mode
if (run_mode == MODE_EIGENVALUE) then
@ -940,7 +882,7 @@ contains
! Receive source sites from other processes
if (i > 0) then
call MPI_RECV(source_bank, int(dims(1)), MPI_BANK, i, i, &
MPI_COMM_WORLD, MPI_STATUS_IGNORE, mpi_err)
mpi_intracomm, MPI_STATUS_IGNORE, mpi_err)
end if
#endif
@ -969,7 +911,7 @@ contains
else
#ifdef MPI
call MPI_SEND(source_bank, int(work), MPI_BANK, 0, rank, &
MPI_COMM_WORLD, mpi_err)
mpi_intracomm, mpi_err)
#endif
end if

159
src/summary.F90
View file

@ -1,5 +1,7 @@
module summary
use hdf5
use constants
use endf, only: reaction_name
use geometry_header, only: BASE_UNIVERSE, Cell, Universe, Lattice, &
@ -8,6 +10,7 @@ module summary
use hdf5_interface
use material_header, only: Material
use mesh_header, only: RegularMesh
use message_passing
use nuclide_header
use output, only: time_stamp
use surface_header
@ -15,8 +18,6 @@ module summary
use tally_header, only: TallyObject
use tally_filter, only: find_offset
use hdf5
implicit none
private
@ -35,45 +36,10 @@ contains
! Create a new file using default properties.
file_id = file_create("summary.h5")
! Write header information
call write_header(file_id)
if (run_CE) then
call write_dataset(file_id, "run_CE", 1)
else
call write_dataset(file_id, "run_CE", 0)
end if
! Write number of particles
call write_dataset(file_id, "n_particles", n_particles)
call write_dataset(file_id, "n_batches", n_batches)
call write_attribute_string(file_id, "n_particles", &
"description", "Number of particles per generation")
call write_attribute_string(file_id, "n_batches", &
"description", "Total number of batches")
! Write eigenvalue information
if (run_mode == MODE_EIGENVALUE) then
! write number of inactive/active batches and generations/batch
call write_dataset(file_id, "n_inactive", n_inactive)
call write_dataset(file_id, "n_active", n_active)
call write_dataset(file_id, "gen_per_batch", gen_per_batch)
! Add description of each variable
call write_attribute_string(file_id, "n_inactive", &
"description", "Number of inactive batches")
call write_attribute_string(file_id, "n_active", &
"description", "Number of active batches")
call write_attribute_string(file_id, "gen_per_batch", &
"description", "Number of generations per batch")
end if
call write_nuclides(file_id)
call write_geometry(file_id)
call write_materials(file_id)
if (n_tallies > 0) then
call write_tallies(file_id)
end if
! Terminate access to the file.
call file_close(file_id)
@ -87,22 +53,16 @@ contains
subroutine write_header(file_id)
integer(HID_T), intent(in) :: file_id
! Write filetype and revision
call write_dataset(file_id, "filetype", "summary")
call write_dataset(file_id, "revision", REVISION_SUMMARY)
! Write version information
call write_dataset(file_id, "version_major", VERSION_MAJOR)
call write_dataset(file_id, "version_minor", VERSION_MINOR)
call write_dataset(file_id, "version_release", VERSION_RELEASE)
! Write filetype and version info
call write_attribute(file_id, "filetype", "summary")
call write_attribute(file_id, "version", VERSION_SUMMARY)
call write_attribute(file_id, "openmc_version", VERSION)
#ifdef GIT_SHA1
call write_attribute(file_id, "git_sha1", GIT_SHA1)
#endif
! Write current date and time
call write_dataset(file_id, "date_and_time", time_stamp())
! Write MPI information
call write_dataset(file_id, "n_procs", n_procs)
call write_attribute_string(file_id, "n_procs", "description", &
"Number of MPI processes")
call write_attribute(file_id, "date_and_time", time_stamp())
end subroutine write_header
@ -119,7 +79,7 @@ contains
! Write useful data from nuclide objects
nuclide_group = create_group(file_id, "nuclides")
call write_dataset(nuclide_group, "n_nuclides_total", n_nuclides_total)
call write_attribute(nuclide_group, "n_nuclides", n_nuclides_total)
! Build array of nuclide names and awrs
allocate(nucnames(n_nuclides_total))
@ -151,9 +111,10 @@ contains
subroutine write_geometry(file_id)
integer(HID_T), intent(in) :: file_id
integer :: i, j, k, m, offset
integer :: i, j, k, m, offset
integer, allocatable :: lattice_universes(:,:,:)
integer, allocatable :: cell_materials(:)
integer, allocatable :: cell_ids(:)
real(8), allocatable :: cell_temperatures(:)
integer(HID_T) :: geom_group
integer(HID_T) :: cells_group, cell_group
@ -171,10 +132,10 @@ contains
! Use H5LT interface to write number of geometry objects
geom_group = create_group(file_id, "geometry")
call write_dataset(geom_group, "n_cells", n_cells)
call write_dataset(geom_group, "n_surfaces", n_surfaces)
call write_dataset(geom_group, "n_universes", n_universes)
call write_dataset(geom_group, "n_lattices", n_lattices)
call write_attribute(geom_group, "n_cells", n_cells)
call write_attribute(geom_group, "n_surfaces", n_surfaces)
call write_attribute(geom_group, "n_universes", n_universes)
call write_attribute(geom_group, "n_lattices", n_lattices)
! ==========================================================================
! WRITE INFORMATION ON CELLS
@ -187,9 +148,6 @@ contains
c => cells(i)
cell_group = create_group(cells_group, "cell " // trim(to_str(c%id)))
! Write internal OpenMC index for this cell
call write_dataset(cell_group, "index", i)
! Write name for this cell
call write_dataset(cell_group, "name", c%name)
@ -198,8 +156,8 @@ contains
! Write information on what fills this cell
select case (c%type)
case (CELL_NORMAL)
call write_dataset(cell_group, "fill_type", "normal")
case (FILL_MATERIAL)
call write_dataset(cell_group, "fill_type", "material")
if (size(c % material) == 1) then
if (c % material(1) == MATERIAL_VOID) then
@ -227,7 +185,7 @@ contains
call write_dataset(cell_group, "temperature", cell_temperatures)
deallocate(cell_temperatures)
case (CELL_FILL)
case (FILL_UNIVERSE)
call write_dataset(cell_group, "fill_type", "universe")
call write_dataset(cell_group, "fill", universes(c%fill)%id)
if (allocated(c%offset)) then
@ -243,7 +201,7 @@ contains
call write_dataset(cell_group, "rotation", c%rotation)
end if
case (CELL_LATTICE)
case (FILL_LATTICE)
call write_dataset(cell_group, "fill_type", "lattice")
call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id)
end select
@ -302,9 +260,6 @@ contains
surface_group = create_group(surfaces_group, "surface " // &
trim(to_str(s%id)))
! Write internal OpenMC index for this surface
call write_dataset(surface_group, "index", i)
! Write name for this surface
call write_dataset(surface_group, "name", s%name)
@ -374,16 +329,16 @@ contains
call write_dataset(surface_group, "coefficients", coeffs)
deallocate(coeffs)
! Write boundary condition
! Write boundary type
select case (s%bc)
case (BC_TRANSMIT)
call write_dataset(surface_group, "boundary_condition", "transmission")
call write_dataset(surface_group, "boundary_type", "transmission")
case (BC_VACUUM)
call write_dataset(surface_group, "boundary_condition", "vacuum")
call write_dataset(surface_group, "boundary_type", "vacuum")
case (BC_REFLECT)
call write_dataset(surface_group, "boundary_condition", "reflective")
call write_dataset(surface_group, "boundary_type", "reflective")
case (BC_PERIODIC)
call write_dataset(surface_group, "boundary_condition", "periodic")
call write_dataset(surface_group, "boundary_type", "periodic")
end select
call close_group(surface_group)
@ -403,11 +358,15 @@ contains
univ_group = create_group(universes_group, "universe " // &
trim(to_str(u%id)))
! Write internal OpenMC index for this universe
call write_dataset(univ_group, "index", i)
! Write list of cells in this universe
if (u%n_cells > 0) call write_dataset(univ_group, "cells", u%cells)
if (u % n_cells > 0) then
allocate(cell_ids(u % n_cells))
do j = 1, u % n_cells
cell_ids(j) = cells(u % cells(j)) % id
end do
call write_dataset(univ_group, "cells", cell_ids)
deallocate(cell_ids)
end if
call close_group(univ_group)
end do UNIVERSE_LOOP
@ -425,13 +384,14 @@ contains
lat => lattices(i)%obj
lattice_group = create_group(lattices_group, "lattice " // trim(to_str(lat%id)))
! Write internal OpenMC index for this lattice
call write_dataset(lattice_group, "index", i)
! Write name, pitch, and outer universe
call write_dataset(lattice_group, "name", lat%name)
call write_dataset(lattice_group, "pitch", lat%pitch)
call write_dataset(lattice_group, "outer", lat%outer)
if (lat % outer > 0) then
call write_dataset(lattice_group, "outer", universes(lat % outer) % id)
else
call write_dataset(lattice_group, "outer", lat % outer)
end if
! Write distribcell offsets if present
if (allocated(lat%offset)) then
@ -535,16 +495,17 @@ contains
material_group = create_group(materials_group, "material " // &
trim(to_str(m%id)))
! Write internal OpenMC index for this material
call write_dataset(material_group, "index", i)
if (m % depletable) then
call write_attribute(material_group, "depletable", 1)
else
call write_attribute(material_group, "depletable", 0)
end if
! Write name for this material
call write_dataset(material_group, "name", m % name)
! Write atom density with units
call write_dataset(material_group, "atom_density", m % density)
call write_attribute_string(material_group, "atom_density", "units", &
"atom/b-cm")
! Copy ZAID for each nuclide to temporary array
allocate(nucnames(m%n_nuclides))
@ -576,34 +537,4 @@ contains
end subroutine write_materials
!===============================================================================
! WRITE_TALLIES
!===============================================================================
subroutine write_tallies(file_id)
integer(HID_T), intent(in) :: file_id
integer :: i
integer(HID_T) :: tallies_group
integer(HID_T) :: tally_group
type(TallyObject), pointer :: t
tallies_group = create_group(file_id, "tallies")
TALLY_METADATA: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
tally_group = create_group(tallies_group, "tally " &
// trim(to_str(t % id)))
! Write the name for this tally
call write_dataset(tally_group, "name", t % name)
call close_group(tally_group)
end do TALLY_METADATA
call close_group(tallies_group)
end subroutine write_tallies
end module summary

19
src/tally.F90
View file

@ -2,10 +2,6 @@ module tally
use, intrinsic :: ISO_C_BINDING
#ifdef MPI
use message_passing
#endif
use algorithm, only: binary_search
use constants
use cross_section, only: multipole_deriv_eval
@ -18,6 +14,7 @@ module tally
mesh_intersects_1d, mesh_intersects_2d, &
mesh_intersects_3d
use mesh_header, only: RegularMesh
use message_passing
use output, only: header
use particle_header, only: LocalCoord, Particle
use string, only: to_str
@ -2030,7 +2027,7 @@ contains
end do SCORE_LOOP
nullify(matxs,nucxs)
nullify(matxs, nucxs)
end subroutine score_general_mg
!===============================================================================
@ -4073,14 +4070,14 @@ contains
! The MPI_IN_PLACE specifier allows the master to copy values into
! a receive buffer without having a temporary variable
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, &
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_SUM, 0, mpi_intracomm, mpi_err)
! Transfer values to value on master
t % results(RESULT_VALUE,:,:) = tally_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, &
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_SUM, 0, mpi_intracomm, mpi_err)
! Reset value on other processors
t % results(RESULT_VALUE,:,:) = ZERO
@ -4094,14 +4091,14 @@ contains
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err)
! Transfer values back to global_tallies on master
global_tallies(RESULT_VALUE, :) = global_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err)
! Reset value on other processors
global_tallies(RESULT_VALUE, :) = ZERO
@ -4111,11 +4108,11 @@ contains
! last realization
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
0, mpi_intracomm, mpi_err)
else
! Receive buffer not significant at other processors
call MPI_REDUCE(total_weight, dummy, 1, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
0, mpi_intracomm, mpi_err)
end if
end subroutine reduce_tally_results

71
src/tally_filter.F90
View file

@ -501,9 +501,17 @@ contains
class(UniverseFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
integer :: i
integer, allocatable :: universe_ids(:)
call write_dataset(filter_group, "type", "universe")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % universes )
allocate(universe_ids(size(this % universes)))
do i = 1, size(this % universes)
universe_ids(i) = universes(this % universes(i)) % id
end do
call write_dataset(filter_group, "bins", universe_ids)
end subroutine to_statepoint_universe
subroutine initialize_universe(this)
@ -562,9 +570,17 @@ contains
class(MaterialFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
integer :: i
integer, allocatable :: material_ids(:)
call write_dataset(filter_group, "type", "material")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % materials )
allocate(material_ids(size(this % materials)))
do i = 1, size(this % materials)
material_ids(i) = materials(this % materials(i)) % id
end do
call write_dataset(filter_group, "bins", material_ids)
end subroutine to_statepoint_material
subroutine initialize_material(this)
@ -639,9 +655,17 @@ contains
class(CellFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
integer :: i
integer, allocatable :: cell_ids(:)
call write_dataset(filter_group, "type", "cell")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % cells )
allocate(cell_ids(size(this % cells)))
do i = 1, size(this % cells)
cell_ids(i) = cells(this % cells(i)) % id
end do
call write_dataset(filter_group, "bins", cell_ids)
end subroutine to_statepoint_cell
subroutine initialize_cell(this)
@ -692,10 +716,10 @@ contains
distribcell_index = cells(this % cell) % distribcell_index
offset = 0
do i = 1, p % n_coord
if (cells(p % coord(i) % cell) % type == CELL_FILL) then
if (cells(p % coord(i) % cell) % type == FILL_UNIVERSE) then
offset = offset + cells(p % coord(i) % cell) % &
offset(distribcell_index)
elseif (cells(p % coord(i) % cell) % type == CELL_LATTICE) then
elseif (cells(p % coord(i) % cell) % type == FILL_LATTICE) then
if (lattices(p % coord(i + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(i + 1) % lattice_x, &
@ -725,7 +749,7 @@ contains
call write_dataset(filter_group, "type", "distribcell")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % cell )
call write_dataset(filter_group, "bins", cells(this % cell) % id)
end subroutine to_statepoint_distribcell
subroutine initialize_distribcell(this)
@ -784,9 +808,16 @@ contains
class(CellbornFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
integer :: i
integer, allocatable :: cell_ids(:)
call write_dataset(filter_group, "type", "cellborn")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % cells )
allocate(cell_ids(size(this % cells)))
do i = 1, size(this % cells)
cell_ids(i) = cells(this % cells(i)) % id
end do
call write_dataset(filter_group, "bins", cell_ids)
end subroutine to_statepoint_cellborn
subroutine initialize_cellborn(this)
@ -852,7 +883,7 @@ contains
call write_dataset(filter_group, "type", "surface")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % surfaces )
call write_dataset(filter_group, "bins", this % surfaces)
end subroutine to_statepoint_surface
subroutine initialize_surface(this)
@ -941,7 +972,7 @@ contains
call write_dataset(filter_group, "type", "energy")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % bins )
call write_dataset(filter_group, "bins", this % bins)
end subroutine to_statepoint_energy
function text_label_energy(this, bin) result(label)
@ -1006,7 +1037,7 @@ contains
call write_dataset(filter_group, "type", "energyout")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % bins )
call write_dataset(filter_group, "bins", this % bins)
end subroutine to_statepoint_energyout
function text_label_energyout(this, bin) result(label)
@ -1048,7 +1079,7 @@ contains
call write_dataset(filter_group, "type", "delayedgroup")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % groups )
call write_dataset(filter_group, "bins", this % groups)
end subroutine to_statepoint_dg
function text_label_dg(this, bin) result(label)
@ -1097,7 +1128,7 @@ contains
call write_dataset(filter_group, "type", "mu")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % bins )
call write_dataset(filter_group, "bins", this % bins)
end subroutine to_statepoint_mu
function text_label_mu(this, bin) result(label)
@ -1160,7 +1191,7 @@ contains
call write_dataset(filter_group, "type", "polar")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % bins )
call write_dataset(filter_group, "bins", this % bins)
end subroutine to_statepoint_polar
function text_label_polar(this, bin) result(label)
@ -1223,7 +1254,7 @@ contains
call write_dataset(filter_group, "type", "azimuthal")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % bins )
call write_dataset(filter_group, "bins", this % bins)
end subroutine to_statepoint_azimuthal
function text_label_azimuthal(this, bin) result(label)
@ -1390,7 +1421,7 @@ contains
c => cells(cell_index)
! Skip normal cells which do not have offsets
if (c % type == CELL_NORMAL) then
if (c % type == FILL_MATERIAL) then
cycle
end if
@ -1399,13 +1430,13 @@ contains
end do
! Ensure we didn't just end the loop by iteration
if (c % type /= CELL_NORMAL) then
if (c % type /= FILL_MATERIAL) then
! There are more cells in this universe that it could be in
later_cell = .true.
! Two cases, lattice or fill cell
if (c % type == CELL_FILL) then
if (c % type == FILL_UNIVERSE) then
temp_offset = c % offset(map)
! Get the offset of the first lattice location
@ -1422,7 +1453,7 @@ contains
end if
end if
if (n == 1 .and. c % type /= CELL_NORMAL) then
if (n == 1 .and. c % type /= FILL_MATERIAL) then
this_cell = .true.
end if
@ -1440,7 +1471,7 @@ contains
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type == CELL_FILL) then
if (c % type == FILL_UNIVERSE) then
! Enter this cell to update the current offset
offset = c % offset(map) + offset
@ -1451,7 +1482,7 @@ contains
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
elseif (c % type == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj

4
src/timer_header.F90
View file

@ -42,11 +42,11 @@ contains
function timer_get_value(self) result(elapsed)
class(Timer), intent(in) :: self ! the timer
real(8) :: elapsed ! total elapsed time
real(8) :: elapsed ! total elapsed time
integer(8) :: end_counts ! current number of counts
integer(8) :: count_rate ! system-dependent counting rate
real :: elapsed_time ! elapsed time since last start
real(8) :: elapsed_time ! elapsed time since last start
if (self % running) then
call system_clock(end_counts, count_rate)

8
src/track_output.F90
View file

@ -114,10 +114,10 @@ contains
!$omp critical (FinalizeParticleTrack)
file_id = file_create(fname)
call write_dataset(file_id, 'filetype', 'track')
call write_dataset(file_id, 'revision', REVISION_TRACK)
call write_dataset(file_id, 'n_particles', n_particle_tracks)
call write_dataset(file_id, 'n_coords', n_coords)
call write_attribute(file_id, 'filetype', 'track')
call write_attribute(file_id, 'version', VERSION_TRACK)
call write_attribute(file_id, 'n_particles', n_particle_tracks)
call write_attribute(file_id, 'n_coords', n_coords)
do i = 1, n_particle_tracks
call write_dataset(file_id, 'coordinates_' // trim(to_str(i)), &
tracks(i)%coords)

3
src/tracking.F90
View file

@ -8,6 +8,7 @@ module tracking
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use output, only: write_message
use message_passing
use particle_header, only: LocalCoord, Particle
use physics, only: collision
use physics_mg, only: collision_mg
@ -104,9 +105,7 @@ contains
p % coord(p % n_coord) % uvw, material_xs)
else
material_xs % total = ZERO
material_xs % elastic = ZERO
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
end if
end if

1
src/trigger.F90
View file

@ -10,6 +10,7 @@ module trigger
use output, only: warning, write_message
use mesh, only: bin_to_mesh_indices
use mesh_header, only: RegularMesh
use message_passing, only: master
use trigger_header, only: TriggerObject
use tally, only: TallyObject
use tally_filter, only: MeshFilter

29
src/volume_calc.F90
View file

@ -9,8 +9,8 @@ module volume_calc
use geometry, only: find_cell
use global
use hdf5_interface, only: file_create, file_close, write_attribute, &
create_group, close_group, write_dataset, write_attribute_string
use output, only: write_message, header
create_group, close_group, write_dataset
use output, only: write_message, header, time_stamp
use message_passing
use particle_header, only: Particle
use random_lcg, only: prn, prn_set_stream, set_particle_seed
@ -272,11 +272,11 @@ contains
if (master) then
#ifdef MPI
do j = 1, n_procs - 1
call MPI_RECV(n, 1, MPI_INTEGER, j, 0, MPI_COMM_WORLD, &
call MPI_RECV(n, 1, MPI_INTEGER, j, 0, mpi_intracomm, &
MPI_STATUS_IGNORE, mpi_err)
allocate(data(2*n))
call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, MPI_COMM_WORLD, &
call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, mpi_intracomm, &
MPI_STATUS_IGNORE, mpi_err)
do k = 0, n - 1
do m = 1, master_indices(i_domain) % size()
@ -340,8 +340,8 @@ contains
data(2*k + 2) = master_hits(i_domain) % data(k + 1)
end do
call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, MPI_COMM_WORLD, mpi_err)
call MPI_SEND(data, 2*n, MPI_INTEGER, 0, 1, MPI_COMM_WORLD, mpi_err)
call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, mpi_intracomm, mpi_err)
call MPI_SEND(data, 2*n, MPI_INTEGER, 0, 1, mpi_intracomm, mpi_err)
deallocate(data)
#endif
end if
@ -426,14 +426,25 @@ contains
! Create HDF5 file
file_id = file_create(filename)
! Write header info
call write_attribute(file_id, "filetype", "volume")
call write_attribute(file_id, "version", VERSION_VOLUME)
call write_attribute(file_id, "openmc_version", VERSION)
#ifdef GIT_SHA1
call write_attribute(file_id, "git_sha1", GIT_SHA1)
#endif
! Write current date and time
call write_attribute(file_id, "date_and_time", time_stamp())
! Write basic metadata
select case (this % domain_type)
case (FILTER_CELL)
call write_attribute_string(file_id, ".", "domain_type", "cell")
call write_attribute(file_id, "domain_type", "cell")
case (FILTER_MATERIAL)
call write_attribute_string(file_id, ".", "domain_type", "material")
call write_attribute(file_id, "domain_type", "material")
case (FILTER_UNIVERSE)
call write_attribute_string(file_id, ".", "domain_type", "universe")
call write_attribute(file_id, "domain_type", "universe")
end select
call write_attribute(file_id, "samples", this % samples)
call write_attribute(file_id, "lower_left", this % lower_left)

9
tests/test_asymmetric_lattice/inputs_true.dat
View file

@ -205,11 +205,10 @@
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<eigenvalue>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="fission">
<parameters>-32 -32 0 32 32 32</parameters>

5
tests/test_asymmetric_lattice/test_asymmetric_lattice.py
View file

@ -89,8 +89,9 @@ class AsymmetricLatticeTestHarness(PyAPITestHarness):
outstr += '\n'.join(map('{:.8e}'.format, tally.std_dev.flatten())) + '\n'
# Extract fuel assembly lattices from the summary
core = sp.summary.get_cell_by_id(1)
fuel = sp.summary.get_cell_by_id(80)
cells = sp.summary.geometry.get_all_cells()
core = cells[1]
fuel = cells[80]
fuel = fuel.fill
core = core.fill

11
tests/test_cmfd_feed/settings.xml
View file

@ -2,11 +2,10 @@
<settings>
<!-- Parameters for criticality calculation -->
<eigenvalue>
<batches>20</batches>
<inactive>10</inactive>
<particles>1000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>20</batches>
<inactive>10</inactive>
<particles>1000</particles>
<!-- How verbose output should be -->
<verbosity value="7" />
@ -19,7 +18,7 @@
</space>
</source>
<!-- Shannon Entropy -->
<!-- Shannon Entropy -->
<entropy>
<dimension> 10 1 1 </dimension>
<lower_left> -10.0 -1.0 -1.0 </lower_left>

11
tests/test_cmfd_nofeed/settings.xml
View file

@ -2,11 +2,10 @@
<settings>
<!-- Parameters for criticality calculation -->
<eigenvalue>
<batches>20</batches>
<inactive>10</inactive>
<particles>1000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>20</batches>
<inactive>10</inactive>
<particles>1000</particles>
<!-- How verbose output should be -->
<verbosity value="7" />
@ -19,7 +18,7 @@
</space>
</source>
<!-- Shannon Entropy -->
<!-- Shannon Entropy -->
<entropy>
<dimension> 10 1 1 </dimension>
<lower_left> -10.0 -1.0 -1.0 </lower_left>

9
tests/test_complex_cell/settings.xml
View file

@ -1,11 +1,10 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<source>
<space type="box">

9
tests/test_confidence_intervals/settings.xml
View file

@ -3,11 +3,10 @@
<confidence_intervals>true</confidence_intervals>
<eigenvalue>
<batches>10</batches>
<inactive>2</inactive>
<particles>100</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>2</inactive>
<particles>100</particles>
<source>
<space type="box">

7
tests/test_create_fission_neutrons/inputs_true.dat
View file

@ -18,10 +18,9 @@
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<fixed_source>
<particles>100</particles>
<batches>10</batches>
</fixed_source>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space type="box">
<parameters>-1 -1 -1 1 1 1</parameters>

9
tests/test_density/settings.xml
View file

@ -1,11 +1,10 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<source>
<space type="box">

9
tests/test_diff_tally/inputs_true.dat
View file

@ -297,11 +297,10 @@
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<eigenvalue>
<particles>100</particles>
<batches>3</batches>
<inactive>0</inactive>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>3</batches>
<inactive>0</inactive>
<source strength="1.0">
<space type="box">
<parameters>-160 -160 -183 160 160 183</parameters>

9
tests/test_distribmat/inputs_true.dat
View file

@ -37,11 +37,10 @@
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<eigenvalue>
<particles>1000</particles>
<batches>5</batches>
<inactive>0</inactive>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>5</batches>
<inactive>0</inactive>
<source strength="1.0">
<space type="box">
<parameters>-1 -1 -1 1 1 1</parameters>

4
tests/test_distribmat/test_distribmat.py
View file

@ -116,7 +116,7 @@ class DistribmatTestHarness(PyAPITestHarness):
def _get_results(self):
outstr = super(DistribmatTestHarness, self)._get_results()
su = openmc.Summary('summary.h5')
outstr += str(su.get_cell_by_id(11))
outstr += str(su.geometry.get_all_cells()[11])
return outstr
def _cleanup(self):
@ -127,5 +127,5 @@ class DistribmatTestHarness(PyAPITestHarness):
if __name__ == '__main__':
harness = DistribmatTestHarness('statepoint.5.*')
harness = DistribmatTestHarness('statepoint.5.h5')
harness.main()

11
tests/test_eigenvalue_genperbatch/settings.xml
View file

@ -1,12 +1,11 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>7</batches>
<inactive>3</inactive>
<particles>1000</particles>
<generations_per_batch>3</generations_per_batch>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>7</batches>
<inactive>3</inactive>
<particles>1000</particles>
<generations_per_batch>3</generations_per_batch>
<source>
<space type="box">

9
tests/test_eigenvalue_no_inactive/settings.xml
View file

@ -1,11 +1,10 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>0</inactive>
<particles>1000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>0</inactive>
<particles>1000</particles>
<source>
<space type="box">

7
tests/test_energy_cutoff/inputs_true.dat
View file

@ -17,10 +17,9 @@
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<fixed_source>
<particles>100</particles>
<batches>10</batches>
</fixed_source>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space type="box">
<parameters>-1 -1 -1 1 1 1</parameters>

9
tests/test_energy_grid/settings.xml
View file

@ -3,11 +3,10 @@
<log_grid_bins>20000</log_grid_bins>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<source>
<space type="box">

Some files were not shown because too many files have changed in this diff Show more