From 6418ad2ebbead041de7db68659aad7c5f2f52c89 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 25 Feb 2012 23:24:14 -0500 Subject: [PATCH 01/14] Moved cross section summary output into cross_sections.out --- src/ace.F90 | 6 +++--- src/constants.F90 | 1 + src/initialize.F90 | 6 ++++-- src/output.F90 | 27 +++++++++++++++++++++++++-- 4 files changed, 33 insertions(+), 7 deletions(-) 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/constants.F90 b/src/constants.F90 index 9168727771..d2fad27b82 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -318,5 +318,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_XS = 14 ! unit # for writing xs summary file end module constants diff --git a/src/initialize.F90 b/src/initialize.F90 index b3196ae4ad..697273533e 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/output.F90 b/src/output.F90 index e1226e1ea3..5da98d1ca7 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -968,7 +968,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 @@ -983,8 +982,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 From ec78616d6944f3d51ee908b36a87dcc62c3e4ef9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 26 Feb 2012 11:26:25 -0500 Subject: [PATCH 02/14] Added hdf5_write_nuclides. --- src/DEPENDENCIES | 2 + src/hdf5_interface.F90 | 118 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 120 insertions(+) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 7b18695de5..121e635e27 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -72,7 +72,9 @@ 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 diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 07289d32a7..006c6cf3c9 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 @@ -92,6 +94,7 @@ contains call hdf5_write_geometry() call hdf5_write_materials() + call hdf5_write_nuclides() if (n_tallies > 0) then call hdf5_write_tallies() end if @@ -598,6 +601,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 !=============================================================================== From 0bb07fcdfc367de396d3fd4cb088ada929f43f5e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 26 Feb 2012 12:45:51 -0500 Subject: [PATCH 03/14] Removed presentations directory. --- presentations/crpg_2011-10-06/Makefile | 7 - presentations/crpg_2011-10-06/annulus.tex | 13 -- .../crpg_2011-10-06/beamercolorthememit.sty | 22 -- .../crpg_2011-10-06/graphical_settings.tex | 47 ----- presentations/crpg_2011-10-06/mit-logo.pdf | Bin 41503 -> 0 bytes presentations/crpg_2011-10-06/slides.tex | 192 ------------------ .../crpg_2011-10-06/union_energy_grid.tex | 51 ----- 7 files changed, 332 deletions(-) delete mode 100644 presentations/crpg_2011-10-06/Makefile delete mode 100644 presentations/crpg_2011-10-06/annulus.tex delete mode 100644 presentations/crpg_2011-10-06/beamercolorthememit.sty delete mode 100644 presentations/crpg_2011-10-06/graphical_settings.tex delete mode 100644 presentations/crpg_2011-10-06/mit-logo.pdf delete mode 100644 presentations/crpg_2011-10-06/slides.tex delete mode 100644 presentations/crpg_2011-10-06/union_energy_grid.tex 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 d64b43adeab78f3df1a0d06fa5430531a1d7a94a..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 41503 zcmeFY1yEeyvM)S@1Og!t+$FfX4I13tJ;*S)3@!l@+3Ze5?dn+L_R)5CPhOq z$bp3Qr9~N$N!85`M8YIu=xAtV3q}+WKm^&C{9*h<{%41O@ch*ckx9nf#DPSIgzZmv z5H3HCzjpnj7Xo@={adjz$idbLV*G;S4<=0oBMXqRwEPp(r{4c2gQ2&Wq7Gwi5u3BX8hZSfBOFs$;R?eul_OOPm2x-(@SVbm{dWoj(-e(iHOLHGym?# z^*7gxWBzAA{wVw>zJCPdZ=WGDsW=%q{uRc*1R*lX8(P1_^u>*cOycHNjv&ZOyqH8q zR7F9?wkDuILDMidd9m-s8bl^ZkU7}wB`R#(e`0B63sJE%GzO8d{4w*zTf*j!4vHX% zh^@7qt<9hCy|Df59^hXbF9`p(`g!(w74Sh?TuK}O0|Nt)eE9*M7XTuE&i-55A0huO z0e?C_{{*1E1suWR!NQ;dUZKLkqQX3P0to*F9N-@t;yZ-5@NjQly?!C9d;kF8UctPA zeT4vz_zE5l9sq!S{RR&H)ms<@L^N~^OjlGaY+P=jOI!-CsL|zofo>A@ZWW5Z=JNLWDzn z1^W+{7xk~<-f-5Ueq`rxaCCC6s)5H)jUAsr!@_12Qc_coi6><-w9kB=2Oz=1yqE@y z3J?Ue`TVQ=|Nr>^J03u5Uz(^$II6u^%-PZ2b=A4Rn459gcTLIZj;q>^*DV~TZVH%_ z@E^-P!(p31Z|c44ai-)FlikC@o6qb^D|G9tR$UBR7vV8T3oWueq+fXkY>*kXEq+G% zPtGw5=ifK!gjU_xR-^JqEHEU@;#4CSXM8jc*Q;;MwP0p#JGE9nY!P6u9^O(I8X6(s zD$`v#WDoAs`8`b|vvZ3m-@!Z@2Y_R|=bBV+zZu#+^(-w&#}FX7WNs0i#Aw+#b8 zP`90Qw`Tyg*TL&g*Rj?S6Bf^ap}t#otI}scIe|QIVI`vNqYFc*vZvim^(~X9?#~Y0o={|q#VSS`fAJMp0U;Bp)UToEsd9Ai*K(~08cf(KO|Kyx;q&mag z*#T^6eWF*Dg6nZW``d5bR|-wXbz%B_RS{$SR=@;siQP%`O8Ci=)I}vFWe7si%${4u zg2*?y>TDi-@+KQ_R(9TNu+AgX5Yv#o(P^ZDGZpR|)tythLs3$Gx^txs)!iJv={}iv z2V*5b6EOM~@d=`+*CVLgBkFt;2;vJv#NRTeWnWH1%}Y4Y?uOX3%C_Y&e&W_K6MW| z_e3qcDZANt^!v%oIhDnvdA_q;!eT@hh=YEn@@1d2O-4yzU4f>vG>{SxV(Uvf@!7up z%4Ykp(rI!y>bf~qUUC<`(iVZxyWDAGpDhNOpekPtZpnc>L~)ZT#<$ZBTbRj(@1h9X zTc8)7sXhZ}UZ=uRm-i2tKLbjB{Peb}W4a4_^V3;VXL>5b^b$bYUqdg06yQqPk~%KC zZA}4cX*Siu~nkU$q63x-AR5aGfJt4>-|6Z>YT&ebAg)q znPNoZ)IgZ@8f!_n(kW-hXQDglqj+$P24n8BD|wNSdHE2?QeS=M8FfLsRGsk3=QS?H zAY6JMI;o~wM>2tE2+6f@qurjAzjVZsZ>UA7dYkRr5g^_o0Hp5hDn@3`iuf7E)|nOB z0jSmw!{@L)5|p);pR+uDM9N#LHDcW{FFvR;1)wL*{V^lmaPrf?c?Jz;x$shR)aw~=X>j8(ABF4tLatu!^Uj~1x&Kc6050}1 zcLu<|{Ze(Oa?sZDX5c223})v90KlJDBQDyTIe7Yb?F2x{0{AfAj!cZ7UmuNJTQI+U z4FJGo&KQ^P)fWWL!!PMjxUg-lW#0OA;ZKLj&j5Ij*RTfvydiy8Fz-rt zSFXNRW;pKwK3q$UCF=0ZHLBEIRmVv>HL=O)1_Dug;7mm^x94oKLGN{=dvNii>BzAO z-rov!Z9HP~sZK-;yk2(JeItFPZgtGeKQQ-2&B?Aez}hc2*-COxrMibzRlzLb5(g_7 z#_fDVzO);~L^c(=nL%K$nNFWD?fUo%gJiDiGlG<* z=>02i%h1L%?X4uI#Fm#SUIVkmL--z0QM4|pkJ$I#uqQCiEcq-C>0g`d2vZuBjIM(XNzI-_;P z8ic7`B+lj*mpAici_URNmhsheC(%2}Ui9|kTcGgT4Doal6nxQqh_3xzu0$Ht21(td zg8PnlT6eMCv!`a=go@pgWbT2LB?U3rrCiNH{CGM(C+`wR8t+C)EgjvD41o_w;>uR6 zL|{ZJZpG^gI-!tiPccv#-sj(+YgRUVKMq_dE7)hVtd1pM(PgQuY5vh?_E96(CQ)Wb z-g-adqh;Cv)%^QP?StO2FdAL{xT{kJ?MQe{&-}DBB>!3(M;x;i8UXMXMhcpK5~|xID>?q z#~h)Dgk$H&190QOHERpltTa}=-@(Q~jnJxoeRsgtr#F>~kzu>0N3wLHb0zwgR$~Vt zsyf@l!w}v+(*80&v+i3s;n)w--MaTx0=zL{NNywAae3@0uGa-Zgxy!Hz422R_*&7WT*r`=oPB`kZ0|WzIx;- zu~iI*CAXsza9}+*XV>Q z#2WWK_LoGSzF zvMC5NMXkQOj?)0PeaM>5_M7`qXm2B>Nj0_|$T2K4OBn28=xClS@#XvayuaZf6a7Ad z1o=%XAxiV*uYCy}N3TpxV(CcE`4k2_f`bXMNjidI8;b!r$1M-{fr|LdjfB@u2SU=@ zHQRM{IPrI8H>=+`dn*fxyQQpT*)nR}pDybbEUEai7mOWy^z^amHSmirzhv;-)|S$U z7O14ika0N9oc{XcUrR)lwnw3$TpuZ~K&wlrx-F6B;|{Ojk6n{FsmLgmH~Bv{HUl%(UQp3=p%Er(k;>q%koXZxPdLaWn9uR zdk7_I1kQOrE~;xVl_@UXeDg~cvdIeXV(SPCsZCC_eBIS(du?OI3GIN0XA{M$U(E%Q z)-sK*Bq*o$L<00U+`Fv=rTZ;y*P0q8aINaD3Y`$s7+NW9r6_LcB2{*ehV{b3Mf6mN ziCq&Uq1yB}wM|leu921BL9XqxhHDTO3fT!>1O!|v&;X6Leb ziAAL5iw7ECTUDb=fC!qaA^Y!Z5>r?=S3YN^irPWh~5Gu$i z?n|+q0*0h*W!~xYSPhL%u!HO#b`c>z=4TUEO7_DqgWcK<2*M;{%<83gYtAx)rC}F0}T8!w!f>hj&L5doJ}0f`}6X(`gco^sKMW^e9@4 zknEePe&rbco}{7}@2XN|ApsTTR5-JSNw0*OhuO`=ZghNOms6IqO7oD_o1~d3$ycAx z$WgDvwi&=fsIYm5+UF_klh5L~_I-D;6j60=IbA;iDxGW1j6+LNpbHv7P%7Y^Y?R10 zua@1XvR0Kez5aaJv)?$$xu}CzaGVz1Ek)7Sr>Y+`);m~4S{a{$1UZlDp5Lvzhu|N5zbW{DO5fKOfwk1*H^JH#+SO9 z4iLD>$DdL(1Y+aPt{ zL611sj)^)mo4CI5b{cKV=yV&>=y=Its|ir6;Fm#t&)D4g8(b^%bwopIMUU-m9?y*eh;Xm6P|SFcXj z3QJmjgwy-=x^k3tQ~K6W(F)ZKaH@FwjJE_6on76Lemd?I%GZXaEjCu&o&uMdz~ z;Cys>weY^5slcn(=~JTcN!JJ_y8Gobz!0S%c#Wo?)p3V>v?{cCVsp``nuUvAF(Fo; zJC_P+;&V9BYj{m}Sr<+x7qtpC{RY{r04D+ound0I!3k!VEGqBlw|G#1WqcQXZ?iF* zfYj_luYHsNkB%HGoFfB`E=G%VRtfd6i`9E>G(3o3U!-DOCsk@{tb(Y$AD)Fs7!%T< zw0U6|kR^a>BwdQRbB91G$#8>wa9{`Nr}y+fGi|`DPrOqC#ds(yPsPswWDjLn=Fb2C zmd6~^#ohcf!0`=n*rnZj#-A_d)EP*BNx1fD^hP!L%=h$?&wEH`U+QT;W^@W1eTTXD z^x@44@XlS|>v8VsP#gx}0|2n|95_9yUw*j%G$YvTeDm%H%*IhB@h9u`>rZpG?MV$t ztS>M;VDrl-jwQaFaAcomVZb{00C0`2V;-$f%N_1)kGp=pQvCG7_1ETdT|2TrJ=V5V zZLoU$*~jwgAI+~mc&xXnwemiVa-RLAZUvGmYDoSR`db59sV3k;@mu!De)3HMc~dHN%hM1g zhklWuRd{iaXpUN}j;gvC{|tPo zTv#uGN@-PdjYc3|4tc)Y)wb~uUSmh=uCEwDzbLs(f`5qKVttGIt*VwQVqzY{C}f@3 z+Q`!cj)}99f;JFUy4@BzDptzle2|)@BxgU^rI71&W~xJOnVAilEUjvuk*eEE2?H{6 z#VB$0T3rw9q4n1q&OX`Gr~vtM6~7O9b63fX2tG6*RpU?U&Ieb7BTH*772RL!9Oi-7 zwQ|Qc8=Kn?@v79L*IE*0JA={cF2SsQxac^Ag>#7AQ?sajw6Vpa&YC!w8H(2#)-Fha z!dJ&~wJ9ba84Hk&<48lPG2EEMzXV#`$`%TL8Wjc=V8Vo@dVES^cnP%G64jWbUTv6f__Oz%Q zPf)UWXK)WAn?h(DqWt0puI5nmVwv!`sM`d4XWb9>Z9A{>d_hObJICW}xUNRDF)zK1}Xp+W$;T@xSc5#v) zEu*0G;dltvQ8MqO1?lPXZGWff&oQ{A&*+2`1j}%A_)+g*k!rQ1^^>7|ngXbc4hGX* zWM-amdx*%|T|o>mxp6L{5opFtArI#gUuOaSh2H8<^G_-Gub03V+q)|itxWStCh-`; z8C?EX39JQz0Gf$%qf;`8<5?RJ_Q9R=I6x1UqA+^RT@M zwGL1*oHY0yKS+IKWiNpCvC^C7QAUSIIgvg~cRL)oOce+p2L^8EK%IqkL`s@F#>nz} z6c2q0oOw-)RB~cNauKh3$j5q&N#EDM)1a&+!~->|Repd9F@gM&znUu5n6Katck%wZ z%BKk3==LCZo$gXRgZU`+ZkIz>Eh;F@-k@{gTCRu*(YI$lYfSUrTrfrf8l2`IPfOW- z@c@E~n`n!6Z6qaa+jz0D4G4&=>m@XWOD$bLZCbrA1_aHbS=$ASv&^@$5Bb+g2GHI* z@@muRZCeEqj+cmjQ(sJ!7S4%Y|o;y~UDiTg_stLI+p_@%b zwD#&eo^8#s5*C=lmty2hvdhO3-HZ&@3YLr>SK4I=3{b0jBL#8Q?4KaA(rlH6bI$37 z9U~0$&MzW99=APOh$8rQbauCymwtSq@C+bxg3-`RX;bx!TzR#AAFBk2(?4b^lQw?| z$e(%moL z_+(T1+X>q;*^3jAPM{AzPJ9A6W296e1H7%}hs^Ko)B(0R_fix!YY=rJ}UptVG zPh{&*KST<{oj!H6O+q({J8yu!-wsCc0^6b zV{*%uygGSN>{bBx3QzWm^1v|rN&6sTAl=CYge5#PPr6On2{Y&-4Z7e%Frr$wyt0| zs*w7%uiU2UU;<*)-^DZI6~A|swCt{aWDM0WHChd|wlcNK&pKXuoCR4r=vYwaS|p`} zdsnLFM|3YDdPTWb#)3Lh29u29-}YUYF_Q}>Os|xSeVG=o&NmKmiASARi^B zJf^r?*wEj(SXl}6dBg8)^Qep_==m=6NKW+CM~|o<&2lee&l$xdKj$0ZBKhlAl=c0- zWM~S8PbVU+GzW4RElx9*;`=y32XAz}ZtJ(3->VF2?hDD=j-cY zjxBw0KT2(ZC{FL<^0TH+8ODGy*ADn$kIGlVv2*AMtyrsOLVu+eY?v9I#DK<<*oA+Y zht`6--&-Zhd_QyUJ)8Wk1Tm7;5csrA9KRC+vSE6sx|+}NRrHb8!~agm;NIexrPA>V z|2uwZWlm-LiYX=^y_Jba{Q^#A#K|p2?a+uXL0ZnJI5Y{nI#FlanJz2w zLNxW&|4a@607-cF{3&Ch`)9zE_u7cXCsQN<;MIr+#Z!ik!I4k)X?htEBR`k}Pkaecvq4<|q2n25Q`9-je9OFy6R7$Ywf;fMCyyeD1k1Oq~| zA2DFg-H%r;2^Ue;vVK1U!~mDhfo9I9JzSBpo2-Ov^V^hg3}rJe3H9}9&h^RjxXE~j zzujl3G&CqG&QQ6yxJXGZixU|5;tJ7XEd1f1^0+%ctuX+w?2LzBl;&lMi%i@NNmjx5#x zD@xDh=IciZ7uXxKCc&{!|0M(e3teFwm+J8B>zy>TOPMmX@k*;OTQCF9La#Q*JtVY7 zZ?c(cIRvgjR#!CLl+xqqh z=4Ik3s8kJRP+)yYe?*1n=V7OPX&YrjZaIrY6|)KgU0mUab@R9(XY}^05-3L5k^ai4 zaM2!Biu+=Bz<4#rt>aruRma97{^jnu2T$!Dg*{D|YE6>fUC4Sp7I5v?BGI~9^bxjN zeB4)|A85?o>I;|<6At09$~BSfGgZG;_pRVq0t~L@^Naj(lnQ*5Q&eV@mH&_~cR%g! zSB=%&4IKEAlv=s-q-_L)X1cgbtU_eH#J zBB6pEw*EA?6=>%UYeDsi+dv_x`_e{tM~A-ow9;zV>E}e%M9_T-2X&qxwWkP?AM<|4 zM1XY+7(OW5(OI&)aUPooDw9ss!WLshG~>Imu~1-mw*=i6l?aBhh7LTnk#TsDSkCJ2 zniU5KgTOaf`b8G1qXY+LQR!p8Md-QCKqxwUR=ASyqy(os4Xh|ZHCpY z;X&JHvU&6n<-4FJfqUXrVQfvUjjUTr=(EKwMPw^GZU1gUiU^C0jCmdV6A0^-P35B? zTd4v>{xSYy#<$zrbW>?@Mz^BP!D_vj9?}z`?N(fwXmS~xF?2Q&GVJiM5*UTmJgL!$ zoujio{&95w*kI#6p@?BPq6B;E$b-P5%SXc*%|?KG5wHsmS?GZx6I{L{)}sJw%hO=uR3odzKL#sdE9bwoQ1`s zk3)qC5=tP%9v2%{3Ej+nfaC2b2Xo$#5vv2zdFd6;9gEdGYEg_?Rd12;H2G9SRV9Ju zlsExZ?YV69&<|X8J$$`xI^$vrmd?ZSki=>Gki6I#)2|8y1R2cODl?`hbKEG_(CiQkW(p`sv@*7yIysdcV8K#RPT0-Jjke~?S&cG?4CGLjc* z{vIB}4-#RkTz^RE%A^L{Fe?6VJz3u(WyyqeV|nIK3C7Y|SfHpli+rKE_CA1Xc@R~G zP2r!#n`6jsYS+0E$r4NbRTqM-rK?D#@|JBJzeANUpfsuKq%DGJ8PkkS@=7`fh>E)k%}+ABq#FpA!PfA!rlutxbv7`WE7GOYveFrH7-B>$ZKy zAwz-DQ_ieNMkl7|k{@|t<=l)o$fG+o91*8*oM(gt(07(N zeoSg8TU+n4)wmzYK2m)7=(>V47k(v6GmJ@U8Pnr7vG_8djmWCmPA)&4t7V29BeaQD z5M1U5WTW&3u0$wL*|c?;aPQPJ z+6+;nV*f^XIUa6r!#vErfoT>hL25B1#AiRkcc#zyJ~U?QDNconce z?5oAIQK`eLw11dSFiQtFOjEX(Sgk&UCdrqDt4m?+94bv=k-TABhSS(%fUhbx+C)fA zy=hXomced0O0A6kaS_xw$$04`QEOvu5s%AXm`rk+YdXfqP}^RZ_jwmF@|G3ebsg+j zpeG+0bJl3g>k-pQBjb9dHC*aRYpgxn`dQ1xTr9z0;tPf~uliJ?MDg8v+uZ%M8hdiT zZfuQ2w$fG58_B}VDdAarZVvkBQrVlH_5$MI8< zxK!qxycpi}Qfuvcv74K#L?mKkG*yg%NOavI±+m*AVIb3BuGc0M36C5mZ@hUgwR;l^74)>2Ws1s!O+0X*hlT5{> zLyQ9%TGrcAWj*?^;Ef_P1X&?$1Ofz&Fyj0hJS$8k?&&?YmvTkbK2}?V2)gYb&j#1p zfVtu+ct>E2&c(9#87cnxfp!&O3PZHqh}CC&2j z$%w0oDLK*-Y8_j#wm})%{dh^9ER_4o<6)&gcSVV$mv58oxg<6%pc0$(+oCB;b`~-B z)I71%T-HsmFV4Sj=IlQ~FA`?cx}w%)I1E&}YC|#l6YW*t!$C~txtWr z^pX7uSnVmErb{6rMzNY5C8g3%v_E8W5GHX-qg#J&qDg%-#jdr^WNqC{1P7|m_5h8c z+)4R0+n{qO-)Xn2Jx@dE4n6D66TRmkTSsc;!wL~_?c@ne!qYR~Iu)CLRcsbgsyq87 z09yXjBo{o944_mDKbrLM=>*qm&7?AbSHIy7{oUIuebH9Hf@58PoyIqj(i3 zgZX;Ig>)^Q(As4o>~sJbl0M^h+O2sW$+rr^RZrB{wI3mH7zVCEfg1_LeNdor!qXE# zd&NdV{YUyxnm~Uo!j`H{P*;TCS@19BTJ&`n zq@3Ao3Yul<^XLXj?QmvK4Hw?1CiZ1~vka*^34<&Ttz6ly#wuCo+;x(W&M=YEbG5jc zlqF*{j$qNo?zYIflg-6bI}JvY{!}R8zNNDiyCExABZ+6(@ho{ zn$MGEmOn~%pMw81ufo*zb}(vKUfe{>RHoG#L1wiUp(S;AWcaM3hac9%l`}sFx;4Vd zv)puf*B>}k+ixyq+M|_9HMMOtZLzK$Jc7bl#3=$1*gIlD<)+2>8mef8BXM%)Xj znqHuJ@#Dh1g_Bbq$MWw?2HKeN6m_-Dy0S`?+86SB|OW_8j{gj>?<{Z|E77BQXU0y2Rar<8Eis z17n8`Z=YHp8-jT9WwqYU}CjK)mA$8Basp{D;0mu^&NKm8vWcnbDIyCWfG52iuL|!KkatwK*1#cq&7g(tB$K?qB1=#^yOp z8p^=t*QIwIo$OjNj3+ z^x&)s-K-sg*1}Og?@+q@`gzuZn6CV13w2Y-TY>)7`#0$oTstG&P+x@BFQ^G(JMRz% zVg^~<6H6_`rOvp1#&VDwe*1Z{^<{2gDogW18P{RBL~ZbaTIN^J97Ds?W6#^c!DV&4 z`MKBSxro#Bl1hR@WO}i}PSUEBWHWhG)AXWPB4gS0zS?&Rn}_E7QW3LDdC)GkPW?kE zBqufuX-pjVE=~w_t{yyxdiijqXk&oFxJh_3~hKRHl69yF`v zm%cHo1I$o2RaOg0TPoO%rM*u3qJC1%D5q$`(#@weeiz*;+%l}08IMnY#~De+Erb&( zrZqt+O>M(UVCx>)@42$3F#WL`+(UJv7Aj^#tEjqLUGkKAfkn%D_!v>OZM+su*vf(D9^WkR}5jKppT-Sqd1Kh$c3ZicZ9RP%B-6?blj&e6?L zR%?&0=-MDnI6515vUkki@`j-md4=|qigL3xnq<%`*jbxvQH#e7B2(aZlyRfkuN_I2 zn_NzG*GsC&T$PoK*bpw0xpvLs1t-Lpj|R2VeG#pgNivM;*f}rM3E`6z9te^2^IhhIgQ)y;J+v^BCY}mo&=x22yT9JX@I4HFK&VJdZB;4;@`hOO zhzMKN(8a?;e%*}$*NdP9LLnuR3aK%~CLfl2JAa$w+Er}QJM4zKz1~O?Yb{h+CHZ+7 z5bSJZWtHyOL%sLTPUF%^aeh@h!M!`$TMbx8u^W@tEnZy!rTZ-Q`fGhd!rDn2XJH19 zkwmHll64W=N4JlQ%k_FKRm0ygxfJ%eq{!pM>(#U+u$v7JM|Vfi9fxjNG&jd7BMzSQgS z%oNa}^evoyy%{?zF_}TJ@htNU_Q0Oy6D2(}2%O8HL|&e}3e4JwG^;y_oO}jIZTdJL zUdXn&^E_!1Gjg^HfuEEFW3RisV`3>A?@=}*x9S)4m)8xj!##;6ZWYY045}Y(3GRUd z{^==qTvUPgNASNtC}>a?U2Wv%v$%Iy7vao}io4Jp*9?>$^mo48v1pp2ja|W+S3>x9 zVNygHM44xlbqRicTcokId{uPqQ85dByH4q@Zm_%5FQzX< zXVEXS0EN|EOK7$~H7a$w1#o4v)FMz3M!YW9p9*kGSM>0jKY-4T6B!p4l(ys$W*27K zpKUQ=H9Nl4s~#FUgk*!HXE!QwJS$sch`$DqNf3)y$zPo&(a39VSEI+>nY7Edk~Og# z<}I<)52t$+cLA&HdG%*AWb^2wH>S`Cjw2460Beom6nH(%9)f=jGI!5m9U0UdRuaB>qVOCqR(WvfIz4-Io)8~oWI_M$zf<#lbOQ2v-(+u^1<{(WlV7p z>&4jMMlA-MCQk>>Vmvorf8Jsa5TaULh+8UJM`i#;V!i_VEw<4h_gqr$)|T-r zU3RVLa}s76QlM2QQP`S@-$;ud<)=l(Z3gAqi-q*@vuvfOo~M|gUX#j$W7tw59DXqs z&>70n5eMUY9-i;&XxFbD(GLl$qhK|N!@Pv&fHEqwUlM|o))Q|XIeYtB!=rnDopw${ zb?;1=_@AZ|!LzZRu4;O_X!3ILM9-b&t@Nu~O6|XK*Zy@AMNg}$*JW%H=+C|G!jhOn zoDsYqhm(RXOP4V4F%S*XTaHS#vxBB-CNwdTZ&PzII0NN(nht6l8Z3;s>YTQF2J}-* zQxvmthrQL~s34Bcv}G8YJqu-^bH^|nW$BrhlpX#xKR=)GdKIWKlkuIvrTA2>v|Mjq z&arGn)xNgEGzlp;$?^cMqbv^C(EH)nHUaxi$A#U^LHyF%9v3Z0`Z7d0Bd0V%DKOHU z11Ry$+)RIHI0tj+P(-rZK zji4oeoiA^{=t(Qzfii-;;xu^}cRS2v+kdWYL+!*?fkwUm+F~|bYlta?u9%lPR2k<$ z+PD)+8@w%r9?#a9qGKd4=BX9O(Tl5Pte==T-|a%d`eEBvj%a*oeTrdm<}hTI&Ov5d zO-+TP>k0fOptp3lIdz(+Fli}%8<;{rAL!pbnQJe;8{;IUr;9AyNK>RjL)AC3I>=8r z=}~Q3)Q|_^)ya?V?kF2p=IGbMvb-#sy*rm9gXkp%7WTgbMnF&5$`f4$Z+{sTlRCd5 zZ=~ghVn%P9yv4CYMkActp9a)FB%LS}@6_AbP@K|H9_n=4l5fs_*8NM${eMNu|1Ws{ONR)6W7WG7 zGePgW@(19dBvQK%0HAej(l)Vg`8bN<#WwYmy%_-GCouPg=E`9C)z*Vpw9oks0Py~h zHgd|z^2+JnIHBzH-JtaUkW6NK* zvoduvV*BG6upqzkAO?O0X!2X^%+uZPWuWfYr_VSX?>*SI?p5z^dkB@@Y@Y2yfd;Bi zK;zhu!G;ln66#}ujRb;>6Ner?hnGLWz7YaXg3u+-?QNcLKQeo5=J3N3vbWP2-vp&w zQkj!cQ@$}rZgg}Jca?9=XmGH|CQB!9q4x>!m-O*fPp6xk#_)g!V`bQ@yUmU1O~OE7 zlzyPF$}lo|P~@2T58f{S%7#r!pIUmlZ zVLE{q!TF7F(aUA@-{1aEPS^=BReDPdEF9>M+{fHdFGAcO-y6Ih%7$N%6=}tJx~c&k z@!cf=0Prgu3cF&ZZMJd8>jaJl=AK*mPWuWhaz6c6Iur-*ja_tWE(mu&{6wns=<}+( znQ46Fx?W#}p(X(UIJgg%?sCYzZm_WzlY`!T{x|(Xo!wxuX%=0zNh?OQ*x#3{mcXy%D8IkJot|GqsVYOaf#XsBh9rVNSf z$RHwsz{>p%h=_d;yP{f8W~{VgE_thp@scmg%o1-jsAfosr=e5VU>$X;?66?d3@3df zl6G!HIk9Sn%F^3R(+~*$=Cm93ZS84%w0{bfr!}axi+;uV0o0F;E=d50nvsTWG6lKoFZ~W(G@OFtFdCbJu<; z@Cr;(tlf?i#h*=~50@A>B_}66L&YHfI$7v8e{h|-^FWF_ zmX(SA0p2iU=LEM{=ziEeQ_3{@?fM?EFrpYq95X9D2ZcdJZV?S#ag``7l1$Qqq&+wF z&+N0_un6+?Q#Jy(!IG03yG{1LW&&Z7(^i)GC)aD681}Y8LdV`%MZFXDMl@AD)UrZ2 zD#FXEsyD^DM3%^&yFfs;1CpZ2YFemL$0p{-vV9~qa{HMIE0!#7v2QYdi5YU2)JnE^ zSK%Y4k5&_x6Oc1pRgf*S|IN8`AG2It#*KC+Rhd$F zc80_yS~y?i)G#Mr97dbJ)6eY#`y9M1qL_kT+HxCT*cG?}9;(Exr5J}>bljr=T^&=KZ_WV#+IoEPvwon&n@XcsIWzQ8x$>*P@{M()!wpC@X^83q# z_g>q{PAv%=+{-d7aBIvnK%un!`jgVlhz4N%008&l9(0jE=k59~;S@_8c4QR*;HoJS zRGTWFdk{I=0}zYC0QAl?@^g2NmoknfKYw`pQrIa{)?+YbtAtYX-r>I9)~9)u7<;xY zd+YFu3H;dSiD3ie_X;TgW;*73nD=ox%-cho^_u95j1K_zI640!YkH-=-{lmUhpFZV zWTNRZKQZqo$b%u~3a+xbOxqp#BIq-mX6t5~j;H0~Q_B6R9U$(d8UQc#->NPCUheg` z#iA%4ZDR%!kE`^=9e%4+C^gT3{n{C(>c^l5loX%3qfPzdd!g}w=;;A?jfK5>_g*SRD+dB7pdd;Ye(MmdAYn`ox<_xE(8V2JpUZ zoDwUp2I7}|4wxU(ab&f|G7jiRau&7d6V0tDMi_G(L!cLJ;HvLdYU;ATs{Y319|_L3 z00yt}@!2{_=#r7&3U1WdSK&3{bDW#nFj6O~xf^U#5<2Nu;ZvQ#g*`rMn& zf+Kt1t@e(U%gWKFBn<3~Q!XRBn}=;)_0InqcR-51=c-}+J*)eTc9}P>X-Afjh~jrW ztNK@b^WbSl3G|%beo8OP;Zz;Iy;3DrP@$Gq<_mucQ5p%PS6o|LjSo==4|_515&2yE z9v{g63e>oFQWbY$rJ{1Th}^EK^siUH$8qp2inBI93#t{0CNHHd(@HiEB=J2p5TBMJ zj6<7E1{2GrWyC+p^(_Tb#%gXMLmZ|z>xA2q&k$Yq!HM7wE zUEfO%$HY2xCd(qgjk}ffHCKz&%yD{5y{+eX-Y82GJH%*3A08UvlN}K7IsTt z&DZy?{$?vzt{ttWWLtwZn#+UWYX9SyQRGJiOMW8^qBq>J{V@s%lRxe&+*>qe(neq+ zWX&Miqn1WpLCnyyXpm|_o!MmTT6{?@ZIwe01ue61P^!EVlzY=Y-PH5+WUB|jo-&z2 zq+oB6UbWyLqLjf>b9t6nnH$JO)(-17Zr){&AJc~CkKhLWh9S-aEb60WJsZI2-1^E5 z&hUYQ;M*>#OqHe052X}jLhk4jowO2y5m6ltFopj#(iY;6pZ@DP&fun6&-C^ zj1o2ih@uxIVLCCgHvlnUFJn5}=+4ElZTGR^B>5_@hT(mhP0rY?Bs~VJqrKxmzWZ6- zxVMv`5OQ?SWTqhSop)R`%>X8u{POo-%v5|QMc3#&bsuU;i7=AgBT(}cTYjmy`wsKQ$aD;I4?bv};h|q(9r4V|tt~v->^^>0 zs=<)1Wz36XwPy>6H!=op@pCCQ0$K2(3Jx0T#2kaO{mM&id)`d_iurM%ypg^Uz;6S6 z#*JPV87Tzqq4KESp{f6ey|<2v>)HMVAtVqyxVyVc;{BmLgPTnyY*_^danh(!0$Q3?szforWY(=Q1Y1pMr z9_FpHjK~v+!=o!2>k2K>klMB{!||aeJ+MW}9W<4G)ynL1C9zN*PsVq&D7s6?ek0!_ z&_9RO!DMc#?Qhh2Vhdt(*5I^~jBz4@MGR@^J=YtjX*JG?kod*h?cnfLW&8ZV891%z z(&BGv5Q;d(hH^@M`hh&W6+VU|o!i%jC{EFo&B0Z-pf)B?2|uB>Qu%0{JWiBc`DpC! z$4-1dW(!X{Z<%~+r>h2vVkzSn6*{5uX=prxKSHV(l{5V$`l14rwI`#Zy=`4(#-(=i z+qdy;3UodOp6~N7Tl58zV;)2lX;_Qt?93cET(P@+hS)re=0_q%A^4V$Z2?<`vX0As zvq}oA$i>NW@+q{jxSQw^dUtb~P>YQkW8rK&O)ONuxmU=CpBd?GPZLh7fyMz~a~!6n zS!BNoa9BrN7;6h;SXbV_Gvo95YT3kOK9`igik!mmyrDo3K~;q(kg)ObIq+GQ*|RTZ ze6^bgdlf0|g$}wiU9eXPRJ1TGw&xM{N;xVJ2bjk&+X{5{JWU;&*!1&YB4jkDLKHUf zKT*j62HkEwpKr-MS)%I;&OiiMwuk@@v29T|FN5>YID^&YnJXSCB@$WEaVbs<%WQ+q zeH~0{$Rt&b?bI;mzPF_R-nErRy2dj{H`|g0)^0z~D^L2@KL(;lvRG%C3_R>`+^qzH4J&?=1fQ@X!!D9f ztQ*zt2YpNT4p!RQj!7@+_+L56%h^-5HV!r(aXB78HXB&f$q_tP&`6WYiA3-1U8?l% z(h2>65qgL8FQUm`g?O}Em9n9b`#&&ce_;Bv=RZ&2Lm_&y?-X7{n!49vWp;%P;1p_L z@FmVuXLBrGwpc$60P)d9Kf}Om%$B_9Xm?v(p|_X71blc4>opaaoto3_8Pu-!S&STd zmPTtQ3DWA66#2n3eX-v4=I3YsPo}SNsm)8BIrlV|f1QD)(VNo3eK}svt;g4VZdqYq z(7SANW~1~2HyJhu_>a0~AIFtX7qU$=VUy3A6qX2wE?RAl1ZOEe{||ytY+aVyZ)Ig) zE?L=Hc|IWZ{54J}0sBE1>P9RTKz9GRCI}7dGgxHq4@}jfLD%#h)uo8zn>fMs*T^T) zpRmTS#B{KI7c&kDxlQ3EgGr}>I%S_;=a5qwv8)95-r={LO} zCwT%RzsV<)1w(xEesr;lQ^;gkS@O=+prv_%ZK|3z6Hoz5k+I*aZQBVAd27m!PoA!+ zr_pJ>do~G7gQDmYDUB&QYX?_oZRw=%nnN-tyaAPmDTP%sP4#tD(qyz};oNJUA>Jbc z{ru!(n9cmN-qUQ<{jxO+TXLnS)pr$jo!%+zp{e6}Xm}}7;ZmEL|58uY z$ee2@z_HM+4Mo<_C`*o1J2-Rp3Mo9)xL8Mx*(kZ|<2pGN`SueG=}u$TGW;?yyRITK zzCZq`To~Yci{i{%5v%8_s9kEJC7RfKvaFewYLT!vRgFv*E8Rh_P}!l!zCW8;Vs{hm z654y!tm0U=x{=brSLPue>1{Xr^&!ti7!ZZI^)Vz?Rhc@&hz9I3Y_c^ibOW1oDE*cm z+*ap);p8Hn-1}y-wE!A5(g={BWl}bW)MLF@PkTn8zE&LG9_X87qr4t|6#~4#@aC9B z4CcL-@66kL5j^dP#r=(%-deoyWQ4xYSN4k`@hL~;1lv_5YiBEm|0hDoa85KkBWy?G zK30_MDd6!Me-7^Ih#X5r%x#>@Ew`j&>FI3ZR~V~epv@x}XVl0v6SF@MFAR#Z(;F@X zkAtI(PhhP*m>x^84SpgC^R|j3Fv}*8&1@Mbu07n|Wp>@ob{Qvbx!jKvMG$?Dx|?qD zVPrNRPsWrZ>|HV4B*i!Z)wz;=#@x(QvAf@X#t^Bf_}{uJ^FfGOE?Z8&`f?_T9UZ^YG&tRK#&|zCqTr zazx17!3=~km5lEd?-oKK9E>7TCWVjtm>nVG%s-j?k+&vXhDDhtF0UTA_TZTXb>$fD zqRbQH)yzLgJ-T4ZQEvxJVXA2pN4*)N@8B4IGtP$MJy3VKOUyA$%F!Rozx$cU)Ye_Q zGvef;Lm?_^ApeBOW2M6N6S*9j&BsPg13O--H3;EWGAUBaF!igaeBBM_$j)=12-g+N zz>RIw8(?E5liyj=wtAs!RNq48va$iwG|Gg;K^}4@tHj1`1joa=qoiU^ZZLLdU+R$lLx2~tUYSdJAG7lPkV_tKFldrIgxow~6o#XvF`#`RAXbZ7xsp^7ad zce2Y?YYntgeBS1i>}o9ZnKaW3iYfGxV2}&V)$uN;EDf*CXb*Skx%kpDT93O=1py6O zxUis|N4|w!f-+`R_v#_NDb4_$i@iR!m z2N+&tl2ASL7Wl?5o%lkbFPh*76qlo>S+PV!>76N?E7!EPNBvC>q~r4%5M{PS+3Lk8BhCowPn zrul+9G|f|GpVH*PwB5le%Gp%pGsL5ZZKg*IO(|AiR>pUngC6xg_bLOs_ z7{)>#CwA*cdDtBchvB6+j~?2kA(eC#86~4oO9L0ANGGAv(DR}A;}iQ>V-@$zarED~ zuur%i5_?|;u^sJ+cK!2K^cVeqZ)oRLK-80n(WdEL&{}G*^lswfk?jMm1HQXzAz?TFg;i@*Pr4TG_W z+EBP{9_pMv?(W8eK@^8>v8BV7KOq!U_dP9Ld-GyUe~t5R1QbuK{RQV&(8inIXMW{3 z&c9V&+YMcRi!R0WJhmVS{q%Zbg1Y^8ybF171|zVy|8xAp+{S%SURT`j9$nFJ#2#;i zI-+Bh_CGMjJcm1!@g(gh(~#8fkl^|RN3h>-!-g{|Cd$KuDAqH3bM<392-7=DmeRyw zKuxZ)Qd1_^>;{r**odK875|aI9`hNf4wr#eHUpnjtWHq-iGD-tv~s&baFetA^qzJ@ zr(t5{;z|B2SJzoHS$+?h)7itxHl9Huk{c}&o6fblR@K#=&g#p6|9AatlkCD7gqQ^# zoJ#$Ed+V3p57lQ81GFsm10w)xqp8@zf-3vyCKMOOka=)ctm#T9cW^}_4&(KY7L4!; ztzlvYb{^H_2^M6sMpg-~=0|TG(-8srmXn0Vub@VV%QKhti zNIRjLm1S0CB_E)|MXT4xG=YGQfh!IccnxQDB530gXQe)}XlFDuLxdDYp4&*>XhGtn zvoI#xRt#Dop30ESdq)vOlOR#ca)H1bzn`IjhZ|~XQzQ6x*cC4#qEJ>F84LLHrs^I-9@8UhHT6o^`~rdh3Kdva13oP3F=Jvu9$i*`tm`Z zP3UnLO~Q2y&z%R(1|q;`Cp~93_zxSn|CFbv6=9SnzMIA(3o-@K|la9~8;%4wf`FuSD;pAz_v)viFcu1m)ql{L8-)D zWG`F1^go1mW@1Vf&I5{0lxi2pS^EMOWGN9M1F2kkan~xp5!qR4r@OG>xS)3|sqH76 z)v0D36f9WKk`bt+IX7J6_92&B4jw(U{HTulXFKgqlAmTF5~h;LcD;qi_6T#7-QW;cU{^yGCs}|IkuntJ*G?eO-6P@JHd*W z>utVP?Vd}6AMCTQxU*3?Slv&UC>|MiQEW15HY0zd1OAGro%L&$JOmgu3l zfpY28hz8$k%bOrW{>1S>kR%$eY1kha^sjRw%m8OAvJ9CXQnm9BCn4e|B!xjIj$DPsG#JmH)!3S^;q^%TcB%cLfCK8(-fwx*np1D8|oCZost`Xa0|TFih7z=876 zx2;CUk{A6(jpC-_<%ROeIRoKU#}9n-^E4X7>9t*xXLg9|g}v0w3wOAD?^F=Ctv5N8 z+q8H&0QUCG$C@`MlV0jm^LzL&wqc!$p5qorS+<{Z4Z5_v^i)@l+i|bGh4$XcZh-Q= zf<&WR_VS3#6rm7M*?5eX`A#-=x}A;|1@fNWVU3+b$C%%bLOV-76e&kwGH$fU#t^D2 zC?`{DE*k2|yN^PluinOkzaS5QVHOU5VLpx zmwn)`!vCi>zwvtx_okKCQf$k*9uM(6v~!AbhY0&FU4tFVSp1|gt27_`lKPLOrWX7x zYM?W~z`*`)tW@eS!5^6Y5tr5PKb`Q@(Q*RhkI!?LN+8m!nteN-r-C~zS!Ih1ek}H5 z^@&pH9Euakun6>Qd{&Yj^jugR+noN;IaD>|0v{4cSMi_~@vg}#yWX- zi2vr0u6N;=gY!UFAk?DID3{dP%!4U9Zu}&#gAz>GzA2B{o0=!ksZ3|@CC6IO zT8(rFS#O}JI@OONrJTnd^EDeXmj~3?txAu*B=7u8(Nw{JPg7hLf_&6ewUk?t0q1Bo z=h9RXq?aUpXF;@HBkuUAT+by0x>U-DPf&TH=rr%XFoGh!I!*-bi}83&qbLqDyA2aV zw=n04j)z@lw|?e8l-hAp@5AVp@mPPrcIv@?-0=fHYzFJYoyWz5`f-k|JMgZ!r38Qy ze0~S}Az4|W>RgsdA|&`!ReVz+@w*CYJRX^BL73;VnSU$8Eas22&RGB<)Daj=Qg9w~ z)Y{R>;jo@q&%488G?2j(Sw?~$#Wln%9nt@aj*1CkvwJBsaBoyR#stNU4`Yq`Y8{Yz z;$ww*KQ3<(Z=V)Jv4EN!msVE;<_?63TJrGDY7b^=Y@{^!8*uwWmqCTC9ZRfgjT+7T z^CDVUNl*_l1yC`oVwmE})=`&GOeloY5)#>;MA*IJ16jGvt7uXZ72D{>B=2&4lkM}Ud7eVwiQR#v@&xi4a;Ya_WBOM{yIf&KrY+0e|28h ztJA=gKkhQ2mHITY$44EI{_={Qd~SoGQ#|w~Gn%1td*)jM(^{s&H06&w~9nL=e6ywg0bZ*j?n4z4=k1ZX}y{?q&L^!6&Cf>V| zAeWvq_Wq$axKEvTY3E9{&3vrmN9NgOX*yo>Co;C(noVy@HQ}8;O&?q!7r1%eA@EGv znOy5CB;_h34*2hGMo2UX@OWw6AvmKq#43Ie#@?G0YFwX7@eME3Mxi zBO13ja{YBew)>N4#;fk1nse_ zwgz={>gS*Zci9Tvg6VhG<;c5XZ{K~l@GNOPuu)$FBjPkN#<*9 zBmEi;(;IqA!{#j8ZS|nYQ*L$7C!Z+<@{c#DnVO)Zwl#%)Q5v>~zCmMVM>{7?^@2l8 zx5r5Yyf{WV-H#SA7jg;VMs3b_-sg>l0&fSYSRaX_nf?tOF`ZnQvZ8Ub!p}K(&9_7$?RK-!{BA zH(~F&$U062sJ@lZ*-1-~=5lBeEcFS+6I ze9Aj4A(sJJ;PPZJP~URO^E2_3wn9p?xWWS~aNiNQn@rRh;;% zg2AlEp5C{3MTOH$9_i$#w;==@+{h8k^2YjCxRM#>oM>wXq zOqe;2`jmAQ9TK@BYuX6zZ8b{F7g<_~)f8xWl_uq&_Kiy?0=I~v0-Kso z`>}6bVu>Hq*@uic>2lx^Ou{U&hD@B-kp+MpZqRE9`nk^)On;Hg^i6vnmX~)pY9I~xr~<# z>8f(EfW~~iR#Vt@l!i59Mv4xrS6X?Dz5^#-omz*shIbjODi@s-_5P6g+4E*Z=$56> zRt{C%ADFN4gjApZ5?KC)_;)qupB;EB&URk1mLMY@V-I=}9#>9G(zP@a#wGHVj8X6( zT|F2yu}CKp=dzW(!Qun^>6Uoe1l!r(7MIs#@Xv{!IZqAt2zNZDe_)DXM5)({Un_}V zUt)B-k9*KZ`JfTo<s}WUA8~fKgibG)b?@EW zho91-os|TBaUz91FtGQN8c%9=POJNZ58)(uE1v*1>%Man8pfQRDWiUHy6cN9!U(>S z>~PFR*Pkl8crJpJULIOtUhtNBk1yD z8&mXzv`mv?<)%3i`?+GL*|%g~n9-WL7kIm0NucVcUt$RVT^TLnSMr6x$zbEUtOtT2 z>-ZQRRas@bLl+&uOk7qo0a9fjNtl&Pm@c?VEQefglTMlheozfncWb$P5DG-O*@92n z`j6?PxoQ1fn!WjJGWXWkql9aL9~(5)Op5{dvbaPo;HiEm{mn+@y~YQS@ziBHY$)3} z6{l)+T@pdu4J|=J5Inh1!s@#PwgBIBqu&zntEVD-!m-*~)FME1*0G*u-~&Vs`4j%BNg*?n?yIv2xK_9il`*oy!m3K4Qu|Q-7{9`bOmA>J+Ci zX}WSAtWW0B?BJAA9aBy@At|k^{18*KFzMJohATM){>aA4g>J#3lVWK=udA1w{(z#9 z4I?n{S4aLQ#J_tE{Co3%$0c|0j_v5tFCc0e>Q4LV9Z{^xr`+P8!kF+4j5_(jKg>$0PqCmNy5+98#Lfz@ihUJO>L~hK&!OX? zhM#0=hhv}Y7>hE)Y85JZ-Is-w_E!`qqFWV+%Nm?)y|&Gj&V%O?XY&0#$m?7UGc-7W z*_OrALWmPJGo}ScezRO!71m9Ikc>DZNaBf)%bfRm4ES|!lF(U5?v0N)#a%%j~GK`QXH*!PFCurPZ|LApYCgK;T;nB${TsCln$JC zrr1)o8DHwoEA~#YE-i(?8tzRpPOB_R;rKcPrh1_ticQD3n7Jl7mQ42caEoW{DvHW) zT`NcnehoR+AUnb+c#ZA9<5T~O3IxA_P-i`zC4J~Y5R#JADY#|zhU#5IUlq?Vp*dTF-_~CpH05~B@ zI}e*{;ac?aC5azg^~V58Z__Qj!tlXHNJz?!m;;e=svBoN4w8xbnd*3=f4)~0_^ulr zpVPG&STP?iTl_NKR_R?)Hj*@TLeI>aUv3KAS}1d3QQ*+KsAx|2oww67Nk~be*~;C& z`p6efd8G<&l3l&-h*tRQp7@%neeLga7rqx>`Ui$~^Iie6`Ui$o^m#MxHuUP}5*b47 zYur=qLx9z<&G1pE7}KM7FBe#$NvDJA3G|ppb^HN{ra5PpKPWPvqF<->0N3h~S%2SZ zcDI?KOX@ebzSOKD<5XUm)!%lw7+K(|6Ri`j;|=F*JGe)g@D;JAk&%_*apmzXthdQK z*Vt3>EH8K`@iBz5O ziqAxJw5*M)+G2|wV*=^3j%9ve(BVj%z!Ad~PTLPz`SOFiZ};k^vnYlMv_0~0xW zqOI(_;aQ-zup4 zTlz&|`Mqmky)uwHv4hk)4XT>jKID4qaYEE zYPP|M2rWh@Z);5mg_U|^!giy{_C9Ncur&``L1^IAsXL=r_F7^Rqc)`N5e|SpTPn9 z|2d&lah{D2)_di!bhLCMg8rl?|Jn!TWti z<3yWE8d2PmSvWMZ$98r^;OKnQ4q{VvE*BkjT_^R>>Pm*kN%cdNjI%x?%v`!J`svh6 z&u_taHVTG3MuS@H3Vk3j2M|Xpom7rKWHLWI%}wK!bDH$;E4K9?D|U7-8EgFSEB3?R zR_x~yvsUEOhE1o(r;caT^3T>TD~Go@*28~bmgBA>YK%n%7qD8UHJ>!P(&fKL3)rVD zS}r6&ECAR+3b?V0)+B00Mng?2IVBDOJG`}5tw7}o7j68r6xrJ9IW})4A}enbZ&yT# zg!?KFdd-W9hm>aPT`Dw7t1;YBV&12u_98m)%^?a-7w;LdFyIQr1eT_t_6a1g^4VwdO((*Eu&)&R!NHwZPg zpH!|SY&b~CWh^+rmj@~K?6^&%Tu!{L_wcz$UNi#VfkX7@+ooA6meO37o0g3IKN(Ix0!O_|bGcmcu$9fg3!>Ji?z6KEkf`F${2qKuX;wvRR>rckh~l*aa3J^#A9fVxkFwS{YMwF0#%zQzL~R5K@LL7HbJUzzNjem z*(g8P$}#yTWpn*d;?!HIBw=nlnjjMQG&6LI2&Mh~l_p9&T=N8a;$&>p($Q^ZKF~tD zJT-976#HUOdQKV}Y|6J!0E$|( z$w-4-)mEYb&Zd}VfgWhuLl-P>2`-we6oD_2`*gNqu^)qp5@5MsP@1NciqxZ)%99V8@j09MSE;zxP6iU|snpX- zgZb$}JOjHdtV1MZ(zRtqg_BvMLABV82s*xrcM+$&)2k4ckV5oY7Di+3Uw)oLk{`Ic zNAp-Q>@~opC|!A9eJ51wrfYN}>{->K3ovAn4iXbZi)b3p(?_`{0Qzg4a|l6!uNz%`P{M2|%5`T)$T+BaK%zJwchB6d@bJRwYKCUZZexum z0q2nO?H>Sy)sWe;{FB@^F8Z(mGIHL~yp#H~#|m54>pY=z(EvkQ^Jq!q6SDYg{d(en zAf`<+GDG}%8oyB?_aaY7OtH}Qc#_3P-#Ns#lz# zC;BT~d)G0M`y=XKQHZ+{f-KaLEK-LT2IDBSMKfc0~1%QZbNbT(B z4dAddfIRezZN!r_1c;Cm8p`Qr`iRTiF_~A{Zg5%XqG=Nq<% z4*wuPw=`IpdZAPgof5mHR(`JK>eIa*TEK?o?wm2iv(C=sTG24}U`Tdn4O!BfEy;Y* z&RsqbTL!LhxX1jIvcmR_XLa*^C>Js{+G?ft7#-}VUu`hFfrqrPw{g`D_ZO2S{*GaC zpSI_gR=mp{VEaFuOifAhK}Ml|%qYOxb1|fGrnkJ};7t zKF&|naoFYUJJpF=vvM&(OhyTi71xqDbsNZ8Axi0{T0fUf z&CcxvQhJP{dlD13GaCReeL}n04D)7RzHp41xh$K?J?yv`1e&xwTupDGYPEf})!C0( z<-rbaYVwnJ(zmv3uxD|LFnFk3K(N z?Z$uJ*ys#<&hQ;L`ZcLOtvJZn@R(QxKSlOjxXQ6+oaRGd`g`&%eF{7 z5U234y-*(il0<(S2>f~}Z$3;)0y%@sV=+wSEQBqH@>!-RH81oQ3~Lg_3SAZzNzx{`YFB z|5bJBhFouV>IC|ovz8Oca7l9U(>z!%n>$FuUVD##DO&u_V@Q}*TZJIa&UAF$f?7nV zP+TfC=WSEAIba5He2Tqbua(it>A_8XFB1mi<%gdXcX!Fb2@@n>w6W4Row%<$K~O#F z#Kt1vgIAh9#{IYJoNr}1OIk0HYZb%%9Eo^ML-Q>3Yk17D*K%mYnF;iGp=o1&JOa|a zyXuJGIV-*tx(^X`t)n!U5hRn9%(X{*0SOpWtUu%8SOr?LI!g;nTLnU=(2$cb6EHub zsTZ6fvO(qi;bPJ+cOy`=@8*7+myHuxEu-Y~&fa9AcQ--{C1z4T@(6hQ=_Rvk{eTt2 zG6lNhH_~H!tA^eD(q|AMHep!RLF(_?-Jh>Mc^P_8g?Yqh%NwcVOGvd!%_NpYX3LLb zrf_i$##G?B1$W$RvG-9Y?`}))#)WfUUHf|n1`=3ew=O3)dd&Zna1fT?-tNTdmoUVP zIKJ|q+$u)MD4LZaQ=jCJ1`c1F1zXdW6sk2(f6N)N7{;%zo`x&<_Pq7|{P#&#;X0RW zk>>e&TG@N4eO@4Xnxq14;x+WfNib<6Flj-VpSZ(Bc zYjW;3K#-Soav-;(>}**-r!zX^zd&kkuK!Cg!X>&n-)xT8YmeRZfE2kX?8eNmCkVRx z+dnYTci+&m?<50o?n5@-bUpK{!lfCWlc`jy-ixU`Ok00Gy*D5JfoCUtsribm)x$Nu z>f5Iz?;t&@$;BO>m#pnoSqmRX{_B8tJQaUp|El2jEfkX8?O5A2b%BZx3x(z*)re3? zEcNW<1h3PvS2LvR+B&e&9@T}y&S(3zBD_vTFTPAE<4CCx-$ZQKq%TJmWGU0w!ZDh_ zZ<L zrRlaa`@fF`t=@roomb*3u9==b7EEfo)sAsc{|qHyOmr=&3h-rqKx4V52=L52{{thc zlM~$NjO=mzPD|EFd3LX%LyDrj@UE%QgR0u}{kKVN_=*&Ts(>1^=)zrW;8wwziof31 zdd%G&784)K_#I{Hv#|*-jrn#cV%N&K<&xH9893s*Ct*fwas%B!g?2S-Ssj&|2Wo+H z@fID9{FTWlgLlbgYDAIQ!85*hk0VfH1+spqCtgJeG8h-BWKJ8`?84%RV`g!acn&yb zMfy=C@8f?n1#mSAlrJbY)f-?!Ru48|NTKDAhU0j>%xWt4L@y}}P}QnF=Un4a&lc#c zBTLtg-wMU4rAzIQEx3%vr4p9sLCRI!8^Y)|1*VnRCNn2U8;mOM^%Q{NHUe4+ZOm4) zkbH+NdlnA+rxu6o$G#NEq!7DAXUhOp>zFg@98|^}(ExB2DBhB|h#4J@7GZpRN=e&p2A+_ZUuD+tQ5q z8W6JvA^m2RM;-EIvDe}eYhc36l^qC;!4G+^>l>^x7c1j9ke=CdeXVNwllHrky)$;y z_(pU?Gs2j1X#UI!)lFrycR1?ec9N%aMcZ`fQF&+dnHzg2SXWQPG%BG0ITDoVY)-)F zjLGtApwH)D-v7Lp{zW5tT}(mFBUf3$^Nns54V7_z%FA9xt(k(7j1T1%oiAO~>r!Wa zn2(f03m|Mpigc=pgwTrHA-o8hN-wb-q31#Bt2C?38-`n_7w~Gp%QvO&u^tf%5oqS; z`b?+pTEyl@U4)5ZNcWwa!~Gov?;jW@Se<|ZJ-Bwggf@LT?&Wpe4r}YARGrwU8re$5 zJ8^Vt?*t3GMaQMO<|;4gbb(@4eIlLeU!nU68_9%;9i(q=cH_G*|IvL6@eQ0pL6??K zd|rFCRXXrJ{#iFoP>peyyzx&rR?|;J-LxuEq#JjfS|XEo_Z$&;buTjSo~BJ(eGCqd zl`p@l`;P_a|ADdM-3&6ta=NZTPgf}qIy+eXdG^uaRU7tSM`e&|O`GXp`F&A;YPG@l z-u{Qep@uDmJ%FSf*(aPnk;>S#;k~c=m_ncXay3Q8q4ew`Va{Nr8y-!+=sq6|t}5I< zUKCeBWKs*&gmtW|Z$P=#**WU*Gl0S9=oriEV=xA>MK-_BNiwOtYs)w*#~~N{F8JXU z3V?e(WQ$W6{;0&mu=fWh?Q#@9=qLKitNkmDhm%>d8D6uNMRQ}rh`V2T#@x#o@7+pr zqb;O(JA4fUPMqfZz4M&gGwpWJdV=$!&Z_bPfO;sn^G`!6-F}^Hj-^h>ZLSpOJ8klJ z-}8Rreys8h4bqHC`~#zG_6KGi;jgDeAi_j}kAr_;TI(K*C*9_mo(c?dp4`Y@k3=f% za*DNIHM>miN{?3Wr`TR5Qzrv54Sa2XfTIiw+fk<-vP|zM zyr;tncP-vxV5F9tQvW^`l=wz|S|lhjZm+Fo_%>2!LCaQGJ0(l7-8OOZ7t|II4bm>= zp4rI>k2DgYcYakhzIKv(jlJOL_!5e;NK|Rm6;ZC1D;-;DCv4kTrfVQ4p9Nr;&n|$_ zWwWQqnHO{4DB&wG*Jeh1|9&Rh3uDE!t}vJB&s$tt{w}j}B_?nK=ZdghOKLiOA-Mts zBw>_>ur+7{IPs4bB0kz1#230dsFt!yJG)Rlm%gq( zSHK-;oNWFEH2wblAu}lX^MALn2)R8g?Hpdp);V*^SrUpX7^V9JtsQ{k?BlGQhgb^vw0Dl&B#?Oc?TbyD}ht~mTAIU?G|4F96Xd#madYNNd7h&Q?Z$ZVI~ zTSG?Yf8IS@4N4PS8&69e`Z=1oGu+%6ovczH&^ci$b=Ixc8(6s*bJD@Q+Swxo{^ zV}>ewn@qhucyqxQu6*L~2#E^Os45ua`m*ZMr5Ek#)RS$rqs|K%HQK@QzO%X-)fJI6 zKd*Eukr#}UZ>I({RWlzzqOgbw#GUBGHcC1LnGbPnN{M)Pt0||YgnODa^mf%F@?j&0 zG|jgr>D*LM8K%tCMR^Io)hA=I$j|jF^x|8NJpYR0El%FPu_11+8oH<@Uj*pqXapq1 zFmqt~mlgN%3iuA}XNs3xDz@p4HH5}oSkz{IkqQ@4s8vhYF-`(@EB(Up zNzvK#aP&9s2cA_0;++YXNGz*}Fw8l6mWiY?q+Va{*7Q+wkxdErlv8oYmwo$q1v`?| z#7Z9)IQMXCKJkEz8EmUqDnV9u4y8oLttKTd<++T0PFglt74~?XlNbXRsKJ@rkD0fstLDzoVwI|@eRHagq^)@Cfj{Ciq_?^Or3 zCM@Yy^|=EI-$0{zs2Pnmv*AvJ^YqACSE6b&eW+GYP1=s62|O6!&7%cZoDm*e5$@f& z@)wyMY}GFl80>IqWl|VbsKIsYd9J%cExYL9+vxUQp0>0<&No_ZvE6j(6_kGmep^or zFO`r2qIsva$gJRr%ZWk5!MOx1XsM(Ntq5pn#Fut};jfltInzCG126x%w-0(bS)-jh zbRFB7&*?UrUUiwORT|N<{(+$^w_ASOo0}~aQ*GLAKDO1GOPF+Q&|>c@);&v*Lg%X$ z5qs&4(jN<#+nMJZ(I5~}rSC_){De7|OIEmOPIAIJW7D5aF@*{Lm0csA#no2vGMys| zqX;;~pscD4zy(56^AWrTk-?z({HJ?RVpyP;U)PpUm`Q`Qyuv$~yI-f2r^a9$tH(iS z?bklw#hBe{vwhr&+~y9=Z@)j;>Wi*v1i(5;D5V=7GC%A=5A*uqxuU?=)#JH=N;oBu zYC*iL_hT&y96X6_lX1-{u$Gk-X@!>5mYo&GWkqcwy-MP32V%AgAJ5R!HwhJ=MCJq6anw)x&<$T2xHvykCRn@70rO-iJfMDv zt|_y^$qh^d1q*0RbV6Stn5A#^J|>`w8`f8-&t3+1?Ab-gAFZ`^?8{Z0bYH4 zKDWHCXQ;0&8pY}bC>ucrh6)=^{=g)m4`n?0j2Hafy9Q#0WOOOJ@^)|cmR~Y z{kC?J<{xCk#EWRUPa4JJO@tmtV`&50Yld4@9ywF3W^qCKz|hklAch*Oxl<%~DjQy(21t;Hzj*LXx}WgHkhkgzS~W8vjqkUUCFMqtNx6Os>*MT3@rJO3;^3_G!J4g4*T7-s<3d9>%9L$#8tgvm~_iFUF zWm6}PCVs)4FO79KII)GFx9Q}3BOwx79w6+?f1&jJ;_JR~l31$cQ0h6+na{Au1n9F| z*Og_q@yija@~kN0;Me$FfPJj6Q=v(O;2zN_{nBw+tDy-ZA?qbVGJZwnS`ME7qR-_d z;bE}StGw94Z?5qxJxqzG!J?%prI^wch#y8BI;dxaO{y$XNU3u?pg=}FA*^YW(=bbrzkSnk@0ztDpzZZ&;Tg-zKPp~y}Dbr-_)xxfyGi#9Kh6t7~g0q0!s z_-TZzW+cC2siLV1nRB#gv0^c3a$6uIBMS#tt|Q(GWa(`+MDQR4!?AWE_v2@vfUO>8 zfqGk5fW=Hq+*p0ponviehpgQ)wO{n3=Fm2k4q5N~Q=5LD0Vl&46BkDF+(}-)3%!;fo+XUO zzAmY9p+b{a-&6P3A^pT(tP+=EvTU|Zo%07Kj1=be;lH*||3CNtOJ)H5&oUC!$nECt zYGLB=H_FM>770MkPX4!lE?!nvxb~PUaTmhDdBWHs?ry3Ut`bfT&Q6XNj_&0ANNkE0j#htb&G{dXh=?HlqrZQa{%pZ~ke8N| zhJl5Jg^`8+VE(MaNc^qzpSu3WVN0sChFh6I}5 z?+qLb+<*Ig2#E0SVcsIa!XZPe0Uuys-@?JbzJYuD4&^@{dI$Fw9s%~vdqfxvOe}09 z99#-KN`5sLLRt=hfVj6$@RtNoLsJX1DL!-x;9&n7Unm7_`5s!shyi5^3kL^__*eCR zyZ{~p69J0?8wVGkl7n0PV}0YF6&MsKi#Hf>7%-wR=YBO$8|mADo;3nd-xUZD)EB&C zEVY$INhuhaS20Je`Bj|?-?9<--ydWrph!6-##v6zYt43zQgj8tbRQV`{vUOnc|4R` z1IH1%B-@lima^|NW(H%Y!X(>>8cPgj#u$tk8M0<6G=xHOWhu)f#aKd%XpEiGT#O}B z8S9O%%98s&-OqjRdtdi`pU>y~@%+y3?B}2Fd7g8Qal{LQ-14=(M}IzdH7tCXV57!4 z`ehwfJCU`-)~#lVKahLrG^~_d>pd>dE07LRP% zesRH$o6R!g8eQw|5T>?xx-3l>r;LanTWv}o|0u_k9I)dy?2#sUzjV(8XDh;thR`w% zyS9D!?KDsIZ-Z>>)yQuB{*wN$TFC)v_3Dy(tBB%HH*T=4)T)x0nJRUssX`7HPl#=z zf-bgm@KSR@jU$H#;GN^#zQ!?;B{P}I!($&vJ)s**8BXpm58Ep)e`vT-Ajo>mqk_t4 zMBndd@xRTK26@811m>FJs_@v_pBsjw1Ndqe><7yzs{tyuPhMQNYmr(@4_i8rKTW)O zpeIaQ_FVo*4v&Cln%KDZk1^mYb)`Va*Bfl-Tz0jlOnhRkVj~njaTU=@zahZ}j4M6-vEwgNN0?N6X9wg`(3| z1@zueWX`8mt&{|$9E!YK@UsgZ6Y!b=5;}Gsn{ln~H^sRCwyhu_gaCNE&$s0sNAg5^F-0;B<6-*W=)zCROH{M81Air$SE$fm&mx6AeI+u4K zbW!y{LZ;%UbY760S2v?IuFj}VP6K<_JbURbalyz@Gicw)k-1@_jpVD`KC4b8QP6%Z z$U&JndY#jLm13#xvc&hxpK{^C;TRn;$N0S(F)IfLJKK7#baE!26`uosdyyAu5czH$ zQ?;!sx06?C!p*=t$NAmwY}&ZuacARa3`gMX9i36)5^&Xf`~|bkOzq7*_XoQ2R=+*e zJyIt7=-a8r219|Z8s91DH_FwB_{b#g-ez*bPw80~xcgWc*(f24Qie5VkqSPEN1vQ) zyG?Fn1b>Mkd<^hG(QI;`8fq@;AM_d2DDx~^Zq4H*X*D;+>{JG`y~-NG57~bHempjo zf&NsBzlT%H=X^|K*3vpKm5b@2`pg9h7vfEg=k-);6%&>i*dUPk*IH9Tc6mY zvQxVy8%jU=zWJ!NDE3={x6iXmX?vLFPBYHf9A|q7ndqHXJb>vqioM@iFZxz_|;Bm2?ZFu#2-`*^$)o;%58=Se^}Rq)Hv(mT%!N;*FIj%4*98j#s#!^icH!nJeFz zCGl$_JC0r0cX{w~=MeqA;l+kT+iLgG{qEjl;GKx7`>H1&%O5=0&(S@TuMt%-0& z{nr#u`JET#Twx~bHDjMyj7z?r^ctBpm1h^O%^lu}wFNfnnNE$a7vf&UtiJsAef4*) zDd(mrT-B%Hn&=8P@y^=9KDq4h@eTh&oiX&~Z{zoX50-&jv3=5Uaa#?*Deh*mCjN6` za+Y!CDue2-SM(XZpG_iJ(XG0z_lXit73c1wDjtD`AKF zrhjV9v#_#RA1i}*CIoFBj})~IIh5Rip@`mWD<5f}%2UR#`W4xNaeGDf@G2U(pH3$0I5j zWlgB~1J=qnF{L5Pfo7eErS~9-4IXsz=qYB}TfWY-AHC>L1dQC;lT$Nkzu0`^mRAVE zP6=6{4;!-^>U=Z(w?}l_R+4!aRXVfttS3~*9Ug8Q@?$X)aCd`FVu)h z;9fsU8HE~4D%0CZ@psS>R}#=MwE<&agi<2~Htb(wvs0Zhlm2Tt*AJ$PQN^-GUEK4x zyoRy_2;Cy8FTF!ICM}EBr1@=K`K#;?8UrhIGjwmVlXHne$|;Uf8L=e=UeHYk z9MtCPS}J+j#LgqneLTgquM`YJBoDsgWka4=wa1m5F5;_)%$ABp*NyS`*LGSyjhgRp zDD3vaO1%BJ-5*nX739Nl2rwFU$*>Li@{8JJ?`PojuIQS@xI)s+HrHCw0S*&*lh6SdHBPm#EKfxF}p05(igDV!@`zQ3BHBguN;Z_Dq@) zTIuw1c?}Cgaf@oYm?0i&GSNJ9ltx$VJ6_v%qMK4Pje6xW^B&Oyp7wpxIJqIzwFLXX zC{$b-y=bw~ya`pD!wa6GN4ZzeF>JeKvul&i>3WCFq0ArdB!D#*o{1RTid?^Cgl_B{okZ0*(W5E4&gPug$lE( zhK`s*uEn1FzA?BthZu+hLJarZkR<1S65AhS^arv*A&9@Q4EiUQ{fTuD*q`svznkEH ziU7L-VD)dS-v6@d|HrCxhCmL)%Buzu$XbCUJV`sm52pt39DtEz0x9?eh73fWFw_M@ zz+gSF9>f3whd=>?9m((i#yN(P{K9)jCkj|K-Ic$WqNir%e3fmHv{41kW`l@y`uj{mQ; zfQHyz7pOi|7Z}U`^1;CR|ELQLYybu(!Y-eI?(T8{A8dD>0Zngr+W^hrZ~2i!F@!*D zsNj#K0B1rJ!1VyXBZF}ydB~6H0{}jXL?X-oSTR8T7^pBZCY1cgFu{Pu0jQvo(n$v^ G!G8ga5UU3O 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} - From 919b2a52fd63fcd1210aad038454d0ae31e53175 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 26 Feb 2012 15:59:48 -0500 Subject: [PATCH 04/14] Changed accumulation of tally scores between cycles. --- src/tally.F90 | 42 ++++++++++++++++++++++-------------------- 1 file changed, 22 insertions(+), 20 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 9920b8ba21..0eb222e1f5 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1048,6 +1048,26 @@ contains 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 history + ! within a cycle to the variables val and val_sq. This will later allow us + ! to calculate a variance on the tallies + + score % val = score % val + score % val_history/n_particles + score % val_sq = score % val_sq + (score % val_history/n_particles)**2 + + ! Reset the single cycle estimate + score % val_history = 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 +1076,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,22 +1087,7 @@ 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 From ed9f8c1133caf5a9e95e86ba61473f6c13b8890e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 26 Feb 2012 16:49:43 -0500 Subject: [PATCH 05/14] Added collision and track-length estimators for keff. Only works on single processor right now. --- src/DEPENDENCIES | 2 +- src/global.F90 | 8 +++++- src/intercycle.F90 | 72 ++++++---------------------------------------- src/physics.F90 | 27 +++++++++++++---- src/tally.F90 | 50 +++++++++++++++++++++++++++++++- 5 files changed, 88 insertions(+), 71 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 121e635e27..e93239e57a 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -122,6 +122,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 @@ -206,7 +207,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/global.F90 b/src/global.F90 index bc414c84c8..4649745890 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -10,7 +10,7 @@ module global use mesh_header, only: StructuredMesh use particle_header, only: Particle use source_header, only: ExtSource - use tally_header, only: TallyObject, TallyMap + use tally_header, only: TallyObject, TallyMap, TallyScore use timing, only: Timer #ifdef MPI @@ -136,6 +136,12 @@ module global real(8) :: keff = ONE real(8) :: keff_std + ! Estimators for the effective neutron multiplication factor + + type(TallyScore) :: k_analog + type(TallyScore) :: k_tracklength + type(TallyScore) :: k_collision + ! Shannon entropy logical :: entropy_on = .false. real(8) :: entropy ! value of shannon entropy diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 3280d65a93..9276aaf7ef 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 @@ -436,18 +388,10 @@ contains 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 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 - #ifdef MPI ! Collect number bank sites onto master process call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & @@ -462,19 +406,21 @@ contains ! 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_analog % val_history = real(total_bank) * keff + k_cycle = k_analog % val_history/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(k_analog) + call accumulate_cycle_estimate(k_tracklength) + call accumulate_cycle_estimate(k_collision) ! Determine mean and standard deviation of mean - keff = k_sum/n - keff_std = sqrt((k_sum_sq/n - keff*keff)/n) + keff = k_analog % val/n + keff_std = sqrt((k_analog % val_sq/n - keff*keff)/n) ! Display output for this cycle if (current_cycle > n_inactive + 1) then diff --git a/src/physics.F90 b/src/physics.F90 index 4eb9852834..4d377f790e 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(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(k_collision, p % wgt * & + material_xs % nu_fission / material_xs % total) + end if + p % surface = NONE call collision() diff --git a/src/tally.F90 b/src/tally.F90 index 0eb222e1f5..654b3849bf 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 @@ -1092,6 +1091,55 @@ contains 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(:,:) % 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 + !=============================================================================== ! WRITE_TALLIES creates an output file and writes out the mean values of all ! tallies and their standard deviations From 3184881e01eb7abbede021b1d5b0c9a8318ecc3f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 26 Feb 2012 17:16:19 -0500 Subject: [PATCH 06/14] Changed name of attributes on TallyScore. --- src/intercycle.F90 | 8 ++--- src/tally.F90 | 84 ++++++++++++++++++++++---------------------- src/tally_header.F90 | 8 ++--- 3 files changed, 50 insertions(+), 50 deletions(-) diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 9276aaf7ef..f57072f545 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -406,8 +406,8 @@ contains ! cycle keff, we need to multiply by that keff to get the current cycle's ! value - k_analog % val_history = real(total_bank) * keff - k_cycle = k_analog % val_history/n_particles + k_analog % value = real(total_bank) * keff + k_cycle = k_analog % value/n_particles if (current_cycle > n_inactive) then ! Active cycle number @@ -419,8 +419,8 @@ contains call accumulate_cycle_estimate(k_collision) ! Determine mean and standard deviation of mean - keff = k_analog % val/n - keff_std = sqrt((k_analog % val_sq/n - keff*keff)/n) + keff = k_analog % sum/n + keff_std = sqrt((k_analog % sum_sq/n - keff*keff)/n) ! Display output for this cycle if (current_cycle > n_inactive + 1) then diff --git a/src/tally.F90 b/src/tally.F90 index 654b3849bf..3ba8d9d167 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1042,8 +1042,8 @@ 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 @@ -1055,15 +1055,15 @@ contains type(TallyScore), intent(inout) :: score - ! 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 + ! 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 % val = score % val + score % val_history/n_particles - score % val_sq = score % val_sq + (score % val_history/n_particles)**2 + 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 % val_history = ZERO + score % value = ZERO end subroutine accumulate_cycle_estimate @@ -1114,7 +1114,7 @@ contains allocate(tally_temp(n,m)) - tally_temp = t % scores(:,:) % val_history + tally_temp = t % scores(:,:) % value if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into a @@ -1122,15 +1122,15 @@ contains 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 + ! 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 val_history on other processors - t % scores(:,:) % val_history = 0 + ! Reset value on other processors + t % scores(:,:) % value = 0 end if deallocate(tally_temp) @@ -1290,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 @@ -1367,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 @@ -1383,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 @@ -1399,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 @@ -1415,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 @@ -1431,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 @@ -1447,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 @@ -1544,16 +1544,16 @@ contains 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 + val = t % scores(j,k) % sum + val2 = t % scores(j,k) % sum_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 + t % scores(j,k) % sum = mean + t % scores(j,k) % sum_sq = std end do end do 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 !=============================================================================== From b633ba42cbd350f7a10cadd1e95d385642de714d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 26 Feb 2012 17:44:32 -0500 Subject: [PATCH 07/14] Created elemental subroutine calculate_statistics. --- src/tally.F90 | 37 +++++++++++++++---------------------- 1 file changed, 15 insertions(+), 22 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 3ba8d9d167..ce4fbd511c 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1521,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 @@ -1530,33 +1544,12 @@ 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() 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) % sum - val2 = t % scores(j,k) % sum_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) % sum = mean - t % scores(j,k) % sum_sq = std - end do - end do - + call calculate_statistics(t % scores) end do end subroutine tally_statistics From fa02077a912fcb7dd81f8f25c4cc664cf3ad6855 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 26 Feb 2012 18:52:12 -0500 Subject: [PATCH 08/14] Moved estimates of keff into global_tallies. --- src/constants.F90 | 8 ++++++++ src/global.F90 | 14 ++++++++------ src/intercycle.F90 | 47 ++++++++++++++++++++++++++++------------------ src/physics.F90 | 8 ++++---- 4 files changed, 49 insertions(+), 28 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index d2fad27b82..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 diff --git a/src/global.F90 b/src/global.F90 index 4649745890..2404172f54 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -96,6 +96,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(:) @@ -136,12 +144,6 @@ module global real(8) :: keff = ONE real(8) :: keff_std - ! Estimators for the effective neutron multiplication factor - - type(TallyScore) :: k_analog - type(TallyScore) :: k_tracklength - type(TallyScore) :: k_collision - ! Shannon entropy logical :: entropy_on = .false. real(8) :: entropy ! value of shannon entropy diff --git a/src/intercycle.F90 b/src/intercycle.F90 index f57072f545..db51dcb31c 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -385,42 +385,53 @@ 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 + integer :: n ! active cycle number + real(8) :: k_cycle ! single cycle estimate of keff + real(8) :: global_temp(N_GLOBAL_TALLIES) message = "Calculate cycle keff..." call write_message(8) + ! 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_analog % value = real(total_bank) * keff - k_cycle = k_analog % value/n_particles + k_cycle = global_tallies(K_ANALOG) % value/n_particles if (current_cycle > n_inactive) then ! Active cycle number n = current_cycle - n_inactive ! Accumulate single cycle realizations of k - call accumulate_cycle_estimate(k_analog) - call accumulate_cycle_estimate(k_tracklength) - call accumulate_cycle_estimate(k_collision) + call accumulate_cycle_estimate(global_tallies) ! Determine mean and standard deviation of mean - keff = k_analog % sum/n - keff_std = sqrt((k_analog % 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/physics.F90 b/src/physics.F90 index 4d377f790e..8ed65c8c7f 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -99,8 +99,8 @@ contains call score_tracklength_tally(distance) ! Score track-length estimate of k-eff - call add_to_score(k_tracklength, p % wgt * distance * & - material_xs % nu_fission) + call add_to_score(global_tallies(K_TRACKLENGTH), & + p % wgt * distance * material_xs % nu_fission) end if if (d_collision > d_boundary) then @@ -126,8 +126,8 @@ contains ! Score collision estimate of keff if (tallies_on) then - call add_to_score(k_collision, p % wgt * & - material_xs % nu_fission / material_xs % total) + call add_to_score(global_tallies(K_COLLISION), & + p % wgt * material_xs % nu_fission / material_xs % total) end if p % surface = NONE From 6a72d02d8d7dee40ac32e26da66175befaeb5f3d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 27 Feb 2012 12:14:15 -0500 Subject: [PATCH 09/14] Added statistics and output for global tallies. --- src/finalize.F90 | 4 ++-- src/geometry.F90 | 17 +++++++++++------ src/hdf5_interface.F90 | 35 +++++++++++++++++++++++++++++++++++ src/intercycle.F90 | 2 ++ src/output.F90 | 21 +++++++++++++++------ src/tally.F90 | 4 ++++ 6 files changed, 69 insertions(+), 14 deletions(-) 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/hdf5_interface.F90 b/src/hdf5_interface.F90 index 006c6cf3c9..fafe30109c 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -101,6 +101,17 @@ contains 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 !=============================================================================== @@ -410,6 +421,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 !=============================================================================== diff --git a/src/intercycle.F90 b/src/intercycle.F90 index db51dcb31c..5ecd79fe84 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -387,7 +387,9 @@ contains 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) diff --git a/src/output.F90 b/src/output.F90 index 5da98d1ca7..0446ba5676 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -942,22 +942,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 diff --git a/src/tally.F90 b/src/tally.F90 index ce4fbd511c..af6d93e738 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1546,12 +1546,16 @@ contains integer :: i ! index in tallies array type(TallyObject), pointer :: t => null() + ! Calculate statistics for user-defined tallies do i = 1, n_tallies t => tallies(i) call calculate_statistics(t % scores) end do + ! Calculate statistics for global tallies + call calculate_statistics(global_tallies) + end subroutine tally_statistics end module tally From 114aebc487174f8428306f8f36e304abec1bc156 Mon Sep 17 00:00:00 2001 From: nhorelik Date: Tue, 28 Feb 2012 13:01:19 -0500 Subject: [PATCH 10/14] finalized ppm plotting --- src/DEPENDENCIES | 7 + src/Makefile | 2 +- src/OBJECTS | 4 +- src/geometry_header.F90 | 1 + src/global.F90 | 8 +- src/input_xml.F90 | 116 ++++++++--- src/main.F90 | 2 +- src/material_header.F90 | 2 + src/output.F90 | 46 ++++- src/plot.F90 | 287 +++++++++++++-------------- src/plot_header.F90 | 31 +++ src/ppmlib.F90 | 101 ++++++++++ src/xml-fortran/templates/plot_t.xml | 16 +- 13 files changed, 429 insertions(+), 194 deletions(-) create mode 100644 src/plot_header.F90 create mode 100644 src/ppmlib.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index e93239e57a..72f98c217d 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -68,6 +68,7 @@ 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 @@ -106,6 +107,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 @@ -155,6 +157,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 @@ -182,8 +185,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 search.o: constants.o search.o: error.o diff --git a/src/Makefile b/src/Makefile index 465845ff32..14b74bd85a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -28,7 +28,7 @@ USE_HDF5 = no ifeq ($(COMPILER),gnu) F90 = gfortran - F90FLAGS := -cpp -fbacktrace + F90FLAGS := -cpp -fbacktrace -DNO_F2008 LDFLAGS = # Debugging diff --git a/src/OBJECTS b/src/OBJECTS index f817f4cd5b..e78fcd0c62 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -28,6 +28,8 @@ output.o \ particle_header.o \ physics.o \ plot.o \ +plot_header.o \ +ppmlib.o \ random_lcg.o \ search.o \ source.o \ @@ -35,4 +37,4 @@ source_header.o \ string.o \ tally.o \ tally_header.o \ -timing.o \ No newline at end of file +timing.o 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 2404172f54..bd4b7471c4 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -9,6 +9,7 @@ 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, TallyScore use timing, only: Timer @@ -38,6 +39,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 @@ -45,6 +47,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 @@ -187,10 +190,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 @@ -235,6 +234,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/input_xml.F90 b/src/input_xml.F90 index 9d83c4b1df..e1a9222854 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -8,6 +8,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 @@ -259,6 +261,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 " // & @@ -529,6 +539,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. @@ -1086,16 +1104,17 @@ 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_PLOT_XML reads data from a plot.xml file !=============================================================================== subroutine read_plot_xml use xml_data_plot_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 plot.xml file exist? + character(MAX_LINE_LEN) :: filename ! absolute path to plot.xml + type(Plot), pointer :: pl => null() ! Check if plot.xml exists filename = trim(path_input) // "plot.xml" @@ -1112,28 +1131,79 @@ contains ! Parse plot.xml file call read_xml_file_plot_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 + end if + + 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_plot_xml diff --git a/src/main.F90 b/src/main.F90 index 4aa5633b9c..af9824a6a0 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -6,7 +6,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 0446ba5676..c4209ace13 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) diff --git a/src/plot.F90 b/src/plot.F90 index ce7b26ed6c..e23ece1c1c 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(3)/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(3) / 2.0 + else if (pl % basis == PLOT_BASIS_YZ) then + in_i = 2 + out_i = 3 + in_pixel = pl % width(2)/dble(pl % pixels(1)) + out_pixel = pl % width(3)/dble(pl % pixels(2)) + xyz(1) = pl % origin(1) + xyz(2) = pl % origin(2) - pl % width(2) / 2.0 + xyz(3) = pl % origin(3) - pl % width(3) / 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 = 0 + g = 0 + b = 0 + 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/xml-fortran/templates/plot_t.xml b/src/xml-fortran/templates/plot_t.xml index b4dc9eaea8..ef0875c47a 100644 --- a/src/xml-fortran/templates/plot_t.xml +++ b/src/xml-fortran/templates/plot_t.xml @@ -3,9 +3,17 @@ - - - - + + + + + + + + + + + + From ef90efae30a51256fb17138c88ef4941451e5762 Mon Sep 17 00:00:00 2001 From: nhorelik Date: Tue, 28 Feb 2012 13:19:41 -0500 Subject: [PATCH 11/14] minor bugfix to xz and yx plotting basis --- src/plot.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/plot.F90 b/src/plot.F90 index e23ece1c1c..97f7c4a200 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -76,18 +76,18 @@ contains in_i = 1 out_i = 3 in_pixel = pl % width(1)/dble(pl % pixels(1)) - out_pixel = pl % width(3)/dble(pl % pixels(2)) + 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(3) / 2.0 + 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(2)/dble(pl % pixels(1)) - out_pixel = pl % width(3)/dble(pl % pixels(2)) + 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(2) / 2.0 - xyz(3) = pl % origin(3) - pl % width(3) / 2.0 + 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 From 7f32a44d49fe8613bf90637e19e6282ac9de6f97 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Feb 2012 15:19:51 -0500 Subject: [PATCH 12/14] Added error message for wrong length and made slice default for plots. --- src/input_xml.F90 | 5 ++++- src/xml-fortran/templates/plot_t.xml | 4 ++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e1a9222854..b72f115e0e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1143,7 +1143,10 @@ contains pl % aspect = plot_(i) % aspect if (size(plot_(i) % pixels) == 2) then - pl % pixels = plot_(i) % pixels + pl % pixels = plot_(i) % pixels + else + message = " must be length 2 in plot " // to_str(i) + call fatal_error() end if select case (plot_(i) % color) diff --git a/src/xml-fortran/templates/plot_t.xml b/src/xml-fortran/templates/plot_t.xml index ef0875c47a..d0bd38cef4 100644 --- a/src/xml-fortran/templates/plot_t.xml +++ b/src/xml-fortran/templates/plot_t.xml @@ -5,12 +5,12 @@ - + - + From c3a4cf23a9ae333478b74dcc7a5fc44b36bbb3bc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Feb 2012 18:37:24 -0500 Subject: [PATCH 13/14] Changed plot.xml to plots.xml and renamed root element accordingly. --- src/input_xml.F90 | 24 +++++++++---------- .../templates/{plot_t.xml => plots_t.xml} | 2 +- 2 files changed, 13 insertions(+), 13 deletions(-) rename src/xml-fortran/templates/{plot_t.xml => plots_t.xml} (95%) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index b72f115e0e..fad53be3ea 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -32,7 +32,7 @@ contains call read_geometry_xml() call read_materials_xml() call read_tallies_xml() - if (plotting) call read_plot_xml() + if (plotting) call read_plots_xml() end subroutine read_input_xml @@ -1104,23 +1104,23 @@ contains end subroutine read_tallies_xml !=============================================================================== -! READ_PLOT_XML reads data from a plot.xml file +! 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 integer i - logical :: file_exists ! does plot.xml file exist? - character(MAX_LINE_LEN) :: filename ! absolute path to plot.xml + 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 @@ -1128,8 +1128,8 @@ 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) ! Allocate plots array n_plots = size(plot_) @@ -1208,7 +1208,7 @@ contains end do - end subroutine read_plot_xml + end subroutine read_plots_xml !=============================================================================== ! READ_CROSS_SECTIONS_XML reads information from a cross_sections.xml file. This diff --git a/src/xml-fortran/templates/plot_t.xml b/src/xml-fortran/templates/plots_t.xml similarity index 95% rename from src/xml-fortran/templates/plot_t.xml rename to src/xml-fortran/templates/plots_t.xml index d0bd38cef4..af8f2c6e27 100644 --- a/src/xml-fortran/templates/plot_t.xml +++ b/src/xml-fortran/templates/plots_t.xml @@ -1,7 +1,7 @@