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

merged master with cmfd and resolved conflicts

Browse source
This commit is contained in:
Bryan Herman 2012-01-19 12:56:54 -05:00
commit cee2ac2eaa
18 changed files with 410 additions and 180 deletions
Download patch file Download diff file Expand all files Collapse all files

2
docs/source/developers.rst
View file

@ -13,8 +13,10 @@ Advisors to the project include:
* `Benoit Forget`_
* `Kord Smith`_
* `Andrew Siegel`_
.. _Paul Romano: mailto:romano7@mit.edu
.. _Bryan Herman: mailto:bherman@mit.edu
.. _Benoit Forget: mailto:bforget@mit.edu
.. _Kord Smith: mailto:kord@mit.edu
.. _Andrew Siegel: mailto:siegela@mcs.anl.gov

10
src/DEPENDENCIES
View file

@ -54,6 +54,12 @@ energy_grid.o: datatypes_header.o
energy_grid.o: global.o
energy_grid.o: output.o
finalize.o: global.o
finalize.o: hdf5_interface.o
finalize.o: output.o
finalize.o: tally.o
finalize.o: timing.o
fission.o: ace_header.o
fission.o: constants.o
fission.o: error.o
@ -84,6 +90,9 @@ global.o: source_header.o
global.o: tally_header.o
global.o: timing.o
hdf5_interface.o: constants.o
hdf5_interface.o: global.o
initialize.o: ace.o
initialize.o: bank_header.o
initialize.o: constants.o
@ -94,6 +103,7 @@ initialize.o: error.o
initialize.o: geometry.o
initialize.o: geometry_header.o
initialize.o: global.o
initialize.o: hdf5_interface.o
initialize.o: input_xml.o
initialize.o: output.o
initialize.o: random_lcg.o

100
src/Makefile
View file

