merged 'develop' into 'xml'

This commit is contained in:
Bryan Herman 2013-09-25 14:43:32 -04:00
commit 02d7679b89
35 changed files with 4092 additions and 1696 deletions

View file

@ -35,17 +35,17 @@ Don't use ``print *`` or ``write(*,*)``. If writing to a file, use a specific
unit. Writing to standard output or standard error should be handled by the
``write_message`` subroutine or functionality in the error module.
-------------------------
Subroutines and Functions
-------------------------
----------
Procedures
----------
Above each subroutine/function, include a comment block giving a brief
description of what the subroutine or function does.
Above each procedure, include a comment block giving a brief description of what
the procedure does.
Arguments to subroutines/functions should be explicitly specified as intent(in),
intent(out), or intent(inout).
Nonpointer dummy arguments to procedures should be explicitly specified as
intent(in), intent(out), or intent(inout).
Include a comment describing what each argument to a subroutine/function is.
Include a comment describing what each argument to a procedure is.
---------
Variables
@ -133,9 +133,11 @@ value of f90-continuation-indent in Emacs.
Whitespace in Expressions
-------------------------
Use a single space between arguments to procedures.
Avoid extraneous whitespace in the following situations:
- In subroutine/function calls::
- In procedure calls::
Yes: call somesub(x, y(2), z)
No: call somesub( x, y( 2 ), z )

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

View file

@ -388,6 +388,14 @@ survival biasing, otherwise known as implicit capture or absorption.
.. _trace:
``<threads>`` Element
---------------------
The ``<threads>`` element indicates the number of OpenMP threads to be used for
a simulation. It has no attributes and accepts a positive integer value.
*Default*: None (Determined by environment variable :envvar:`OMP_NUM_THREADS`)
``<trace>`` Element
-------------------

View file

@ -141,6 +141,10 @@ MPI
Enables parallel runs using the Message Passing Interface. The MPI_DIR
variable should be set to the base directory of the MPI implementation.
OPENMP
Enables shared-memory parallelism using the OpenMP API. The Fortran compiler
being used must support OpenMP.
HDF5
Enables HDF5 output in addition to normal screen and text file output. The
HDF5_DIR variable should be set to the base directory of the HDF5
@ -326,6 +330,20 @@ Alternatively, you could run from any directory:
Note that in the latter case, any output files will be placed in the present
working directory which may be different from ``/home/username/somemodel``.
Command-Line Flags
------------------
OpenMC accepts the following command line flags:
-g, --geometry-debug Run in geometry debugging mode, where cell overlaps are
checked for after each move of a particle
-n, --particles N Use *N* particles per generation or batch
-p, --plot Run in plotting mode
-r, --restart file Restart a previous run from a state point or a particle
restart file
-s, --threads N Run with *N* OpenMP threads
-v, --version Show version information
-----------------------------------------------------
Configuring Input Validation with GNU Emacs nXML mode
-----------------------------------------------------

View file

@ -24,6 +24,9 @@ Run in plotting mode
Restart a previous run from a state point or a particle restart file named
\fIbinaryFile\fP.
.TP
.BI \-s " N" "\fR,\fP \-\-threads" " N"
Use \fIN\fP OpenMP threads.
.TP
.B "\-v\fR, \fP\-\-version"
Show version information.
.TP

View file

@ -167,7 +167,6 @@ 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
@ -198,6 +197,7 @@ initialize.o: global.o
initialize.o: hdf5_interface.o
initialize.o: hdf5_summary.o
initialize.o: input_xml.o
initialize.o: omp_lib.o
initialize.o: output.o
initialize.o: output_interface.o
initialize.o: random_lcg.o
@ -278,21 +278,20 @@ 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: output_interface.o
particle_restart.o: particle_header.o
particle_restart.o: random_lcg.o
particle_restart.o: tracking.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: output_interface.o
particle_restart_write.o: particle_header.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
@ -381,9 +380,11 @@ tracking.o: global.o
tracking.o: output.o
tracking.o: particle_header.o
tracking.o: physics.o
tracking.o: random_lcg.o
tracking.o: string.o
tracking.o: tally.o
xml_interface.o: constants.o
xml_interface.o: error.o
xml_interface.o: global.o

View file

@ -14,6 +14,7 @@ DEBUG = no
PROFILE = no
OPTIMIZE = no
MPI = no
OPENMP = no
HDF5 = no
PETSC = no
@ -176,6 +177,25 @@ else
endif
endif
# OpenMP for shared-memory parallelism
ifeq ($(OPENMP),yes)
ifeq ($(COMPILER),intel)
F90FLAGS += -openmp -DOPENMP
LDFLAGS += -openmp
endif
ifeq ($(COMPILER),gnu)
F90FLAGS += -fopenmp -DOPENMP
LDFLAGS += -fopenmp
endif
ifeq ($(COMPILER),ibm)
F90FLAGS += -qsmp=omp -WF,-DOPENMP
LDFLAGS += -qsmp=omp
endif
endif
# PETSC for CMFD functionality
ifeq ($(PETSC),yes)
@ -222,7 +242,7 @@ endif
all: xml $(program)
xml:
cd xml; make MACHINE=$(MACHINE) F90=$(F90) F90FLAGS="$(F90FLAGS)"
cd xml; make MACHINE=$(MACHINE) F90=$(F90) F90FLAGS="$(F90FLAGS)" LDFLAGS="$(LDFLAGS)"
$(program): $(objects)
$(F90) $(objects) $(xml_lib) $(LDFLAGS) -o $@
install:

View file

@ -51,7 +51,9 @@ contains
! allocate arrays for ACE table storage and cross section cache
allocate(nuclides(n_nuclides_total))
allocate(sab_tables(n_sab_tables))
!$omp parallel
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
! ==========================================================================
! READ ALL ACE CROSS SECTION TABLES

View file

