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

Cleaned up source, added comments, removed unnecessary variables.

Browse source
This commit is contained in:
Paul Romano 2012-03-24 11:19:39 -04:00
parent 2c24c3d405
commit 0adb96d042
18 changed files with 245 additions and 270 deletions
Download patch file Download diff file Expand all files Collapse all files

4
src/DEPENDENCIES
View file

@ -99,6 +99,7 @@ initialize.o: random_lcg.o
initialize.o: source.o
initialize.o: string.o
initialize.o: tally.o
initialize.o: tally_header.o
initialize.o: timing.o
input_xml.o: constants.o
@ -143,7 +144,6 @@ main.o: intercycle.o
main.o: output.o
main.o: physics.o
main.o: plot.o
main.o: random_lcg.o
main.o: source.o
main.o: string.o
main.o: tally.o
@ -177,6 +177,7 @@ physics.o: geometry.o
physics.o: geometry_header.o
physics.o: global.o
physics.o: interpolation.o
physics.o: material_header.o
physics.o: output.o
physics.o: particle_header.o
physics.o: random_lcg.o
@ -198,7 +199,6 @@ plot.o: ppmlib.o
plot.o: source.o
plot.o: string.o
search.o: constants.o
search.o: error.o
search.o: global.o

7
src/ace.F90
View file

@ -22,8 +22,6 @@ module ace
real(8), allocatable :: XSS(:) ! Cross section data
integer :: XSS_index ! current index in XSS data
type(DictionaryCI), pointer :: already_read => null()
private :: NXS
private :: JXS
private :: XSS
@ -47,6 +45,7 @@ contains
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: nuc => null()
type(SAB_Table), pointer :: sab => null()
type(DictionaryCI), pointer :: already_read => null()
! ==========================================================================
! COUNT NUMBER OF TABLES AND CREATE DICTIONARIES
@ -154,6 +153,8 @@ contains
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
! Read the ACE table into the appropriate entry on the nuclides
! array
call read_ace_table(index_nuclides, index_list)
! Print out information on table to cross_sections.out file
@ -172,6 +173,8 @@ contains
index_list = dict_get_key(xs_listing_dict, name)
index_sab = dict_get_key(sab_dict, name)
! Read the ACE table into the appropriate entry on the sab_tables
! array
call read_ace_table(index_sab, index_list)
call dict_add_key(already_read, name, 0)

4
src/bank_header.F90
View file

@ -9,7 +9,11 @@ module bank_header
!===============================================================================
type Bank
! The 'sequence' attribute is used here to ensure that the data listed
! appears in the given order. This is important for MPI purposes when bank
! sites are sent from one processor to another.
sequence
integer(8) :: id ! Unique ID
real(8) :: xyz(3) ! location of bank particle
real(8) :: uvw(3) ! diretional cosines

2
src/energy_grid.F90
View file

@ -37,7 +37,7 @@ contains
! Create dictionary for keeping track of cross sections already added
call dict_create(already_added)
! Loop over all files
! Loop over all materials
do i = 1, n_materials
mat => materials(i)

14
src/initialize.F90
View file

@ -3,13 +3,13 @@ module initialize
use ace, only: read_xs
use bank_header, only: Bank
use constants
use datatypes, only: dict_create, dict_add_key, dict_get_key, &
dict_has_key, dict_keys
use datatypes, only: dict_create, dict_get_key, dict_has_key, &
dict_keys
use datatypes_header, only: ListKeyValueII, DictionaryII
use energy_grid, only: unionized_grid
use error, only: fatal_error
use geometry, only: neighbor_lists
use geometry_header, only: Cell, Surface, Universe, Lattice, BASE_UNIVERSE
use geometry_header, only: Cell, Universe, Lattice, BASE_UNIVERSE
use global
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict, read_plots_xml
@ -19,7 +19,8 @@ module initialize
use random_lcg, only: initialize_prng
use source, only: initialize_source
use string, only: to_str, starts_with, ends_with, lower_case
use tally, only: create_tally_map, TallyObject
use tally, only: create_tally_map
use tally_header, only: TallyObject
use timing, only: timer_start, timer_stop
#ifdef MPI
@ -33,8 +34,6 @@ module initialize
implicit none
type(DictionaryII), pointer :: build_dict => null()
contains
!===============================================================================
@ -286,9 +285,6 @@ contains
! Create special dictionary used in input_xml
call dict_create(cells_in_univ_dict)
! Create special dictionary for building universes
call dict_create(build_dict)
end subroutine create_dictionaries

