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

merged Paul's latest commits and propagated changes from Tally_header file

Browse source
This commit is contained in:
Bryan Herman 2012-03-06 10:01:34 -08:00
commit 925fe5015e
33 changed files with 877 additions and 749 deletions
Download patch file Download diff file Expand all files Collapse all files

7
presentations/crpg_2011-10-06/Makefile
View file

@ -1,7 +0,0 @@
all: slides.tex
mkdir -p build
pdflatex --output-directory=build slides.tex > /dev/null
pdflatex --output-directory=build slides.tex > /dev/null
clean:
rm -rf build

13
presentations/crpg_2011-10-06/annulus.tex
View file

@ -1,13 +0,0 @@
\begin{figure}
\begin{center}
\begin{tikzpicture}
[
transform canvas = {scale=1.0},
]
\fill[mitred] (0,0) circle (1cm);
\fill[white] (0,0) circle (0.50cm);
\end{tikzpicture}
\end{center}
\end{figure}

22
presentations/crpg_2011-10-06/beamercolorthememit.sty
View file

@ -1,22 +0,0 @@
% Copyright 2004 by Madhusudan Singh <madhusudan.singh@gmail.com>
%
% This file may be distributed and/or modified
%
% 1. under the LaTeX Project Public License and/or
% 2. under the GNU Public License.
%
% See the file doc/licenses/LICENSE for more details.
%\DefineNamedColor{named}{mitred} {cmyk}{0,1,0.65,0.34}
\DefineNamedColor{named}{mitred} {rgb}{0.6,0.2,0.2}
\DefineNamedColor{named}{mitgray} {rgb}{0.4,0.4,0.4}
\DefineNamedColor{named}{darkgray} {cmyk}{0,0,0,0.90}
\mode<presentation>
\setbeamercolor{alerted text}{fg=green!80!yellow}
\setbeamercolor*{palette primary}{bg=mitred,fg=white}
\setbeamercolor*{palette secondary}{fg=white,bg=mitgray}
\setbeamercolor*{palette tertiary}{fg=white,bg=darkgray}
\setbeamercolor*{palette quaternary}{fg=white,bg=yellow}
\setbeamercolor*{structure}{fg=mitred,bg=white}
\setbeamercolor{frametitle}{bg=mitred,fg=white}
\mode<all>

47
presentations/crpg_2011-10-06/graphical_settings.tex
View file

@ -1,47 +0,0 @@
\tikzstyle{sNormalBlockStyle} =
[
draw,
rectangle,
rounded corners = 0.1cm,
fill = blue!20,
minimum height = 3em,
minimum width = 6em,
]
\tikzstyle{RectBlue} =
[
draw,
rectangle,
fill = blue!50,
minimum height = 2em,
minimum width = 2em,
]
\tikzstyle{RectGreen} =
[
draw,
rectangle,
fill = green!75,
minimum height = 2em,
minimum width = 2em,
]
\tikzstyle{RectWhite} =
[
draw,
rectangle,
minimum height = 2em,
minimum width = 2em,
]
\tikzstyle{sTextBlockStyle} =
[
draw,
rectangle,
drop shadow,
rounded corners = 0.1cm,
fill = blue!10,
thick,
inner xsep = 0.2cm, % minimum distance between text and borders along x dimension
inner ysep = 0.2cm % minimum distance between text and borders along y dimension
]

BIN
presentations/crpg_2011-10-06/mit-logo.pdf
View file

Binary file not shown.

192
presentations/crpg_2011-10-06/slides.tex
View file