@ -57,7 +57,8 @@ contains
FILTER_SURFACE, IN_RIGHT, OUT_RIGHT, IN_FRONT, &
OUT_FRONT, IN_TOP, OUT_TOP, CMFD_NOACCEL, ZERO, ONE
use error, only: fatal_error
use global, only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes
use global, only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes,&
matching_bins
use mesh, only: mesh_indices_to_bin
use mesh_header, only: StructuredMesh
use tally_header, only: TallyObject
@ -137,21 +138,21 @@ contains
TALLY: if (ital == 1) then
! reset all bins to 1
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
! set ijk as mesh indices
ijk = (/ i, j, k /)
! get bin number for mesh indices
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
! apply energy in filter
if (i_filter_ein > 0) then
t % matching_bins(i_filter_ein) = ng - h + 1
matching_bins(i_filter_ein) = ng - h + 1
end if
! calculate score index from bins
score_index = sum((t % matching_bins - 1) * t%stride) + 1
score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
! get flux
flux = t % results(1,score_index) % sum
@ -180,24 +181,24 @@ contains
INGROUP: do g = 1, ng
! reset all bins to 1
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
! set ijk as mesh indices
ijk = (/ i, j, k /)
! get bin number for mesh indices
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
if (i_filter_ein > 0) then
! apply energy in filter
t % matching_bins(i_filter_ein) = ng - h + 1
matching_bins(i_filter_ein) = ng - h + 1
! set energy out bin
t % matching_bins(i_filter_eout) = ng - g + 1
matching_bins(i_filter_eout) = ng - g + 1
end if
! calculate score index from bins
score_index = sum((t % matching_bins - 1) * t%stride) + 1
score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
! get scattering
cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /&
@ -216,69 +217,69 @@ contains
else if (ital == 3) then
! initialize and filter for energy
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
if (i_filter_ein > 0) then
t % matching_bins(i_filter_ein) = ng - h + 1
matching_bins(i_filter_ein) = ng - h + 1
end if
! left surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i-1, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum
! right surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum
! back surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j-1, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum
! front surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum
! bottom surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k-1 /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum
! top surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum
end if TALLY

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)

View file

@ -14,6 +14,7 @@ module cross_section
save
integer :: union_grid_index
!$omp threadprivate(union_grid_index)
contains
@ -32,7 +33,8 @@ contains
integer :: j ! index in mat % i_sab_nuclides
real(8) :: atom_density ! atom density of a nuclide
logical :: check_sab ! should we check for S(a,b) table?
type(Material), pointer :: mat => null() ! current material
type(Material), pointer, save :: mat => null() ! current material
!$omp threadprivate(mat)
! Set all material macroscopic cross sections to zero
material_xs % total = ZERO
@ -139,7 +141,8 @@ contains
integer :: i_grid ! index on nuclide energy grid
real(8) :: f ! interp factor on nuclide energy grid
type(Nuclide), pointer :: nuc => null()
type(Nuclide), pointer, save :: nuc => null()
!$omp threadprivate(nuc)
! Set pointer to nuclide
nuc => nuclides(i_nuclide)
@ -258,7 +261,8 @@ contains
real(8) :: f ! interp factor on S(a,b) energy grid
real(8) :: inelastic ! S(a,b) inelastic cross section
real(8) :: elastic ! S(a,b) elastic cross section
type(SAlphaBeta), pointer :: sab => null()
type(SAlphaBeta), pointer, save :: sab => null()
!$omp threadprivate(sab)
! Set flag that S(a,b) treatment should be used for scattering
micro_xs(i_nuclide) % index_sab = i_sab
@ -346,9 +350,10 @@ contains
real(8) :: capture ! (n,gamma) cross section
real(8) :: fission ! fission cross section
real(8) :: inelastic ! inelastic cross section
type(UrrData), pointer :: urr => null()
type(Nuclide), pointer :: nuc => null()
type(Reaction), pointer :: rxn => null()
type(UrrData), pointer, save :: urr => null()
type(Nuclide), pointer, save :: nuc => null()
type(Reaction), pointer, save :: rxn => null()
!$omp threadprivate(urr, nuc, rxn)
micro_xs(i_nuclide) % use_ptable = .true.

View file

@ -39,6 +39,7 @@ contains
subroutine run_eigenvalue()
type(Particle) :: p
integer :: i_work
if (master) call header("K EIGENVALUE SIMULATION", level=1)
@ -71,7 +72,9 @@ contains
! ====================================================================
! LOOP OVER PARTICLES
PARTICLE_LOOP: do current_work = 1, work
!$omp parallel do schedule(static) firstprivate(p)
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work
! grab source particle from bank
call get_source_particle(p, current_work)
@ -80,6 +83,7 @@ contains
call transport(p)
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()
@ -126,7 +130,9 @@ contains
tallies_on = .true.
! Add user tallies to active tallies list
!$omp parallel
call setup_active_usertallies()
!$omp end parallel
end if
! check CMFD initialize batch
@ -160,6 +166,11 @@ contains
subroutine finalize_generation()
#ifdef OPENMP
! Join the fission bank from each thread into one global fission bank
call join_bank_from_threads()
#endif
! Distribute fission bank across processors evenly
call time_bank % start()
call synchronize_bank()
@ -378,7 +389,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()
@ -394,11 +405,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
@ -423,7 +434,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
@ -435,13 +446,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
@ -811,4 +821,45 @@ contains
end subroutine replay_batch_history
#ifdef OPENMP
!===============================================================================
! JOIN_BANK_FROM_THREADS
!===============================================================================
subroutine join_bank_from_threads()
integer :: total ! total number of fission bank sites
integer :: i ! loop index for threads
! Initialize the total number of fission bank sites
total = 0
!$omp parallel
! Copy thread fission bank sites to one shared copy
!$omp do ordered schedule(static)
do i = 1, n_threads
!$omp ordered
master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank)
total = total + n_bank
!$omp end ordered
end do
!$omp end do
! Make sure all threads have made it to this point
!$omp barrier
! Now copy the shared fission bank sites back to the master thread's copy.
if (thread_id == 0) then
n_bank = total
fission_bank(1:n_bank) = master_fission_bank(1:n_bank)
else
n_bank = 0
end if
!$omp end parallel
end subroutine join_bank_from_threads
#endif
end module eigenvalue

View file

@ -12,6 +12,7 @@ module fixed_source
use tracking, only: transport
type(Bank), pointer :: source_site => null()
!$omp threadprivate(source_site)
contains
@ -23,11 +24,15 @@ contains
if (master) call header("FIXED SOURCE TRANSPORT SIMULATION", level=1)
! Allocate particle and dummy source site
!$omp parallel
allocate(source_site)
!$omp end parallel
! Turn timer and tallies on
tallies_on = .true.
!$omp parallel
call setup_active_usertallies()
!$omp end parallel
call time_active % start()
! ==========================================================================
@ -47,15 +52,16 @@ contains
! =======================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(static) firstprivate(p)
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)
@ -67,6 +73,7 @@ contains
call transport(p)
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()

