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

Merge branch 'develop' into lost_particle_plotting01

Browse source 
Conflicts:
	src/DEPENDENCIES
	src/OBJECTS
	src/particle_restart.F90
	src/physics.F90
	src/source.F90
This commit is contained in:
Sterling Harper 2013-09-21 16:28:05 -04:00
commit 315ab9716d
262 changed files with 233962 additions and 3234 deletions
Download patch file Download diff file Expand all files Collapse all files

3
.gitignore vendored
View file

@ -24,3 +24,6 @@ src/xml-fortran/xmlreader
# Modules built from XML templates
src/templates/*.f90
# Test results error file
results_error.dat

5
docs/source/publications.rst
View file

@ -4,6 +4,11 @@
Publications
============
- Paul K. Romano, Andrew R. Siegel, Benoit Forget, and Kord Smith, "Data
decomposition of Monte Carlo particle transport simulations via tally
servers," *J. Comput. Phys.*, **252**, 20--36
(2013). `<http://dx.doi.org/10.1016/j.jcp.2013.06.011>`_
- Paul K. Romano, Bryan R. Herman, Nicholas E. Horelik, Benoit Forget, Kord
Smith, and Andrew R. Siegel, "Progress and Status of the OpenMC Monte Carlo
Code," *Proc. Int. Conf. Mathematics and Computational Methods Applied to

1
docs/source/quickinstall.rst
View file

@ -43,6 +43,7 @@ following commands in a terminal:
git clone git://github.com/mit-crpg/openmc.git
cd openmc/src
git checkout master
make
sudo make install

45
docs/source/releasenotes/notes_0.5.3.rst Normal file
View file

@ -0,0 +1,45 @@
.. _notes_0.5.3:
==============================
Release Notes for OpenMC 0.5.3
==============================
.. note::
These release notes are for an upcoming release of OpenMC and are still
subject to change.
-------------------
System Requirements
-------------------
There are no special requirements for running the OpenMC code. As of this
release, OpenMC has been tested on a variety of Linux distributions, Mac OS X,
and Microsoft Windows 7. Memory requirements will vary depending on the size of
the problem at hand (mostly on the number of nuclides in the problem).
------------
New Features
------------
- Special run mode --tallies removed.
- Particle restarts and state point restarts are both identified with the -r
command line flag.
- New regression test suite.
- All memory leaks fixed.
- Shared-memory parallelism with OpenMP.
---------
Bug Fixes
---------
- 2b1e8a_: Normalize direction vector after reflecting particle.
- 5853d2_: Set blank default for cross section listing alias.
- e178c7_: Fix infinite loop with words greater than 80 characters in write_message.
- c18a6e_: Chcek for valid secondary mode on S(a,b) tables.
- 82c456_: Fix bug where last process could have zero particles.
.. _2b1e8a: https://github.com/mit-crpg/openmc/commit/2b1e8a
.. _5853d2: https://github.com/mit-crpg/openmc/commit/5853d2
.. _e178c7: https://github.com/mit-crpg/openmc/commit/e178c7
.. _c18a6e: https://github.com/mit-crpg/openmc/commit/c18a6e
.. _82c456: https://github.com/mit-crpg/openmc/commit/82c456

6
docs/source/usersguide/install.rst
View file

@ -102,6 +102,12 @@ repository::
git clone git://github.com/mit-crpg/openmc.git
By default, the cloned repository will be set to the development branch. To
switch to the source of the latest stable release, run the following commands::
cd openmc/src
git checkout master
.. _GitHub: https://github.com/mit-crpg/openmc
.. _git: http://git-scm.com
.. _ssh: http://en.wikipedia.org/wiki/Secure_Shell

576
src/DEPENDENCIES
View file

@ -1,24 +1,32 @@
ace.o: ace_header.o
ace.o: constants.o
ace.o: endf.o
ace.o: error.o
ace.o: fission.o
ace.o: global.o
ace.o: material_header.o
ace.o: output.o
ace.o: set_header.o
ace.o: string.o
set_header.o: constants.o
set_header.o: list_header.o
ace_header.o: constants.o
ace_header.o: endf_header.o
energy_grid.o: constants.o
energy_grid.o: global.o
energy_grid.o: list_header.o
energy_grid.o: output.o
cmfd_data.o: cmfd_header.o
cmfd_data.o: constants.o
cmfd_data.o: error.o
cmfd_data.o: global.o
cmfd_data.o: mesh.o
cmfd_data.o: mesh_header.o
cmfd_data.o: tally_header.o
list_header.o: constants.o
cmfd_slepc_solver.o: cmfd_loss_operator.o
cmfd_slepc_solver.o: cmfd_prod_operator.o
cmfd_slepc_solver.o: constants.o
cmfd_slepc_solver.o: global.o
cmfd_loss_operator.o: constants.o
cmfd_loss_operator.o: global.o
particle_restart.o: bank_header.o
particle_restart.o: constants.o
particle_restart.o: geometry_header.o
particle_restart.o: global.o
particle_restart.o: output.o
particle_restart.o: output_interface.o
particle_restart.o: particle_header.o
particle_restart.o: random_lcg.o
particle_restart.o: tracking.o
doppler.o: constants.o
cmfd_execute.o: cmfd_data.o
cmfd_execute.o: cmfd_message_passing.o
@ -33,7 +41,75 @@ cmfd_execute.o: output.o
cmfd_execute.o: search.o
cmfd_execute.o: tally.o
cmfd_header.o: constants.o
cmfd_power_solver.o: cmfd_loss_operator.o
cmfd_power_solver.o: cmfd_prod_operator.o
cmfd_power_solver.o: constants.o
cmfd_power_solver.o: global.o
cmfd_power_solver.o: string.o
cmfd_message_passing.o: cmfd_header.o
cmfd_message_passing.o: global.o
random_lcg.o: global.o
plot.o: constants.o
plot.o: error.o
plot.o: geometry.o
plot.o: geometry_header.o
plot.o: global.o
plot.o: output.o
plot.o: particle_header.o
plot.o: plot_header.o
plot.o: ppmlib.o
plot.o: string.o
fixed_source.o: constants.o
fixed_source.o: global.o
fixed_source.o: output.o
fixed_source.o: particle_header.o
fixed_source.o: random_lcg.o
fixed_source.o: source.o
fixed_source.o: state_point.o
fixed_source.o: string.o
fixed_source.o: tally.o
fixed_source.o: tracking.o
source.o: bank_header.o
source.o: constants.o
source.o: error.o
source.o: geometry_header.o
source.o: global.o
source.o: math.o
source.o: output.o
source.o: particle_header.o
source.o: random_lcg.o
source.o: string.o
cmfd_prod_operator.o: constants.o
cmfd_prod_operator.o: global.o
ace_header.o: constants.o
ace_header.o: endf_header.o
fission.o: ace_header.o
fission.o: constants.o
fission.o: error.o
fission.o: global.o
fission.o: interpolation.o
fission.o: search.o
cmfd_jacobian_operator.o: cmfd_loss_operator.o
cmfd_jacobian_operator.o: cmfd_prod_operator.o
cmfd_jacobian_operator.o: constants.o
cmfd_jacobian_operator.o: global.o
cmfd_snes_solver.o: cmfd_jacobian_operator.o
cmfd_snes_solver.o: cmfd_loss_operator.o
cmfd_snes_solver.o: cmfd_power_solver.o
cmfd_snes_solver.o: cmfd_prod_operator.o
cmfd_snes_solver.o: constants.o
cmfd_snes_solver.o: global.o
cmfd_snes_solver.o: string.o
cmfd_input.o: cmfd_message_passing.o
cmfd_input.o: error.o
@ -46,163 +122,6 @@ cmfd_input.o: tally_header.o
cmfd_input.o: tally_initialize.o
cmfd_input.o: templates/cmfd_t.o
cmfd_jacobian_operator.o: cmfd_loss_operator.o
cmfd_jacobian_operator.o: cmfd_prod_operator.o
cmfd_jacobian_operator.o: constants.o
cmfd_jacobian_operator.o: global.o
cmfd_loss_operator.o: constants.o
cmfd_loss_operator.o: global.o
cmfd_message_passing.o: cmfd_header.o
cmfd_message_passing.o: global.o
cmfd_output.o: cmfd_data.o
cmfd_output.o: cmfd_header.o
cmfd_output.o: constants.o
cmfd_output.o: global.o
cmfd_power_solver.o: cmfd_loss_operator.o
cmfd_power_solver.o: cmfd_prod_operator.o
cmfd_power_solver.o: constants.o
cmfd_power_solver.o: global.o
cmfd_power_solver.o: string.o
cmfd_prod_operator.o: constants.o
cmfd_prod_operator.o: global.o
cmfd_slepc_solver.o: cmfd_loss_operator.o
cmfd_slepc_solver.o: cmfd_prod_operator.o
cmfd_slepc_solver.o: constants.o
cmfd_slepc_solver.o: global.o
cmfd_snes_solver.o: cmfd_jacobian_operator.o
cmfd_snes_solver.o: cmfd_loss_operator.o
cmfd_snes_solver.o: cmfd_power_solver.o
cmfd_snes_solver.o: cmfd_prod_operator.o
cmfd_snes_solver.o: constants.o
cmfd_snes_solver.o: global.o
cmfd_snes_solver.o: string.o
cross_section.o: ace_header.o
cross_section.o: constants.o
cross_section.o: error.o
cross_section.o: fission.o
cross_section.o: global.o
cross_section.o: material_header.o
cross_section.o: random_lcg.o
cross_section.o: search.o
doppler.o: constants.o
eigenvalue.o: cmfd_execute.o
eigenvalue.o: constants.o
eigenvalue.o: error.o
eigenvalue.o: global.o
eigenvalue.o: math.o
eigenvalue.o: mesh.o
eigenvalue.o: mesh_header.o
eigenvalue.o: output.o
eigenvalue.o: physics.o
eigenvalue.o: random_lcg.o
eigenvalue.o: search.o
eigenvalue.o: source.o
eigenvalue.o: state_point.o
eigenvalue.o: string.o
eigenvalue.o: tally.o
endf.o: constants.o
endf.o: string.o
energy_grid.o: constants.o
energy_grid.o: global.o
energy_grid.o: list_header.o
energy_grid.o: output.o
error.o: global.o
finalize.o: cmfd_output.o
finalize.o: global.o
finalize.o: hdf5_interface.o
finalize.o: output.o
finalize.o: tally.o
fission.o: ace_header.o
fission.o: constants.o
fission.o: error.o
fission.o: global.o
fission.o: interpolation.o
fission.o: search.o
fixed_source.o: constants.o
fixed_source.o: global.o
fixed_source.o: output.o
fixed_source.o: physics.o
fixed_source.o: random_lcg.o
fixed_source.o: source.o
fixed_source.o: state_point.o
fixed_source.o: string.o
fixed_source.o: tally.o
geometry.o: constants.o
geometry.o: error.o
geometry.o: geometry_header.o
geometry.o: global.o
geometry.o: output.o
geometry.o: particle_header.o
geometry.o: particle_restart_write.o
geometry.o: string.o
geometry.o: tally.o
global.o: ace_header.o
global.o: bank_header.o
global.o: cmfd_header.o
global.o: constants.o
global.o: dict_header.o
global.o: geometry_header.o
global.o: hdf5_interface.o
global.o: material_header.o
global.o: mesh_header.o
global.o: particle_header.o
global.o: plot_header.o
global.o: set_header.o
global.o: source_header.o
global.o: tally_header.o
global.o: timer_header.o
hdf5_summary.o: ace_header.o
hdf5_summary.o: constants.o
hdf5_summary.o: endf.o
hdf5_summary.o: geometry_header.o
hdf5_summary.o: global.o
hdf5_summary.o: material_header.o
hdf5_summary.o: mesh_header.o
hdf5_summary.o: output_interface.o
hdf5_summary.o: output.o
hdf5_summary.o: string.o
hdf5_summary.o: tally_header.o
initialize.o: ace.o
initialize.o: bank_header.o
initialize.o: constants.o
initialize.o: dict_header.o
initialize.o: energy_grid.o
initialize.o: error.o
initialize.o: geometry.o
initialize.o: geometry_header.o
initialize.o: global.o
initialize.o: hdf5_interface.o
initialize.o: hdf5_summary.o
initialize.o: input_xml.o
initialize.o: output_interface.o
initialize.o: output.o
initialize.o: random_lcg.o
initialize.o: source.o
initialize.o: state_point.o
initialize.o: string.o
initialize.o: tally_header.o
initialize.o: tally_initialize.o
input_xml.o: cmfd_input.o
input_xml.o: constants.o
input_xml.o: dict_header.o
@ -224,15 +143,6 @@ input_xml.o: templates/plots_t.o
input_xml.o: templates/settings_t.o
input_xml.o: templates/tallies_t.o
interpolation.o: constants.o
interpolation.o: endf_header.o
interpolation.o: error.o
interpolation.o: global.o
interpolation.o: search.o
interpolation.o: string.o
list_header.o: constants.o
main.o: constants.o
main.o: eigenvalue.o
main.o: finalize.o
@ -242,7 +152,88 @@ main.o: initialize.o
main.o: particle_restart.o
main.o: plot.o
particle_restart_write.o: bank_header.o
particle_restart_write.o: global.o
particle_restart_write.o: output_interface.o
particle_restart_write.o: particle_header.o
particle_restart_write.o: string.o
timer_header.o: constants.o
math.o: constants.o
math.o: random_lcg.o
interpolation.o: constants.o
interpolation.o: endf_header.o
interpolation.o: error.o
interpolation.o: global.o
interpolation.o: search.o
interpolation.o: string.o
cmfd_data.o: cmfd_header.o
cmfd_data.o: constants.o
cmfd_data.o: error.o
cmfd_data.o: global.o
cmfd_data.o: mesh.o
cmfd_data.o: mesh_header.o
cmfd_data.o: tally_header.o
cmfd_output.o: cmfd_data.o
cmfd_output.o: cmfd_header.o
cmfd_output.o: constants.o
cmfd_output.o: global.o
global.o: ace_header.o
global.o: bank_header.o
global.o: cmfd_header.o
global.o: constants.o
global.o: dict_header.o
global.o: geometry_header.o
global.o: hdf5_interface.o
global.o: material_header.o
global.o: mesh_header.o
global.o: plot_header.o
global.o: set_header.o
global.o: source_header.o
global.o: tally_header.o
global.o: timer_header.o
string.o: constants.o
string.o: error.o
string.o: global.o
finalize.o: cmfd_output.o
finalize.o: global.o
finalize.o: hdf5_interface.o
finalize.o: output.o
finalize.o: tally.o
tally_initialize.o: constants.o
tally_initialize.o: global.o
tally_initialize.o: tally_header.o
tally.o: ace_header.o
tally.o: constants.o
tally.o: error.o
tally.o: global.o
tally.o: math.o
tally.o: mesh.o
tally.o: mesh_header.o
tally.o: output.o
tally.o: particle_header.o
tally.o: search.o
tally.o: string.o
tally.o: tally_header.o
particle_header.o: constants.o
particle_header.o: geometry_header.o
output_interface.o: constants.o
output_interface.o: error.o
output_interface.o: global.o
output_interface.o: hdf5_interface.o
output_interface.o: mpiio_interface.o
output_interface.o: tally_header.o
mesh.o: constants.o
mesh.o: global.o
@ -250,6 +241,123 @@ mesh.o: mesh_header.o
mesh.o: particle_header.o
mesh.o: search.o
endf.o: constants.o
endf.o: string.o
initialize.o: ace.o
initialize.o: bank_header.o
initialize.o: constants.o
initialize.o: dict_header.o
initialize.o: energy_grid.o
initialize.o: error.o
initialize.o: geometry.o
initialize.o: geometry_header.o
initialize.o: global.o
initialize.o: hdf5_interface.o
initialize.o: hdf5_summary.o
initialize.o: input_xml.o
initialize.o: output.o
initialize.o: output_interface.o
initialize.o: random_lcg.o
initialize.o: source.o
initialize.o: state_point.o
initialize.o: string.o
initialize.o: tally_header.o
initialize.o: tally_initialize.o
cross_section.o: ace_header.o
cross_section.o: constants.o
cross_section.o: error.o
cross_section.o: fission.o
cross_section.o: global.o
cross_section.o: material_header.o
cross_section.o: particle_header.o
cross_section.o: random_lcg.o
cross_section.o: search.o
state_point.o: constants.o
state_point.o: error.o
state_point.o: global.o
state_point.o: output.o
state_point.o: output_interface.o
state_point.o: string.o
state_point.o: tally_header.o
eigenvalue.o: cmfd_execute.o
eigenvalue.o: constants.o
eigenvalue.o: error.o
eigenvalue.o: global.o
eigenvalue.o: math.o
eigenvalue.o: mesh.o
eigenvalue.o: mesh_header.o
eigenvalue.o: output.o
eigenvalue.o: particle_header.o
eigenvalue.o: random_lcg.o
eigenvalue.o: search.o
eigenvalue.o: source.o
eigenvalue.o: state_point.o
eigenvalue.o: string.o
eigenvalue.o: tally.o
eigenvalue.o: tracking.o
search.o: error.o
search.o: global.o
tracking.o: cross_section.o
tracking.o: error.o
tracking.o: geometry.o
tracking.o: geometry_header.o
tracking.o: global.o
tracking.o: output.o
tracking.o: particle_header.o
tracking.o: particle_track.o
tracking.o: physics.o
tracking.o: random_lcg.o
tracking.o: string.o
tracking.o: tally.o
ace.o: ace_header.o
ace.o: constants.o
ace.o: endf.o
ace.o: error.o
ace.o: fission.o
ace.o: global.o
ace.o: material_header.o
ace.o: output.o
ace.o: set_header.o
ace.o: string.o
geometry.o: constants.o
geometry.o: error.o
geometry.o: geometry_header.o
geometry.o: global.o
geometry.o: output.o
geometry.o: particle_header.o
geometry.o: particle_restart_write.o
geometry.o: string.o
geometry.o: tally.o
plot_header.o: constants.o
cmfd_header.o: constants.o
tally_header.o: constants.o
hdf5_summary.o: ace_header.o
hdf5_summary.o: constants.o
hdf5_summary.o: endf.o
hdf5_summary.o: geometry_header.o
hdf5_summary.o: global.o
hdf5_summary.o: hdf5_interface.o
hdf5_summary.o: material_header.o
hdf5_summary.o: mesh_header.o
hdf5_summary.o: output.o
hdf5_summary.o: output_interface.o
hdf5_summary.o: string.o
hdf5_summary.o: tally_header.o
error.o: global.o
output.o: ace_header.o
output.o: constants.o
output.o: endf.o
@ -264,119 +372,25 @@ output.o: plot_header.o
output.o: string.o
output.o: tally_header.o
output_interface.o: constants.o
output_interface.o: error.o
output_interface.o: global.o
output_interface.o: hdf5_interface.o
output_interface.o: mpiio_interface.o
output_interface.o: tally_header.o
particle_header.o: constants.o
particle_track.o: constants.o
particle_track.o: global.o
particle_track.o: output_interface.o
particle_track.o: string.o
particle_restart.o: bank_header.o
particle_restart.o: constants.o
particle_restart.o: geometry_header.o
particle_restart.o: global.o
particle_restart.o: hdf5_interface.o
particle_restart.o: output.o
particle_restart.o: particle_header.o
particle_restart.o: physics.o
particle_restart.o: random_lcg.o
particle_restart.o: source.o
particle_restart_write.o: bank_header.o
particle_restart_write.o: global.o
particle_restart_write.o: hdf5_interface.o
particle_restart_write.o: string.o
physics.o: ace_header.o
physics.o: constants.o
physics.o: cross_section.o
physics.o: endf.o
physics.o: error.o
physics.o: fission.o
physics.o: geometry.o
physics.o: geometry_header.o
physics.o: global.o
physics.o: interpolation.o
physics.o: material_header.o
physics.o: math.o
physics.o: mesh.o
physics.o: output.o
physics.o: particle_header.o
physics.o: particle_track.o
physics.o: particle_restart_write.o
physics.o: random_lcg.o
physics.o: search.o
physics.o: string.o
physics.o: tally.o
plot.o: constants.o
plot.o: error.o
plot.o: geometry.o
plot.o: geometry_header.o
plot.o: global.o
plot.o: output.o
plot.o: particle_header.o
plot.o: plot_header.o
plot.o: ppmlib.o
plot.o: source.o
plot.o: string.o
plot_header.o: constants.o
random_lcg.o: global.o
search.o: error.o
search.o: global.o
set_header.o: constants.o
set_header.o: list_header.o
source.o: bank_header.o
source.o: constants.o
source.o: error.o
source.o: geometry_header.o
source.o: global.o
source.o: output.o
source.o: particle_header.o
source.o: physics.o
source.o: random_lcg.o
source.o: string.o
state_point.o: constants.o
state_point.o: error.o
state_point.o: global.o
state_point.o: output.o
state_point.o: output_interface.o
state_point.o: string.o
state_point.o: tally_header.o
string.o: constants.o
string.o: error.o
string.o: global.o
tally.o: ace_header.o
tally.o: constants.o
tally.o: error.o
tally.o: global.o
tally.o: math.o
tally.o: mesh.o
tally.o: mesh_header.o
tally.o: output.o
tally.o: particle_header.o
tally.o: search.o
tally.o: string.o
tally.o: tally_header.o
tally_header.o: constants.o
tally_initialize.o: constants.o
tally_initialize.o: global.o
tally_initialize.o: tally_header.o
timer_header.o: constants.o

10
src/Makefile
View file

@ -1,17 +1,13 @@
program = openmc
prefix = /usr/local
source = $(wildcard *.F90)
objects = $(source:.F90=.o)
templates = $(wildcard templates/*.o)
xml_fort = xml-fortran/xmlparse.o \
xml-fortran/read_xml_primitives.o \
xml-fortran/write_xml_primitives.o
#===============================================================================
# Object Files
#===============================================================================
include OBJECTS
#===============================================================================
# User Options
#===============================================================================
@ -31,7 +27,7 @@ PETSC = no
MPI_DIR = /opt/mpich/3.0.4-$(COMPILER)
HDF5_DIR = /opt/hdf5/1.8.11-$(COMPILER)
PHDF5_DIR = /opt/phdf5/1.8.11-$(COMPILER)
PETSC_DIR = /opt/petsc/3.3-p6-$(COMPILER)
PETSC_DIR = /opt/petsc/3.4.2-$(COMPILER)
#===============================================================================
# Add git SHA-1 hash

63
src/OBJECTS
View file

@ -1,63 +0,0 @@
objects = \
ace.o \
ace_header.o \
bank_header.o \
cmfd_data.o \
cmfd_execute.o \
cmfd_header.o \
cmfd_input.o \
cmfd_jacobian_operator.o \
cmfd_loss_operator.o \
cmfd_message_passing.o \
cmfd_output.o \
cmfd_power_solver.o \
cmfd_prod_operator.o \
cmfd_snes_solver.o \
constants.o \
cross_section.o \
dict_header.o \
doppler.o \
eigenvalue.o \
endf.o \
endf_header.o \
energy_grid.o \
error.o \
finalize.o \
fission.o \
fixed_source.o \
geometry.o \
geometry_header.o \
global.o \
hdf5_interface.o \
hdf5_summary.o \
initialize.o \
interpolation.o \
input_xml.o \
list_header.o \
main.o \
material_header.o \
math.o \
mesh_header.o \
mesh.o \
mpiio_interface.o \
output.o \
output_interface.o \
particle_header.o \
particle_track.o \
particle_restart.o \
particle_restart_write.o \
physics.o \
plot.o \
plot_header.o \
ppmlib.o \
random_lcg.o \
search.o \
set_header.o \
source.o \
source_header.o \
state_point.o \
string.o \
tally.o \
tally_header.o \
tally_initialize.o \
timer_header.o

10
src/ace.F90
View file

@ -10,7 +10,7 @@ module ace
use material_header, only: Material
use output, only: write_message
use set_header, only: SetChar
use string, only: str_to_int, str_to_real, lower_case, to_str
use string, only: to_str
implicit none
@ -1178,8 +1178,14 @@ contains
integer :: NMU ! number of outgoing angles
integer :: JXS4 ! location of elastic energy table
! read secondary energy mode for inelastic scattering
! read secondary energy mode for inelastic scattering and check
table % secondary_mode = NXS(7)
if (table % secondary_mode /= SAB_SECONDARY_EQUAL .and. &
table % secondary_mode /= SAB_SECONDARY_SKEWED) then
message = "Unsupported secondary mode on S(a,b) table " // &
trim(adjustl(table % name)) // ": " // to_str(table % secondary_mode)
call fatal_error()
end if
! read number of inelastic energies and allocate arrays
NE_in = int(XSS(JXS(1)))

11
src/cmfd_execute.F90
View file

@ -281,8 +281,7 @@ contains
use constants, only: ZERO, ONE
use error, only: warning, fatal_error
use global, only: n_particles, meshes, source_bank, work, &
n_user_meshes, message, cmfd, master, mpi_err, &
bank_first, bank_last
n_user_meshes, message, cmfd, master, mpi_err
use mesh_header, only: StructuredMesh
use mesh, only: count_bank_sites, get_mesh_indices
use search, only: binary_search
@ -296,7 +295,6 @@ contains
integer :: ijk(3) ! spatial bin location
integer :: e_bin ! energy bin of source particle
integer :: n_groups ! number of energy groups
integer(8) :: size_bank ! size of source bank
logical :: outside ! any source sites outside mesh
logical :: in_mesh ! source site is inside mesh
logical :: new_weights ! calcualte new weights
@ -312,9 +310,6 @@ contains
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! compute size of source bank
size_bank = bank_last - bank_first + 1_8
! allocate arrays in cmfd object (can take out later extend to multigroup)
if (.not.allocated(cmfd%sourcecounts)) then
allocate(cmfd%sourcecounts(ng,nx,ny,nz))
@ -337,7 +332,7 @@ contains
! count bank sites in mesh
call count_bank_sites(m, source_bank, cmfd%sourcecounts, egrid, &
sites_outside=outside, size_bank = size_bank)
sites_outside=outside, size_bank=work)
! check for sites outside of the mesh
if (master .and. outside) then
@ -360,7 +355,7 @@ contains
end if
! begin loop over source bank
do i = 1, int(size_bank, 4)
do i = 1, int(work,4)
! determine spatial bin
call get_mesh_indices(m, source_bank(i)%xyz, ijk, in_mesh)

2
src/cmfd_slepc_solver.F90
View file

@ -1,6 +1,6 @@
module cmfd_slepc_solver
#ifdef PETSC
#ifdef SLEPC
use cmfd_loss_operator, only: loss_operator, init_M_operator, &
build_loss_matrix, destroy_M_operator
use cmfd_prod_operator, only: prod_operator, init_F_operator, &

3
src/constants.F90
View file

@ -256,8 +256,7 @@ module constants
EVENT_SURFACE = -2, &
EVENT_LATTICE = -1, &
EVENT_SCATTER = 1, &
EVENT_ABSORB = 2, &
EVENT_FISSION = 3
EVENT_ABSORB = 2
! Tally score type
integer, parameter :: N_SCORE_TYPES = 14

96
src/cross_section.F90
View file

@ -6,10 +6,14 @@ module cross_section
use fission, only: nu_total
use global
use material_header, only: Material
use particle_header, only: Particle
use random_lcg, only: prn
use search, only: binary_search
implicit none
save
integer :: union_grid_index
contains
@ -18,7 +22,9 @@ contains
! particle is currently traveling through.
!===============================================================================
subroutine calculate_xs()
subroutine calculate_xs(p)
type(Particle), intent(inout) :: p
integer :: i ! loop index over nuclides
integer :: i_nuclide ! index into nuclides array
@ -42,7 +48,7 @@ contains
mat => materials(p % material)
! Find energy index on unionized grid
if (grid_method == GRID_UNION) call find_energy_index()
if (grid_method == GRID_UNION) call find_energy_index(p % E)
! Determine if this material has S(a,b) tables
check_sab = (mat % n_sab > 0)
@ -84,7 +90,7 @@ contains
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E) then
call calculate_nuclide_xs(i_nuclide, i_sab)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E)
end if
! ========================================================================
@ -125,10 +131,11 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab)
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer :: i_grid ! index on nuclide energy grid
real(8) :: f ! interp factor on nuclide energy grid
@ -143,19 +150,19 @@ contains
! If we're using the unionized grid with pointers, finding the index on
! the nuclide energy grid is as simple as looking up the pointer
i_grid = nuc % grid_index(p % index_grid)
i_grid = nuc % grid_index(union_grid_index)
case (GRID_NUCLIDE)
! If we're not using the unionized grid, we have to do a binary search on
! the nuclide energy grid in order to determine which points to
! interpolate between
if (p % E < nuc % energy(1)) then
if (E < nuc % energy(1)) then
i_grid = 1
elseif (p % E > nuc % energy(nuc % n_grid)) then
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
else
i_grid = binary_search(nuc % energy, nuc % n_grid, p % E)
i_grid = binary_search(nuc % energy, nuc % n_grid, E)
end if
end select
@ -164,7 +171,7 @@ contains
if (nuc % energy(i_grid) == nuc % energy(i_grid+1)) i_grid = i_grid + 1
! calculate interpolation factor
f = (p%E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid))
f = (E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid))
micro_xs(i_nuclide) % index_grid = i_grid
micro_xs(i_nuclide) % interp_factor = f
@ -213,22 +220,22 @@ contains
! need to correct it by subtracting the non-S(a,b) elastic cross section and
! then add back in the calculated S(a,b) elastic+inelastic cross section.
if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab)
if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab, E)
! if the particle is in the unresolved resonance range and there are
! probability tables, we need to determine cross sections from the table
if (urr_ptables_on .and. nuc % urr_present) then
if (p % E > nuc % urr_data % energy(1) .and. &
p % E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then
call calculate_urr_xs(i_nuclide)
if (E > nuc % urr_data % energy(1) .and. &
E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then
call calculate_urr_xs(i_nuclide, E)
end if
end if
! Set last evaluated energy -- if we're in S(a,b) region, force
! re-calculation of cross-section
if (i_sab == 0) then
micro_xs(i_nuclide) % last_E = p % E
micro_xs(i_nuclide) % last_E = E
else
micro_xs(i_nuclide) % last_E = ZERO
end if
@ -241,10 +248,11 @@ contains
! whatever data were taken from the normal Nuclide table.
!===============================================================================
subroutine calculate_sab_xs(i_nuclide, i_sab)
subroutine calculate_sab_xs(i_nuclide, i_sab, E)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer :: i_grid ! index on S(a,b) energy grid
real(8) :: f ! interp factor on S(a,b) energy grid
@ -259,12 +267,12 @@ contains
sab => sab_tables(i_sab)
! Get index and interpolation factor for inelastic grid
if (p%E < sab % inelastic_e_in(1)) then
if (E < sab % inelastic_e_in(1)) then
i_grid = 1
f = ZERO
else
i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E)
f = (p%E - sab%inelastic_e_in(i_grid)) / &
i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
f = (E - sab%inelastic_e_in(i_grid)) / &
(sab%inelastic_e_in(i_grid+1) - sab%inelastic_e_in(i_grid))
end if
@ -273,32 +281,32 @@ contains
f * sab % inelastic_sigma(i_grid + 1)
! Check for elastic data
if (p % E < sab % threshold_elastic) then
if (E < sab % threshold_elastic) then
! Determine whether elastic scattering is given in the coherent or
! incoherent approximation. For coherent, the cross section is
! represented as P/E whereas for incoherent, it is simply P
if (sab % elastic_mode == SAB_ELASTIC_EXACT) then
if (p % E < sab % elastic_e_in(1)) then
if (E < sab % elastic_e_in(1)) then
! If energy is below that of the lowest Bragg peak, the elastic
! cross section will be zero
elastic = ZERO
else
i_grid = binary_search(sab % elastic_e_in, &
sab % n_elastic_e_in, p % E)
elastic = sab % elastic_P(i_grid) / p % E
sab % n_elastic_e_in, E)
elastic = sab % elastic_P(i_grid) / E
end if
else
! Determine index on elastic energy grid
if (p % E < sab % elastic_e_in(1)) then
if (E < sab % elastic_e_in(1)) then
i_grid = 1
else
i_grid = binary_search(sab % elastic_e_in, &
sab % n_elastic_e_in, p%E)
sab % n_elastic_e_in, E)
end if
! Get interpolation factor for elastic grid
f = (p%E - sab%elastic_e_in(i_grid))/(sab%elastic_e_in(i_grid+1) - &
f = (E - sab%elastic_e_in(i_grid))/(sab%elastic_e_in(i_grid+1) - &
sab%elastic_e_in(i_grid))
! Calculate S(a,b) elastic scattering cross section
@ -325,9 +333,10 @@ contains
! from probability tables
!===============================================================================
subroutine calculate_urr_xs(i_nuclide)
subroutine calculate_urr_xs(i_nuclide, E)
integer, intent(in) :: i_nuclide ! index into nuclides array
real(8), intent(in) :: E ! energy
integer :: i_energy ! index for energy
integer :: i_table ! index for table
@ -350,12 +359,12 @@ contains
! determine energy table
i_energy = 1
do
if (p % E < urr % energy(i_energy + 1)) exit
if (E < urr % energy(i_energy + 1)) exit
i_energy = i_energy + 1
end do
! determine interpolation factor on table
f = (p % E - urr % energy(i_energy)) / &
f = (E - urr % energy(i_energy)) / &
(urr % energy(i_energy + 1) - urr % energy(i_energy))
! sample probability table using the cumulative distribution
@ -377,7 +386,7 @@ contains
f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table)
elseif (urr % interp == LOG_LOG) then
! Get logarithmic interpolation factor
f = log(p % E / urr % energy(i_energy)) / &
f = log(E / urr % energy(i_energy)) / &
log(urr % energy(i_energy + 1) / urr % energy(i_energy))
! Calculate elastic cross section/factor
@ -440,44 +449,31 @@ contains
! Determine nu-fission cross section
if (nuc % fissionable) then
micro_xs(i_nuclide) % nu_fission = nu_total(nuc, p % E) * &
micro_xs(i_nuclide) % nu_fission = nu_total(nuc, E) * &
micro_xs(i_nuclide) % fission
end if
end subroutine calculate_urr_xs
!===============================================================================
! FIND_ENERGY_INDEX determines the index on the union energy grid and the
! interpolation factor for a particle at a certain energy
! FIND_ENERGY_INDEX determines the index on the union energy grid at a certain
! energy
!===============================================================================
subroutine find_energy_index()
subroutine find_energy_index(E)
integer :: i_grid ! index on union energy grid
real(8) :: E ! energy of particle
real(8) :: interp ! interpolation factor
! copy particle's energy
E = p % E
real(8), intent(in) :: E ! energy of particle
! if particle's energy is outside of energy grid range, set to first or last
! index. Otherwise, do a binary search through the union energy grid.
if (E < e_grid(1)) then
i_grid = 1
union_grid_index = 1
elseif (E > e_grid(n_grid)) then
i_grid = n_grid - 1
union_grid_index = n_grid - 1
else
i_grid = binary_search(e_grid, n_grid, E)
union_grid_index = binary_search(e_grid, n_grid, E)
end if
! calculate the interpolation factor -- note this will be outside of [0,1)
! for a particle outside the energy range of the union grid
interp = (E - e_grid(i_grid))/(e_grid(i_grid+1) - e_grid(i_grid))
! set particle attributes
p % index_grid = i_grid
p % interp = interp
end subroutine find_energy_index
end module cross_section

32
src/eigenvalue.F90
View file

@ -13,7 +13,7 @@ module eigenvalue
use mesh_header, only: StructuredMesh
use output, only: write_message, header, print_columns, &
print_batch_keff, print_generation
use physics, only: transport
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use source, only: get_source_particle
@ -21,6 +21,7 @@ module eigenvalue
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
use tracking, only: transport
private
public :: run_eigenvalue
@ -37,10 +38,9 @@ contains
subroutine run_eigenvalue()
if (master) call header("K EIGENVALUE SIMULATION", level=1)
type(Particle) :: p
! Allocate particle
allocate(p)
if (master) call header("K EIGENVALUE SIMULATION", level=1)
! Display column titles
call print_columns()
@ -74,10 +74,10 @@ contains
PARTICLE_LOOP: do current_work = 1, work
! grab source particle from bank
call get_source_particle(current_work)
call get_source_particle(p, current_work)
! transport particle
call transport()
call transport(p)
end do PARTICLE_LOOP
@ -98,6 +98,9 @@ contains
! END OF RUN WRAPUP
if (master) call header("SIMULATION FINISHED", level=1)
! Clear particle
call p % clear()
end subroutine run_eigenvalue
@ -375,7 +378,7 @@ contains
end if
! the last processor should not be sending sites to right
finish = bank_last
finish = work_index(rank + 1)
end if
call time_bank_sample % stop()
@ -391,11 +394,11 @@ contains
if (start < n_particles) then
! Determine the index of the processor which has the first part of the
! source_bank for the local processor
neighbor = start / maxwork
neighbor = binary_search(work_index, n_procs + 1, start) - 1
SEND_SITES: do while (start < finish)
! Determine the number of sites to send
n = min((neighbor + 1)*maxwork, finish) - start
n = min(work_index(neighbor + 1), finish) - start
! Initiate an asynchronous send of source sites to the neighboring
! process
@ -420,7 +423,7 @@ contains
! ==========================================================================
! RECEIVE BANK SITES FROM NEIGHBORS OR TEMPORARY BANK
start = bank_first - 1
start = work_index(rank)
index_local = 1
! Determine what process has the source sites that will need to be stored at
@ -432,13 +435,12 @@ contains
neighbor = binary_search(bank_position, n_procs, start) - 1
end if
RECV_SITES: do while (start < bank_last)
RECV_SITES: do while (start < work_index(rank + 1))
! Determine how many sites need to be received
if (neighbor == n_procs - 1) then
n = min(n_particles, (rank+1)*maxwork) - start
n = work_index(rank + 1) - start
else
n = min(bank_position(neighbor+2), min(n_particles, &
(rank+1)*maxwork)) - start
n = min(bank_position(neighbor + 2), work_index(rank + 1)) - start
end if
if (neighbor /= rank) then
@ -749,7 +751,7 @@ contains
else
! count number of source sites in each ufs mesh cell
call count_bank_sites(ufs_mesh, source_bank, source_frac, &
sites_outside=sites_outside)
sites_outside=sites_outside, size_bank=work)
! Check for sites outside of the mesh
if (master .and. sites_outside) then

94
src/energy_grid.F90
View file

@ -2,7 +2,7 @@ module energy_grid
use constants, only: MAX_LINE_LEN
use global
use list_header, only: ListElemReal
use list_header, only: ListReal
use output, only: write_message
contains
@ -18,9 +18,8 @@ contains
subroutine unionized_grid()
integer :: i ! index in nuclides array
type(ListElemReal), pointer :: list => null()
type(ListElemReal), pointer :: current => null()
type(Nuclide), pointer :: nuc => null()
type(ListReal), pointer :: list => null()
type(Nuclide), pointer :: nuc => null()
message = "Creating unionized energy grid..."
call write_message(5)
@ -31,24 +30,18 @@ contains
call add_grid_points(list, nuc % energy)
end do
! determine size of list
n_grid = 0
current => list
do while (associated(current))
n_grid = n_grid + 1
current => current % next
end do
! Set size of unionized energy grid
n_grid = list % size()
! create allocated array from linked list
allocate(e_grid(n_grid))
current => list
do i = 1, n_grid
e_grid(i) = current % data
current => current % next
e_grid(i) = list % get_item(i)
end do
! delete linked list and dictionary
! call list_delete(list)
call list % clear()
deallocate(list)
! Set pointers to unionized energy grid for each nuclide
call grid_pointers()
@ -62,89 +55,64 @@ contains
subroutine add_grid_points(list, energy)
type(ListElemReal), pointer :: list
type(ListReal), pointer :: list
real(8), intent(in) :: energy(:)
integer :: i ! index in energy array
integer :: n ! size of energy array
real(8) :: E ! actual energy value
type(ListElemReal), pointer :: current => null()
type(ListElemReal), pointer :: previous => null()
type(ListElemReal), pointer :: head => null()
type(ListElemReal), pointer :: tmp => null()
integer :: i ! index in energy array
integer :: n ! size of energy array
integer :: current ! current index
real(8) :: E ! actual energy value
i = 1
n = size(energy)
! if the original list is empty, we need to allocate the first element and
! If the original list is empty, we need to allocate the first element and
! store first energy point
if (.not. associated(list)) then
allocate(list)
current => list
do i = 1, n
current % data = energy(i)
if (i == n) then
current % next => null()
return
end if
allocate(current % next)
current => current % next
call list % append(energy(i))
end do
return
end if
current => list
head => list
! Set current index to beginning of the list
current = 1
do while (i <= n)
E = energy(i)
! If we've reached the end of the grid energy list, add the remaining
! energy points to the end
if (.not. associated(current)) then
! finish remaining energies
if (current > list % size()) then
! Finish remaining energies
do while (i <= n)
allocate(previous % next)
current => previous % next
current % data = energy(i)
previous => current
call list % append(energy(i))
i = i + 1
end do
current%next => null()
exit
end if
if (E < current % data) then
! create new element and insert it in energy grid list
allocate(tmp)
tmp % data = E
tmp % next => current
if (associated(previous)) then
previous % next => tmp
previous => tmp
else
previous => tmp
head => previous
end if
nullify(tmp)
if (E < list % get_item(current)) then
! advance index
! Insert new energy in this position
call list % insert(current, E)
! Advance index in linked list and in new energy grid
i = i + 1
current = current + 1
elseif (E == current % data) then
! found the exact same energy, no need to store duplicates so just
elseif (E == list % get_item(current)) then
! Found the exact same energy, no need to store duplicates so just
! skip and move to next index
i = i + 1
current = current + 1
else
previous => current
current => current % next
current = current + 1
end if
end do
! It's possible that an element was inserted at the front of the list, so we
! need to move the list pointer back to the start of the list
list => head
end subroutine add_grid_points
!===============================================================================

9
src/error.F90
View file

@ -17,7 +17,9 @@ contains
! stream.
!===============================================================================
subroutine warning()
subroutine warning(force)
logical, optional :: force ! force write from proc other than master
integer :: i_start ! starting position
integer :: i_end ! ending position
@ -26,7 +28,7 @@ contains
integer :: indent ! length of indentation
! Only allow master to print to screen
if (.not. master) return
if (.not. master .and. .not. present(force)) return
! Write warning at beginning
write(ERROR_UNIT, fmt='(1X,A)', advance='no') 'WARNING: '
@ -139,9 +141,6 @@ contains
if (current_batch > 0) then
write(ERROR_UNIT,'(1X,A,I12) ') 'Batch: ', current_batch
write(ERROR_UNIT,'(1X,A,I12) ') 'Generation:', current_gen
write(ERROR_UNIT,'(1X,A,I12)') 'Particle: ', p % id
write(ERROR_UNIT,'(1X,A,3ES12.4)') 'Location: ', p % coord0 % xyz
write(ERROR_UNIT,'(1X,A,3ES12.4)') 'Direction: ', p % coord0 % uvw
write(ERROR_UNIT,*)
end if

4
src/finalize.F90
View file

@ -69,6 +69,10 @@ contains
#endif
#ifdef MPI
! Free all MPI types
call MPI_TYPE_FREE(MPI_BANK, mpi_err)
call MPI_TYPE_FREE(MPI_TALLYRESULT, mpi_err)
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(mpi_err)
#endif

38
src/fixed_source.F90
View file

@ -1,15 +1,15 @@
module fixed_source
use constants, only: ZERO
use constants, only: ZERO
use global
use output, only: write_message, header
use physics, only: transport
use random_lcg, only: set_particle_seed
use source, only: initialize_particle, sample_external_source, &
copy_source_attributes
use state_point, only: write_state_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies
use output, only: write_message, header
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use source, only: sample_external_source, copy_source_attributes
use state_point, only: write_state_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies
use tracking, only: transport
type(Bank), pointer :: source_site => null()
@ -17,12 +17,12 @@ contains
subroutine run_fixedsource()
integer(8) :: i ! index over histories in single cycle
integer(8) :: i ! index over histories in single cycle
type(Particle) :: p
if (master) call header("FIXED SOURCE TRANSPORT SIMULATION", level=1)
! Allocate particle and dummy source site
allocate(p)
allocate(source_site)
! Turn timer and tallies on
@ -50,21 +50,21 @@ contains
PARTICLE_LOOP: do i = 1, work
! Set unique particle ID
p % id = (current_batch - 1)*n_particles + bank_first + i - 1
p % id = (current_batch - 1)*n_particles + work_index(rank) + i
! set particle trace
trace = .false.
if (current_batch == trace_batch .and. current_gen == trace_gen .and. &
bank_first + i - 1 == trace_particle) trace = .true.
work_index(rank) + i == trace_particle) trace = .true.
! set random number seed
call set_particle_seed(p % id)
! grab source particle from bank
call sample_source_particle()
call sample_source_particle(p)
! transport particle
call transport()
call transport(p)
end do PARTICLE_LOOP
@ -120,16 +120,18 @@ contains
! SAMPLE_SOURCE_PARTICLE
!===============================================================================
subroutine sample_source_particle()
subroutine sample_source_particle(p)
type(Particle), intent(inout) :: p
! Set particle
call initialize_particle()
call p % initialize()
! Sample the external source distribution
call sample_external_source(source_site)
! Copy source attributes to the particle
call copy_source_attributes(source_site)
call copy_source_attributes(p, source_site)
end subroutine sample_source_particle

96
src/geometry.F90
View file

@ -5,7 +5,7 @@ module geometry
use geometry_header, only: Cell, Surface, Universe, Lattice
use global
use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord
use particle_header, only: LocalCoord, deallocate_coord, Particle
use particle_restart_write, only: write_particle_restart
use string, only: to_str
use tally, only: score_surface_current
@ -19,10 +19,11 @@ contains
! particle are inside a cell defined as the intersection of a series of surfaces
!===============================================================================
function simple_cell_contains(c) result(in_cell)
function simple_cell_contains(c, p) result(in_cell)
type(Cell), pointer :: c
logical :: in_cell
type(Cell), pointer :: c
type(Particle), intent(inout) :: p
logical :: in_cell
integer :: i ! index of surfaces in cell
integer :: i_surface ! index in surfaces array (with sign)
@ -47,7 +48,7 @@ contains
! Determine the specified sense of the surface in the cell and the actual
! sense of the particle with respect to the surface
s => surfaces(abs(i_surface))
actual_sense = sense(s)
actual_sense = sense(p, s)
specified_sense = (c % surfaces(i) > 0)
! Compare sense of point to specified sense
@ -68,7 +69,9 @@ contains
! position using simple_cell_contains and the LocalCoord's built up by find_cell
!===============================================================================
subroutine check_cell_overlap()
subroutine check_cell_overlap(p)
type(Particle), intent(inout) :: p
integer :: i ! cell loop index on a level
integer :: n ! number of cells to search on a level
@ -92,7 +95,7 @@ contains
index_cell = univ % cells(i)
c => cells(index_cell)
if (simple_cell_contains(c)) then
if (simple_cell_contains(c, p)) then
! the particle should only be contained in one cell per level
if (index_cell /= coord % cell) then
message = "Overlapping cells detected: " // &
@ -120,10 +123,11 @@ contains
! as it's within the geometry
!===============================================================================
recursive subroutine find_cell(found, search_cells)
recursive subroutine find_cell(p, found, search_cells)
logical, intent(inout) :: found
integer, optional :: search_cells(:)
type(Particle), intent(inout) :: p
logical, intent(inout) :: found
integer, optional :: search_cells(:)
integer :: i ! index over cells
integer :: i_x, i_y, i_z ! indices in lattice
@ -166,7 +170,7 @@ contains
! get pointer to cell
c => cells(index_cell)
if (simple_cell_contains(c)) then
if (simple_cell_contains(c, p)) then
! Set cell on this level
p % coord % cell = index_cell
@ -209,7 +213,7 @@ contains
p % coord % rotated = .true.
end if
call find_cell(found)
call find_cell(p, found)
if (.not. found) exit
elseif (c % type == CELL_LATTICE) then
@ -331,7 +335,7 @@ contains
end if
if (.not. outside_lattice) then
call find_cell(found)
call find_cell(p, found)
if (.not. found) exit
end if
@ -352,9 +356,10 @@ contains
! the geometry, is reflected, or crosses into a new lattice or cell
!===============================================================================
subroutine cross_surface(last_cell)
subroutine cross_surface(p, last_cell)
integer, intent(in) :: last_cell ! last cell particle was in
type(Particle), intent(inout) :: p
integer, intent(in) :: last_cell ! last cell particle was in
real(8) :: x ! x-x0 for sphere
real(8) :: y ! y-y0 for sphere
@ -395,7 +400,7 @@ contains
! physically moving the particle forward slightly
p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
call score_surface_current()
call score_surface_current(p)
end if
! Score to global leakage tally
@ -417,7 +422,7 @@ contains
if (.not. associated(p % coord, p % coord0)) then
message = "Cannot reflect particle " // trim(to_str(p % id)) // &
" off surface in a lower universe."
call handle_lost_particle()
call handle_lost_particle(p)
return
end if
@ -427,7 +432,7 @@ contains
if (active_current_tallies % size() > 0) then
p % coord0 % xyz = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw
call score_surface_current()
call score_surface_current(p)
p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
end if
@ -568,10 +573,10 @@ contains
if (associated(p % coord0 % next)) then
call deallocate_coord(p % coord0 % next)
call find_cell(found)
call find_cell(p, found)
if (.not. found) then
message = "Couldn't find particle after reflecting from surface."
call handle_lost_particle()
call handle_lost_particle(p)
return
end if
end if
@ -594,14 +599,14 @@ contains
! If coming from negative side of surface, search all the neighboring
! cells on the positive side
call find_cell(found, surf % neighbor_pos)
call find_cell(p, found, surf % neighbor_pos)
if (found) return
elseif (p % surface < 0 .and. allocated(surf % neighbor_neg)) then
! If coming from positive side of surface, search all the neighboring
! cells on the negative side
call find_cell(found, surf % neighbor_neg)
call find_cell(p, found, surf % neighbor_neg)
if (found) return
end if
@ -613,7 +618,7 @@ contains
p % surface = NONE
p % coord => p % coord0
call deallocate_coord(p % coord % next)
call find_cell(found)
call find_cell(p, found)
if (run_mode /= MODE_PLOTTING .and. (.not. found)) then
! If a cell is still not found, there are two possible causes: 1) there is
@ -624,7 +629,7 @@ contains
p % coord => p % coord0
call deallocate_coord(p % coord % next)
p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
call find_cell(found)
call find_cell(p, found)
! Couldn't find next cell anywhere! This probably means there is an actual
! undefined region in the geometry.
@ -633,7 +638,7 @@ contains
message = "After particle " // trim(to_str(p % id)) // " crossed surface " &
// trim(to_str(surfaces(i_surface) % id)) // " it could not be &
&located in any cell and it did not leak."
call handle_lost_particle()
call handle_lost_particle(p)
return
end if
end if
@ -644,9 +649,10 @@ contains
! CROSS_LATTICE moves a particle into a new lattice element
!===============================================================================
subroutine cross_lattice(lattice_crossed)
subroutine cross_lattice(p, lattice_crossed)
integer, intent(in) :: lattice_crossed
type(Particle), intent(inout) :: p
integer, intent(in) :: lattice_crossed
integer :: i_x, i_y, i_z ! indices in lattice
integer :: n_x, n_y, n_z ! size of lattice
@ -722,11 +728,11 @@ contains
p % coord => p % coord0
! Search for particle
call find_cell(found)
call find_cell(p, found)
if (.not. found) then
message = "Could not locate particle " // trim(to_str(p % id)) // &
" after crossing a lattice boundary."
call handle_lost_particle()
call handle_lost_particle(p)
return
end if
else
@ -734,7 +740,7 @@ contains
p % coord % universe = lat % universes(i_x, i_y, i_z)
! Find cell in next lattice element
call find_cell(found)
call find_cell(p, found)
if (.not. found) then
! In some circumstances, a particle crossing the corner of a cell may not
! be able to be found in the next universe. In this scenario we cut off
@ -745,11 +751,11 @@ contains
p % coord => p % coord0
! Search for particle
call find_cell(found)
call find_cell(p, found)
if (.not. found) then
message = "Could not locate particle " // trim(to_str(p % id)) // &
" after crossing a lattice boundary."
call handle_lost_particle()
call handle_lost_particle(p)
return
end if
end if
@ -763,11 +769,12 @@ contains
! that has a parent cell, also include the surfaces of the edge of the universe.
!===============================================================================
subroutine distance_to_boundary(dist, surface_crossed, lattice_crossed)
subroutine distance_to_boundary(p, dist, surface_crossed, lattice_crossed)
real(8), intent(out) :: dist
integer, intent(out) :: surface_crossed
integer, intent(out) :: lattice_crossed
type(Particle), intent(inout) :: p
real(8), intent(out) :: dist
integer, intent(out) :: surface_crossed
integer, intent(out) :: lattice_crossed
integer :: i ! index for surface in cell
integer :: index_surf ! index in surfaces array (with sign)
@ -1342,10 +1349,11 @@ contains
! is in.
!===============================================================================
recursive function sense(surf) result(s)
recursive function sense(p, surf) result(s)
type(Surface), pointer :: surf ! surface
logical :: s ! sense of particle
type(Particle), intent(inout) :: p
type(Surface), pointer :: surf ! surface
logical :: s ! sense of particle
real(8) :: x,y,z ! coordinates of particle
real(8) :: func ! surface function evaluated at point
@ -1451,7 +1459,7 @@ contains
! Particle may be coincident with this surface. Artifically move the
! particle forward a tiny bit.
p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
s = sense(surf)
s = sense(p, surf)
elseif (func > 0) then
s = .true.
else
@ -1544,11 +1552,13 @@ contains
! HANDLE_LOST_PARTICLE
!===============================================================================
subroutine handle_lost_particle()
subroutine handle_lost_particle(p)
type(Particle), intent(inout) :: p
! Print warning and write lost particle file
call warning()
call write_particle_restart()
call warning(force = .true.)
call write_particle_restart(p)
! Increment number of lost particles
p % alive = .false.

29
src/global.F90
View file

@ -9,7 +9,6 @@ module global
use geometry_header, only: Cell, Universe, Lattice, Surface
use material_header, only: Material
use mesh_header, only: StructuredMesh
use particle_header, only: Particle
use plot_header, only: ObjectPlot
use set_header, only: SetInt
use source_header, only: ExtSource
@ -27,11 +26,6 @@ module global
implicit none
save
! ============================================================================
! THE PARTICLE
type(Particle), pointer :: p => null()
! ============================================================================
! GEOMETRY-RELATED VARIABLES
@ -166,10 +160,8 @@ module global
type(Bank), allocatable, target :: source_bank(:)
type(Bank), allocatable, target :: fission_bank(:)
integer(8) :: n_bank ! # of sites in fission bank
integer(8) :: bank_first ! index of first particle in bank
integer(8) :: bank_last ! index of last particle in bank
integer(8) :: work ! number of particles per processor
integer(8) :: maxwork ! maximum number of particles per processor
integer(8), allocatable :: work_index(:) ! starting index in source bank for each process
integer(8) :: current_work ! index in source bank of current history simulated
! Temporary k-effective values
@ -420,6 +412,19 @@ contains
if (allocated(xs_listings)) deallocate(xs_listings)
if (allocated(micro_xs)) deallocate(micro_xs)
! Deallocate external source
if (allocated(external_source % params_space)) &
deallocate(external_source % params_space)
if (allocated(external_source % params_angle)) &
deallocate(external_source % params_angle)
if (allocated(external_source % params_energy)) &
deallocate(external_source % params_energy)
! Deallocate k and entropy
if (allocated(k_generation)) deallocate(k_generation)
if (allocated(entropy)) deallocate(entropy)
if (allocated(entropy_p)) deallocate(entropy_p)
! Deallocate tally-related arrays
if (allocated(meshes)) deallocate(meshes)
if (allocated(tallies)) then
@ -440,6 +445,9 @@ contains
if (allocated(source_bank)) deallocate(source_bank)
if (allocated(entropy_p)) deallocate(entropy_p)
! Deallocate array of work indices
if (allocated(work_index)) deallocate(work_index)
! Deallocate cmfd
call deallocate_cmfd(cmfd)
@ -464,6 +472,9 @@ contains
call nuclide_dict % clear()
call sab_dict % clear()
call xs_listing_dict % clear()
! Clear statepoint batch set
call statepoint_batch % clear()
end subroutine free_memory

1642
src/hdf5_interface.F90
View file

File diff suppressed because it is too large Load diff

296
src/hdf5_summary.F90
View file

@ -14,6 +14,10 @@ module hdf5_summary
use string, only: to_str
use tally_header, only: TallyObject
implicit none
type(BinaryOutput) :: su
contains
!===============================================================================
@ -25,7 +29,7 @@ contains
character(MAX_FILE_LEN) :: filename = "summary.h5"
! Create a new file using default properties.
call file_create(filename, "serial")
call su % file_create(filename)
! Write header information
call hdf5_write_header()
@ -34,24 +38,24 @@ contains
if (run_mode == MODE_EIGENVALUE) then
! Write number of particles
call write_data(n_particles, "n_particles")
call su % write_data(n_particles, "n_particles")
! Use H5LT interface to write n_batches, n_inactive, and n_active
call write_data(n_batches, "n_batches")
call write_data(n_inactive, "n_inactive")
call write_data(n_active, "n_active")
call write_data(gen_per_batch, "gen_per_batch")
call su % write_data(n_batches, "n_batches")
call su % write_data(n_inactive, "n_inactive")
call su % write_data(n_active, "n_active")
call su % write_data(gen_per_batch, "gen_per_batch")
! Add description of each variable
call write_attribute_string("n_particles", &
call su % write_attribute_string("n_particles", &
"description", "Number of particles per generation")
call write_attribute_string("n_batches", &
call su % write_attribute_string("n_batches", &
"description", "Total number of batches")
call write_attribute_string("n_inactive", &
call su % write_attribute_string("n_inactive", &
"description", "Number of inactive batches")
call write_attribute_string("n_active", &
call su % write_attribute_string("n_active", &
"description", "Number of active batches")
call write_attribute_string("gen_per_batch", &
call su % write_attribute_string("gen_per_batch", &
"description", "Number of generations per batch")
end if
@ -63,7 +67,7 @@ contains
end if
! Terminate access to the file.
call file_close("serial")
call su % file_close()
end subroutine hdf5_write_summary
@ -74,16 +78,16 @@ contains
subroutine hdf5_write_header()
! Write version information
call write_data(VERSION_MAJOR, "version_major")
call write_data(VERSION_MINOR, "version_minor")
call write_data(VERSION_RELEASE, "version_release")
call su % write_data(VERSION_MAJOR, "version_major")
call su % write_data(VERSION_MINOR, "version_minor")
call su % write_data(VERSION_RELEASE, "version_release")
! Write current date and time
call write_data(time_stamp(), "date_and_time")
call su % write_data(time_stamp(), "date_and_time")
! Write MPI information
call write_data(n_procs, "n_procs")
call write_attribute_string("n_procs", "description", &
call su % write_data(n_procs, "n_procs")
call su % write_attribute_string("n_procs", "description", &
"Number of MPI processes")
end subroutine hdf5_write_header
@ -103,53 +107,53 @@ contains
type(Lattice), pointer :: lat => null()
! Use H5LT interface to write number of geometry objects
call write_data(n_cells, "n_cells", group="geometry")
call write_data(n_surfaces, "n_surfaces", group="geometry")
call write_data(n_universes, "n_universes", group="geometry")
call write_data(n_lattices, "n_lattices", group="geometry")
call su % write_data(n_cells, "n_cells", group="geometry")
call su % write_data(n_surfaces, "n_surfaces", group="geometry")
call su % write_data(n_universes, "n_universes", group="geometry")
call su % write_data(n_lattices, "n_lattices", group="geometry")
! ==========================================================================
! WRITE INFORMATION ON CELLS
! Create a cell group (nothing directly written in this group) then close
call hdf5_open_group("geometry/cells")
call hdf5_close_group()
call su % open_group("geometry/cells")
call su % close_group()
! Write information on each cell
CELL_LOOP: do i = 1, n_cells
c => cells(i)
! Write universe for this cell
call write_data(universes(c % universe) % id, "universe", &
call su % write_data(universes(c % universe) % id, "universe", &
group="geometry/cells/cell " // trim(to_str(c % id)))
! Write information on what fills this cell
select case (c % type)
case (CELL_NORMAL)
call write_data("normal", "fill_type", &
call su % write_data("normal", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
if (c % material == MATERIAL_VOID) then
call write_data(-1, "material", &
call su % write_data(-1, "material", &
group="geometry/cells/cell " // trim(to_str(c % id)))
else
call write_data(materials(c % material) % id, "material", &
call su % write_data(materials(c % material) % id, "material", &
group="geometry/cells/cell " // trim(to_str(c % id)))
end if
case (CELL_FILL)
call write_data("universe", "fill_type", &
call su % write_data("universe", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call write_data(universes(c % fill) % id, "material", &
call su % write_data(universes(c % fill) % id, "material", &
group="geometry/cells/cell " // trim(to_str(c % id)))
case (CELL_LATTICE)
call write_data("lattice", "fill_type", &
call su % write_data("lattice", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call write_data(lattices(c % fill) % id, "lattice", &
call su % write_data(lattices(c % fill) % id, "lattice", &
group="geometry/cells/cell " // trim(to_str(c % id)))
end select
! Write list of bounding surfaces
if (c % n_surfaces > 0) then
call write_data(c % surfaces, "surfaces", length= c % n_surfaces, &
call su % write_data(c % surfaces, "surfaces", length= c % n_surfaces, &
group="geometry/cells/cell " // trim(to_str(c % id)))
end if
@ -159,8 +163,8 @@ contains
! WRITE INFORMATION ON SURFACES
! Create surfaces group (nothing directly written here) then close
call hdf5_open_group("geometry/surfaces")
call hdf5_close_group()
call su % open_group("geometry/surfaces")
call su % close_group()
! Write information on each surface
SURFACE_LOOP: do i = 1, n_surfaces
@ -169,54 +173,54 @@ contains
! Write surface type
select case (s % type)
case (SURF_PX)
call write_data("X Plane", "type", &
call su % write_data("X Plane", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_PY)
call write_data("Y Plane", "type", &
call su % write_data("Y Plane", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_PZ)
call write_data("Z Plane", "type", &
call su % write_data("Z Plane", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_PLANE)
call write_data("Plane", "type", &
call su % write_data("Plane", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CYL_X)
call write_data("X Cylinder", "type", &
call su % write_data("X Cylinder", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CYL_Y)
call write_data("Y Cylinder", "type", &
call su % write_data("Y Cylinder", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CYL_Z)
call write_data("Z Cylinder", "type", &
call su % write_data("Z Cylinder", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_SPHERE)
call write_data("Sphere", "type", &
call su % write_data("Sphere", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CONE_X)
call write_data("X Cone", "type", &
call su % write_data("X Cone", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CONE_Y)
call write_data("Y Cone", "type", &
call su % write_data("Y Cone", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (SURF_CONE_Z)
call write_data("Z Cone", "type", &
call su % write_data("Z Cone", "type", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
end select
! Write coefficients for surface
call write_data(s % coeffs, "coefficients", length=size(s % coeffs), &
call su % write_data(s % coeffs, "coefficients", length=size(s % coeffs), &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
! Write positive neighbors
if (allocated(s % neighbor_pos)) then
call write_data(s % neighbor_pos, "neighbors_positive", &
call su % write_data(s % neighbor_pos, "neighbors_positive", &
length=size(s % neighbor_pos), &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
end if
! Write negative neighbors
if (allocated(s % neighbor_neg)) then
call write_data(s % neighbor_neg, "neighbors_negative", &
call su % write_data(s % neighbor_neg, "neighbors_negative", &
length=size(s % neighbor_neg), &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
end if
@ -224,16 +228,16 @@ contains
! Write boundary condition
select case (s % bc)
case (BC_TRANSMIT)
call write_data("transmission", "boundary_condition", &
call su % write_data("transmission", "boundary_condition", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (BC_VACUUM)
call write_data("vacuum", "boundary_condition", &
call su % write_data("vacuum", "boundary_condition", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (BC_REFLECT)
call write_data("reflective", "boundary_condition", &
call su % write_data("reflective", "boundary_condition", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
case (BC_PERIODIC)
call write_data("periodic", "boundary_condition", &
call su % write_data("periodic", "boundary_condition", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
end select
@ -243,8 +247,8 @@ contains
! WRITE INFORMATION ON UNIVERSES
! Create universes group (nothing directly written here) then close
call hdf5_open_group("geometry/universes")
call hdf5_close_group()
call su % open_group("geometry/universes")
call su % close_group()
! Write information on each universe
UNIVERSE_LOOP: do i = 1, n_universes
@ -252,7 +256,7 @@ contains
! Write list of cells in this universe
if (u % n_cells > 0) then
call write_data(u % cells, "cells", length=u % n_cells, &
call su % write_data(u % cells, "cells", length=u % n_cells, &
group="geometry/universes/universe " // trim(to_str(u % id)))
end if
@ -262,8 +266,8 @@ contains
! WRITE INFORMATION ON LATTICES
! Create lattices group (nothing directly written here) then close
call hdf5_open_group("geometry/lattices")
call hdf5_close_group()
call su % open_group("geometry/lattices")
call su % close_group()
! Write information on each lattice
LATTICE_LOOP: do i = 1, n_lattices
@ -272,21 +276,21 @@ contains
! Write lattice type
select case(lat % type)
case (LATTICE_RECT)
call write_data("rectangular", "type", &
call su % write_data("rectangular", "type", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
case (LATTICE_HEX)
call write_data("hexagonal", "type", &
call su % write_data("hexagonal", "type", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end select
! Write lattice dimensions, lower left corner, and width of element
call write_data(lat % dimension, "dimension", &
call su % write_data(lat % dimension, "dimension", &
length=lat % n_dimension, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
call write_data(lat % lower_left, "lower_left", &
call su % write_data(lat % lower_left, "lower_left", &
length=lat % n_dimension, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
call write_data(lat % width, "width", &
call su % write_data(lat % width, "width", &
length=lat % n_dimension, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
@ -308,7 +312,7 @@ contains
end do
end do
end do
call write_data(lattice_universes, "universes", &
call su % write_data(lattice_universes, "universes", &
length=(/n_x, n_y, n_z/), &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
deallocate(lattice_universes)
@ -329,16 +333,16 @@ contains
type(Material), pointer :: m => null()
! Use H5LT interface to write number of materials
call write_data(n_materials, "n_materials", group="materials")
call su % write_data(n_materials, "n_materials", group="materials")
! Write information on each material
do i = 1, n_materials
m => materials(i)
! Write atom density with units
call write_data(m % density, "atom_density", &
call su % write_data(m % density, "atom_density", &
group="materials/material " // trim(to_str(m % id)))
call write_attribute_string("atom_density", "units", "atom/b-cm", &
call su % write_attribute_string("atom_density", "units", "atom/b-cm", &
group="materials/material " // trim(to_str(m % id)))
! Copy ZAID for each nuclide to temporary array
@ -348,23 +352,23 @@ contains
end do
! Write temporary array to 'nuclides'
call write_data(zaids, "nuclides", length=m % n_nuclides, &
call su % write_data(zaids, "nuclides", length=m % n_nuclides, &
group="materials/material " // trim(to_str(m % id)))
! Deallocate temporary array
deallocate(zaids)
! Write atom densities
call write_data(m % atom_density, "nuclide_densities", &
call su % write_data(m % atom_density, "nuclide_densities", &
length=m % n_nuclides, &
group="materials/material " // trim(to_str(m % id)))
! Write S(a,b) information if present
if (m % n_sab > 0) then
call write_data(m % i_sab_nuclides, "i_sab_nuclides", &
call su % write_data(m % i_sab_nuclides, "i_sab_nuclides", &
length=m % n_sab, &
group="materials/material " // trim(to_str(m % id)))
call write_data(m % i_sab_tables, "i_sab_tables", &
call su % write_data(m % i_sab_tables, "i_sab_tables", &
length=m % n_sab, &
group="materials/material " // trim(to_str(m % id)))
end if
@ -385,72 +389,72 @@ contains
type(TallyObject), pointer :: t => null()
! Write total number of meshes
call write_data(n_meshes, "n_meshes", group="tallies")
call su % write_data(n_meshes, "n_meshes", group="tallies")
! Write information for meshes
MESH_LOOP: do i = 1, n_meshes
m => meshes(i)
! Write type and number of dimensions
call write_data(m % type, "type", &
call su % write_data(m % type, "type", &
group="tallies/mesh " // trim(to_str(m % id)))
call write_data(m % n_dimension, "n_dimension", &
call su % write_data(m % n_dimension, "n_dimension", &
group="tallies/mesh " // trim(to_str(m % id)))
! Write mesh information
call write_data(m % dimension, "dimension", &
call su % write_data(m % dimension, "dimension", &
length=m % n_dimension, &
group="tallies/mesh " // trim(to_str(m % id)))
call write_data(m % lower_left, "lower_left", &
call su % write_data(m % lower_left, "lower_left", &
length=m % n_dimension, &
group="tallies/mesh " // trim(to_str(m % id)))
call write_data(m % upper_right, "upper_right", &
call su % write_data(m % upper_right, "upper_right", &
length=m % n_dimension, &
group="tallies/mesh " // trim(to_str(m % id)))
call write_data(m % width, "width", &
call su % write_data(m % width, "width", &
length=m % n_dimension, &
group="tallies/mesh " // trim(to_str(m % id)))
end do MESH_LOOP
! Write number of tallies
call write_data(n_tallies, "n_tallies", group="tallies")
call su % write_data(n_tallies, "n_tallies", group="tallies")
TALLY_METADATA: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
! Write size of each tally
call write_data(t % total_score_bins, "total_score_bins", &
call su % write_data(t % total_score_bins, "total_score_bins", &
group="tallies/tally " // trim(to_str(t % id)))
call write_data(t % total_filter_bins, "total_filter_bins", &
call su % write_data(t % total_filter_bins, "total_filter_bins", &
group="tallies/tally " // trim(to_str(t % id)))
! Write number of filters
call write_data(t % n_filters, "n_filters", &
call su % write_data(t % n_filters, "n_filters", &
group="tallies/tally " // trim(to_str(t % id)))
FILTER_LOOP: do j = 1, t % n_filters
! Write type of filter
call write_data(t % filters(j) % type, "type", &
call su % write_data(t % filters(j) % type, "type", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
! Write number of bins for this filter
call write_data(t % filters(j) % n_bins, "n_bins", &
call su % write_data(t % filters(j) % n_bins, "n_bins", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
! Write filter bins
if (t % filters(j) % type == FILTER_ENERGYIN .or. &
t % filters(j) % type == FILTER_ENERGYOUT) then
call write_data(t % filters(j) % real_bins, "bins", &
call su % write_data(t % filters(j) % real_bins, "bins", &
length=size(t % filters(j) % real_bins), &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
else
call write_data(t % filters(j) % int_bins, "bins", &
call su % write_data(t % filters(j) % int_bins, "bins", &
length=size(t % filters(j) % int_bins), &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
@ -459,35 +463,35 @@ contains
! Write name of type
select case (t % filters(j) % type)
case(FILTER_UNIVERSE)
call write_data("universe", "type_name", &
call su % write_data("universe", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_MATERIAL)
call write_data("material", "type_name", &
call su % write_data("material", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_CELL)
call write_data("cell", "type_name", &
call su % write_data("cell", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_CELLBORN)
call write_data("cellborn", "type_name", &
call su % write_data("cellborn", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_SURFACE)
call write_data("surface", "type_name", &
call su % write_data("surface", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_MESH)
call write_data("mesh", "type_name", &
call su % write_data("mesh", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_ENERGYIN)
call write_data("energy", "type_name", &
call su % write_data("energy", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_ENERGYOUT)
call write_data("energyout", "type_name", &
call su % write_data("energyout", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
end select
@ -495,7 +499,7 @@ contains
end do FILTER_LOOP
! Write number of nuclide bins
call write_data(t % n_nuclide_bins, "n_nuclide_bins", &
call su % write_data(t % n_nuclide_bins, "n_nuclide_bins", &
group="tallies/tally " // trim(to_str(t % id)))
! Create temporary array for nuclide bins
@ -509,14 +513,14 @@ contains
end do NUCLIDE_LOOP
! Write and deallocate nuclide bins
call write_data(temp_array, "nuclide_bins", length=t % n_nuclide_bins, &
call su % write_data(temp_array, "nuclide_bins", length=t % n_nuclide_bins, &
group="tallies/tally " // trim(to_str(t % id)))
deallocate(temp_array)
! Write number of score bins
call write_data(t % n_score_bins, "n_score_bins", &
call su % write_data(t % n_score_bins, "n_score_bins", &
group="tallies/tally " // trim(to_str(t % id)))
call write_data(t % score_bins, "score_bins", length=t % n_score_bins, &
call su % write_data(t % score_bins, "score_bins", length=t % n_score_bins, &
group="tallies/tally " // trim(to_str(t % id)))
end do TALLY_METADATA
@ -539,7 +543,7 @@ contains
type(UrrData), pointer :: urr => null()
! Use H5LT interface to write number of nuclides
call write_data(n_nuclides_total, "n_nuclides", group="nuclides")
call su % write_data(n_nuclides_total, "n_nuclides", group="nuclides")
! Write information on each nuclide
NUCLIDE_LOOP: do i = 1, n_nuclides_total
@ -550,27 +554,27 @@ contains
size_total = size_xs
! Write some basic attributes
call write_data(nuc % zaid, "zaid", &
call su % write_data(nuc % zaid, "zaid", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % awr, "awr", &
call su % write_data(nuc % awr, "awr", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % kT, "kT", &
call su % write_data(nuc % kT, "kT", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % n_grid, "n_grid", &
call su % write_data(nuc % n_grid, "n_grid", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % n_reaction, "n_reactions", &
call su % write_data(nuc % n_reaction, "n_reactions", &
group="nuclides/" // trim(nuc % name))
call write_data(nuc % n_fission, "n_fission", &
call su % write_data(nuc % n_fission, "n_fission", &
group="nuclides/" // trim(nuc % name))
call write_data(size_xs, "size_xs", &
call su % write_data(size_xs, "size_xs", &
group="nuclides/" // trim(nuc % name))
! =======================================================================
! WRITE INFORMATION ON EACH REACTION
! Create overall group for reactions and close it
call hdf5_open_group("nuclides/" // trim(nuc % name) // "/reactions")
call hdf5_close_group()
call su % open_group("nuclides/" // trim(nuc % name) // "/reactions")
call su % close_group()
RXN_LOOP: do j = 1, nuc % n_reaction
! Information on each reaction
@ -591,19 +595,19 @@ contains
end if
! Write information on reaction
call write_data(rxn % Q_value, "Q_value", &
call su % write_data(rxn % Q_value, "Q_value", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
call write_data(rxn % multiplicity, "multiplicity", &
call su % write_data(rxn % multiplicity, "multiplicity", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
call write_data(rxn % threshold, "threshold", &
call su % write_data(rxn % threshold, "threshold", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
call write_data(size_angle, "size_angle", &
call su % write_data(size_angle, "size_angle", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
call write_data(size_energy, "size_energy", &
call su % write_data(size_energy, "size_energy", &
group="nuclides/" // trim(nuc % name) // "/reactions/" // &
trim(reaction_name(rxn % MT)))
@ -616,24 +620,24 @@ contains
if (nuc % urr_present) then
urr => nuc % urr_data
call write_data(urr % n_energy, "urr_n_energy", &
call su % write_data(urr % n_energy, "urr_n_energy", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % n_prob, "urr_n_prob", &
call su % write_data(urr % n_prob, "urr_n_prob", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % interp, "urr_interp", &
call su % write_data(urr % interp, "urr_interp", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % inelastic_flag, "urr_inelastic", &
call su % write_data(urr % inelastic_flag, "urr_inelastic", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % absorption_flag, "urr_absorption", &
call su % write_data(urr % absorption_flag, "urr_absorption", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % energy(1), "urr_min_E", &
call su % write_data(urr % energy(1), "urr_min_E", &
group="nuclides/" // trim(nuc % name))
call write_data(urr % energy(urr % n_energy), "urr_max_E", &
call su % write_data(urr % energy(urr % n_energy), "urr_max_E", &
group="nuclides/" // trim(nuc % name))
end if
! Write total memory used
call write_data(size_total, "size_total", &
call su % write_data(size_total, "size_total", &
group="nuclides/" // trim(nuc % name))
end do NUCLIDE_LOOP
@ -650,63 +654,63 @@ contains
real(8) :: speed
! Write timing data
call write_data(time_initialize % elapsed, "time_initialize", &
call su % write_data(time_initialize % elapsed, "time_initialize", &
group="timing")
call write_data(time_read_xs % elapsed, "time_read_xs", &
call su % write_data(time_read_xs % elapsed, "time_read_xs", &
group="timing")
call write_data(time_unionize % elapsed, "time_unionize", &
call su % write_data(time_unionize % elapsed, "time_unionize", &
group="timing")
call write_data(time_transport % elapsed, "time_transport", &
call su % write_data(time_transport % elapsed, "time_transport", &
group="timing")
call write_data(time_bank % elapsed, "time_bank", &
call su % write_data(time_bank % elapsed, "time_bank", &
group="timing")
call write_data(time_bank_sample % elapsed, "time_bank_sample", &
call su % write_data(time_bank_sample % elapsed, "time_bank_sample", &
group="timing")
call write_data(time_bank_sendrecv % elapsed, "time_bank_sendrecv", &
call su % write_data(time_bank_sendrecv % elapsed, "time_bank_sendrecv", &
group="timing")
call write_data(time_tallies % elapsed, "time_tallies", &
call su % write_data(time_tallies % elapsed, "time_tallies", &
group="timing")
call write_data(time_inactive % elapsed, "time_inactive", &
call su % write_data(time_inactive % elapsed, "time_inactive", &
group="timing")
call write_data(time_active % elapsed, "time_active", &
call su % write_data(time_active % elapsed, "time_active", &
group="timing")
call write_data(time_finalize % elapsed, "time_finalize", &
call su % write_data(time_finalize % elapsed, "time_finalize", &
group="timing")
call write_data(time_total % elapsed, "time_total", &
call su % write_data(time_total % elapsed, "time_total", &
group="timing")
! Add descriptions to timing data
call write_attribute_string("time_initialize", "description", &
call su % write_attribute_string("time_initialize", "description", &
"Total time elapsed for initialization (s)", group="timing")
call write_attribute_string("time_read_xs", "description", &
call su % write_attribute_string("time_read_xs", "description", &
"Time reading cross-section libraries (s)", group="timing")
call write_attribute_string("time_unionize", "description", &
call su % write_attribute_string("time_unionize", "description", &
"Time unionizing energy grid (s)", group="timing")
call write_attribute_string("time_transport", "description", &
call su % write_attribute_string("time_transport", "description", &
"Time in transport only (s)", group="timing")
call write_attribute_string("time_bank", "description", &
call su % write_attribute_string("time_bank", "description", &
"Total time synchronizing fission bank (s)", group="timing")
call write_attribute_string("time_bank_sample", "description", &
call su % write_attribute_string("time_bank_sample", "description", &
"Time between generations sampling source sites (s)", group="timing")
call write_attribute_string("time_bank_sendrecv", "description", &
call su % write_attribute_string("time_bank_sendrecv", "description", &
"Time between generations SEND/RECVing source sites (s)", &
group="timing")
call write_attribute_string("time_tallies", "description", &
call su % write_attribute_string("time_tallies", "description", &
"Time between batches accumulating tallies (s)", group="timing")
call write_attribute_string("time_inactive", "description", &
call su % write_attribute_string("time_inactive", "description", &
"Total time in inactive batches (s)", group="timing")
call write_attribute_string("time_active", "description", &
call su % write_attribute_string("time_active", "description", &
"Total time in active batches (s)", group="timing")
call write_attribute_string("time_finalize", "description", &
call su % write_attribute_string("time_finalize", "description", &
"Total time for finalization (s)", group="timing")
call write_attribute_string("time_total", "description", &
call su % write_attribute_string("time_total", "description", &
"Total time elapsed (s)", group="timing")
! Write calculation rate
total_particles = n_particles * n_batches * gen_per_batch
speed = real(total_particles) / (time_inactive % elapsed + &
time_active % elapsed)
call write_data(speed, "neutrons_per_second", group="timing")
call su % write_data(speed, "neutrons_per_second", group="timing")
end subroutine hdf5_write_timing

84
src/initialize.F90
View file

@ -14,7 +14,7 @@ module initialize
use output, only: title, header, write_summary, print_version, &
print_usage, write_xs_summary, print_plot, &
write_message
use output_interface, only: file_open, file_close, read_data
use output_interface
use random_lcg, only: initialize_prng
use source, only: initialize_source
use state_point, only: load_state_point
@ -242,6 +242,9 @@ contains
! Commit derived type for tally scores
call MPI_TYPE_COMMIT(MPI_TALLYRESULT, mpi_err)
! Free temporary MPI type
call MPI_TYPE_FREE(temp_type, mpi_err)
end subroutine initialize_mpi
#endif
@ -303,6 +306,7 @@ contains
integer :: filetype
character(MAX_FILE_LEN) :: pwd ! present working directory
character(MAX_WORD_LEN), allocatable :: argv(:) ! command line arguments
type(BinaryOutput) :: sp
! Get working directory
call GET_ENVIRONMENT_VARIABLE("PWD", pwd)
@ -343,9 +347,9 @@ contains
i = i + 1
! Check what type of file this is
call file_open(argv(i), 'parallel', 'r')
call read_data(filetype, 'filetype')
call file_close('parallel')
call sp % file_open(argv(i), 'r', serial = .false.)
call sp % read_data(filetype, 'filetype')
call sp % file_close()
! Set path and flag for type of run
select case (filetype)
@ -420,6 +424,8 @@ contains
integer, allocatable :: index_cell_in_univ(:) ! the index in the univ%cells
! array for each universe
type(ElemKeyValueII), pointer :: pair_list => null()
type(ElemKeyValueII), pointer :: current => null()
type(ElemKeyValueII), pointer :: next => null()
type(Universe), pointer :: univ => null()
type(Cell), pointer :: c => null()
@ -431,24 +437,27 @@ contains
! cells_in_univ_dict, it's the id of the universe and the number of cells.
pair_list => universe_dict % keys()
do while (associated(pair_list))
! find index of universe in universes array
i_univ = pair_list % value
current => pair_list
do while (associated(current))
! Find index of universe in universes array
i_univ = current % value
univ => universes(i_univ)
univ % id = pair_list % key
univ % id = current % key
! check for lowest level universe
! Check for lowest level universe
if (univ % id == 0) BASE_UNIVERSE = i_univ
! find cell count for this universe
! Find cell count for this universe
n_cells_in_univ = cells_in_univ_dict % get_key(univ % id)
! allocate cell list for universe
! Allocate cell list for universe
allocate(univ % cells(n_cells_in_univ))
univ % n_cells = n_cells_in_univ
! move to next universe
pair_list => pair_list % next
! Move to next universe
next => current % next
deallocate(current)
current => next
end do
! Also allocate a list for keeping track of where cells have been assigned
@ -460,16 +469,19 @@ contains
do i = 1, n_cells
c => cells(i)
! get pointer to corresponding universe
! Get pointer to corresponding universe
i_univ = universe_dict % get_key(c % universe)
univ => universes(i_univ)
! increment the index for the cells array within the Universe object and
! Increment the index for the cells array within the Universe object and
! then store the index of the Cell object in that array
index_cell_in_univ(i_univ) = index_cell_in_univ(i_univ) + 1
univ % cells(index_cell_in_univ(i_univ)) = i
end do
! Clear dictionary
call cells_in_univ_dict % clear()
end subroutine prepare_universes
!===============================================================================
@ -746,15 +758,37 @@ contains
subroutine calculate_work()
! Determine maximum amount of particles to simulate on each processor
maxwork = ceiling(real(n_particles)/n_procs,8)
integer :: i ! loop index
integer :: remainder ! Number of processors with one extra particle
integer(8) :: i_bank ! Running count of number of particles
integer(8) :: min_work ! Minimum number of particles on each proc
integer(8) :: work_i ! Number of particles on rank i
! ID's of first and last source particles
bank_first = rank*maxwork + 1
bank_last = min((rank+1)*maxwork, n_particles)
allocate(work_index(0:n_procs))
! number of particles for this processor
work = bank_last - bank_first + 1
! Determine minimum amount of particles to simulate on each processor
min_work = n_particles/n_procs
! Determine number of processors that have one extra particle
remainder = int(mod(n_particles, int(n_procs,8)), 4)
i_bank = 0
work_index(0) = 0
do i = 0, n_procs - 1
! Number of particles for rank i
if (i < remainder) then
work_i = min_work + 1
else
work_i = min_work
end if
! Set number of particles
if (rank == i) work = work_i
! Set index into source bank for rank i
i_bank = i_bank + work_i
work_index(i+1) = i_bank
end do
end subroutine calculate_work
@ -767,7 +801,7 @@ contains
integer :: alloc_err ! allocation error code
! Allocate source bank
allocate(source_bank(maxwork), STAT=alloc_err)
allocate(source_bank(work), STAT=alloc_err)
! Check for allocation errors
if (alloc_err /= 0) then
@ -776,7 +810,7 @@ contains
end if
! Allocate fission bank
allocate(fission_bank(3*maxwork), STAT=alloc_err)
allocate(fission_bank(3*work), STAT=alloc_err)
! Check for allocation errors
if (alloc_err /= 0) then

93
src/list_header.F90
View file

@ -19,16 +19,19 @@ module list_header
type :: ListElemInt
integer :: data
type(ListElemInt), pointer :: next => null()
type(ListElemInt), pointer :: prev => null()
end type ListElemInt
type :: ListElemReal
real(8) :: data
type(ListElemReal), pointer :: next => null()
type(ListElemReal), pointer :: prev => null()
end type ListElemReal
type :: ListElemChar
character(MAX_WORD_LEN) :: data
type(ListElemChar), pointer :: next => null()
type(ListElemChar), pointer :: prev => null()
end type ListElemChar
!===============================================================================
@ -128,6 +131,7 @@ contains
else
! Otherwise append element at end of list
this % tail % next => elem
elem % prev => this % tail
this % tail => this % tail % next
end if
@ -152,6 +156,7 @@ contains
else
! Otherwise append element at end of list
this % tail % next => elem
elem % prev => this % tail
this % tail => this % tail % next
end if
@ -176,6 +181,7 @@ contains
else
! Otherwise append element at end of list
this % tail % next => elem
elem % prev => this % tail
this % tail => this % tail % next
end if
@ -499,8 +505,13 @@ contains
else
! Default case with new element somewhere in middle of list
i = 0
elem => this % head
if (i_list >= this % last_index) then
i = this % last_index
elem => this % last_elem
else
i = 0
elem => this % head
end if
do while (associated(elem))
i = i + 1
if (i == i_list - 1) then
@ -509,11 +520,17 @@ contains
new_elem % data = data
! Put it before the i-th element
new_elem % next => elem % next
elem % next => new_elem
new_elem % prev => elem % prev
new_elem % next => elem
new_elem % prev % next => new_elem
new_elem % next % prev => new_elem
this % count = this % count + 1
this % last_index = i_list
this % last_elem => new_elem
exit
end if
i = i + 1
elem => elem % next
end do
end if
@ -544,21 +561,31 @@ contains
else
! Default case with new element somewhere in middle of list
i = 0
elem => this % head
if (i_list >= this % last_index) then
i = this % last_index
elem => this % last_elem
else
i = 0
elem => this % head
end if
do while (associated(elem))
i = i + 1
if (i == i_list - 1) then
if (i == i_list) then
! Allocate new element
allocate(new_elem)
new_elem % data = data
! Put it before the i-th element
new_elem % next => elem % next
elem % next => new_elem
new_elem % prev => elem % prev
new_elem % next => elem
new_elem % prev % next => new_elem
new_elem % next % prev => new_elem
this % count = this % count + 1
this % last_index = i_list
this % last_elem => new_elem
exit
end if
i = i + 1
elem => elem % next
end do
end if
@ -589,21 +616,31 @@ contains
else
! Default case with new element somewhere in middle of list
i = 0
elem => this % head
if (i_list >= this % last_index) then
i = this % last_index
elem => this % last_elem
else
i = 0
elem => this % head
end if
do while (associated(elem))
i = i + 1
if (i == i_list - 1) then
if (i == i_list) then
! Allocate new element
allocate(new_elem)
new_elem % data = data
! Put it before the i-th element
new_elem % next => elem % next
elem % next => new_elem
new_elem % prev => elem % prev
new_elem % next => elem
new_elem % prev % next => new_elem
new_elem % next % prev => new_elem
this % count = this % count + 1
this % last_index = i_list
this % last_elem => new_elem
exit
end if
i = i + 1
elem => elem % next
end do
end if
@ -620,7 +657,6 @@ contains
integer :: data
type(ListElemInt), pointer :: elem => null()
type(ListElemInt), pointer :: prev => null()
elem => this % head
do while (associated(elem))
@ -632,12 +668,14 @@ contains
if (associated(elem, this % head)) then
this % head => elem % next
if (associated(elem, this % tail)) nullify(this % tail)
if (associated(this % head)) nullify(this % head % prev)
deallocate(elem)
else if (associated(elem, this % tail)) then
this % tail => prev
this % tail => elem % prev
deallocate(this % tail % next)
else
prev % next => elem % next
elem % prev % next => elem % next
elem % next % prev => elem % prev
deallocate(elem)
end if
@ -647,7 +685,6 @@ contains
end if
! Advance pointers
prev => elem
elem => elem % next
end do
@ -659,7 +696,6 @@ contains
real(8) :: data
type(ListElemReal), pointer :: elem => null()
type(ListElemReal), pointer :: prev => null()
elem => this % head
do while (associated(elem))
@ -671,12 +707,14 @@ contains
if (associated(elem, this % head)) then
this % head => elem % next
if (associated(elem, this % tail)) nullify(this % tail)
if (associated(this % head)) nullify(this % head % prev)
deallocate(elem)
else if (associated(elem, this % tail)) then
this % tail => prev
this % tail => elem % prev
deallocate(this % tail % next)
else
prev % next => elem % next
elem % prev % next => elem % next
elem % next % prev => elem % prev
deallocate(elem)
end if
@ -686,7 +724,6 @@ contains
end if
! Advance pointers
prev => elem
elem => elem % next
end do
@ -698,7 +735,6 @@ contains
character(*) :: data
type(ListElemChar), pointer :: elem => null()
type(ListElemChar), pointer :: prev => null()
elem => this % head
do while (associated(elem))
@ -710,12 +746,14 @@ contains
if (associated(elem, this % head)) then
this % head => elem % next
if (associated(elem, this % tail)) nullify(this % tail)
if (associated(this % head)) nullify(this % head % prev)
deallocate(elem)
else if (associated(elem, this % tail)) then
this % tail => prev
this % tail => elem % prev
deallocate(this % tail % next)
else
prev % next => elem % next
elem % prev % next => elem % next
elem % next % prev => elem % prev
deallocate(elem)
end if
@ -725,7 +763,6 @@ contains
end if
! Advance pointers
prev => elem
elem => elem % next
end do

51
src/math.F90
View file

@ -1,6 +1,7 @@
module math
use constants, only: PI, ONE, TWO
use constants, only: PI, ONE, TWO
use random_lcg, only: prn
implicit none
@ -162,4 +163,52 @@ contains
end function calc_pn
!===============================================================================
! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based
! on a direct sampling scheme. The probability distribution function for a
! Maxwellian is given as p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can
! be sampled using rule C64 in the Monte Carlo Sampler LA-9721-MS.
!===============================================================================
function maxwell_spectrum(T) result(E_out)
real(8), intent(in) :: T ! tabulated function of incoming E
real(8) :: E_out ! sampled energy
real(8) :: r1, r2, r3 ! random numbers
real(8) :: c ! cosine of pi/2*r3
r1 = prn()
r2 = prn()
r3 = prn()
! determine cosine of pi/2*r
c = cos(PI/TWO*r3)
! determine outgoing energy
E_out = -T*(log(r1) + log(r2)*c*c)
end function maxwell_spectrum
!===============================================================================
! WATT_SPECTRUM samples the outgoing energy from a Watt energy-dependent fission
! spectrum. Although fitted parameters exist for many nuclides, generally the
! continuous tabular distributions (LAW 4) should be used in lieu of the Watt
! spectrum. This direct sampling scheme is an unpublished scheme based on the
! original Watt spectrum derivation (See F. Brown's MC lectures).
!===============================================================================
function watt_spectrum(a, b) result(E_out)
real(8), intent(in) :: a ! Watt parameter a
real(8), intent(in) :: b ! Watt parameter b
real(8) :: E_out ! energy of emitted neutron
real(8) :: w ! sampled from Maxwellian
w = maxwell_spectrum(a)
E_out = w + a*a*b/4. + (TWO*prn() - ONE)*sqrt(a*a*b*w)
end function watt_spectrum
end module math

605
src/mpiio_interface.F90
View file

@ -5,9 +5,40 @@ module mpiio_interface
implicit none
integer :: mpi_fh ! MPI file handle
integer :: mpiio_err ! MPI error code
! Generic HDF5 write procedure interface
interface mpi_write_data
module procedure mpi_write_double
module procedure mpi_write_double_1Darray
module procedure mpi_write_double_2Darray
module procedure mpi_write_double_3Darray
module procedure mpi_write_double_4Darray
module procedure mpi_write_integer
module procedure mpi_write_integer_1Darray
module procedure mpi_write_integer_2Darray
module procedure mpi_write_integer_3Darray
module procedure mpi_write_integer_4Darray
module procedure mpi_write_long
module procedure mpi_write_string
end interface mpi_write_data
! Generic HDF5 read procedure interface
interface mpi_read_data
module procedure mpi_read_double
module procedure mpi_read_double_1Darray
module procedure mpi_read_double_2Darray
module procedure mpi_read_double_3Darray
module procedure mpi_read_double_4Darray
module procedure mpi_read_integer
module procedure mpi_read_integer_1Darray
module procedure mpi_read_integer_2Darray
module procedure mpi_read_integer_3Darray
module procedure mpi_read_integer_4Darray
module procedure mpi_read_long
module procedure mpi_read_string
end interface mpi_read_data
contains
!===============================================================================
@ -40,7 +71,7 @@ contains
! Determine access mode
open_mode = MPI_MODE_RDONLY
if (mode == 'w') then
open_mode = MPI_MODE_WRONLY
open_mode = ior(MPI_MODE_APPEND, MPI_MODE_WRONLY)
end if
! Create the file
@ -65,173 +96,505 @@ contains
! MPI_WRITE_INTEGER writes integer scalar data using MPI File I/O
!===============================================================================
subroutine mpi_write_integer(fh, buffer)
subroutine mpi_write_integer(fh, buffer, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: buffer ! data to write
integer, intent(in) :: fh ! file handle
integer, intent(in) :: buffer ! data to write
logical, intent(in) :: collect ! collective I/O
call MPI_FILE_WRITE(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_integer
!===============================================================================
! MPI_WRITE_INTEGER_1DARRAY writes integer 1-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_integer_1Darray(fh, buffer, length)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
integer, intent(in) :: buffer(:) ! data to write
call MPI_FILE_WRITE(fh, buffer, length, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_integer_1Darray
!===============================================================================
! MPI_WRITE_LONG writes long integer scalar data using MPI file I/O
!===============================================================================
subroutine mpi_write_long(fh, buffer)
integer, intent(in) :: fh ! file handle
integer(8), intent(in) :: buffer ! data to write
call MPI_FILE_WRITE(fh, buffer, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_long
!===============================================================================
! MPI_WRITE_DOUBLE writes double precision scalar data using MPI file I/O
!===============================================================================
subroutine mpi_write_double(fh, buffer)
integer, intent(in) :: fh ! file handle
real(8), intent(in) :: buffer ! data to write
call MPI_FILE_WRITE(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_double
!===============================================================================
! MPI_WRITE_DOUBLE_1DARRAY writes double precision 1-D array using MPI file I/O
!===============================================================================
subroutine mpi_write_double_1Darray(fh, buffer, length)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
real(8), intent(in) :: buffer(:) ! data to write
call MPI_FILE_WRITE(fh, buffer, length, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_double_1Darray
!===============================================================================
! MPI_WRITE_STRING writes string data using MPI file I/O
!===============================================================================
subroutine mpi_write_string(fh, buffer, length)
character(*), intent(in) :: buffer ! data to write
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of data
call MPI_FILE_WRITE(fh, buffer, length, MPI_CHARACTER, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_write_string
!===============================================================================
! MPI_READ_INTEGER reads integer scalar data using MPI file I/O
!===============================================================================
subroutine mpi_read_integer(fh, buffer)
subroutine mpi_read_integer(fh, buffer, collect)
integer, intent(in) :: fh ! file handle
integer, intent(inout) :: buffer ! read data to here
integer, intent(in) :: fh ! file handle
integer, intent(inout) :: buffer ! read data to here
logical, intent(in) :: collect ! collective I/O
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_integer
!===============================================================================
! MPI_WRITE_INTEGER_1DARRAY writes integer 1-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_integer_1Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
integer, intent(in) :: buffer(:) ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, length, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, length, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_integer_1Darray
!===============================================================================
! MPI_READ_INTEGER_1DARRAY reads integer 1-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_integer_1Darray(fh, buffer, length)
subroutine mpi_read_integer_1Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
integer, intent(inout) :: buffer(:) ! read data to here
logical, intent(in) :: collect ! collective I/O
call MPI_FILE_READ(fh, buffer, length, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, length, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, length, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_integer_1Darray
!===============================================================================
! MPI_WRITE_INTEGER_2DARRAY writes integer 2-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_integer_2Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(2) ! length of array
integer, intent(in) :: buffer(length(1),length(2)) ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_integer_2Darray
!===============================================================================
! MPI_READ_INTEGER_2DARRAY reads integer 2-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_integer_2Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(2) ! length of array
integer, intent(inout) :: buffer(length(1),length(2)) ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_integer_2Darray
!===============================================================================
! MPI_WRITE_INTEGER_3DARRAY writes integer 3-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_integer_3Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(3) ! length of array
integer, intent(in) :: buffer(length(1),length(2),&
length(3)) ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_integer_3Darray
!===============================================================================
! MPI_READ_INTEGER_3DARRAY reads integer 3-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_integer_3Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(3) ! length of array
integer, intent(inout) :: buffer(length(1),length(2), &
length(3)) ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_integer_3Darray
!===============================================================================
! MPI_WRITE_INTEGER_4DARRAY writes integer 4-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_integer_4Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(4) ! length of array
integer, intent(in) :: buffer(length(1),length(2),&
length(3),length(4)) ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_integer_4Darray
!===============================================================================
! MPI_READ_INTEGER_4DARRAY reads integer 4-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_integer_4Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(4) ! length of array
integer, intent(inout) :: buffer(length(1),length(2), &
length(3),length(4)) ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, product(length), MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_integer_4Darray
!===============================================================================
! MPI_WRITE_DOUBLE writes integer scalar data using MPI File I/O
!===============================================================================
subroutine mpi_write_double(fh, buffer, collect)
integer, intent(in) :: fh ! file handle
real(8), intent(in) :: buffer ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_double
!===============================================================================
! MPI_READ_DOUBLE reads integer scalar data using MPI file I/O
!===============================================================================
subroutine mpi_read_double(fh, buffer, collect)
integer, intent(in) :: fh ! file handle
real(8), intent(inout) :: buffer ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_double
!===============================================================================
! MPI_WRITE_DOUBLE_1DARRAY writes integer 1-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_double_1Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
real(8), intent(in) :: buffer(:) ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, length, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, length, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_double_1Darray
!===============================================================================
! MPI_READ_DOUBLE_1DARRAY reads integer 1-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_double_1Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
real(8), intent(inout) :: buffer(:) ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, length, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, length, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_double_1Darray
!===============================================================================
! MPI_WRITE_DOUBLE_2DARRAY writes integer 2-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_double_2Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(2) ! length of array
real(8), intent(in) :: buffer(length(1),length(2)) ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_double_2Darray
!===============================================================================
! MPI_READ_DOUBLE_2DARRAY reads integer 2-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_double_2Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(2) ! length of array
real(8), intent(inout) :: buffer(length(1),length(2)) ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_double_2Darray
!===============================================================================
! MPI_WRITE_DOUBLE_3DARRAY writes integer 3-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_double_3Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(3) ! length of array
real(8), intent(in) :: buffer(length(1),length(2),&
length(3)) ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_double_3Darray
!===============================================================================
! MPI_READ_DOUBLE_3DARRAY reads integer 3-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_double_3Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(3) ! length of array
real(8), intent(inout) :: buffer(length(1),length(2), &
length(3)) ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_double_3Darray
!===============================================================================
! MPI_WRITE_DOUBLE_4DARRAY writes integer 4-D array data using MPI File I/O
!===============================================================================
subroutine mpi_write_double_4Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(4) ! length of array
real(8), intent(in) :: buffer(length(1),length(2),&
length(3),length(4)) ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_double_4Darray
!===============================================================================
! MPI_READ_DOUBLE_4DARRAY reads integer 4-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_double_4Darray(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length(4) ! length of array
real(8), intent(inout) :: buffer(length(1),length(2), &
length(3),length(4)) ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, product(length), MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_double_4Darray
!===============================================================================
! MPI_WRITE_LONG writes long integer scalar data using MPI file I/O
!===============================================================================
subroutine mpi_write_long(fh, buffer, collect)
integer, intent(in) :: fh ! file handle
integer(8), intent(in) :: buffer ! data to write
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_write_long
!===============================================================================
! MPI_READ_LONG reads long integer scalar data using MPI file I/O
!===============================================================================
subroutine mpi_read_long(fh, buffer)
subroutine mpi_read_long(fh, buffer, collect)
integer, intent(in) :: fh ! file handle
integer(8), intent(inout) :: buffer ! read data to here
integer, intent(in) :: fh ! file handle
integer(8), intent(inout) :: buffer ! read data to here
logical, intent(in) :: collect ! collective I/O
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpiio_err)
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_long
!===============================================================================
! MPI_READ_DOUBLE reads double precision scalar data using MPI file I/O
! MPI_WRITE_STRING writes string data using MPI file I/O
!===============================================================================
subroutine mpi_read_double(fh, buffer)
subroutine mpi_write_string(fh, buffer, length, collect)
integer, intent(in) :: fh ! file handle
real(8), intent(inout) :: buffer ! read data to here
character(*), intent(in) :: buffer ! data to write
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of data
logical, intent(in) :: collect ! collective I/O
call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
if (collect) then
call MPI_FILE_WRITE_ALL(fh, buffer, length, MPI_CHARACTER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_WRITE(fh, buffer, length, MPI_CHARACTER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_double
!===============================================================================
! MPI_READ_DOUBLE_1DARRAY reads double precision 1-D array using MPI file I/O
!===============================================================================
subroutine mpi_read_double_1Darray(fh, buffer, length)
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of array
real(8), intent(inout) :: buffer(:) ! read data to here
call MPI_FILE_READ(fh, buffer, length, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end subroutine mpi_read_double_1Darray
end subroutine mpi_write_string
!===============================================================================
! MPI_READ_STRING reads string data using MPI file I/O
!===============================================================================
subroutine mpi_read_string(fh, buffer, length)
subroutine mpi_read_string(fh, buffer, length, collect)
character(*), intent(inout) :: buffer ! read data to here
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of string
character(*), intent(inout) :: buffer ! read data to here
integer, intent(in) :: fh ! file handle
integer, intent(in) :: length ! length of string
logical, intent(in) :: collect ! collective I/O
call MPI_FILE_READ(fh, buffer, length, MPI_CHARACTER, &
MPI_STATUS_IGNORE, mpiio_err)
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, length, MPI_CHARACTER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, length, MPI_CHARACTER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_string

8
src/output.F90
View file

@ -11,7 +11,7 @@ module output
use math, only: t_percentile
use mesh_header, only: StructuredMesh
use mesh, only: mesh_indices_to_bin, bin_to_mesh_indices
use particle_header, only: LocalCoord
use particle_header, only: LocalCoord, Particle
use plot_header
use string, only: upper_case, to_str
use tally_header, only: TallyObject
@ -236,7 +236,9 @@ contains
! PRINT_PARTICLE displays the attributes of a particle
!===============================================================================
subroutine print_particle()
subroutine print_particle(p)
type(Particle), intent(in) :: p
integer :: i ! index for coordinate levels
type(Cell), pointer :: c => null()
@ -302,8 +304,6 @@ contains
! Display weight, energy, grid index, and interpolation factor
write(ou,*) ' Weight = ' // to_str(p % wgt)
write(ou,*) ' Energy = ' // to_str(p % E)
write(ou,*) ' Energy grid index = ' // to_str(p % index_grid)
write(ou,*) ' Interpolation factor = ' // to_str(p % interp)
write(ou,*)
end subroutine print_particle

1576
src/output_interface.F90
View file

File diff suppressed because it is too large Load diff

65
src/particle_header.F90
View file

@ -1,6 +1,7 @@
module particle_header
use constants, only: NEUTRON, ONE, NONE, ZERO
use constants, only: NEUTRON, ONE, NONE, ZERO
use geometry_header, only: BASE_UNIVERSE
implicit none
@ -57,7 +58,8 @@ module particle_header
real(8) :: absorb_wgt ! weight absorbed for survival biasing
! What event last took place
integer :: event ! scatter, absorption, fission
logical :: fission ! did the particle cause implicit fission
integer :: event ! scatter, absorption
integer :: event_nuclide ! index in nuclides array
integer :: event_MT ! reaction MT
@ -65,10 +67,6 @@ module particle_header
integer :: n_bank ! number of fission sites banked
real(8) :: wgt_bank ! weight of fission sites banked
! Energy grid data
integer :: index_grid ! index on unionized energy grid
real(8) :: interp ! interpolation factor for energy grid
! Indices for various arrays
integer :: surface ! index for surface particle is on
integer :: cell_born ! index for cell particle was born in
@ -78,6 +76,9 @@ module particle_header
! Statistical data
integer :: n_collision ! # of collisions
contains
procedure :: initialize => initialize_particle
procedure :: clear => clear_particle
end type Particle
contains
@ -102,4 +103,56 @@ contains
end subroutine deallocate_coord
!===============================================================================
! INITIALIZE_PARTICLE sets default attributes for a particle from the source
! bank
!===============================================================================
subroutine initialize_particle(this)
class(Particle) :: this
! Clear coordinate lists
call this % clear()
! Set particle to neutron that's alive
this % type = NEUTRON
this % alive = .true.
! clear attributes
this % surface = NONE
this % cell_born = NONE
this % material = NONE
this % last_material = NONE
this % wgt = ONE
this % last_wgt = ONE
this % absorb_wgt = ZERO
this % n_bank = 0
this % wgt_bank = ZERO
this % n_collision = 0
this % fission = .false.
! Set up base level coordinates
allocate(this % coord0)
this % coord0 % universe = BASE_UNIVERSE
this % coord => this % coord0
end subroutine initialize_particle
!===============================================================================
! CLEAR_PARTICLE
!===============================================================================
subroutine clear_particle(this)
class(Particle) :: this
! remove any coordinate levels
call deallocate_coord(this % coord0)
! Make sure coord pointer is nullified
nullify(this % coord)
end subroutine clear_particle
end module particle_header

183
src/particle_restart.F90
View file

@ -2,152 +2,97 @@ module particle_restart
use, intrinsic :: ISO_FORTRAN_ENV
use bank_header, only: Bank
use bank_header, only: Bank
use constants
use geometry_header, only: BASE_UNIVERSE
use geometry_header, only: BASE_UNIVERSE
use global
use particle_header, only: deallocate_coord
use output, only: write_message, print_particle
use physics, only: transport
use random_lcg, only: set_particle_seed
use source, only: initialize_particle
#ifdef HDF5
use hdf5_interface
#endif
use output, only: write_message, print_particle
use output_interface, only: BinaryOutput
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use tracking, only: transport
implicit none
private
public :: run_particle_restart
#ifdef HDF5
integer(HID_T) :: hdf5_particle_file
#endif
! Short names for output and error units
integer :: ou = OUTPUT_UNIT
integer :: eu = ERROR_UNIT
! Binary file
type(BinaryOutput) :: pr
contains
#ifdef HDF5
!===============================================================================
! READ_HDF5_PARTICLE_RESTART
!===============================================================================
subroutine read_hdf5_particle_restart()
! write meessage
message = "Loading particle restart file " // trim(path_particle_restart) &
// "..."
call write_message(1)
! open hdf5 file
call h5fopen_f(path_particle_restart, H5F_ACC_RDONLY_F, hdf5_particle_file,&
hdf5_err)
! read data from file
call hdf5_read_integer(hdf5_particle_file, 'current_batch', current_batch)
call hdf5_read_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch)
call hdf5_read_integer(hdf5_particle_file, 'current_gen', current_gen)
call hdf5_read_long(hdf5_particle_file, 'n_particles', n_particles, hdf5_integer8_t)
call hdf5_read_long(hdf5_particle_file, 'id', p % id, hdf5_integer8_t)
call hdf5_read_double(hdf5_particle_file, 'weight', p % wgt)
call hdf5_read_double(hdf5_particle_file, 'energy', p % E)
dims1 = (/3/)
call h5ltread_dataset_double_f(hdf5_particle_file, 'xyz', p % coord % xyz, &
dims1, hdf5_err)
call h5ltread_dataset_double_f(hdf5_particle_file, 'uvw', p % coord % uvw, &
dims1, hdf5_err)
! set particle last attributes
p % last_wgt = p % wgt
p % last_xyz = p % coord % xyz
p % last_E = p % E
! close hdf5 file
call h5fclose_f(hdf5_particle_file, hdf5_err)
end subroutine read_hdf5_particle_restart
#endif
!===============================================================================
! READ_BINARY_PARTICLE_RESTART
!===============================================================================
subroutine read_binary_particle_restart()
integer :: filetype
integer :: revision
! write meessage
message = "Loading particle restart file " // trim(path_particle_restart) &
// "..."
call write_message(1)
! open file
open(UNIT=UNIT_PARTICLE, FILE=path_particle_restart, STATUS='old', &
ACCESS='stream')
! read data from file
read(UNIT_PARTICLE) filetype
read(UNIT_PARTICLE) revision
read(UNIT_PARTICLE) current_batch
read(UNIT_PARTICLE) gen_per_batch
read(UNIT_PARTICLE) current_gen
read(UNIT_PARTICLE) n_particles
read(UNIT_PARTICLE) p % id
read(UNIT_PARTICLE) p % wgt
read(UNIT_PARTICLE) p % E
read(UNIT_PARTICLE) p % coord % xyz
read(UNIT_PARTICLE) p % coord % uvw
! set particle last attributes
p % last_wgt = p % wgt
p % last_xyz = p % coord % xyz
p % last_E = p % E
! close hdf5 file
close(UNIT_PARTICLE)
end subroutine read_binary_particle_restart
!===============================================================================
! RUN_PARTICLE_RESTART
! RUN_PARTICLE_RESTART is the main routine that runs the particle restart
!===============================================================================
subroutine run_particle_restart()
integer(8) :: particle_seed
character(MAX_FILE_LEN) :: filename
type(Particle) :: p
! initialize the particle to be tracked
allocate(p)
call initialize_particle()
! Set verbosity high
verbosity = 10
! read in the restart information
#ifdef HDF5
call read_hdf5_particle_restart()
#else
call read_binary_particle_restart()
#endif
! Initialize the particle to be tracked
call p % initialize()
! set all tallies to 0 for now (just tracking errors)
! Read in the restart information
call read_particle_restart(p)
! Set all tallies to 0 for now (just tracking errors)
n_tallies = 0
! compute random number seed
! Compute random number seed
particle_seed = ((current_batch - 1)*gen_per_batch + &
current_gen - 1)*n_particles + p % id
call set_particle_seed(particle_seed)
! transport neutron
call transport()
! Transport neutron
call transport(p)
! write output if particle made it
call print_particle()
! Write output if particle made it
call print_particle(p)
end subroutine run_particle_restart
!===============================================================================
! READ_PARTICLE_RESTART reads the particle restart file
!===============================================================================
subroutine read_particle_restart(p)
integer :: int_scalar
type(Particle), intent(inout) :: p
! Write meessage
message = "Loading particle restart file " // trim(path_particle_restart) &
// "..."
call write_message(1)
! Open file
call pr % file_open(path_particle_restart, 'r')
! Read data from file
call pr % read_data(int_scalar, 'filetype')
call pr % read_data(int_scalar, 'revision')
call pr % read_data(current_batch, 'current_batch')
call pr % read_data(gen_per_batch, 'gen_per_batch')
call pr % read_data(current_gen, 'current_gen')
call pr % read_data(n_particles, 'n_particles')
call pr % read_data(p % id, 'id')
call pr % read_data(p % wgt, 'weight')
call pr % read_data(p % E, 'energy')
call pr % read_data(p % coord % xyz, 'xyz', length=3)
call pr % read_data(p % coord % uvw, 'uvw', length=3)
! Set particle last attributes
p % last_wgt = p % wgt
p % last_xyz = p % coord % xyz
p % last_E = p % E
! Close hdf5 file
call pr % file_close()
end subroutine read_particle_restart
end module particle_restart

140
src/particle_restart_write.F90
View file

@ -1,125 +1,65 @@
module particle_restart_write
use, intrinsic :: ISO_FORTRAN_ENV
use bank_header, only: Bank
use bank_header, only: Bank
use global
use string, only: to_str
#ifdef HDF5
use hdf5_interface
#endif
use output_interface, only: BinaryOutput
use particle_header, only: Particle
use string, only: to_str
implicit none
private
public :: write_particle_restart
#ifdef HDF5
integer(HID_T) :: hdf5_particle_file
#endif
! Binary output file
type(BinaryOutput) :: pr
contains
!===============================================================================
! WRITE_PARTICLE_RESTART
! WRITE_PARTICLE_RESTART is the main routine that writes out the particle file
!===============================================================================
subroutine write_particle_restart()
subroutine write_particle_restart(p)
type(Particle), intent(in) :: p
character(MAX_FILE_LEN) :: filename
type(Bank), pointer :: src => null()
! Dont write another restart file if in particle restart mode
if (run_mode == MODE_PARTICLE) return
! write out binary or HDF5 file
! Set up file name
filename = trim(path_output) // 'particle_' // trim(to_str(current_batch)) &
// '_' // trim(to_str(p % id))
#ifdef HDF5
call write_particle_restart_hdf5()
filename = trim(filename) // '.h5'
#else
call write_particle_restart_binary()
filename = trim(filename) // '.binary'
#endif
! Create file
call pr % file_create(filename)
! Get information about source particle
src => source_bank(current_work)
! Write data to file
call pr % write_data(FILETYPE_PARTICLE_RESTART, 'filetype')
call pr % write_data(REVISION_PARTICLE_RESTART, 'revision')
call pr % write_data(current_batch, 'current_batch')
call pr % write_data(gen_per_batch, 'gen_per_batch')
call pr % write_data(current_gen, 'current_gen')
call pr % write_data(n_particles, 'n_particles')
call pr % write_data(p % id, 'id')
call pr % write_data(src % wgt, 'weight')
call pr % write_data(src % E, 'energy')
call pr % write_data(src % xyz, 'xyz', length = 3)
call pr % write_data(src % uvw, 'uvw', length = 3)
! Close file
call pr % file_close()
end subroutine write_particle_restart
#ifdef HDF5
!===============================================================================
! WRITE_PARTICLE_RESTART_HDF5
!===============================================================================
subroutine write_particle_restart_hdf5()
character(MAX_FILE_LEN) :: filename
type(Bank), pointer :: src => null()
! set up file name
filename = trim(path_output) // 'particle_' // trim(to_str(current_batch)) &
// '_' // trim(to_str(current_work)) // '.h5'
! create hdf5 file
call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_particle_file, hdf5_err)
! get information about source particle
src => source_bank(current_work)
! write data to file
call hdf5_write_integer(hdf5_particle_file, 'filetype', &
FILETYPE_PARTICLE_RESTART)
call hdf5_write_integer(hdf5_particle_file, 'revision', &
REVISION_PARTICLE_RESTART)
call hdf5_write_integer(hdf5_particle_file, 'current_batch', current_batch)
call hdf5_write_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch)
call hdf5_write_integer(hdf5_particle_file, 'current_gen', current_gen)
call hdf5_write_long(hdf5_particle_file, 'n_particles', n_particles, hdf5_integer8_t)
call hdf5_write_long(hdf5_particle_file, 'id', p % id, hdf5_integer8_t)
call hdf5_write_double(hdf5_particle_file, 'weight', src % wgt)
call hdf5_write_double(hdf5_particle_file, 'energy', src % E)
dims1 = (/3/)
call h5ltmake_dataset_double_f(hdf5_particle_file, 'xyz', 1, dims1, &
src % xyz, hdf5_err)
call h5ltmake_dataset_double_f(hdf5_particle_file, 'uvw', 1, dims1, &
src % uvw, hdf5_err)
! close hdf5 file
call h5fclose_f(hdf5_particle_file, hdf5_err)
end subroutine write_particle_restart_hdf5
#endif
!===============================================================================
! WRITE_PARTICLE_RESTART_BINARY
!===============================================================================
subroutine write_particle_restart_binary()
character(MAX_FILE_LEN) :: filename
type(Bank), pointer :: src => null()
! set up file name
filename = trim(path_output) // 'particle_' // trim(to_str(current_batch)) &
// '_' // trim(to_str(current_work)) // '.binary'
! create hdf5 file
open(UNIT=UNIT_PARTICLE, FILE=filename, STATUS='replace', &
ACCESS='stream')
! get information about source particle
src => source_bank(current_work)
! write data to file
write(UNIT_PARTICLE) FILETYPE_PARTICLE_RESTART
write(UNIT_PARTICLE) REVISION_PARTICLE_RESTART
write(UNIT_PARTICLE) current_batch
write(UNIT_PARTICLE) gen_per_batch
write(UNIT_PARTICLE) current_gen
write(UNIT_PARTICLE) n_particles
write(UNIT_PARTICLE) p % id
write(UNIT_PARTICLE) src % wgt
write(UNIT_PARTICLE) src % E
write(UNIT_PARTICLE) src % xyz
write(UNIT_PARTICLE) src % uvw
! close hdf5 file
close(UNIT_PARTICLE)
end subroutine write_particle_restart_binary
end module particle_restart_write

27
src/particle_track.F90
View file

@ -7,7 +7,8 @@ module particle_track
use constants
use global
use output_interface, only: file_create, file_close, write_data
use output_interface, only: BinaryOutput
use particle_header, only: Particle
use string, only: to_str
implicit none
@ -30,7 +31,9 @@ contains
! WRITE_PARTICLE_TRACK copies particle position to an array.
!===============================================================================
subroutine write_particle_track()
subroutine write_particle_track(p)
type(Particle), intent(in) :: p
real(8) :: old_coords(3, n_tracks)
! Save the coordinates gathered in previous calls.
@ -59,34 +62,40 @@ contains
#ifdef HDF5
subroutine finalize_particle_track()
subroutine finalize_particle_track(p)
type(Particle), intent(in) :: p
character(MAX_FILE_LEN) :: fname
type(BinaryOutput) :: binout
fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) &
// '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) &
// '.h5'
call file_create(fname, 'serial')
call file_create(fname)
call write_data(coords, 'coordinates', length=(/3, n_tracks/))
call file_close('serial')
call file_close()
deallocate(coords)
end subroutine finalize_particle_track
#else
subroutine finalize_particle_track()
subroutine finalize_particle_track(p)
type(Particle), intent(in) :: p
character(MAX_FILE_LEN) :: fname
type(BinaryOutput) :: binout
integer :: i
real(8) :: flat_coords(3*n_tracks)
fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) &
// '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) &
// '.binary'
call file_create(fname, 'serial')
call binout % file_create(fname)
do i=1, n_tracks
flat_coords(3*i-2 : 3*i) = coords(:,i)
end do
call write_data(flat_coords, 'coordinates', length=3*n_tracks)
call file_close('serial')
call binout % write_data(flat_coords, 'coordinates', length=3*n_tracks)
call binout % file_close()
deallocate(coords)
end subroutine finalize_particle_track

1075
src/physics.F90
View file

File diff suppressed because it is too large Load diff

24
src/plot.F90
View file

@ -6,11 +6,10 @@ module plot
use geometry_header, only: Cell, BASE_UNIVERSE
use global
use output, only: write_message
use particle_header, only: deallocate_coord
use particle_header, only: deallocate_coord, Particle
use plot_header
use ppmlib, only: Image, init_image, allocate_image, &
deallocate_image, set_pixel
use source, only: initialize_particle
use string, only: to_str
implicit none
@ -49,8 +48,9 @@ contains
! current particle's position
!===============================================================================
subroutine position_rgb(pl, rgb, id)
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
@ -61,8 +61,8 @@ contains
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
call find_cell(found_cell)
if (check_overlaps) call check_cell_overlap()
call find_cell(p, found_cell)
if (check_overlaps) call check_cell_overlap(p)
if (.not. found_cell) then
! If no cell, revert to default color
@ -110,6 +110,7 @@ contains
real(8) :: out_pixel
real(8) :: xyz(3)
type(Image) :: img
type(Particle) :: p
! Initialize and allocate space for image
call init_image(img)
@ -142,8 +143,7 @@ contains
end if
! allocate and initialize particle
allocate(p)
call initialize_particle()
call p % initialize()
p % coord % xyz = xyz
p % coord % uvw = (/ 0.5, 0.5, 0.5 /)
p % coord % universe = BASE_UNIVERSE
@ -152,7 +152,7 @@ contains
do x = 1, img % width
! get pixel color
call position_rgb(pl, rgb, id)
call position_rgb(p, pl, rgb, id)
! Create a pixel at (x,y) with color (r,g,b)
call set_pixel(img, x, y, rgb(1), rgb(2), rgb(3))
@ -228,6 +228,7 @@ contains
integer :: id ! id of cell or material
real(8) :: vox(3) ! x, y, and z voxel widths
real(8) :: ll(3) ! lower left starting point for each sweep direction
type(Particle) :: p
! compute voxel widths in each direction
vox = pl % width/dble(pl % pixels)
@ -236,8 +237,7 @@ contains
ll = pl % origin - pl % width / 2.0
! allocate and initialize particle
allocate(p)
call initialize_particle()
call p % initialize()
p % coord0 % xyz = ll
p % coord0 % uvw = (/ 0.5, 0.5, 0.5 /)
p % coord0 % universe = BASE_UNIVERSE
@ -257,7 +257,7 @@ contains
do z = 1, pl % pixels(3)
! get voxel color
call position_rgb(pl, rgb, id)
call position_rgb(p, pl, rgb, id)
! write to plot file
write(UNIT_PLOT) id

57
src/source.F90
View file

@ -5,9 +5,9 @@ module source
use error, only: fatal_error
use geometry_header, only: BASE_UNIVERSE
use global
use math, only: maxwell_spectrum, watt_spectrum
use output, only: write_message
use particle_header, only: deallocate_coord
use physics, only: maxwell_spectrum, watt_spectrum
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed
use string, only: to_str
@ -47,7 +47,7 @@ contains
src => source_bank(i)
! initialize random number seed
id = bank_first + i - 1
id = work_index(rank) + i
call set_particle_seed(id)
! sample external source distribution
@ -149,23 +149,24 @@ contains
! GET_SOURCE_PARTICLE returns the next source particle
!===============================================================================
subroutine get_source_particle(index_source)
subroutine get_source_particle(p, index_source)
integer(8), intent(in) :: index_source
integer :: i ! iterator index
type(Particle), intent(inout) :: p
integer(8), intent(in) :: index_source
integer(8) :: particle_seed ! unique index for particle
integer :: i
type(Bank), pointer :: src => null()
! set defaults
call initialize_particle()
call p % initialize()
! Copy attributes from source to particle
src => source_bank(index_source)
call copy_source_attributes(src)
call copy_source_attributes(p, src)
! set identifier for particle
p % id = bank_first + index_source - 1
p % id = work_index(rank) + index_source
! set random number seed
particle_seed = (overall_gen - 1)*n_particles + p % id
@ -197,9 +198,10 @@ contains
! COPY_SOURCE_ATTRIBUTES
!===============================================================================
subroutine copy_source_attributes(src)
subroutine copy_source_attributes(p, src)
type(Bank), pointer :: src
type(Particle), intent(inout) :: p
type(Bank), pointer :: src
! copy attributes from source bank site
p % wgt = src % wgt
@ -212,37 +214,4 @@ contains
end subroutine copy_source_attributes
!===============================================================================
! INITIALIZE_PARTICLE sets default attributes for a particle from the source
! bank
!===============================================================================
subroutine initialize_particle()
! Set particle to neutron that's alive
p % type = NEUTRON
p % alive = .true.
! clear attributes
p % surface = NONE
p % cell_born = NONE
p % material = NONE
p % last_material = NONE
p % wgt = ONE
p % last_wgt = ONE
p % absorb_wgt = ZERO
p % n_bank = 0
p % wgt_bank = ZERO
p % n_collision = 0
! remove any original coordinates
call deallocate_coord(p % coord0)
! Set up base level coordinates
allocate(p % coord0)
p % coord0 % universe = BASE_UNIVERSE
p % coord => p % coord0
end subroutine initialize_particle
end module source

341
src/state_point.F90
View file

@ -22,6 +22,8 @@ module state_point
implicit none
type(BinaryOutput) :: sp ! statepoint/source output file
contains
!===============================================================================
@ -51,78 +53,78 @@ contains
message = "Creating state point " // trim(filename) // "..."
call write_message(1)
! Create statepoint file
call file_create(filename, 'serial')
if (master) then
! Create statepoint file
call sp % file_create(filename)
! Write file type
call write_data(FILETYPE_STATEPOINT, "filetype")
call sp % write_data(FILETYPE_STATEPOINT, "filetype")
! Write revision number for state point file
call write_data(REVISION_STATEPOINT, "revision")
call sp % write_data(REVISION_STATEPOINT, "revision")
! Write OpenMC version
call write_data(VERSION_MAJOR, "version_major")
call write_data(VERSION_MINOR, "version_minor")
call write_data(VERSION_RELEASE, "version_release")
call sp % write_data(VERSION_MAJOR, "version_major")
call sp % write_data(VERSION_MINOR, "version_minor")
call sp % write_data(VERSION_RELEASE, "version_release")
! Write current date and time
call write_data(time_stamp(), "date_and_time")
call sp % write_data(time_stamp(), "date_and_time")
! Write path to input
call write_data(path_input, "path")
call sp % write_data(path_input, "path")
! Write out random number seed
call write_data(seed, "seed")
call sp % write_data(seed, "seed")
! Write run information
call write_data(run_mode, "run_mode")
call write_data(n_particles, "n_particles")
call write_data(n_batches, "n_batches")
call sp % write_data(run_mode, "run_mode")
call sp % write_data(n_particles, "n_particles")
call sp % write_data(n_batches, "n_batches")
! Write out current batch number
call write_data(current_batch, "current_batch")
call sp % write_data(current_batch, "current_batch")
! Write out information for eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call write_data(n_inactive, "n_inactive")
call write_data(gen_per_batch, "gen_per_batch")
call write_data(k_generation, "k_generation", &
call sp % write_data(n_inactive, "n_inactive")
call sp % write_data(gen_per_batch, "gen_per_batch")
call sp % write_data(k_generation, "k_generation", &
length=current_batch*gen_per_batch)
call write_data(entropy, "entropy", length=current_batch*gen_per_batch)
call write_data(k_col_abs, "k_col_abs")
call write_data(k_col_tra, "k_col_tra")
call write_data(k_abs_tra, "k_abs_tra")
call write_data(k_combined, "k_combined", length=2)
call sp % write_data(entropy, "entropy", length=current_batch*gen_per_batch)
call sp % write_data(k_col_abs, "k_col_abs")
call sp % write_data(k_col_tra, "k_col_tra")
call sp % write_data(k_abs_tra, "k_abs_tra")
call sp % write_data(k_combined, "k_combined", length=2)
end if
! Write number of meshes
call write_data(n_meshes, "n_meshes", group="tallies")
call sp % write_data(n_meshes, "n_meshes", group="tallies")
! Write information for meshes
MESH_LOOP: do i = 1, n_meshes
call write_data(meshes(i) % id, "id", &
call sp % write_data(meshes(i) % id, "id", &
group="tallies/mesh" // to_str(i))
call write_data(meshes(i) % type, "type", &
call sp % write_data(meshes(i) % type, "type", &
group="tallies/mesh" // to_str(i))
call write_data(meshes(i) % n_dimension, "n_dimension", &
call sp % write_data(meshes(i) % n_dimension, "n_dimension", &
group="tallies/mesh" // to_str(i))
call write_data(meshes(i) % dimension, "dimension", &
call sp % write_data(meshes(i) % dimension, "dimension", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
call write_data(meshes(i) % lower_left, "lower_left", &
call sp % write_data(meshes(i) % lower_left, "lower_left", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
call write_data(meshes(i) % upper_right, "upper_right", &
call sp % write_data(meshes(i) % upper_right, "upper_right", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
call write_data(meshes(i) % width, "width", &
call sp % write_data(meshes(i) % width, "width", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
end do MESH_LOOP
! Write number of tallies
call write_data(n_tallies, "n_tallies", group="tallies")
call sp % write_data(n_tallies, "n_tallies", group="tallies")
! Write all tally information except results
TALLY_METADATA: do i = 1, n_tallies
@ -130,41 +132,41 @@ contains
t => tallies(i)
! Write id
call write_data(t % id, "id", group="tallies/tally" // to_str(i))
call sp % write_data(t % id, "id", group="tallies/tally" // to_str(i))
! Write number of realizations
call write_data(t % n_realizations, "n_realizations", &
call sp % write_data(t % n_realizations, "n_realizations", &
group="tallies/tally" // to_str(i))
! Write size of each tally
call write_data(t % total_score_bins, "total_score_bins", &
call sp % write_data(t % total_score_bins, "total_score_bins", &
group="tallies/tally" // to_str(i))
call write_data(t % total_filter_bins, "total_filter_bins", &
call sp % write_data(t % total_filter_bins, "total_filter_bins", &
group="tallies/tally" // to_str(i))
! Write number of filters
call write_data(t % n_filters, "n_filters", &
call sp % write_data(t % n_filters, "n_filters", &
group="tallies/tally" // to_str(i))
! Write filter information
FILTER_LOOP: do j = 1, t % n_filters
! Write type of filter
call write_data(t % filters(j) % type, "type", &
call sp % write_data(t % filters(j) % type, "type", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j))
! Write number of bins for this filter
call write_data(t % filters(j) % n_bins, "n_bins", &
call sp % write_data(t % filters(j) % n_bins, "n_bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j))
! Write bins
if (t % filters(j) % type == FILTER_ENERGYIN .or. &
t % filters(j) % type == FILTER_ENERGYOUT) then
call write_data(t % filters(j) % real_bins, "bins", &
call sp % write_data(t % filters(j) % real_bins, "bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
length=size(t % filters(j) % real_bins))
else
call write_data(t % filters(j) % int_bins, "bins", &
call sp % write_data(t % filters(j) % int_bins, "bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
length=size(t % filters(j) % int_bins))
end if
@ -172,7 +174,7 @@ contains
end do FILTER_LOOP
! Write number of nuclide bins
call write_data(t % n_nuclide_bins, "n_nuclide_bins", &
call sp % write_data(t % n_nuclide_bins, "n_nuclide_bins", &
group="tallies/tally" // to_str(i))
! Set up nuclide bin array and then write
@ -184,20 +186,20 @@ contains
temp_array(j) = t % nuclide_bins(j)
end if
end do NUCLIDE_LOOP
call write_data(temp_array, "nuclide_bins", &
call sp % write_data(temp_array, "nuclide_bins", &
group="tallies/tally" // to_str(i), length=t % n_nuclide_bins)
deallocate(temp_array)
! Write number of score bins, score bins, and scatt order
call write_data(t % n_score_bins, "n_score_bins", &
call sp % write_data(t % n_score_bins, "n_score_bins", &
group="tallies/tally" // to_str(i))
call write_data(t % score_bins, "score_bins", &
call sp % write_data(t % score_bins, "score_bins", &
group="tallies/tally" // to_str(i), length=t % n_score_bins)
call write_data(t % scatt_order, "scatt_order", &
call sp % write_data(t % scatt_order, "scatt_order", &
group="tallies/tally" // to_str(i), length=t % n_score_bins)
! Write number of user score bins
call write_data(t % n_user_score_bins, "n_user_score_bins", &
call sp % write_data(t % n_user_score_bins, "n_user_score_bins", &
group="tallies/tally" // to_str(i))
end do TALLY_METADATA
@ -214,18 +216,18 @@ contains
elseif (master) then
! Write number of global realizations
call write_data(n_realizations, "n_realizations")
call sp % write_data(n_realizations, "n_realizations")
! Write global tallies
call write_data(N_GLOBAL_TALLIES, "n_global_tallies")
call write_tally_result(global_tallies, "global_tallies", &
call sp % write_data(N_GLOBAL_TALLIES, "n_global_tallies")
call sp % write_tally_result(global_tallies, "global_tallies", &
n1=N_GLOBAL_TALLIES, n2=1)
! Write tallies
if (tallies_on) then
! Indicate that tallies are on
call write_data(1, "tallies_present", group="tallies")
call sp % write_data(1, "tallies_present", group="tallies")
! Write all tally results
TALLY_RESULTS: do i = 1, n_tallies
@ -234,7 +236,7 @@ contains
t => tallies(i)
! Write sum and sum_sq for each bin
call write_tally_result(t % results, "results", &
call sp % write_tally_result(t % results, "results", &
group="tallies/tally" // to_str(i), &
n1=size(t % results, 1), n2=size(t % results, 2))
@ -243,10 +245,13 @@ contains
else
! Indicate tallies are off
call write_data(0, "tallies_present", group="tallies")
call sp % write_data(0, "tallies_present", group="tallies")
end if
! Close the file for serial writing
call sp % file_close()
end if
! Check for eigenvalue calculation
@ -255,9 +260,6 @@ contains
! Check for writing source out separately
if (source_separate) then
! Close statepoint file
call file_close('serial')
! Set filename for source
filename = trim(path_output) // 'source.' // &
trim(to_str(current_batch))
@ -271,30 +273,21 @@ contains
message = "Creating source file " // trim(filename) // "..."
call write_message(1)
! Create statepoint file
call file_create(filename, 'parallel')
! Create source file
call sp % file_create(filename, serial = .false.)
#ifdef HDF5
# ifdef MPI
else
! Close HDF5 serial file and reopen in parallel
call file_close('serial')
call file_open(filename, 'parallel', 'w')
# endif
#endif
! Reopen state point file in parallel
call sp % file_open(filename, 'w', serial = .false.)
end if
! Write out source
call write_source_bank()
call sp % write_source_bank()
! Close file, all files in parallel mode
call file_close('parallel') ! even if no MPI, this will work for HDF5
else
! Close file if not in eigenvalue mode or no source writing
call file_close('serial')
! Close file
call sp % file_close()
end if
@ -323,10 +316,10 @@ contains
if (master) then
! Write number of realizations
call write_data(n_realizations, "n_realizations")
call sp % write_data(n_realizations, "n_realizations")
! Write number of global tallies
call write_data(N_GLOBAL_TALLIES, "n_global_tallies")
call sp % write_data(N_GLOBAL_TALLIES, "n_global_tallies")
end if
! Copy global tallies into temporary array for reducing
@ -355,7 +348,7 @@ contains
! Write out global tallies sum and sum_sq
call write_tally_result(tallyresult_temp, "global_tallies", &
call sp % write_tally_result(tallyresult_temp, "global_tallies", &
n1=N_GLOBAL_TALLIES, n2=1)
! Deallocate temporary tally result
@ -371,7 +364,7 @@ contains
if (tallies_on) then
! Indicate that tallies are on
if (master) then
call write_data(1, "tallies_present", group="tallies")
call sp % write_data(1, "tallies_present", group="tallies")
end if
! Write all tally results
@ -409,7 +402,7 @@ contains
tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:)
! Write reduced tally results to file
call write_tally_result(t % results, "results", &
call sp % write_tally_result(t % results, "results", &
group="tallies/tally" // to_str(i), n1=m, n2=n)
! Deallocate temporary tally result
@ -428,7 +421,7 @@ contains
else
if (master) then
! Indicate that tallies are off
call write_data(0, "tallies_present", group="tallies")
call sp % write_data(0, "tallies_present", group="tallies")
end if
end if
@ -455,14 +448,14 @@ contains
call write_message(1)
! Open file for reading
call file_open(path_state_point, 'parallel', 'r')
call sp % file_open(path_state_point, 'r')
! Read filetype
call read_data(int_array(1), "filetype", option="collective")
call sp % read_data(int_array(1), "filetype")
! Read revision number for state point file and make sure it matches with
! current version
call read_data(int_array(1), "revision", option="collective")
call sp % read_data(int_array(1), "revision")
if (int_array(1) /= REVISION_STATEPOINT) then
message = "State point version does not match current version " &
// "in OpenMC."
@ -470,9 +463,9 @@ contains
end if
! Read OpenMC version
call read_data(int_array(1), "version_major", option="collective")
call read_data(int_array(2), "version_minor", option="collective")
call read_data(int_array(3), "version_release", option="collective")
call sp % read_data(int_array(1), "version_major")
call sp % read_data(int_array(2), "version_minor")
call sp % read_data(int_array(3), "version_release")
if (int_array(1) /= VERSION_MAJOR .or. int_array(2) /= VERSION_MINOR &
.or. int_array(3) /= VERSION_RELEASE) then
message = "State point file was created with a different version " &
@ -481,70 +474,68 @@ contains
end if
! Read date and time
call read_data(current_time, "date_and_time", option="collective")
call sp % read_data(current_time, "date_and_time")
! Read path to input
call read_data(path_temp, "path", option="collective")
call sp % read_data(path_temp, "path")
! Read and overwrite random number seed
call read_data(seed, "seed", option="collective")
call sp % read_data(seed, "seed")
! Read and overwrite run information except number of batches
call read_data(run_mode, "run_mode", option="collective")
call read_data(n_particles, "n_particles", option="collective")
call read_data(int_array(1), "n_batches", option="collective")
call sp % read_data(run_mode, "run_mode")
call sp % read_data(n_particles, "n_particles")
call sp % read_data(int_array(1), "n_batches")
! Take maximum of statepoint n_batches and input n_batches
n_batches = max(n_batches, int_array(1))
! Read batch number to restart at
call read_data(restart_batch, "current_batch", option="collective")
call sp % read_data(restart_batch, "current_batch")
! Read information specific to eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call read_data(int_array(1), "n_inactive", option="collective")
call read_data(gen_per_batch, "gen_per_batch", option="collective")
call read_data(k_generation, "k_generation", &
length=restart_batch*gen_per_batch, option="collective")
call read_data(entropy, "entropy", length=restart_batch*gen_per_batch, &
option="collective")
call read_data(k_col_abs, "k_col_abs", option="collective")
call read_data(k_col_tra, "k_col_tra", option="collective")
call read_data(k_abs_tra, "k_abs_tra", option="collective")
call read_data(real_array(1:2), "k_combined", length=2, &
option="collective")
call sp % read_data(int_array(1), "n_inactive")
call sp % read_data(gen_per_batch, "gen_per_batch")
call sp % read_data(k_generation, "k_generation", &
length=restart_batch*gen_per_batch)
call sp % read_data(entropy, "entropy", length=restart_batch*gen_per_batch)
call sp % read_data(k_col_abs, "k_col_abs")
call sp % read_data(k_col_tra, "k_col_tra")
call sp % read_data(k_abs_tra, "k_abs_tra")
call sp % read_data(real_array(1:2), "k_combined", length=2)
! Take maximum of statepoint n_inactive and input n_inactive
n_inactive = max(n_inactive, int_array(1))
end if
! Read number of meshes
call read_data(n_meshes, "n_meshes", group="tallies", option="collective")
call sp % read_data(n_meshes, "n_meshes", group="tallies")
! Read and overwrite mesh information
MESH_LOOP: do i = 1, n_meshes
call read_data(meshes(i) % id, "id", &
group="tallies/mesh" // to_str(i), option="collective")
call read_data(meshes(i) % type, "type", &
group="tallies/mesh" // to_str(i), option="collective")
call read_data(meshes(i) % n_dimension, "n_dimension", &
group="tallies/mesh" // to_str(i), option="collective")
call read_data(meshes(i) % dimension, "dimension", &
call sp % read_data(meshes(i) % id, "id", &
group="tallies/mesh" // to_str(i))
call sp % read_data(meshes(i) % type, "type", &
group="tallies/mesh" // to_str(i))
call sp % read_data(meshes(i) % n_dimension, "n_dimension", &
group="tallies/mesh" // to_str(i))
call sp % read_data(meshes(i) % dimension, "dimension", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension, option="collective")
call read_data(meshes(i) % lower_left, "lower_left", &
length=meshes(i) % n_dimension)
call sp % read_data(meshes(i) % lower_left, "lower_left", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension, option="collective")
call read_data(meshes(i) % upper_right, "upper_right", &
length=meshes(i) % n_dimension)
call sp % read_data(meshes(i) % upper_right, "upper_right", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension, option="collective")
call read_data(meshes(i) % width, "width", &
length=meshes(i) % n_dimension)
call sp % read_data(meshes(i) % width, "width", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension, option="collective")
length=meshes(i) % n_dimension)
end do MESH_LOOP
! Read and overwrite number of tallies
call read_data(n_tallies, "n_tallies", group="tallies", option="collective")
call sp % read_data(n_tallies, "n_tallies", group="tallies")
! Read in tally metadata
TALLY_METADATA: do i = 1, n_tallies
@ -553,18 +544,17 @@ contains
t => tallies(i)
! Read tally id
call read_data(t % id, "id", group="tallies/tally" // to_str(i), &
option="collective")
call sp % read_data(t % id, "id", group="tallies/tally" // to_str(i))
! Read number of realizations
call read_data(t % n_realizations, "n_realizations", &
group="tallies/tally" // to_str(i), option="collective")
call sp % read_data(t % n_realizations, "n_realizations", &
group="tallies/tally" // to_str(i))
! Read size of tally results
call read_data(int_array(1), "total_score_bins", &
group="tallies/tally" // to_str(i), option="collective")
call read_data(int_array(2), "total_filter_bins", &
group="tallies/tally" // to_str(i), option="collective")
call sp % read_data(int_array(1), "total_score_bins", &
group="tallies/tally" // to_str(i))
call sp % read_data(int_array(2), "total_filter_bins", &
group="tallies/tally" // to_str(i))
! Check size of tally results array
if (int_array(1) /= t % total_score_bins .and. &
@ -574,45 +564,42 @@ contains
end if
! Read number of filters
call read_data(t % n_filters, "n_filters", &
group="tallies/tally" // to_str(i), option="collective")
call sp % read_data(t % n_filters, "n_filters", &
group="tallies/tally" // to_str(i))
! Read filter information
FILTER_LOOP: do j = 1, t % n_filters
! Read type of filter
call read_data(t % filters(j) % type, "type", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
option="collective")
call sp % read_data(t % filters(j) % type, "type", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j))
! Read number of bins for this filter
call read_data(t % filters(j) % n_bins, "n_bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
option="collective")
call sp % read_data(t % filters(j) % n_bins, "n_bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j))
! Read bins
if (t % filters(j) % type == FILTER_ENERGYIN .or. &
t % filters(j) % type == FILTER_ENERGYOUT) then
call read_data(t % filters(j) % real_bins, "bins", &
call sp % read_data(t % filters(j) % real_bins, "bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
length=size(t % filters(j) % real_bins), option="collective")
length=size(t % filters(j) % real_bins))
else
call read_data(t % filters(j) % int_bins, "bins", &
call sp % read_data(t % filters(j) % int_bins, "bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
length=size(t % filters(j) % int_bins), option="collective")
length=size(t % filters(j) % int_bins))
end if
end do FILTER_LOOP
! Read number of nuclide bins
call read_data(t % n_nuclide_bins, "n_nuclide_bins", &
group="tallies/tally" // to_str(i), option="collective")
call sp % read_data(t % n_nuclide_bins, "n_nuclide_bins", &
group="tallies/tally" // to_str(i))
! Set up nuclide bin array and then write
allocate(temp_array(t % n_nuclide_bins))
call read_data(temp_array, "nuclide_bins", &
group="tallies/tally" // to_str(i), length=t % n_nuclide_bins, &
option="collective")
call sp % read_data(temp_array, "nuclide_bins", &
group="tallies/tally" // to_str(i), length=t % n_nuclide_bins)
NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins
if (temp_array(j) > 0) then
nuclides(t % nuclide_bins(j)) % zaid = temp_array(j)
@ -623,18 +610,16 @@ contains
deallocate(temp_array)
! Write number of score bins, score bins, and scatt order
call read_data(t % n_score_bins, "n_score_bins", &
group="tallies/tally" // to_str(i), option="collective")
call read_data(t % score_bins, "score_bins", &
group="tallies/tally" // to_str(i), length=t % n_score_bins, &
option="collective")
call read_data(t % scatt_order, "scatt_order", &
group="tallies/tally" // to_str(i), length=t % n_score_bins, &
option="collective")
call sp % read_data(t % n_score_bins, "n_score_bins", &
group="tallies/tally" // to_str(i))
call sp % read_data(t % score_bins, "score_bins", &
group="tallies/tally" // to_str(i), length=t % n_score_bins)
call sp % read_data(t % scatt_order, "scatt_order", &
group="tallies/tally" // to_str(i), length=t % n_score_bins)
! Write number of user score bins
call read_data(t % n_user_score_bins, "n_user_score_bins", &
group="tallies/tally" // to_str(i), option="collective")
call sp % read_data(t % n_user_score_bins, "n_user_score_bins", &
group="tallies/tally" // to_str(i))
end do TALLY_METADATA
@ -642,21 +627,21 @@ contains
if (master) then
! Read number of realizations for global tallies
call read_data(n_realizations, "n_realizations")
call sp % read_data(n_realizations, "n_realizations", collect=.false.)
! Read number of global tallies
call read_data(int_array(1), "n_global_tallies")
call sp % read_data(int_array(1), "n_global_tallies", collect=.false.)
if (int_array(1) /= N_GLOBAL_TALLIES) then
message = "Number of global tallies does not match in state point."
call fatal_error()
end if
! Read global tally data
call read_tally_result(global_tallies, "global_tallies", &
call sp % read_tally_result(global_tallies, "global_tallies", &
n1=N_GLOBAL_TALLIES, n2=1)
! Check if tally results are present
call read_data(int_array(1), "tallies_present", group="tallies")
call sp % read_data(int_array(1), "tallies_present", group="tallies", collect=.false.)
! Read in sum and sum squared
if (int_array(1) == 1) then
@ -666,12 +651,19 @@ contains
t => tallies(i)
! Read sum and sum_sq for each bin
call read_tally_result(t % results, "results", &
call sp % read_tally_result(t % results, "results", &
group="tallies/tally" // to_str(i), &
n1=size(t % results, 1), n2=size(t % results, 2))
end do TALLY_RESULTS
end if
#ifdef MPI
! If using MPI, file needs to be closed and reopened in parallel
! If serial, we cannot close the file or we will lose our file position
call sp % file_close()
# endif
end if
! Read source if in eigenvalue mode
@ -681,7 +673,7 @@ contains
if (source_separate) then
! Close statepoint file
call file_close('parallel')
call sp % file_close()
! Set filename for source
filename = trim(path_output) // 'source.' // &
@ -696,28 +688,27 @@ contains
message = "Loading source file " // trim(filename) // "..."
call write_message(1)
! Create statepoint file
call file_open(filename, 'parallel', 'r')
! Open source file
call sp % file_open(filename, 'r', serial = .false.)
else
#ifdef MPI
! Reopen statepoint file in parallel, but only if MPI
! We will compute the position where the source begins
call sp % file_open(path_state_point, 'r', serial = .false.)
#endif
end if
! Write out source
call read_source_bank()
! Close file
if (source_separate) then
call file_close('parallel')
else
call file_close('parallel')
end if
else
! Close file if not in eigenvalue mode
call file_close('parallel')
call sp % read_source_bank()
end if
! Close file
call sp % file_close()
end subroutine load_state_point
subroutine read_source

39
src/string.F90
View file

@ -322,49 +322,18 @@ contains
end function str_to_int
!===============================================================================
! STR_TO_REAL converts an arbitrary string to a real(8). Generally this function
! is intended for strings for which the exact format is not known. If the format
! of the number is known a priori, the appropriate format descriptor should be
! used in lieu of this routine because of the extra overhead.
!
! Arguments:
! string = character(*) containing number to convert
! STR_TO_REAL converts an arbitrary string to a real(8)
!===============================================================================
function str_to_real(string) result(num)
character(*), intent(in) :: string
real(8) :: num
real(8) :: num
integer :: index_decimal ! index of decimal point
integer :: index_exponent ! index of exponent character
integer :: w ! total field width
integer :: d ! number of digits to right of decimal point
integer :: ioError
character(8) :: fmt ! format for reading string
! Determine total field width
w = len_trim(string)
! Determine number of digits to right of decimal point
index_decimal = index(string, '.')
index_exponent = max(index(string, 'd'), index(string, 'D'), &
index(string, 'e'), index(string, 'E'))
if (index_decimal > 0) then
if (index_exponent > 0) then
d = index_exponent - index_decimal - 1
else
d = w - index_decimal
end if
else
d = 0
end if
! Create format specifier for reading string
write(fmt, '("(E",I2,".",I2,")")') w, d
integer :: ioError
! Read string
read(UNIT=string, FMT=fmt, IOSTAT=ioError) num
read(UNIT=string, FMT=*, IOSTAT=ioError) num
if (ioError > 0) num = ERROR_REAL
end function str_to_real

59
src/tally.F90
View file

@ -10,7 +10,7 @@ module tally
mesh_intersects_2d, mesh_intersects_3d
use mesh_header, only: StructuredMesh
use output, only: header
use particle_header, only: LocalCoord
use particle_header, only: LocalCoord, Particle
use search, only: binary_search
use string, only: to_str
use tally_header, only: TallyResult, TallyMapItem, TallyMapElement
@ -32,7 +32,9 @@ contains
! triggered at every collision, not every event
!===============================================================================
subroutine score_analog_tally()
subroutine score_analog_tally(p)
type(Particle), intent(in) :: p
integer :: i
integer :: i_tally
@ -70,7 +72,7 @@ contains
! =======================================================================
! DETERMINE SCORING BIN COMBINATION
call get_scoring_bins(i_tally, found_bin)
call get_scoring_bins(p, i_tally, found_bin)
if (.not. found_bin) cycle
! =======================================================================
@ -256,7 +258,7 @@ contains
else
! Skip any non-fission events
if (p % event /= EVENT_FISSION) cycle SCORE_LOOP
if (.not. p % fission) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for the
@ -276,7 +278,7 @@ contains
else
! Skip any non-fission events
if (p % event /= EVENT_FISSION) cycle SCORE_LOOP
if (.not. p % fission) cycle SCORE_LOOP
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Normally, we only need to make contributions to one scoring
@ -285,7 +287,7 @@ contains
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout(t, score_index)
call score_fission_eout(p, t, score_index)
cycle SCORE_LOOP
else
@ -312,7 +314,7 @@ contains
else
! Skip any non-fission events
if (p % event /= EVENT_FISSION) cycle SCORE_LOOP
if (.not. p % fission) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
@ -364,8 +366,9 @@ contains
! neutrons produced with different energies.
!===============================================================================
subroutine score_fission_eout(t, i_score)
subroutine score_fission_eout(p, t, i_score)
type(Particle), intent(in) :: p
type(TallyObject), pointer :: t
integer, intent(in) :: i_score ! index for score
@ -424,9 +427,10 @@ contains
! that require post-collision information.
!===============================================================================
subroutine score_tracklength_tally(distance)
subroutine score_tracklength_tally(p, distance)
real(8), intent(in) :: distance
type(Particle), intent(in) :: p
real(8), intent(in) :: distance
integer :: i
integer :: i_tally
@ -464,14 +468,14 @@ contains
! since multiple bins can be scored to with a single track
if (t % find_filter(FILTER_MESH) > 0) then
call score_tl_on_mesh(i_tally, distance)
call score_tl_on_mesh(p, i_tally, distance)
cycle
end if
! =======================================================================
! DETERMINE SCORING BIN COMBINATION
call get_scoring_bins(i_tally, found_bin)
call get_scoring_bins(p, i_tally, found_bin)
if (.not. found_bin) cycle
! =======================================================================
@ -485,7 +489,7 @@ contains
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
if (t % all_nuclides) then
call score_all_nuclides(i_tally, flux, filter_index)
call score_all_nuclides(p, i_tally, flux, filter_index)
else
NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins
@ -721,11 +725,12 @@ contains
! the user requests <nuclides>all</nuclides>.
!===============================================================================
subroutine score_all_nuclides(i_tally, flux, filter_index)
subroutine score_all_nuclides(p, i_tally, flux, filter_index)
integer, intent(in) :: i_tally
real(8), intent(in) :: flux
integer, intent(in) :: filter_index
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
real(8), intent(in) :: flux
integer, intent(in) :: filter_index
integer :: i ! loop index for nuclides in material
integer :: j ! loop index for scoring bin types
@ -946,10 +951,11 @@ contains
! these tallies, it is possible to score to multiple mesh cells for each track.
!===============================================================================
subroutine score_tl_on_mesh(i_tally, d_track)
subroutine score_tl_on_mesh(p, i_tally, d_track)
integer, intent(in) :: i_tally
real(8), intent(in) :: d_track
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
real(8), intent(in) :: d_track
integer :: i ! loop index for filter/score bins
integer :: j ! loop index for direction
@ -1145,7 +1151,7 @@ contains
if (t % all_nuclides) then
! Score reaction rates for each nuclide in material
call score_all_nuclides(i_tally, flux, filter_index)
call score_all_nuclides(p, i_tally, flux, filter_index)
else
NUCLIDE_BIN_LOOP: do b = 1, t % n_nuclide_bins
@ -1260,10 +1266,11 @@ contains
! for a tally based on the particle's current attributes.
!===============================================================================
subroutine get_scoring_bins(i_tally, found_bin)
subroutine get_scoring_bins(p, i_tally, found_bin)
integer, intent(in) :: i_tally
logical, intent(out) :: found_bin
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
logical, intent(out) :: found_bin
integer :: i ! loop index for filters
integer :: n ! number of bins for single filter
@ -1370,7 +1377,9 @@ contains
! determining which mesh surfaces were crossed
!===============================================================================
subroutine score_surface_current()
subroutine score_surface_current(p)
type(Particle), intent(in) :: p
integer :: i
integer :: i_tally

210
src/tracking.F90 Normal file
View file

@ -0,0 +1,210 @@
module tracking
use cross_section, only: calculate_xs
use error, only: fatal_error, warning
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
cross_lattice, check_cell_overlap
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use output, only: write_message
use particle_header, only: LocalCoord, Particle
use particle_track, only: initialize_particle_track, write_particle_track, &
finalize_particle_track
use physics, only: collision
use random_lcg, only: prn
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_surface_current
contains
!===============================================================================
! TRANSPORT encompasses the main logic for moving a particle through geometry.
!===============================================================================
subroutine transport(p)
type(Particle), intent(inout) :: p
integer :: surface_crossed ! surface which particle is on
integer :: lattice_crossed ! lattice boundary which particle crossed
integer :: last_cell ! most recent cell particle was in
integer :: n_event ! number of collisions/crossings
real(8) :: d_boundary ! distance to nearest boundary
real(8) :: d_collision ! sampled distance to collision
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
type(LocalCoord), pointer :: coord => null()
! Display message if high verbosity or trace is on
if (verbosity >= 9 .or. trace) then
message = "Simulating Particle " // trim(to_str(p % id))
call write_message()
end if
! If the cell hasn't been determined based on the particle's location,
! initiate a search for the current cell
if (p % coord % cell == NONE) then
call find_cell(p, found_cell)
! Particle couldn't be located
if (.not. found_cell) then
message = "Could not locate particle " // trim(to_str(p % id))
call fatal_error()
end if
! set birth cell attribute
p % cell_born = p % coord % cell
end if
! Initialize number of events to zero
n_event = 0
! Add paricle's starting weight to count for normalizing tallies later
total_weight = total_weight + p % wgt
! Force calculation of cross-sections by setting last energy to zero
micro_xs % last_E = ZERO
! Prepare to write out particle track.
if (write_track) then
call initialize_particle_track()
endif
do while (p % alive)
! Write particle track.
if (write_track) call write_particle_track(p)
if (check_overlaps) call check_cell_overlap(p)
! Calculate microscopic and macroscopic cross sections -- note: if the
! material is the same as the last material and the energy of the
! particle hasn't changed, we don't need to lookup cross sections again.
if (p % material /= p % last_material) call calculate_xs(p)
! Find the distance to the nearest boundary
call distance_to_boundary(p, d_boundary, surface_crossed, lattice_crossed)
! Sample a distance to collision
if (material_xs % total == ZERO) then
d_collision = INFINITY
else
d_collision = -log(prn()) / material_xs % total
end if
! Select smaller of the two distances
distance = min(d_boundary, d_collision)
! Advance particle
coord => p % coord0
do while (associated(coord))
coord % xyz = coord % xyz + distance * coord % uvw
coord => coord % next
end do
! Score track-length tallies
if (active_tracklength_tallies % size() > 0) &
call score_tracklength_tally(p, distance)
! Score track-length estimate of k-eff
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + p % wgt * distance * &
material_xs % nu_fission
if (d_collision > d_boundary) then
! ====================================================================
! PARTICLE CROSSES SURFACE
last_cell = p % coord % cell
p % coord % cell = NONE
if (lattice_crossed /= NONE) then
! Particle crosses lattice boundary
p % surface = NONE
call cross_lattice(p, lattice_crossed)
p % event = EVENT_LATTICE
else
! Particle crosses surface
p % surface = surface_crossed
call cross_surface(p, last_cell)
p % event = EVENT_SURFACE
end if
else
! ====================================================================
! PARTICLE HAS COLLISION
! Score collision estimate of keff
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + p % wgt * &
material_xs % nu_fission / material_xs % total
! score surface current tallies -- this has to be done before the collision
! since the direction of the particle will change and we need to use the
! pre-collision direction to figure out what mesh surfaces were crossed
if (active_current_tallies % size() > 0) call score_surface_current(p)
! Clear surface component
p % surface = NONE
call collision(p)
! Score collision estimator tallies -- this is done after a collision
! has occurred rather than before because we need information on the
! outgoing energy for any tallies with an outgoing energy filter
if (active_analog_tallies % size() > 0) call score_analog_tally(p)
! Reset banked weight during collision
p % n_bank = 0
p % wgt_bank = ZERO
! Reset fission logical
p % fission = .false.
! Save coordinates for tallying purposes
p % last_xyz = p % coord0 % xyz
! Set last material to none since cross sections will need to be
! re-evaluated
p % last_material = NONE
! Set all uvws to base level -- right now, after a collision, only the
! base level uvws are changed
coord => p % coord0
do while(associated(coord % next))
if (coord % next % rotated) then
! If next level is rotated, apply rotation matrix
coord % next % uvw = matmul(cells(coord % cell) % &
rotation, coord % uvw)
else
! Otherwise, copy this level's direction
coord % next % uvw = coord % uvw
end if
! Advance coordinate level
coord => coord % next
end do
end if
! If particle has too many events, display warning and kill it
n_event = n_event + 1
if (n_event == MAX_EVENTS) then
message = "Particle " // trim(to_str(p%id)) // " underwent maximum &
&number of events."
call warning()
p % alive = .false.
end if
end do
! Finish particle track output.
if (write_track) then
call write_particle_track(p)
call finalize_particle_track(p)
endif
end subroutine transport
end module tracking

4
src/utils/build_dependencies.py
View file

@ -9,9 +9,9 @@ for src in glob.iglob('*.F90'):
module = src.strip('.F90')
deps = set()
d = re.findall(r'\n\s*use\s+(\w+)',
open(src,'r').read())
open(src, 'r').read())
for name in d:
if name in ['mpi','hdf5','h5lt']:
if name in ['mpi', 'hdf5', 'h5lt']:
continue
if name.startswith('xml_data_'):
name = name.replace('xml_data_', 'templates/')

39
src/utils/convert_xsdata.py
View file

@ -6,6 +6,7 @@ from xml.dom.minidom import getDOMImplementation
types = {1: "neutron", 2: "dosimetry", 3: "thermal"}
class Xsdata(object):
def __init__(self, filename):
@ -94,7 +95,7 @@ class XsdataTable(object):
continue
string = " ".join(self.alias)
else:
string = "{0}".format(getattr(self,attribute))
string = "{0}".format(getattr(self, attribute))
# Skip metastable and binary if 0
if attribute == "metastable" and self.metastable == 0:
@ -109,7 +110,7 @@ class XsdataTable(object):
# node.appendChild(nodeAttr)
node.setAttribute(attribute, string)
return node
if __name__ == '__main__':
# Read command line arguments
@ -121,25 +122,25 @@ if __name__ == '__main__':
# Read xsdata and create XML document object
xsdataObject = Xsdata(xsdataFile)
doc = xsdataObject.to_xml()
# Reduce number of lines
lines = doc.toprettyxml(indent=' ')
lines = lines.replace('<alias>\n ','<alias>')
lines = lines.replace('\n </alias>','</alias>')
lines = lines.replace('<zaid>\n ','<zaid>')
lines = lines.replace('\n </zaid>','</zaid>')
lines = lines.replace('<type>\n ','<type>')
lines = lines.replace('\n </type>','</type>')
lines = lines.replace('<awr>\n ','<awr>')
lines = lines.replace('\n </awr>','</awr>')
lines = lines.replace('<temperature>\n ','<temperature>')
lines = lines.replace('\n </temperature>','</temperature>')
lines = lines.replace('<path>\n ','<path>')
lines = lines.replace('\n </path>','</path>')
lines = lines.replace('<metastable>\n ','<metastable>')
lines = lines.replace('\n </metastable>','</metastable>')
lines = lines.replace('<binary>\n ','<binary>')
lines = lines.replace('\n </binary>','</binary>')
lines = lines.replace('<alias>\n ', '<alias>')
lines = lines.replace('\n </alias>', '</alias>')
lines = lines.replace('<zaid>\n ', '<zaid>')
lines = lines.replace('\n </zaid>', '</zaid>')
lines = lines.replace('<type>\n ', '<type>')
lines = lines.replace('\n </type>', '</type>')
lines = lines.replace('<awr>\n ', '<awr>')
lines = lines.replace('\n </awr>', '</awr>')
lines = lines.replace('<temperature>\n ', '<temperature>')
lines = lines.replace('\n </temperature>', '</temperature>')
lines = lines.replace('<path>\n ', '<path>')
lines = lines.replace('\n </path>', '</path>')
lines = lines.replace('<metastable>\n ', '<metastable>')
lines = lines.replace('\n </metastable>', '</metastable>')
lines = lines.replace('<binary>\n ', '<binary>')
lines = lines.replace('\n </binary>', '</binary>')
# Write document in pretty XML to specified file
f = open(xmlFile, 'w')

18
src/utils/convert_xsdir.py
View file

@ -16,6 +16,7 @@ elements = [None, "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na",
"Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg",
"Cn"]
class Xsdir(object):
def __init__(self, filename):
@ -54,14 +55,14 @@ class Xsdir(object):
# Check for end of second section
if len(words) % 2 != 0 or words[0] == 'directory':
break
for zaid, awr in zip(words[::2], words[1::2]):
self.awr[zaid] = awr
# Read third section
while words[0] != 'directory':
words = self.f.readline().split()
while True:
words = self.f.readline().split()
if not words:
@ -76,7 +77,7 @@ class Xsdir(object):
# Create XsdirTable object and add to line
table = XsdirTable(self.directory)
self.tables.append(table)
# All tables have at least 7 attributes
table.name = words[0]
table.awr = float(words[1])
@ -115,7 +116,6 @@ class Xsdir(object):
self.entries = list(self.entries)[0]
else:
self.recordlength = None
def to_xml(self):
# Create XML document
@ -154,7 +154,7 @@ class Xsdir(object):
# Add a node for each table
for table in self.tables:
if table.name[-1] in ['e', 'p', 'u', 'h', 'g' ,'m', 'd']:
if table.name[-1] in ['e', 'p', 'u', 'h', 'g', 'm', 'd']:
continue
node = table.to_xml_node(doc)
root.appendChild(node)
@ -230,7 +230,7 @@ class XsdirTable(object):
return "{0}-{1}.{2}".format(elements[Z], s, self.xs)
else:
return None
@property
def zaid(self):
if self.name.endswith('c'):
@ -245,10 +245,10 @@ class XsdirTable(object):
def to_xml_node(self, doc):
node = doc.createElement("ace_table")
node.setAttribute("name", self.name)
for attribute in ["alias", "zaid", "type", "metastable", "awr",
for attribute in ["alias", "zaid", "type", "metastable", "awr",
"temperature", "path", "location"]:
if hasattr(self, attribute):
string = str(getattr(self,attribute))
string = str(getattr(self, attribute))
# Skip metastable and binary if 0
if attribute == "metastable" and self.metastable == 0:
@ -274,7 +274,7 @@ if __name__ == '__main__':
# Read xsdata and create XML document object
xsdirObject = Xsdir(xsdirFile)
doc = xsdirObject.to_xml()
# Reduce number of lines
lines = doc.toprettyxml(indent=' ')

1
src/utils/geometry.py
View file

@ -2,6 +2,7 @@
from xml.dom.minidom import parse
class Geometry(object):
def __init__(self, filename):

4
src/utils/memory_usage.py
View file

@ -1,7 +1,7 @@
#!/usr/bin/env python2
# This script reads a cross_sections.out file, adds up the memory usage for each
# nuclide and S(a,b) table, and displays the total memory usage
# This script reads a cross_sections.out file, adds up the memory usage for
# each nuclide and S(a,b) table, and displays the total memory usage
import sys
import os

73
src/utils/particle_restart.py Normal file
View file

@ -0,0 +1,73 @@
#!/usr/bin/env python2
import struct
class Particle(object):
def __init__(self, filename):
if filename.endswith('.h5'):
import h5py
self._f = h5py.File(filename, 'r')
self._hdf5 = True
else:
self._f = open(filename, 'rb')
self._hdf5 = False
# Read all metadata
self._read_data()
def _read_data(self):
# Read filetype
self.filetype = self._get_int(path='filetype')[0]
# Read statepoint revision
self.revision = self._get_int(path='revision')[0]
# Read current batch
self.current_batch = self._get_int(path='current_batch')[0]
# Read run information
self.gen_per_batch = self._get_int(path='gen_per_batch')[0]
self.current_gen = self._get_int(path='current_gen')[0]
self.n_particles = self._get_long(path='n_particles')[0]
# Read particle properties
self.id = self._get_long(path='id')[0]
self.weight = self._get_double(path='weight')[0]
self.energy = self._get_double(path='energy')[0]
self.xyz = self._get_double(3, path='xyz')
self.uvw = self._get_double(3, path='uvw')
def _get_data(self, n, typeCode, size):
return list(struct.unpack('={0}{1}'.format(n, typeCode),
self._f.read(n*size)))
def _get_int(self, n=1, path=None):
if self._hdf5:
return [int(v) for v in self._f[path].value]
else:
return [int(v) for v in self._get_data(n, 'i', 4)]
def _get_long(self, n=1, path=None):
if self._hdf5:
return [long(v) for v in self._f[path].value]
else:
return [long(v) for v in self._get_data(n, 'q', 8)]
def _get_float(self, n=1, path=None):
if self._hdf5:
return [float(v) for v in self._f[path].value]
else:
return [float(v) for v in self._get_data(n, 'f', 4)]
def _get_double(self, n=1, path=None):
if self._hdf5:
return [float(v) for v in self._f[path].value]
else:
return [float(v) for v in self._get_data(n, 'd', 8)]
def _get_string(self, n=1, path=None):
if self._hdf5:
return str(self._f[path].value)
else:
return str(self._get_data(n, 's', 1)[0])

8
src/utils/split_nuclide.py
View file

@ -4,15 +4,13 @@ import sys
if len(sys.argv) != 3:
print("Must supply element and atom/b-cm")
exit
sys.exit(1)
element = sys.argv[1]
ao = float(sys.argv[2])
for line in open('abundances_modified.txt','r'):
for line in open('abundances_modified.txt', 'r'):
words = line.split()
if words[1] == element:
print('<nuclide name="{0}-{1}" ao="{2:10.4e}" />'.format(
element, words[2], float(words[3])*ao))
element, words[2], float(words[3])*ao))

9
src/utils/statepoint_histogram.py
View file

@ -56,14 +56,15 @@ for i, t in enumerate(sp.tallies):
relative_error = t_value*sqrt((s2/n - s*s)/(n-1))/s
uncertainties.append(relative_error)
# Display number of tallies with less than 1% CIs
fraction = (sum([1.0 if re < 0.01 else 0.0 for re in uncertainties])
/n_bins)
# Display number of tallies with less than 1% CIs
fraction = (sum([1.0 if re < 0.01 else 0.0 for re in uncertainties])
/ n_bins)
print("Tally " + str(i+1))
print(" Fraction under 1% = {0}".format(fraction))
print(" Min relative error = {0}".format(min(uncertainties)))
print(" Max relative error = {0}".format(max(uncertainties)))
print(" Non-scoring bins = {0}".format(1.0 - float(len(uncertainties))/n_bins))
print(" Non-scoring bins = {0}".format(
1.0 - float(len(uncertainties))/n_bins))
# Plot histogram
plt.hist(uncertainties, 100)

24
src/utils/statepoint_meshplot.py
View file

@ -60,9 +60,9 @@ for t in sp.tallies:
assert n_bins == nx*ny*nz*ns
# Create lists for tallies
mean = np.zeros((nx,ny))
error = np.zeros((nx,ny))
criteria = np.zeros((nx,ny))
mean = np.zeros((nx, ny))
error = np.zeros((nx, ny))
criteria = np.zeros((nx, ny))
# Determine starting position for data
start = sp._f.tell() + (axial_level-1)*ns*16 + (score - 1)*16
@ -75,19 +75,19 @@ for t in sp.tallies:
# Read sum and sum-squared
s, s2 = sp._get_double(2)
s /= n
mean[x,y] = s
mean[x, y] = s
if s != 0.0:
error[x,y] = t_value*sqrt((s2/n - s*s)/(n-1))/s
criteria[x,y] = 1.0 if error[x,y] < 0.05 else 0.0
error[x, y] = t_value*sqrt((s2/n - s*s)/(n-1))/s
criteria[x, y] = 1.0 if error[x, y] < 0.05 else 0.0
# Make figure
print("Making colorplot...")
plt.imshow(mean, interpolation="nearest")
plt.colorbar()
plt.xlim((0,nx))
plt.ylim((0,ny))
plt.xticks(np.linspace(0,nx,5),
np.linspace(m.lower_left[0],m.upper_right[0],5))
plt.yticks(np.linspace(0,ny,5),
np.linspace(m.lower_left[1],m.upper_right[1],5))
plt.xlim((0, nx))
plt.ylim((0, ny))
plt.xticks(np.linspace(0, nx, 5),
np.linspace(m.lower_left[0], m.upper_right[0], 5))
plt.yticks(np.linspace(0, ny, 5),
np.linspace(m.lower_left[1], m.upper_right[1], 5))
plt.show()

23
tests/convert_precision.py Executable file
View file

@ -0,0 +1,23 @@
#!/usr/bin/env python
import os
import glob
dirs = glob.glob('test_*')
for adir in dirs:
os.chdir(adir)
files = glob.glob('results.py')
if len(files) > 0:
files = files[0]
with open(files, 'r') as fh:
intxt = fh.read()
intxt = intxt.replace('14.8E', '12.6E')
with open(files, 'w') as fh:
fh.write(intxt)
os.chdir('..')

41
tests/nose_mpi.py Normal file
View file

@ -0,0 +1,41 @@
import os
import logging
from nose.plugins import Plugin
log = logging.getLogger('nose.plugins.nose_mpi')
class NoseMPI(Plugin):
mpi_np = "0"
mpi_exec = "/opt/mpich/3.0.4-gnu/bin/mpiexec"
def options(self, parser, env=os.environ):
"""Define the command line options for plugin."""
super(NoseMPI, self).options(parser, env)
parser.add_option(
"--mpi-np", dest="mpi_np", default=0,
help="Number of MPI processors to execute OpenMC executable.")
parser.add_option(
"--mpi-exec", dest="mpi_exec",
default="/opt/mpich/3.0.4-gnu/bin/mpiexec",
help="Absolute path to mpiexec file.")
def configure(self, options, conf):
"""Configure plugin based on command line options"""
super(NoseMPI, self).configure(options, conf)
try:
mpi_np = int(options.mpi_np)
self.enabled = True
NoseMPI.mpi_np = options.mpi_np
except:
self.enabled = False
return
if not os.path.exists(options.mpi_exec):
print 'Need to specify valid mpiexec path.'
exit()
NoseMPI.mpi_exec = options.mpi_exec

101
tests/run_tests.py Executable file
View file

@ -0,0 +1,101 @@
#!/usr/bin/env python
from __future__ import print_function
import os
import sys
import nose
import glob
from subprocess import call
from nose_mpi import NoseMPI
def run_compile():
print('-'*17)
print('Compilation tests')
print('-'*17)
# clean up all previous executables
openmc_exe = glob.glob(pwd + '/../src/openmc*')
for exe in openmc_exe:
os.remove(exe)
# run compile test
result = nose.run(argv=['run_tests.py', 'test_compile'] + flags)
if not result:
print('Did not pass compile tests.')
results.append(('compile', result))
def run_suite(name=None, mpi=False):
print('-'*(len(name) + 6))
print(name + ' tests')
print('-'*(len(name) + 6))
# Set arguments list. Note that the first argument is a dummy argument (the
# script name). It's not actually recursively calling run_tests.py
argv = ['run_tests.py', '--exclude', 'test_compile'] + flags
# Add MPI plugin if set
if mpi:
plugins = [NoseMPI()]
argv += ['--mpi-np', '3', '--mpi-exec', mpiexec]
else:
plugins = None
try:
os.chdir(pwd)
os.rename(pwd + '/../src/openmc-' + name, pwd + '/../src/openmc')
result = nose.run(argv=argv, addplugins=plugins)
finally:
os.rename(pwd + '/../src/openmc', pwd + '/../src/openmc-' + name)
if not result:
print('Did not pass ' + name + ' tests')
results.append((name, result))
# set mpiexec path
if os.environ.has_key('COMPILER'):
compiler = os.environ['COMPILER']
else:
compiler = 'gnu'
mpiexec = '/opt/mpich/3.0.4-{0}/bin/mpiexec'.format(compiler)
# get current working directory
pwd = os.getcwd()
sys.path.append(pwd)
# Set list of tests, either default or from command line
flags = []
tests = ['compile', 'normal', 'debug', 'optimize', 'hdf5', 'mpi', 'phdf5',
'petsc', 'phdf5-petsc', 'phdf5-petsc-optimize']
if len(sys.argv) > 1:
flags = [i for i in sys.argv[1:] if i.startswith('-')]
tests_ = [i for i in sys.argv[1:] if not i.startswith('-')]
tests = tests_ if tests_ else tests
# Run tests
results = []
for name in tests:
if name == 'compile':
run_compile()
elif name in ['normal', 'debug', 'optimize', 'hdf5']:
run_suite(name=name)
elif name in ['mpi', 'phdf5', 'petsc', 'phdf5-petsc',
'phdf5-petsc-optimize']:
run_suite(name=name, mpi=True)
# print out summary of results
print('\n' + '='*54)
print('Summary of Compilation Option Testing:\n')
OK = '\033[92m'
FAIL = '\033[91m'
ENDC = '\033[0m'
BOLD = '\033[1m'
for name, result in results:
print(name + '.'*(50 - len(name)), end='')
if result:
print(BOLD + OK + '[OK]' + ENDC)
else:
print(BOLD + FAIL + '[FAILED]' + ENDC)

25
tests/test_basic/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

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

@ -0,0 +1,2 @@
k-combined:
3.011353E-01 2.854556E-03

27
tests/test_basic/test_basic.py
View file

@ -1,7 +1,10 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
@ -10,16 +13,32 @@ def setup():
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
assert os.path.exists(pwd + '/statepoint.10.binary')
statepoint = glob.glob(pwd + '/statepoint.10.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_results():
statepoint = glob.glob(pwd + '/statepoint.10.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.10.binary']
output = glob.glob(pwd + '/statepoint.10.*')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

16
tests/test_cmfd_feed/cmfd.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<cmfd>
<mesh>
<lower_left> -10 -1 -1 </lower_left>
<upper_right> 10 1 1 </upper_right>
<dimension> 10 1 1 </dimension>
<albedo> 0.0 0.0 1.0 1.0 1.0 1.0 </albedo>
</mesh>
<begin> 5 </begin>
<display> dominance </display>
<solver> power </solver>
<feedback> true </feedback>
</cmfd>

43
tests/test_cmfd_feed/geometry.xml Normal file
View file

@ -0,0 +1,43 @@
<?xml version="1.0" encoding="UTF-8"?>
<geometry>
<!-- Definition of Cells -->
<cell id="1">
<universe>0</universe>
<surfaces>-1 2 -3 4 -5 6</surfaces>
<material>1</material>
</cell>
<!-- Defition of Surfaces -->
<surface id="1">
<type>x-plane</type>
<coeffs>10</coeffs>
<boundary> vacuum </boundary>
</surface>
<surface id="2">
<type>x-plane</type>
<coeffs>-10</coeffs>
<boundary> vacuum </boundary>
</surface>
<surface id="3">
<type>y-plane</type>
<coeffs>1</coeffs>
<boundary>reflective</boundary>
</surface>
<surface id="4">
<type>y-plane</type>
<coeffs>-1</coeffs>
<boundary>reflective</boundary>
</surface>
<surface id="5">
<type>z-plane</type>
<coeffs>1</coeffs>
<boundary>reflective</boundary>
</surface>
<surface id="6">
<type>z-plane</type>
<coeffs>-1</coeffs>
<boundary>reflective</boundary>
</surface>
</geometry>

12
tests/test_cmfd_feed/materials.xml Normal file
View file

@ -0,0 +1,12 @@
<?xml version="1.0"?>
<materials>
<!-- Definition of materials -->
<material id="1">
<density value="19" units="g/cc" />
<nuclide name="U-235" xs="70c" wo="0.21" />
<nuclide name="U-238" xs="70c" wo="0.68" />
<nuclide name="O-16" xs="70c" wo="0.11" />
</material>
</materials>

58
tests/test_cmfd_feed/results.py Normal file
View file

@ -0,0 +1,58 @@
#!/usr/bin/env python
import sys
import numpy as np
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.20.binary')
sp.read_results()
# extract tally results and convert to vector
results1 = sp.tallies[0].results
shape1 = results1.shape
size1 = (np.product(shape1))
results1 = np.reshape(results1, size1)
results2 = sp.tallies[1].results
shape2 = results2.shape
size2 = (np.product(shape2))
results2 = np.reshape(results2, size2)
results3 = sp.tallies[2].results
shape3 = results3.shape
size3 = (np.product(shape3))
results3 = np.reshape(results3, size3)
results4 = sp.tallies[3].results
shape4 = results4.shape
size4 = (np.product(shape4))
results4 = np.reshape(results4, size4)
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write out tally results
outstr += 'tally 1:\n'
for item in results1:
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 2:\n'
for item in results2:
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 3:\n'
for item in results3:
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 4:\n'
for item in results4:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

654
tests/test_cmfd_feed/results_true.dat Normal file
View file

@ -0,0 +1,654 @@
k-combined:
1.167124E+00 1.217343E-02
tally 1:
1.126891E+01
1.275755E+01
2.086556E+01
4.369501E+01
2.853629E+01
8.169559E+01
3.404506E+01
1.165019E+02
3.723908E+01
1.389881E+02
3.760048E+01
1.416691E+02
3.455469E+01
1.197186E+02
2.837503E+01
8.086291E+01
2.151789E+01
4.648315E+01
1.194541E+01
1.432331E+01
tally 2:
2.292173E+01
2.650792E+01
1.598183E+01
1.288692E+01
2.194143E+00
2.510557E-01
4.217816E+01
8.938583E+01
2.974466E+01
4.447038E+01
3.939269E+00
7.882413E-01
5.742289E+01
1.654599E+02
4.067581E+01
8.308454E+01
5.566955E+00
1.570072E+00
6.810586E+01
2.331019E+02
4.831543E+01
1.174673E+02
6.244999E+00
1.967938E+00
7.258131E+01
2.646756E+02
5.185146E+01
1.351425E+02
6.622488E+00
2.231695E+00
7.246115E+01
2.641295E+02
5.159542E+01
1.339441E+02
6.702239E+00
2.277825E+00
6.716397E+01
2.266815E+02
4.765063E+01
1.141627E+02
6.384098E+00
2.063760E+00
5.549542E+01
1.545361E+02
3.940819E+01
7.794880E+01
5.183375E+00
1.369483E+00
4.155685E+01
8.679439E+01
2.926744E+01
4.307145E+01
3.938852E+00
7.877949E-01
2.335180E+01
2.758241E+01
1.622683E+01
1.331996E+01
2.248921E+00
2.594934E-01
tally 3:
1.537807E+01
1.193746E+01
1.044189E+00
5.627797E-02
2.861614E+01
4.118322E+01
1.846208E+00
1.729007E-01
3.921462E+01
7.723576E+01
2.427723E+00
2.983153E-01
4.653207E+01
1.089818E+02
3.128394E+00
4.945409E-01
4.989896E+01
1.251503E+02
3.242191E+00
5.330743E-01
4.967378E+01
1.241657E+02
3.184751E+00
5.164524E-01
4.592989E+01
1.061019E+02
2.936424E+00
4.356307E-01
3.783939E+01
7.188489E+01
2.485680E+00
3.117618E-01
2.819358E+01
3.997588E+01
1.875908E+00
1.772000E-01
1.562421E+01
1.236020E+01
1.013140E+00
5.260420E-02
tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.126365E+00
4.930682E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.456410E+00
1.499416E+00
2.727755E+00
3.788867E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.453480E+00
2.790159E+00
5.240103E+00
1.386895E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.741563E+00
3.830225E+00
7.137968E+00
2.555042E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.221731E+00
4.266573E+00
8.474116E+00
3.604647E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.054937E+00
4.124537E+00
9.066119E+00
4.133833E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.362253E+00
3.517245E+00
9.086548E+00
4.156118E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.871573E+00
2.369190E+00
8.399584E+00
3.549784E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.998103E+00
1.258194E+00
7.145296E+00
2.568115E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.716076E+00
3.753568E-01
5.427593E+00
1.484154E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.056645E+00
4.702349E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00

32
tests/test_cmfd_feed/settings.xml Normal file
View file

@ -0,0 +1,32 @@
<?xml version="1.0" encoding="UTF-8"?>
<settings>
<!-- Parameters for criticality calculation -->
<eigenvalue>
<batches>20</batches>
<inactive>10</inactive>
<particles>1000</particles>
</eigenvalue>
<!-- How verbose output should be -->
<verbosity value="7" />
<!-- Starting source -->
<source>
<space>
<type>box</type>
<parameters>-10 -1 -1 10 1 1</parameters>
</space>
</source>
<!-- Shannon Entropy -->
<entropy>
<dimension> 10 1 1 </dimension>
<lower_left> -10.0 -1.0 -1.0 </lower_left>
<upper_right> 10.0 1.0 1.0 </upper_right>
</entropy>
<!-- Run CMFD -->
<run_cmfd> true </run_cmfd>
</settings>

16
tests/test_cmfd_feed/tallies.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1">
<type>rectangular</type>
<lower_left>-10 -1 -1 </lower_left>
<upper_right>10 1 1</upper_right>
<dimension>10 1 1</dimension>
</mesh>
<tally id="1">
<filter type="mesh" bins="1" />
<scores>flux</scores>
</tally>
</tallies>

68
tests/test_cmfd_feed/test_cmfd_feed.py Normal file
View file

@ -0,0 +1,68 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from nose.plugins.skip import SkipTest
from nose_mpi import NoseMPI
pwd = os.path.dirname(__file__)
skipAll = False
def setup():
os.putenv('PWD', pwd)
os.chdir(pwd)
global skipAll
skipAll = False
def test_run():
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
output = proc.communicate()[0]
print(output)
if 'CMFD is not available' in output:
global skipAll
skipAll = True
raise SkipTest
assert returncode == 0
def test_created_statepoint():
if skipAll:
raise SkipTest
statepoint = glob.glob(pwd + '/statepoint.20.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_output_exists():
if skipAll:
raise SkipTest
assert os.path.exists(pwd + '/tallies.out')
def test_results():
if skipAll:
raise SkipTest
statepoint = glob.glob(pwd + '/statepoint.20.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = glob.glob(pwd + '/statepoint.20.*')
output.append(pwd + '/tallies.out')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

16
tests/test_cmfd_nofeed/cmfd.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<cmfd>
<mesh>
<lower_left> -10 -1 -1 </lower_left>
<upper_right> 10 1 1 </upper_right>
<dimension> 10 1 1 </dimension>
<albedo> 0.0 0.0 1.0 1.0 1.0 1.0 </albedo>
</mesh>
<begin> 5 </begin>
<display> dominance </display>
<solver> power </solver>
<feedback> false </feedback>
</cmfd>

43
tests/test_cmfd_nofeed/geometry.xml Normal file
View file

@ -0,0 +1,43 @@
<?xml version="1.0" encoding="UTF-8"?>
<geometry>
<!-- Definition of Cells -->
<cell id="1">
<universe>0</universe>
<surfaces>-1 2 -3 4 -5 6</surfaces>
<material>1</material>
</cell>
<!-- Defition of Surfaces -->
<surface id="1">
<type>x-plane</type>
<coeffs>10</coeffs>
<boundary> vacuum </boundary>
</surface>
<surface id="2">
<type>x-plane</type>
<coeffs>-10</coeffs>
<boundary> vacuum </boundary>
</surface>
<surface id="3">
<type>y-plane</type>
<coeffs>1</coeffs>
<boundary>reflective</boundary>
</surface>
<surface id="4">
<type>y-plane</type>
<coeffs>-1</coeffs>
<boundary>reflective</boundary>
</surface>
<surface id="5">
<type>z-plane</type>
<coeffs>1</coeffs>
<boundary>reflective</boundary>
</surface>
<surface id="6">
<type>z-plane</type>
<coeffs>-1</coeffs>
<boundary>reflective</boundary>
</surface>
</geometry>

12
tests/test_cmfd_nofeed/materials.xml Normal file
View file

@ -0,0 +1,12 @@
<?xml version="1.0"?>
<materials>
<!-- Definition of materials -->
<material id="1">
<density value="19" units="g/cc" />
<nuclide name="U-235" xs="70c" wo="0.21" />
<nuclide name="U-238" xs="70c" wo="0.68" />
<nuclide name="O-16" xs="70c" wo="0.11" />
</material>
</materials>

58
tests/test_cmfd_nofeed/results.py Normal file
View file

@ -0,0 +1,58 @@
#!/usr/bin/env python
import sys
import numpy as np
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.20.binary')
sp.read_results()
# extract tally results and convert to vector
results1 = sp.tallies[0].results
shape1 = results1.shape
size1 = (np.product(shape1))
results1 = np.reshape(results1, size1)
results2 = sp.tallies[1].results
shape2 = results2.shape
size2 = (np.product(shape2))
results2 = np.reshape(results2, size2)
results3 = sp.tallies[2].results
shape3 = results3.shape
size3 = (np.product(shape3))
results3 = np.reshape(results3, size3)
results4 = sp.tallies[3].results
shape4 = results4.shape
size4 = (np.product(shape4))
results4 = np.reshape(results4, size4)
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write out tally results
outstr += 'tally 1:\n'
for item in results1:
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 2:\n'
for item in results2:
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 3:\n'
for item in results3:
outstr += "{0:12.6E}\n".format(item)
outstr += 'tally 4:\n'
for item in results4:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

654
tests/test_cmfd_nofeed/results_true.dat Normal file
View file

@ -0,0 +1,654 @@
k-combined:
1.137460E+00 1.215114E-02
tally 1:
1.168490E+01
1.400522E+01
2.192157E+01
4.824543E+01
2.936949E+01
8.648246E+01
3.435607E+01
1.185184E+02
3.612841E+01
1.309721E+02
3.552232E+01
1.266988E+02
3.314322E+01
1.107274E+02
2.765588E+01
7.705156E+01
2.044445E+01
4.210167E+01
1.125932E+01
1.276689E+01
tally 2:
2.372613E+01
2.862691E+01
1.660400E+01
1.404387E+01
2.237703E+00
2.562122E-01
4.397041E+01
9.722086E+01
3.107800E+01
4.859826E+01
4.052018E+00
8.364586E-01
5.850026E+01
1.715533E+02
4.157400E+01
8.666187E+01
5.420325E+00
1.492090E+00
6.875182E+01
2.377104E+02
4.875600E+01
1.196878E+02
6.167355E+00
1.922001E+00
7.240156E+01
2.639189E+02
5.161800E+01
1.342800E+02
6.657325E+00
2.248300E+00
7.074222E+01
2.516854E+02
5.016500E+01
1.266005E+02
6.598272E+00
2.224189E+00
6.625191E+01
2.209754E+02
4.694400E+01
1.110192E+02
6.197674E+00
1.947375E+00
5.574790E+01
1.561211E+02
3.963100E+01
7.893339E+01
5.104628E+00
1.314367E+00
4.068463E+01
8.316576E+01
2.878400E+01
4.165629E+01
4.000895E+00
8.100229E-01
2.238335E+01
2.532998E+01
1.555300E+01
1.221694E+01
2.312206E+00
2.751529E-01
tally 3:
1.599200E+01
1.304060E+01
1.075546E+00
5.989825E-02
2.991600E+01
4.504778E+01
2.016558E+00
2.060921E-01
4.008500E+01
8.057644E+01
2.465241E+00
3.073562E-01
4.690900E+01
1.108192E+02
3.050094E+00
4.696181E-01
4.972600E+01
1.246434E+02
3.127136E+00
4.939882E-01
4.836400E+01
1.176964E+02
3.102574E+00
4.868820E-01
4.524800E+01
1.031819E+02
2.966738E+00
4.482948E-01
3.809700E+01
7.295318E+01
2.526411E+00
3.226394E-01
2.772600E+01
3.866039E+01
1.852304E+00
1.725388E-01
1.494500E+01
1.129101E+01
9.995136E-01
5.115569E-02
tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.251000E+00
5.355770E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.759000E+00
1.676451E+00
2.883000E+00
4.269210E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.633000E+00
2.922951E+00
5.354000E+00
1.442474E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.787000E+00
3.871605E+00
7.137000E+00
2.553987E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.343000E+00
4.385923E+00
8.470000E+00
3.603838E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.076000E+00
4.142704E+00
8.977000E+00
4.056459E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.252000E+00
3.428360E+00
8.957000E+00
4.041545E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.868000E+00
2.375736E+00
8.404000E+00
3.562144E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.914000E+00
1.211724E+00
7.074000E+00
2.511368E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.634000E+00
3.523220E-01
5.349000E+00
1.437693E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.050000E+00
4.687600E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00

32
tests/test_cmfd_nofeed/settings.xml Normal file
View file

@ -0,0 +1,32 @@
<?xml version="1.0" encoding="UTF-8"?>
<settings>
<!-- Parameters for criticality calculation -->
<eigenvalue>
<batches>20</batches>
<inactive>10</inactive>
<particles>1000</particles>
</eigenvalue>
<!-- How verbose output should be -->
<verbosity value="7" />
<!-- Starting source -->
<source>
<space>
<type>box</type>
<parameters>-10 -1 -1 10 1 1</parameters>
</space>
</source>
<!-- Shannon Entropy -->
<entropy>
<dimension> 10 1 1 </dimension>
<lower_left> -10.0 -1.0 -1.0 </lower_left>
<upper_right> 10.0 1.0 1.0 </upper_right>
</entropy>
<!-- Run CMFD -->
<run_cmfd> true </run_cmfd>
</settings>

16
tests/test_cmfd_nofeed/tallies.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1">
<type>rectangular</type>
<lower_left>-10 -1 -1 </lower_left>
<upper_right>10 1 1</upper_right>
<dimension>10 1 1</dimension>
</mesh>
<tally id="1">
<filter type="mesh" bins="1" />
<scores>flux</scores>
</tally>
</tallies>

67
tests/test_cmfd_nofeed/test_cmfd_nofeed.py Normal file
View file

@ -0,0 +1,67 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from nose.plugins.skip import SkipTest
from nose_mpi import NoseMPI
pwd = os.path.dirname(__file__)
skipAll = False
def setup():
os.putenv('PWD', pwd)
os.chdir(pwd)
global skipAll
skipAll = False
def test_run():
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
output = proc.communicate()[0]
print(output)
if 'CMFD is not available' in output:
global skipAll
skipAll = True
raise SkipTest
assert returncode == 0
def test_created_statepoint():
if skipAll:
raise SkipTest
statepoint = glob.glob(pwd + '/statepoint.20.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_output_exists():
if skipAll:
raise SkipTest
assert os.path.exists(pwd + '/tallies.out')
def test_results():
if skipAll:
raise SkipTest
statepoint = glob.glob(pwd + '/statepoint.20.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = glob.glob(pwd + '/statepoint.20.*')
output.append(pwd + '/tallies.out')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

103
tests/test_compile/test_compile.py
View file

@ -15,100 +15,73 @@ import shutil
pwd = os.path.dirname(__file__)
if os.environ.has_key('COMPILER'):
compiler = 'COMPILER=' + os.environ['COMPILER']
else:
compiler = 'COMPILER=gnu'
def setup():
# Change to source directory
os.chdir(pwd + '/../../src')
def test_gfortran():
def test_normal():
returncode = run(['make','distclean'])
returncode = run(['make'])
returncode = run(['make', compiler])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran')
shutil.move('openmc', 'openmc-normal')
def test_gfortran_debug():
def test_debug():
returncode = run(['make','distclean'])
returncode = run(['make', 'DEBUG=yes'])
returncode = run(['make', compiler, 'DEBUG=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-debug')
def test_profile():
returncode = run(['make','distclean'])
returncode = run(['make', compiler, 'PROFILE=yes'])
assert returncode == 0
def test_gfortran_profile():
def test_optimize():
returncode = run(['make','distclean'])
returncode = run(['make', 'PROFILE=yes'])
returncode = run(['make', compiler, 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-optimize')
def test_gfortran_optimize():
def test_mpi():
returncode = run(['make','distclean'])
returncode = run(['make', 'OPTIMIZE=yes'])
returncode = run(['make', compiler, 'MPI=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-mpi')
def test_gfortran_mpi():
def test_hdf5():
returncode = run(['make','distclean'])
returncode = run(['make', 'MPI=yes'])
returncode = run(['make', compiler, 'HDF5=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-hdf5')
def test_gfortran_hdf5():
def test_petsc():
returncode = run(['make','distclean'])
returncode = run(['make', 'HDF5=yes'])
returncode = run(['make', compiler, 'MPI=yes', 'PETSC=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-petsc')
def test_gfortran_petsc():
def test_mpi_hdf5():
returncode = run(['make','distclean'])
returncode = run(['make', 'MPI=yes', 'PETSC=yes'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5')
def test_gfortran_mpi_hdf5():
def test_mpi_hdf5_petsc():
returncode = run(['make','distclean'])
returncode = run(['make', 'MPI=yes', 'HDF5=yes'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-petsc')
def test_gfortran_mpi_hdf5_petsc():
def test_mpi_hdf5_petsc_optimize():
returncode = run(['make','distclean'])
returncode = run(['make', 'MPI=yes', 'HDF5=yes', 'PETSC=yes'])
assert returncode == 0
def test_intel():
returncode = run(['make','distclean'])
returncode = run(['make', 'COMPILER=intel'])
assert returncode == 0
def test_intel_debug():
returncode = run(['make','distclean'])
returncode = run(['make', 'COMPILER=intel', 'DEBUG=yes'])
assert returncode == 0
def test_intel_profile():
returncode = run(['make','distclean'])
returncode = run(['make', 'COMPILER=intel', 'PROFILE=yes'])
assert returncode == 0
def test_intel_optimize():
returncode = run(['make','distclean'])
returncode = run(['make', 'COMPILER=intel', 'OPTIMIZE=yes'])
assert returncode == 0
def test_intel_mpi():
returncode = run(['make','distclean'])
returncode = run(['make', 'COMPILER=intel', 'MPI=yes'])
assert returncode == 0
def test_intel_hdf5():
returncode = run(['make','distclean'])
returncode = run(['make', 'COMPILER=intel', 'HDF5=yes'])
assert returncode == 0
def test_intel_petsc():
returncode = run(['make','distclean'])
returncode = run(['make', 'COMPILER=intel', 'MPI=yes', 'PETSC=yes'])
assert returncode == 0
def test_intel_mpi_hdf5():
returncode = run(['make','distclean'])
returncode = run(['make', 'COMPILER=intel', 'MPI=yes', 'HDF5=yes'])
assert returncode == 0
def test_intel_mpi_hdf5_petsc():
returncode = run(['make','distclean'])
returncode = run(['make', 'COMPILER=intel', 'MPI=yes', 'HDF5=yes', 'PETSC=yes'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes', 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-petsc-optimize')
def run(commands):
proc = Popen(commands, stderr=STDOUT, stdout=PIPE)
@ -118,4 +91,4 @@ def run(commands):
def teardown(commands):
returncode = run(['make','distclean'])
shutil.copy('openmc-gfortran','openmc')
shutil.copy('openmc-normal','openmc')

37
tests/test_confidence_intervals/results.py Normal file
View file

@ -0,0 +1,37 @@
#!/usr/bin/env python
import sys
import numpy as np
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write out tally results
outstr += 'tallies:\n'
for item in results:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

5
tests/test_confidence_intervals/results_true.dat Normal file
View file

@ -0,0 +1,5 @@
k-combined:
2.937103E-01 7.391597E-03
tallies:
6.267771E+01
4.928030E+02

28
tests/test_confidence_intervals/test_confidence_intervals.py
View file

@ -1,7 +1,10 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
@ -10,19 +13,36 @@ def setup():
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
assert os.path.exists(pwd + '/statepoint.10.binary')
statepoint = glob.glob(pwd + '/statepoint.10.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_created_output():
assert os.path.exists(pwd + '/tallies.out')
def test_results():
statepoint = glob.glob(pwd + '/statepoint.10.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.10.binary', pwd + '/tallies.out']
output = glob.glob(pwd + '/statepoint.10.*')
output.append(pwd + '/results_test.dat')
output.append(pwd + '/tallies.out')
for f in output:
if os.path.exists(f):
os.remove(f)

25
tests/test_density_atombcm/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

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

@ -0,0 +1,2 @@
k-combined:
1.739614E+00 1.281892E-02

30
tests/test_density_atombcm/test_density_atombcm.py
View file

@ -1,22 +1,44 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
def setup():
def setup():
os.putenv('PWD', pwd)
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
statepoint = glob.glob(pwd + '/statepoint.10.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_results():
statepoint = glob.glob(pwd + '/statepoint.10.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.10.binary']
output = glob.glob(pwd + '/statepoint.10.*')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

25
tests/test_density_atomcm3/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

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

@ -0,0 +1,2 @@
k-combined:
1.081464E+00 2.283652E-02

28
tests/test_density_atomcm3/test_density_atomcm3.py
View file

@ -1,7 +1,10 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
@ -10,13 +13,32 @@ def setup():
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
statepoint = glob.glob(pwd + '/statepoint.10.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_results():
statepoint = glob.glob(pwd + '/statepoint.10.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.10.binary']
output = glob.glob('statepoint.10.*')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

25
tests/test_density_kgm3/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

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

@ -0,0 +1,2 @@
k-combined:
8.058880E-01 5.121476E-03

28
tests/test_density_kgm3/test_density_kgm3.py
View file

@ -1,7 +1,10 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
@ -10,13 +13,32 @@ def setup():
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
statepoint = glob.glob(pwd + '/statepoint.10.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_results():
statepoint = glob.glob(pwd + '/statepoint.10.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.10.binary']
output = glob.glob(pwd + '/statepoint.10.*')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

25
tests/test_density_sum/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

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

@ -0,0 +1,2 @@
k-combined:
3.166569E-01 6.367411E-03

27
tests/test_density_sum/test_density_sum.py
View file

@ -1,7 +1,10 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
@ -10,16 +13,32 @@ def setup():
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
assert os.path.exists(pwd + '/statepoint.10.binary')
statepoint = glob.glob(pwd + '/statepoint.10.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_results():
statepoint = glob.glob('statepoint.10.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.10.binary']
output = glob.glob(pwd + '/statepoint.10.*')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

25
tests/test_eigenvalue_genperbatch/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.7.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

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

@ -0,0 +1,2 @@
k-combined:
3.004375E-01 2.308696E-03

27
tests/test_eigenvalue_genperbatch/test_eigenvalue_genperbatch.py
View file

@ -1,7 +1,10 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
@ -10,16 +13,32 @@ def setup():
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
assert os.path.exists(pwd + '/statepoint.7.binary')
statepoint = glob.glob(pwd + '/statepoint.7.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_results():
statepoint = glob.glob(pwd + '/statepoint.7.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.7.binary']
output = glob.glob(pwd + '/statepoint.7.*')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

25
tests/test_eigenvalue_no_inactive/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

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

@ -0,0 +1,2 @@
k-combined:
3.118508E-01 6.764530E-03

27
tests/test_eigenvalue_no_inactive/test_eigenvalue_no_inactive.py
View file

@ -1,7 +1,10 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
@ -10,16 +13,32 @@ def setup():
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
assert os.path.exists(pwd + '/statepoint.10.binary')
statepoint = glob.glob(pwd + '/statepoint.10.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_results():
statepoint = glob.glob(pwd + '/statepoint.10.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.10.binary']
output = glob.glob(pwd + '/statepoint.10.*')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

25
tests/test_energy_grid/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

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

@ -0,0 +1,2 @@
k-combined:
3.243018E-01 9.673636E-04

28
tests/test_energy_grid/test_energy_grid.py
View file

@ -1,7 +1,10 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
@ -10,13 +13,32 @@ def setup():
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
statepoint = glob.glob(pwd + '/statepoint.10.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_results():
statepoint = glob.glob(pwd + '/statepoint.10.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.10.binary']
output = glob.glob(pwd + '/statepoint.10.*')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

31
tests/test_entropy/results.py Normal file
View file

@ -0,0 +1,31 @@
#!/usr/bin/env python
import sys
import numpy as np
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write out tally results
outstr += 'entropy:\n'
for item in sp.entropy:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

13
tests/test_entropy/results_true.dat Normal file
View file

@ -0,0 +1,13 @@
k-combined:
3.011353E-01 2.854556E-03
entropy:
7.114513E+00
8.055602E+00
8.233973E+00
8.235572E+00
8.284484E+00
8.289848E+00
8.420080E+00
8.279431E+00
8.299997E+00
8.363968E+00

27
tests/test_entropy/test_entropy.py
View file

@ -1,7 +1,10 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
from nose_mpi import NoseMPI
import glob
pwd = os.path.dirname(__file__)
@ -10,16 +13,32 @@ def setup():
os.chdir(pwd)
def test_run():
proc = Popen([pwd + '/../../src/openmc'], stderr=STDOUT, stdout=PIPE)
openmc_path = pwd + '/../../src/openmc'
if int(NoseMPI.mpi_np) > 0:
proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
print(proc.communicate()[0])
assert returncode == 0
def test_created_statepoint():
assert os.path.exists(pwd + '/statepoint.10.binary')
statepoint = glob.glob(pwd + '/statepoint.10.*')
assert len(statepoint) == 1
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
def test_results():
statepoint = glob.glob(pwd + '/statepoint.10.*')
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
output = [pwd + '/statepoint.10.binary']
output = glob.glob(pwd + '/statepoint.10.*')
output.append(pwd + '/results_test.dat')
for f in output:
if os.path.exists(f):
os.remove(f)

37
tests/test_filter_cell/results.py Normal file
View file

@ -0,0 +1,37 @@
#!/usr/bin/env python
import sys
import numpy as np
# import statepoint
sys.path.append('../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write out tally results
outstr += 'tallies:\n'
for item in results:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

11
tests/test_filter_cell/results_true.dat Normal file
View file

@ -0,0 +1,11 @@
k-combined:
9.377383E-01 6.215432E-02
tallies:
0.000000E+00
0.000000E+00
1.882286E+01
7.203128E+01
4.035920E+00
3.309607E+00
5.615336E+01
6.414295E+02

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