@ -1,192 +0,0 @@
\documentclass{beamer}
%\documentclass[handout,t]{beamer}
\batchmode
% \usepackage{pgfpages}
% \pgfpagesuselayout{4 on 1}[letterpaper,landscape,border shrink=5mm]
\usepackage{amsmath,amssymb,enumerate,epsfig,bbm,calc,color,ifthen,capt-of}
% TikZ packages
\usepackage{tikz}
\usetikzlibrary{arrows,decorations.pathmorphing,decorations.footprints,fadings,calc,trees,mindmap,shadows,decorations.text,patterns,positioning,shapes,matrix,fit}
\usetheme{Berlin}
\usecolortheme{mit}
\title{The OpenMC Monte Carlo Code}
\author{Paul K. Romano}
\date{October 6, 2011}
\pgfdeclareimage[height=0.5cm]{mit-logo}{mit-logo.pdf}
\logo{\pgfuseimage{mit-logo}\hspace*{0.3cm}}
\AtBeginSection[]
{
\begin{frame}<beamer>
\frametitle{Outline}
\tableofcontents[currentsection]
\end{frame}
}
\beamerdefaultoverlayspecification{<+->}
\include{graphical_settings}
% -----------------------------------------------------------------------------
\begin{document}
% -----------------------------------------------------------------------------
\frame{\titlepage}
\section[Outline]{}
\begin{frame}{Outline}
\tableofcontents
\end{frame}
% -----------------------------------------------------------------------------
\section{Introduction}
\begin{frame}{Background}
My background:
\begin{itemize}
\item<1-> B.S. Nuclear Engineering, RPI (2007)
\item<1-> M.S. Nuclear Scince and Engineering, MIT (2009)
\end{itemize}
Origin of OpenMC:
\begin{itemize}
\item<1-> Working on advanced parallelization for Monte Carlo
\item<1-> Needed testing platform - MC21?
\item<1-> Ultimately decided to start from scratch
\end{itemize}
\end{frame}
\begin{frame}{Goals}
Overall objectives of OpenMC:
\begin{itemize}
\item<1-> Fully featured, capable of realistic physics
\item<1-> Written in a modern programming language (F2003)
\item<1-> Easy to understand inner workings
\item<1-> High performance
\item<1-> Extensible for research purposes
\item<1-> Open source and freely available
\end{itemize}
\end{frame}
\begin{frame}{Monte Carlo}
Rather than solving the transport equation by discretizing the spatial,
energy, and angular variables, Monte Carlo solves the transport equation by
following individual particles as they move stochastically through a medium.
\begin{itemize}
\item<1-> Standard deviation $\propto \frac{1}{\sqrt{N}}$
\end{itemize}
\end{frame}
% -----------------------------------------------------------------------------
\section{Methods and Theory}
\begin{frame}{Solution Algorithm}
\begin{itemize}
\item<1-> Loop over particles
\begin{enumerate}
\item<1-> Sample particle from distribution
\item<1-> Track particle to next collision
\item<1-> Sample nuclide within material
\item<1-> Sample reaction within nuclide
\item<1-> Repeat from step 2 until particle dies
\end{enumerate}
\end{itemize}
While particle is being tracked, keep track of collision rate, fission rate,
etc.
\begin{equation*}
\phi = \frac{1}{NV} \sum\limits_{\text{all collisions} \atop \text{in cell}} \frac{w_i}{\Sigma_t (E_i)}
\end{equation*}
\end{frame}
\begin{frame}{Geometry}
\begin{itemize}
\item<1-> All geometry is described as constructive solid geometry (unions and
intersections of second-order surfaces)
\item<1-> $Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Iz + J = 0$
\item<1-> e.g. to construct an annulus
\end{itemize}
\include{annulus}
\end{frame}
\begin{frame}{Cross-Sections}
\begin{itemize}
\item<1-> No need to reinvent the wheel
\item<1-> ACE (A Compact ENDF) Format
\begin{itemize}
\item<1-> MCNP
\item<1-> SERPENT
\end{itemize}
\item<1-> Three arrays: NXS, JXS, XSS
\item<1-> Parse arrays into internal derived types
\end{itemize}
\end{frame}
\begin{frame}{Union Energy Grid}
\begin{itemize}
\item<1-> Each nuclide has cross-sections tabulated at different energy points
\item<1-> To determine collision type, need to look up microscopic reaction
cross sections
\end{itemize}
\include{union_energy_grid}
\end{frame}
\begin{frame}{Tallies}{}
{\bf Remember:} In Monte Carlo, the only answer you get is the one you ask
for. In general, you don't get the global solution like you do in
deterministic land.
\begin{itemize}
\item<1-> Implemented a robust tally system to calculate user-specified
quantities of interest
\item<1-> {\bf Filters:} Spatial location, incoming/outgoing energy, birth region, mesh
\item<1-> {\bf Responses:} Flux, reaction rates, currents, etc.
\end{itemize}
\end{frame}
\begin{frame}{Parallel Fission Bank}
\end{frame}
% -----------------------------------------------------------------------------
\section{Using OpenMC}
\begin{frame}{Compiling}
\end{frame}
\begin{frame}{XML Input Format}
Unlike many nuclear codes, OpenMC uses a modular XML input format
\begin{itemize}
\item<1-> Building a geometric model -- {\bf geometry.xml}
\item<1-> Assigning materials to volumes -- {\bf materials.xml}
\item<1-> Setting parameters for the simulation -- {\bf settings.xml}
\item<1-> Determining which quantities to score -- {\bf tallies.xml}
\end{itemize}
\end{frame}
\begin{frame}{Running}
\end{frame}
\begin{frame}{Post-run Analysis}
\end{frame}
% -----------------------------------------------------------------------------
\section{Development}
\begin{frame}{Version Control}
\end{frame}
% -----------------------------------------------------------------------------
\section{Conclusions}
\begin{frame}{Questions?}
References:
\begin{itemize}
\item<1-> Github: \url{http://github.com/paulromano/openmc}.
\item<1-> Documentation: \url{http://paulromano.github.com/openmc/}.
\end{itemize}
\end{frame}
% -----------------------------------------------------------------------------
\end{document}

51
presentations/crpg_2011-10-06/union_energy_grid.tex
View file

@ -1,51 +0,0 @@
\begin{figure}
\begin{center}
\begin{tikzpicture}
[
transform canvas = {scale=0.7},
xshift = -5.0cm,
yshift = 1.0cm
]
\node (LabelOne) [sTextBlockStyle] {Union Energy Grid};
\node (LabelTwo) [sTextBlockStyle, below of = LabelOne, node distance = 1.5cm] {Nuclide Pointers};
\node (LabelThree) [sTextBlockStyle, below of = LabelTwo, node distance = 1.5cm] {Nuclide Energy Grid};
\node (EOne) [RectBlue, right of = LabelOne, node distance = 3cm] {$E_1$};
\node (ETwo) [RectBlue, right of = EOne] {$E_2$};
\node (EThree) [RectBlue, right of = ETwo] {$E_3$};
\node (EFour) [RectBlue, right of = EThree] {$E_4$};
\node (EFive) [RectBlue, right of = EFour] {$E_5$};
\node (ESix) [RectBlue, right of = EFive] {$E_6$};
\node (ESeven) [RectBlue, right of = ESix] {$E_7$};
\node (EEight) [RectBlue, right of = ESeven] {$E_8$};
\node (POne) [RectWhite, below of = EOne, node distance = 1.5cm] {0};
\node (PTwo) [RectWhite, below of = ETwo, node distance = 1.5cm] {0};
\node (PThree) [RectWhite, below of = EThree, node distance = 1.5cm] {1};
\node (PFour) [RectWhite, below of = EFour, node distance = 1.5cm] {1};
\node (PFive) [RectWhite, below of = EFive, node distance = 1.5cm] {1};
\node (PSix) [RectWhite, below of = ESix, node distance = 1.5cm] {2};
\node (PSeven) [RectWhite, below of = ESeven, node distance = 1.5cm] {3};
\node (PEight) [RectWhite, below of = EEight, node distance = 1.5cm] {3};
\node (ENucOne) [RectGreen, below of = PThree, node distance = 1.5cm] {$E_1$};
\node (ENucTwo) [RectGreen, below of = PSix, node distance = 1.5cm] {$E_2$};
\node (ENucThree) [RectGreen, below of = PSeven, node distance = 1.5cm] {$E_3$};
\draw [-triangle 45, thick] (EOne) -- (POne);
\draw [-triangle 45, thick] (ETwo) -- (PTwo);
\draw [-triangle 45, thick] (EThree) -- (PThree);
\draw [-triangle 45, thick] (EFour) -- (PFour);
\draw [-triangle 45, thick] (EFive) -- (PFive);
\draw [-triangle 45, thick] (ESix) -- (PSix);
\draw [-triangle 45, thick] (ESeven) -- (PSeven);
\draw [-triangle 45, thick] (EEight) -- (PEight);
\draw [dashed] (-2,-0.70) -- (10.5,-0.70);
\draw [dashed] (-2,-2.20) -- (10.5,-2.20);
\end{tikzpicture}
\end{center}
\end{figure}

11
src/DEPENDENCIES
View file

@ -114,11 +114,14 @@ global.o: geometry_header.o
global.o: material_header.o
global.o: mesh_header.o
global.o: particle_header.o
global.o: plot_header.o
global.o: source_header.o
global.o: tally_header.o
global.o: timing.o
hdf5_interface.o: ace_header.o
hdf5_interface.o: constants.o
hdf5_interface.o: endf.o
hdf5_interface.o: geometry_header.o
hdf5_interface.o: global.o
hdf5_interface.o: material_header.o
@ -151,6 +154,7 @@ input_xml.o: geometry_header.o
input_xml.o: global.o
input_xml.o: mesh_header.o
input_xml.o: output.o
input_xml.o: plot_header.o
input_xml.o: string.o
input_xml.o: tally_header.o
input_xml.o: xml-fortran/templates/cross_sections_t.o
@ -167,6 +171,7 @@ intercycle.o: output.o
intercycle.o: random_lcg.o
intercycle.o: search.o
intercycle.o: string.o
intercycle.o: tally.o
intercycle.o: tally_header.o
intercycle.o: timing.o
@ -204,6 +209,7 @@ output.o: geometry_header.o
output.o: global.o
output.o: mesh_header.o
output.o: particle_header.o
output.o: plot_header.o
output.o: string.o
output.o: tally_header.o
@ -231,8 +237,12 @@ plot.o: error.o
plot.o: geometry.o
plot.o: geometry_header.o
plot.o: global.o
plot.o: output.o
plot.o: particle_header.o
plot.o: plot_header.o
plot.o: ppmlib.o
plot.o: source.o
plot.o: string.o
random_lcg.o: global.o
@ -258,7 +268,6 @@ string.o: global.o
tally.o: constants.o
tally.o: error.o
tally.o: global.o
tally.o: intercycle.o
tally.o: mesh.o
tally.o: mesh_header.o
tally.o: output.o

3
src/Makefile
View file

@ -16,7 +16,6 @@ include OBJECTS
#===============================================================================
F90X = /home/bherman/petsc/release-ifort/bin/mpif90
F90XFLAGS = -Wall -Wno-unused-variable
COMPILER = petsc
DEBUG = no
PROFILE = no
@ -30,7 +29,7 @@ USE_HDF5 = yes
ifeq ($(COMPILER),gnu)
F90 = gfortran
F90FLAGS := -cpp -fbacktrace -DN2008
F90FLAGS := -cpp -fbacktrace -DNO_F2008
LDFLAGS =
# Debugging

2
src/OBJECTS
View file

@ -39,6 +39,8 @@ output.o \
particle_header.o \
physics.o \
plot.o \
plot_header.o \
ppmlib.o \
random_lcg.o \
search.o \
source.o \

6
src/ace.F90
View file

@ -137,7 +137,7 @@ contains
! READ ALL ACE CROSS SECTION TABLES
! display header in summary.out
if (master) call header("CROSS SECTION TABLES", unit=UNIT_SUMMARY)
if (master) call header("CROSS SECTION TABLES", unit=UNIT_XS)
call dict_create(already_read)
@ -156,9 +156,9 @@ contains
call read_ace_table(index_nuclides, index_list)
! Print out information on table to summary.out file
! Print out information on table to cross_sections.out file
nuc => nuclides(index_nuclides)
if (master) call print_nuclide(nuc, unit=UNIT_SUMMARY)
if (master) call print_nuclide(nuc, unit=UNIT_XS)
call dict_add_key(already_read, name, 0)
call dict_add_key(already_read, alias, 0)

36
src/cmfd_data.F90
View file

@ -136,7 +136,7 @@ contains
score_index = sum((bins - 1) * t%stride) + 1
! get flux
flux = t % scores(score_index,1) % val
flux = t % scores(score_index,1) % sum
cmfd % flux(h,i,j,k) = flux
! detect zero flux
@ -146,10 +146,10 @@ contains
end if
! get total rr and convert to total xs
cmfd % totalxs(h,i,j,k) = t % scores(score_index,2) % val / flux
cmfd % totalxs(h,i,j,k) = t % scores(score_index,2) % sum / flux
! get p1 scatter rr and convert to p1 scatter xs
cmfd % p1scattxs(h,i,j,k) = t % scores(score_index,3) % val / flux
cmfd % p1scattxs(h,i,j,k) = t % scores(score_index,3) % sum / flux
! calculate diffusion coefficient
cmfd % diffcof(h,i,j,k) = 1.0_8/(3.0_8*(cmfd % totalxs(h,i,j,k) -&
@ -180,16 +180,16 @@ contains
score_index = sum((bins - 1) * t%stride) + 1
! get scattering
cmfd % scattxs(h,g,i,j,k) = t % scores(score_index,1) % val /&
cmfd % scattxs(h,g,i,j,k) = t % scores(score_index,1) % sum /&
& cmfd % flux(h,i,j,k)
! get nu-fission
cmfd % nfissxs(h,g,i,j,k) = t % scores(score_index,2) % val /&
cmfd % nfissxs(h,g,i,j,k) = t % scores(score_index,2) % sum /&
& cmfd % flux(h,i,j,k)
! bank source
cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + &
& t % scores(score_index,2) % val
& t % scores(score_index,2) % sum
end do INGROUP
@ -203,55 +203,55 @@ contains
bins(1:3) = (/ i-1, j, k /) + 1
bins(SURF_FILTER_SURFACE) = IN_RIGHT
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(1,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(1,h,i,j,k) = t % scores(score_index,1) % sum
bins(SURF_FILTER_SURFACE) = OUT_RIGHT
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
cmfd % current(2,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(2,h,i,j,k) = t % scores(score_index,1) % sum
! right surface
bins(1:3) = (/ i, j, k /) + 1
bins(SURF_FILTER_SURFACE) = IN_RIGHT
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
cmfd % current(3,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(3,h,i,j,k) = t % scores(score_index,1) % sum
bins(SURF_FILTER_SURFACE) = OUT_RIGHT
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(4,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(4,h,i,j,k) = t % scores(score_index,1) % sum
! back surface
bins(1:3) = (/ i, j-1, k /) + 1
bins(SURF_FILTER_SURFACE) = IN_FRONT
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(5,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(5,h,i,j,k) = t % scores(score_index,1) % sum
bins(SURF_FILTER_SURFACE) = OUT_FRONT
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
cmfd % current(6,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(6,h,i,j,k) = t % scores(score_index,1) % sum
! front surface
bins(1:3) = (/ i, j, k /) + 1
bins(SURF_FILTER_SURFACE) = IN_FRONT
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
cmfd % current(7,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(7,h,i,j,k) = t % scores(score_index,1) % sum
bins(SURF_FILTER_SURFACE) = OUT_FRONT
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(8,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(8,h,i,j,k) = t % scores(score_index,1) % sum
! bottom surface
bins(1:3) = (/ i, j, k-1 /) + 1
bins(SURF_FILTER_SURFACE) = IN_TOP
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(9,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(9,h,i,j,k) = t % scores(score_index,1) % sum
bins(SURF_FILTER_SURFACE) = OUT_TOP
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
cmfd % current(10,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(10,h,i,j,k) = t % scores(score_index,1) % sum
! top surface
bins(1:3) = (/ i, j, k /) + 1
bins(SURF_FILTER_SURFACE) = IN_TOP
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
cmfd % current(11,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(11,h,i,j,k) = t % scores(score_index,1) % sum
bins(SURF_FILTER_SURFACE) = OUT_TOP
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(12,h,i,j,k) = t % scores(score_index,1) % val
cmfd % current(12,h,i,j,k) = t % scores(score_index,1) % sum
end if TALLY

12
src/cmfd_output.F90
View file

@ -18,7 +18,7 @@ contains
subroutine neutron_balance()
use constants, only: ONE
use global, only: cmfd,UNIT_CMFD,keff,MAX_FILE_LEN,current_cycle
use global, only: cmfd,keff,MAX_FILE_LEN,current_cycle
use string
integer :: i ! iteration counter for x
@ -95,13 +95,6 @@ contains
! bank res in cmfd object
cmfd%resnb(g,i,j,k) = res
! write output
! label = "MESH (" // trim(int4_to_str(i)) // ". " // &
! & trim(int4_to_str(j)) // ", " // trim(int4_to_str(k)) // &
! & ") GROUP " // trim(int4_to_str(g))
! write(UNIT=UNIT_CMFD, FMT='(A,T35,A)') label, &
! & trim(real_to_str(res))
end do GROUPG
end do XLOOP
@ -110,9 +103,6 @@ contains
end do ZLOOP
! close file
close(UNIT=UNIT_CMFD)
end subroutine neutron_balance
!===============================================================================

10
src/constants.F90
View file

@ -293,6 +293,14 @@ module constants
IN_TOP = 5, &
OUT_TOP = 6
! Global tallY parameters
integer, parameter :: N_GLOBAL_TALLIES = 4
integer, parameter :: &
K_ANALOG = 1, &
K_COLLISION = 2, &
K_TRACKLENGTH = 3, &
LEAKAGE = 4
! ============================================================================
! MISCELLANEOUS CONSTANTS
@ -318,6 +326,6 @@ module constants
integer, parameter :: UNIT_SUMMARY = 11 ! unit # for writing summary file
integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file
integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file
integer, parameter :: UNIT_CMFD = 14 ! unit # for writing cmfd file
integer, parameter :: UNIT_XS = 14 ! unit # for writing xs summary file
end module constants

4
src/finalize.F90
View file

@ -6,7 +6,7 @@ module finalize
use timing, only: timer_start, timer_stop
#ifdef HDF5
use hdf5_interface, only: hdf5_write_timing, hdf5_close_output
use hdf5_interface, only: hdf5_write_results, hdf5_close_output
#endif
implicit none
@ -37,7 +37,7 @@ contains
#ifdef HDF5
! Write time statistics to HDF5 output
if (master) then
call hdf5_write_timing()
call hdf5_write_results()
call hdf5_close_output()
end if
#endif

17
src/geometry.F90
View file

@ -8,7 +8,7 @@ module geometry
use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord
use string, only: to_str
use tally, only: score_surface_current
use tally, only: score_surface_current, add_to_score
implicit none
@ -263,12 +263,17 @@ contains
! forward slightly so that if the mesh boundary is on the surface, it is
! still processed
if (tallies_on .and. n_current_tallies > 0) then
! TODO: Find a better solution to score surface currents than
! physically moving the particle forward slightly
if (tallies_on) then
if (n_current_tallies > 0) then
! TODO: Find a better solution to score surface currents than
! physically moving the particle forward slightly
p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
call score_surface_current()
p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
call score_surface_current()
end if
! Score to global leakage tally
call add_to_score(global_tallies(LEAKAGE), p % wgt)
end if
! Display message

1
src/geometry_header.F90
View file

@ -64,6 +64,7 @@ module geometry_header
& surfaces(:) ! List of surfaces bounding cell -- note that
! parentheses, union, etc operators will be listed
! here too
integer, allocatable :: rgb(:) ! plotting color
end type Cell
! array index of universe 0

18
src/global.F90
View file

@ -10,8 +10,9 @@ module global
use material_header, only: Material
use mesh_header, only: StructuredMesh
use particle_header, only: Particle
use plot_header, only: Plot
use source_header, only: ExtSource
use tally_header, only: TallyObject, TallyMap
use tally_header, only: TallyObject, TallyMap, TallyScore
use timing, only: Timer
#ifdef MPI
@ -39,6 +40,7 @@ module global
type(Lattice), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(Plot), allocatable, target :: plots(:)
! Size of main arrays
integer :: n_cells ! # of cells
@ -46,6 +48,7 @@ module global
integer :: n_lattices ! # of lattices
integer :: n_surfaces ! # of surfaces
integer :: n_materials ! # of materials
integer :: n_plots ! # of plots
! These dictionaries provide a fast lookup mechanism -- the key is the
! user-specified identifier and the value is the index in the corresponding
@ -97,6 +100,14 @@ module global
integer, allocatable :: tracklength_tallies(:)
integer, allocatable :: current_tallies(:)
! Global tallies
! 1) analog estimate of k-eff
! 2) collision estimate of k-eff
! 3) track-length estimate of k-eff
! 4) leakage fraction
type(TallyScore) :: global_tallies(N_GLOBAL_TALLIES)
! Tally map structure
type(TallyMap), allocatable :: tally_maps(:)
@ -183,10 +194,6 @@ module global
! PLOTTING VARIABLES
logical :: plotting = .false.
real(8) :: plot_origin(3)
real(8) :: plot_width(2)
real(8) :: plot_basis(6)
real(8) :: pixel
! ============================================================================
! HDF5 VARIABLES
@ -253,6 +260,7 @@ contains
if (allocated(lattices)) deallocate(lattices)
if (allocated(surfaces)) deallocate(surfaces)
if (allocated(materials)) deallocate(materials)
if (allocated(plots)) deallocate(plots)
! Deallocate cross section data, listings, and cache
if (allocated(nuclides)) deallocate(nuclides)

153
src/hdf5_interface.F90
View file

@ -1,6 +1,8 @@
module hdf5_interface
use ace_header, only: Reaction, UrrData
use constants
use endf, only: reaction_name
use geometry_header, only: Cell, Surface, Universe, Lattice
use global
use material_header, only: Material
@ -116,12 +118,24 @@ contains
call hdf5_write_geometry()
call hdf5_write_materials()
call hdf5_write_nuclides()
if (n_tallies > 0) then
call hdf5_write_tallies()
end if
end subroutine hdf5_write_summary
!===============================================================================
! HDF5_WRITE_RESULTS
!===============================================================================
subroutine hdf5_write_results()
call hdf5_write_timing()
call hdf5_write_global_tallies()
end subroutine hdf5_write_results
!===============================================================================
! HDF5_WRITE_GEOMETRY
!===============================================================================
@ -431,6 +445,30 @@ contains
end subroutine hdf5_write_materials
!===============================================================================
! HDF5_WRITE_GLOBAL_TALLIES
!===============================================================================
subroutine hdf5_write_global_tallies()
integer :: rank = 1
integer(HSIZE_T) :: dims(1) = (/ 2 /)
call h5ltmake_dataset_double_f(hdf5_output_file, "k_analog", &
rank, dims, (/ global_tallies(K_ANALOG) % sum, &
global_tallies(K_ANALOG) % sum_sq /), hdf5_err)
call h5ltmake_dataset_double_f(hdf5_output_file, "k_collision", &
rank, dims, (/ global_tallies(K_COLLISION) % sum, &
global_tallies(K_COLLISION) % sum_sq /), hdf5_err)
call h5ltmake_dataset_double_f(hdf5_output_file, "k_tracklength", &
rank, dims, (/ global_tallies(K_TRACKLENGTH) % sum, &
global_tallies(K_TRACKLENGTH) % sum_sq /), hdf5_err)
call h5ltmake_dataset_double_f(hdf5_output_file, "leakage", &
rank, dims, (/ global_tallies(LEAKAGE) % sum, &
global_tallies(LEAKAGE) % sum_sq /), hdf5_err)
end subroutine hdf5_write_global_tallies
!===============================================================================
! HDF5_WRITE_TALLIES
!===============================================================================
@ -622,6 +660,121 @@ contains
end subroutine hdf5_write_tallies
!===============================================================================
! HDF5_WRITE_NUCLIDES
!===============================================================================
subroutine hdf5_write_nuclides()
integer :: i, j
integer :: size_total
integer :: size_xs
integer :: size_angle
integer :: size_energy
integer(HID_T) :: group
integer(HID_T) :: nuclide_group
integer(HID_T) :: reactions_group
integer(HID_T) :: rxn_group
type(Nuclide), pointer :: nuc => null()
type(Reaction), pointer :: rxn => null()
type(UrrData), pointer :: urr => null()
! Create group for nuclides
call h5gcreate_f(hdf5_output_file, "/nuclides", group, hdf5_err)
! Use H5LT interface to write number of nuclides
call hdf5_make_integer(group, "n_nuclides", n_nuclides_total)
! Write information on each nuclide
do i = 1, n_nuclides_total
nuc => nuclides(i)
! Determine size of cross-sections
size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8
size_total = size_xs
! Create group for i-th nuclide
call h5gcreate_f(group, trim(nuc % name), nuclide_group, hdf5_err)
! Write some basic attributes
call hdf5_make_integer(nuclide_group, "zaid", nuc % zaid)
call hdf5_make_double(nuclide_group, "awr", nuc % awr)
call hdf5_make_double(nuclide_group, "kT", nuc % kT)
call hdf5_make_integer(nuclide_group, "n_grid", nuc % n_grid)
call hdf5_make_integer(nuclide_group, "n_reactions", nuc % n_reaction)
call hdf5_make_integer(nuclide_group, "n_fission", nuc % n_fission)
call hdf5_make_integer(nuclide_group, "size_xs", size_xs)
! =======================================================================
! WRITE INFORMATION ON EACH REACTION
! Create overall group
call h5gcreate_f(nuclide_group, "reactions", reactions_group, hdf5_err)
do j = 1, nuc % n_reaction
! Information on each reaction
rxn => nuc % reactions(j)
! Determine size of angle distribution
if (rxn % has_angle_dist) then
size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8
else
size_angle = 0
end if
! Determine size of energy distribution
if (rxn % has_energy_dist) then
size_energy = size(rxn % edist % data) * 8
else
size_energy = 0
end if
! Create reaction group
call h5gcreate_f(reactions_group, reaction_name(rxn % MT), &
rxn_group, hdf5_err)
! Write information on reaction
call hdf5_make_double(rxn_group, "Q_value", rxn % Q_value)
call hdf5_make_integer(rxn_group, "n_neutrons", rxn % TY)
call hdf5_make_integer(rxn_group, "grid_index", rxn % IE)
call hdf5_make_integer(rxn_group, "size_angle", size_angle)
call hdf5_make_integer(rxn_group, "size_energy", size_energy)
call h5gclose_f(rxn_group, hdf5_err)
! Accumulate data size
size_total = size_total + size_angle + size_energy
end do
! Close overall group for reactions
call h5gclose_f(reactions_group, hdf5_err)
! =======================================================================
! WRITE INFORMATION ON URR PROBABILITY TABLES
if (nuc % urr_present) then
urr => nuc % urr_data
call hdf5_make_integer(nuclide_group, "urr_n_energy", urr % n_energy)
call hdf5_make_integer(nuclide_group, "urr_n_prob", urr % n_prob)
call hdf5_make_integer(nuclide_group, "urr_interp", urr % interp)
call hdf5_make_integer(nuclide_group, "urr_inelastic", urr % inelastic_flag)
call hdf5_make_integer(nuclide_group, "urr_absorption", urr % absorption_flag)
call hdf5_make_double(nuclide_group, "urr_min_E", urr % energy(1))
call hdf5_make_double(nuclide_group, "urr_max_E", urr % energy(urr % n_energy))
end if
! Write total memory used
call hdf5_make_integer(nuclide_group, "size_total", size_total)
! Close group for i-th nuclide
call h5gclose_f(nuclide_group, hdf5_err)
end do
! Close group for nuclides
call h5gclose_f(group, hdf5_err)
end subroutine hdf5_write_nuclides
!===============================================================================
! HDF5_WRITE_TIMING
!===============================================================================

6
src/initialize.F90
View file

@ -14,7 +14,8 @@ module initialize
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict
use output, only: title, header, print_summary, print_geometry, &
print_plot, create_summary_file
print_plot, create_summary_file, &
create_xs_summary_file
use random_lcg, only: initialize_prng
use source, only: initialize_source
use string, only: to_str, starts_with, ends_with, lower_case
@ -56,8 +57,9 @@ contains
call read_command_line()
if (master) then
! Create summary.out file
! Create output files
call create_summary_file()
call create_xs_summary_file()
#ifdef HDF5
! Open HDF5 output file for writing and write header information

137
src/input_xml.F90
View file

@ -9,6 +9,8 @@ module input_xml
use global
use mesh_header, only: StructuredMesh
use output, only: write_message
use plot_header
use random_lcg, only: prn
use string, only: lower_case, to_str, str_to_int, str_to_real, &
split_string, starts_with, ends_with
use tally_header, only: TallyObject
@ -32,7 +34,7 @@ contains
call read_materials_xml()
call read_tallies_xml()
call read_cmfd_xml()
if (plotting) call read_plot_xml()
if (plotting) call read_plots_xml()
end subroutine read_input_xml
@ -264,6 +266,14 @@ contains
c % material = cell_(i) % material
c % fill = cell_(i) % fill
! Set plot color
if (plotting) then
allocate(c % rgb(3))
c % rgb(1) = prn()*255
c % rgb(2) = prn()*255
c % rgb(3) = prn()*255
end if
! Check to make sure that either material or fill was specified
if (c % material == 0 .and. c % fill == 0) then
message = "Neither material nor fill was specified for cell " // &
@ -534,6 +544,14 @@ contains
! Copy material id
m % id = material_(i) % id
! Set plot color
if (plotting) then
allocate(m % rgb(3))
m % rgb(1) = prn()*255
m % rgb(2) = prn()*255
m % rgb(3) = prn()*255
end if
! 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.
@ -1091,22 +1109,23 @@ contains
end subroutine read_tallies_xml
!===============================================================================
! READ_TALLIES_XML reads data from a tallies.xml file and parses it, checking
! for errors and placing properly-formatted data in the right data structures
! READ_PLOTS_XML reads data from a plots.xml file
!===============================================================================
subroutine read_plot_xml
subroutine read_plots_xml
use xml_data_plot_t
use xml_data_plots_t
logical :: file_exists ! does tallies.xml file exist?
character(MAX_LINE_LEN) :: filename ! absolute path to tallies.xml
integer i
logical :: file_exists ! does plots.xml file exist?
character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml
type(Plot), pointer :: pl => null()
! Check if plot.xml exists
filename = trim(path_input) // "plot.xml"
! Check if plots.xml exists
filename = trim(path_input) // "plots.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
message = "Plot XML file '" // trim(filename) // "' does not exist!"
message = "Plots XML file '" // trim(filename) // "' does not exist!"
call fatal_error()
end if
@ -1114,33 +1133,87 @@ contains
message = "Reading plot XML file..."
call write_message(5)
! Parse plot.xml file
call read_xml_file_plot_t(filename)
! Parse plots.xml file
call read_xml_file_plots_t(filename)
! Copy plotting origin
if (size(origin_) == 3) then
plot_origin = origin_
end if
! Allocate plots array
n_plots = size(plot_)
allocate(plots(n_plots))
! Copy plotting width
if (size(width_) == 2) then
plot_width = width_
end if
do i = 1, n_plots
pl => plots(i)
! Read basis
select case (basis_)
case ("xy")
plot_basis = (/ 1, 0, 0, 0, 1, 0 /)
case ("yz")
plot_basis = (/ 0, 1, 0, 0, 0, 1 /)
case ("xz")
plot_basis = (/ 1, 0, 0, 0, 0, 1 /)
end select
! Copy data into plots
pl % id = plot_(i) % id
pl % aspect = plot_(i) % aspect
! Read pixel width
pixel = pixel_
if (size(plot_(i) % pixels) == 2) then
pl % pixels = plot_(i) % pixels
else
message = "<pixels> must be length 2 in plot " // to_str(i)
call fatal_error()
end if
end subroutine read_plot_xml
select case (plot_(i) % color)
case ("cell")
pl % color = PLOT_COLOR_CELLS
case ("mat")
pl % color = PLOT_COLOR_MATS
case default
message = "Unsupported plot color '" // plot_(i) % color &
// "' in plot " // trim(to_str(i))
call fatal_error()
end select
select case (plot_(i) % type)
case ("slice")
pl % type = PLOT_TYPE_SLICE
!case ("points")
! pl % type = PLOT_TYPE_POINTS
case default
message = "Unsupported plot type '" // plot_(i) % type &
// "' in plot " // trim(to_str(i))
call fatal_error()
end select
select case (plot_(i) % basis)
case ("xy")
pl % basis = PLOT_BASIS_XY
case ("xz")
pl % basis = PLOT_BASIS_XZ
case ("yz")
pl % basis = PLOT_BASIS_YZ
case default
message = "Unsupported plot basis '" // plot_(i) % basis &
// "' in plot " // trim(to_str(i))
call fatal_error()
end select
! Copy plotting origin
if (size(plot_(i) % origin) == 3) then
pl % origin = plot_(i) % origin
else
message = "Origin must be length 3 " &
// "in plot " // trim(to_str(i))
call fatal_error()
end if
! Copy plotting width
if (size(plot_(i) % width) == 3) then
pl % width = plot_(i) % width
else if (size(plot_(i) % width) == 2) then
pl % width(1) = plot_(i) % width(1)
pl % width(2) = plot_(i) % width(2)
else
message = "Bad plot width " &
// "in plot " // trim(to_str(i))
call fatal_error()
end if
end do
end subroutine read_plots_xml
!===============================================================================
! READ_CROSS_SECTIONS_XML reads information from a cross_sections.xml file. This

107
src/intercycle.F90
View file

@ -10,6 +10,7 @@ module intercycle
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use string, only: to_str
use tally, only: accumulate_cycle_estimate
use tally_header, only: TallyObject
use timing, only: timer_start, timer_stop
@ -273,55 +274,6 @@ contains
end subroutine synchronize_bank
!===============================================================================
! REDUCE_TALLIES collects all the results from tallies onto one processor
!===============================================================================
#ifdef MPI
subroutine reduce_tallies()
integer :: i ! loop index for tallies
integer :: n ! number of filter bins
integer :: m ! number of score bins
integer :: n_bins ! total number of bins
real(8), allocatable :: tally_temp(:,:) ! contiguous array of scores
type(TallyObject), pointer :: t => null()
do i = 1, n_tallies
t => tallies(i)
n = t % n_total_bins
m = t % n_score_bins
n_bins = n*m
allocate(tally_temp(n,m))
tally_temp = t % scores(:,:) % val_history
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into a
! receive buffer without having a temporary variable
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
! Transfer values to val_history on master
t % scores(:,:) % val_history = tally_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(tally_temp, tally_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
! Reset val_history on other processors
t % scores(:,:) % val_history = 0
end if
deallocate(tally_temp)
end do
end subroutine reduce_tallies
#endif
!===============================================================================
! SHANNON_ENTROPY calculates the Shannon entropy of the fission source
! distribution to assess source convergence
@ -433,48 +385,55 @@ contains
subroutine calculate_keff()
integer(8) :: total_bank ! total number of source sites
integer :: n ! active cycle number
real(8) :: k_cycle ! single cycle estimate of keff
real(8), save :: k_sum ! accumulated keff
real(8), save :: k_sum_sq ! accumulated keff**2
integer :: n ! active cycle number
real(8) :: k_cycle ! single cycle estimate of keff
#ifdef MPI
real(8) :: global_temp(N_GLOBAL_TALLIES)
#endif
message = "Calculate cycle keff..."
call write_message(8)
! initialize sum and square of sum at beginning of run
if (current_cycle == 1) then
k_sum = ZERO
k_sum_sq = ZERO
end if
! Since the creation of bank sites was originally weighted by the last
! cycle keff, we need to multiply by that keff to get the current cycle's
! value
#ifdef MPI
! Collect number bank sites onto master process
call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
global_tallies(K_ANALOG) % value = n_bank * keff
! Copy global tallies into array to be reduced
global_temp = global_tallies(:) % value
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Transfer values back to global_tallies on master
global_tallies(:) % value = global_temp
else
call MPI_REDUCE(global_temp, global_temp, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Reset value on other processors
global_tallies(:) % value = ZERO
end if
#else
total_bank = n_bank
global_tallies(K_ANALOG) % value = n_bank * keff
#endif
! Collect statistics and print output
if (master) then
! Since the creation of bank sites was originally weighted by the last
! cycle keff, we need to multiply by that keff to get the current cycle's
! value
k_cycle = real(total_bank)/real(n_particles)*keff
k_cycle = global_tallies(K_ANALOG) % value/n_particles
if (current_cycle > n_inactive) then
! Active cycle number
n = current_cycle - n_inactive
! Accumulate cycle estimate of k
k_sum = k_sum + k_cycle
k_sum_sq = k_sum_sq + k_cycle*k_cycle
! Accumulate single cycle realizations of k
call accumulate_cycle_estimate(global_tallies)
! Determine mean and standard deviation of mean
keff = k_sum/n
keff_std = sqrt((k_sum_sq/n - keff*keff)/n)
keff = global_tallies(K_ANALOG) % sum/n
keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - keff*keff)/n)
! Display output for this cycle
if (current_cycle > n_inactive + 1) then

2
src/main.F90
View file

@ -7,7 +7,7 @@ program main
use initialize, only: initialize_run
use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank
use output, only: write_message, header
use plot, only: run_plot
use plotter, only: run_plot
use physics, only: transport
use random_lcg, only: set_particle_seed
use source, only: get_source_particle

2
src/material_header.F90
View file

@ -21,6 +21,8 @@ module material_header
character(12) :: sab_name ! name of S(a,b) table
integer :: sab_table = 0 ! index in sab_tables
integer :: sab_nuclide = 0 ! index of nuclide which has S(a,b) table
integer, allocatable :: rgb(:) ! plotting color
end type Material
end module material_header

94
src/output.F90
View file

@ -10,6 +10,7 @@ module output
use global
use mesh_header, only: StructuredMesh
use particle_header, only: LocalCoord
use plot_header
use string, only: upper_case, to_str
use tally_header, only: TallyObject
@ -891,21 +892,46 @@ contains
subroutine print_plot()
integer i
type(Plot), pointer :: pl => null()
! Display header for plotting
call header("PLOTTING SUMMARY")
! Print plotting origin
write(ou,100) "Plotting Origin:", trim(to_str(plot_origin(1))) // &
" " // trim(to_str(plot_origin(2))) // " " // &
trim(to_str(plot_origin(3)))
do i=1,n_plots
pl => plots(i)
! Print plotting width
write(ou,100) "Plotting Width:", trim(to_str(plot_width(1))) // &
" " // trim(to_str(plot_width(2)))
! Print plot id
write(ou,100) "Plot ID:", trim(to_str(pl % id))
! Print pixel width
write(ou,100) "Pixel Width:", trim(to_str(pixel))
write(ou,*)
! Print 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
if (pl % type == PLOT_TYPE_SLICE) then
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
" " // trim(to_str(pl % width(2)))
write(ou,100) "Coloring:", trim(to_str(pl % color))
write(ou,100) "Basis:", trim(to_str(pl % basis))
write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // &
trim(to_str(pl % pixels(2)))
else if (pl % type == PLOT_TYPE_POINTS) then
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
" " // trim(to_str(pl % width(2))) // " " &
// trim(to_str(pl % width(3)))
write(ou,100) "Coloring:", trim(to_str(pl % color))
write(ou,100) "Ray Spacing:", trim(to_str(pl % aspect))
end if
write(ou,*)
end do
! Format descriptor for columns
100 format (1X,A,T25,A)
@ -943,22 +969,31 @@ contains
write(ou,100) "Total time for finalization", time_finalize % elapsed
write(ou,100) "Total time elapsed", time_total % elapsed
! display header block
call header("Run Statistics")
! display calculate rate and final keff
total_particles = n_particles * n_cycles
speed = real(total_particles) / (time_inactive % elapsed + &
time_active % elapsed)
string = to_str(speed)
write(ou,101) "Calculation Rate", trim(string)
write(ou,102) "Final Keff", keff, keff_std
! display header block for results
call header("Results")
! write global tallies
write(ou,102) "k-effective (Analog)", global_tallies(K_ANALOG) % sum, &
global_tallies(K_ANALOG) % sum_sq
write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) % sum, &
global_tallies(K_COLLISION) % sum_sq
write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum, &
global_tallies(K_TRACKLENGTH) % sum_sq
write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, &
global_tallies(LEAKAGE) % sum_sq
write(ou,*)
! format for write statements
100 format (1X,A,T35,"= ",ES11.4," seconds")
101 format (1X,A,T20,"= ",A," neutrons/second")
102 format (1X,A,T20,"= ",F8.5," +/- ",F8.5)
101 format (1X,A,T35,"= ",A," neutrons/second")
102 format (1X,A,T30,"= ",F8.5," +/- ",F8.5)
end subroutine print_runtime
@ -969,7 +1004,6 @@ contains
subroutine create_summary_file()
integer :: io_error
logical :: file_exists ! does log file already exist?
character(MAX_FILE_LEN) :: path ! path of summary file
@ -984,8 +1018,32 @@ contains
! Open log file for writing
open(UNIT=UNIT_SUMMARY, FILE=path, STATUS='replace', &
ACTION='write', IOSTAT=io_error)
ACTION='write')
end subroutine create_summary_file
!===============================================================================
! CREATE_XS_SUMMARY_FILE
!===============================================================================
subroutine create_xs_summary_file()
logical :: file_exists ! does log file already exist?
character(MAX_FILE_LEN) :: path ! path of summary file
! Create filename for log file
path = "cross_sections.out"
! Check if log file already exists
inquire(FILE=path, EXIST=file_exists)
if (file_exists) then
! Possibly copy old log file
end if
! Open log file for writing
open(UNIT=UNIT_XS, FILE=path, STATUS='replace', &
ACTION='write')
end subroutine create_xs_summary_file
end module output

27
src/physics.F90
View file

@ -17,7 +17,7 @@ module physics
use search, only: binary_search
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_surface_current
score_surface_current, add_to_score
implicit none
@ -93,11 +93,20 @@ contains
coord => coord % next
end do
! Score track-length tallies
if (tallies_on .and. n_tracklength_tallies > 0) &
call score_tracklength_tally(distance)
if (tallies_on) then
! Score track-length tallies
if (n_tracklength_tallies > 0) &
call score_tracklength_tally(distance)
! Score track-length estimate of k-eff
call add_to_score(global_tallies(K_TRACKLENGTH), &
p % wgt * distance * material_xs % nu_fission)
end if
if (d_collision > d_boundary) then
! ====================================================================
! PARTICLE CROSSES SURFACE
last_cell = p % coord % cell
p % coord % cell = NONE
if (lattice_crossed /= NONE) then
@ -112,7 +121,15 @@ contains
p % event = EVENT_SURFACE
end if
else
! collision
! ====================================================================
! PARTICLE HAS COLLISION
! Score collision estimate of keff
if (tallies_on) then
call add_to_score(global_tallies(K_COLLISION), &
p % wgt * material_xs % nu_fission / material_xs % total)
end if
p % surface = NONE
call collision()

287
src/plot.F90
View file

@ -1,4 +1,4 @@
module plot
module plotter
use constants
use error, only: fatal_error
@ -6,185 +6,172 @@ module plot
cross_lattice, cell_contains
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord
use plot_header
use ppmlib, only: Image, init_image, allocate_image, &
deallocate_image, set_pixel
use source, only: initialize_particle
use string, only: to_str
implicit none
contains
!===============================================================================
! RUN_PLOT generates a binary stream file containing a list of surface/lattice
! crossings and what cell was traveled through. A Python script can then be used
! to generate a plot based on the recorded crossings and cells
! RUN_PLOT
!===============================================================================
subroutine run_plot()
integer :: i ! loop index
integer :: surface_crossed ! surface which particle is on
integer :: lattice_crossed ! is surface crossing in lattice?
integer :: last_cell ! most recent cell particle was in
integer :: enter_surface ! entrance surface
real(8) :: xyz(3) ! starting coordinates
real(8) :: last_x_coord ! bounding x coordinate
real(8) :: last_y_coord ! bounding y coordinate
real(8) :: d ! distance to boundary
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
type(Cell), pointer :: c => null()
type(Universe), pointer :: univ => null()
type(LocalCoord), pointer :: coord => null()
integer :: i
type(Plot), pointer :: pl => null()
! Open plot file for binary writing
path_plot = trim(path_input) // "plot.out"
open(UNIT=UNIT_PLOT, FILE=path_plot, STATUS="replace", ACCESS="stream")
do i=1,n_plots
pl => plots(i)
! Write origin, width, basis, and pixel width to file
write(UNIT=UNIT_PLOT) plot_origin
write(UNIT=UNIT_PLOT) plot_width
write(UNIT=UNIT_PLOT) plot_basis
write(UNIT=UNIT_PLOT) pixel
! Display output message
message = "Processing plot " // trim(to_str(pl % id)) // "..."
call write_message(5)
! Determine coordinates of the upper-left corner of the plot
xyz(1) = plot_origin(1) - plot_width(1) / 2.0
xyz(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0
xyz(3) = plot_origin(3)
if (pl % type == PLOT_TYPE_SLICE) then
! Determine bounding x and y coordinates for plot
last_x_coord = plot_origin(1) + plot_width(1) / 2.0
last_y_coord = plot_origin(2) - plot_width(2) / 2.0
! create 2d image
call create_ppm(pl)
end if
end do
end subroutine run_plot
!===============================================================================
! create_ppm
!===============================================================================
subroutine create_ppm(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()
call init_image(img)
call allocate_image(img, pl % pixels(1), pl % pixels(2))
if (pl % basis == PLOT_BASIS_XY) then
in_i = 1
out_i = 2
in_pixel = pl % width(1)/dble(pl % pixels(1))
out_pixel = pl % width(2)/dble(pl % pixels(2))
xyz(1) = pl % origin(1) - pl % width(1) / 2.0
xyz(2) = pl % origin(2) - pl % width(2) / 2.0
xyz(3) = pl % origin(3)
else if (pl % basis == PLOT_BASIS_XZ) then
in_i = 1
out_i = 3
in_pixel = pl % width(1)/dble(pl % pixels(1))
out_pixel = pl % width(2)/dble(pl % pixels(2))
xyz(1) = pl % origin(1) - pl % width(1) / 2.0
xyz(2) = pl % origin(2)
xyz(3) = pl % origin(3) - pl % width(2) / 2.0
else if (pl % basis == PLOT_BASIS_YZ) then
in_i = 2
out_i = 3
in_pixel = pl % width(1)/dble(pl % pixels(1))
out_pixel = pl % width(2)/dble(pl % pixels(2))
xyz(1) = pl % origin(1)
xyz(2) = pl % origin(2) - pl % width(1) / 2.0
xyz(3) = pl % origin(3) - pl % width(2) / 2.0
end if
! allocate and initialize particle
allocate(p)
call initialize_particle()
p % coord % xyz = xyz
p % coord % uvw = (/ 1, 0, 0 /)
p % coord % universe = BASE_UNIVERSE
! loop over horizontal rays
do while(xyz(2) > last_y_coord)
do y=1, img%height
do x=1, img%width
! initialize the particle and set starting coordinate and direction
call initialize_particle()
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
p % coord % xyz = xyz
p % coord % uvw = (/ 1, 0, 0 /)
call find_cell(found_cell)
! write starting coordinate to file
write(UNIT=UNIT_PLOT) p % coord % xyz
! Find cell that particle is currently in
call find_cell(found_cell)
! =======================================================================
! MOVE PARTICLE FORWARD TO NEXT CELL
if (.not. found_cell) then
! Clear any coordinates beyond first level
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
distance = INFINITY
univ => universes(BASE_UNIVERSE)
do i = 1, univ % n_cells
p % coord0 % xyz = xyz
p % coord0 % cell = univ % cells(i)
call distance_to_boundary(d, surface_crossed, lattice_crossed)
if (d < distance) then
! Check to make sure particle is actually going into this cell
! by moving it slightly forward and seeing if the cell contains
! that coordinate
p % coord0 % xyz = p % coord0 % xyz + (d + TINY_BIT) * p % coord0 % uvw
c => cells(p % coord0 % cell)
if (.not. cell_contains(c)) cycle
! Set new distance and retain pointer to this cell
distance = d
enter_surface = surface_crossed
end if
end do
! No cell was found on this horizontal ray
if (distance == INFINITY) then
p % coord0 % xyz(1) = last_x_coord
write(UNIT_PLOT) p % coord0 % xyz, 0
! Move to next horizontal ray
xyz(2) = xyz(2) - pixel
cycle
end if
! Write coordinate where next cell begins
p % coord0 % xyz = xyz + distance * p % coord0 % uvw
write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
! Process surface crossing for next cell
p % coord0 % cell = NONE
p % surface = -enter_surface
call cross_surface(enter_surface)
end if
! =======================================================================
! MOVE PARTICLE ACROSS HORIZONTAL TRACK
do while (p % alive)
! save particle's current cell
last_cell = p % coord % cell
! Calculate distance to next boundary
call distance_to_boundary(distance, surface_crossed, lattice_crossed)
! Advance particle
coord => p % coord0
do while (associated(coord))
coord % xyz = coord % xyz + distance * coord % uvw
coord => coord % next
end do
! If next boundary crossing is out of range of the plot, only include
! the visible portion and move to next horizontal ray
if (p % coord0 % xyz(1) >= last_x_coord) then
p % alive = .false.
p % coord0 % xyz(1) = last_x_coord
! If there is no cell beyond this boundary, mark it as cell 0
if (distance == INFINITY) p % coord % cell = 0
! Write ending coordinates to file
write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
cycle
end if
! Write boundary crossing coordinates to file
write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
p % coord % cell = 0
if (lattice_crossed /= NONE) then
p % surface = NONE
call cross_lattice(lattice_crossed)
if (.not. found_cell) then
r = 255
g = 255
b = 255
else
c => cells(p % coord % cell)
if (pl % color == PLOT_COLOR_MATS) then
r = materials(c % material) % rgb(1)
g = materials(c % material) % rgb(2)
b = materials(c % material) % rgb(3)
else if (pl % color == PLOT_COLOR_CELLS) then
r = c % rgb(1)
g = c % rgb(2)
b = c % rgb(3)
else
p % surface = surface_crossed
call cross_surface(last_cell)
! Since boundary conditions are disabled in plotting mode, we need
! to manually add the last segment
if (surfaces(abs(surface_crossed)) % bc /= BC_TRANSMIT) then
p % coord0 % xyz(1) = last_x_coord
write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
exit
end if
r = 0
g = 0
b = 0
end if
end if
call set_pixel(img, x, y, r, g, b)
p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel
end do
p % coord0 % xyz(in_i) = xyz(in_i)
p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) + out_pixel
end do
call output_ppm(pl,img)
call deallocate_image(img)
end subroutine create_ppm
!===============================================================================
! output_ppm
!===============================================================================
subroutine output_ppm(pl,img)
type(Plot), pointer :: pl
type(Image), intent(in) :: img
integer :: i, j
character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
path_plot = trim(path_input) // "slice" // trim(to_str(pl % id)) // ".ppm"
open(UNIT=UNIT_PLOT, FILE=path_plot)
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(UNIT_PLOT, '(3A1)', advance='no') achar(img%red(i,j)), &
achar(img%green(i,j)), &
achar(img%blue(i,j))
end do
! Move y-coordinate to next position
xyz(2) = xyz(2) - pixel
end do
! Close plot file
close(UNIT=UNIT_PLOT)
end subroutine output_ppm
end subroutine run_plot
end module plot
end module plotter

31
src/plot_header.F90 Normal file
View file

@ -0,0 +1,31 @@
module plot_header
implicit none
!===============================================================================
! PLOT hold plot information
!===============================================================================
type Plot
integer :: id ! Unique ID
integer :: type ! Type
integer :: color ! quantity to color regions by
real(8) :: origin(3) ! xyz center of plot location
real(8) :: aspect ! spacing between rays in raytracer
real(8) :: width(3) ! xyz widths of plot
integer :: basis ! direction of plot slice
integer :: pixels(2) ! pixel width/height of plot slice
end type Plot
integer :: PLOT_TYPE_SLICE = 1
integer :: PLOT_TYPE_POINTS = 2
integer :: PLOT_BASIS_XY = 1
integer :: PLOT_BASIS_XZ = 2
integer :: PLOT_BASIS_YZ = 3
integer :: PLOT_COLOR_CELLS = 1
integer :: PLOT_COLOR_MATS = 2
end module plot_header

101
src/ppmlib.F90 Normal file
View file

@ -0,0 +1,101 @@
module ppmlib
implicit none
!===============================================================================
! Image holds RGB information for output PPM image
!===============================================================================
type Image
integer, dimension(:,:), pointer :: red, green, blue
integer :: width, height
end type Image
contains
subroutine init_image(img)
type(Image), intent(out) :: img
nullify(img % red)
nullify(img % green)
nullify(img % blue)
img % width = 0
img % height = 0
end subroutine init_image
subroutine allocate_image(img, w, h)
type(Image), intent(inout) :: img
integer, intent(in) :: w, h
allocate(img % red(w, h))
allocate(img % green(w, h))
allocate(img % blue(w, h))
img % width = w
img % height = h
end subroutine allocate_image
subroutine deallocate_image(img)
type(Image) :: img
if ( associated(img % red) ) deallocate(img % red)
if ( associated(img % green) ) deallocate(img % green)
if ( associated(img % blue) ) deallocate(img % blue)
end subroutine deallocate_image
function inside_image(img, x, y) result(inside)
type(Image), intent(in) :: img
integer, intent(in) :: x, y
logical :: inside
inside = .false.
if ( ( x < img % width) .and. &
( y < img % height ) .and. &
( x >= 0 ) .and. &
( y >= 0 ) ) inside = .true.
end function inside_image
function valid_image(img) result(valid)
type(Image), intent(in) :: img
logical :: valid
valid = .false.
if ( img % width == 0 ) return
if ( img % height == 0 ) return
if ( .not. associated(img % red) .or. &
.not. associated(img % green) .or. &
.not. associated(img % blue) ) return
valid = .true.
end function valid_image
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
if ( inside_image(img, x, y) .and. valid_image(img)) then
img % red(x+1,y+1) = mod(abs(r), 256)
img % green(x+1, y+1) = mod(abs(g), 256)
img % blue(x+1, y+1) = mod(abs(b), 256)
end if
end subroutine set_pixel
end module ppmlib

189
src/tally.F90
View file

@ -11,7 +11,6 @@ module tally
use tally_header, only: TallyScore, TallyMapItem, TallyMapElement
#ifdef MPI
use intercycle, only: reduce_tallies
use mpi
#endif
@ -1043,11 +1042,31 @@ contains
type(TallyScore), intent(inout) :: score
real(8), intent(in) :: val
score % n_events = score % n_events + 1
score % val_history = score % val_history + val
score % n_events = score % n_events + 1
score % value = score % value + val
end subroutine add_to_score
!===============================================================================
! ACCUMULATE_CYCLE_ESTIMATE
!===============================================================================
elemental subroutine accumulate_cycle_estimate(score)
type(TallyScore), intent(inout) :: score
! Add the sum and square of the sum of contributions from each cycle
! within a cycle to the variables sum and sum_sq. This will later allow us
! to calculate a variance on the tallies
score % sum = score % sum + score % value/n_particles
score % sum_sq = score % sum_sq + (score % value/n_particles)**2
! Reset the single cycle estimate
score % value = ZERO
end subroutine accumulate_cycle_estimate
!===============================================================================
! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history
! within the cycle to a new random variable
@ -1056,10 +1075,7 @@ contains
subroutine synchronize_tallies()
integer :: i ! index in tallies array
integer :: j ! index over filter bins
integer :: k ! index over scoring bins
real(8) :: val ! value of accumulated tally
type(TallyObject), pointer :: t
type(TallyObject), pointer :: t => null()
#ifdef MPI
call reduce_tallies()
@ -1070,26 +1086,60 @@ contains
t => tallies(i)
! Loop over all filter and scoring bins
do k = 1, t % n_score_bins
do j = 1, t % n_total_bins
! Add the sum and square of the sum of contributions from each
! history within a cycle to the variables val and val_sq. This will
! later allow us to calculate a variance on the tallies
val = t % scores(j,k) % val_history / n_particles
t % scores(j,k) % val = t % scores(j,k) % val + val
t % scores(j,k) % val_sq = t % scores(j,k) % val_sq + val*val
! Reset the within-cycle accumulation variable
t % scores(j,k) % val_history = ZERO
end do
end do
call accumulate_cycle_estimate(t % scores)
end do
end subroutine synchronize_tallies
!===============================================================================
! REDUCE_TALLIES collects all the results from tallies onto one processor
!===============================================================================
#ifdef MPI
subroutine reduce_tallies()
integer :: i ! loop index for tallies
integer :: n ! number of filter bins
integer :: m ! number of score bins
integer :: n_bins ! total number of bins
real(8), allocatable :: tally_temp(:,:) ! contiguous array of scores
type(TallyObject), pointer :: t => null()
do i = 1, n_tallies
t => tallies(i)
n = t % n_total_bins
m = t % n_score_bins
n_bins = n*m
allocate(tally_temp(n,m))
tally_temp = t % scores(:,:) % value
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into a
! receive buffer without having a temporary variable
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
! Transfer values to value on master
t % scores(:,:) % value = tally_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(tally_temp, tally_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
! Reset value on other processors
t % scores(:,:) % value = 0
end if
deallocate(tally_temp)
end do
end subroutine reduce_tallies
#endif
!===============================================================================
! WRITE_TALLIES creates an output file and writes out the mean values of all
! tallies and their standard deviations
@ -1240,8 +1290,8 @@ contains
do k = 1, t % n_score_bins
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name(abs(t % score_bins(k) % scalar)), &
to_str(t % scores(score_index,k) % val), &
trim(to_str(t % scores(score_index,k) % val_sq))
to_str(t % scores(score_index,k) % sum), &
trim(to_str(t % scores(score_index,k) % sum_sq))
end do
indent = indent - 2
@ -1317,15 +1367,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
bins(SURF_FILTER_SURFACE) = OUT_RIGHT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
! Right Surface
bins(1:3) = (/ i, j, k /) + 1
@ -1333,15 +1383,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
bins(SURF_FILTER_SURFACE) = OUT_RIGHT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
! Back Surface
bins(1:3) = (/ i, j-1, k /) + 1
@ -1349,15 +1399,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
bins(SURF_FILTER_SURFACE) = OUT_FRONT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
! Front Surface
bins(1:3) = (/ i, j, k /) + 1
@ -1365,15 +1415,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
bins(SURF_FILTER_SURFACE) = OUT_FRONT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
! Bottom Surface
bins(1:3) = (/ i, j, k-1 /) + 1
@ -1381,15 +1431,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
bins(SURF_FILTER_SURFACE) = OUT_TOP
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
! Top Surface
bins(1:3) = (/ i, j, k /) + 1
@ -1397,15 +1447,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
bins(SURF_FILTER_SURFACE) = OUT_TOP
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % sum), &
trim(to_str(t % scores(score_index,1) % sum_sq))
end do
end do
@ -1471,6 +1521,20 @@ contains
end function get_label
!===============================================================================
! CALCULATE_STATISTICS
!===============================================================================
elemental subroutine calculate_statistics(score)
type(TallyScore), intent(inout) :: score
! Calculate mean and standard deviation of the mean
score % sum = score % sum/n_active
score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2)/n_active)
end subroutine calculate_statistics
!===============================================================================
! TALLY_STATISTICS computes the mean and standard deviation of the mean of each
! tally and stores them in the val and val_sq attributes of the TallyScores
@ -1480,35 +1544,18 @@ contains
subroutine tally_statistics()
integer :: i ! index in tallies array
integer :: j ! loop index for filter bins
integer :: k ! loop index for scoring bins
real(8) :: val ! sum(x)
real(8) :: val2 ! sum(x*x)
real(8) :: mean ! mean value
real(8) :: std ! standard deviation of the mean
type(TallyObject), pointer :: t => null()
! Calculate statistics for user-defined tallies
do i = 1, n_tallies
t => tallies(i)
do k = 1, t % n_score_bins
do j = 1, t % n_total_bins
! Copy values from tallies
val = t % scores(j,k) % val
val2 = t % scores(j,k) % val_sq
! Calculate mean and standard deviation
mean = val/n_active
std = sqrt((val2/n_active - mean*mean)/n_active)
! Copy back into TallyScore
t % scores(j,k) % val = mean
t % scores(j,k) % val_sq = std
end do
end do
call calculate_statistics(t % scores)
end do
! Calculate statistics for global tallies
call calculate_statistics(global_tallies)
end subroutine tally_statistics
end module tally

8
src/tally_header.F90
View file

@ -34,10 +34,10 @@ module tally_header
!===============================================================================
type TallyScore
integer :: n_events = 0
real(8) :: val_history = 0.
real(8) :: val = 0.
real(8) :: val_sq = 0.
integer :: n_events = 0
real(8) :: value = 0.
real(8) :: sum = 0.
real(8) :: sum_sq = 0.
end type TallyScore
!===============================================================================

11
src/xml-fortran/templates/plot_t.xml
View file

@ -1,11 +0,0 @@
<?xml version="1.0"?>
<template>
<options rootname="plot" />
<variable name="origin_" tag="origin" type="double-array" />
<variable name="width_" tag="width" type="double-array" />
<variable name="basis_" tag="basis" type="word" length="3" default="'xy'" />
<variable name="pixel_" tag="pixel" type="double" default="0.01" />
</template>

19
src/xml-fortran/templates/plots_t.xml Normal file
View file

@ -0,0 +1,19 @@
<?xml version="1.0"?>
<template>
<options rootname="plots" />
<typedef name="plot_xml">
<component name="id" type="integer" />
<component name="type" type="word" length="10" default="'slice'"/>
<component name="color" type="word" length="10" default="'cell'"/>
<component name="origin" type="double-array" />
<component name="width" type="double-array" />
<component name="basis" type="word" length="3" default="'xy'" />
<component name="pixels" type="integer-array" />
<component name="aspect" type="double" default="0.1" />
</typedef>
<variable name="plot_" tag="plot" type="plot_xml" dimension="1" />
</template>