View file

@ -29,7 +29,8 @@ contains
integer :: i_surface ! index in surfaces array (with sign)
logical :: specified_sense ! specified sense of surface in list
logical :: actual_sense ! sense of particle wrt surface
type(Surface), pointer :: s => null()
type(Surface), pointer, save :: s => null()
!$omp threadprivate(s)
SURFACE_LOOP: do i = 1, c % n_surfaces
! Lookup surface
@ -76,9 +77,10 @@ contains
integer :: i ! cell loop index on a level
integer :: n ! number of cells to search on a level
integer :: index_cell ! index in cells array
type(Cell), pointer :: c ! pointer to cell
type(Universe), pointer :: univ ! universe to search in
type(LocalCoord), pointer :: coord ! particle coordinate to search on
type(Cell), pointer, save :: c => null() ! pointer to cell
type(Universe), pointer, save :: univ => null() ! universe to search in
type(LocalCoord), pointer, save :: coord => null() ! particle coordinate to search on
!$omp threadprivate(c, univ, coord)
coord => p % coord0
@ -139,9 +141,10 @@ contains
logical :: use_search_cells ! use cells provided as argument
logical :: outside_lattice ! if particle is not inside lattice bounds
logical :: lattice_edge ! if particle is on a lattice edge
type(Cell), pointer :: c ! pointer to cell
type(Lattice), pointer :: lat ! pointer to lattice
type(Universe), pointer :: univ ! universe to search in
type(Cell), pointer, save :: c => null() ! pointer to cell
type(Lattice), pointer, save :: lat => null() ! pointer to lattice
type(Universe), pointer, save :: univ => null() ! universe to search in
!$omp threadprivate(c, lat, univ)
! Remove coordinates for any lower levels
call deallocate_coord(p % coord % next)
@ -375,7 +378,8 @@ contains
real(8) :: norm ! "norm" of surface normal
integer :: i_surface ! index in surfaces
logical :: found ! particle found in universe?
type(Surface), pointer :: surf => null()
type(Surface), pointer, save :: surf => null()
!$omp threadprivate(surf)
i_surface = abs(p % surface)
surf => surfaces(i_surface)
@ -404,8 +408,12 @@ contains
end if
! Score to global leakage tally
if (tallies_on) global_tallies(LEAKAGE) % value = &
if (tallies_on) then
!$omp critical
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + p % wgt
!$omp end critical
end if
! Display message
if (verbosity >= 10 .or. trace) then
@ -658,7 +666,8 @@ contains
integer :: n_x, n_y, n_z ! size of lattice
real(8) :: x0, y0, z0 ! half width of lattice element
logical :: found ! particle found in cell?
type(Lattice), pointer :: lat => null()
type(Lattice), pointer, save :: lat => null()
!$omp threadprivate(lat)
lat => lattices(p % coord % lattice)
@ -787,11 +796,12 @@ contains
real(8) :: a,b,c,k ! quadratic equation coefficients
real(8) :: quad ! discriminant of quadratic equation
logical :: on_surface ! is particle on surface?
type(Cell), pointer :: cl => null()
type(Surface), pointer :: surf => null()
type(Lattice), pointer :: lat => null()
type(LocalCoord), pointer :: coord => null()
type(LocalCoord), pointer :: final_coord => null()
type(Cell), pointer, save :: cl => null()
type(Surface), pointer, save :: surf => null()
type(Lattice), pointer, save :: lat => null()
type(LocalCoord), pointer, save :: coord => null()
type(LocalCoord), pointer, save :: final_coord => null()
!$omp threadprivate(cl, surf, lat, coord, final_coord)
! inialize distance to infinity (huge)
dist = INFINITY
@ -1562,7 +1572,9 @@ contains
! Increment number of lost particles
p % alive = .false.
!$omp critical
n_lost_particles = n_lost_particles + 1
!$omp end critical
! Abort the simulation if the maximum number of lost particles has been
! reached

View file

@ -97,6 +97,7 @@ module global
type(StructuredMesh), allocatable, target :: meshes(:)
type(TallyObject), allocatable, target :: tallies(:)
integer, allocatable :: matching_bins(:)
! Pointers for different tallies
type(TallyObject), pointer :: user_tallies(:) => null()
@ -111,6 +112,8 @@ module global
type(SetInt) :: active_tracklength_tallies
type(SetInt) :: active_current_tallies
type(SetInt) :: active_tallies
!$omp threadprivate(active_analog_tallies, active_tracklength_tallies, &
!$omp& active_current_tallies, active_tallies)
! Global tallies
! 1) collision estimate of k-eff
@ -159,11 +162,12 @@ module global
! Source and fission bank
type(Bank), allocatable, target :: source_bank(:)
type(Bank), allocatable, target :: fission_bank(:)
#ifdef OPENMP
type(Bank), allocatable, target :: master_fission_bank(:)
#endif
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
@ -205,6 +209,11 @@ module global
integer :: MPI_BANK ! MPI datatype for fission bank
integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult
#ifdef OPENMP
integer :: n_threads = NONE ! number of OpenMP threads
integer :: thread_id ! ID of a given thread
#endif
! No reduction at end of batch
logical :: reduce_tallies = .true.
@ -375,6 +384,9 @@ module global
logical :: output_xs = .false.
logical :: output_tallies = .true.
!$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, &
!$omp& trace, thread_id, current_work, matching_bins)
contains
!===============================================================================
@ -432,16 +444,25 @@ contains
! Now deallocate the tally array
deallocate(tallies)
end if
if (allocated(matching_bins)) deallocate(matching_bins)
if (allocated(tally_maps)) deallocate(tally_maps)
! Deallocate energy grid
if (allocated(e_grid)) deallocate(e_grid)
! Deallocate fission and source bank and entropy
!$omp parallel
if (allocated(fission_bank)) deallocate(fission_bank)
!$omp end parallel
#ifdef OPENMP
if (allocated(master_fission_bank)) deallocate(master_fission_bank)
#endif
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)