@ -15,52 +15,25 @@ include OBJECTS
# User Options
#===============================================================================
COMPILER = petsc
DEBUG = yes
PROFILE = no
OPTIMIZE = no
USE_MPI = no
USE_HDF5 = yes
#===============================================================================
# Intel Fortran compiler options
#===============================================================================
ifeq ($(COMPILER),intel)
F90 = ifort
F90FLAGS := -fpp -warn -assume byterecl
LDFLAGS =
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -traceback -ftrapuv -fp-stack-check -check all -fpe0
LDFLAGS += -g
endif
# Profiling
ifeq ($(PROFILE),yes)
F90FLAGS += -pg
LDFLAGS += -pg
endif
# Optimization
ifeq ($(OPTIMIZE),yes)
F90FLAGS += -O3 -ipo
endif
endif
COMPILER = yes
DEBUG = yes
PROFILE = no
OPTIMIZE = no
USE_MPI = no
USE_HDF5 = no
#===============================================================================
# GNU Fortran compiler options
#===============================================================================
ifeq ($(COMPILER),gfortran)
ifeq ($(COMPILER),gnu)
F90 = gfortran
F90FLAGS := -cpp -Wall
F90FLAGS := -cpp -fbacktrace
LDFLAGS =
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -pedantic -std=f2008 -fbacktrace -fbounds-check \
F90FLAGS += -g -Wall -pedantic -std=f2008 -fbounds-check \
-ffpe-trap=invalid,zero,overflow,underflow
LDFLAGS += -g
endif
@ -77,6 +50,33 @@ ifeq ($(COMPILER),gfortran)
endif
endif
#===============================================================================
# Intel Fortran compiler options
#===============================================================================
ifeq ($(COMPILER),intel)
F90 = ifort
F90FLAGS := -cpp -warn -assume byterecl -traceback
LDFLAGS =
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -ftrapuv -fp-stack-check -check all -fpe0
LDFLAGS += -g
endif
# Profiling
ifeq ($(PROFILE),yes)
F90FLAGS += -pg
LDFLAGS += -pg
endif
# Optimization
ifeq ($(OPTIMIZE),yes)
F90FLAGS += -O3 -ipo
endif
endif
# use petsc compiler
ifeq ($(COMPILER),petsc)
F90 = ${FLINKER}
@ -88,11 +88,6 @@ ifeq ($(COMPILER),petsc)
# LDFLAGS = ${PETSC_SYS_LIB}
LDFLAGS = ${SLEPC_LIB}
ifeq ($(USE_HDF5),yes)
F90FLAGS += -I${HDF5_PREFIX}/include
LDFLAGS += -L${HDF5_PREFIX}/lib ${HDF5_PREFIX}/lib/libhdf5hl_fortran.a ${HDF5_PREFIX}/lib/libhdf5_hl.a ${HDF5_PREFIX}/lib/libhdf5_fortran.a ${HDF5_PREFIX}/lib/libhdf5.a -lz -lrt -lm -Wl,-rpath -Wl,${HDF5_PREFIX}/lib
endif
ifeq ($(DEBUG),yes)
F90FLAGS += -g -fbacktrace -fbounds-check
LDFLAGS += -g
@ -106,12 +101,12 @@ endif
ifeq ($(COMPILER),pgi)
F90 = pgf90
F90FLAGS := -Mpreprocess -DNO_F2008 -Minform=inform
F90FLAGS := -Mpreprocess -DNO_F2008 -Minform=inform -traceback
LDFLAGS =
# Debugging
ifeq ($(DEBUG),yes)
F90FLAGS += -g -Mbounds -Mchkptr -Mchkstk -traceback
F90FLAGS += -g -Mbounds -Mchkptr -Mchkstk
LDFLAGS += -g
endif
@ -175,11 +170,22 @@ endif
# Use MPI for parallelism
ifeq ($(USE_MPI),yes)
MPI = /opt/mpich2/1.4.1-intel
F90 = $(MPI)/bin/mpif90
MPI_ROOT = /opt/mpich2/1.4.1-$(COMPILER)
F90 = $(MPI_ROOT)/bin/mpif90
F90FLAGS += -DMPI
endif
# Use HDF5 for output
ifeq ($(USE_HDF5),yes)
HDF5_ROOT = /opt/hdf5/1.8.8-$(COMPILER)
F90FLAGS += -DHDF5 -I${HDF5_ROOT}/include
LDFLAGS += -L${HDF5_ROOT}/lib ${HDF5_ROOT}/lib/libhdf5hl_fortran.a \
${HDF5_ROOT}/lib/libhdf5_hl.a ${HDF5_ROOT}/lib/libhdf5_fortran.a \
${HDF5_ROOT}/lib/libhdf5.a -lz -lrt -lm -Wl,-rpath -Wl,${HDF5_ROOT}/lib
endif
#===============================================================================
# Options for IBM Blue Gene/P ANL supercomputer
#===============================================================================
@ -190,7 +196,7 @@ ifeq ($(MACHINE),bluegene)
endif
#===============================================================================
# Options for Cray XK6 ORNL Jaguar supercomputer
# Options for Cray XK6 ORNL Titan supercomputer
#===============================================================================
ifeq ($(MACHINE),crayxk6)
@ -224,7 +230,7 @@ neat:
.PHONY: all xml-fortran clean neat distclean
%.o: %.F90
$(F90) -Ixml-fortran -Ixml-fortran/templates $(F90FLAGS) -c $<
$(F90) $(F90FLAGS) -Ixml-fortran -Ixml-fortran/templates -c $<
#===============================================================================
# Dependencies

2
src/OBJECTS
View file

@ -13,10 +13,12 @@ endf.o \
endf_header.o \
energy_grid.o \
error.o \
finalize.o \
fission.o \
geometry.o \
geometry_header.o \
global.o \
hdf5_interface.o \
initialize.o \
intercycle.o \
interpolation.o \

3
src/constants.F90
View file

@ -22,7 +22,8 @@ module constants
real(8), parameter :: TINY_BIT = 1e-8_8
! User for precision in geometry
real(8), parameter :: FP_PRECISION = 1e-5_8
real(8), parameter :: FP_PRECISION = 1e-14_8
real(8), parameter :: FP_REL_PRECISION = 1e-5_8
! Maximum number of collisions/crossings
integer, parameter :: MAX_EVENTS = 10000

49
src/finalize.F90 Normal file
View file