3
src/input_xml.F90
View file

@ -1,8 +1,7 @@
module input_xml
use constants
use datatypes, only: dict_create, dict_add_key, dict_has_key, &
dict_get_key
use datatypes, only: dict_add_key, dict_has_key, dict_get_key
use error, only: fatal_error, warning
use geometry_header, only: Cell, Surface, Lattice
use global

1
src/main.F90
View file

@ -8,7 +8,6 @@ program main
use output, only: write_message, header
use plotter, 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, add_to_score

14
src/mesh_header.F90
View file

@ -8,13 +8,13 @@ module mesh_header
!===============================================================================
type StructuredMesh
integer :: id
integer :: type
integer :: n_dimension
integer, allocatable :: dimension(:)
real(8), allocatable :: lower_left(:)
real(8), allocatable :: upper_right(:)
real(8), allocatable :: width(:)
integer :: id ! user-specified id
integer :: type ! rectangular, hexagonal
integer :: n_dimension ! rank of mesh
integer, allocatable :: dimension(:) ! number of cells in each direction
real(8), allocatable :: lower_left(:) ! lower-left corner of mesh
real(8), allocatable :: upper_right(:) ! upper-right corner of mesh
real(8), allocatable :: width(:) ! width of each mesh cell
end type StructuredMesh
end module mesh_header

190
src/output.F90
View file

