diff --git a/presentations/crpg_2011-10-06/Makefile b/presentations/crpg_2011-10-06/Makefile deleted file mode 100644 index a0bf7d5f9e..0000000000 --- a/presentations/crpg_2011-10-06/Makefile +++ /dev/null @@ -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 diff --git a/presentations/crpg_2011-10-06/annulus.tex b/presentations/crpg_2011-10-06/annulus.tex deleted file mode 100644 index dc603ec3c2..0000000000 --- a/presentations/crpg_2011-10-06/annulus.tex +++ /dev/null @@ -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} diff --git a/presentations/crpg_2011-10-06/beamercolorthememit.sty b/presentations/crpg_2011-10-06/beamercolorthememit.sty deleted file mode 100644 index dd1943f312..0000000000 --- a/presentations/crpg_2011-10-06/beamercolorthememit.sty +++ /dev/null @@ -1,22 +0,0 @@ -% Copyright 2004 by Madhusudan Singh -% -% 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 -\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 diff --git a/presentations/crpg_2011-10-06/graphical_settings.tex b/presentations/crpg_2011-10-06/graphical_settings.tex deleted file mode 100644 index 3cd943c5fa..0000000000 --- a/presentations/crpg_2011-10-06/graphical_settings.tex +++ /dev/null @@ -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 -] diff --git a/presentations/crpg_2011-10-06/mit-logo.pdf b/presentations/crpg_2011-10-06/mit-logo.pdf deleted file mode 100644 index d64b43adea..0000000000 Binary files a/presentations/crpg_2011-10-06/mit-logo.pdf and /dev/null differ diff --git a/presentations/crpg_2011-10-06/slides.tex b/presentations/crpg_2011-10-06/slides.tex deleted file mode 100644 index aa46a38931..0000000000 --- a/presentations/crpg_2011-10-06/slides.tex +++ /dev/null @@ -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} - \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} diff --git a/presentations/crpg_2011-10-06/union_energy_grid.tex b/presentations/crpg_2011-10-06/union_energy_grid.tex deleted file mode 100644 index 171040b68f..0000000000 --- a/presentations/crpg_2011-10-06/union_energy_grid.tex +++ /dev/null @@ -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} - diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index dea486de22..a1b474f12f 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -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 diff --git a/src/Makefile b/src/Makefile index 734a4af2a3..bb25277980 100644 --- a/src/Makefile +++ b/src/Makefile @@ -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 diff --git a/src/OBJECTS b/src/OBJECTS index 1a9f2b8cc7..734b6da54f 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -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 \ diff --git a/src/ace.F90 b/src/ace.F90 index 13e4caff7b..0406c7dbca 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -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) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index a26c607026..8de22310f8 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -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 diff --git a/src/cmfd_output.F90 b/src/cmfd_output.F90 index 45b2a7167a..3d67b9ed08 100644 --- a/src/cmfd_output.F90 +++ b/src/cmfd_output.F90 @@ -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 !=============================================================================== diff --git a/src/constants.F90 b/src/constants.F90 index dce69efc85..89490eb183 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 diff --git a/src/finalize.F90 b/src/finalize.F90 index 93a3c7b086..c4dcaf0623 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -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 diff --git a/src/geometry.F90 b/src/geometry.F90 index 856b149c51..3457719556 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index fbcf0353a9..aa1fcca5f6 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index 49c303ce29..5ed20ccd12 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 2732c334cc..4c5091f05f 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -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 !=============================================================================== diff --git a/src/initialize.F90 b/src/initialize.F90 index bf07944a2f..db2ca4b585 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 18083ee45b..c7adcae12f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 = " 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 diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 2861f0ab76..adfeaa0288 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -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 diff --git a/src/main.F90 b/src/main.F90 index f2fa99055b..5f75d57f92 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -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 diff --git a/src/material_header.F90 b/src/material_header.F90 index 916578923b..d8e9ca4ea2 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -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 diff --git a/src/output.F90 b/src/output.F90 index bd5fec3381..dc3a9f4fbb 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/physics.F90 b/src/physics.F90 index 4eb9852834..8ed65c8c7f 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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() diff --git a/src/plot.F90 b/src/plot.F90 index ce7b26ed6c..cd2aec0caa 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -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 diff --git a/src/plot_header.F90 b/src/plot_header.F90 new file mode 100644 index 0000000000..bfe6f9b071 --- /dev/null +++ b/src/plot_header.F90 @@ -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 diff --git a/src/ppmlib.F90 b/src/ppmlib.F90 new file mode 100644 index 0000000000..2c9687617a --- /dev/null +++ b/src/ppmlib.F90 @@ -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 diff --git a/src/tally.F90 b/src/tally.F90 index 9920b8ba21..af6d93e738 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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 diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 8e04a334d4..fce21943e7 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -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 !=============================================================================== diff --git a/src/xml-fortran/templates/plot_t.xml b/src/xml-fortran/templates/plot_t.xml deleted file mode 100644 index b4dc9eaea8..0000000000 --- a/src/xml-fortran/templates/plot_t.xml +++ /dev/null @@ -1,11 +0,0 @@ - - diff --git a/src/xml-fortran/templates/plots_t.xml b/src/xml-fortran/templates/plots_t.xml new file mode 100644 index 0000000000..af8f2c6e27 --- /dev/null +++ b/src/xml-fortran/templates/plots_t.xml @@ -0,0 +1,19 @@ + +