File diff suppressed because it is too large Load diff

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

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
@ -26,6 +26,10 @@ module initialize
use mpi
#endif
#ifdef OPENMP
use omp_lib
#endif
#ifdef HDF5
use hdf5_interface
use hdf5_summary, only: hdf5_write_summary
@ -306,6 +310,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)
@ -330,7 +335,7 @@ contains
run_mode = MODE_PLOTTING
check_overlaps = .true.
case ('-n', '-n_particles', '--n_particles')
case ('-n', '-particles', '--particles')
! Read number of particles per cycle
i = i + 1
n_particles = str_to_int(argv(i))
@ -346,9 +351,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)
@ -363,6 +368,23 @@ contains
case ('-g', '-geometry-debug', '--geometry-debug')
check_overlaps = .true.
case ('-s', '--threads')
! Read number of threads
i = i + 1
#ifdef OPENMP
! Read and set number of OpenMP threads
n_threads = str_to_int(argv(i))
if (n_threads < 1) then
message = "Invalid number of threads specified on command line."
call fatal_error()
end if
call omp_set_num_threads(n_threads)
#else
message = "Ignoring number of threads specified on command line."
call warning()
#endif
case ('-?', '-help', '--help')
call print_usage()
stop
@ -754,15 +776,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
@ -772,10 +816,10 @@ contains
subroutine allocate_banks()
integer :: alloc_err ! allocation error code
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
@ -783,8 +827,26 @@ contains
call fatal_error()
end if
! Allocate fission bank
allocate(fission_bank(3*maxwork), STAT=alloc_err)
#ifdef OPENMP
! If OpenMP is being used, each thread needs its own private fission
! bank. Since the private fission banks need to be combined at the end of a
! generation, there is also a 'master_fission_bank' that is used to collect
! the sites from each thread.
!$omp parallel
n_threads = omp_get_num_threads()
thread_id = omp_get_thread_num()
if (thread_id == 0) then
allocate(fission_bank(3*work))
else
allocate(fission_bank(3*work/n_threads))
end if
!$omp end parallel
allocate(master_fission_bank(3*work), STAT=alloc_err)
#else
allocate(fission_bank(3*work), STAT=alloc_err)
#endif
! Check for allocation errors
if (alloc_err /= 0) then

View file

@ -225,6 +225,23 @@ contains
call get_node_value(node_verb, "value", verbosity)
end if
! Number of OpenMP threads
if (threads_ /= NONE) then
#ifdef OPENMP
if (n_threads == NONE) then
n_threads = threads_
if (n_threads < 1) then
message = "Invalid number of threads: " // to_str(n_threads)
call fatal_error()
end if
call omp_set_num_threads(n_threads)
end if
#else
message = "Ignoring number of threads."
call warning()
#endif
end if
! ==========================================================================
! EXTERNAL SOURCE

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

View file

@ -167,9 +167,11 @@ contains
write(OUTPUT_UNIT,*)
write(OUTPUT_UNIT,*) 'Options:'
write(OUTPUT_UNIT,*) ' -g, --geometry-debug Run in geometry debugging mode'
write(OUTPUT_UNIT,*) ' -n, --particles Number of particles per generation'
write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode'
write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run from a state point'
write(OUTPUT_UNIT,*) ' or a particle restart file'
write(OUTPUT_UNIT,*) ' -s, --threads Number of OpenMP threads'
write(OUTPUT_UNIT,*) ' -v, --version Show version information'
write(OUTPUT_UNIT,*) ' -?, --help Show this message'
end if
@ -1644,7 +1646,7 @@ contains
! to be used for a given tally.
! Initialize bins, filter level, and indentation
t % matching_bins = 0
matching_bins(1:t%n_filters) = 0
j = 1
indent = 0
@ -1654,16 +1656,16 @@ contains
if (t % n_filters == 0) exit find_bin
! Increment bin combination
t % matching_bins(j) = t % matching_bins(j) + 1
matching_bins(j) = matching_bins(j) + 1
! =================================================================
! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER
if (t % matching_bins(j) > t % filters(j) % n_bins) then
if (matching_bins(j) > t % filters(j) % n_bins) then
! If this is the first filter, then exit
if (j == 1) exit print_bin
t % matching_bins(j) = 0
matching_bins(j) = 0
j = j - 1
indent = indent - 2
@ -1696,7 +1698,7 @@ contains
! bins below the lowest filter level will be zeros
if (t % n_filters > 0) then
filter_index = sum((max(t % matching_bins,1) - 1) * t % stride) + 1
filter_index = sum((max(matching_bins(1:t%n_filters),1) - 1) * t % stride) + 1
else
filter_index = 1
end if
@ -1805,7 +1807,7 @@ contains
m => meshes(t % filters(i_filter_mesh) % int_bins(1))
! initialize bins array
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
! determine how many energy in bins there are
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
@ -1831,7 +1833,7 @@ contains
do l = 1, n
if (print_ebin) then
! Set incoming energy bin
t % matching_bins(i_filter_ein) = l
matching_bins(i_filter_ein) = l
! Write incoming energy bin
write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') &
@ -1839,102 +1841,102 @@ contains
end if
! Left Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Right Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Back Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Front Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Bottom Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Top Surface
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
to_str(t % results(1,filter_index) % sum), &
@ -1965,7 +1967,7 @@ contains
real(8) :: E1 ! upper bound for energy bin
type(StructuredMesh), pointer :: m => null()
bin = t % matching_bins(i_filter)
bin = matching_bins(i_filter)
select case(t % filters(i_filter) % type)
case (FILTER_UNIVERSE)

File diff suppressed because it is too large Load diff

View file

@ -2,123 +2,27 @@ 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 output, only: write_message, print_particle
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use tracking, only: transport
#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(p)
type(Particle), intent(inout) :: p
! 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(p)
type(Particle), intent(inout) :: p
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()
@ -126,30 +30,69 @@ contains
integer(8) :: particle_seed
type(Particle) :: p
! initialize the particle to be tracked
! Set verbosity high
verbosity = 10
! Initialize the particle to be tracked
call p % initialize()
! read in the restart information
#ifdef HDF5
call read_hdf5_particle_restart(p)
#else
call read_binary_particle_restart(p)
#endif
! Read in the restart information
call read_particle_restart(p)
! set all tallies to 0 for now (just tracking errors)
! 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
! Transport neutron
call transport(p)
! write output if particle made it
! 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

View file

@ -1,132 +1,65 @@
module particle_restart_write
use, intrinsic :: ISO_FORTRAN_ENV
use bank_header, only: Bank
use bank_header, only: Bank
use global
use particle_header, only: Particle
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(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(p)
filename = trim(filename) // '.h5'
#else
call write_particle_restart_binary(p)
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(p)
type(Particle), intent(in) :: p
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(p % id)) // '.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(p)
type(Particle), intent(in) :: p
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(p % id)) // '.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

