From 8f5723ca23066c76876e4a6544465f7ef72d5024 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 13 Dec 2024 15:19:01 -0500 Subject: [PATCH 1/6] update the "Reading and plotting manually" portion of docs (#304) this adds how to exclude ghost cells. --- docs/source/manual_plot.png | Bin 20382 -> 21175 bytes docs/source/output.rst | 35 ++++++++++++++++++++++++++++------- 2 files changed, 28 insertions(+), 7 deletions(-) diff --git a/docs/source/manual_plot.png b/docs/source/manual_plot.png index 61619fc2be8cb0102da576283cd6cc3accec7307..358368709f5fc5ad32f70be36a50d178aa99ec4c 100644 GIT binary patch literal 21175 zcma&ObySpV7(F^PNDd9s0}4uaBLYKrt27J}($Xa$If#^`pnx<2BGR3bA|(w1BHdm0 z9nbpZ{o~%XjtdX-&D-Cz_p_h<8l|P7NO*_#4g>-rR92GLfk4n<5D40PTx{^0R&;j` z@I%Z)!N^0`8R6k&;bsF-xA1UrboOwxvt;tLadWqGb`szgEs300#m|i~9abkVb4{BlHZ`}guZ!ypA6DzfJ1oD@KYD;Q`Savy=jxr42UskE zI+(?r+Sa-g71qq0X^sN!h(A|1sbcpN1Kzp2A%+^W;&6z2t?!+$bgT^c&1avkT|-%w z6CTOS%jc0L&SNmKu;h$&n}4;5Mag)neSi*E#WF9 zgiTLBIoh6~x~F@29+!s;AIg@AySnmkijW0e@`<|tA-Fi*UD#R6y#qHWd-7Vp#9(Qx zNEZ^!XVty9-^19`=yN#Tc$+Td;4&pRS)j<_dIC%z%Z3wLM?{`7;p)4v~CdptMuNwlJNfH}8bAGFTXF@|S3v zuP|`5w6q{KGqFe{L8(~_L#|?+(_u3Vj%>|vZjRbt?#a#1#~zlx}k~d zd*6F;-#aPJgqw;)MJg(Xdi$+DEXY}^$6PQ|_l1k;c#qgP_L_b`lwd*p8`-z4%F%!D z{hCaZl9NAzwsxCN8BBw&U-}<>zrT&m%7r6jKgK@m+Cy*D?C;ZX>IZ&Z={WVppr+iU zfpUF)eS{zdbx?l(_z^6z+ejuPBt$>q52r_d6?B;3B%l{3FK;|hkv}!7bz$bD#Gm)w zEc4tQH9)Y+$S{tNkJnRGP*PD@oor3v>5?aHi55&RHWtUZ7ZVPZT#(vnfGPN%?#$Wo zl6qywKVtg9>GC&nVYdgY3UAQO716G3&5T9O%+ooi>L)L=f63EE#Q&c=)1C>xxCY1S z9k)re&*2}7)?l>4va+t}Y6oo{oshu{5ex_=*ORco>x++tg}9KtFK|%3%Bo`-NEiL4|gcn@T>@nEaBIaWAl~|Gazk7#lY;1$lrj?dh z2Bk*OnO?At>}-2urYtNhKPf_?|L2blw`wQot4GhaUqto>H#S2|Z+q}E>&8BuY8A%t zQGeuq`89_UOJQd3hr}me;@6R+0}IU$dV&HysBdlO0-GWpn;_e8UluK-(8;~;iMX~g zJYu*5WlL2_xeKXi*}BtQR66b_I`*0xsg)w;UWCRNXiOCL&OS1-eY`%PuVrD*&tS4r zYtoYUvB{rU`uj`7WM<38>A_lP2QJ>0)0I> zW+6R2_=Hd)*Vi^`C7X@xV!pDI`&yG0prPNx4Vt;yM^+wLWFaJtTvsL}(0G6Kzf*(Q z5O3SzGz|t>c$;~c;2Pz%QPm+iWbwX01WN~{1dGP7mCGN`kNCz zdSevS$hUuPY!{556`$q=#R(2DX2hr|@=_Na@(MhbS6>8nXCaYF0n-4ytil1sO1rpCemly&#@=Y7A-B1}O?4u2>p9@ziR;fH|$+6&gOz(_xFE8l|G(gQW1?Qx__T6m=O3+_>)I<>?Hpoe_IiGD!>cqBV} z<5HU1y!vSGoM1R=u_ywel@Z*WX{9vGIh5(1qk=Z@6%k$rY}mg~*?2~9(anj=ZV}S;Q-p%7-Bu?4D3nd_i=cp+=WP=2xcP7RiLGqUudqNeUI`Fkj`5 z-RGlbR{u@~t@GCJ<}Q^dfAiE&_b%8Fz1&~o0wEP02lg~b5@_%f&AaqQQ%9UO;54aM zJ9^)951+rTz<&t+_}+%LZen2Yu~&x`@8fiPAMRLAIph2#g>d%qbM)Z9#7>9K?rGs^ z(EmObudF%Oa>aV)3eKxc?Nt1<~qk_Z8Uvh0q zkL)+(9(R-f2ga5+SLY~%|6ITEK=S%%w70jnb!V<&rtvn0r#TV;C=Id96)!xN_VsFi ziq*l&cklDnlyO=aHqlSP)i#JuhGH5EL>T*Oog7x_kbw|Pok!czAK~~8cW3k-fN`&d z0IK8TV|I4-=ZDe_p4*QtEG(*K{tUg3A-pT>aaJwGKrG1+cTm~CiVlg)!PO?WqIa8J z*)3b6PEZdO5j8ZBwHJ1_`)Hf>thG-}pDMd#V>#1aj$DCFAi+sUD8o!PfeAO56|I_( zAuny_R(y;#KC2}D4fBTiO6w~O&};C_HXIJeo+vkq8_89y@mYwDzO%BjBG|xn|2_`D z9{RwzoE#R;r{7A)x1WA0pJNzEj`1zw<*$?-$1zo<{9a_#TjvCfK1_NV`W-(xx zz5}WC?86+${#7%zmVNbmcTYHyJ=Ikze9ii=VBboD(y+ddO|F184oo!{7$yoO3Isy2Q}sJ3Yx0yNb2_h47PVT5NNA7Y_U_vYwjFZvNgv4f6*8IR5kVQ z?*Qkq+McdH8!r!hZFX~zf_jOB9e)K^TJ;p)CTcy_nQp5|ronq@@!|-tOj=RCNNxxz zN?)dXAz3FPBu{~C-0!|iAK09impA?bpNhM){WX3rBMqCBazZ;ltjrKopiWIL=|DaR zV##zCo;NsvJ6xcQCz%m0xO|Mmd5?R5OLJQx=G~t4Gk`7v zPR^sW&QJRV2aqK3{zi!UL4ugNV^tkpL>NFUkhciJct%py(Mw-EvtBU=zmqRg-D%I9 zf*W8$wn}M~@lb8hQq8pP|PdHu~oy^WAU2W}GaZ$`R)(5w)f65E_4m~z4a zJg>@zHf;L`Y#Fe}^!2SVQ^9^D=vzf-sZF5)tMpxd3J3Q%HHurQ*%q9HY_`dBa9!;+ z?0tY`b~%QO(@Cl&zTW8uDXk#P?70zOJ1^oFDpIlVDBW@T>nL{$|1kv`F{~L2L?^(* z96dH$`I#aU{a>phT%}wBoM0i+p7WzfmWtFt##@yzn{iOfK^SOz0(}3iut1y|6w^995oOqu9AAeiRX`3<89UE zIZ5Nx76{f@t|h`aFn^YW8#HqOKc9rKVDw%F)m2!*ODIATw)@7eJP6W0ly5Nz&TqT~ zkC#DfU&jFw^`*ywI95XV^rytAZ`yFA#p61_f1$V5g7bNgWxYAN!k5(4I#Ne-1w|JB zO*l$v=o1?1_Q?oXJW6Fh9TA@LyoJqtJooaRjPD=I1z(n9R!lY{+PVWD1uJp+;3AX+%7~-_LMyY}tgo)O3lfqXL<(g_{$RtzVV7b02b= zUd7T1-SL`r#d*hLHvcCe%l|<6KL&z=Gr&MNJ3eurtU0Uw(tXkX0^_>2$Qc|Os%HLW z>$2e6-3K!Fjx5WeIePoM4TF2i=vk>-a+-dZrzQ57S=`mF7Vq}VMpNmP?-3Fb)*cQ? z0Lp<;Fgmq!ZP|z9H2#9g7c`8GW^mf2;(oLA9V1Ey{nnx6rss)adTB?BKdPA@qe%oELZ~T(wfYz7>cDO zB#;s?NQ!Q}I=>F26LoFGfrK$3X-?j833Ng;GKrCvLfH&Ft9)RFwqaIazEny5fCyF1`?j_O)dL0OJl zZ4Q%kFTPo%`@te4=BV`&qIUK)@X`{Ieg-q#xQaZ~XYC7EP1M4%A6>lqb@LKqHfntULK|HM%D-NIH8-n-&7#J~vRZ6}Dg zO9na>@tmWW0#2wXP9Ljn4TrQ^Sg3l#P^nU6K1Kd>-KI5?@B^jY16?WQY13jt_v-r^@RF#5cKMM4z z(QG9DMIHk+*u3W{b}Ly52O0GIlpqL;$I&6>y;GB~x+6hCSOnX!p+&$6P~-f}AwUY! zQd9<+D<~DS*ozxhjR!XR-#AB(sAY{t1|!AQ@&0V4Rl~v{`TD0B_yG z)$WV@h6ro0L;jofU7Nv3W89rav`v6|@KPa|Wc6e!fS`@Wo5e1cp97Ihf+lZM!EzD0 zC3Yyyw25HB29W<2)8uQvpt+YqPSWxpe!|xP=!;`K9eCPl8TTtBbH)}=6l&D0YWS)7 z_)*gVm|7^Hj#V1JQ@`v)i`v5CL>|f!Od!0Rd%Iu*hyp~7oenIS$+I*l4=Q*+3S5Tz zNO)uaV~TqB%{aFz$B;lDvi8oAfs8OpG9{`HM`;mn!gzY%{`Tc!-m!a;V09zG zN^$3>$7}Jn=iAD6-)~tb31b%*5)iO_j@r)=Pjf5*aC*AviwRY$(wT7|%g0%ul~e-I z9`~^{`VX!nJ?|5?-!YncHU4iwEb-whb0|>+`I)cpz2NhAG)qG5k*KBm&0!8wc$VeF z1&`kA*DP)Mw%vf->J^Yzz^ke!yVF6CTzX@b7A^$q1}9tZeu4J3whi1{m{7Yo@)JUS z=qomBYtL5X5@1N7^yXOSuIpoAr?xowLS8U?N7{vmTRmA;H8fYZoRZvfh?=6Qe z^*28&<^#%0UVXkaEWTxPX~HyT45~U92^bW|SVEXWX?F+cjcP&G2Q_H!cXNUe^BT z9%vu;*V(`i)PpPLtm)FQr8KB=g5ce|&~=`Z^!OkAgaU=(|N0%CFm;hi?cG}eYE~&K z2cXu_Yg5wJhJJ3ou(w7aSO5y+Q*yR~gjgbAwqdt9CLX5PzT{#)TkHC(*40YVDj#$5 zYqdzUIV-0NU>EU~-(PQl1u)iY>%E`tAJP>=Wd#&xc5d!)_UKdrTcVB4&6$9Kq4yF% z0{2>+_uK6PU|MP3hR$hJ0Wk|WPNrQhz|7+ul^^-hYbngph96@Th2;8!ZXgf}&HW

iG}JxcIdpn|L@Vx}pU17)Pr+f@?I1F)$=f76n&>-^KVeHZS1Y+C=rxWf3QQVF18*Z}O?+u3r6Ixkkxr&uA1O--83N4kNljR!3}N(g z18qu7B(0&u(vdU*2y+1appKD~gDxf}Ht5v<{kyD!ffzo45b z&nE%kXV~aIqw4d)zG{nFi)%-uBN)lNG&75eGrv8`Vo4k1`u-`s&I;U*S zje4(2;IZ9TA+Jn?=dJ$-?DJv4RNry7MKA_NZqL-@2HoD!&x?qNP$efP4?2~UmdXQe z*brAavt2Rsv|zo=pk#+ELupw_9*sFQusTh>V#*OJGh|=QGRWXB9*y?uu!Qu#CCP5g z%t&)G{KKE%^U^4E5|^RSO7LPjnsH}f6l=y90TG3d5 zwT%rk7Z)L<;`s>6(dHyE*j7j|4udyDWYz^E%Wp?))~%n*sK$vN9KTRQvjF{=>jONgsA=&uD>W52(%TZpY1&GcxO0}83-OR{dZ~J6)1MX2RG7y7{ zYv0+e4B&!NP)Bq@aCOl-5E61xmq)Ng-#cvM5+6$@O@z7q{bIfAz4u!Maif##HYN4_n2 zIBIcj25|3YLuvGxn7Sry!_Rekv|ec(z#~7Mjfgs( z5X|#oC=`a;`h&tTqUQO{h`0y~C{1)GP1H*f!gc{zqcrs$)U+D=1(&Fbo}AD!^OFnR zF-m^72LKQCNpcV40)Ond0!{ ztgWP)yuaVx&8`HEv#*FyAre237xjbyc8KDq_vXWG`N3u*3yQqQ(Pc?2Z^Zw{xX%}h zy*CtE!4k^pEO(lQHBIfGd240Gqk8Y03DBSTnh_3|EBxcvkhok^OrT((L>b!gcXohG z^3w@DX@mKvg*Vv!H_g~^Tgclk)PA@Y{kX%hC17a;PKD5cMDkgR-LF^4z0Z@-Cq$pL zg<_S412IOAY?s-G-v(8WVd>DiW3)VotXNNDR=6&%47l{11(oZ+Gc`PzvaR8ScL7uY zX%FyJ|2`Io z6pGiPDh&27wh3E|>Sm>nwh%FXj`DfsPpa=2ae=AH@4qh8P8D@!kz^q5enTGv$PB3f zV&M5{^Z(5OHd5+ST57Val&oX378{M!A9uQs&1ZU6rD^XTFYCO3@xQ3_a!}L%hkCTx z>DWAw4;j_x$IR0{La?*^9@xcs(8!zrG=Ra-co4=*O`B+S3p8H$EVpfzx4cOcu>GBV zdv*VwX4dQC1!ZDt1cb_-UjvIlZ6w3!eCB7#W-`O#4PK`{p;W%FDdD-o_x^+9 z;R6gdNMowNp*V72P|qpO)3kxHCDJf%OW{ck!9eyxV7OE=+fR3W)2Y8|DO{sr8Nx4U z0YMqrS=;uM*z)u%nfag4FBr&CoB6xKI43LnoqKvR!~HKh3u5d%9_l<+bic96S#Uoa_U6u zX)cskm`(BW$a6w)|^*{;XVrIYBuk5TAdPvlOI=E4rfW# z0^vV7Zgu(%VBU8fCQ4;}e8fPZC21HKOs+0lZY5BC9UXZf8D$t4gRWew%4NGh7~mK{ zT*OzpZ5{D6leUr=V4nDD`DdOz-_xMP0