@ -29,8 +29,8 @@ contains
subroutine title()
character(10) :: today_date
character(8) :: today_time
character(10) :: date_
character(10) :: time_
write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') &
' .d88888b. 888b d888 .d8888b.', &
@ -55,18 +55,23 @@ contains
#endif
! Write the date and time
call get_today(today_date, today_time)
call date_and_time(DATE=date_, TIME=time_)
date_ = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8)
time_ = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A)') &
trim(today_date), trim(today_time)
trim(date_), trim(time_)
! Write information to summary file
call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1)
write(UNIT=UNIT_SUMMARY, FMT=*) &
"Copyright: 2011 Massachusetts Institute of Technology"
"Copyright: 2011-2012 Massachusetts Institute of Technology"
write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') &
"Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1
#endif
write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A,1X,A)') &
trim(today_date), trim(today_time)
trim(date_), trim(time_)
! Write information on number of processors
#ifdef MPI
@ -89,8 +94,8 @@ contains
integer, optional :: unit ! unit to write to
integer, optional :: level ! specified header level
integer :: n
integer :: m
integer :: n ! number of = signs on left
integer :: m ! number of = signs on right
integer :: unit_ ! unit to write to
integer :: header_level ! actual header level
character(MAX_LINE_LEN) :: line
@ -140,14 +145,15 @@ contains
subroutine write_message(level)
integer, optional :: level
integer, optional :: level ! verbosity level
integer :: n_lines
integer :: i
integer :: n_lines ! number of lines needed
integer :: i ! index for lines
! Only allow master to print to screen
if (.not. master .and. present(level)) return
! TODO: Take care of line wrapping so words don't get cut off
if (.not. present(level) .or. level <= verbosity) then
n_lines = (len_trim(message)-1)/79 + 1
do i = 1, n_lines
@ -157,55 +163,13 @@ contains
end subroutine write_message
!===============================================================================
! GET_TODAY determines the date and time at which the program began execution
! and returns it in a readable format
!===============================================================================
subroutine get_today(today_date, today_time)
character(10), intent(out) :: today_date
character(8), intent(out) :: today_time
integer :: val(8)
character(8) :: date_
character(10) :: time_
character(5) :: zone
call date_and_time(date_, time_, zone, val)
! val(1) = year (YYYY)
! val(2) = month (MM)
! val(3) = day (DD)
! val(4) = timezone
! val(5) = hours (HH)
! val(6) = minutes (MM)
! val(7) = seconds (SS)
! val(8) = milliseconds
if (val(2) < 10) then
if (val(3) < 10) then
today_date = date_(6:6) // "/" // date_(8:8) // "/" // date_(1:4)
else
today_date = date_(6:6) // "/" // date_(7:8) // "/" // date_(1:4)
end if
else
if (val(3) < 10) then
today_date = date_(5:6) // "/" // date_(8:8) // "/" // date_(1:4)
else
today_date = date_(5:6) // "/" // date_(7:8) // "/" // date_(1:4)
end if
end if
today_time = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
end subroutine get_today
!===============================================================================
! PRINT_PARTICLE displays the attributes of a particle
!===============================================================================
subroutine print_particle()
integer :: i
integer :: i ! index for coordinate levels
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
@ -266,6 +230,7 @@ contains
write(ou,*) ' Surface = ' // to_str(sign(s % id, p % surface))
end if
! Display weight, energy, grid index, and interpolation factor
write(ou,*) ' Weight = ' // to_str(p % wgt)
write(ou,*) ' Energy = ' // to_str(p % E)
write(ou,*) ' IE = ' // to_str(p % IE)
@ -301,27 +266,35 @@ contains
subroutine print_cell(c, unit)
type(Cell), pointer :: c
integer, optional :: unit
integer, optional :: unit ! specified unit to write to
integer :: temp
integer :: i
integer :: unit_
integer :: index_cell ! index in cells array
integer :: i ! loop index for surfaces
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
type(Material), pointer :: m => null()
! Set unit to stdout if not already set
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write user-specified id for cell
write(unit_,*) 'Cell ' // to_str(c % id)
temp = dict_get_key(cell_dict, c % id)
write(unit_,*) ' Array Index = ' // to_str(temp)
! Find index in cells array and write
index_cell = dict_get_key(cell_dict, c % id)
write(unit_,*) ' Array Index = ' // to_str(index_cell)
! Write what universe this cell is in
u => universes(c % universe)
write(unit_,*) ' Universe = ' // to_str(u % id)
! Write information on fill for cell
select case (c % type)
case (CELL_NORMAL)
write(unit_,*) ' Fill = NONE'
@ -332,12 +305,16 @@ contains
l => lattices(c % fill)
write(unit_,*) ' Fill = Lattice ' // to_str(l % id)
end select
! Write information on material
if (c % material == 0) then
write(unit_,*) ' Material = NONE'
else
m => materials(c % material)
write(unit_,*) ' Material = ' // to_str(m % id)
end if
! Write surface specification
string = ""
do i = 1, c % n_surfaces
select case (c % surfaces(i))
@ -367,24 +344,31 @@ contains
type(Universe), pointer :: univ
integer, optional :: unit
integer :: i
integer :: unit_
integer :: i ! loop index for cells in this universe
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
type(Universe), pointer :: base_u
type(Universe), pointer :: base_u => null()
! Set default unit to stdout if not specified
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Get a pointer to the base universe
base_u => universes(BASE_UNIVERSE)
! Write user-specified id for this universe
write(unit_,*) 'Universe ' // to_str(univ % id)
! If this is the base universe, indicate so
if (associated(univ, base_u)) then
write(unit_,*) ' Base Universe'
end if
! Write list of cells in this universe
string = ""
do i = 1, univ % n_cells
c => cells(univ % cells(i))
@ -404,14 +388,16 @@ contains
type(Lattice), pointer :: lat
integer, optional :: unit
integer :: unit_
integer :: unit_ ! unit to write to
! set default unit if not specified
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write information about lattice
write(unit_,*) 'Lattice ' // to_str(lat % id)
write(unit_,*) ' n_x = ' // to_str(lat % n_x)
write(unit_,*) ' n_y = ' // to_str(lat % n_y)
@ -430,19 +416,23 @@ contains
subroutine print_surface(surf, unit)
type(Surface), pointer :: surf
integer, optional :: unit
integer, optional :: unit ! specified unit to write to
integer :: i
integer :: unit_
integer :: i ! loop index for coefficients
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
! set default unit if not specified
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write user-specified id of surface
write(unit_,*) 'Surface ' // to_str(surf % id)
! Write type of surface
select case (surf % type)
case (SURF_PX)
string = "X Plane"
@ -469,12 +459,14 @@ contains
end select
write(unit_,*) ' Type = ' // trim(string)
! Write coefficients for this surface
string = ""
do i = 1, size(surf % coeffs)
string = trim(string) // ' ' // to_str(surf % coeffs(i), 4)
end do
write(unit_,*) ' Coefficients = ' // trim(string)
! Write neighboring cells on positive side of this surface
string = ""
if (allocated(surf % neighbor_pos)) then
do i = 1, size(surf % neighbor_pos)
@ -483,6 +475,7 @@ contains
end if
write(unit_,*) ' Positive Neighbors = ' // trim(string)
! Write neighboring cells on negative side of this surface
string = ""
if (allocated(surf % neighbor_neg)) then
do i = 1, size(surf % neighbor_neg)
@ -490,6 +483,8 @@ contains
end do
end if
write(unit_,*) ' Negative Neighbors =' // trim(string)
! Write boundary condition for this surface
select case (surf % bc)
case (BC_TRANSMIT)
write(unit_,*) ' Boundary Condition = Transmission'
@ -513,12 +508,13 @@ contains
type(Material), pointer :: mat
integer, optional :: unit
integer :: i
integer :: unit_
real(8) :: density
integer :: i ! loop index for nuclides
integer :: unit_ ! unit to write to
real(8) :: density ! density in atom/b-cm
character(MAX_LINE_LEN) :: string
type(Nuclide), pointer :: nuc => null()
! set default unit to stdout if not specified
if (present(unit)) then
unit_ = unit
else
@ -559,9 +555,9 @@ contains
type(TallyObject), pointer :: t
integer, optional :: unit
integer :: i
integer :: id
integer :: unit_
integer :: i ! index for filter or score bins
integer :: id ! user-specified id
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
@ -569,14 +565,17 @@ contains
type(Material), pointer :: m => null()
type(StructuredMesh), pointer :: sm => null()
! set default unit to stdout if not specified
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write user-specified id of tally
write(unit_,*) 'Tally ' // to_str(t % id)
! Write any cells bins if present
if (t % n_filter_bins(FILTER_CELL) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_CELL)
@ -587,6 +586,7 @@ contains
write(unit_, *) ' Cell Bins:' // trim(string)
end if
! Write any surface bins if present
if (t % n_filter_bins(FILTER_SURFACE) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_SURFACE)
@ -597,6 +597,7 @@ contains
write(unit_, *) ' Surface Bins:' // trim(string)
end if
! Write any universe bins if present
if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_UNIVERSE)
@ -607,6 +608,7 @@ contains
write(unit_, *) ' Material Bins:' // trim(string)
end if
! Write any material bins if present
if (t % n_filter_bins(FILTER_MATERIAL) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_MATERIAL)
@ -617,6 +619,7 @@ contains
write(unit_, *) ' Material Bins:' // trim(string)
end if
! Write any mesh bins if present
if (t % n_filter_bins(FILTER_MESH) > 0) then
string = ""
id = t % mesh
@ -628,6 +631,7 @@ contains
write(unit_, *) ' Mesh Bins:' // trim(string)
end if
! Write any birth region bins if present
if (t % n_filter_bins(FILTER_CELLBORN) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_CELLBORN)
@ -638,6 +642,7 @@ contains
write(unit_, *) ' Birth Region Bins:' // trim(string)
end if
! Write any incoming energy bins if present
if (t % n_filter_bins(FILTER_ENERGYIN) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_ENERGYIN) + 1
@ -647,6 +652,7 @@ contains
write(unit_,*) ' Incoming Energy Bins:' // trim(string)
end if
! Write any outgoing energy bins if present
if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then
string = ""
do i = 1, t % n_filter_bins(FILTER_ENERGYOUT) + 1
@ -656,6 +662,7 @@ contains
write(unit_,*) ' Outgoing Energy Bins:' // trim(string)
end if
! Write any score bins if present
if (t % n_score_bins > 0) then
string = ""
do i = 1, t % n_score_bins
@ -674,7 +681,7 @@ contains
string = trim(string) // ' nu-fission'
end select
end do
write(unit_,*) ' Macro Reactions:' // trim(string)
write(unit_,*) ' Scores:' // trim(string)
end if
write(unit_,*)
@ -687,7 +694,7 @@ contains
subroutine print_geometry()
integer :: i
integer :: i ! loop index for various arrays
type(Surface), pointer :: s => null()
type(Cell), pointer :: c => null()
type(Universe), pointer :: u => null()
@ -735,12 +742,12 @@ contains
type(Nuclide), pointer :: nuc
integer, optional :: unit
integer :: i
integer :: unit_
integer :: size_total
integer :: size_xs
integer :: size_angle
integer :: size_energy
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_total ! memory used by nuclide (bytes)
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_angle ! memory used for angle distributions (bytes)
integer :: size_energy ! memory used for energy distributions (bytes)
type(Reaction), pointer :: rxn => null()
type(UrrData), pointer :: urr => null()
@ -824,7 +831,7 @@ contains
subroutine print_summary()
integer :: i
integer :: i ! loop index
character(15) :: string
type(Material), pointer :: m => null()
type(TallyObject), pointer :: t => null()
@ -896,24 +903,24 @@ contains
subroutine print_plot()
integer i
type(Plot), pointer :: pl => null()
integer :: i ! loop index for plots
type(Plot), pointer :: pl => null()
! Display header for plotting
call header("PLOTTING SUMMARY")
do i=1,n_plots
do i = 1, n_plots
pl => plots(i)
! Print plot id
! Write plot id
write(ou,100) "Plot ID:", trim(to_str(pl % id))
! Print plotting origin
! Write plotting origin
write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // &
" " // trim(to_str(pl % origin(2))) // " " // &
trim(to_str(pl % origin(3)))
! Print plotting width
! Write plotting width
if (pl % type == PLOT_TYPE_SLICE) then
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
@ -940,8 +947,8 @@ contains
subroutine print_runtime()
integer(8) :: total_particles
real(8) :: speed
integer(8) :: total_particles ! total # of particles simulated
real(8) :: speed ! # of neutrons/second
character(15) :: string
! display header block
@ -1017,7 +1024,8 @@ contains
end subroutine create_summary_file
!===============================================================================
! CREATE_XS_SUMMARY_FILE
! CREATE_XS_SUMMARY_FILE creates an output file to write information about the
! cross section tables used in the simulation.
!===============================================================================
subroutine create_xs_summary_file()