View file

@ -74,7 +74,8 @@ contains
integer :: i_nuclide ! index in nuclides array
integer :: i_reaction ! index in nuc % reactions array
type(Nuclide), pointer :: nuc => null()
type(Nuclide), pointer, save :: nuc => null()
!$omp threadprivate(nuc)
i_nuclide = sample_nuclide(p, 'total ')
@ -129,7 +130,8 @@ contains
real(8) :: cutoff
real(8) :: atom_density ! atom density of nuclide in atom/b-cm
real(8) :: sigma ! microscopic total xs for nuclide
type(Material), pointer :: mat => null()
type(Material), pointer, save :: mat => null()
!$omp threadprivate(mat)
! Get pointer to current material
mat => materials(p % material)
@ -191,8 +193,9 @@ contains
real(8) :: f
real(8) :: prob
real(8) :: cutoff
type(Nuclide), pointer :: nuc => null()
type(Reaction), pointer :: rxn => null()
type(Nuclide), pointer, save :: nuc => null()
type(Reaction), pointer, save :: rxn => null()
!$omp threadprivate(nuc, rxn)
! Get pointer to nuclide
nuc => nuclides(i_nuclide)
@ -251,18 +254,22 @@ contains
p % last_wgt = p % wgt
! Score implicit absorption estimate of keff
!$omp critical
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % absorb_wgt * &
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
!$omp end critical
else
! See if disappearance reaction happens
if (micro_xs(i_nuclide) % absorption > &
prn() * micro_xs(i_nuclide) % total) then
! Score absorption estimate of keff
!$omp critical
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % wgt * &
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
!$omp end critical
p % alive = .false.
p % event = EVENT_ABSORB
@ -306,8 +313,9 @@ contains
real(8) :: f
real(8) :: prob
real(8) :: cutoff
type(Nuclide), pointer :: nuc => null()
type(Reaction), pointer :: rxn => null()
type(Nuclide), pointer, save :: nuc => null()
type(Reaction), pointer, save :: rxn => null()
!$omp threadprivate(nuc, rxn)
! Get pointer to nuclide and grid index/interpolation factor
nuc => nuclides(i_nuclide)
@ -410,7 +418,8 @@ contains
real(8) :: v_cm(3) ! velocity of center-of-mass
real(8) :: v_t(3) ! velocity of target nucleus
real(8) :: uvw_cm(3) ! directional cosines in center-of-mass
type(Nuclide), pointer :: nuc => null()
type(Nuclide), pointer, save :: nuc => null()
!$omp threadprivate(nuc)
! get pointer to nuclide
nuc => nuclides(i_nuclide)
@ -487,7 +496,8 @@ contains
real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j)
real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j)
real(8) :: prob ! probability for sampling Bragg edge
type(SAlphaBeta), pointer :: sab => null()
type(SAlphaBeta), pointer, save :: sab => null()
!$omp threadprivate(sab)
! Get pointer to S(a,b) table
sab => sab_tables(i_sab)
@ -715,8 +725,9 @@ contains
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
logical :: in_mesh ! source site in ufs mesh?
type(Nuclide), pointer :: nuc => null()
type(Reaction), pointer :: rxn => null()
type(Nuclide), pointer, save :: nuc => null()
type(Reaction), pointer, save :: rxn => null()
!$omp threadprivate(nuc, rxn)
! Get pointers
nuc => nuclides(i_nuclide)
@ -816,7 +827,8 @@ contains
real(8) :: xi ! random number
real(8) :: yield ! delayed neutron precursor yield
real(8) :: prob ! cumulative probability
type(DistEnergy), pointer :: edist => null()
type(DistEnergy), pointer, save :: edist => null()
!$omp threadprivate(edist)
! Determine total nu
nu_t = nu_total(nuc, E)

View file

@ -15,6 +15,8 @@ module random_lcg
integer(8) :: prn_stride ! stride between particles
real(8) :: prn_norm ! 2^(-M)
!$omp threadprivate(prn_seed)
public :: prn
public :: initialize_prng
public :: set_particle_seed

View file