m%RM$f7Fx`#&*^wON~djmK7O%FykZ- zlpJ$&^Ky=8_%fDno}j{<;%J@V{v?m!F%3{&SwTHPFV_iG6#xNaBC6^6$gApHk?y2t zs5muNMuW#ST8qlAx9s)pv4|lj6$2%IZFw!pQ88j(%MSrU6#KJSY>4B-#kGiLfoG#0)xczzvrrgg?eZjpPLk)9)!pxtrS1|Xx=$hU**ZUMJX&TK#K}& zLJ}1xN;U2FM9r6VT5Ib^H}Xlf+++1<$^5YA-+5pEz$*4o^u`=o_Ku;-J ziWqO^&d8!cajw_XWs5dbfKnxWP5P}*v8h19fCSkxFF2hKfJ_n?)S$m=&YCpc)C%(E zyivWzV$&wvcCS&C5D@|<*Pg45U{&Ncw*f-`I=d8HuapYSYm}N4G~Lyu#rRe&LSel8 zp;;#&T}l35K`{cY1hrB6blhHO5crC?{&l^ZZ%zFWeo*0@B$N$j`<-1#SH=HWm59e} zy?>iGuwhLCjwG0^jBj0vvP_Jk5~2+H1?RTf%=7Z&CoywT$KtI3EccKfv_@Hl*HrIX z6wXfbJq>I=&nS|>#+>-Q7oU|q2?3nliE2}_(8GrbcY|k?$+77i9VE&LeOM{N{em4| z1l~Fcy%i--Tcaz;d4>y^RTyGTW$dM-ABEb` zRw=jE1{y~Z`4k2rqmA5z^B@raq(8lt^AiN)CfNBHwzUxzww~%Vp?#O+ZHy_zGDYr;udH^Gf zPTUtYa40Y>Uc$^WZFm#5fScXcin7E8?k6gW2uf!NAXiL$6X15xGb^Pf|t$*oHd%^C}rT_;EtOqplHJP)`X$?mDO z_9Z`r)!zaHa_#4;S?s>+EN^GW!Qj6qI}1z$zZ!kqIhQ|Sx;1wbrp(~PD7RK@5vaFJ z24}V5KbiW3t=>Wm?JEB=L0B#WiBbMhT45s0oJ4;fqEM>LweijS@XufeaEh`t-)@YT z@Okfb(-wEU#{*n@Ep0=;i%&c5md*O`O4L)Q-c~8~+~S=MDs1 zUilt3Qt64CBJe-fomaOkGM}q@PW$grR%R7{9z7O&%gOTIFq&>n(5;|2bWQ(wcj>Q3 z^Y8bQd_HZZTSLp^XOo7W?X0-KD3%(wcw$z?^e8tr5KkGhH@7eLC`1xSG@*amloS-O z+}zyK=iVhIVgTOL-3`zC^r_lc4X9=BxSmdIoNzt;#w2CZaNMKh@2cE-K}J#khI!mT+yIs4`t z>Tx2}ZX;e=>**4gs#}@Eo>N;(U1d8QK5qZwG|-AiXz?_{YccwH_n^q|tM5sVonKm# ztT*bch!(`poXn#wI*R3MQolN%B;MI;wjbE+48f+8+ygi7%ZUE6NZtzRRkt9k$Z7%_GwR*%f2aW$%* zlX_gnG2_R>4{GgtP9x=+s<}3CPgc=x|r7!-%P`)3OwzE^Wwg>p-{q1k_`5@}Xf9^Jg$w%#OtuF4` z=Ff%ftUnvtP+2NhwM<}B%=t`{CgzTah)Av;H+lJIAPxAmgN|}lAGw{VZZ^}7(fY-b zut}c}J`NoA7$h3b`}Vqp`y81?qx^CJ`8NVki_ZHzxqnU_{X)_0dhCf#3}s0X@Uv)I zwh8btU@N~R$ix%FyGl$nNibKdquL_dqWbs=3qF^8r>ZG_LQLCqqp_}HETd3P_k@_0NcPH13*MZE2)GWW0PF&a`AStt9#1WSm5q9d zKUO-GMj2X4GXm@hIh`+y=3_2?klMHGaGH`<0_cGvG+C3PA& z3l-~~4C$xV4()x0W9(5|0yz(X6rl3V5u;p&do&OC1T+U^h$z*-#2?vGaE}i333BG} z`e{CaoVY#nyj+b|`LSc8#$uU>Od91c{jkt@cBO3H#ipAsxSey_N@ z1T3$gdT@z*M`PcoL6;B1^)V|&RT+h2@O+Hh$ayt)*;F`fA72MyA6fF|P~kYhk;7yR zK;VFjl-^}lk<5_^J~W(#q=3jLfI-Qd21#u=3Q7qu-;I&M@MdEoNNcX3XWBy6B`ogSVP@wyxycdJ0R9&i!pReq=O42VIX)6^sHi(s?xMM5~c8 z-yQbO4w$bUCTLP3>lwU$vJ&m&E2FvL7yt@1cf1=PuGWABoM%178%GR8+A7N$>75^s zXqIGsx>t#|^DVNW*LWep+HTN?DfW>b%aqrm+9t%j0tqDGw4rSYBp7wGOV|S!)%q$} zTi|>4^{ir+NPMw)+${~`c(g9iZqCZw%dsT43V3K<${rgnL*hr7Z=srKjyTyy@CX|? zE5~UhH)Ud4a`~$iKH_8SOC0ee3EKM!>w9pU z68SG%)_o)b`Y_&U=$X@__dPh^Tg(CV>r!;~+VDk)er+vf>_*kKb;JZFxmxa$(H2Q~(Wx_cLS4 zbJD0Q+AbT~(1X3P`7i$R&Q7TPEP^xcG;;5JMzCa?Vsc?ePv6a=?I}}r(Xbk zd*eZaCT0~2%#D8};=Li63b3j53x=~6G|~avnAVE9^=VX$a&>tYds=`mC-oH-!78bG zjvdN)+W8hptOH2aB-iUNZ`UhYS5KBo+#Lxu$0$7Mi4njdjgce@nBfGHI&r4Qb1$mb zku(Xtp1!A6PcEQrl&QtS*qKZ}wL#co@2Sxk5M+-92{FtHN=@@@=SlJsj=QM!Gkv4|TV zi36N;7-`Wq(7RT^eP4fsQk(xmko;)ayh-VC!m#Ec>qs>(Fsg3XEWj-+n6tcp<895< z5TiKYY@NHdzW(?Q>WWQYe}A$JO0DslH5p=j%XAUi6%Xv?teGDJ_}?yn?byU;k4=qK z3P@n(NP{}9tE;(jyqf7vi(ktQVmMt0k^vyXfMk%O1vN9cudn#L|KXq~BMcYa&D}j& z22`W=-!%05blm6y6ZLN>CYr!ws%}kVk;i#yiQ~tGE_D_H8N4sGh?^qhUMO{=8>k9; zo%lBWYVhI#a#8UtfXNx0a971*vGo8 zyw46eO@d*;RdP$MY2*JuB&QnpyhYU6}~_g=ejD3k=2We-l#D z4^?C{M?Y+v;Z!XqVq3Q%A)TUTd%#_qxWFUgedITor108)3lZTvmGH zg&74bJ8>J&|71|h1WMZ5+wW~v4O?Cut&4i^F`+I}{l7K05dvc&LrIAFE9JL)r)5uH z(uBrSP*L?<3(bDK&Ba)=-G);D<6ASJiLC9dk=E&=s@2JSh$M==Jc`+ zh4eG=Zx{mtu2_YgBtLH&(bCeg)Mx)m*#z!yP*`GT)XbKq&xOafKP+zy0d{Io;{2qV zl2cKkg{nUOFZ=J9Z|@87fJ(GqH5dkDynZK+MrEZ;hN0}ypyX=sLgUu|k1kz<7P0${ z{6O}YVKU!6SoK=6XaPiPcYkR_iVcEZyZdgvnC@+lkN!gZ@AAs+D~^?u{+_Gz?efPS zSpS^ASg9nxpi$JH=+5TNC`b%c#$mKZxi^;?S)Sa!g|(!8PBUC>#93NWXu1E$Xre zIXF1TIiILc#DDe6lC^S9m@k#^{K!bC*nSQ~!(J_Eb}1al*-mF*F*ZvcH?UojQdGdN ze?BoXikvj{BMiD)W^_HPruuvnZvkRDS{^5Dx| zTK|=z+jzvEHGNkXZ8L#Q71S>igk=nMbbf9vQ#>><2Lf(oEcYvKfM&Ec$l(rMMuzqj zTL+h`$Migz@(<2F97(p8D(Sd@O4_o0A9649LVi-B24nKTxM<12f&Y8yC^Us}El37wDxc7~v9mQ)H3>Yw;^#@~Iu z89j)fY4eMH)|y^yB#`h75jP}TKZyLIO86`2Pc*7+`yddlbzIrgpx~}bM|%89f>V|G z&|h=yn4cnL?tryp8 zE;}l@s0LM9c<{cF-OjHB_Nv(Gez6AH4~A8C{E7fjQpgF`vCJN|mt_C-lg7~I7x!+) zo_V9!#tM;XX09gAR=an?x+Tz_GJ(s{hMPfjFPTjQc;8DeuL!uue8rn&phCtJebKZ0 z%{y9xwF^e%EOoHfMeOEMHk;G~%Sa^J-`B4Bb1R2M@%=Btea5g&m4GAtuNm z;Y8zj>)*K5p~OnTJb>3|B@IN1iQQOnD3M@t*72p%pSLt<NW;M3l zazYkuKbmJHL~wCZE5EA;D+{$mb2uQyM2nH58k@XzMNj6=`lC7)zu+W(XM|@dy-kV_ ztPU8kHh9K>Ikn3c(P=)*;s8&Wx^neNO+(WNk@@zXl7rg%L7Ifmf?GCrNPC-E%-zs}6 zo`SBsWWS?yL8C^MY^c+=I)bQ#rsaHI6w-Qoy=HfH{XBVIj^(FG=ww8@mkM~XQbmNR zs^r$hwGRLl0VF~zm4@%lx}4uK&Lk1IOmb_x9Y|MFJPP)yKAsitc7KOTR7X2zo<>;m ze2&kTI7;uwS2dI=oIL<8Y8AeB=GDEpAuFY|xxb_sG~I!Cp7h0Z+LPcev0ke^WyJ-K z?%QRLzf!C^uh^3Qt^3N~pWQS9D8i>NlNtiVA(xzb2G#j#oNO7fK)a~1FY6_DB+agb zAm-;UtoUAPX@YfpP{{`zaZXw4R-#P;Ti!26$n%T`&w510y^pRj_U?ZdBa+{r4*1P- zXMqC?GHb$1ax4BK7|fc&HLW|?|0ei=rOR`#n?K)~(qL8^ibxDTJ;~|{bxwI~9Uexo z!w4{!_;!zmwlVOTNTW#L@%s_fG!ZC0w6^0$AdQAEE2HPn?^TI=ZI-B@nAww1eVLH2 zO``qv8+}n{Q_w6%4Rz}PG+=K!r>G*`$3lSigof|Ko*)&l-rLRBN;cL%zXK?wFk@+F zlf%iPTa)3g1+-|RsmXSKS!Vb3BvtFyh?w{4XO2GCLu<|nEy{)Wtj?upKWWj#25qWR z&42U9TFeRZ<2cU;2>mQYY>67g?lykKFw1JnIjxBN{b&DtzyBC#k8w8;4V-t0<-@(p zpbUbB*RLYAg@mL)<&}43ATrbr+HpCtH}8FGJIYsRC_ZFj*LmAuPc?s0N|4amZlI#E z_F3@vWr>%?woE_BWaO$ow+G#BvT_enXU={|OoL*?!&D_gwVi zx4$&p^7^i#^`?~I28#-sLxMHjn^Xt0RxL4I?H-y;odr80wy*8yX-g+yoXbVT#1?aF zHs0PlN&~dXgtChgYASTe+#;s47liXSR|HpRcyo>jNw-6#tHtS8yaIOW?i6>zk1Vww*VhAV+?3Gx zO`M&j_swmH;eDvh0)1%n+~^Sl!#dN{1-O6^6|9J=_21J^Q$xSmPJeYaRNHLSbS-+b zx30V;7C__Ge4>fm&3=7ZemBI&%H6_3o9X4pWvY$6?l*~N*-f^|RPQ-3{1&N1;yS;! zIc{h5C6r;HB4w>{s|vI6y1~`&Kk;WvN3z~NGklFzYXr@bPeZd?su{2A3zp2vhkDwJ z&U*)cTUg~0jxe@PbR-*{p3@}JHGPFIlgmnOP$&EQ6$>+MT9-WM*K*vxra{|Y8_Zr0yzXzhy(!|92K>hGX-OqhYa=V^BG z(xGn3{{GSffB)J$G@kf>kV9>#$jnzN^3&0uNZT6Rzz1T#qqNbh+v@{KfI5cKC>VdO zd%1k6d&pPoC1TQ}9Hk6Wzkx^;*aCZu1f79Evq+r`>dzbf@!xacWm}?J74`ZAMG@`O z)9la=`*!5VCfCtIM}I$Sq7=76ikiPYw)c4~X5t@^GZuH&J&f&J1LomosRws!)VCr~ z)3u{{83|qm9ea4HQvohP|4`>opWMBeI#^}-7XQs@%ptq9UeuwGlcSW*>oUdSJ?iH6 z)MgFPzHqZpBk{cz+e#s?J9jWy={~zyfB=RHVrN@pY~CH7GQYU)f%rr0M3WsemFpgF zvBs9c(vc4Ox&YxlEXSN5FIBwoem>Dg)cqC-Y*1~e+*YFa#_+lGM!?WzuMr=p9(xNj zZ`v$`Ec(C^1sr9PX#Ms6bwfwDTILcEU;m$+V_~VBcqp97sQa>m*OsSd-JJM z;y9CpyB(jTqVYh9=$QFp-&aj+d~FV=cRzG{m+>DZ=Tt9KvPna7U-8R zyn28j<6$)|ZGnkSs-ZWz++St5w>m}Z?<1UgZ9nDF2Nnxpwc6b=veAg+?iDntR4mx| zRC|^VgBjD3*KwBz+tWHXw!qVdg^wtf3(`|no~AIn(fn$^Th#T=mYd^p*Nrg_0!C@N z+heoa*XM_8Jpw!2N(^6;#8dt<3pV_tkCaQO9_*PaH)7 z4A<&=Fn<~NlS2r_E(KFBd~{*ObmbI5EE#h6C9%e z(G}(=i9fSrEZ>0_Cx)t@Og)L?WlW(wOA;$fb=%dA{eQwUzTyL*Wmx09)4T72> zW55v zt@ZQ|sljl;3>p9+c{7}CK{JC&RZq;TaCtM-04X#YK;*zP7e;%tEH(l>a>g~*Xt@?Nr_NJbgjoE zVT%q1p_AF<681TF;Ekp{sN2&Z(|%Bsrgm>2T?ijsM7BTV`@gRvf4vL<#vsY8d@lxW z@NX3`MYG6jJD$3ZTyzYL4xa4Sthx9QLRwjr`b3rJ51zScXYJ;2)F_k0UPN7XB-KBYz+Fk!Lu4oiWjh zjo8~nyU>CSa;K&xd*ncWNL8#?x=t%G-&?u{DpIU+w2j2Ip9tZPXn~R~tGHi&N(>{f zD$GA@^f$)ew7|hXvtel1v3m-;e;=$Z&~ilGR+KX0iV6##AKrk=%`fH|A^*3xXoR3n z&Ve}9wV31Z5``_A#cNvzc}^^=W)T}(8A+8 z&yO$|WEKf~gYf%$ZXhI&_MBycgRg2Z&)Arl)*CP_LNH}<5l`B0gmEf!KYNSw|?J*Ni^+*g*K#31i?MV{raZK`&ZyoPRj|3$q<|qDJh(VIm#~+TCXM zxW6}24w~j#`Sn)7v9zVG`S{&$9@PD4;s5)7^w^PD=)uEy{XB%k_E}1l*xBVoS*nzu zaOp(n|I^Bu1~rwXVHjvz1`vb>QCtYhq6Rby7#3M<8kJxvpe%+hiV8*~0U=>ngLW&j z$*{{VX)+OMKq0b42sYS&hy(-?NI+1;giV@F!;<+fPIZ<2b87zNR;q4p&OP^>b5FhB z`@EjC6RK-vue*T{%1d)xaLG5?w4XPh&-aaN#*)dAw4(L_$%gH&va+Q&rNvXvCu74) z*K%+|PKlRZ7VR;i7CpmwEvdTmQ4#FoY_%k)ZXtO%>|pZUc5g>#XDqv}y!?lnni@X~ z9k47vSo$uCPm1P~$P{+Z6{q>u+RLJU-JilRv%A%$|6;|F7@U*p-EzGf$l}M9y4E7v zRwqW=;gWBR#F9j6d}0P=6v477N9O zaIMrJAZkdQ6{>r+zht82sxcakJ`gZseYKy3+k1MSuI|sR-KqYWfW^Ho^@t;^2wgsp zVrVhLjnNkXg~GkMfTSL=HvL_V5LSW~;x+vH%}3$75X+)(OdW$fpIsi{z4@OG&$_Qa zF2)s9G-;a0p4EltSx_LkKO>HL;=wOpvv0lz;b4+JI7af*3C#m3Do<_v{Psy`SXvUG zBB7+@kb1?}Hz_TxXKL*Gq6&X)$#qhaIhHi_%*;&1bw5@GRIk-o-+|H1-!vC=NSnV$ z>uBSw;DVxs$$<$)etPmQ126qdE|k^kq6?5;h^4p5e~mYShYVRq1o{i&aWy41xwA31 zZ`-&)!Tjriey(a@rb*eZfmAhauG0ASo6w~{Oisp1IU^mmAIwzW{O;Rtu4mybyk36m zy6DF1zNIzm^CK4U(TlW2e(*t<^~T#3;-(?H5jT3og??7BH}Qtf1X!C>2P{J$!)PP6I%~0M{FKS`1lGw%!$Z5yWHY z?l29q5dl%JjW2N)K!SsZ#e(&w$(+HEvD_{&Mlc}c;zaaDRC*-PbhdarzN@dohaEl< zqt_pjmTMMrTG(5bAH=8h;;Q&mH)yAt^+pbU#xELK}{ zUJE_9m$N*IN#>I>fa2%a?6%8Gg**X0aNFbHhV3&F%^m!B)efN(kD3cXP$x)fln^-e zjZ68)0amtlb{VpAa*Us%6>ylkx#Wx3YzDwAl)k%%8_#LIU=d#x0x_VAT3U%LA`jJ{ zc5raWUKFn`axi?r0@_LJ4_Us2kgV05m3RmbZQb`J++weR!5Pszl=%Meo)kO1S>f!L zLV#>flyO#TgWzc!Y$5LKH}QCUCC%_+i(>2eJQX>QHCy}787~3EgIuAyoj+e-2*_)D zb5{H=nPci9@0Exk#KXU>EDSJe0>`eH)wTp=>NJhddrc44!N^Fzmfd>s_uudp4}pxt zqAK9idr>+HPOS_RJ1V#3RCaT@5)LYz6+U+I$FtnMa(kJ1p|X8_ebta=!M6MVJ{{<%rO$T}>EcN=z!!F|;U zI%=N9$s_lkl-hfE$U|sIPEL+@)!6&@Z(EWy{Ve`Xr_gS zdh0hIm-^f*vvRSklIzl88Hz(-VH0d*UpUP>d^c(m@7B7&2Xwn2d`&RzSgpcA; zB|Ebj>&*#fkDlq51rLVohRVTA2{}S8S?lT8>@ZU6kNna`pmyF*VAl1b|8j%~uvj$% z&D4i^As^I`D2egjvgpQPZkREy7*Z&L`IOC`VUvgEruc5t!&{DNC>{aw+B0I_35nIH z8j3uu9yW0hHfgTY1=^2c>@9yV;;b17kb=)NM{I#q{q+d-?$ZTHx9+OR2>5KH7v2VF zUEOTHCmtWe6%_p5MT z%>ML{mH9GS0(+`f4h9-Z0s+GAmz6aqM$T!PExfsKAoH`jsq18&vq|pDUs__LGLv*5 zKLnArVTz{>NL`lzjs}Ln4J_)6y*+ngRl?rgoj0}CW?2c&$|U(w_?K(&o{@o8PSQJd z5-*QYpj?A!{t%9mW;wJn6&&y`!HCXcsQ~5ev8=3>cQ&sWJlYMwD5o|kNYt;2eBY0b z9#oWxx#YQR>sHtAw_j^${13UxCPv%v9I%i!niS9~}g+(M~O|s9=bfJH*{HOdGS# zm*dKIW;-u=d)t!h=F@{lnq-J=>HE641wS#sM6zVMsuIJftaR?@&^RobI8r8RZh5xj z#M0t)29zZB;gYq?VJqL|5zJ@s2!s7gpBgu=t$k2@Y%QsanPOykCa>SHMFW;rVn2)6 z1^)N64!(b*d>21>h(ED9{bcq|P>+)f4sK+tfp`2|=W&$afUrrmJx5`?vAq;m*a7V; zc=b`)h1JpHzu{QF>SWFoujWiGEJ|8z{(gBied|z>fT^9|!EA&XY&eUKrO&-jv9pEc zE^WI)bZUM+24xqZ*ZNCAz0e5WhsA}4Z=R9OWugb6&{4LNa|H=HnLa=Hqt)- zeCk(nw#LEfqp%Y(^=khI1jY%#jr#vP2*uKY&pdJJ7VWkWA}N78ddj-g>fDY00Mq43 A00000 literal 20382 zcmeGEbySq!_dg6@Go*x|fFe=~f`EjSNGm9ygmjl64T3ZX3zuGV>PplUOcV%$P^%~_ zY9R<3gCM9$a#Hw4toQp*_#k$_uA)N@e;$%sKZE~2=BjMyjvzFaguke4xh#7GVMkOH zuj%+CEdKC+IPfd+=wQp6>;3yD%*Z)za_(#Nv^VIfanBEWwJLO$_ssC4dQSRFqb=h` zn>c>wf`ZDkRRuTL-ugM1){d1L7(X?2tWds@XJH)WgkpHXK1+N)>It(=^!t^A)}K>Rj)AXO6%b+sS!G1CBM7FI7C|A1bO=I@Ad{E>-&g)`ON@(qjL5ar;8LbS=<%H5 z;yWRh*8UO*lES9iucNK~nwOXN=bCTXshbG%Bf=AadV(hu%%zziqpA(i%;^n({_kG=`&%!!3 z|39^k!Rn5V#~0d$Hm&^`5hpw6&T9>SmwTOPh?L(aIME=*kCl~Fy=4#dZ+y5M@Zjcq z-jLaQ{>q5cw6azFvQO?!- zpnBa*o(jbek{*wS$f}7@v{B`=cJ9;Np5arMB8V|R2E-~a9M-QG+&-m*PCrigIZO5o z66xVvLWXIkZX0^CC194*wR*@_pgalf^5k=5`zuGJ-S3X~ zAVe)q%~!d(x$AWk2ZcjpMJ6>`A(kF>Pp+p@MKF4m?aj4)PR3^==(!J-sg6iGn%vv9 zMPRbTs&(OTZ)GTkN836sGw^UN86W3Z*o!ML77N&2XI;_O(MeB~^5$w3r}FajHQeH{ zAdJ8v2l@2WgGF!hK5`rWbsmGmU~_a^f&4t zD36UP9*cfa9wzJ`lglLEbw!y57Y6eNzkkn%A)DRz{`*rUoK4Q`&6_v1zvCRsty&}c z16N!n^1pof;&NJIVsGO|mIti-rA<c(jXiT!gE1qaA06g zBWU+^*2V?^IKPGjVe&XfPgmg?KzuZVTL*nu&8Tm5P$dYByT?W#ty+I+Vp_S4j1LI1Ox zX`_|S@$KKf^-J>(zp5be>C3j{<$qa4beDx4;n^{|qsQF$q>JXM=;M&H)gRE+P9wsm zzQ5iuU2#c>i;F9F802*xF1{w%xc!0C{?Bv^lZegB+`POCDqG|JMFKTHPpTzJ8*-~B z^w0#Z5-t6x)4wwD&|N%;+B(=iE=&wVr^Jw*e$i~|mxKfotusa&#D9^G!hZkE9}4!? zj~}T)&mx(xerS1d8oxT#%qt*(>v%7+G5tbtVTSFfn$v!wG2rh`bW7JvYb5@yv5Tg& zbHUA*jHx4KwvIDj!W&$eg8c(}m^HPHL#Y&)f^Ih6dNIu)AMu!oU8N$tAbMIGiz|(y z@ji|ud}Tnb+RE)2CI>ZHHix*F3z|2Ye0wLTG4pJa9=E;{JLOC9Ygace8##%fIERq8 zr@4m^@k-A*#oLO_YeXZ=OhTN6rbOiHp&-EN>}_uUJoHz3Aw(s2=^_6`{^LA?QTtKK ze2fw-DB7gTx^JYkD!e7)!9Ap;m|7?F{JyK58c{z5sk##KB+HQBBA+6q;@Fcl4C;y0 z2T2_Bs~9Q_T~`VN)%Re8T%!z!s57(3k$jFr%br=j<6h0-ID?J~z>$5FkI-mdbH8wb z4^#dqxc>w)W!M}0-k{Y=4RxB>UG&oAJ5H9VdJSy9L!IgoQ^+}$Fg+n1WXgkB-UpXT zSMc29GyUqG&Shyqw}1y6`Xa^ag{`~OMiaN(=K=?rM!C&HRry?6n8;KY z!q_o?D0-1+F#&BuM?PsEMq%N$g~YNwWBV{T%G$Xdi_tlUd@B2tC^())EDx3{CLehl zbP*wSpD8+nyM&%pgn8;ofB`zNzcXrxlDzbw-ux7Funm8 zcY;PWh&_56-7i~kVib#*4Q$aDgQ8NfC%k~Qq}b{sH;TRv7GDDo>d!@6jUuRmF^63a zo?b|B6?57UAGC(PnzhdoFF#32kaP7Kovbz%k53ymSr*u4P16-SC#aUskVP^bnCl=~^61cDwdIG%R>fL2C`rFjfdive`X zPmuq#c0OPo)nQA)xGIXTnb1GW2-6fMOtYvAL+5M(J~1gK8`&&B#`z8*pM~yUtM)Ih zU9wO|qw4uiiM)11{NvD@WygrSU~6($D0O~L_LrutF55cIAdmI0)!XzdM7%ap;8F== ze_b2+U+;5&B8VQ0#1~%vdBek)nnE`erL_W7~ z+Y)VCL}&yiyR}bvCyLEU#v>BP;4%97KJ*tM^?=~+jk95=55;jPbW$I#byqHDFCT?R zrgT7nTIRFJ?PyEJ_N+6alPWlO*{0BWbHeXZ!TVt>iyjrzflqypgLsY9`3a`+&i;$#pQc8zNpFU$`6m%l`cPW zZ@f5V-}~XlXk&nXh*E^p)gSc51HVLRc`=x*@*93&*P@+BJMlybZ44;asZUR97jpAN z@3qq>dg2Us$F~}@9z8-g{P}WLX8bTUl(mmrZb#*Z?`)#vY=XPTM8MYkBT9ju$e{fh zsz~OMlP)U<%oewwQLqz))0lNi=#=9ohQ~tU4IkPJDPrpdd-MtLr5Wo7)OtpB? zpu>HaA%25m`FL^Xg*E@sW=k<~aj!b;`}4j1MW}Mw17rSkh>yZ;IbIB^izuYIh6`~S zt5hQ)Apt*JyMZ&_SS;&y7%jJ#8LzCY)XT`o_)*Bm&#&w3?3^r@R$FV-($W%k``*2K z?O!50N_jB{hD6F|8L`HryJe(NYV58OL&Vq|?X29)%$q^G{mJXU1cMB(jQgaG`{4^> z4qx8;_BJOq^_<&CQ&7;YCOKj5@4b;clZdkgA_(j2o2YbslhY9Oa|6{iu-D`gdW;wy zC+%CxX}`HRvir^a@Ixr8pK_$k<}gncKOdjAsi|odeqmujM^6to+fDv2mdH6PTFTXz ztj3;?B9_kaL>S~9#Zj(ILpq4urBg3>n2&F5d+X}z8Z`w5RN*019Ky$4tJe5HW8t!* zfx^>c#Q&`$#HtNVvYg_A25aQGu!`ZPXKmn)3+j)dz}%wx*gxbxln2-;3BTw462{Xq7v3Q|W6A*!#J5 zVPbB!D=`mxpO2r|&ul#0C~|C%;V~++3AgEpeN|wNI&&)Xgeb;ry@*=XmJle1KDdGG zk8}s(W+(0s7f;3D3mZ0#O#PQC?3qq13YNAc2rDCyLu{bM~XWb+X==${V*R@ehQB;}(fD08)K&#))NA=kf*+5r za4+f>7g|J0ZCZ)28GA`GPc5TvbMEprCZxI*Hji;NJXeI9R^tl;Tsd3L5sWE>u7(K3g@+aUB? z#ugT_nc6ZSS?tIs_mgTVdvu8V4CO~dB6uWd2>HEDWk^Rpb+L%Xm6b4+CK-^6J+cm} zlqz;IC_uuLq(|uG&zoRcUD*i_$ju_g`Fga(^4PORV>8;=EMuaJ7vvO6Nomy;kteh3 zTg(<=FuOZ|3 z;5|wE$ZOw#tv=>VM2G~(bo7Eg#{V$JvEg0m{>ZDd-B0+S!?_HsQ3b1bY+j9pl8F!kuY!F%bFj*3dyDSjhD zQoQM#{Dx1zUFfCDyE_+SQcrNiM~A0q9CgdZ&b;n^c2pDYiNWMu5XY=r9VCCI7fR@a z{NjziS#R`(+k{lv?nkYs!~Oe1+_Ha;iCwv}<=gb(!+93jhZ1fhb>8EPcKgoyi7KIDN-GI`RZQz+aS~o;LGI9UJ1)FC+}9(rM&FWZ-&Q!+5t46=cyF z{yxR;zaqSt=$O-AHQhPeeX*Et)a7WqL@q_xxj&y%OoobyiYiCnj9?y_T~bey&iOOb znpaw?^sNJ>iyj`pqa&@|6{fBN+Emb-CM-i61RC3^4ypQehL&YrC^osxx?UboBM{QWHDsw_0Cx+-4>jiUh6PrY}}Rh`)E`BEP{>uukSuOgdNn$ec~`R~KyB zFZ~k1B9*fD`~RCc7gZTnI=yFNVmkTmU~h}Y^Xu0~s3O!N)pLHZsf9SEuqN>*!^CJB z*Bc}~qG|!Z_{@9Dx@GgI&x+eke-0x$WQEuP(1(6qebYg_GXM~eK2F{Vp*UMCTd-vd zd*_$89%4^KRBil1jO2Sa0`WaN_7-V&kdc!)wvy_a|g`9=1L2hK2V0#$S~U zu?fe`mpuS3RJhIDRErnq*KWS~R`CN-Ji0RTu^kLOeuH^sg}zu9_0m<9lNs{ko^th- zWpLVOYXce{?mXZ(ZVgE4&;tK zfvqb5on$yR3*w+*^IXA@W+W9S`A?JSWh6n!OP8I{RFIv$; zWMP>dgz7O;Vh(@5z)t+h5{490D21U;3qlnO#qjYMvLKV{IQRPOWK^*uA5X;N$Ni9f z@r1Ep4-yXFL=R#zD>utu9+TgEp9vlm7uwiAKdQoh9eL9H?=g4GUJXj|7Af`%!O<84 zJox;}bk?m6Zwl~;F*7f^Iq#sJ{A-cO8O3oS5A`a7k=h*ft`v6!U}Mu5R@qei5ws)( zhOS+{8}&B*kwqy}k0JxD>jlR!CCiIh6Zpp$)QN%#sTvQF44J`med2uD;o){kd2Bcn z;E==7b0o?nw8s@dqEzb;8$+fKnN0Xe2CT)AIL09QpSgSrzVyu8RqWdjU7i?l7m1{! zL_H@7D*XR_Rq8C~RRP$NTkxM@!^s`GZL>NJBW$>A2NB}Z*2PSvxF1Tn4c1$|F zWKBoTOo*S(bbxb*CznWnkQJ{D5WE@_tnq~sQddMtMy2iP6u5}*JXcAiv9cHlHWUal zc?Nb^A}FMRHZ~RIf}S5G{UX{K6Gnr5!o1N2*elvMkhTdi}{l)??3SKwd|xYuJ3sAIciF(#pZJ`Q$hxJ%h=zR#1&%; z5X}Q?#y@iXR7|d?^T;f&2`PnaD6V8#BR?OXOZc}SP9#0EB099zEbNmvA1|RU0$O6k zCV(PD2il~@sF?H*v<(RA3X>C!e9kuk3y=@gI;8Iv1Gs=!x5&hxu9ToHdf9m~xDet0 z3a`cQ)$!wyH}f%K6Qi}UkFya72L5M}Q?KbF_N3I+y>~KUtA5HMY=+M&>5QH>4)^v zN@jO1&7Ms}xZCK+mQf1}p#-Md`fTh+eT1Vo_QPfEzzMASNmOYIKW zd+Ag^i}&}`+ky=`A~qcmNmo=<@bdF7d70YV=Qubz`Zs7uFhicD{_w}}u*5ow!wOX6(vjN6Ws-gB_CQ#RRLnMjUy;}D64A?c&iAy<;kv_Ec7_betK z-A75QKK&fAi_3JWTM{SW_q)w(hbdun+)10V|2xUdA`}W%=%`3_b#w-n(6hz3S*i;3 z@Br>y{{lA=w)Ydu{S50cWX7mV z?JxAnB=2U9j#{--QZ4Tx#Waq)Lqh3~u#@8_VRwX?XrdlMA!MMZCk4E3PJO+!o9{wS z*2$`os{6^u*<=&8{g%Ji5XQqAwC|SaJ?^`^Q4~}l=7?LbaU<>$AE1~l#+d7cl4S;? zyo-~*kpwgEkdAz6ku-c})8l$lUkBL}VW|ReTD!1rbJWpusbVP8Qn5I#@|+%_>;m>8 z+?mC}34iApAA+==Z6~F@{WTZqqhdNf{hWfmhfSdXx#dKR`fhWP-{O)UNpUy@J2QY} z-Nn3^kJp%qY*PV25@=(=Dik3Pn3)s)`%~&QalrdJQ_d31c2xAgEknp9I3T;lhKe39 zAxnw*PVoNf;N)z}FW`{lsp!~NCV8xymLbaew9vnXz~xRhYn2}3>?Ml=m<}C{PfBdTB6{A%p3jV4-0Q%YF3fEnEcK zPe^&1jZ#={-5@C@Bn_mv(3(e+lYE$K@FQV}kHLS?+-?!D$V1zS#gipqV6w^(&p8k~ ztDQx;R%X(#qsnoJQRieJ0Fb|7A0YIX zC~zkr1WHVA>@BcuOp@Mf0LY*0#(z|HK~E{1OL$CJ2N+5+j;vXTD*S7t1_O&@M7p&J zQkvGtb4i8Cl##h_tHezE{5tB%V+a7wl8nF?Sqyvl3ld z5BiIXzajIlsR9($_Uj@(OC0LcIYLK@l(FYPH7%*xCm;64=FmtIS+vpeVPun2F@scb zpumJi$g?IHC0;Nf@t;^6!dLCtTt8TKo@h~R6H0&8j8-@dOdv&P!CKIPK3<~w&3P(7 zn}j%FURk?7a`sH~qxL^E2c92|$QXZ1TSySYT1MHS6CY8eSTeTRPre2WBH^O$S@}z0 zW*9FZAd1~ywOR6B4vUEv=OPvU&7798zcTih(@XwN_^dsr+v+cZ{b9 zDWzj{Enn9~ppHTcI+e=*GyusE$O;xBp_0(5#UUY(E&OXCU*ng84B2qB1oacZYF=rv zA_Z#%@B2?4I!VOmFT)$AxkCskd}&D#l8!z?r;SxB4@M=-*>Nz)`w$?Hs6sI)$Z$Mw z>4^wwF`3A@%QV9mWoowk_2~bS5 zaZQckY@&CPUcSNQl)a_zHQJ#xZ0aVrZr$n=y>6`yBmTV`5Voc;LXMXB_DbNaJXN!D z+g7=(i5BXKcKYwMNNN)(E7tyP=OcUA6diq%noU-k8uAOyF0{HPhaFF&&%S}Xq7qZH zWZk8E6Q|_kGW?*7vhS(_wr47iV$y<8edQGt93tR-4Trn|Qza78?_IxPfK*TxcwL2Ve4V6i_N-)r1A#{F09 zcmDp`?LOM?W?d_<5D=Ut3J?Bn3BH9M1pQ5zb_kAKHP5?Ba$!%Aikj1uy ztM2Qgjzt^0TOTd%e6(0F6NNPgb1R0^%;;y+iUI$$DzW{L;jnKnrdp{knHEU&3v&3$1LfynsPWJ zcc>7(KKYo&0W#V+*??MD@%8@XL!F@QX14Xci6dF#CUoD$5mK%V1}ujY8pYf~fw^!a zH6}m<-+T|k3c{kh0yJAY{haTlxv}GhOZ}gl1wcK3b;kkRdNPPDOmUqc144g)*Ve?` zf2%gE_uc-{NcNrMGD$+aDK>$ON+QBOuFVw7lZdtfZ_xs6g7k$u$0L6l)jYUYRkTuz z7x>{hWNOdKun5YZ&_`iW@7FVcORe0rYga#mijt)+SaOK`GQoF~fmmxZV4jFP3k-(z zITi-&hdUAIM8RDh6w0_G^#n)0k$-x| zLm#2(r=%^e@Yq1g4cLUiRvshAxXR#LBw{itQp6HHpxA;b$iOp7+>ZKBAT)H_o5WRh zT0#P3mB^y1ipr`&e^p_3MH7a+#DMMEP*jThqvV1KSV4Z}iP~5oVL%@2<#CZi8w(ze`JX7i6E~|-@$X0!jfQR!@ zPuvMwqemhX`$Y@`a4HBwtfC@oX%=Cu2EY4M>aH+VPPs&EWwKSnmtsmrb%I8~s-=n+ zqL@#0Xu6pda)Jr=vf}S#)TLF$tIFS{W`|rDfqu=PY^3$T?W$BT%D?2TjZV57(KDz) z{6}l(!0Gy8_bs7FUKt(6ia80D5<|i?Ug{zR;p`B+7+;N%AmMC45sD!Y-&7|JzUKP@ zaf+MkA~C7cg%LDYCrtOXJUOOV4un@wQef(QMuI~9R^>Yhp9F+_ScS_Exn-(B2u3%d z+?#hWl)Tv{n$dQFgvm6jag!!m)Jt){^A{ANi{Ava8gX7oPMa5&mP=tyPIXmMHQr>^ zz;~lV0IL|H@K)UTiTCGJIoj@35vJlOkh?#fZ5cgXhk1`iFQ!7z;GCnQyL&x&$bZir z)mY?n&Zlk!-6p)>E{0l#%GcnhrYqm~TZT+GY2vB1vC8kwG+gaUI4@<_oY>i?c+K%(`9MGu}64b>$_l zt5=dm@#>rTW)#-{Z4aC{>{VhL{2WHNfhD3*j_$fT%iI8Yj8Wsmnp2kyq6#g~E8v>k z6<1iC^(=*yhPuzZD@m5J9}pT`luP}w0U*?IFYiFc_?T4WVoM&T$#+fu1BHbD{)2^ta?(P*ZB7UBYI_s3b6yz?Y@b7k~s#qUf+6mF!oKdk3RFp+i)LoL{C?D9;gI> zeAIg4b#!$H@d5vE^L24L)=PMmH*Bl6)~$Sg;+JC6N8_e&S`GF6@p{8Do6f#`!#_XtZod;4%r`7sFSs9YJKAG=bxM&jCLI8a zg{r250sqJtFvopE$@J>s&kO zSL=P(c1>)1`+3+PmgN8|m#Fh>dCCM3x)Jj7+PS_wse(C_^B}ZHxT98rVUazm@3@a141!nQ?L2`q^x7XRoC>g5fLbS8Rwz&h@YmkMt?s+7F3kIZ{kV5@3AnKIK`4Q3( ziWHU}O4_die#j`Lv&CX4%tjS8MyZXJ_i9LvMgOL!l~Ptmm67%ckYfZa19Cr5zm!p3 zc-RjNts6_4W~&%gFub0V68sj6oWAtTAmM*AHN{qH8PV$mVv8bTa~jynfw7FqsyCnb z5X-|r>Mokg9RnVBb;~~YEe%zFZ|~!_`*Onb^Jhq* zZERd31q?zQ2JuRoBH}2+%WXOKu1a4e5wSe`>q_syxj0SDgkP_)pXiGZ{;#Q-s;*FB z#sI@Ob59rPW-*R|2cqqWfayhzPPyGf`Xw-*Vlv`!`12MCI>xV0AmfTq2#*HQserq$ zgn@#i=+tb;fcbP57=VtFZ*MPe=(zNO4irhYN%vpTq77&bu+J*?q5Y3{F4j>yuz zkP5nOi^r<3!`bB}A&4Gd_>MeSHsHm0-R)GbP*U@Bm5k;Gp`_qH81Z5@g)`FcpqexQ zyBcl1Ny}ca$?B-7U9IquIhm^r^;xgB@k)g!EwY3`=Fq062T9#Qfw+U2le+{WlgUdD zEUextL8+?ZaAooD`Cr?<{N=9;SJW1XE6I}130TOye{@wgLZ&qJ4AfDD&)1-yYy^4D ztev6wY0v99RvdT5AJQ4J9-CAbOYG*)3|jgYm*Eh82USjlPNC$?#Ow?>8lTpXR){YD z)3AFflGv3SU*p-9{t904vtrXW<#b-17?;7SWYO(I29?-qC**&A8&l_HEM@(f765y8J#UgFa6eVd56u8IqRSPAE=*-V^y+o4qtm}mCHjUlO_kK! zg@%#4-UM3!10%ZggFVwGf z^+PGZ{yVA}xcwv6~vQ-y0d2ZtPNXgHU}Ja~>5XlqcG|yN6zD$d2FxsQ%INr?rfX;y#|p zI$noKeCj16?;v!52P9aigGt`l^#<}%@%(7%y%)c`las5vX4mWw_jiEDrd_lp<@-1N zOC(!Sd;2X@_fE0V5-TOHkoX7AqK~{{B2~l4u=-xCnoj}^uqowAR(`*Jq(!IKydFWO ziE%e?fGz#V5SH55Xx3}D9z$XF--<bO$lGp%8K8Rg(?NG$kk1F`S?U6rpFDFGyjrfA||)a z)cV=faLxxEZruTkBJmn!%2sx zMkN=QLWs{ZOzXTh{Nn!?Kkq2&OCoPR$|LFNtgdo2u0dBxTgL)3I-0h75V4JK*f4o? zP(nFu1{E!$V2n+7Qo#m(h`aDb%Uz(>H9);%cSGhwZ+4~gSmxiqf6eayPx0d}8Sar0 z4v}olZ=Yfkk=E#Zm2EaYvZwBeWQ`WDWvhos>1ZTs`*rGJu6yT`z)$eU^)0(-ZxXOYDBHk*I$TSSk^gBv>ft8oR2^q~i={H{> z$sZ@iHKu@eTFI!w58zn)pNRYwvwA7>b}A((-!K>WLJ@YVR|Ch(O2d>UrfLN;GthIF z{@73YkUai@l|Pc+8fuQ7j+di)Ja5iDzsov!lV%{F1CcZtYAL_gh4xB#nz&?Pb}*(1 zh(5vWHi(s%m@lM{ST ziAC|`9oB2SvQwRnIX$1>_X#X7(a*IrVYJFfhb%g!`GPSP_OZ^DAnRyO`(Ixg>Mx=j zo_Ty;UnyF7#rb3XRT!K{B7LDJvCpt{Jn|L0d5n#|#b%F*Y)GzXR9? zFF(aXzv_|Xsm>HaAn9tFg_%;V2~0km>hrlV6DGH;neVQD1tZZfy78IK`*w7r&%7du z!OEfmfaB$l+9Z>19B;ksy1csO&a2DV(ObyDT`^^!#s z34|qr$6{;y)AwV9G!{J4!3H)A8Z3UgM0OUOAkELMGrQ;G<(LIyn(f)Wt~t0$L zKA-Xp8@D;y{3Io;hC_cEQ|%vUYFLfnoYxllCT8Pf_Ezah>fO`D0kjB9*Ft8OfQ9$= z!hRY3*FW{uQPeaL$oT^jaP zRr{a<(rujb>Fb+%--DmmwIk?MWhA|dyD6_uY_lSB?bo;N0XSVE zIna}t>@@slX$mncnf_%e>GMI?Dr9#4+Zvaxzl8D$$4XwHKNu>(8jfPCHO#O@b#h{A z(S)6`5g6e9R*Bdex3|$5me}MS-K@V^TQ(9 zZye}8{GfNF!_MuMThWbmb|GENrcu96oEx+wvSI@zD3iVMRN`#^8(DYDj6;`AgBI2K zql2P?e^Edbrd$G~QB^d~+aFVytPLk6@;f%8qDX-_OQZF;`Kck%oI{8u^k_4;SbBXh z`Q==t?)?YPBu%mkWO6SwJ=AVB9u>1veFfu{8B@3(2jLS9Y@Ev^i-EqVCf1#OCkA5o z^>b*lDTZ3x%{V92aH-<{05O*(Rszt_a~$T7!(P;x(3!XcKm z&ax6h%!YNZL>Kc0T^7yj>AzD5@ z)GlDrUG4GrqKsFw6Un~pfs21Z<8I#31>#R7=UFZv^b3v;gd{34v~SF}37Vw8;Ck7y zSsZ2au?rc#g>*>5oHbGFZ+*3IuYa8tyi+*QbF{y5E?{1+D!ZgWW~>=anto_)kbE?s zTz;6yLGgQg{b%^;`DhheYkz~U7{#-UIK+kp)iY&4G$(U3o1j&&y+}fMD%3k-?#Sp* ztrwfj`z1@%LFtnGW7&gM$KAEPRvp@o-fowAWq0n#b{NY!h!6=fG^DhDkKTN%p(@y4 zbU!ig?{w6-)M&n=UVY)l{h@pQ2J%I))B5AL@GGo8vK6SE<3t@|ahGaz*xp~?UONc^ zCmE;=aF3XIF=p<|!*>fBmb5xq<*utav((Pd#(U28ALWJUWj?VK3{@{*VADIi;eEJU zMB)9}qh+Ai&#NsX&g6QUZ4vNVfG>h0SO+o8xSkecU#_tLPR>M9Gvu61!!L|%$v#8F z`l&&AYjaVXh12&J2f5!i<-YS~)y1y#CsJlEWL+^!pRWJzL_=d31&nTLRopLYe}u7q z(@vG;vbNR(N11x=*qNcKQ|kS?vv21M&5v~J@bR8&K~6X&uE@TjLTO_#+P41BRvG(6 z39f4Q@Wsm9kJ|0s0m?7bLt95KnGj6d_at^Ld65j*E{L67>NNX- zS-1h_zR|94`|h{aGa{36K*ypiS8u)?2m{w4cjT|EvF}zjp5q)%AEdS$8mT%^@GiN{ zL-1s6qkct&k}fT0h&5$o{bNI*jYULW@$k^TDr`=3cyQ^*!x2uI^27Un^GlKq>zgMC zf7OaE_Dz`UJSq z55|yGI1)1z!B}8X`VyJ5Xgej8f27cD^O6Tg_d;;MUrA#jE9iULM7^)P)N!$e7Ohj zWMvC<*0SC#aA@Dmx%m04#5Cx*RE64Zvy{?$-Bs6eOU!IsXE2vnqsZB%m_UJ zd`%Bcs!1}Fkv(amGqF5bS(3Cb3eCQ2`?DBB08|(b+$!d2ki3XzLLAAyLQgnX(-u?9Uo?yXgLG!Z@+!2)(3ihFQG*U!Ywa-_X;=Y7+Bu5Ij9j?J@D;o1A-RqdX3|9W6RnfRh~pCL*zlO~QMaf+Rl{6JR0Hi9+N{@0B-JEPtvZoYr8@F`NnElbg(~A{$n|7}_tth!yKz%}VmS zMZh4_B4u+4DJhNYb_xOR21FyCsh_m46ZBHhMajrn&&q6QYPNRF*UL3ztKJTg!Z#M2* zKDPDs+lJ$(Y|C8cSAY!q@(Lc$dFYeGD|uvW#Q0MbgA!8H1_(#;drxjN_ozp9s&1(I z=c9-J(qYG|_b96JK3S-Lcp%mALAk$`eja*unw0+{MB+pfVjpKyN%o53wE4plf!ewA z2ivW3Z#i4;UWPjS0-yQL^Ty4JqQnP#gk%D06@G9c?c|W@gSHEjNBe6Vb7_&Acw2FR z->a()sHv#|-@er;D2yEKd@=jcU{fc1wLvvN#H@VeD;OB}+9^qei zbN|||HJ{MmC@HG7HZKqhH(Nd6+j~dnt^>k7tjpK z$2=)dT_TV<{0^>8xf)ElN)U(#Gqg(MWv46IaFk12wdG}kO^Rl#Xv0!q7W~N5hMUAau=~f_t{sRDX{vC%CjRvfS3c%JnPXz7%R5ANvfn4GGG}Zwx1Z;p5~~yM&`SSj&1O-zfnwGp z*;69A6U^$wSUUa_iFVHV3b*S=;iPe4<)+WmOMxnS-$v6)d>Hs`a+`CrJ$EEei43y7 z$f>#-A?s)tu%9vV%{(sirwG-~!{~`;`m6`C96Pqnxyrd?o{;1GCHg~QYv!`w1R>k? znm~rxUk+fI0ZtMvEt8alW#aftR37xm%5cSS?6|ny*vF9m@j%!Q#8^>)PC8;B@%PMQ zA}AtIA%Drug0@m>L@rinl}@vGM`G#tUW2&(*}gSgHS$%c>2DA{fL9DXE1zGq3u zU2$|giu(Dy#VS}{zRb28x4VI-8-l%4VBTZ~{7IQ{0&#L|>4H7fTpGq%5d_N*G{u_aO?j*IdW|XF7{BC0wU4>G1H$b`A??*Xcl-e z6>#(+U@d-PV5sTPlVD*TarthlH1cpA%f|<5Y*?q1NVu(~t9h^LHgJgHPLZSkKkoZ5 zai=#CHwev!`$TMfIiZoJ{4Hlewh+RFB$=69Dq28KeX9hmlji2+sK9B=5!Rb}uBTRM zKi|E~BI%(F=Q5yP>o{ZQHhfh3AOGj0Ss*gcKF-0xK?|r; z4`uHchNPh-h-H&C;O>gfO`(&|pFj6Wy+y=2G%#~njpn}}QqqNH1^b5eXZ#5_&%mK& z!p$>>_k;GvCb~A+g**USf(nqf`th>eoX<0He$tuqAOGjWfd&XTgv)Dy#r=P}wx;JC zP_I&t(Q#d*rJSI-dzsIom=y+$|yC3o1Z_R-?Xa>1m$~n<(pdzc|{Izgf}-kJIay~PFZ%^ zuasFflbX7IBpjEGgQKZ>g(mM0gU@$8St>AV7^Nac{|3a;etEd0x&GH%;|@5=84q|H zob}xpGQV;V*t7%i)Wk!4K~)EEwc>36e5-LDtE7tjK2Ew@9k^MRET*4tVCnLscJLwa zdOHOlPEp#pxs~`%VuRtKM$+9_BRC6d6>&!7OZ{T$w+_d`8ZL!2&xl~z!YS9``qfXA zo#67@9|F6!L&kqI_N=&*kfVcxu=ImHra3FZ8QMN7(#O9AgLZDaL$Ulf_+Y!Y@7~1| zN{S-zkRQN~o%bKkoc#^Bqtz|pY%5i6AMr!%v&UR=cVt)puTtAGOVfi(I2*>%LOwY zaWG&uTyL-CsZ*!opBXQFSrtY*vr9Ok1SYBv4o2x18Rec4v6(YJ z+AwDbzbPW-WPG;6N)qmEaULnXa9OP$WGDkZpKu;rSlH&ubOUdlU9+?M#>U3;n`!^} z@Bx4>kHcATxVpyJcx+|~O7#8W>RE)2ERX&Q9({m)scnX>cMyu4cPwhiOg>n&1uM|@ znU=*~C=$Q?qi!QqZvR^H-uKyc04sq#XS$mew98#lRfSLTUvB$c?CcD$67C5Z92&xN zuv-zkIn_QCSHTJ;B=?YoxpwEc^GMqXVM`zlx-Q zO~a?xp?uQ%3T{vFXr{h07sj2`ldYLnD#NTF)^eAZho^TYQmzoLFu~X23;VA${-Gv- z`6H}xR)mXIj&_TVa+CLeauW!Ogo|e82$kpMx|K#J8fr#*UX1t^hyTy6X#D(J2VRfA z53(K~SP`sQljE@c`=Wz3>!!*7`_R7Y|G(c~FL~={0?!Y6DYI-b&lkhm_j|v80iFkx z4a~8fuYu>osA#IIyIYsPyD~F<-pqTB)&~;)tp%oaZFBSRy8ZTrcP*dKS$ykOl&g{; zD^uf^zMJpZKF_#U`~B{&hi%d&K-YOsI%oSm2DrWPwRillN&7zb=AU@Jc%D(f-o4?X zqN0I+Pc!^GdMidj0K5yW?(eT7rNE7audlBU2UcDE-&d|$^{TAeGirNY?0Sg9Cw&2C z{`l{A%jd^cJYq(VPG8%W~NI72sybnKJuwd7s^X z1x(+ofMfKkOj}j<0{3wM z>*I3SnjMRuo|^g$c*;*|-AnIyBjD1;yldC4&03Z@%luA3^Utm8Yjq2M1KZno&zwAI z2wFKgD}VprZI$P3zh7Wm#?D+D^M3#TI%5Hj8T)~I&Ci18c6@z(Z~hX^Nng zj&z`+?{9CXgEBH`{2Dk5;`bR?P#b%Dd7XJ*_kA}^vF-Oe$sk3*{P92a-MziH{{Vv{ ze_rLYna|SgKF`eg%eC$Dl`AQrvhNOX>x3{14l$#PBU| zLI3Mxz0&{xm%Y7l%u#{k@CrRfbH`044$m$uoe<30`Feqf*CK`q$1b{-pH#K;{$Ksm zYcXfo{IKZnzxV#yxBgm`?fmn**1bL(T)%Dqn=dx!qpoi~R=DQ&+pxn1S2pRC${e@% z`8V@h?f` method. The following sequence (done in a python session) reads in stored data (from the compressible Sedov problem) and plots data falling on a line in the x -direction through the y-center of the domain (note: this will include -the ghost cells). +direction through the y-center of the domain. The return value of +``read`` is a ``Simulation`` object. + .. code-block:: python import matplotlib.pyplot as plt import pyro.util.io_pyro as io - sim = io.read("sedov_unsplit_0000.h5") + + sim = io.read("sedov_unsplit_0290.h5") dens = sim.cc_data.get_var("density") - plt.plot(dens.g.x, dens[:,dens.g.ny//2]) - plt.show() + + fig, ax = plt.subplots() + ax.plot(dens.g.x, dens[:,dens.g.qy//2]) + ax.grid() .. image:: manual_plot.png :align: center -Note: this includes the ghost cells, by default, seen as the small -regions of zeros on the left and right. +.. note:: + + This includes the ghost cells, by default, seen as the small + regions of zeros on the left and right. The total number of cells, + including ghost cells in the y-direction is ``qy``, which is why + we use that in our slice. + +If we wanted to exclude the ghost cells, then we could use the ``.v()`` method +on the density array to exclude the ghost cells, and then manually index ``g.x`` +to just include the valid part of the domain: + +.. code:: python + + ax.plot(dens.g.x[g.ilo:g.ihi+1], dens.v()[:, dens.g.ny//2]) + +.. note:: + + In this case, we are using ``ny`` since that is the width of the domain + excluding ghost cells. From 81873bc4f77369978aecf83a3771614d446273f5 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 13 Dec 2024 17:53:28 -0500 Subject: [PATCH 2/6] add a demonstration of doing a lateral average (#305) --- docs/source/index.rst | 1 + docs/source/rt_average.ipynb | 155 +++++++++++++++++++++++++++++++++++ 2 files changed, 156 insertions(+) create mode 100644 docs/source/rt_average.ipynb diff --git a/docs/source/index.rst b/docs/source/index.rst index 688308baa..5dfa8a232 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -72,6 +72,7 @@ new ideas. compressible-rt-compare.ipynb adding_a_problem_jupyter.ipynb + rt_average.ipynb advection-error.ipynb compressible-convergence.ipynb diff --git a/docs/source/rt_average.ipynb b/docs/source/rt_average.ipynb new file mode 100644 index 000000000..14fdaa201 --- /dev/null +++ b/docs/source/rt_average.ipynb @@ -0,0 +1,155 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e405b84e-09c0-4c8f-b178-02102b2782be", + "metadata": {}, + "source": [ + "# Horizontal Averages of Rayleigh-Taylor" + ] + }, + { + "cell_type": "markdown", + "id": "6395d3eb-9572-454b-823a-b5c2f8c4267d", + "metadata": {}, + "source": [ + "Here we show how to do a lateral / horizonal average of the compressible Rayleigh-Taylor\n", + "simulation to see what the average vertical profile of the mixing looks like." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "571afb49-3a5d-49d1-a7e3-8e51bb6c7d11", + "metadata": {}, + "outputs": [], + "source": [ + "from pyro import Pyro" + ] + }, + { + "cell_type": "markdown", + "id": "97989792-3858-43f4-8e59-dca0a2686777", + "metadata": {}, + "source": [ + "First we'll run the `rt` problem using the problem defaults." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7e5a5084-c403-41d4-aec7-0bd48f0be310", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/zingale/development/pyro2/pyro/compressible/problems/rt.py:71: RuntimeWarning: invalid value encountered in divide\n", + " 0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]\n" + ] + } + ], + "source": [ + "p = Pyro(\"compressible\")\n", + "p.initialize_problem(\"rt\")\n", + "p.run_sim()" + ] + }, + { + "cell_type": "markdown", + "id": "9ab0b191-8a7a-409f-8773-ba0f1b1dd49e", + "metadata": {}, + "source": [ + "Now we'll get the density variable and use `np.average()` to compute the average in the x-direction.\n", + "Note that we operate only on the valid data (`dens.v()`)." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "bdbb2ac1-54ae-4f5a-b96a-fdc3eb602221", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "1a0909a9-5a61-4606-9e87-83d4559213c6", + "metadata": {}, + "outputs": [], + "source": [ + "dens = p.get_var(\"density\")\n", + "dens_avg = np.average(dens.v(), axis=0)" + ] + }, + { + "cell_type": "markdown", + "id": "81ae06fb-6125-4a67-a4f1-fbd157f79b03", + "metadata": {}, + "source": [ + "Finally, we can plot the profile." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "950658f0-5c24-4392-94a9-17dbfc769a63", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'average density')" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMe0lEQVR4nO3deXwU9f0/8NfsbnZzX4RcJJBAuI8QzoIHiBwCxaIVVGyB4lEVK15Y6E+l1CPiVyu21Vq0hWK9L7QeSAQJilAECVfkDiSE3NdmN8me8/tjDwgkkE12d3ZnXs/HIw+yk5ndd4awefE5BVEURRARERHJhErqAoiIiIi8ieGGiIiIZIXhhoiIiGSF4YaIiIhkheGGiIiIZIXhhoiIiGSF4YaIiIhkRSN1Af5mt9tx9uxZREVFQRAEqcshIiKiDhBFEY2NjUhNTYVKdem2GcWFm7NnzyI9PV3qMoiIiKgTSkpKkJaWdslzFBduoqKiAOfNiY6OlrocIiIi6gC9Xo/09HT37/FLUVy4cXVFRUdHM9wQEREFmY4MKeGAYiIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhXFbZxJRESkFKIowmYXYbGJMNvssNjssNpEWO12qFWC40MQoFGpoFIBapXjc60muNs+GG6IiCjgiKIIk9UOk8UOs83uPi4IgOD+XIBKALQaFbRqFTTq4P6F7CKKIprMNhhMVjS2WGEwWWFoscJgsqCx5bxjrb5ucT92HWs221rdO09oVALCtWqEazUI16kRodUgTKtGhFaN6LAQRIeGICYsBNFhGkSHhrR5LC5C6/V70+H6JXtlIiKSPVEUUW0wo7S+GZX6FtQ1mVFjNKPOaEat0YJaowm1TRY0NlvQbLGhxWJDi8WOFqsNoujZa6kEQKdRO8KORgWd80+tWgVdiBo6tQq6kHPHdRr1eZ+r3NfqXB8haue1KoSoVRCcgQo4P2A5PhzHBFhsdlhsIiw2O8xWu7u1xGy1o9lig9FkhdHs/NNkQ5PZ2uqYweR4bPfwe/eERiVApRJgt4uwiWKb99lqF6FvsULfYu3Ua8SGh6DgialdL7aTGG6IiKhLRFFEWUMLTlQZcLzSgBNVBhRVG3G2vgWl9c0wWzvXeuApuwg0W2xottj88nq+plYJiNRpEKnTICr0vD9DQ1odi9RpEBmqQZTzT9fXwrQaaNUqhKgFhKgdAc0VbM7nCjk2uwi7KMJiFdFksaLJbEOTM4A1mW0wOoNYY4sV+maLI/w0W9DQbIG+xQJ9s9X5pwXRoSGS3Tcw3BARkacq9C3YW1yPgpJ67Cupx4HSBhhM7f8PXyUASdGhSIoORbcILeIitIiP0CIuXIv4iBDER+gQHero9ggNUSNUo0ZoiAqhWsfnIWrB3WICZ5hysdkdY0nMVseHyf1hcz9u/eeFx22tznE9bnXMYoPFZofofn3nn+c9EJ2fapxBQqdROQOFAK1GjRCVgDCtGpE6DcK1GkTo1IjQaRB+3rFInaMLKCpUgyhdCEJDVK2+b19RqQSoICBE7TygBWLQtXBi82XTUwcw3BAR0SWZrDZsPVKFrUeqsO1oFUrrmy86R6MS0KtbOLISI9GneyR6d49EWlwYesSGITkmFCFeHA9z/i98jVqARq1CuHTDO6gNapXvQ9mlSBpucnNz8dFHH+Hw4cMICwvD+PHjsWrVKvTv3/+S173//vt4/PHHcerUKfTt2xerVq3CjBkz/FY3EZESlNQ24a1dxXjvhxLUGM3u4yoB6JcUhZyesRieHovs9Fj06R7p1QBD1BWShpv8/HwsXrwYo0ePhtVqxR/+8AdMnToVhYWFiIiIaPOa77//Hrfeeityc3Px85//HG+99RZmz56NH3/8EUOGDPH790BEJCc2u4htR6vwn52nseVIpbsLJjk6FNcNScaEft0xJjMeETo2/FPgEkTR0/HovlNVVYXExETk5+fj6quvbvOcm2++GUajEZ999pn72M9+9jMMHz4cr7766kXnm0wmmEwm92O9Xo/09HQ0NDQgOjraR98JEVFwaWiy4O0fivHm/06jpPZct9NVfRNw29hemDwwUTZTrSk46fV6xMTEdOj3d0BF74aGBgBAfHx8u+fs2LEDDz30UKtj06ZNw4YNG9o8Pzc3FytXrvRypURE8tDQZMHr353Euu2n0OgcFBwdqsGcUem4bWxP9O4eKXWJRB4LmHBjt9vxwAMP4Iorrrhk91J5eTmSkpJaHUtKSkJ5eXmb5y9fvrxVGHK13BAREfCbdbvwY3E9AKB/UhRuvyoTs4alIkyrvuy1RIEqYMLN4sWLcfDgQXz33XdefV6dTgedTufV5yQikoPTNUb8WFwPjUrAX2/NwbTByRetgUIUjAIi3Nx333347LPPsG3bNqSlpV3y3OTkZFRUVLQ6VlFRgeTkZB9XSUQkL9uOVQMARvaKw/ShKVKXQ+Q1ko4OE0UR9913Hz7++GNs2bIFmZmZl71m3Lhx2Lx5c6tjeXl5GDdunA8rJSKSn21HqwAAV/frLnUpRF4lacvN4sWL8dZbb+GTTz5BVFSUe9xMTEwMwsLCAADz589Hjx49kJubCwBYsmQJJkyYgBdeeAEzZ87EO++8g927d2PNmjVSfitEREHFYrNjx4kaAMDVfRluSF4kbbn5+9//joaGBkycOBEpKSnuj3fffdd9TnFxMcrKytyPx48fj7feegtr1qxBdnY2PvjgA2zYsIFr3BAReWBvcT0MJiviI7QYnMplMUheJG256cgSO1u3br3o2Jw5czBnzhwfVUVEJH/fHnN0SV2ZlcBBxCQ7XJGJiEiBXONtruqbIHUpRF7HcENEpDB1RjP2lzoWTb2K421IhhhuiIgUZn9pA0QR6N09AskxoVKXQ+R1DDdERApT59zhOymKwYbkieGGiEhhGpotAICYsBCpSyHyCYYbIiKF0TPckMwx3BARKYy75Sac4YbkieGGiEhh2C1FcsdwQ0SkMAw3JHcMN0RECsNwQ3LHcENEpDAMNyR3DDdERArD2VIkdww3REQKw5YbkjuGGyIiBbHY7DCabQDDDckYww0RkYK4uqQAIJrhhmSK4YaISEFcXVJROg3UKkHqcoh8guGGiEhBXOGGrTYkZww3REQKwsHEpAQMN0RECsJwQ0rAcENEpCBc44aUgOGGiEhB2HJDSsBwQ0SkIO5wE85wQ/LFcENEpCBsuSElYLghIlIQTgUnJWC4ISJSELbckBIw3BARKUhDsxVguCGZY7ghIlIQTgUnJWC4ISJSEHZLkRIw3BARKYTVZofBxG4pkj+GGyIihdC3WN2fR4dqJK2FyJcYboiIFMLVJRWp00Cj5ts/yRd/uomIFILjbUgpGG6IiBSCC/iRUjDcEBEpxLmWG463IXljuCEiUgh2S5FSMNwQESkEF/AjpWC4ISJSCLbckFIw3BARKURDE8MNKQPDDRGRQuhbOFuKlIHhhohIIRqdKxRHcXVikjmGGyIihWh0ttxE6dhyQ/LGcENEpBBsuSGlYLghIlIIvTvcsOWG5I3hhohIIdzdUmy5IZljuCEiUgCz1Q6T1Q4AiGbLDckcww0RkQK4Wm0AIJItNyRzDDdERArgGkwcoVVDrRKkLofIpxhuiIgUoJGDiUlBGG6IiBSg0eTolmKXFCkBww0RkQJwjRtSEoYbIiIFYLcUKQnDDRGRAnCNG1IShhsiIgVwtdxEM9yQAjDcEBEpwLmWG3ZLkfwx3BARKYB7zI2OLTckf5KGm23btmHWrFlITU2FIAjYsGHDZa958803kZ2djfDwcKSkpGDRokWoqanxS71ERMGKs6VISSQNN0ajEdnZ2Xj55Zc7dP727dsxf/583H777Th06BDef/997Nq1C3feeafPayUiCmZ6dkuRgkga4adPn47p06d3+PwdO3YgIyMD999/PwAgMzMTv/3tb7Fq1ap2rzGZTDCZTO7Her2+i1UTEQUfttyQkgTVmJtx48ahpKQEX3zxBURRREVFBT744APMmDGj3Wtyc3MRExPj/khPT/drzUREgYADiklJgircXHHFFXjzzTdx8803Q6vVIjk5GTExMZfs1lq+fDkaGhrcHyUlJX6tmYgoELDlhpQkqMJNYWEhlixZgieeeAJ79uzBxo0bcerUKdx9993tXqPT6RAdHd3qg4hIac6tc8OWG5K/oIrwubm5uOKKK7B06VIAwLBhwxAREYGrrroKTz31FFJSUqQukYgo4FhsdjRbbABbbkghgqrlpqmpCSpV65LVajUAQBRFiaoiIgpsBmerDbgrOCmEpOHGYDCgoKAABQUFAICioiIUFBSguLgYcI6XmT9/vvv8WbNm4aOPPsLf//53nDx5Etu3b8f999+PMWPGIDU1VbLvg4gokLm6pMJC1AhRB9X/aYk6RdIIv3v3blxzzTXuxw899BAAYMGCBVi3bh3KysrcQQcAFi5ciMbGRvztb3/Dww8/jNjYWEyaNOmSU8GJiJROz00zSWEEUWH9OXq9HjExMWhoaODgYiJShB0nanDrazvRu3sEtjw8UepyiDrFk9/fbJ8kIpI5rnFDSsNwQ0Qkc+emgbNbipSB4YaISOYMJi7gR8rCcENEJHPubikdu6VIGRhuiIhkjlsvkNIw3BARyZzeHW7YckPKwHBDRCRzjVznhhSG4YaISObYLUVKw3BDRCRzXOeGlIbhhohI5rjODSkNww0Rkcw1ckAxKQzDDRGRzHFAMSkNww0RkYxZbHYYzTYAQHQYW25IGRhuiIhkrMZgBgCoVQJiGW5IIRhuiIhkrNpgAgDER2ihUglSl0PkFww3REQyVuUMNwmROqlLIfIbhhsiIhmrbnSEm+5RDDekHAw3REQyVu0cc5MQqZW6FCK/YbghIpIx15ib7uyWIgVhuCEikrFqjrkhBWK4ISKSMXe4iWK3FCkHww0RkYxVNbLlhpSH4YaISMbODShmuCHlYLghIpIpq82OuiaGG1IehhsiIpmqNZohioBKcKxQTKQUDDdERDJV5d56QQc1t14gBWG4ISKSKS7gR0rFcEMkExv2lmLkk3nYdKhc6lIoQHDrBVIqhhsiGThdY8QfPj6AGqMZT3xyCE1mq9QlUQDgAn6kVAw3REHOZhex9P39aDLbAADl+ha8tq1I6rIoAJwLN+yWImVhuCEKcmu3F2HXqVpEaNX4/XUDAACv5p9Ahb5F6tJIYlzAj5SK4YYoiFXqW/DCpqMAgMd+Pgh3T+iNET1j0Wyx4c/O46RcXMCPlIrhhiiIvfj1MTRbbBjRMxa3jE6HIAh4ZFp/AMCmQg4sVrpz+0ox3JCyeBxuJkyYgPXr16O5udk3FRFRhxyvNOC93SUAgOUzBkIQHOuYDE6JAQDUNVnQ7ByHQ8rkCjfd2XJDCuNxuMnJycEjjzyC5ORk3Hnnndi5c6dvKiOiS3pu42HY7CImD0zC6Ix49/HoMA3CtWoAQFkD/xOiVDa7iFqjs1uKO4KTwngcblavXo2zZ89i7dq1qKysxNVXX41Bgwbh+eefR0VFhW+qJKJWfiyuw6bCCqgE4PfX9W/1NUEQkBobBgAoa+CgYqWqNZphFwFBAOLDGW5IWTo15kaj0eDGG2/EJ598gjNnzmDevHl4/PHHkZ6ejtmzZ2PLli3er5SI3P66+RgA4Jcj0tA3Keqir6fEhAIASuvZcqNUri6p+HAtNGoOryRl6dJP/K5du7BixQq88MILSExMxPLly5GQkICf//zneOSRR7xXJRG5HSxtwDdHqqASgHuvyWrznNQYZ8tNPVtulIoL+JGSaTy9oLKyEm+88QbWrl2LY8eOYdasWXj77bcxbdo094DGhQsX4rrrrsPzzz/vi5qJFO1vW44DAGZlpyIzIaLNc1JiHS03HHOjXJV610wpdkmR8ngcbtLS0tCnTx8sWrQICxcuRPfu3S86Z9iwYRg9erS3aiQip6MVjdjo3DvqvnZabXBey81ZjrlRrJK6JgBAely41KUQ+Z3H4Wbz5s246qqrLnlOdHQ0vvnmm67URURteHXrCQDA9CHJbY61cXEPKOaYG8UqrnGGm3iGG1Iej8fcrFixAvX19Rcd1+v1mDRpkrfqIqILVDWa8N/9ZwEAv53Q55LnurqlztY3QxRFv9RHgeV0rSPc9OrGcEPK43G4yc/Ph9lsvuh4S0sLvv32W2/VRUQXeHtXMSw2EcPTYzE8PfaS57q6pYxmG/Qt3CFciYqd4aYnW25IgTrcLbV//34AgCiKKCwsRHn5uaXdbTYbNm7ciB49evimSiKFs9jsePN/pwEAC8dnXPb8MK0aseEhqG+yoKyhGTFhIX6okgJFk9nq3jSzV3zbg86J5KzD4Wb48OEQBAGCILTZ/RQWFoa//vWv3q6PiABsPFiOCr0J3aN0mDE0pUPXpMSEOcJNfQsGJEf7ukQKICW1jrFW0aEaxIQz2JLydDjcFBUVQRRF9O7dG7t27Wo1S0qr1SIxMRFqtdpXdRIp2r+/PwUAmDemJ7SajvUm94gNxU9lepzldHDFOV1jBAD06sZWG1KmDoebXr16AQDsdrsv6yGiCxyvNGD36TpoVAJuG9uzw9eluKaDc8aU4nC8DSldh8LNp59+iunTpyMkJASffvrpJc+9/vrrvVUbEQH4pKAUADChX3ckRod2+Dr3Qn5cpVhx3OGGM6VIoToUbmbPno3y8nIkJiZi9uzZ7Z4nCAJsNps36yNSNFEU8fFeR7iZnePZgP1zC/mx5UZp2HJDStehcHN+VxS7pYj8Z/fpOpypa0akToPJA5M8upY7gyuXawG/Xgw3pFBe2Sq2rUX9iKjrXK021w1JRpjWswH7rp3ByxpaYLdzIT+lsNlFnKlztNZxdWJSKo/DzapVq/Duu++6H8+ZMwfx8fHo0aMH9u3b5+36iBTLZLXh8/1lAIAbPeySAoDkmFAIAmC22lFjvHjhTZKncn0LzDY7NCrB3XpHpDQeh5tXX30V6enpAIC8vDx8/fXX2LhxI6ZPn46lS5f6okYiRco/UoWGZguSo0Mxtnc3j68PUavQPVIHcHdwRXF1SaXFhUGtEqQuh0gSHm+cWV5e7g43n332GebOnYupU6ciIyMDY8eO9UWNRIr01aEKAMCMoSmd/iWVGK1DZaMJNQa23ChFca1jjZueXOOGFMzjlpu4uDiUlJQAADZu3IjJkycDzlkdnClF5B1Wmx2bDzvCzbTBng0kPl9cuBYAUMtuKcU4N1OKXVKkXB6HmxtvvBHz5s3DlClTUFNTg+nTpwMA9u7di6ysLI+ea9u2bZg1axZSU1MhCAI2bNhw2WtMJhP+3//7f+jVqxd0Oh0yMjLwr3/9y9Nvgyig/XCqDvVNFsSFh2Bkr7hOP098BMON0px2z5Riyw0pl8fdUi+++CIyMjJQUlKC5557DpGRkQCAsrIy3HvvvR49l9FoRHZ2NhYtWoQbb7yxQ9fMnTsXFRUV+Oc//4msrCyUlZVxejrJTl6ho9Xm2oFJ0Kg7P6nRHW6aGG6U4kSVa+sFzpQi5fI43ISEhOCRRx656PiDDz7o8YtPnz7d3fLTERs3bkR+fj5OnjyJ+Ph4AEBGxqV3SDaZTDCZTO7Her3e4zqJ/EkURWwqLAcATB3U+S4pAIh3dkvVseVGESw2O05UGgCAm6WSonkcbgDg2LFj+Oabb1BZWXlRq8kTTzzhrdou8umnn2LUqFF47rnn8MYbbyAiIgLXX389nnzySYSFtd2/nJubi5UrV/qsJiJv+6msEWfqmhEaosJVfbt34Ir2xbFbSlFO1xhhttkRrlUjLY5jbki5PA43r732Gu655x4kJCQgOTkZgnBuFocgCD4NNydPnsR3332H0NBQfPzxx6iursa9996LmpoarF27ts1rli9fjoceesj9WK/Xu2d7EQUiV6vNVX27e7xw34Vc3VJ17JZShCPljlabvklRUHEaOCmYx+HmqaeewtNPP43f//73vqnoEux2OwRBwJtvvomYmBgAwJ///GfcdNNNeOWVV9psvdHpdNDpdH6vlaizvjlcCQCY0sUuKXC2lOIcKXd0u/dPipS6FCJJeTxSsa6uDnPmzPFNNZeRkpKCHj16uIMNAAwcOBCiKOLMmTOS1ETkTXVGM/aXNgDOXcC76lzLjaXLz0WB70hFIwCgP8fbkMJ5HG7mzJmDTZs2+aaay7jiiitw9uxZGAwG97GjR49CpVIhLS1NkpqIvOn7EzUQRaBfUiSSokO7/HxxESEAgPomM2zcX0r2jlY43hv7J0VJXQqRpDzulsrKysLjjz+OnTt3YujQoQgJCWn19fvvv7/Dz2UwGHD8+HH346KiIhQUFCA+Ph49e/bE8uXLUVpaivXr1wMA5s2bhyeffBK/+c1vsHLlSlRXV2Pp0qVYtGhRuwOKiYLJd8erAABXZnW91QbndUvZRUDfbHEPMCb5abHYcKrGMQ28fzLDDSmbx+FmzZo1iIyMRH5+PvLz81t9TRAEj8LN7t27cc0117gfuwb+LliwAOvWrUNZWRmKi4vdX4+MjEReXh5+97vfYdSoUejWrRvmzp2Lp556ytNvgyjgiKKIb49VAwCu6pvglecMUasQFapBY4sVtU1mhhsZO1ZhgCg6uiITIvn3TMrmcbgpKiry2otPnDgRoth+U/m6desuOjZgwADk5eV5rQaiQHG6pgln6poRohYwtne81563W4QWjS1Wx1o33mkQogDkGm/TLymy1SxWIiXq9NKnZrMZR44cgdVq9W5FRAr17XFHq82InnEI13ZqCao2uVprajhjStaOOsMNF+8j6kS4aWpqwu23347w8HAMHjzY3W30u9/9Ds8++6wvaiRShG+POsbbeKtLyoWrFCvD4XJXyw3H2xB5HG6WL1+Offv2YevWrQgNPTebY/LkyXj33Xe9XR+RIlhtduw4WQMAuLKLqxJfKI77SynC0XLXNHCucUPkcdv3hg0b8O677+JnP/tZq37dwYMH48SJE96uj0gRfiprRGOLFVGhGgztEdOBKzrOvdYNW25kq77JjHJ9C+BcnZhI6TxuuamqqkJiYuJFx41GIwexEXXSrlO1AIBRveKg9vKy+edWKeZCfnJ1wLnwY69u4YgODbns+URy53G4GTVqFD7//HP3Y1egef311zFu3DjvVkekELud4WZ0pvdmSbnEOxfy4/5S8rX/jCPcDEuLlboUooDgcbfUM888g+nTp6OwsBBWqxUvvfQSCgsL8f3331+07g0RXZ4oivjBFW4yvB9uuL+U/B1whRsvd2kSBSuPW26uvPJKFBQUwGq1YujQodi0aRMSExOxY8cOjBw50jdVEsnYqZomVBvM0GpUGJbm/V9O3Blc/vafqQcADPXBzw9RMOrUYhp9+vTBa6+95v1qiBTohyJHq012Wgx0GrXXn98VbthyI09VjSacbWiBIABD2HJDBHQ03Oj1+g4/YXQ0F5Ai8oQvu6RwXrhpbLHCYrMjRN3ptTspAB10Dibu0z0SkTrvLf5IFMw69C8hNja2wzOhbDZbV2siUhRfh5vo0BCoBMfmmXVNZiRGdX23cQoc+5xdUhxvQ3ROh8LNN9984/781KlTWLZsGRYuXOieHbVjxw78+9//Rm5uru8qJZKhysYWnKppgiAAI3rF+eQ1VCoBceFa1BjNqDUy3MiNazAxx9sQndOhcDNhwgT353/605/w5z//Gbfeeqv72PXXX4+hQ4dizZo1WLBggW8qJZKhPafqAAD9k6IQE+a79UniIs6FG5IPURSxv5TTwIku5HHn+44dOzBq1KiLjo8aNQq7du3yVl1EirC3xNGlMNJHrTYu5/aX4kJ+clKub0FVowlqlYBBKRzvSOTicbhJT09vc6bU66+/jvT0dG/VRaQIrim82em+/V93nHMhP+4vJS+uxfv6JkYiTOv9mXZEwcrjofUvvvgifvnLX+LLL7/E2LFjAQC7du3CsWPH8OGHH/qiRiJZsttFHCx1zET0xfo25+P+UvK0t9gZjtklRdSKxy03M2bMwLFjx3D99dejtrYWtbW1mDVrFo4ePYoZM2b4pkoiGTpZbYTBZEVoiApZ3X27kzPXupGnPacdM+1GZvi2W5Mo2HRqUYS0tDQ8/fTT3q+GSEEOlDr+1z0kNQYaH689wy0Y5MdktWGfs1vKV8sIEAUrruZFJJF9Jf6bwpsQqQMA1BhNPn8t8o+DpQ0wW+3oFqFFRrdwqcshCigMN0QSOeCewuv7cOPqlqoxsOVGLn5wLiMwKiOuw4usEikFww2RBKw2Ow6d9d/6JN0ineGG3VKysdsVbnqxS4roQgw3RBI4VmlAi8WOSJ0Gmd0ifP563SIc3VJ1RjPsdtHnr0e+JYqiezDxKA4mJrpIp8KN1WrF119/jX/84x9obGwEAJw9exYGg8Hb9RHJkmvJ/CE9oqFS+b5LwdUtZbWL0LdwIb9gd6LKiLomC3QaFQanctsFogt5PFvq9OnTuO6661BcXAyTyYQpU6YgKioKq1atgslkwquvvuqbSolkxLXZob/WJ9FqVIgK1aCxxYpqgxmxztlTFJx2OzdbHZ4eC62GDfBEF/L4X8WSJUswatQo1NXVISwszH38hhtuwObNm71dH5EsHSx1tdz473/drhlTnA4e/HafPjeYmIgu5nHLzbfffovvv/8eWm3r//llZGSgtLTUm7URyZLNLuJIhaM7d3Cq//YDio/QoqjaiBoDp4MHu11FrvE2HExM1BaPW27sdjtsNttFx8+cOYOoqChv1UUkWyW1TWix2KHTqNDLD4OJXbpFcMaUHJTUNqG4tglqlcDF+4ja4XG4mTp1KlavXu1+LAgCDAYDVqxYwe0XiDrgcLmj1aZvUiTUfhhM7OKeDs61boLajhM1AIDstBhE6jq1yDyR7Hn8L+OFF17AtGnTMGjQILS0tGDevHk4duwYEhIS8Pbbb/umSiIZOeIMN/2T/NclhfOmg9dyleKg9v2JagDAFVkJUpdCFLA8DjdpaWnYt28f3nnnHezfvx8GgwG33347brvttlYDjImobUed4236J/t2s8wLuaaDV7NbKmiJoojtzpabcX26SV0OUcDqVJumRqPBr371K+9XQ6QAh8v1AID+yX5uuXF2S9WyWyponagyoKrRBJ1GhRE9OVOKqD0eh5tPP/20zeOCICA0NBRZWVnIzMz0Rm1EstNiseFUTRMAYECyfwfgu7qluHlm8Np+3NFqMyojDqEhaqnLIQpYHoeb2bNnQxAEiGLrJdxdxwRBwJVXXokNGzYgLo7/syA634kqA2x2ETFhIUiM0vn1tTmgOPi5xtuM78PxNkSX4vFsqby8PIwePRp5eXloaGhAQ0MD8vLyMHbsWHz22WfYtm0bampq8Mgjj/imYqIg5h5MnBzl952cXeGmrskMG/eXCjo2u4idJx3r24zneBuiS/K45WbJkiVYs2YNxo8f7z527bXXIjQ0FHfddRcOHTqE1atXY9GiRd6ulSjoucKNv7ukACDOueWCXQTqm8zoFunfliPqmoOlDWhotiBSp8FQP65sTRSMPG65OXHiBKKjLx4IGR0djZMnTwIA+vbti+rqau9USCQjrjVu+iX5P9yEqFWIDQ8BuAVDUNp6pAoAcGVWAjRq7idFdCke/wsZOXIkli5diqqqKvexqqoqPProoxg9ejQA4NixY0hPT/dupUQy4JoGLkXLDc6fDs5xN0HnmyOVAICJ/btLXQpRwPM43Pzzn/9EUVER0tLSkJWVhaysLKSlpeHUqVN4/fXXAQAGgwGPPfaYL+olCloNTRaUNbQAAPpJFG4SIrh5ZjCqNZrdO8lP7J8odTlEAc/jMTf9+/dHYWEhNm3ahKNHj7qPTZkyBSqVIyvNnj3b+5USBbnjVQYAQHJ0KKJDQySpId69vxSngweTbUerIIqOFr/kmFCpyyEKeJ1axE+lUuG6667Ddddd5/2KiGTqVLURAJCZ4L/NMi/kmjHFbqngstXZJXXNALbaEHVEp8KN0WhEfn4+iouLYTa3fpO8//77vVUbkaycrnGEmwwpw42z5Yb7SwUPm11E/lHHGMeJ/TjehqgjPA43e/fuxYwZM9DU1ASj0Yj4+HhUV1cjPDwciYmJDDdE7Shyrkyc0S1cshpc07+5kF/w2HemHnVNFkSFajCiFxdGJeoIjwcUP/jgg5g1axbq6uoQFhaGnTt34vTp0xg5ciSef/5531RJJAOubikpW27OjblhuAkWWw87uqSu6puAEE4BJ+oQj/+lFBQU4OGHH4ZKpYJarYbJZEJ6ejqee+45/OEPf/BNlURBThRFnKoJnDE3NQZ2SwWLTYUVAIBrByRJXQpR0PA43ISEhLhnRSUmJqK4uBgAEBMTg5KSEu9XSCQDtUYzGlusAICe8RJ2S7k3z2TLTTA4XWPE4fJGqFUCrh3IwcREHeXxmJucnBz88MMP6Nu3LyZMmIAnnngC1dXVeOONNzBkyBDfVEkU5Fw7gafGhEq6m7Nrs876JgtaLDbuLB3gNh1ytNqMzYxHrHP7DCK6PI9bbp555hmkpKQAAJ5++mnExcXhnnvuQVVVFdasWeOLGomCXiCMtwGA2PAQRGgdgaa0vlnSWujyNhWWAwCmDU6WuhSioOJRy40oikhMTHS30CQmJmLjxo2+qo1INlzjbXp1kzbcCIKAHnFhOFphQGldM/p0j5S0HmpftcGE3afrAABTBnG8DZEnPGq5EUURWVlZHFtD5CFXt1RmgnTjbVx6xIYBAM7UseUmkH1dWAFRBIb2iEGq8++MiDrGo3CjUqnQt29f1NTU+K4iIhlyd0tJ3HIDAD3iHL8oS+ubpC6FLsE1S2raYLbaEHnK4zE3zz77LJYuXYqDBw/6piIimRFFMWDG3ABAWpyj9aiULTcBq6HZgu+OVQMApnK8DZHHPJ4tNX/+fDQ1NSE7OxtarRZhYa2bS2tra71ZH1HQqzWa0WiyQhCknQbu4uqW4oDiwJVXWAGzzY6+iZHolyTNDvJEwczjcLN69WrfVEIkU67BxKkxYQEx9drVLcUxN4Hrs/1nAQA/H5YqdSlEQcnjcLNgwQLfVEIkU6eqHWNbekm4p9T50pzhpkLfAovNziX9A0yd0ezukpo5LEXqcoiCUqfe1U6cOIHHHnsMt956KyorHfuefPnllzh06JBHz7Nt2zbMmjULqampEAQBGzZs6PC127dvh0ajwfDhwz2un8ifTgfINHCXhAgdtBoV7CJQ3tAidTl0gU2F5bDaRQxIjkJWIqfqE3WGx+EmPz8fQ4cOxf/+9z989NFHMBgMAIB9+/ZhxYoVHj2X0WhEdnY2Xn75ZY+uq6+vx/z583Httdd6dB2RFM46A4SrxURqKpXgHndTUscZU4Hms/1lAIBZ2eySIuosj8PNsmXL8NRTTyEvLw9a7bnlwCdNmoSdO3d69FzTp0/HU089hRtuuMGj6+6++27MmzcP48aN8+g6Iim4WkdSYkKlLsXNPaiY424CSo3BhO9POJbamDmUXVJEneVxuDlw4ECbYSQxMRHV1dXeqqtda9euxcmTJzvcSmQymaDX61t9EPlTWYMjQCQHULhJi+OMqUD0331nYbOLGNojJiCWDSAKVh6Hm9jYWJSVlV10fO/evejRo4e36mrTsWPHsGzZMvznP/+BRtOxsdC5ubmIiYlxf6Snp/u0RqLziaKIMnfLTWB0S4GrFAesD38sBQDcOMK376VEcudxuLnlllvw+9//HuXl5RAEAXa7Hdu3b8cjjzyC+fPn+6ZKADabDfPmzcPKlSvRr1+/Dl+3fPlyNDQ0uD+4dQT5U6PJiiazDQCQHB04LTfuVYoZbgLG0YpGHChtgEYl4HqOtyHqEo+ngj/zzDNYvHgx0tPTYbPZMGjQIHfweOyxx3xTJYDGxkbs3r0be/fuxX333QcAsNvtEEURGo0GmzZtwqRJky66TqfTQafT+awuoktxjbeJDQ9BmFb6NW5c3KsUs1sqYHy45wwA4JoBiegWyfcsoq7wONxotVq89tprePzxx3Hw4EEYDAbk5OSgb9++vqnQKTo6GgcOHGh17JVXXsGWLVvwwQcfIDMz06evT9QZri6pQGq1wXktN2UNzbDZRahVgtQlKZrVZsfHex1dUr8ckSZ1OURBz+Nw89133+HKK69Ez5490bNnzy69uMFgwPHjx92Pi4qKUFBQgPj4ePTs2RPLly9HaWkp1q9fD5VKhSFDhrS6PjExEaGhoRcdJwoU5QE4mBgAkqJ0UKsEWGwiKhtbEEjjgZTo2+PVqGw0IS48BJMGJEpdDlHQ83jMzaRJk5CZmYk//OEPKCws7NKL7969Gzk5OcjJyQEAPPTQQ8jJycETTzwBACgrK0NxcXGXXoNISmUBOA0cADRqlbs1ieNupPfBbkeX1PXZqdBquGI0UVd5/K/o7NmzePjhh5Gfn48hQ4Zg+PDh+L//+z+cOXPG4xefOHEiRFG86GPdunUAgHXr1mHr1q3tXv/HP/4RBQUFHr8ukb9U6F3dUoHXMuLaxPNUDRfyk1K1wYRNheUAgJtHd601nIgcPA43CQkJuO+++7B9+3acOHECc+bMwb///W9kZGS0OaCXSMkCteUGAPonO3abPlzGtZ+k9OGeM7DYRGSnx2JQarTU5RDJQpfaPzMzM7Fs2TI8++yzGDp0KPLz871XGZEMuGZLBdqYGwAY4Ao35Y1Sl6JYoiji7V2Orvd5Y7gGF5G3dDrcbN++Hffeey9SUlIwb948DBkyBJ9//rl3qyMKcoHccjMwxdFKcLicLTdS2XGyBqdqmhCp0+Dnw7i2DZG3eDxbavny5XjnnXdw9uxZTJkyBS+99BJ+8YtfIDw83DcVEgWpJrMVDc0WIEBbbvolRUEQgGqDGVWNJnSP4toq/vb2Lseior8YnooIncdvx0TUDo//NW3btg1Lly7F3LlzkZCQ4JuqiGTA1SUVqdMgKjRE6nIuEqZVI7NbBE5WG/FTmR7do7pLXZKiVDWasPGgYyubW8dwIDGRN3kcbrZv3+6bSohkxhVukqIDt0VkQEoUTlYbcbhcj6v7Mdz401v/K4bFJiKnZyyG9IiRuhwiWel0O2hhYSGKi4thNptbHb/++uu9URdR0AvEDTMvNDA5Gl8cKMfhMg4q9iez1Y43/3caALBwfIbU5RDJjsfh5uTJk7jhhhtw4MABCIIAURQBAILgWL7dZrN5v0qiIFSuD9yZUi4DnIOKf+KMKb/aeKgclc5xTtOHpEhdDpHseDxbasmSJcjMzERlZSXCw8Nx6NAhbNu2DaNGjbrkgntESlMewDOlXFzTwY9XNsJis0tdjmL8+/tTAIDbxvbkisREPuDxv6odO3bgT3/6ExISEqBSqaBSqXDllVciNzcX999/v2+qJApCZQG8xo1LWlwYonQaWGwiTlQZpC5HEfafqcee03UIUQuYN5YDiYl8weNwY7PZEBXl+N9eQkICzp49CwDo1asXjhw54v0KiYJUud6xZ1Mgt9wIgoABKa6Vitvumjpa0Yh/f38KD7+3Dyv/ewgmK7ueu+If204CAGYOTUFiVOD+bBAFM4/H3AwZMgT79u1DZmYmxo4di+eeew5arRZr1qxB7969fVMlURCq0JsAIOB/gQ1IjsYPp+rwU7kes9Gj1df+u+8s7n9nL5xD6wAAJbVNeOW2kexO6YSiaiO+POCY/n33xD5Sl0MkWx6/Oz322GOw2x1983/6059QVFSEq666Cl988QX+8pe/+KJGoqAjiiLqjI6ZhPERWqnLuaTBzv2M8o9UuScIAMDpGiOWf3QAogiMyYjHnVdlQqdR4eufKrHknb2wcoyOx9ZsOwm7CEwakIgBydxHishXPG65mTZtmvvzrKwsHD58GLW1tYiLi3PPmCJSOoPJCqvdERTiwgM73Fw3JBkr/1uIw+WN2HmyFuP6dIPZasf9b++FwWTF6Iw4vHXnWGjUKlyRlYC71u/BlwfL8Wr+Cdw3qa/U5QeNSn0LPtxzBgBwL1ttiHzKK+3K8fHxDDZE56lvcmy7oNOoEKZVS13OJcWGa3HjCEd31NrtRQCAZ788jH1nGhATFoLVt+RAo3a8VUzsn4inbhgCAHj9uyIYTFYJKw8u//yuCGabHaMz4jAqI17qcohkjZ3mRD5Q1+Tokgr0VhuX31zhWEgu76cKPP/VEfzLGXKeu2kYesS2XoTwlyPS0DshAvVNFryx47Qk9QabWqMZb+x03Kt72GpD5HMMN0Q+UOdsuYkND7w9pdqSlRiFq/t1hygCf/vmOADgwcn9MG1w8kXnqlUCFl+TBQB47duTaDJ3rfWmzmjG+h2ncOf63dhyuKJLzxWoXv/2JJrMNgztEYNr+idKXQ6R7DHcEPlAfZC13OC81hsAmJWdivuvzWr33F8MT0XP+HDUGs1463/FnX7N93eXYMwzX+OJTw4hr7ACL2w62unnClR1RrN70b77r+3LLnwiP2C4IfIB10ypuIjgaLkBgAl9u2P28FRMH5KM/7tp2CV/CWvUKtznbL15Nf8EjJ0Ye2Ox2bFq42FYbCKyEiMB55o6Zqu8ZmH9a3sRjGYbBqZEY/JAttoQ+QPDDZEPnOuWCp6WG5VKwOpbcvD3X41EaMjlB0HfMKIHMrqFo9pgdg9E9sSWw5WoNpjRPUqHL5dchahQx0rJxyvls1JyrdGMddsdrTZLrs1iqw2RnzDcEPnAuW6p4Gm58VSIWoUHp/QDAPwj/6S7taqj3vuhBABw44geCFGrMMi5ieehsw0+qFYaL319FI0mKwalRGPqoIvHLxGRbzDcEPmAq+UmmMbcdMasYakYmBKNRpMVr+af6PB1FfoWfHOkEgAwZ2Q6AGBwagwA4NBZvY+q9a8TVQa86RyP9NjMgVCp2GpD5C8MN0Q+4JoKHkzdUp2hUgl4dFp/AMC670+hUt/Soes++rEUdhEY1SvOPd7GtVJyoUzCzbNfHobVLmLSgESMz0qQuhwiRWG4IfKBenfLjXy7pVwm9u+O4emxMFnt+ODHM5c9XxRFvL/b0SU1d1S6+/jgHs5wU6aH3S62e30w2HmyBnmFFVCrBCyfPkDqcogUh+GGyAeU0nID587i88b0BAC8v/tMq/2p2vLVoQqcrDYiUqfBzGEp7uN9ukdCq1HBYLKiuLbJ53X7it0u4pkvfgIA3DI6HX2ToqQuiUhxGG6IfEBJLTcAMHNYCsK1ahRVG7H7dF2759ntIlZ/7VjLZuH4DETozm1vF6JWYUCyIwgE87ibT/edxf4zDYjQqvHA5H5Sl0OkSAw3RF5mttrdey7JfUCxS4ROg587W2Hedc6CastXh8pxuLwRkToN7rgq86Kvu8bdBOuMqRaLDf/31RHAuc1C9yid1CURKRLDDZGX1Tc7uqQEAYgOU0bLDQDcPNoxfubz/WVtbqjpaLU5BjhXQ26ry25QkM+Y+tf2IpTWNyMlJhS3X9lb6nKIFIvhhsjLXF1SMWEhUCto+u+InnHo3T0CzRYbNuwtvejrH+0txZGKRkTpNLijnV/851pugi/clDU0429bHPtyPTy1f8DvBk8kZww3RF7m3npBIV1SLucPLH5h0xFUNZrcXztdY8QfPz0EALh7Yh/EtDMWaWByNAQBqDaYUG0wtXlOoHrys0I0mW0Y0TMWN+b0kLocIkVjuCHysmDbEdyb5o/LwMCUaNQ1WfDYhgMQRRFmqx33v70XBpMVozPi8Nur2++uCdOq0T3SMU6lrL5ja+YEgvyjVfjiQDlUAvDU7KFcsI9IYgw3RF4WjDuCe4tWo8Lzc4ZBoxLw1aEKPPrBftz62k7sO9OAmLAQrL4lBxr1pd92kqJDAecqxsGgxWLDik8OAgAWjs/EIGfXGhFJh+GGyMuU3HID5zYKv5vUFwDw/p4z2HO6DioBWPXLYegRG3bZ613hpjxIws3qr4/hVE0TEqN0eHBKX6nLISIAmg6cQ0QeUHLLjcu91/RBXZMZ+hYLctJjMa5PN2Qldmwxu+QYR7dUR7dykNK+knqs2ebYU+vpG4YiKlSZgZYo0DDcEHlZnQJ2BL+cELUKf7x+cKeuTYoKjpYbk9WGRz/YD7sIXJ+diimDkqQuiYic2C1F5GXnuqWU23LTFUkxrnAT2LOl/rL5GI5UNKJbhLbTQY6IfIPhhsjLlDoV3FuSnWNuArlbaldRLV7Z6uiOemr2EMRH8O+aKJAw3BB5GbuluibQBxTrWyx48N0CiCJw08g0TB+a0oGriMifGG6IvKye3VJd4mq5qW+yoMVik7qci6z45BBK65vRMz6c3VFEAYrhhsiLRFFEfbNzR/AIttx0RnSYBjqN462pMsDG3Xy67yw+3lsKlQC8eHM2InWck0EUiBhuiLxI32KFzS4CHHPTaYIgIDkm8LqmSuub8f8+PgAAuG9SX4zsFS91SUTUDoYbIi9yrXETGqJCaAg3TuysQBt3Y7OLePi9AjS2WDE8PRb3T8qSuiQiugSGGyIvck0DZ6tN1yQF2Iypv2w+hp0naxGuVWP1zcMvu4UEEUmL/0KJvEjvHG8TE8bxNl2RHO1Ypbi8Qfpws/VIJf6y5RgA4OkbhiAjIULqkojoMhhuiLzIYLICAKJCOdC0KwKlW6q0vtk97fu2sT1xQ06apPUQUccw3BB5kaHFEW44i6ZrznVLSTdbqtlsw91v7EFdkwVDe8Tg8Z8PkqwWIvIMww2RFzU6W24iuYFil0g9W0oURTz64X4cKG1AfIQWr9w2ggPEiYIIww2RFxld4UbHX4RdkXxet5Qoin5//Ve2nsB/952FRiXg77eNQHp8uN9rIKLOY7gh8iKDid1S3tA9yjGg2Gy1o8E5SNtf8gor8PymIwCAlb8YjLG9u/n19Ymo6xhuiLyo0T3mht1SXREaonbvzeXPrqmjFY144J29EEXg1z/rhdvG9vLbaxOR9zDcEHmRu+WGs6W6zD1jyk/TweuMZty5fjeMZht+1jseT8ziAGKiYMVwQ+RFhhZHF0oUu6W6zJ8zpkxWG+56YzdO1zQhLS4Mr9w2EiFcqI8oaPFfL5EXGU2OXazZctN1Kc4ZU2cbmn36OqIo4tEP9uOHU3WI0mnwr4WjER/BFaaJghnDDZEXuaaCR7DlpsvS4sIAACW1vg03L359DJ8UOGdG/Wok+iVF+fT1iMj3GG6IvMhgcnRLcbZU17mmX5fUNfnsNT7ccwZ/2ezYWuGp2UNwZd8En70WEfmPpOFm27ZtmDVrFlJTUyEIAjZs2HDJ8z/66CNMmTIF3bt3R3R0NMaNG4evvvrKb/USXY5rhWJuv9B17nBT65tws+NEDZZ9tB8AcM/EPrhlTE+fvA4R+Z+k4cZoNCI7Oxsvv/xyh87ftm0bpkyZgi+++AJ79uzBNddcg1mzZmHv3r0+r5WoI7jOjfekxznCTbm+BSarzavPfbzSgN++sRsWm4iZQ1OwdGp/rz4/EUlL0nfg6dOnY/r06R0+f/Xq1a0eP/PMM/jkk0/w3//+Fzk5OT6okKjjTFYbLDbHarocUNx1CZFahIWo0Wyx4Wx9CzK9tBt3jcGERet+gL7FipyesXhhbjZUKsErz01EgSGox9zY7XY0NjYiPj6+3XNMJhP0en2rDyJfcHVJAUCEluGmqwRBOG9QsXe6plosNtz1xh4U1zYhPT4Mr80fxT2jiGQoqMPN888/D4PBgLlz57Z7Tm5uLmJiYtwf6enpfq2RlMPVJRWuVUPNlgCv8OagYteU7z2n6xAdqsHahaOREKnzQpVEFGiCNty89dZbWLlyJd577z0kJia2e97y5cvR0NDg/igpKfFrnaQc57ZeYKuNt/R0Dyru+nTwlzYfw6fOzTBf/dVIZCVyyjeRXAXlu/A777yDO+64A++//z4mT558yXN1Oh10Ov7vjHzPyK0XvM5b3VKfFJRi9dfnpnyPz+KUbyI5C7qWm7fffhu/+c1v8Pbbb2PmzJlSl0Pk5uqW4tYL3uONbqkfi+uw9APHlO+7ru7NKd9ECiDpu7DBYMDx48fdj4uKilBQUID4+Hj07NkTy5cvR2lpKdavXw84u6IWLFiAl156CWPHjkV5eTkAICwsDDExMZJ9H0Tgppk+4ZoO3tmWmzN1Tbhr/W6YrXZMHpiE3183wMsVElEgkrTlZvfu3cjJyXFP437ooYeQk5ODJ554AgBQVlaG4uJi9/lr1qyB1WrF4sWLkZKS4v5YsmSJZN8DkYtrzA1nSnlPeryjW6quyeIOjx3V2GLB7et2o9pgxsCUaLx0y3AO9CZSCEnfhSdOnAhRFNv9+rp161o93rp1qx+qIuocttx4X1RoCGLDQ1DfZEFJbRMGpkR36DqrzY7fvb0XRyoakRilwz8XjOJ+X0QKEnRjbogClZFjbnyiZye2YXj6i5+w9UgVQkNUeH3BKKTGhvmwQiIKNAw3RF7ingrOlhuvco+7qevYdPA3dp7G2u2nAAB/njscw9JifVofEQUehhsiLzm3r1SI1KXISlp8x6eDbztahT9+eggAsHRaf8wYmuLz+ogo8DDcEHmJwb2IH5fz96aOzpg6VtGIxW/+CJtdxI0jeuDeiX38VCERBRqGGyIv4YBi3+jVzRFuDpc3tjsBocZgwqJ//4BGkxWjM+KQe+NQCAJnRhEpFcMNkZewW8o3RvaKg1ajQml9M45WGC76umszzJLaZvSMD8ervxoJnYatZ0RKxnBD5CXnwg1bbrwpXKvBlc7tEvIKy1t9zW4XsfS8zTD/tXA0unEzTCLFY7gh8hLXmJsodkt53ZRBSQCAvJ8qWx1f/fVR/LfVZpiRElVIRIGE4YbIS1wtN1wszvuuHZgIANhXUo8KfQsA4MM9Z/CXLY7tW565cSg3wyQiN4YbIi+w20V2S/lQYlQohqc71qvZ/FMlvjlSiWUfOTbDvGdiH8wdlS5xhUQUSPguTOQFTRab+3N2S/nGlEFJKCipx1+3HEOFvgV2EZgxNBlLp/aXujQiCjBsuSHyAtd4G41KgE7Df1a+4Bp3U9bgCDZzR6XhxZuHQ8XNMInoAvwvJpEXGEwWwLnGDddX8Y2+iZEYlBKNIxWNeGzmQCwcn8F7TURtYrgh8gL3vlIcb+MzgiDg3d/+DI0tVm6ESUSXxHdiIi/gYGL/iAoNQVQoF0kkokvj4AAiLzAy3BARBQyGGyIvcHdLcaYUEZHkGG6IvIDdUkREgYPhhsgLuPUCEVHgYLgh8gL31gtahhsiIqkx3BB5gbtbii03RESSY7gh8oL6ZsciftGcpkxEJDmGGyIvqG8yAwDiI7RSl0JEpHgMN0ReUGd0tNzEhrPlhohIagw3RF7garmJC2fLDRGR1BhuiLygluGGiChgMNwQdVGLxYYWix0AEBvBbikiIqkx3BB1UZ2z1UajEhDFFYqJiCTHcEPURecGE2shCILU5RARKR7DDVEXnRtMzC4pIqJAwHBD1EUcTExEFFgYboi6qK6Ja9wQEQUShhuiLqo3suWGiCiQMNwQdZGr5SaOWy8QEQUEhhuiLuKAYiKiwMJwQ9RFdRxQTEQUUBhuiLqolgOKiYgCCsMNURe5u6U45oaIKCAw3BB1UZ2RY26IiAIJww1RF1htduhbrADH3BARBQyGG6IuaGi2uD+PCWPLDRFRIGC4IeoC10yp6FANNGr+cyIiCgR8NybqAi7gR0QUeBhuiLrANZg4luNtiIgCBsMNURfUu1puOFOKiChgMNwQdYFrzE08W26IiAIGww1RF9Q2sVuKiCjQMNwQdUG9kd1SRESBhuGGqAtc3VKxnC1FRBQwGG6IuoADiomIAg/DDVEXcEAxEVHgYbgh6oI6DigmIgo4DDdEnXS80oAa5yJ+idE6qcshIiInhhuiTnox7yhEEZgyKAkJkQw3RESBguGGqBMOljbg8wNlEATg4an9pC6HiIjOI2m42bZtG2bNmoXU1FQIgoANGzZc9pqtW7dixIgR0Ol0yMrKwrp16/xSK5GLxWbH85uOAACuz07FgORoqUsiIqLzaKR8caPRiOzsbCxatAg33njjZc8vKirCzJkzcffdd+PNN9/E5s2bcccddyAlJQXTpk3zS83tsdlFlDU0QxQdj0URECE6/wREUXT+CaDV8XPn2UWxc9ef9zX3Ne1cj4vOafu57c7zztXiOM9FvPDrznNwwWs5n7FVvRcecz1Jq69f+Lid12/99dbH3M/Txjnnf19mq93xYbPBZLHDbHM8NlntMJqsMDg/jCYrGlusMFntAAC1SsCDk9lqQ0QUaCQNN9OnT8f06dM7fP6rr76KzMxMvPDCCwCAgQMH4rvvvsOLL74oebipMZhw5apvJK2B/EcQgLuu7o2MhAipSyEiogtIGm48tWPHDkyePLnVsWnTpuGBBx5o9xqTyQSTyeR+rNfrfVKbIAjQaVQQBECA4PzTcVwAgPMfX/A1Qbjg8wuuh/ucyz+36rxzcNFrtL4erV7z3Lkq5wHX1xwvIbif4/zH53//gvvzC2p3P49wwdfbf27X61/4PZx7vnPPDef31ZHnxvnfv/M8rVoNrUYFnUbl/tP1eYROgwidBlE6DSJDNYjQahAV6jgWouaQNSKiQBRU4aa8vBxJSUmtjiUlJUGv16O5uRlhYWEXXZObm4uVK1f6vLbuUTocearjrVBERETkG7L/r+fy5cvR0NDg/igpKZG6JCIiIvKhoGq5SU5ORkVFRatjFRUViI6ObrPVBgB0Oh10Oq5BQkREpBRB1XIzbtw4bN68udWxvLw8jBs3TrKaiIiIKLBIGm4MBgMKCgpQUFAAOKd6FxQUoLi4GHB2Kc2fP999/t13342TJ0/i0UcfxeHDh/HKK6/gvffew4MPPijZ90BERESBRdJws3v3buTk5CAnJwcA8NBDDyEnJwdPPPEEAKCsrMwddAAgMzMTn3/+OfLy8pCdnY0XXngBr7/+uuTTwImIiChwCOL5K6MpgF6vR0xMDBoaGhAdzZVliYiIgoEnv7+DaswNERER0eUw3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsBNXGmd7gWrNQr9dLXQoRERF1kOv3dkfWHlZcuGlsbAQApKenS10KEREReaixsRExMTGXPEdx2y/Y7XacPXsWUVFREATBa8+r1+uRnp6OkpISbutwGbxXnuH96jjeK8/wfnUc75VnfHG/RFFEY2MjUlNToVJdelSN4lpuVCoV0tLSfPb80dHR/MHvIN4rz/B+dRzvlWd4vzqO98oz3r5fl2uxceGAYiIiIpIVhhsiIiKSFYYbL9HpdFixYgV0Op3UpQQ83ivP8H51HO+VZ3i/Oo73yjNS3y/FDSgmIiIieWPLDREREckKww0RERHJCsMNERERyQrDDREREckKw40HXn75ZWRkZCA0NBRjx47Frl27Lnn++++/jwEDBiA0NBRDhw7FF1984bdapebJvVq3bh0EQWj1ERoa6td6pbJt2zbMmjULqampEAQBGzZsuOw1W7duxYgRI6DT6ZCVlYV169b5pdZA4On92rp160U/W4IgoLy83G81SyU3NxejR49GVFQUEhMTMXv2bBw5cuSy1yn1fasz90up711///vfMWzYMPcCfePGjcOXX355yWv8/XPFcNNB7777Lh566CGsWLECP/74I7KzszFt2jRUVla2ef7333+PW2+9Fbfffjv27t2L2bNnY/bs2Th48KDfa/c3T+8VnKtYlpWVuT9Onz7t15qlYjQakZ2djZdffrlD5xcVFWHmzJm45pprUFBQgAceeAB33HEHvvrqK5/XGgg8vV8uR44cafXzlZiY6LMaA0V+fj4WL16MnTt3Ii8vDxaLBVOnToXRaGz3GiW/b3XmfkGh711paWl49tlnsWfPHuzevRuTJk3CL37xCxw6dKjN8yX5uRKpQ8aMGSMuXrzY/dhms4mpqalibm5um+fPnTtXnDlzZqtjY8eOFX/729/6vFapeXqv1q5dK8bExPixwsAEQPz4448vec6jjz4qDh48uNWxm2++WZw2bZqPqws8Hblf33zzjQhArKur81tdgaqyslIEIObn57d7jpLfty7UkfvF965z4uLixNdff73Nr0nxc8WWmw4wm83Ys2cPJk+e7D6mUqkwefJk7Nixo81rduzY0ep8AJg2bVq758tFZ+4VABgMBvTq1Qvp6emX/B+A0in156qrhg8fjpSUFEyZMgXbt2+XuhxJNDQ0AADi4+PbPYc/X+d05H6B712w2Wx45513YDQaMW7cuDbPkeLniuGmA6qrq2Gz2ZCUlNTqeFJSUrt99+Xl5R6dLxeduVf9+/fHv/71L3zyySf4z3/+A7vdjvHjx+PMmTN+qjp4tPdzpdfr0dzcLFldgSolJQWvvvoqPvzwQ3z44YdIT0/HxIkT8eOPP0pdml/Z7XY88MADuOKKKzBkyJB2z1Pq+9aFOnq/lPzedeDAAURGRkKn0+Huu+/Gxx9/jEGDBrV5rhQ/V4rbFZwCz7hx41ol/vHjx2PgwIH4xz/+gSeffFLS2ii49e/fH/3793c/Hj9+PE6cOIEXX3wRb7zxhqS1+dPixYtx8OBBfPfdd1KXEhQ6er+U/N7Vv39/FBQUoKGhAR988AEWLFiA/Pz8dgOOv7HlpgMSEhKgVqtRUVHR6nhFRQWSk5PbvCY5Odmj8+WiM/fqQiEhIcjJycHx48d9VGXwau/nKjo6GmFhYZLVFUzGjBmjqJ+t++67D5999hm++eYbpKWlXfJcpb5vnc+T+3UhJb13abVaZGVlYeTIkcjNzUV2djZeeumlNs+V4ueK4aYDtFotRo4cic2bN7uP2e12bN68ud0+xnHjxrU6HwDy8vLaPV8uOnOvLmSz2XDgwAGkpKT4sNLgpNSfK28qKChQxM+WKIq477778PHHH2PLli3IzMy87DVK/vnqzP26kJLfu+x2O0wmU5tfk+TnymdDlWXmnXfeEXU6nbhu3TqxsLBQvOuuu8TY2FixvLxcFEVR/PWvfy0uW7bMff727dtFjUYjPv/88+JPP/0krlixQgwJCREPHDgg4XfhH57eq5UrV4pfffWVeOLECXHPnj3iLbfcIoaGhoqHDh2S8Lvwj8bGRnHv3r3i3r17RQDin//8Z3Hv3r3i6dOnRVEUxWXLlom//vWv3eefPHlSDA8PF5cuXSr+9NNP4ssvvyyq1Wpx48aNEn4X/uPp/XrxxRfFDRs2iMeOHRMPHDggLlmyRFSpVOLXX38t4XfhH/fcc48YExMjbt26VSwrK3N/NDU1uc/h+9Y5nblfSn3vWrZsmZifny8WFRWJ+/fvF5ctWyYKgiBu2rRJFAPk54rhxgN//etfxZ49e4parVYcM2aMuHPnTvfXJkyYIC5YsKDV+e+9957Yr18/UavVioMHDxY///xzCaqWhif36oEHHnCfm5SUJM6YMUP88ccfJarcv1xTlS/8cN2fBQsWiBMmTLjomuHDh4tarVbs3bu3uHbtWomq9z9P79eqVavEPn36iKGhoWJ8fLw4ceJEccuWLRJ+B/7T1n0C0Ornhe9b53Tmfin1vWvRokVir169RK1WK3bv3l289tpr3cFGDJCfK0F0/KUSERERyQLH3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsMNwQERGRrDDcEBERkaww3BAREZGsMNwQERGRrDDcEFHQW79+Pbp16waTydTq+OzZs/HrX/9asrqISBoMN0QU9ObMmQObzYZPP/3UfayyshKff/45Fi1aJGltROR/DDdEFPTCwsIwb948rF271n3sP//5D3r27ImJEydKWhsR+R/DDRHJwp133olNmzahtLQUALBu3TosXLgQgiBIXRoR+ZkgiqIodRFERN4wcuRI3HTTTZg6dSrGjBmDU6dOIT09XeqyiMjPNFIXQETkLXfccQdWr16N0tJSTJ48mcGGSKHYckNEstHQ0IDU1FRYrVasX78eN998s9QlEZEEOOaGiGQjJiYGv/zlLxEZGYnZs2dLXQ4RSYThhohkpbS0FLfddht0Op3UpRCRRNgtRUSyUFdXh61bt+Kmm25CYWEh+vfvL3VJRCQRDigmIlnIyclBXV0dVq1axWBDpHBsuSEiIiJZ4ZgbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpIVhhsiIiKSFYYbIiIikhWGGyIiIpKV/w/5qXuLOBy97QAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "g = dens.g\n", + "ax.plot(g.y[g.jlo:g.jhi+1], dens_avg)\n", + "ax.set_xlabel(\"y\")\n", + "ax.set_ylabel(\"average density\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From f73c3e043a4feb3017ab7181aa3a8e1af2bb144a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 7 Jan 2025 11:23:35 -0500 Subject: [PATCH 3/6] add a multimode RT problem (#307) This seems to work with both compressible and compressible_fv4. --- .../compressible/problems/inputs.rt_multimode | 33 +++++++ pyro/compressible/problems/rt_multimode.py | 85 +++++++++++++++++++ 2 files changed, 118 insertions(+) create mode 100644 pyro/compressible/problems/inputs.rt_multimode create mode 100644 pyro/compressible/problems/rt_multimode.py diff --git a/pyro/compressible/problems/inputs.rt_multimode b/pyro/compressible/problems/inputs.rt_multimode new file mode 100644 index 000000000..5483e39cf --- /dev/null +++ b/pyro/compressible/problems/inputs.rt_multimode @@ -0,0 +1,33 @@ +# simple inputs files for the four-corner problem. + +[driver] +max_steps = 10000 +tmax = 3.0 + + +[io] +basename = rt_ +n_out = 100 + + +[mesh] +nx = 64 +ny = 192 +xmax = 1.0 +ymax = 3.0 + +xlboundary = periodic +xrboundary = periodic + +ylboundary = hse +yrboundary = hse + + +[rt_multimode] +amp = 0.25 +nmodes = 12 + +[compressible] +grav = -1.0 + +limiter = 2 diff --git a/pyro/compressible/problems/rt_multimode.py b/pyro/compressible/problems/rt_multimode.py new file mode 100644 index 000000000..37df4839a --- /dev/null +++ b/pyro/compressible/problems/rt_multimode.py @@ -0,0 +1,85 @@ +"""A multi-mode Rayleigh-Taylor instability.""" + +import numpy as np + +from pyro.util import msg + +DEFAULT_INPUTS = "inputs.rt_multimode" + +PROBLEM_PARAMS = {"rt_multimode.dens1": 1.0, + "rt_multimode.dens2": 2.0, + "rt_multimode.amp": 1.0, + "rt_multimode.sigma": 0.1, + "rt_multimode.nmodes": 10, + "rt_multimode.p0": 10.0} + + +def init_data(my_data, rp): + """ initialize the rt problem """ + + # see the random number generator + rng = np.random.default_rng(12345) + + if rp.get_param("driver.verbose"): + msg.bold("initializing the rt problem...") + + # get the density, momenta, and energy as separate variables + dens = my_data.get_var("density") + xmom = my_data.get_var("x-momentum") + ymom = my_data.get_var("y-momentum") + ener = my_data.get_var("energy") + + gamma = rp.get_param("eos.gamma") + + grav = rp.get_param("compressible.grav") + + dens1 = rp.get_param("rt_multimode.dens1") + dens2 = rp.get_param("rt_multimode.dens2") + p0 = rp.get_param("rt_multimode.p0") + amp = rp.get_param("rt_multimode.amp") + sigma = rp.get_param("rt_multimode.sigma") + nmodes = rp.get_param("rt_multimode.nmodes") + + # initialize the components, remember, that ener here is + # rho*eint + 0.5*rho*v**2, where eint is the specific + # internal energy (erg/g) + xmom[:, :] = 0.0 + ymom[:, :] = 0.0 + dens[:, :] = 0.0 + + # set the density to be stratified in the y-direction + myg = my_data.grid + + ycenter = 0.5*(myg.ymin + myg.ymax) + + p = myg.scratch_array() + + j = myg.jlo + while j <= myg.jhi: + if myg.y[j] < ycenter: + dens[:, j] = dens1 + p[:, j] = p0 + dens1*grav*myg.y[j] + + else: + dens[:, j] = dens2 + p[:, j] = p0 + dens1*grav*ycenter + dens2*grav*(myg.y[j] - ycenter) + + j += 1 + + # add multiple modes to the vertical velocity + L = myg.xmax - myg.xmin + for k in range(1, nmodes+1): + phase = rng.random() * 2 * np.pi + mode_amp = amp * rng.random() + ymom[:, :] += (mode_amp * np.cos(2.0 * np.pi * k*myg.x2d / L + phase) * + np.exp(-(myg.y2d - ycenter)**2 / sigma**2)) + ymom /= nmodes + ymom *= dens + + # set the energy (P = cs2*dens) + ener[:, :] = p[:, :]/(gamma - 1.0) + \ + 0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :] + + +def finalize(): + """ print out any information to the user at the end of the run """ From 128987909d19203c792e27d827d93fd76f7d5eaa Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 8 Jan 2025 15:57:16 -0500 Subject: [PATCH 4/6] use a seeded random number generator for convection (#310) we use a random velocity field at initialization this makes the problem reproducible --- pyro/compressible/problems/convection.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyro/compressible/problems/convection.py b/pyro/compressible/problems/convection.py index 90768cd74..8b8b53872 100644 --- a/pyro/compressible/problems/convection.py +++ b/pyro/compressible/problems/convection.py @@ -42,6 +42,9 @@ def init_data(my_data, rp): ymom[:, :] = 0.0 dens[:, :] = dens_cutoff + # create a seeded random number generator + rng = np.random.default_rng(12345) + # set the density to be stratified in the y-direction myg = my_data.grid @@ -75,7 +78,7 @@ def init_data(my_data, rp): ener[:, :] = p[:, :]/(gamma - 1.0) # pairs of random numbers between [-1, 1] - vel_pert = 2.0 * np.random.random_sample((myg.qx, myg.qy, 2)) - 1 + vel_pert = 2.0 * rng.random(size=(myg.qx, myg.qy, 2)) - 1 cs = np.sqrt(gamma * p / dens) From c1b0584ab959b5a75943b2671b7e161d850b6e9a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 9 Jan 2025 13:39:40 -0500 Subject: [PATCH 5/6] add some asserts to the compressible solver (#312) this catches negative density and internal energy if these become a problem, then the solution is to do a dual energy formulation or to introduce some floors. --- pyro/compressible/simulation.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index af152d319..993f10b8f 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -65,6 +65,9 @@ def cons_to_prim(U, gamma, ivars, myg): out=np.zeros_like(U[:, :, ivars.iener]), where=(U[:, :, ivars.idens] != 0.0)) + assert e.v().min() > 0.0 + assert q.v(n=ivars.irho).min() > 0.0 + q[:, :, ivars.ip] = eos.pres(gamma, q[:, :, ivars.irho], e) if ivars.naux > 0: From 6cc755151f40e53ffa2c8626a795ef517e9dc97c Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 9 Jan 2025 16:04:32 -0500 Subject: [PATCH 6/6] implement a sponge term in compressible (#313) for the CTU solver, the sponge update is implicit in time for the RK/SDC solvers, it is an explicit source this damps the velocity field over a timescale "sponge_timescale" The sponge is essentially the same as in Castro --- pyro/compressible/__init__.py | 2 +- pyro/compressible/_defaults | 8 ++++++++ pyro/compressible/simulation.py | 31 +++++++++++++++++++++++++++++ pyro/compressible_fv4/_defaults | 8 ++++++++ pyro/compressible_fv4/simulation.py | 9 ++++++++- pyro/compressible_rk/_defaults | 6 ++++++ pyro/compressible_rk/simulation.py | 7 +++++++ pyro/compressible_sdc/_defaults | 8 ++++++++ 8 files changed, 77 insertions(+), 2 deletions(-) diff --git a/pyro/compressible/__init__.py b/pyro/compressible/__init__.py index d78b69b9b..9aacc0ca0 100644 --- a/pyro/compressible/__init__.py +++ b/pyro/compressible/__init__.py @@ -7,4 +7,4 @@ __all__ = ["simulation"] from .simulation import (Simulation, Variables, cons_to_prim, - get_external_sources, prim_to_cons) + get_external_sources, get_sponge_factor, prim_to_cons) diff --git a/pyro/compressible/_defaults b/pyro/compressible/_defaults index 27edf68e4..850d3d729 100644 --- a/pyro/compressible/_defaults +++ b/pyro/compressible/_defaults @@ -21,6 +21,14 @@ grav = 0.0 ; gravitational acceleration (in y-direction) riemann = HLLC ; HLLC or CGF + +[sponge] +do_sponge = 0 ; do we include a sponge source term + +sponge_rho_begin = 1.e-2 ; density below which to begin the sponge +sponge_rho_full = 1.e-3 ; density below which the sponge is fully enabled +sponge_timescale = 1.e-2 ; the timescale over which the sponge should act + [particles] do_particles = 0 particle_generator = grid diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 993f10b8f..dd9e65e71 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -159,6 +159,29 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None, problem_source return S +def get_sponge_factor(U, ivars, rp, myg): + """compute the sponge factor, f / tau, that goes into a + sponge damping term of the form S = - (f / tau) (rho U)""" + + rho = U[:, :, ivars.idens] + rho_begin = rp.get_param("sponge.sponge_rho_begin") + rho_full = rp.get_param("sponge.sponge_rho_full") + + assert rho_begin > rho_full + + f = myg.scratch_array() + + f[:, :] = np.where(rho > rho_begin, + 0.0, + np.where(rho < rho_full, + 1.0, + 0.5 * (1.0 - np.cos(np.pi * (rho - rho_begin) / + (rho_full - rho_begin))))) + + tau = rp.get_param("sponge.sponge_timescale") + return f / tau + + class Simulation(NullSimulation): """The main simulation class for the corner transport upwind compressible hydrodynamics solver @@ -395,6 +418,14 @@ def evolve(self): var = self.cc_data.get_var_by_index(n) var.v()[:, :] += 0.5 * self.dt * (S_new.v(n=n) - S_old.v(n=n)) + # finally, do the sponge, if desired -- this is formulated as an + # implicit update to the velocity + if self.rp.get_param("sponge.do_sponge"): + kappa_f = get_sponge_factor(self.cc_data.data, self.ivars, self.rp, myg) + + self.cc_data.data[:, :, self.ivars.ixmom] /= (1.0 + self.dt * kappa_f) + self.cc_data.data[:, :, self.ivars.iymom] /= (1.0 + self.dt * kappa_f) + if self.particles is not None: self.particles.update_particles(self.dt) diff --git a/pyro/compressible_fv4/_defaults b/pyro/compressible_fv4/_defaults index c9f936e0d..18d1c3707 100644 --- a/pyro/compressible_fv4/_defaults +++ b/pyro/compressible_fv4/_defaults @@ -24,4 +24,12 @@ grav = 0.0 ; gravitational acceleration (in y-direction) riemann = CGF +[sponge] +do_sponge = 0 ; do we include a sponge source term + +sponge_rho_begin = 1.e-2 ; density below which to begin the sponge +sponge_rho_full = 1.e-3 ; density below which the sponge is fully enabled +sponge_timescale = 1.e-2 ; the timescale over which the sponge should act + + diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index df662ae4a..0097832b6 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -1,6 +1,6 @@ import pyro.compressible_fv4.fluxes as flx from pyro import compressible_rk -from pyro.compressible import get_external_sources +from pyro.compressible import get_external_sources, get_sponge_factor from pyro.mesh import fv, integration @@ -48,6 +48,13 @@ def substep(self, myd): (flux_x.v(n=n) - flux_x.ip(1, n=n))/myg.dx + \ (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + S.v(n=n) + # finally, add the sponge source, if desired + if self.rp.get_param("sponge.do_sponge"): + kappa_f = get_sponge_factor(myd.data, self.ivars, self.rp, myg) + + k.v(n=self.ivars.ixmom)[:, :] -= kappa_f.v() * myd.data.v(n=self.ivars.ixmom) + k.v(n=self.ivars.iymom)[:, :] -= kappa_f.v() * myd.data.v(n=self.ivars.iymom) + return k def preevolve(self): diff --git a/pyro/compressible_rk/_defaults b/pyro/compressible_rk/_defaults index 293a4ec6b..8f1aa0467 100644 --- a/pyro/compressible_rk/_defaults +++ b/pyro/compressible_rk/_defaults @@ -26,3 +26,9 @@ riemann = HLLC ; HLLC or CGF well_balanced = 0 ; use a well-balanced scheme to keep the model in hydrostatic equilibrium +[sponge] +do_sponge = 0 ; do we include a sponge source term + +sponge_rho_begin = 1.e-2 ; density below which to begin the sponge +sponge_rho_full = 1.e-3 ; density below which the sponge is fully enabled +sponge_timescale = 1.e-2 ; the timescale over which the sponge should act diff --git a/pyro/compressible_rk/simulation.py b/pyro/compressible_rk/simulation.py index f0306403a..7378fd478 100644 --- a/pyro/compressible_rk/simulation.py +++ b/pyro/compressible_rk/simulation.py @@ -33,6 +33,13 @@ def substep(self, myd): (flux_x.v(n=n) - flux_x.ip(1, n=n))/myg.dx + \ (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + S.v(n=n) + # finally, add the sponge source, if desired + if self.rp.get_param("sponge.do_sponge"): + kappa_f = compressible.get_sponge_factor(myd.data, self.ivars, self.rp, myg) + + k.v(n=self.ivars.ixmom)[:, :] -= kappa_f.v() * myd.data.v(n=self.ivars.ixmom) + k.v(n=self.ivars.iymom)[:, :] -= kappa_f.v() * myd.data.v(n=self.ivars.iymom) + return k def method_compute_timestep(self): diff --git a/pyro/compressible_sdc/_defaults b/pyro/compressible_sdc/_defaults index 438576e9d..1fb4e9e13 100644 --- a/pyro/compressible_sdc/_defaults +++ b/pyro/compressible_sdc/_defaults @@ -22,3 +22,11 @@ temporal_method = RK4 ; integration method (see mesh/integration.py) grav = 0.0 ; gravitational acceleration (in y-direction) riemann = CGF + + +[sponge] +do_sponge = 0 ; do we include a sponge source term + +sponge_rho_begin = 1.e-2 ; density below which to begin the sponge +sponge_rho_full = 1.e-3 ; density below which the sponge is fully enabled +sponge_timescale = 1.e-2 ; the timescale over which the sponge should act