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

Workingto integrate consistent scattering matrices with consolidation made by @nelsonag

Browse source
This commit is contained in:
Will Boyd 2017-02-27 19:49:22 -05:00
commit 0495ff7eab
71 changed files with 3445 additions and 5009 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -318,7 +318,6 @@ set(LIBOPENMC_FORTRAN_SRC
src/physics_mg.F90
src/plot.F90
src/plot_header.F90
src/ppmlib.F90
src/product_header.F90
src/progress_header.F90
src/random_lcg.F90

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

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

26
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,7 +56,7 @@ def plot_geometry(output=True, openmc_exec='openmc', cwd='.'):
def run(particles=None, threads=None, geometry_debug=False,
restart_file=None, tracks=False, output=True, cwd='.',
restart_file=None, tracks=False, output='full', cwd='.',
openmc_exec='openmc', mpi_args=None):
"""Run an OpenMC simulation.
@ -63,8 +75,10 @@ def run(particles=None, threads=None, geometry_debug=False,
Path to restart file to use
tracks : bool, optional
Write tracks for all particles. Defaults to False.
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.

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

43
openmc/material.py
View file

@ -259,6 +259,49 @@ class Material(object):
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

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):
@ -598,7 +584,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':
@ -1332,6 +1318,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()
@ -1353,11 +1340,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):

1342
openmc/mgxs/mgxs.py
View file

File diff suppressed because it is too large Load diff

35
openmc/mgxs_library.py
View file

@ -351,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)
@ -1074,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
@ -1102,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,
@ -1124,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
@ -1154,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,
@ -1307,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
@ -1334,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,
@ -1598,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])

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

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

646
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,231 +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, group in self._f['materials'].items():
if key == 'n_materials':
continue
material_id = int(key.lstrip('material '))
index = group['index'].value
name = group['name'].value.decode()
density = group['atom_density'].value
nuc_densities = group['nuclide_densities'][...]
nuclides = group['nuclides'].value
# Create the Material
material = openmc.Material(material_id=material_id, name=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()
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 = []
@ -263,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
@ -550,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

29
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

167
src/eigenvalue.F90
View file

@ -446,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
@ -458,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

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

79
src/hdf5_interface.F90
View file

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

26
src/initialize.F90
View file

@ -13,18 +13,18 @@ module initialize
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
@ -307,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.
@ -319,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
@ -333,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 &
@ -560,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 "&
@ -575,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)) &
@ -934,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

10
src/input_xml.F90
View file

@ -4261,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
@ -5412,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

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

5
src/output.F90
View file

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

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"/>

296
src/state_point.F90
View file

@ -16,7 +16,6 @@ module state_point
use hdf5
use constants
use dict_header, only: ElemKeyValueII, ElemKeyValueCI
use endf, only: reaction_name
use error, only: fatal_error, warning
use global
@ -43,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.' // &
@ -68,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)
@ -107,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)
@ -125,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)
@ -137,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
@ -145,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)
@ -228,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
@ -256,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")
@ -372,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
@ -393,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")
@ -521,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
! ==========================================================================
@ -570,92 +539,72 @@ contains
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_intracomm, 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_intracomm, 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
@ -667,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) &
@ -685,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
@ -703,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
@ -717,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")
@ -730,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
@ -745,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), &
@ -759,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
@ -800,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

154
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, &
@ -16,8 +18,6 @@ module summary
use tally_header, only: TallyObject
use tally_filter, only: find_offset
use hdf5
implicit none
private
@ -36,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)
@ -88,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
@ -120,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))
@ -152,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
@ -172,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
@ -188,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)
@ -199,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
@ -228,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
@ -244,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
@ -303,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)
@ -375,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)
@ -404,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
@ -426,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
@ -542,16 +501,11 @@ contains
call write_attribute(material_group, "depletable", 0)
end if
! Write internal OpenMC index for this material
call write_dataset(material_group, "index", i)
! 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))
@ -583,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

2
src/tally.F90
View file

@ -2027,7 +2027,7 @@ contains
end do SCORE_LOOP
nullify(matxs,nucxs)
nullify(matxs, nucxs)
end subroutine score_general_mg
!===============================================================================

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

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)

2
src/tracking.F90
View file

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

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

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

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

1
tests/test_mg_convert/test_mg_convert.py
View file

@ -153,7 +153,6 @@ class MGXSTestHarness(PyAPITestHarness):
outstr += 'k-combined:\n'
form = '{0:12.6E} {1:12.6E}\n'
outstr += form.format(sp.k_combined[0], sp.k_combined[1])
sp.close()
return outstr

98
tests/test_mg_survival_biasing/inputs_true.dat Normal file
View file

@ -0,0 +1,98 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="10000" material="10000" region="10000 -10001 10002 -10003 10004 -10005" universe="0" />
<cell id="10001" material="10001" region="10000 -10001 10002 -10003 10005 -10006" universe="0" />
<cell id="10002" material="10002" region="10000 -10001 10002 -10003 10006 -10007" universe="0" />
<cell id="10003" material="10003" region="10000 -10001 10002 -10003 10007 -10008" universe="0" />
<cell id="10004" material="10004" region="10000 -10001 10002 -10003 10008 -10009" universe="0" />
<cell id="10005" material="10005" region="10000 -10001 10002 -10003 10009 -10010" universe="0" />
<cell id="10006" material="10006" region="10000 -10001 10002 -10003 10010 -10011" universe="0" />
<cell id="10007" material="10007" region="10000 -10001 10002 -10003 10011 -10012" universe="0" />
<cell id="10008" material="10008" region="10000 -10001 10002 -10003 10012 -10013" universe="0" />
<cell id="10009" material="10009" region="10000 -10001 10002 -10003 10013 -10014" universe="0" />
<cell id="10010" material="10010" region="10000 -10001 10002 -10003 10014 -10015" universe="0" />
<cell id="10011" material="10011" region="10000 -10001 10002 -10003 10015 -10016" universe="0" />
<surface boundary="reflective" coeffs="0.0" id="10000" type="x-plane" />
<surface boundary="reflective" coeffs="10.0" id="10001" type="x-plane" />
<surface boundary="reflective" coeffs="0.0" id="10002" type="y-plane" />
<surface boundary="reflective" coeffs="10.0" id="10003" type="y-plane" />
<surface boundary="reflective" coeffs="0.0" id="10004" type="z-plane" />
<surface coeffs="0.4167" id="10005" type="z-plane" />
<surface coeffs="0.8334" id="10006" type="z-plane" />
<surface coeffs="1.2501" id="10007" type="z-plane" />
<surface coeffs="1.6668" id="10008" type="z-plane" />
<surface coeffs="2.0835" id="10009" type="z-plane" />
<surface coeffs="2.5002" id="10010" type="z-plane" />
<surface coeffs="2.9169" id="10011" type="z-plane" />
<surface coeffs="3.3336" id="10012" type="z-plane" />
<surface coeffs="3.7503" id="10013" type="z-plane" />
<surface coeffs="4.167" id="10014" type="z-plane" />
<surface coeffs="4.5837" id="10015" type="z-plane" />
<surface boundary="reflective" coeffs="5.0" id="10016" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>../1d_mgxs.h5</cross_sections>
<material id="10000" name="1">
<density units="macro" value="1.0" />
<macroscopic name="uo2_ang" />
</material>
<material id="10001" name="2">
<density units="macro" value="1.0" />
<macroscopic name="uo2_ang_mu" />
</material>
<material id="10002" name="3">
<density units="macro" value="1.0" />
<macroscopic name="uo2_iso" />
</material>
<material id="10003" name="4">
<density units="macro" value="1.0" />
<macroscopic name="uo2_iso_mu" />
</material>
<material id="10004" name="5">
<density units="macro" value="1.0" />
<macroscopic name="clad_ang" />
</material>
<material id="10005" name="6">
<density units="macro" value="1.0" />
<macroscopic name="clad_ang_mu" />
</material>
<material id="10006" name="7">
<density units="macro" value="1.0" />
<macroscopic name="clad_iso" />
</material>
<material id="10007" name="8">
<density units="macro" value="1.0" />
<macroscopic name="clad_iso_mu" />
</material>
<material id="10008" name="9">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_ang" />
</material>
<material id="10009" name="10">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_ang_mu" />
</material>
<material id="10010" name="11">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_iso" />
</material>
<material id="10011" name="12">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_iso_mu" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>0.0 0.0 0.0 10.0 10.0 5.0</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
<survival_biasing>true</survival_biasing>
</settings>

2
tests/test_mg_survival_biasing/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
1.080832E+00 1.336780E-02

21
tests/test_mg_survival_biasing/test_mg_survival_biasing.py Normal file
View file

@ -0,0 +1,21 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class MGBasicTestHarness(PyAPITestHarness):
def _build_inputs(self):
"""Write input XML files."""
self._input_set.build_default_materials_and_geometry()
self._input_set.build_default_settings()
self._input_set.settings.survival_biasing = True
self._input_set.export()
if __name__ == '__main__':
harness = MGBasicTestHarness('statepoint.10.h5', False, mg=True)
harness.main()

3
tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py
View file

@ -71,9 +71,6 @@ class MGXSTestHarness(PyAPITestHarness):
if os.path.exists('./tallies.xml'):
os.remove('./tallies.xml')
# Close the statepoint to allow writing
sp.close()
# Re-run MG mode.
if self._opts.mpi_exec is not None:
mpi_args = [self._opts.mpi_exec, '-n', self._opts.mpi_np]

2
tests/test_mgxs_library_distribcell/test_mgxs_library_distribcell.py
View file

@ -34,7 +34,7 @@ class MGXSTestHarness(PyAPITestHarness):
self.mgxs_lib.num_delayed_groups = 6
self.mgxs_lib.legendre_order = 3
self.mgxs_lib.domain_type = 'distribcell'
cells = self.mgxs_lib.openmc_geometry.get_all_material_cells()
cells = self.mgxs_lib.geometry.get_all_material_cells().values()
self.mgxs_lib.domains = [c for c in cells if c.name == 'fuel']
self.mgxs_lib.build_library()

4
tests/test_multipole/test_multipole.py
View file

@ -107,7 +107,7 @@ class MultipoleTestHarness(PyAPITestHarness):
def _get_results(self):
outstr = super(MultipoleTestHarness, 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):
@ -118,5 +118,5 @@ class MultipoleTestHarness(PyAPITestHarness):
if __name__ == '__main__':
harness = MultipoleTestHarness('statepoint.5.*')
harness = MultipoleTestHarness('statepoint.5.h5')
harness.main()

4
tests/test_particle_restart_eigval/results_true.dat
View file

@ -1,11 +1,11 @@
current batch:
1.000000E+01
current gen:
current generation:
1.000000E+00
particle id:
1.030000E+03
run mode:
k-eigenvalue
eigenvalue
particle weight:
1.000000E+00
particle energy:

2
tests/test_particle_restart_fixed/results_true.dat
View file

@ -1,6 +1,6 @@
current batch:
7.000000E+00
current gen:
current generation:
1.000000E+00
particle id:
1.440000E+02

11
tests/test_plot/settings.xml
View file

@ -1,16 +1,7 @@
<?xml version="1.0"?>
<settings>
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
<run_mode>plot</run_mode>
<entropy>
<dimension>5 4 3</dimension>

35
tests/test_plot/test_plot.py
View file

@ -33,28 +33,25 @@ class PlotTestHarness(TestHarness):
for fname in self._plot_names:
path = os.path.join(os.getcwd(), fname)
if os.path.exists(path):
#os.remove(path)
pass
os.remove(path)
def _get_results(self):
"""Return a string hash of the plot files."""
outstr = bytes()
# Add PPM output to results
ppm_files = glob.glob(os.path.join(os.getcwd(), '*.ppm'))
for fname in sorted(ppm_files):
with open(fname, 'rb') as fh:
outstr += fh.read()
# Add voxel data to results
voxel_files = glob.glob(os.path.join(os.getcwd(), '*.voxel'))
for fname in sorted(voxel_files):
with h5py.File(fname, 'r') as fh:
outstr += fh['filetype'].value
outstr += fh['num_voxels'].value.tostring()
outstr += fh['lower_left'].value.tostring()
outstr += fh['voxel_width'].value.tostring()
outstr += fh['data'].value.tostring()
for fname in self._plot_names:
if fname.endswith('.ppm'):
# Add PPM output to results
with open(fname, 'rb') as fh:
outstr += fh.read()
elif fname.endswith('.h5'):
# Add voxel data to results
with h5py.File(fname, 'r') as fh:
outstr += fh.attrs['filetype']
outstr += fh.attrs['num_voxels'].tostring()
outstr += fh.attrs['lower_left'].tostring()
outstr += fh.attrs['voxel_width'].tostring()
outstr += fh['data'].value.tostring()
# Hash the information and return.
sha512 = hashlib.sha512()
@ -65,6 +62,6 @@ class PlotTestHarness(TestHarness):
if __name__ == '__main__':
harness = PlotTestHarness(('1_plot.ppm', '2_plot.ppm', '3_plot.ppm',
'4_plot.voxel'))
harness = PlotTestHarness(('plot_1.ppm', 'plot_2.ppm', 'plot_3.ppm',
'plot_4.h5'))
harness.main()

4
tests/testing_harness.py
View file

@ -219,8 +219,8 @@ class ParticleRestartTestHarness(TestHarness):
outstr = ''
outstr += 'current batch:\n'
outstr += "{0:12.6E}\n".format(p.current_batch)
outstr += 'current gen:\n'
outstr += "{0:12.6E}\n".format(p.current_gen)
outstr += 'current generation:\n'
outstr += "{0:12.6E}\n".format(p.current_generation)
outstr += 'particle id:\n'
outstr += "{0:12.6E}\n".format(p.id)
outstr += 'run mode:\n'