@ -0,0 +1,49 @@
module finalize
use global
use output, only: print_runtime
use tally, only: write_tallies, tally_statistics
use timing, only: timer_stop
#ifdef HDF5
use hdf5_interface, only: hdf5_write_timing, hdf5_close_output
#endif
implicit none
contains
!===============================================================================
! FINALIZE_RUN
!===============================================================================
subroutine finalize_run()
! Calculate statistics for tallies and write to tallies.out
if (.not. plotting) then
call tally_statistics()
if (master) call write_tallies()
end if
! show timing statistics
call timer_stop(time_total)
if (master) call print_runtime()
#ifdef HDF5
if (master) then
call hdf5_write_timing()
call hdf5_close_output()
end if
#endif
! deallocate arrays
call free_memory()
#ifdef MPI
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(mpi_err)
#endif
end subroutine finalize_run
end module finalize

120
src/geometry.F90
View file

@ -149,7 +149,6 @@ contains
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
! Create new level of coordinates
p % in_lower_universe = .true.
allocate(p % coord % next)
! Move particle to next level and set universe
@ -187,7 +186,6 @@ contains
else
! Create new level of coordinates
p % in_lower_universe = .true.
allocate(p % coord % next)
! adjust local position of particle
@ -284,7 +282,7 @@ contains
! PARTICLE REFLECTS FROM SURFACE
! Do not handle reflective boundary conditions on lower universes
if (p % in_lower_universe) then
if (.not. associated(p % coord, p % coord0)) then
message = "Cannot reflect particle " // trim(to_str(p % id)) // &
" off surface in a lower universe."
call fatal_error()
@ -307,11 +305,14 @@ contains
select case (surf%type)
case (SURF_PX)
p % coord0 % uvw = (/ -u, v, w /)
u = -u
case (SURF_PY)
p % coord0 % uvw = (/ u, -v, w /)
v = -v
case (SURF_PZ)
p % coord0 % uvw = (/ u, v, -w /)
w = -w
case (SURF_PLANE)
! Find surface coefficients and norm of vector normal to surface
n1 = surf % coeffs(1)
@ -325,8 +326,6 @@ contains
v = v - 2*dot_prod*n2/norm
w = w - 2*dot_prod*n3/norm
! Set vector
p % coord0 % uvw = (/ u, v, w /)
case (SURF_CYL_X)
! Find y-y0, z-z0 and dot product of direction and surface normal
y = p % coord0 % xyz(2) - surf % coeffs(1)
@ -338,8 +337,6 @@ contains
v = v - 2*dot_prod*y/(R*R)
w = w - 2*dot_prod*z/(R*R)
! Set vector
p % coord0 % uvw = (/ u, v, w /)
case (SURF_CYL_Y)
! Find x-x0, z-z0 and dot product of direction and surface normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
@ -351,8 +348,6 @@ contains
u = u - 2*dot_prod*x/(R*R)
w = w - 2*dot_prod*z/(R*R)
! Set vector
p % coord0 % uvw = (/ u, v, w /)
case (SURF_CYL_Z)
! Find x-x0, y-y0 and dot product of direction and surface normal
x = p % coord0 % xyz(1) - surf % coeffs(1)
@ -364,8 +359,6 @@ contains
u = u - 2*dot_prod*x/(R*R)
v = v - 2*dot_prod*y/(R*R)
! Set vector
p % coord0 % uvw = (/ u, v, w /)
case (SURF_SPHERE)
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
! normal
@ -380,18 +373,32 @@ contains
v = v - 2*dot_prod*y/(R*R)
w = w - 2*dot_prod*z/(R*R)
! Set vector
p % coord0 % uvw = (/ u, v, w /)
case default
message = "Reflection not supported for surface " // &
trim(to_str(surf % id))
call fatal_error()
end select
! Set new particle direction
p % coord0 % uvw = (/ u, v, w /)
! Reassign particle's cell and surface
p % coord0 % cell = last_cell
p % surface = -p % surface
! If a reflective surface is coincident with a lattice or universe
! boundary, it is necessary to redetermine the particle's coordinates in
! the lower universes.
if (associated(p % coord0 % next)) then
call deallocate_coord(p % coord0 % next)
call find_cell(p, found)
if (.not. found) then
message = "Couldn't find particle after reflecting from surface."
call fatal_error()
end if
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
@ -495,29 +502,7 @@ contains
! Check to make sure still in lattice
i_x = p % coord % lattice_x
i_y = p % coord % lattice_y
if (i_x < 1 .or. i_x > lat % n_x) then
message = "Particle " // trim(to_str(p % id)) // " reached edge of &
&lattice " // trim(to_str(lat % id)) // " at position (" // &
trim(to_str(i_x)) // "," // trim(to_str(i_y)) // ")."
call fatal_error()
elseif (i_y < 1 .or. i_y > lat % n_y) then
message = "Particle " // trim(to_str(p % id)) // " reached edge of &
&lattice " // trim(to_str(lat % id)) // " at position (" // &
trim(to_str(i_x)) // "," // trim(to_str(i_y)) // ")."
call fatal_error()
end if
! Find universe for next lattice element
p % coord % universe = lat % element(i_x, i_y)
! Find cell in next lattice element
call find_cell(p, found)
if (.not. found) then
! In some circumstances, a particle crossing the corner of a cell may not
! be able to be found in the next universe. In this scenario we cut off
! all lower-level coordinates and search from universe zero
! Remove lower coordinates
if (i_x < 1 .or. i_x > lat % n_x .or. i_y < 1 .or. i_y > lat % n_y) then
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
@ -528,6 +513,29 @@ contains
" after crossing a lattice boundary."
call fatal_error()
end if
else
! Find universe for next lattice element
p % coord % universe = lat % element(i_x, i_y)
! Find cell in next lattice element
call find_cell(p, found)
if (.not. found) then
! In some circumstances, a particle crossing the corner of a cell may not
! be able to be found in the next universe. In this scenario we cut off
! all lower-level coordinates and search from universe zero
! Remove lower coordinates
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
! Search for particle
call find_cell(p, found)
if (.not. found) then
message = "Could not locate particle " // trim(to_str(p % id)) // &
" after crossing a lattice boundary."
call fatal_error()
end if
end if
end if
end subroutine cross_lattice
@ -852,10 +860,12 @@ contains
! Check is calculated distance is new minimum
if (d < dist) then
dist = d
surface_crossed = -cl % surfaces(i)
lattice_crossed = NONE
final_coord => coord
if (abs(d - dist)/dist >= FP_PRECISION) then
dist = d
surface_crossed = -cl % surfaces(i)
lattice_crossed = NONE
final_coord => coord
end if
end if
end do SURFACE_LOOP
@ -872,16 +882,16 @@ contains
z = coord % xyz(3)
! determine oncoming edge
x0 = lat % width_x * 0.5_8
y0 = lat % width_y * 0.5_8
x0 = sign(lat % width_x * 0.5_8, u)
y0 = sign(lat % width_y * 0.5_8, v)
! left and right sides
if (u == ZERO) then
if (abs(x - x0) < FP_PRECISION) then
d = INFINITY
elseif (u == ZERO) then
d = INFINITY
elseif (u > 0) then
d = (x0 - x)/u
else
d = -(x + x0)/u
d = (x0 - x)/u
end if
! If the lattice boundary is coincident with the parent cell boundary,
@ -892,7 +902,7 @@ contains
! point precision.
if (d < dist) then
if (abs(d - dist)/dist >= FP_PRECISION) then
if (abs(d - dist)/dist >= FP_REL_PRECISION) then
dist = d
if (u > 0) then
lattice_crossed = LATTICE_RIGHT
@ -904,16 +914,16 @@ contains
end if
! top and bottom sides
if (v == ZERO) then
if (abs(y - y0) < FP_PRECISION) then
d = INFINITY
elseif (v == ZERO) then
d = INFINITY
elseif (v > 0) then
d = (y0 - y)/v
else
d = -(y + y0)/v
d = (y0 - y)/v
end if
if (d < dist) then
if (abs(d - dist)/dist >= FP_PRECISION) then
if (abs(d - dist)/dist >= FP_REL_PRECISION) then
dist = d
if (v > 0) then
lattice_crossed = LATTICE_TOP

11
src/global.F90
View file

@ -18,6 +18,10 @@ module global
use mpi
#endif
#ifdef HDF5
use hdf5
#endif
implicit none
save
@ -164,6 +168,13 @@ module global
real(8) :: plot_basis(6)
real(8) :: pixel
! ============================================================================
! HDF5 VARIABLES
#ifdef HDF5
integer(HID_T) :: hdf5_output_file
#endif
! ============================================================================
! MISCELLANEOUS VARIABLES

159
src/hdf5_interface.F90 Normal file
View file

@ -0,0 +1,159 @@
module hdf5_interface
use constants
use global
#ifdef HDF5
use hdf5
use h5lt
#endif
implicit none
contains
#ifdef HDF5
!===============================================================================
! HDF5_OPEN_OUTPUT
!===============================================================================
subroutine hdf5_open_output()
character(9), parameter :: filename = "output.h5" ! File name
integer :: error ! Error flag
! Initialize FORTRAN interface.
call h5open_f (error)
! Create a new file using default properties.
call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_output_file, error)
end subroutine hdf5_open_output
!===============================================================================
! HDF5_WRITE_SUMMARY
!===============================================================================
subroutine hdf5_write_summary()
integer :: error
integer :: rank = 1
integer(HSIZE_T) :: dims(1) = (/1/)
character(8) :: date_
character(10) :: time_
character(19) :: current_time
! Write version information
call h5ltmake_dataset_int_f(hdf5_output_file, "version_major", &
rank, dims, (/ VERSION_MAJOR /), error)
call h5ltmake_dataset_int_f(hdf5_output_file, "version_minor", &
rank, dims, (/ VERSION_MINOR /), error)
call h5ltmake_dataset_int_f(hdf5_output_file, "version_release", &
rank, dims, (/ VERSION_RELEASE /), error)
! Write current date and time
call date_and_time(DATE=date_, TIME=time_)
current_time = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) // &
" " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
call h5ltmake_dataset_string_f(hdf5_output_file, "/date_and_time", &
current_time, error)
! Write MPI information
call h5ltmake_dataset_int_f(hdf5_output_file, "n_procs", &
rank, dims, (/ n_procs /), error)
call h5ltset_attribute_string_f(hdf5_output_file, "n_procs", &
"description", "Number of MPI processes", error)
end subroutine hdf5_write_summary
!===============================================================================
! HDF5_WRITE_TIMING
!===============================================================================
subroutine hdf5_write_timing()
integer :: error
integer :: rank = 1
integer(HSIZE_T) :: dims(1) = (/1/)
integer(HID_T) :: timing_group
! Create group for timing
call h5gcreate_f(hdf5_output_file, "/timing", timing_group, error)
! Write timing data
call h5ltmake_dataset_double_f(timing_group, "time_initialize", &
rank, dims, (/ time_initialize % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_read_xs", &
rank, dims, (/ time_read_xs % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_unionize", &
rank, dims, (/ time_unionize % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_compute", &
rank, dims, (/ time_compute % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_intercycle", &
rank, dims, (/ time_intercycle % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_tallies", &
rank, dims, (/ time_ic_tallies % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_sample", &
rank, dims, (/ time_ic_sample % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_sendrecv", &
rank, dims, (/ time_ic_sendrecv % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_rebuild", &
rank, dims, (/ time_ic_rebuild % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_inactive", &
rank, dims, (/ time_inactive % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_active", &
rank, dims, (/ time_active % elapsed /), error)
call h5ltmake_dataset_double_f(timing_group, "time_total", &
rank, dims, (/ time_total % elapsed /), error)
! Add descriptions to timing data
call h5ltset_attribute_string_f(timing_group, "time_initialize", &
"description", "Total time elapsed for initialization (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_read_xs", &
"description", "Time reading cross-section libraries (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_unionize", &
"description", "Time unionizing energy grid (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_compute", &
"description", "Total time in computation (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_intercycle", &
"description", "Total time between cycles (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_tallies", &
"description", "Time between cycles accumulating tallies (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_sample", &
"description", "Time between cycles sampling source sites (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_sendrecv", &
"description", "Time between cycles SEND/RECVing source sites (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_rebuild", &
"description", "Time between cycles reconstructing source bank (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_inactive", &
"description", "Total time in inactive cycles (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_active", &
"description", "Total time in active cycles (s)", error)
call h5ltset_attribute_string_f(timing_group, "time_total", &
"description", "Total time elapsed (s)", error)
! Close timing group
call h5gclose_f(timing_group, error)
end subroutine hdf5_write_timing
!===============================================================================
! HDF5_CLOSE_OUTPUT
!===============================================================================
subroutine hdf5_close_output()
integer :: error ! Error flag
! Terminate access to the file.
call h5fclose_f(hdf5_output_file, error)
! Close FORTRAN interface.
call h5close_f(error)
end subroutine hdf5_close_output
#endif
end module hdf5_interface

25
src/initialize.F90
View file

@ -25,6 +25,10 @@ module initialize
use mpi
#endif
#ifdef HDF5
use hdf5_interface, only: hdf5_open_output, hdf5_write_summary
#endif
implicit none
type(DictionaryII), pointer :: build_dict => null()
@ -42,7 +46,8 @@ contains
type(Universe), pointer :: univ
! Start initialization timer
! Start total and initialization timer
call timer_start(time_total)
call timer_start(time_initialize)
! Setup MPI
@ -50,13 +55,21 @@ contains
! Read command line arguments
call read_command_line()
if (master) call create_summary_file()
! Print the OpenMC title and version/date/time information
if (master) call title()
if (master) then
! Create summary.out file
call create_summary_file()
! Print initialization header block
if (master) call header("INITIALIZATION", level=1)
#ifdef HDF5
! Open HDF5 output file for writing
call hdf5_open_output()
call hdf5_write_summary()
#endif
! Display title and initialization header
call title()
call header("INITIALIZATION", level=1)
end if
! Initialize random number generator
call initialize_prng()

26
src/input_xml.F90
View file

@ -2,14 +2,14 @@ module input_xml
use cmfd_utils, only: read_cmfd_xml
use constants
use datatypes, only: dict_create, dict_add_key, dict_has_key, &
use datatypes, only: dict_create, dict_add_key, dict_has_key, &
dict_get_key
use error, only: fatal_error, warning
use geometry_header, only: Cell, Surface, Lattice
use global
use mesh_header, only: StructuredMesh
use output, only: write_message
use string, only: lower_case, to_str, str_to_int, str_to_real, &
use string, only: lower_case, to_str, str_to_int, str_to_real, &
split_string, starts_with, ends_with
use tally_header, only: TallyObject
@ -131,6 +131,9 @@ contains
! Survival biasing
if (trim(survival_) == 'on') survival_biasing = .true.
! Probability tables
if (ptables_ == 'off') urr_ptables_on = .false.
! Cutoffs
if (size(cutoff_) > 0) then
weight_cutoff = cutoff_(1) % weight
@ -495,7 +498,13 @@ contains
! Copy density -- the default value for the units is given in the
! material_t.xml file and doesn't need to be specified here, hence case
! default results in an error.
val = material_(i) % density % value
val = material_(i) % density % value
if (val <= ZERO) then
message = "Need to specify a positive density on material " // &
trim(to_str(m % id)) // "."
call fatal_error()
end if
units = material_(i) % density % units
call lower_case(units)
select case(trim(units))
@ -1025,8 +1034,15 @@ contains
! Parse cross_sections.xml file
call read_xml_file_cross_sections_t(path_cross_sections)
! Copy directory information if present
directory = trim(directory_)
if (len_trim(directory_) > 0) then
! Copy directory information if present
directory = trim(directory_)
else
! If no directory is listed in cross_sections.xml, by default select the
! directory in which the cross_sections.xml file resides
i = index(path_cross_sections, "/", BACK=.true.)
directory = path_cross_sections(1:i)
end if
! determine whether binary/ascii
if (filetype_ == 'ascii') then

29
src/main.F90
View file

@ -4,17 +4,17 @@ program main
use cmfd_utils, only: read_hdf5
use constants
use global
use finalize, only: finalize_run
use initialize, only: initialize_run
use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank
use output, only: write_message, header, print_runtime
use output, only: write_message, header
use particle_header, only: Particle
use plot, only: run_plot
use physics, only: transport
use random_lcg, only: set_particle_seed
use source, only: get_source_particle
use string, only: to_str
use tally, only: synchronize_tallies, write_tallies, &
tally_statistics
use tally, only: synchronize_tallies
use timing, only: timer_start, timer_stop
#ifdef MPI
@ -23,6 +23,7 @@ program main
implicit none
<<<<<<< HEAD
integer :: nargs ! number of arguments
character(len=32) :: arg ! the argument
@ -57,32 +58,18 @@ program main
case default
! start timer for total run time
call timer_start(time_total)
! set up problem
call initialize_run()
! start problem
if (plotting) then
call run_plot()
call run_plot()
else
call run_problem()
! Calculate statistics for tallies and write to tallies.out
call tally_statistics()
if (master) call write_tallies()
! show timing statistics
call timer_stop(time_total)
if (master) call print_runtime()
call run_problem()
end if
! call cmfd calculation
call execute_cmfd()
! deallocate arrays
call free_memory()
! finalize run
call finalize_run()
end select

2
src/output.F90
View file

@ -261,7 +261,7 @@ contains
! Print surface
if (p % surface /= NONE) then
s => surfaces(p % surface)
write(ou,*) ' Surface = ' // to_str(s % id)
write(ou,*) ' Surface = ' // to_str(s % id)
end if
write(ou,*) ' Weight = ' // to_str(p % wgt)

2
src/particle_header.F90
View file

@ -39,7 +39,6 @@ module particle_header
integer :: type ! Particle type (n, p, e, etc)
! Particle coordinates
logical :: in_lower_universe ! is particle in lower universe?
type(LocalCoord), pointer :: coord0 => null() ! coordinates on universe 0
type(LocalCoord), pointer :: coord => null() ! coordinates on lowest universe
@ -107,7 +106,6 @@ contains
allocate(p % coord0)
p % coord0 % universe = BASE_UNIVERSE
p % coord => p % coord0
p % in_lower_universe = .false.
end subroutine initialize_particle

7
src/physics.F90
View file

@ -45,6 +45,7 @@ contains
! Particle couldn't be located
if (.not. found_cell) then
message = "Could not locate particle " // trim(to_str(p % id))
print *, p % coord0 % xyz
call fatal_error()
end if
@ -79,7 +80,11 @@ contains
call distance_to_boundary(p, d_boundary, surface_crossed, lattice_crossed)
! Sample a distance to collision
d_collision = -log(prn()) / material_xs % total
if (material_xs % total == ZERO) then
d_collision = INFINITY
else
d_collision = -log(prn()) / material_xs % total
end if
! Select smaller of the two distances
distance = min(d_boundary, d_collision)

40
src/utils/core_map.py
View file

@ -1,40 +0,0 @@
#!/usr/bin/env python
import sys
import numpy as np
import matplotlib.pyplot as pyplot
# Get filename and file object
filename = sys.argv[1]
axial_level = int(sys.argv[2])
if len(sys.argv) > 3:
stage = int(sys.argv[3])
else:
stage = 1
fh = open(filename,'r')
# Read size of mesh
words = fh.readline().split()
nx, ny, nz, nstage = [int(item) for item in words]
# Read values
value_dict = {}
for line in fh:
words = line.split()
i, j, k, m = [int(item) for item in words[:4]]
value = float(words[-1])
value_dict[i,j,k,m] = value
# Set up matrix
matrix = np.array([[value_dict[i+1,j+1,axial_level,stage] for i in range(nx)]
for j in range(ny)])
# Make figure
pyplot.pcolor(matrix)
pyplot.colorbar()
pyplot.xlim([0,nx])
pyplot.ylim([0,ny])
pyplot.xticks([])
pyplot.yticks([])
pyplot.show()

2
src/xml-fortran/templates/materials_t.xml
View file

@ -8,7 +8,7 @@
<!-- Type for specifying density of a material -->
<typedef name="density_xml">
<component name="value" type="double" />
<component name="value" type="double" default="0.0" />
<component name="units" type="word" length="10" default="'atom/b-cm'" />
</typedef>

1
src/xml-fortran/templates/settings_t.xml
View file

@ -23,6 +23,7 @@
<variable name="verbosity_" tag="verbosity" type="integer" />
<variable name="source_" tag="source" type="source_xml" />
<variable name="survival_" tag="survival_biasing" type="word" length="3" />
<variable name="ptables_" tag="ptables" type="word" length="3" />
<variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
<variable name="cross_sections_" tag="cross_sections" type="word" length="255" />
<variable name="trace_" tag="trace" type="integer-array" />