76
src/physics.F90
View file

@ -3,7 +3,7 @@ module physics
use ace_header, only: Nuclide, Reaction, DistEnergy
use constants
use cross_section, only: calculate_xs
use endf, only: reaction_name, is_fission, is_scatter
use endf, only: reaction_name
use error, only: fatal_error, warning
use fission, only: nu_total, nu_prompt, nu_delayed
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
@ -11,6 +11,7 @@ module physics
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use interpolation, only: interpolate_tab1
use material_header, only: Material
use output, only: write_message
use particle_header, only: LocalCoord
use random_lcg, only: prn
@ -2061,77 +2062,4 @@ contains
end function watt_spectrum
!===============================================================================
! WIGNER samples a Wigner distribution of energy level spacings. Note that this
! scheme is C50 in the Monte Carlo Sampler from Los Alamos (LA-9721-MS).
!===============================================================================
function wigner(D_avg) result (D)
real(8), intent(in) :: D_avg ! average level spacing
real(8) :: D ! sampled level spacing
real(8) :: c
c = -4.*D_avg*D_avg/PI * log(prn())
D = sqrt(c)
end function wigner
!===============================================================================
! CHI_SQUARED samples a chi-squared distribution with n degrees of freedom. The
! distribution of resonance widths in the unresolved region is given by a
! chi-squared distribution. For the special case of n=1, this is a Porter-Thomas
! distribution. For cases with n odd, rule C64 is used whereas for cases with n
! even, rule C45 is used.
!===============================================================================
function chi_squared(n, G_avg) result(G)
integer, intent(in) :: n ! number of degrees of freedom
real(8), intent(in), optional :: G_avg ! average resonance width
integer :: i ! loop index
real(8) :: G ! sampled random variable (or resonance width)
real(8) :: x, y, c ! dummy variables
real(8) :: r1, r2 ! psuedorandom numbers
select case (mod(n,2))
case (0)
! Even number of degrees of freedom can be sampled via rule C45. We can
! sample x as -2/n*log(product(r_i, i = 1 to n/2))
x = ONE
do i = 1, n/2
x = x * prn()
end do
x = -2./n * log(x)
case (1)
! Odd number of degrees of freedom can be sampled via rule C64. We can
! sample x as -2/n*(log(r)*cos^2(pi/2*r) + log(product(r_i, i = 1 to
! floor(n/2)))
! Note that we take advantage of integer division on n/2
y = ONE
do i = 1, n/2
y = y * prn()
end do
r1 = prn()
r2 = prn()
c = cos(PI/2.*r2)
x = -2./n * (log(y) + log(r1)*c*c)
end select
! If sampling a chi-squared distribution for a resonance width and the
! average resonance width has been given, return the sampled resonance
! width.
if (present(G_avg)) then
G = x * G_avg
else
G = x
end if
end function chi_squared
end module physics

111
src/plot.F90
View file

@ -2,12 +2,11 @@ module plotter
use constants
use error, only: fatal_error
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
cross_lattice, cell_contains
use geometry, only: find_cell
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord
use particle_header, only: deallocate_coord
use plot_header
use ppmlib, only: Image, init_image, allocate_image, &
deallocate_image, set_pixel
@ -18,17 +17,16 @@ module plotter
contains
!===============================================================================
! RUN_PLOT
! RUN_PLOT controls the logic for making one or many plots
!===============================================================================
subroutine run_plot()
integer :: i
type(Plot), pointer :: pl => null()
integer :: i ! loop index for plots
type(Plot), pointer :: pl => null()
do i=1,n_plots
do i = 1, n_plots
pl => plots(i)
! Display output message
@ -36,31 +34,34 @@ contains
call write_message(5)
if (pl % type == PLOT_TYPE_SLICE) then
! create 2d image
call create_ppm(pl)
end if
end do
end subroutine run_plot
!===============================================================================
! create_ppm
! CREATE_PPM creates an image based on user input from a plots.xml <plot>
! specification in the portable pixmap format (PPM)
!===============================================================================
subroutine create_ppm(pl)
type(Plot), pointer :: pl
type(Plot), pointer :: pl
type(Image) :: img
integer :: in_i, out_i
integer :: x, y
integer :: r, g, b
real(8) :: in_pixel, out_pixel
real(8) :: xyz(3)
logical :: found_cell
type(Cell), pointer :: c => null()
integer :: in_i
integer :: out_i
integer :: x, y ! pixel location
integer :: r, g, b ! colors (red, green, blue) from 0-255
real(8) :: in_pixel
real(8) :: out_pixel
real(8) :: xyz(3)
logical :: found_cell
type(Image) :: img
type(Cell), pointer :: c => null()
! Initialize and allocate space for image
call init_image(img)
call allocate_image(img, pl % pixels(1), pl % pixels(2))
@ -97,8 +98,8 @@ contains
p % coord % uvw = (/ 1, 0, 0 /)
p % coord % universe = BASE_UNIVERSE
do y=1, img%height
do x=1, img%width
do y = 1, img % height
do x = 1, img % width
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
@ -106,63 +107,74 @@ contains
call find_cell(found_cell)
if (.not. found_cell) then
r = pl % not_found % rgb(1)
g = pl % not_found % rgb(2)
b = pl % not_found % rgb(3)
! If no cell, revert to default color
r = pl % not_found % rgb(1)
g = pl % not_found % rgb(2)
b = pl % not_found % rgb(3)
else
if (pl % color_by == PLOT_COLOR_MATS) then
c => cells(p % coord % cell)
r = pl % colors(c % material) % rgb(1)
g = pl % colors(c % material) % rgb(2)
b = pl % colors(c % material) % rgb(3)
else if (pl % color_by == PLOT_COLOR_CELLS) then
r = pl % colors(p % coord % cell) % rgb(1)
g = pl % colors(p % coord % cell) % rgb(2)
b = pl % colors(p % coord % cell) % rgb(3)
else
r = 0
g = 0
b = 0
end if
if (pl % color_by == PLOT_COLOR_MATS) then
! Assign color based on material
c => cells(p % coord % cell)
r = pl % colors(c % material) % rgb(1)
g = pl % colors(c % material) % rgb(2)
b = pl % colors(c % material) % rgb(3)
else if (pl % color_by == PLOT_COLOR_CELLS) then
! Assign color based on cell
r = pl % colors(p % coord % cell) % rgb(1)
g = pl % colors(p % coord % cell) % rgb(2)
b = pl % colors(p % coord % cell) % rgb(3)
else
r = 0
g = 0
b = 0
end if
end if
! Create a pixel at (x,y) with color (r,g,b)
call set_pixel(img, x, y, r, g, b)
! Advance pixel in first direction
p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel
end do
! Advance pixel in second direction
p % coord0 % xyz(in_i) = xyz(in_i)
p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel
end do
! Write out the ppm to a file
call output_ppm(pl,img)
! Free up space
call deallocate_image(img)
end subroutine create_ppm
!===============================================================================
! output_ppm
! OUTPUT_PPM writes out a previously generated image to a PPM file
!===============================================================================
subroutine output_ppm(pl,img)
type(Plot), pointer :: pl
type(Image), intent(in) :: img
subroutine output_ppm(pl, img)
integer :: i, j
type(Plot), pointer :: pl
type(Image), intent(in) :: img
integer :: i ! loop index for height
integer :: j ! loop index for width
! Open PPM file for writing
open(UNIT=UNIT_PLOT, FILE=pl % path_plot)
! Write header
write(UNIT_PLOT, '(A2)') 'P6'
write(UNIT_PLOT, '(I0,'' '',I0)') img%width, img%height
write(UNIT_PLOT, '(A)') '255'
do j=1, img%height
do i=1, img%width
! Write color for each pixel
do j = 1, img % height
do i = 1, img % width
write(UNIT_PLOT, '(3A1)', advance='no') achar(img%red(i,j)), &
achar(img%green(i,j)), &
achar(img%blue(i,j))
achar(img%green(i,j)), achar(img%blue(i,j))
end do
end do
@ -171,5 +183,4 @@ contains
end subroutine output_ppm
end module plotter

2
src/plot_header.F90
View file

@ -26,7 +26,7 @@ module plot_header
integer :: basis ! direction of plot slice
integer :: pixels(2) ! pixel width/height of plot slice
type(ObjectColor) :: not_found ! color for positions where no cell found
type(ObjectColor),allocatable :: colors(:) ! colors of cells/mats
type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats
end type Plot
integer :: PLOT_TYPE_SLICE = 1

66
src/ppmlib.F90
View file

@ -13,6 +13,10 @@ module ppmlib
contains
!===============================================================================
! INIT_IMAGE initializes the Image derived type
!===============================================================================
subroutine init_image(img)
type(Image), intent(out) :: img
@ -26,29 +30,46 @@ contains
end subroutine init_image
!===============================================================================
! ALLOCATE_IMAGE sets the width and height of an image and allocates color
! arrays
!===============================================================================
subroutine allocate_image(img, w, h)
type(Image), intent(inout) :: img
integer, intent(in) :: w, h
integer, intent(in) :: w ! width of image
integer, intent(in) :: h ! height of image
! allocate red, green, and blue array
allocate(img % red(w, h))
allocate(img % green(w, h))
allocate(img % blue(w, h))
! set width and height
img % width = w
img % height = h
end subroutine allocate_image
!===============================================================================
! DEALLOCATE_IMAGE
!===============================================================================
subroutine deallocate_image(img)
type(Image) :: img
if ( associated(img % red) ) deallocate(img % red)
if ( associated(img % green) ) deallocate(img % green)
if ( associated(img % blue) ) deallocate(img % blue)
if (associated(img % red)) deallocate(img % red)
if (associated(img % green)) deallocate(img % green)
if (associated(img % blue)) deallocate(img % blue)
end subroutine deallocate_image
!===============================================================================
! INSIDE_IMAGE determines whether a point (x,y) is inside the image
!===============================================================================
function inside_image(img, x, y) result(inside)
@ -58,13 +79,15 @@ contains
inside = .false.
if ( ( x < img % width) .and. &
( y < img % height ) .and. &
( x >= 0 ) .and. &
( y >= 0 ) ) inside = .true.
if ((x < img % width) .and. (y < img % height) .and. &
(x >= 0) .and. (y >= 0)) inside = .true.
end function inside_image
!===============================================================================
! VALID_IMAGE checks whether the image has a width and height and if its color
! arrays are allocated
!===============================================================================
function valid_image(img) result(valid)
@ -73,29 +96,32 @@ contains
valid = .false.
if ( img % width == 0 ) return
if ( img % height == 0 ) return
if ( .not. associated(img % red) .or. &
if (img % width == 0) return
if (img % height == 0) return
if (.not. associated(img % red) .or. &
.not. associated(img % green) .or. &
.not. associated(img % blue) ) return
.not. associated(img % blue)) return
valid = .true.
end function valid_image
!===============================================================================
! SET_PIXEL sets the
!===============================================================================
subroutine set_pixel(img, x, y, r, g, b)
type(Image), intent(inout) :: img
integer, intent(in) :: x, y
integer, intent(in) :: r, g, b
integer, intent(in) :: x, y ! coordinates
integer, intent(in) :: r, g, b ! red, green, and blue
if ( inside_image(img, x, y) .and. valid_image(img)) then
img % red(x+1,y+1) = mod(abs(r), 256)
img % green(x+1, y+1) = mod(abs(g), 256)
img % blue(x+1, y+1) = mod(abs(b), 256)
if (inside_image(img, x, y) .and. valid_image(img)) then
img % red(x+1,y+1) = mod(abs(r), 256)
img % green(x+1, y+1) = mod(abs(g), 256)
img % blue(x+1, y+1) = mod(abs(b), 256)
end if
end subroutine set_pixel
end subroutine set_pixel
end module ppmlib

1
src/search.F90
View file

@ -1,6 +1,5 @@
module search
use constants, only: ONE, MAX_LINE_LEN
use error, only: fatal_error
use global, only: message

2
src/source.F90
View file

@ -1,7 +1,7 @@
module source
use bank_header, only: Bank
use constants, only: ONE
use constants
use error, only: fatal_error
use geometry_header, only: BASE_UNIVERSE
use global

2
src/source_header.F90
View file

@ -8,7 +8,7 @@ module source_header
!===============================================================================
type ExtSource
integer :: type ! type of source, e.g. 'box' or 'cell'
integer :: type ! type of source, e.g. 'box' or 'point'
real(8), allocatable :: values(:) ! values for particular source type
end type ExtSource

6
src/tally.F90
View file

@ -1062,7 +1062,8 @@ contains
end subroutine add_to_score
!===============================================================================
! ACCUMULATE_BATCH_ESTIMATE
! ACCUMULATE_BATCH_ESTIMATE accumulates scores from many histories (or many
! generations) into a single realization of a random variable.
!===============================================================================
elemental subroutine accumulate_batch_estimate(score)
@ -1540,7 +1541,8 @@ contains
end function get_label
!===============================================================================
! CALCULATE_STATISTICS
! CALCULATE_STATISTICS determines the sample mean and the standard deviation of
! the mean for a TallyScore.
!===============================================================================
elemental subroutine calculate_statistics(score)

10
src/timing.F90
View file

@ -13,7 +13,7 @@ module timing
type Timer
logical :: running = .false. ! is timer running?
integer :: start_counts = 0 ! counts when started
real(8) :: elapsed = 0. ! total time elapsed in seconds
real(8) :: elapsed = ZERO ! total time elapsed in seconds
!!$ contains
!!$ procedure :: start => timer_start
!!$ procedure :: get_value => timer_get_value
@ -24,7 +24,7 @@ module timing
contains
!===============================================================================
! TIMER_START
! TIMER_START starts running a timer and measures the current time
!===============================================================================
subroutine timer_start(self)
@ -38,7 +38,7 @@ contains
end subroutine timer_start
!===============================================================================
! TIMER_GET_VALUE
! TIMER_GET_VALUE returns the current value of the timer
!===============================================================================
function timer_get_value(self) result(elapsed)
@ -61,7 +61,7 @@ contains
end function timer_get_value
!===============================================================================
! TIMER_STOP
! TIMER_STOP stops the timer and sets the elapsed time
!===============================================================================
subroutine timer_stop(self)
@ -78,7 +78,7 @@ contains
end subroutine timer_stop
!===============================================================================
! TIMER_RESET
! TIMER_RESET resets a timer to have a zero value
!===============================================================================
subroutine timer_reset(self)