@ -101,6 +101,8 @@ element settings {
element survival_biasing { xsd:boolean }? &
element threads { xsd:positiveInteger }? &
element trace { list { xsd:positiveInteger+ } }? &
element verbosity { xsd:positiveInteger }? &

View file

@ -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
@ -155,7 +155,8 @@ contains
integer(8), intent(in) :: index_source
integer(8) :: particle_seed ! unique index for particle
type(Bank), pointer :: src => null()
type(Bank), pointer, save :: src => null()
!$omp threadprivate(src)
! set defaults
call p % initialize()
@ -165,7 +166,7 @@ contains
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

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

View file

@ -23,6 +23,7 @@ module tally
! Tally map positioning array
integer :: position(N_FILTER_TYPES - 3) = 0
!$omp threadprivate(position)
contains
@ -54,7 +55,8 @@ contains
real(8) :: macro_total ! material macro total xs
real(8) :: macro_scatt ! material macro scatt xs
logical :: found_bin ! scoring bin found?
type(TallyObject), pointer :: t => null()
type(TallyObject), pointer, save :: t => null()
!$omp threadprivate(t)
! Copy particle's pre- and post-collision weight and angle
last_wgt = p % last_wgt
@ -83,7 +85,7 @@ contains
! be accumulating the tally values
! Determine scoring index for this filter combination
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
! Check for nuclide bins
k = 0
@ -200,8 +202,10 @@ contains
! get the score and tally it
score = last_wgt * calc_pn(n, mu)
!$omp critical
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do
j = j + t % scatt_order(j)
cycle SCORE_LOOP
@ -338,8 +342,10 @@ contains
end select
! Add score to tally
!$omp critical
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do SCORE_LOOP
@ -382,7 +388,7 @@ contains
! save original outgoing energy bin and score index
i = t % find_filter(FILTER_ENERGYOUT)
bin_energyout = t % matching_bins(i)
bin_energyout = matching_bins(i)
! Get number of energies on filter
n = size(t % filters(i) % real_bins)
@ -405,18 +411,20 @@ contains
E_out > t % filters(i) % real_bins(n)) cycle
! change outgoing energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, n, E_out)
matching_bins(i) = binary_search(t % filters(i) % real_bins, n, E_out)
! determine scoring index
i_filter = sum((t % matching_bins - 1) * t % stride) + 1
i_filter = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
! Add score to tally
!$omp critical
t % results(i_score, i_filter) % value = &
t % results(i_score, i_filter) % value + score
!$omp end critical
end do
! reset outgoing energy bin and score index
t % matching_bins(i) = bin_energyout
matching_bins(i) = bin_energyout
end subroutine score_fission_eout
@ -449,9 +457,10 @@ contains
real(8) :: score ! actual score (e.g., flux*xs)
real(8) :: atom_density ! atom density of single nuclide in atom/b-cm
logical :: found_bin ! scoring bin found?
type(TallyObject), pointer :: t => null()
type(Material), pointer :: mat => null()
type(Reaction), pointer :: rxn => null()
type(TallyObject), pointer, save :: t => null()
type(Material), pointer, save :: mat => null()
type(Reaction), pointer, save :: rxn => null()
!$omp threadprivate(t, mat, rxn)
! Determine track-length estimate of flux
flux = p % wgt * distance
@ -486,7 +495,7 @@ contains
! be accumulating the tally values
! Determine scoring index for this filter combination
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
if (t % all_nuclides) then
call score_all_nuclides(p, i_tally, flux, filter_index)
@ -698,8 +707,10 @@ contains
score_index = (k - 1)*t % n_score_bins + j
! Add score to tally
!$omp critical
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do SCORE_LOOP
@ -742,9 +753,10 @@ contains
real(8) :: f ! interpolation factor
real(8) :: score ! actual scoring tally value
real(8) :: atom_density ! atom density of single nuclide in atom/b-cm
type(TallyObject), pointer :: t => null()
type(Material), pointer :: mat => null()
type(Reaction), pointer :: rxn => null()
type(TallyObject), pointer, save :: t => null()
type(Material), pointer, save :: mat => null()
type(Reaction), pointer, save :: rxn => null()
!$omp threadprivate(t, mat, rxn)
! Get pointer to tally
t => tallies(i_tally)
@ -840,8 +852,10 @@ contains
score_index = (i_nuclide - 1)*t % n_score_bins + j
! Add score to tally
!$omp critical
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do SCORE_LOOP
@ -938,8 +952,10 @@ contains
score_index = n_nuclides_total*t % n_score_bins + j
! Add score to tally
!$omp critical
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do MATERIAL_SCORE_LOOP
@ -982,13 +998,14 @@ contains
logical :: found_bin ! was a scoring bin found?
logical :: start_in_mesh ! starting coordinates inside mesh?
logical :: end_in_mesh ! ending coordinates inside mesh?
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
type(Material), pointer :: mat => null()
type(LocalCoord), pointer :: coord => null()
type(TallyObject), pointer, save :: t => null()
type(StructuredMesh), pointer, save :: m => null()
type(Material), pointer, save :: mat => null()
type(LocalCoord), pointer, save :: coord => null()
!$omp threadprivate(t, m, mat, coord)
t => tallies(i_tally)
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
! ==========================================================================
! CHECK IF THIS TRACK INTERSECTS THE MESH
@ -1028,11 +1045,11 @@ contains
case (FILTER_UNIVERSE)
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
case (FILTER_MATERIAL)
t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
p % material, i_tally)
case (FILTER_CELL)
@ -1040,21 +1057,21 @@ contains
coord => p % coord0
do while(associated(coord))
position(FILTER_CELL) = 0
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
matching_bins(i) = get_next_bin(FILTER_CELL, &
coord % cell, i_tally)
if (t % matching_bins(i) /= NO_BIN_FOUND) exit
if (matching_bins(i) /= NO_BIN_FOUND) exit
coord => coord % next
end do
nullify(coord)
case (FILTER_CELLBORN)
! determine next cellborn bin
t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
p % cell_born, i_tally)
case (FILTER_SURFACE)
! determine next surface bin
t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
matching_bins(i) = get_next_bin(FILTER_SURFACE, &
p % surface, i_tally)
case (FILTER_ENERGYIN)
@ -1064,17 +1081,17 @@ contains
! check if energy of the particle is within energy bins
if (p % E < t % filters(i) % real_bins(1) .or. &
p % E > t % filters(i) % real_bins(k + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
matching_bins(i) = binary_search(t % filters(i) % real_bins, &
k + 1, p % E)
end if
end select
! Check if no matching bin was found
if (t % matching_bins(i) == NO_BIN_FOUND) return
if (matching_bins(i) == NO_BIN_FOUND) return
end do FILTER_LOOP
@ -1144,10 +1161,10 @@ contains
flux = p % wgt * distance
! Determine mesh bin
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, ijk_cross)
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, ijk_cross)
! Determining scoring index
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
if (t % all_nuclides) then
! Score reaction rates for each nuclide in material
@ -1245,8 +1262,10 @@ contains
score_index = (b - 1)*t % n_score_bins + j
! Add score to tally
!$omp critical
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do SCORE_LOOP
@ -1275,13 +1294,14 @@ contains
integer :: i ! loop index for filters
integer :: n ! number of bins for single filter
real(8) :: E ! particle energy
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
type(LocalCoord), pointer :: coord => null()
type(TallyObject), pointer, save :: t => null()
type(StructuredMesh), pointer, save :: m => null()
type(LocalCoord), pointer, save :: coord => null()
!$omp threadprivate(t, m, coord)
found_bin = .true.
t => tallies(i_tally)
t % matching_bins = 1
matching_bins(1:t%n_filters) = 1
FILTER_LOOP: do i = 1, t % n_filters
@ -1291,16 +1311,16 @@ contains
m => meshes(t % filters(i) % int_bins(1))
! Determine if we're in the mesh first
call get_mesh_bin(m, p % coord0 % xyz, t % matching_bins(i))
call get_mesh_bin(m, p % coord0 % xyz, matching_bins(i))
case (FILTER_UNIVERSE)
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
case (FILTER_MATERIAL)
t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
p % material, i_tally)
case (FILTER_CELL)
@ -1308,21 +1328,21 @@ contains
coord => p % coord0
do while(associated(coord))
position(FILTER_CELL) = 0
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
matching_bins(i) = get_next_bin(FILTER_CELL, &
coord % cell, i_tally)
if (t % matching_bins(i) /= NO_BIN_FOUND) exit
if (matching_bins(i) /= NO_BIN_FOUND) exit
coord => coord % next
end do
nullify(coord)
case (FILTER_CELLBORN)
! determine next cellborn bin
t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
p % cell_born, i_tally)
case (FILTER_SURFACE)
! determine next surface bin
t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
matching_bins(i) = get_next_bin(FILTER_SURFACE, &
p % surface, i_tally)
case (FILTER_ENERGYIN)
@ -1339,10 +1359,10 @@ contains
! check if energy of the particle is within energy bins
if (E < t % filters(i) % real_bins(1) .or. &
E > t % filters(i) % real_bins(n + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
matching_bins(i) = binary_search(t % filters(i) % real_bins, &
n + 1, E)
end if
@ -1353,17 +1373,17 @@ contains
! check if energy of the particle is within energy bins
if (p % E < t % filters(i) % real_bins(1) .or. &
p % E > t % filters(i) % real_bins(n + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
matching_bins(i) = binary_search(t % filters(i) % real_bins, &
n + 1, p % E)
end if
end select
! If the current filter didn't match, exit this subroutine
if (t % matching_bins(i) == NO_BIN_FOUND) then
if (matching_bins(i) == NO_BIN_FOUND) then
found_bin = .false.
return
end if
@ -1403,8 +1423,9 @@ contains
logical :: x_same ! same starting/ending x index (i)
logical :: y_same ! same starting/ending y index (j)
logical :: z_same ! same starting/ending z index (k)
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
type(TallyObject), pointer, save :: t => null()
type(StructuredMesh), pointer, save :: m => null()
!$omp threadprivate(t, m)
TALLY_LOOP: do i = 1, active_current_tallies % size()
! Copy starting and ending location of particle
@ -1454,7 +1475,7 @@ contains
end if
! search to find incoming energy bin
t % matching_bins(j) = binary_search(t % filters(j) % real_bins, &
matching_bins(j) = binary_search(t % filters(j) % real_bins, &
n + 1, p % E)
end if
@ -1471,24 +1492,28 @@ contains
do j = ijk0(3), ijk1(3) - 1
ijk0(3) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_TOP
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_TOP
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
else
do j = ijk0(3) - 1, ijk1(3), -1
ijk0(3) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_TOP
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_TOP
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
end if
@ -1499,24 +1524,28 @@ contains
do j = ijk0(2), ijk1(2) - 1
ijk0(2) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_FRONT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_FRONT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
else
do j = ijk0(2) - 1, ijk1(2), -1
ijk0(2) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_FRONT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_FRONT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
end if
@ -1527,24 +1556,28 @@ contains
do j = ijk0(1), ijk1(1) - 1
ijk0(1) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_RIGHT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_RIGHT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
else
do j = ijk0(1) - 1, ijk1(1), -1
ijk0(1) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_RIGHT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_RIGHT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
end if
@ -1565,7 +1598,7 @@ contains
do k = 1, n_cross
! Reset scoring bin index
t % matching_bins(i_filter_surf) = 0
matching_bins(i_filter_surf) = 0
! Calculate distance to each bounding surface. We need to treat
! special case where the cosine of the angle is zero since this would
@ -1592,8 +1625,8 @@ contains
! Crossing into right mesh cell -- this is treated as outgoing
! current from (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_RIGHT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_RIGHT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(1) = ijk0(1) + 1
@ -1604,8 +1637,8 @@ contains
ijk0(1) = ijk0(1) - 1
xyz_cross(1) = xyz_cross(1) - m % width(1)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_RIGHT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_RIGHT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
@ -1614,8 +1647,8 @@ contains
! Crossing into front mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_FRONT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_FRONT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(2) = ijk0(2) + 1
@ -1626,8 +1659,8 @@ contains
ijk0(2) = ijk0(2) - 1
xyz_cross(2) = xyz_cross(2) - m % width(2)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_FRONT
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_FRONT
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
@ -1636,8 +1669,8 @@ contains
! Crossing into top mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_TOP
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = OUT_TOP
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(3) = ijk0(3) + 1
@ -1648,16 +1681,16 @@ contains
ijk0(3) = ijk0(3) - 1
xyz_cross(3) = xyz_cross(3) - m % width(3)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_TOP
t % matching_bins(i_filter_mesh) = &
matching_bins(i_filter_surf) = IN_TOP
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
end if
! Determine scoring index
if (t % matching_bins(i_filter_surf) > 0) then
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
if (matching_bins(i_filter_surf) > 0) then
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
! Check for errors
if (filter_index <= 0 .or. filter_index > &
@ -1667,8 +1700,10 @@ contains
end if
! Add to surface current tally
!$omp critical
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
! Calculate new coordinates

View file

@ -86,7 +86,6 @@ module tally_header
! mapped onto one dimension in the results array, the stride attribute gives
! the stride for a given filter type within the results array
integer, allocatable :: matching_bins(:)
integer, allocatable :: stride(:)
! This array provides a way to lookup what index in the filters array a
@ -169,8 +168,6 @@ module tally_header
deallocate(this % filters)
end if
if (allocated(this % matching_bins)) &
deallocate(this % matching_bins)
if (allocated(this % stride)) &
deallocate(this % stride)

View file

@ -31,9 +31,10 @@ contains
subroutine setup_tally_arrays()
integer :: i ! loop index for tallies
integer :: j ! loop index for filters
integer :: n ! temporary stride
integer :: i ! loop index for tallies
integer :: j ! loop index for filters
integer :: n ! temporary stride
integer :: max_n_filters = 0 ! maximum number of filters
type(TallyObject), pointer :: t => null()
TALLY_LOOP: do i = 1, n_tallies
@ -42,7 +43,7 @@ contains
! Allocate stride and matching_bins arrays
allocate(t % stride(t % n_filters))
allocate(t % matching_bins(t % n_filters))
max_n_filters = max(max_n_filters, t % n_filters)
! The filters are traversed in opposite order so that the last filter has
! the shortest stride in memory and the first filter has the largest
@ -63,6 +64,11 @@ contains
end do TALLY_LOOP
! Allocate array for matching filter bins
!$omp parallel
allocate(matching_bins(max_n_filters))
!$omp end parallel
end subroutine setup_tally_arrays
!===============================================================================

View file

@ -32,7 +32,8 @@ contains
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()
type(LocalCoord), pointer, save :: coord => null()
!$omp threadprivate(coord)
! Display message if high verbosity or trace is on
if (verbosity >= 9 .or. trace) then
@ -59,7 +60,9 @@ contains
n_event = 0
! Add paricle's starting weight to count for normalizing tallies later
!$omp critical
total_weight = total_weight + p % wgt
!$omp end critical
! Force calculation of cross-sections by setting last energy to zero
micro_xs % last_E = ZERO
@ -99,9 +102,11 @@ contains
call score_tracklength_tally(p, distance)
! Score track-length estimate of k-eff
!$omp critical
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + p % wgt * distance * &
material_xs % nu_fission
!$omp end critical
if (d_collision > d_boundary) then
! ====================================================================
@ -125,9 +130,11 @@ contains
! PARTICLE HAS COLLISION
! Score collision estimate of keff
!$omp critical
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + p % wgt * &
material_xs % nu_fission / material_xs % total
!$omp end critical
! 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

View file

@ -55,7 +55,11 @@ def run_suite(name=None, mpi=False):
results.append((name, result))
# set mpiexec path
mpiexec = '/opt/mpich/3.0.4-gnu/bin/mpiexec'
if 'COMPILER' in os.environ:
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()
@ -63,23 +67,23 @@ sys.path.append(pwd)
# Set list of tests, either default or from command line
flags = []
tests = ['compile', 'gfortran', 'gfortran-dbg', 'gfortran-opt',
'gfortran-hdf5', 'gfortran-mpi', 'gfortran-phdf5',
'gfortran-petsc', 'gfortran-phdf5-petsc',
'gfortran-phdf5-petsc-opt']
tests = ['compile', 'normal', 'debug', 'optimize', 'mpi', 'omp', 'hdf5',
'petsc', 'mpi-omp', 'omp-hdf5', 'phdf5', 'phdf5-petsc',
'phdf5-petsc-optimize', 'omp-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_ = [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 ['gfortran', 'gfortran-dbg', 'gfortran-opt', 'gfortran-hdf5']:
elif name in ['normal', 'debug', 'optimize', 'omp', 'hdf5', 'omp-hdf5']:
run_suite(name=name)
elif name in ['gfortran-mpi', 'gfortran-phdf5', 'gfortran-petsc',
'gfortran-phdf5-petsc', 'gfortran-phdf5-petsc-opt']:
elif name in ['mpi', 'mpi-omp' 'phdf5', 'petsc', 'phdf5-petsc',
'phdf5-petsc-optimize', 'omp-phdf5-petsc-optimize']:
run_suite(name=name, mpi=True)
# print out summary of results

View file

@ -3,8 +3,8 @@
"""Compilation tests
This set of tests makes sure that OpenMC can compile for many different
combinations of compilation flags and options. By default, only the gfortran and
ifort compilers are tested. This requires that the MPI, HDF5, and PETSC
combinations of compilation flags and options. By default, only the gfortran
and ifort compilers are tested. This requires that the MPI, HDF5, and PETSC
libraries are already set up appropriately in the Makefile.
"""
@ -15,68 +15,115 @@ import shutil
pwd = os.path.dirname(__file__)
if 'COMPILER' in os.environ:
compiler = 'COMPILER=' + os.environ['COMPILER']
else:
compiler = 'COMPILER=gnu'
def setup():
# Change to source directory
os.chdir(pwd + '/../../src')
def test_gfortran():
returncode = run(['make','distclean'])
returncode = run(['make'])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran')
def test_gfortran_debug():
returncode = run(['make','distclean'])
returncode = run(['make', 'DEBUG=yes'])
def test_normal():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran-dbg')
shutil.move('openmc', 'openmc-normal')
def test_gfortran_profile():
returncode = run(['make','distclean'])
returncode = run(['make', 'PROFILE=yes'])
def test_debug():
returncode = run(['make', 'distclean'])
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_optimize():
returncode = run(['make','distclean'])
returncode = run(['make', 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran-opt')
def test_gfortran_mpi():
returncode = run(['make','distclean'])
returncode = run(['make', 'MPI=yes'])
def test_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran-mpi')
shutil.move('openmc', 'openmc-optimize')
def test_gfortran_hdf5():
returncode = run(['make','distclean'])
returncode = run(['make', 'HDF5=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran-hdf5')
def test_gfortran_petsc():
returncode = run(['make','distclean'])
returncode = run(['make', 'MPI=yes', 'PETSC=yes'])
def test_mpi():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran-petsc')
shutil.move('openmc', 'openmc-mpi')
def test_gfortran_mpi_hdf5():
returncode = run(['make','distclean'])
returncode = run(['make', 'MPI=yes', 'HDF5=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran-phdf5')
def test_gfortran_mpi_hdf5_petsc():
returncode = run(['make','distclean'])
returncode = run(['make', 'MPI=yes', 'HDF5=yes', 'PETSC=yes'])
def test_omp():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'OPENMP=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran-phdf5-petsc')
shutil.move('openmc', 'openmc-omp')
def test_gfortran_mpi_hdf5_petsc_optimize():
returncode = run(['make','distclean'])
returncode = run(['make', 'MPI=yes', 'HDF5=yes', 'PETSC=yes', 'OPTIMIZE=yes'])
def test_hdf5():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'HDF5=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-gfortran-phdf5-petsc-opt')
shutil.move('openmc', 'openmc-hdf5')
def test_petsc():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'PETSC=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-petsc')
def test_mpi_omp():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-mpi-omp')
def test_omp_hdf5():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OPENMP=yes', 'HDF5=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-hdf5')
def test_mpi_hdf5():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5')
def test_mpi_hdf5_petsc():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-petsc')
def test_mpi_hdf5_petsc_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'HDF5=yes', 'PETSC=yes',
'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-phdf5-petsc-optimize')
def test_mpi_omp_hdf5_petsc_optimize():
returncode = run(['make', 'distclean'])
returncode = run(['make', compiler, 'MPI=yes', 'OMP=yes', 'HDF5=yes',
'PETSC=yes', 'OPTIMIZE=yes'])
assert returncode == 0
shutil.move('openmc', 'openmc-omp-phdf5-petsc-optimize')
def run(commands):
proc = Popen(commands, stderr=STDOUT, stdout=PIPE)
@ -84,6 +131,7 @@ def run(commands):
print(proc.communicate()[0])
return returncode
def teardown(commands):
returncode = run(['make','distclean'])
shutil.copy('openmc-gfortran','openmc')
returncode = run(['make', 'distclean'])
shutil.copy('openmc-normal', 'openmc')