From d8771b8a4095ba1ce521be51415b4c149d3f3d7c Mon Sep 17 00:00:00 2001 From: Areef Waeming Date: Fri, 11 Oct 2024 12:29:09 +0100 Subject: [PATCH] Add CosmoHamTagging --- Examples/Cosmo/CosmoLevel.cpp | 42 ++++++++++++-- Examples/Cosmo/SimulationParameters.hpp | 8 ++- Examples/Cosmo/params_cheap.txt | 10 +++- Examples/Cosmo/plot.py | 18 ++++-- Examples/Cosmo/plot_efold.pdf | Bin 103746 -> 103746 bytes Examples/Cosmo/plot_lineout.pdf | Bin 101788 -> 134978 bytes Examples/Cosmo/plot_rho_mean.pdf | Bin 144305 -> 144305 bytes .../CosmoHamTaggingCriterion.hpp | 52 ++++++++++++++++++ 8 files changed, 116 insertions(+), 14 deletions(-) create mode 100644 Source/TaggingCriteria/CosmoHamTaggingCriterion.hpp diff --git a/Examples/Cosmo/CosmoLevel.cpp b/Examples/Cosmo/CosmoLevel.cpp index 0b501e38b..9ec48ba82 100644 --- a/Examples/Cosmo/CosmoLevel.cpp +++ b/Examples/Cosmo/CosmoLevel.cpp @@ -18,7 +18,7 @@ #include "NewMatterConstraints.hpp" // For tag cells -#include "FixedGridsTaggingCriterion.hpp" +#include "CosmoHamTaggingCriterion.hpp" // Problem specific includes #include "AMRReductions.hpp" @@ -78,6 +78,25 @@ void CosmoLevel::initialData() InitialK my_initial_K(scalar_field, m_dx, m_p.G_Newton); BoxLoops::loop(my_initial_K, m_state_new, m_state_new, EXCLUDE_GHOST_CELLS); + // Calculate constraints and some diagnostics as we need it in tagging + // criterion + BoxLoops::loop(MatterConstraints( + scalar_field, m_dx, m_p.G_Newton, c_Ham, + Interval(c_Mom, c_Mom), c_Ham_abs_sum, + Interval(c_Mom_abs_sum, c_Mom_abs_sum)), + m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS); + CosmoDiagnostics cosmo_diagnostics( + scalar_field, m_dx, m_p.G_Newton); + BoxLoops::loop(cosmo_diagnostics, m_state_new, m_state_diagnostics, + EXCLUDE_GHOST_CELLS); + // Assign initial rho_mean here + // from rho = 1/2 m^2 phi^2; + // phi0 = A sin(2 pi n x /L); + // mean of phi0^2 = A^2 sin^2 = 0.5*A^2; + // m = m_mode + m_cosmo_amr.set_rho_mean( + 0.5 * m_p.scalar_field_mode * m_p.scalar_field_mode * + (0.5 * m_p.initial_params.amplitude * m_p.initial_params.amplitude)); } #ifdef CH_USE_HDF5 @@ -141,17 +160,28 @@ void CosmoLevel::specificUpdateODE(GRLevelData &a_soln, void CosmoLevel::preTagCells() { - // we don't need any ghosts filled for the fixed grids tagging criterion - // used here so don't fill any + // Pre tagging - fill ghost cells and calculate Ham terms + fillAllGhosts(); + Potential potential(m_p.potential_params, m_p.L, m_p.scalar_field_mode); + ScalarFieldWithPotential scalar_field(potential); + BoxLoops::loop(MatterConstraints( + scalar_field, m_dx, m_p.G_Newton, c_Ham, + Interval(c_Mom, c_Mom), c_Ham_abs_sum, + Interval(c_Mom_abs_sum, c_Mom_abs_sum)), + m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS); + CosmoDiagnostics cosmo_diagnostics( + scalar_field, m_dx, m_p.G_Newton); + BoxLoops::loop(cosmo_diagnostics, m_state_new, m_state_diagnostics, + EXCLUDE_GHOST_CELLS); } void CosmoLevel::computeTaggingCriterion( FArrayBox &tagging_criterion, const FArrayBox ¤t_state, const FArrayBox ¤t_state_diagnostics) { - BoxLoops::loop( - FixedGridsTaggingCriterion(m_dx, m_level, 2.0 * m_p.L, m_p.center), - current_state, tagging_criterion); + BoxLoops::loop(CosmoHamTaggingCriterion(m_dx, m_p.center_tag, m_p.rad, + m_cosmo_amr.get_rho_mean()), + current_state_diagnostics, tagging_criterion); } void CosmoLevel::specificPostTimeStep() { diff --git a/Examples/Cosmo/SimulationParameters.hpp b/Examples/Cosmo/SimulationParameters.hpp index 36e9d4867..ad3c169d6 100644 --- a/Examples/Cosmo/SimulationParameters.hpp +++ b/Examples/Cosmo/SimulationParameters.hpp @@ -36,8 +36,13 @@ class SimulationParameters : public SimulationParametersBase pp.load("scalar_mass", potential_params.scalar_mass, 0.1); pp.load("scalar_field_mode", scalar_field_mode, 1.0); + // Lineout params pp.load("lineout_num_points", lineout_num_points, 10); + // Tagging params + pp.load("center_tag", center_tag, center); + pp.load("rad", rad, L); + #ifdef USE_AHFINDER double AH_guess = 8. * initial_params.amplitude * initial_params.amplitude; @@ -59,8 +64,9 @@ class SimulationParameters : public SimulationParametersBase // Initial data for matter and potential and BH double G_Newton; - double scalar_field_mode; + double scalar_field_mode, rad; int lineout_num_points; + std::array center_tag; InitialScalarData::params_t initial_params; Potential::params_t potential_params; KerrBH::params_t kerr_params; diff --git a/Examples/Cosmo/params_cheap.txt b/Examples/Cosmo/params_cheap.txt index e289e8dba..a94b5eaad 100644 --- a/Examples/Cosmo/params_cheap.txt +++ b/Examples/Cosmo/params_cheap.txt @@ -52,6 +52,10 @@ scalar_field_mode = 1.0 # lineout params lineout_num_points = 10 +# Tagging +# center_tag = 0 0 0 +# rad = 7 + ################################################# # Grid parameters @@ -66,15 +70,15 @@ N_full = 16 L_full = 1 # Maximum number of times you can regrid above coarsest level -max_level = 0 # There are (max_level+1) grids, so min is zero +max_level = 1 # There are (max_level+1) grids, so min is zero # Frequency of regridding at each level and thresholds on the tagging # Need one for each level except the top one, ie max_level items # Generally you do not need to regrid frequently on every level # in this example turn off regridding on all levels # Level Regridding: 0 1 2 3 4 5 -regrid_interval = 0 0 0 0 0 0 -# regrid_threshold = 0.5 +regrid_interval = 10 0 0 0 0 0 +regrid_threshold = 500 regrid_interval = 0 0 0 0 0 0 diff --git a/Examples/Cosmo/plot.py b/Examples/Cosmo/plot.py index 2c201e6ef..b918c485c 100644 --- a/Examples/Cosmo/plot.py +++ b/Examples/Cosmo/plot.py @@ -12,8 +12,15 @@ plt.rcParams.update({'figure.figsize' : '6, 4.2'}) plt.rcParams.update({'figure.autolayout': True}) -data = np.loadtxt('cheap_ex_data/data_out.dat') -rho_dat = np.loadtxt('cheap_ex_data/lineout.dat') +# Read header +with open('cheap_data/lineout.dat', 'r') as file: + header = file.readline().strip() +coord_interp = header.split() +coord_interp = coord_interp[2:] + +# Read data +data = np.loadtxt('cheap_data/data_out.dat') +rho_dat = np.loadtxt('cheap_data/lineout.dat') t_data = data[:,0] chi_mean = data[:,3] rho_mean = data[:,4] @@ -47,11 +54,13 @@ dt_multiplier = 0.25 dt = dx*dt_multiplier # dt = dx * dt_multiplier = 0.5 * 0.25 t = 1/dt -num_t_step = 25 # lineout every time = num_t_step +num_t_step = 50 # lineout every time = num_t_step t_plot = np.arange(0,len(rho_interp),t*num_t_step) +x_tick = np.arange(len(coord_interp)) + for it_plot, t_plot in enumerate(t_plot): - plt.plot(rho_interp[int(t_plot)],color=cm.Reds(8./(1+it_plot*10.),1.), + plt.plot(x_tick, rho_interp[int(t_plot)],color=cm.Reds(8./(1+it_plot*10.),1.), label='t = '+ str(t_plot*dt), marker = "o") plt.title('lineout of ' + r'$\rho$' + ' along x-axis') @@ -62,5 +71,6 @@ plt.ylabel(r'$\rho$', fontsize=14) #plt.yscale('log') #plt.ylim(0.,1.3e-3) +plt.xticks(ticks=x_tick ,labels=coord_interp, rotation='vertical', fontsize=10) plt.savefig('plot_lineout.pdf',dpi=256, bbox_inches='tight',pad_inches = 0.1) plt.close() \ No newline at end of file diff --git a/Examples/Cosmo/plot_efold.pdf b/Examples/Cosmo/plot_efold.pdf index c98cbe567b0d55afda3279a691366728e2d68672..e3b829054f1dc35843d78e339df2f80de23d79bd 100644 GIT binary patch delta 28 kcmX@KitW%UwuUW?9Fw^W4Gawpjf@P8jkXI-W;9>|0FEXII{*Lx delta 28 kcmX@KitW%UwuUW?9Fw^WEDbFU&5SKf4YmtSW;9>|0FQ+TSO5S3 diff --git a/Examples/Cosmo/plot_lineout.pdf b/Examples/Cosmo/plot_lineout.pdf index 97318f5061735d0abfe0f3ebcfd2ed902bd3d8b1..c4e9fa84c04ae41f22af75c1efc5418c582a4d26 100644 GIT binary patch delta 36857 zcmZ7dV{~3o*F6r$wr$&J)Yxts+l}p1qyOE6r6ZSsU)h#27s148bvF~*RyQ8v8{z4sg z{o@m_>LBtZHh3gKrRl`;DX)m*^UyEaIKz>?~dnfO@YM@T8Gln6Qw5)r1};xMh?5_Po$j?Io{cvnta$V!^|*s!vn2tBoZ5+ z74%Si;KM;nZHKO`mzP4<(JaPH$`&lS1E7FBY?;gw!AM zpQUVASWJo_IoVs8!W8U7NMbqZ-C^WgiguxnE%sx4huev3G=Z|)5wwmk2(>}i8tN3@ z(4M$^WC78ICMkywg_5$jOKBLg-^tseH~a}Cjln?`Pve;eVuF{xDztZ(-eH-8zw6Xe z+mAr{?1KZ&ySojMA;w$Qy1N?Eq_8=l65%&I5>)5}Ntpx(EI)tFjtnCA+)Kh+Pbb11 zyJVhtgLKr%KPt`kgfFKsDnQsjFRG$xE;s@cPT!O!p!}6UAbKo|}M`JiHug58>%Kc&woGHD;nLIs3 z!#i5dmrhBJVlTh+yDO20#2X$oiwrx|#_?88+o#e~(Yj*Ihl)y~3 zxF#GdCnXdw(l5AB6~1fwb=t8S&YlMlBqV=1~w_mU47XP{?CxI9OjmW z*gW-T_z~qkyS(JWE%H(_{w@7v;oI=D#1khHjNt^Ro3?}xI|@{`K7*WzQwSZo?=PE@ zSjmw&eWUu+D?$Pv^ieGa;s#B5hIK7Q%?na2Ks(<1Pko)n@tEl<-$LpD-{yoXAn=_O zT+mykZN*n5a9*-=KzelFplOiyX2H%uRqtrjqh357s~2?ncW$ndq;~UQCXwt_^V~0| z=y_B#yPp0$q`I!7!1+jZ?8T2)pOTmW7aTuX$BorIK|>>G5!_1lQA@(Oi*CAuf^sF55&A`P_v0b#NC6 zrrsvmez~cknknWp%Gm*>$T4*rJ?z4fZSl6#%IH}%G)xdD;7inZoQPhAL_Ynx;2+P8 zskHB?;i~xc~MUsC_ z<~zcrSZ5yzztyg>=?8__dJ-a#(KeX%=;)9eCsA)zDoIg8)CRK31RLrj^s-`f0y*j{ zN0#ZD#fS8cjuEBd? zv}_nR#Dt6B;DmjwiHsqnCJ@u(wCq0S!hg%{r!%t_{OYQGN8~{UKUEp@)I1)CWBJ9 z(FU_Tk{m_s#UniSdH7r=KruF_*KeyFeND$w5q4n~70Z{v2rDh#!~~VJbxHeNCXA{N6Gg8JBEvmlZ1zz zorIg4n}my(gM^EVi-en#g@l`p?H~64KXUT0lW=l!{)@kR{s!Bhee5ipaGWeGB>!*x z<>BDq0RAofkLYg^FAEz9D+}v?PW{X9H~udCVf@GRKSlr8d3Z?J|HR+*zbgM>|2_Ht zJ@NmU|3?+@Uupk4=z_6MI{j_ zG53#Q|5*Nxe~kTC^Vjshb^qpU|FPh&MSt-BwD|W&`Gd3bvi^xb`~KSXUyPfLiS6Iw z--3VSe-343`L8GcGwlCv8ynmI?iL;m3kkD?wV}Ba2?rdrgq^Lkh=`p#i8cu%2MY-! z+rQ>mc>ng65{CqepAeaZoYIB_N(aUM@3ti{A`+&+AcJCovL_fKV$iekaFa00ncA8= zTd@6Ug!hkt|4#qb|2HLJXJ=vKNl^|5B}?%`2894Hj_(B#;=$r)IHkBn-Q3)`5m=mM zBAo+Ca-6wgXr!dvq?%@QXBe-DCj=h+4&E1@JNvc9_+Ki`yDl~+H$9)AREQ6QhHk-7RZcsZkKX9=c5n^ouYP$-yFtQ+MDXh&Vdm6rf`) zk_S1xf;?U4uKc0a@EKzCo~y;FpP{ln@6^!y9n=!LIo(gpLCO z=p%s)LA;u`b*T$kLo|L^{R~e_3;dBE-a);zfo1^O2nm`t*=a$jgzUzZ=mW5^_irbz?7G(-*m z>)XvZ2%5n@(tTrA%I{&5Eeq%IGfB32kAeX#k@-v$HL zz*q41kLdTVuI9V%A??a40l&tS{rnDH2JNU?J@ zAK5KhIUfbJU_%6_W)DKWY*heN>@CXAruw%FQ?S%uOz{TDu`O&%^JB+CthWz-1ETlK z0cXM^NK9MM9iS{CMCY4k5|kM5Hl(-X!Rj7kUm*`zu-XjnckgRkkmrC(^Ei>r~<(C12J=m ztvt#1$Ly}Kh4eHs0B0G^SV~J5(pmx?TiF6op#3zTP#wTh=ev^?k+uTd4`HO8!e&1> z9Ynus2r&Kj|C}NT6Ygq8$rrna9*q)@g%yJ$GQBP{DmnDd@-)mc2=2>)Thv$mpt&btrSu;LFq+s>AcBs;VUX;HF_Bw-N)6CLls#`-H?s_TTOLaTNUz07 z342|J21vAZF=jA!u?ZQ&##s)vg8bM+d&1}MnIrc+mw-jkqrtUQ4ElD!5P859wvQOf zu~7_6kM=97Fpp3O7j{tjE_r}bqNf;&=gc2mdLLDIXh@I@2D?QP;y2kx7>N!B5o%E6 zHhF-(vgaH>4EEk1Jn#As{)PAtzVa6b#88;<{%jV%K^1Pc5csqC>~GNxO*0C>)EaZh z{ELmX4z%SxPK3WAX8G*^WK{_MeT*q9v>lW$qh(ohDS&7vE_8?4A^0|0F(EVb=Hl#Q z4S&1Oc;)JICd(hm)MHtRuwmHyid&)(g8s;sHxME`+JgQlMPL8ihu^bc5V`jS2sNL1 zmihfrR1`;G3wjNG6UDYf?`i=6jI@Fv>oW}fAoV3F=HKMtPZE>;Ct0LM=@~hL4XW{c zgAmbZ5h_hZ`AYPMFzEmt#&GWrLXf7T<4KheTyTy$cD>6QIq%TS&T{+I%-m}GK4?`m z+!1m$UC1LQl+6?gN>~Q&Fe~cs+dI-VNb>!D~W?<;qNomKl5GI#z z=Q*~UXRmcZQcL1O^MIPwK=N<_V%_<>VgU#IJBiyCY?`)}fcrj53{d6-2TM_$zY*2|E|%7mDrE<`y2=ST+A`l?41vM@@9$+ulPG$fnpEK|F4W62&PJWj9O7B-28yfkHA(3Z}eCCfBV>6`9 zgpq-1@~3n0x;?KGrk?{k@}A=2=NQ zj6{QSohA4^rVGTLZ(W`y!<{o0YD8PZh>}ghKMQ&XD`;wyd2tNcmX#)yHV&39F=EcJ z!cK5*5)+2MDejxuBxXn6Wi8D7AS2*6HA17L!NO?OYRmGGwZtCYqemvjQeO2RmR|SNr$bEva2^SJFGD zG)VrQ{>*h21s_i1aoa^Z1HMZCs-42_otGiOhf``L3qB z`Nyiea)kFYPPy5J-W zJH@874_jEx8fbg-?M)MoFG?NPqq8q4;jK0 ze*$(Z#?x&jv0ornd3IZ5mK;qs-dx*Fcjc(hj^%&yB2dxsqL9yeL3Wrw>U@uT?0a=r zhi8O?BP73guA9eOQO6Q9Ye0M5R6aeFWs^*S{}tve@8ZSD{gj!4LD`ZfJx%b-Y75lH zTv1uUt;c%-n|*WAXzossF7$AJn`G{efCxYsCYw#6l_r*(qu4BQUC`H#S|H{gW5RYI4HXs@F>muyIJ5&=boxDW>`upQ^g(`39n@541$k*zyM7Y zrq8(K_;h=F?uNEfM1*OvlII=gl{K}CFFbS?%NO#hON#l|x3}sn1%Rp@X?4qo&<^Ob z@>5^Tz~vXJcsAMWjVZF%5*G0thUHH0i^g;OHM!dCtN)Ds94wVris%qj`aIIS9uYcxMhRY-_cUR)5FWhuY zKZq#ii-)hGV*>^1dxu5$L zgmrt@jVi1!n&|zKI2U9_EhE}|>SXw}3!D|kR27NvIsxHL@jsgvRX3OIXfuEw0}^n= zAY>Dl8h0D>?^6n2h8zfwbKohG6qwq*z{vz~RLC#ItgS;A8+?hRt<3her!;xtOZus} zl65$}oCZ`e*zrRt-L^M+QILjSj@!SC3M|fGUVf>^ggTAJJ!o&ID%O zaDR3hW#EVZwt8}C?q1lNq-$g=2mQ%RVn$sK`z@3}Ak#a;2q_&dNQMyTRN9w{y zgNHmHX_ga(R2@x~-S>5#gt4ePdti`Ln=f8{yIek=ay#T+MAh^AN<^2!&b`|8sFs&W zM|+vqqBOtpLO_7A!5cUGi3mH4bD@1Bt;nDdFGXrKFEuT2Y3}WcBb}#qJz8#p(5Gs3f zF=_{m{5&GX{QP(o6{9jps=nFiE10;a$Dp7*J4rG&F4jVYfLSVF+3-DJpGxI?t6lk7 z$I+I{Ej=?jZ?+))y%5EwH%dRqk1%(J^-|16YEUrZ^KmB^#!oPy>hvf%WpmQZ1N7UM zvP+UL9=zxkMKdhOa&4%Q#H^EF;OJ7}D8Tf5lRYw6M=La9O{n>Ejj`EWN<z&08}n4X4kkm6kg)-|<&{eZLmq+iYhL zEnUM;S5zUIm!qa2vQ*2qh(Ic5QBNrLXvq)$6wb5?$XD9}rwI+8v46(Vva%^l4tzo@ zYOfnfx;s>&gWYhal9~AEpyxO3Jxp8Y@<+=;D=VecJRrGPUpzZKRzWoZtkPcRrBJKX zGWwFd=9pXLjdV1krL8V+XY*g_bRMs@TbX~pFlnGY;Ja#BSY{Z5R)q5wCST>Vh(*ds zKTSh#yUgeVq7Mh+9&o-pk;wifYQc`<4jbwX(-_UV;+l<mPh84a7;t*TA z{#Mf2S(P}`Qj7Y;YYzhFn_FFBIZ=hMq2=mAY?`GXO$kaucdDaUSDhebiS6K(zn{wR zGzP#0v42?Hd~{YiKCHZeWuvg;5GFD4oGT42n9)>-77dP95}7#p<{Bdi@hN`!R0mfO zMk7=WScR%5l*jR?9cx;R?XoQ+$!Ld*jI&6YajrPGg3X8la_GGg?&X$W+es<<3A4+f z1BcD)!h#B8mL^w!XpCa$8~;A=ThI=R{mFO4KbNDp|GL%(@xywylHb$JxUNneT27Rj zoHkIcnAwF@Z?}jjBQojhjnk*byf!Rh< z9W2YL&bWxrmi1g|Y>+$;6Mk`1dX1nG7ZbN(#|D@%<7N#vvucCD!l4gi@C0??Dvk|-ukDlW7ZJ@W?VJ$F3N2 zB{|~H=<{c!O4l2PMUzJ*5%{QY(sjno0oWa%;4x&!=_E?xP$in54z)x%_OHH2`6~w5 z4{7O&+5Euy>a~8%c#tv6cam3uT!Xean(AS3J#%LiudYMg+NJJIBL>=oWQUW{rm)1?oA7 z3BzDPK$|R*f~M`uFJ?cTTaSHMs&j-&Z+i(2*2HsA&&rMEmEAED<*dUlH+~fx2c5o~ zTn3T1x7;TEO_-NZbBw^crR7$QJjGGiZ)`q}n=krcNquU*d;53H`~F+&Ji(5aKOW z;C8`HAYNrM&^_IUy)cq_X2^gRbn=Aj6XqJq6>bBz7<~ecq!qLrr$Se&60mRC7y$j{ zR|4*H7*^x1RQ35o6$K*z^mG9BX~;_a&|lm$ng#uw|PajYQv z99B0Rt=lxHS=`t6;Nyqv_fo8EmqKl@lVtC^hqkU4P^%#hWUapJFo;av5<41+Mg}MD z-VJxl3FC5b=<791*;DIC@&eOgo%%Z*3kw7nr?3k5IX^ulBA{GRMmJ*aB>T@{jh`Kc zs*5R8u=R_E7uvJEs)}VIJR+y7c_zZF4a*jojnzd!#pLXh4reyi)lAvzG22+5`-F23 z6_1gVO39d>15cqljL?jn>Z2iY;If#E?y!4c@DP4BcgxCM4n7j!CIZ#*59Q`=MvBM9 zf#r)E?#=PlFJizT#0fVz7UW$}&vM+euox8Xs14DD?77X#c8T?ecDh zAG%F$eL3b0GHXMYRt_#B`^uZJk+ShE@mZr&Hz_gtYGx6n^``I;+sH`=H_`v}yES1& zyQhKdk8H-EUk7VF?M65{xk*n}>V~NS>uuG>rn*tSbf$0i1bl$x=A#&gSftVI>Q_qx zT};i2-VK^zVL4-`;|9yS*Juv96;K$+L>CKxYXZX=`OX|CLylBKiuo@$SV;;JYoJe* zNqid8ENn3|{X1NBCES}WZ59W%&jNeekZBLxQgpBVS1IV?49~l;JTF8I$S(!}{s-j4)!iF7kP`g;m zS$LyUFVR~kVsBZtp{9P!VDc{H_C~=j>m46M)3*U?^3Lx-;(f zZ7qJDn{x~&Sy}>Zivw7E=f|~Gj1)(xESos03>*^=Y4qX4-m>#ch{YGpFpVfRa%Vv}#okrJ0ak<;0 z1k}Ntk!*LOXy@vU<-Y6q2?u{n<5joJQVe4Y4aue>x@IX~;mRV`dCNuBXM)~1q&WF* zzh0(NCwxg@V+#3ZBEoh+&!bg;s5=h4%YL>$>5#kK;d4vy>R5ogvdjwDs2o3R)yboo z@H9vIcMSe}QR=md@%&jn`dK9lB_jPdyiw+$ll6Vd;qO^vat-yEfjbn0b*yDD1;R7Z zEpJgO-SG21+3n#N>fBf|Y&WIgF3Am)a79mdVm|r6`26Q=w%_{BNaE^t(g~skQr|XW z9g#1gcA=Wl=%NqWT&Q@1t&7<&jnPI1?-a;X5xp%GO5g8gY^$(j#mN+qI_-o#6)H(8 z(8h>XFqZ6=Rlrs%y~(CDcO*kRo&uPd)^+G~4l0Ph`(7RRHHmMI5A z8W@tfWG}PN%6uxpbdN?7+Lzj6CNpItBk!TNMqMy%Z=WP5;VKwDvv_8vZeluAj8_U& zbj0D$8)lW(_7C}4cI}TUJYwWK4>P4G|!tDEInR24(3o^DFsRr+9@2QoI(5xu`Pwv-@W`6xtv-+9xTeD#j}7ZRH7dk zF~~)UwfDy)9ISXFD+Hr@KuHZJE$ics7lHpFjxpqJjPEbf;)Zy_X-CAf#7^dlzAlt< zR_)AoOnCI{sdccD}Cx=2{gXbqj>M0Y2$X2OwX!Q zOv;ux(e%;sB;uAU)=KjlGp`Tgyro#gfcJXEHQ=(^v6Q$q6fn~Ea!>%Or^+pXlU0-C zq)GWX#VXDA`StZZo8VxEXbOwwK?J?}jnpl@CV}R1SW9^4Tl=p14XVEfiTZk?ZEtAYokIL zHtO(&o zo5kgHxsIe7+@fANq3K1%cCr8qoB``;Pqd~eCfJGZ8GeetWZ6oo(K-jn6y_5(nC0Ft zm(}p6ykpp#mW(Ri#61%yqq;%$``pWNwr5Z|KOMV_1aVa=+)4*X)l4^)AA*#?*xkSt zG`rDp7At+<2v8yFk~=Q1>kK!kHn64di!>gUEigZSRLRryDxe-B4YMa$?tw?L=@J>tNVYeYH8tBo&o1X!)suI#e9mmBj<}w(Pw$Qld zPb}7gKT)HAZ&eG7IXi87-u49dH8T?aWI#NisAD>+NOw_D8Bzr1&b;mJS|012Q{hm( zL?ty#J&;50FrZYs*n^8A8OZA!P&$S%;!mMJXfoAnqrLK9RYhAfk+i!b?HAh735rgZx2y6(gq|X`&?gDR#Xv!>b}cP_GZA+%R3Z z#!mjlmu=#ma7VhF&lnw^W1g<~tORyVBE3+T)l@-1(|d}+!8Sqecjpe*cT_6H5w6sx z5cB2dWWMRd4#wHwgzEuAdf!Uj5wxbZsl8siT5M7uJG>R(b4n0}1A|9+JvK+Vn6&30 z72U49s*<*VZ6*@Lz!`s2CL20|EWsK#4SK?_5D|F)D5EbGg|DK0%Z^mmOW*vV_mu+i zbF#!IH#a-_Xo0t!kV++d+X?>$u8=d;B?fcrk_xKTtpv{N#bd=%^(a`j&63fU=(!B& zHRL>(13qcs*B~XsuXfj9I3BQ0)-$c&YBg*EYm?y&1WAr83i(XFN9uE=xTNIsFUe5V zxjmo5EA8;BY3k7V%72M{0}py$BS`A@E1g8F4mPJ3(}TLDw-qm}CR*7JoP2RB>^-g@ z z>}`lKFP?Z9IutKL>;hU?W-SfDG55olmaBSq``2j9!pRM9??>X$^-x6JFSTbF*T%pS z!Vc;*-HUkl;Ci$ZgFxdv--ta~tDTqAhE_$e93$oDvR7YI^X18zwAX#~zA}8Y+?%LV z&O=`UNQR9_65m!Wi<10}AVX;#sW)uS4vTeNGyN%E-dsXgzO2&2M*XZS^q`jQMnFP8#YAJP}#wL~Bf$yjl08u9xThc}D#V)0x+? zI9WCbgDa^cu~fKOnk0dFZswl#aKy2^SaXL2Qh9|CMWzP%!c_HJj7UG?UfhN{V ze(PpnNv|EQc;b=ftB|~NOkw(D)T|R;jA_=SFbwCtbz7DidSe=nH1r{*)5A%=*@Zeh z&?++H{!AmxOeSY~DtpK1CRP0uxT&hbG!`$P#O7nV^c?3j5|!WJ=w8AlbV(MfNO2F- zUTuK&Vfq0!{3g#+m&KjKaS4};t>tSh6vYegvke0eDvNE*5&=8jsbO$`2@=XH`E!K& zsjpI7-V^at#u#s2Jnn3UruoeQ)lI;mF`zsj;`Q9lUiM+-y+4T$J5&Sxl zt=j6sZdUek72-1lPMoI;1A{i(krW#0yvdou{<^#yF>~>s zGW&L4y9&OkHnCEAA|`_wpA?Y6GLn;l;*DLXU)ZhT&v!KB>r}7c&dVDP#T&{1#}2Mo z7|viK)cW71>*GQvnvQ`R^xyCG;>{h^3J(tK%cw1XpxR zs&)^QJ4Y`s?&!rMnE6Ct>QxXF%(ai@dB&MuWsk0UVnG7!7FYEs`gZB1%U-8q#X%jY z!ZH5O6nbA)^Aw%R>CCcKP=40}=B6fynM_q*qT^c-sLM0#dY@ra7oA&%!bB3&gLMEkH@~3wVIK=gR&#y0xT`1H+C=TbHaSfB zi^r=L1mCZ6bjf?0jhGG_lhZewbg}ZL5?z@0WH@c?N*asiSp*gr?5xMYl{xD>3;#u< z6wzY3b_q99T_v|4Kh(>=7uwH1BpH2|binT3c98aE9V5JEpyU$(z!D&?g5pnG)xWlA z`$9{ZT7Rvi;;lV(*{aF0^0_J#)o zo2B*%(MW1$HBqg9@io4q=*y>zGEDHVDvSOt5>V0e;VJm8aIVs$!Z_pU z{$4VJOfW00J#C|Cs+LW~3@MqfRNdy1`AuVhDL{lfwssw zO=T_m#^}~Y<6quTw_}+$89#@Ue*Jk7fLjqMcTi!^MU-vgUQ-^*)q(wMfcwftryDk6 zwyPDUc;7r0AJQ-)^Si3Yon+N&)Hseyg{IBZgN?=P_;mE5D0RC0OcVgGc-KAID?{o-}Y9X!3gNQ6t3s@emavCYcPHd z+l+oAY#Rm@RdY)Ud~u$KZ3_?Q^l=Z>7_P`g3YBcfCkEOZ0KVO)PY|NWoYGpaGozwH zXvmI_3t~EAKD2UaH^Ve%Rs^?u!o!jUtA)q)@f_qOZj37MT*WkBc`p`kKS7x5v$jyl z$srMHBfW27ke}M18B#CFmqryb)P;SOQYT-2Ifpq7>BmQhUgk_zh&d~vHE%EJyOJ{! z1c`2E_yV|ckrK6i!zkE5uxSfxi}gPfZgPN;z0y%~4n-HdQ<9J^y>7UHocSe$p6Owg z3Xo>d0+*S(%N;%<-rtSn=J#yRLc{f))TCG@QL>c%935~c9&$hgl6i!dW@30&-Y%p$ z$Gdytu$a)%29&tlsc(|BtEoqRpYWJhOoKu1Py(Cw5zJWMhjsQ+`Vm#*n|9C17uQa9 zmDTJX<(F1n^K(=^m>aO*)DjLbkbQAr{l$Gz9tugy&}pnNq0mm{0VJJ<2lB!O=dqi%%xm>k)RYm(8?7bY)=n zRsk=`tqf4b6Y9A!GEZ%)0Id}LJfDgkZt8;}b(2dw{1hg0UQKvZZ8rv6Q=yYH-t@|_ ze^hQ#iUyMwx(6CRlR|DE&21aL;JBJ*{xbVY%ySEgp=_5DVpdhN<+IGQd|I)PY&i%n z?}=?%b(k8bLQ2H~!Kry)e7`XiD9;LiByg^Co2TC*te(ix#yQH4Bz*`Wa(6(lCeIzM zRO#rmv@TL5oAj;5L6tEV%dH^T5-Uw!N5^=#p#K#$Bt0;}sx|J&d^GqRId{?jRC z!kLjX7Y?qpsHp8dQw^KJhq_pyIX+(tW2BD|m%lj}<}`<$lG05cpws)^)kSB|6$tnY z#Xv7qV*J%6b5ia+t#x1dM_e*x*BM$P8B#;Sci7R5;Y8|`glxht#bhY%dQlVMmRm4= zg(NNG13~CHy%9AndcC#Jh$}fNZaU^Plqktg@KyL8ZN7m zEg}n^3io@a7Y7W{$vulKIfNQnoK*+SyE=Ld=D&{Q2fp+LLD*)en)i?mZ&b)5*__KW zT319z}+khO=Bp3Q!(yUfpm_(RZD^}EU+nTj2AXQkBG5wWV>ckO!x@wzfZ=4 z3DrP0<+&qh;~k!)I(8wfBLkOme3AmF@0qfdZsi7_@5{hG9y~pHU#SMuckA0QKKZZE z*)$sQ%KIE!CCj%<%}_9#Fly1so?+joQY_ac`gCjaPDjrqw~_7@%rFKCba4El zJIt0pEdVj&P%cR3t*9(Bsk2%)p*vR?$SaNPUCWMz&2%R0dABAOpEG_NMG` z69*5-T6dA$Dn6IhOyT>3S43CiqSY9WvJ{>9u-fr={^llskB-Mj-N%SFyb)}%;boSb zPmc@Rq-Xzzm!EM_$+Rp;rD!6|9zHdnr{f<#qUD}cwC=b-7|~93l^k6tve5XM`D6J> zIX_eE){)VuPn#-DzaJ=z>$GS~?^9>1dW(< z&!ScKon?Dz^j`i69$afC{Cy%sx^F}hbp6+hNRp0P<63DltDW|YbxBv(4s^0SA1Fm% zDdDKog*QXsLvrh{%8FX?iTShV-gau)4Akt|$8p?)zDfQ?fEM_XGS@Pf^wi3en>@)n z8rG~l;&G?!E2_LEMtetVM_Tl#wEa2Ky7Qoj?cm0;xdBQFDEDEf-SU`IGVA!%@->B{ z)EpN_HjlA`V^n5FvGPROKI~Se8(_7gvdt1g);oi45wV$!QYlW3B%>1Ty>u(~Sj^HDMF|F#EQ&oYQY9b3PH{TQ_z8_xwH zlTb?;&c%kfJbZGQ7%shVq>0!s7vb-heUS`-!dIKEY|W*U&3<}7p!2aguqK|Tu(bAw z8gHVZJJrOvv8}(w{2`mkGUJ z`4!p0!Yu)0Xn?(Fb~uorsFQ7bA(zRMbDkAn+z+t3}9K z4bp2J#0I|DrPiceLb#%4%pAwY`fS(d@r@s+H#9{~g&8QKf>Vv<4y4HL< z=m@Hq5@+0TpG*%CqQAmjmvsAfsOKE+D9c8$Yv(uHSmQG+|c#H+9GA8iZPuqG7jAgYJJsFz@Q{pST9vSzDN|;>TjXyPqu){h%{NIYCLpz{Gkf_AvnuG0KgGSVR!9J}vae4n)~Gh*-1KGR0k1px3Zdd>)K-e~ z1%wS_oucFbRD^(xI0kpQHa_pb9X*>b!|+Icb*OBZkU5e1bA)p}g-C|{gZc~5)rFye zyY>s;L`VSvQRR!K4#URS&l%G%Of}bBzFae1ZsT?#-|i52K2_uN64DAUm(9lJYR*M; z#4!}T&wNRfmRLkkuw|{#%1kDg6v9Hmx|-c^|0<->fL8Zg2i@IT=ru&}?p z%3Zu`k33*&v(+(o#s)sR*39%0DAzMvpJ1%?`eCshN2wFPrI{ncAqsvz%Y9rUJd7Bz z_(dZj3@IwGQ1Y~c-;g=X+*QHq2lTpc$Lr)=+P;L=QsFTloP=^N1sR2Qk{l3f)Z8TEvz%7l2T$J|d`01&bZKuk zI}`}jP6Zwj?dbPwZ!tRR6)&W_y?Lrn&3K+)=B>{Xf?R1^6wiFtfJT7{sEltr^_jqq zg&a5EVjL?0kXbNZMcd;!97RqT9pQ1p<@5(rZ zw<@}>f=oXdT{sHyO=cer?~HMFiop~@L+ks?(LwtRRggIK)PAm-7Ul6^1Tv4$8P$)p z8_b$D814j@_WY`i2da{jP5j!7ULORY=k4hgcN)~s_3raW*&bDI)F_zpz4w|M?KbvD zuBF7Y8_YY3j)%f31XVMq>}6aemakHR!;=u2IE&i3os-OKk;QtNnh_qMgf90Z$(aC0 zwzWz_Wv8|!O>w_XarY~N`TaNtr*9Ra!s5{IPXmKHfV*vyyS&Y2?two<;Uc4u zUVeF0&2-eIl$LX9DUnhhQ33&>?v77%BrA?tiEg*=!kw)H-JfT8Q)VY6n4tM5C?=}Z zB{nT5X$BTE${cFgHyI=;iyP_uZ6W4(_$gQ~3MAccY<)0`2gi}MIYP!dbasYF4@VZ= z9w=GE=3c?sz&Gi-^hp!+NluPw$9|q(V{M=8I9pqbZ)nl_qtWX0 z&!>Vl=O%E@+gejG-6#Z9J~%hoDWcvBN@P#0Bbg6IfKT-;_C=a>0C5`4O))WISJIiG zd(_l*YnTB@el#J{PcFD!wMUGwt6DL!jeXqCxIGi&Qyzw z-FYlI!gz?^yRzQG7r*`b?$5B%)|obW2mf ze&{2v4Ovy`QJpkHqOj@wlEPPwS0*O)e5qtf1FV=t{q@GDn$1bJSGB5$y#nU67k%3B zE83w+6+|7I1k|;WnfNB&Op5X!e)n9ICF`L^nKj1J?93&2!};vcwS!liI9rZe2VUVf z4l%K69M2`#9FIzrC2DK76`K;n>RbSzVGp(Xt+vB>MXlfLH-}aJn-Zy!pMoS6Hve0Z zy`Y;x#}P}$Z}>|evJKAD$8A3+jm5Id$`dq(cm5UJq>TMc;r^UTUO~CZa#m1N3W*;y z`*w@vs$a;iyywo@+n8Ue@#g2+`670cKAmGnUhTa-?Cf{wTwyCCUe82Mw>3flymW9d zfB+vZOOrjf71c+(Q-_RePLeAEog9m=z7dj;r!z7K+rDt&_K<#j(`iJkxb*27S@ z2^-9%oc{~Ic1GxiJ!_igYy(pexb_hb(e(=5=w$HLITOs{i;u#Dj{bdMK^iOPWuBzg zz3r&UQ2HK7MA+pQD_NAyMUzpQ-XScPp^d{dqYLJ{)?;RMIcQj*e4P{H(cht+@uXC2 ziafFr7wr#Lk1+#={+$4uKR@qg!lT!Wn|W;RJCj0{x9^~1Xvu=e&d69A;9p*Pk3!-& z(&kJd;eNp$p^!jEfTrT9L^peoueazuCcn%5{{ccky}!5_DAf+jWirWV#7yk{1f zJ9^JVbUmR5kQXN3XfU@;IG!bWekRboWMboU*WCG5(IBgU0?YK_K(gb4(CrIcwbjsi z<~RBrWSRRR#F6%>cv!Xum(-7QBWZt|gMObC0x87nzXp&3e^RSYd^qRxwdv_GC}xcs z2wVH%^EpzEOh7NTUv@jiV2}nrI@ixy9@FE^Y*sikHG?*3DQ!X>Gp9!&Tc>MEi~is?Q zOflCqg)xdeujaddX0h{ytN2#>IZgw%;4#`!@4O|2Cy8I|qWFH?wHyMH>MBv@UBlb@Vp+(o z3lL9Pk6p56M)vEyYRB3bLp?6!qy(W$g~o0oc&Pkv9L`Ztx;rU?%J6@+M;z}b!+cBf z^R5CpVZK10jO86q+(E;eOq}?f8?&}BL)GdpL`XaMOd*JQH0cF%dBF=SzYm6$U^o|G z^fE#2Mtt}33^J>5dfyN8PHu6HRpk)WXVOZm3^KGT>IHf5;x)d1O`Ub?=}RIOLBu6S@sB%Pw!|m z2F2BsQ}q26ZIdltFHtUJR@@@n@;`rBde_PMRXN^4+LC!}jO6p>ez!ukUUZaoZi!-= zh`alTX&I+?E442b6f^sK%`&Fty;;Zuf9719TLz%>r8KkzU!s4-ENNXG6F`w+e_0{m z(;R)#VtI4`OPrB(rMkAE>fbT?j9Tbq40`jmH~ZSH{MoWFgK@;SXpBlHkFvR$g%mQH zpSxEp&0V2O(+vbxrof^>gX9E1MM)B$=eRF%2zlW9==J_~UlEb+ua42+txOl8cv{fOWc^BnHpZFE7pCiL-5{2}$ zq?OP`l3hlVSQ4qP(T4E;T*#X|cDTh6hN^LrV9b(sUwMZ*zn=2OgPE~QGUOIh)`>{` z>-Z0!@(Kx-+5@EHYz5s#pmFj$D=Dy*^kQDW0BO`pL1TY_8M2%tE-sw*^R!sWWZAd0 z)wn|gcN$nO;*?dOG9QIe6^$989n;Wp(@Zt}%H<-O>_m1FpxoOK*lum#oY+;oU$Io$ z&FLBR^E(3SrV%SG=1yEWN_n)@{_rdljUREt+%!GXJ zR|HraCXat-qu5(qZ`d^gd9lrH;;>Il@v6 zdzl8S>8DD_DRzB}TP!^pAKOO~+{B+EDWwmPP51s*z`B z4a}HxMzvTkaO3yROR~JYyrB#&A`8d-ghqeJWLM(mOjOB^*s;0~#+F1$2Dq9M%O@k! zRYpScM*11rdJ3+Cg-q(&%FkV5Q(1;ImB_XybageR?(SCneAJ@3^gIv5Fp3DF6@fpf z?GTn2SZ2)4$(s4Qqqx-)EFAS?k*aYoIZ@eJrhH`=Od@|C zpC%?a_-guX0FC{%%|Vv-5}l9Gb)di{-$d`yOxjS{x(Yd@YR z+?0D)DcZ8-9^uz)tG2VLD=a&T+!ZgGp{Q*KH%allB(ZgAQaJh9^=^rd&2z`O`y~{Q zBfTEDq2u1=LagR?RI4`)`2)_TTP^KERBbKq=mX$C>&~#$)Xn%5W|Er2-hY3=;nZS4Z7xMtwSuQgprm|r6H%!AEjfq46i+ncHt;*1%Ek5PS zScH77$iim3v>cVioy^1!wIqK&2och$J$Nghu(oG9fw0t{g6r$?;~1{h`-q6cqu{N> zO^RU)EHL1X!p*^;fillSl!D|E?WUzpid?CCB37LqiSgdL5PS#mgL4VPkV6IBX1~Kc zKlYDk_$Bv?R=V}#E<%IqT}TnHo=GeV&hZOet5@SSV>3ptDw!5_@!-%mo>JE>jwxPv6ddGTWV?@O)ndyJcdVqo94>ha51+MYF((mc04N z$U4&9-lEMQpfxaq&5vFN1%+zo=pfuGJQ_ID`!S*VVUILE76+;?PA=6#M1vb;D z(?|EzGJv~s8uiW}V|jlO>0PxcPHE9<>c{eu$wca~n8q30)w9hW9QG?wmV#MRlWZ9? zn0QkAm_Yy(O;Oe?Qs|YIf5OYGlhxH!#9>w`!L&00{4r?ixYA%3%=gTW1H)xe4a>cO zyE9j42@}}A2?4OyE(A;;%W=`}-E z&R-kUnOb>f$zXpbyZ~V6sMhPbX^;^&-y&*p<(UvZe{I(%{hk~jtBBdm1Iw<}3q>yA zKn6X?o~C(UZHa#{1p}$1j);FXQJtpvbcAhV-1OdU9m=xmbl9-_MR6oB)_0|++_0zv zDYYSxOlURxmZv@Km?B5({VK-9K}<+jQ(T% zz{8AIJpI@^WKSf!8E5aiK%EKa&i5W9!Us-`u0H2TmI9pv1!2Q$ME&3nuWbJhgGx!e zO0XGaKf`}65!g5NLBphcEiufU3;iBH2YoF4NO`*6O+kYm!CNdSd zfamVg?a*x``-A49lL6z)cB4l#q3Fb9=8GO)Sf`SxYEM$6FyeropD`Wu3}qccLO4^#_{uO3XY|K5I7&3Q;It5Vqtrkt*jakKrP-2TtG4Jv^T=jR)UHASj+2BxHTO7Y+ zWHgL&x!^~ObmgiZVkMfOcj6_H%vqEM@5y__DTjt{WfMeL)SkN zZoi_k_DLkoT(2xnu8`OCc zVgSR_LKiIjBvF}mB&-!1osJWDWwlrF#g(&8T~c(|pk1x1-7|x1JOjML|C>*>h5vtH zJ-d|HM5))Z)2S2f-cjZB+Vs>pG%{YVu&*uUB)UFd-IjlkQiswBwt*j;R<;@N>Iy@N zcJFJ4$kX5TvMAbu*|k!(S)>D|Q&BuijpV8oa3V8&AW$9>AcJzPbQ^kX3Se zw~k;_cYH!F2aaIZ%k6_VO_=!V5Nv-yDLDnAWUMRPD`#)LV+Ok&c(m&b0Ngj&sA344>1h) z!8{IAS?-J>P<<|s%0Lu-goNP00`1)RA6tsQPT7?mRSGb-x?qgJ_-rc*0b+l)@rWKr zkO?(yQP3VXv}92n6v55D8gF<40g|gVj)7-~mlwfbBIwtL0ZHY{3+Tr?7r|hg?aXlI z+aB7aHk;!L#9t*|l#^|UY`Yb#7x7#s3BP&g}+{8bjH!Bpvz0r-9oOk7{ey-2wCuMcQw;-qAh zRLjy?Yyd$9>|6q#KXNnwO5&CQ1u6HWyEo>82dMP2Ihb&nxm+IHRSbVhaK%uPr_Jmy z7diiB6VUN`_-@yPlEWP2PrOWx3|L|-P*+}ksC=}^d@`Vdn`WBOH^^Yp;tv&!z~=P% zUG9Y(LbRBVNWS1kH-reg7TNqtv#ReKwpK!-2V)#9jBwq|n1{!*qQS~H*M4qa-76D{ z4kgNnt#*G+mwy}TKJtI-Xgqewu2rKCI-^-Mj}R!iXHx)txjCsKv{Ov3drV=7Wdsg8 z8y0^+)@BsRE=p8r!Tw{5K)9qarT2$;Qn)P~5r86Z;Ah>7?3nc;X`pI%<~0t^Z;gpOfw352^X|r0tdH!(ke^mZX%b|(8g>BM**DRNurdR2RvjV@k9nDsMzv; z(rLMdYz5xCSa$rF3v==}%KT}*I!i`mg@sZx(hT~gXa9fb&Gmp_ErPC0| zG7XE($`^nG^52)c_K*mtRmo!0B6E`La2sdAF#&{ncG2s5g}#i%QltVZkt5@hO;ukNv}(LpQ(_ z)Hr_^cxWm39j0h;Mp>TlSYh!9!lh`lvu6;SgwHolNt%8v;3ZFTS_tF zQ@G>}w$x%<%UWS}^!b_C_=;ueZ>G-s9a6hHdexSnj=Owy)_h&G; z-<~|^tD58Gsn8So1f9uf6Isx8Bv!yKT@{4~2OK8XRHC|$>1rPntV5e?dz6%(2eAZH z0CAzfa!W87)^h78JEe3ImCv5AX&(7R;n_A<-HsOVG(>qq|B+Mr&Uh@W00rmTgtLDq zP9pPq!~QA^gvi*ySox~FL?>(0!SG|+Mh#Yp--VZVKe23Du%lf@pyy{8C6?KdSy`>0 zcA4Y^Ev$W z3t6PL%|1yn_`e$HU>=JzecEIs>0aYz)*!dHxi=K%L~GlobB8(>oIib7_AItR3F(Ea z%u`^mdhs4AP@mcQdSBy!#g?kU!MwgqouT?zG@xz-*0{It3>QJ=fW*BPB-p z_=|-MP$cO<>80MDnr?W155>WuTSoe&nmx2D@#>D}|?oi5J zj~Sdegd~0;x6ME$29h=^gr+bQvpIuW&3VR*mnV-%)`c@ej9G?HiL%ZNlgw0~Pfiwu zoNa@`a<3fTzt&lGvKxOAo&-=c1!vt+;E8Dd1VTMe`+XfxkSe+%^uQk@xq%!)w3v@b zzTjjg&wP#;O#^3s%Ni$b`K;MYG3i^+J%dN0PuzZ1Pr1<_yVxz$xEYZ<4mvvp4$%Cy ze&+yO$KQLKSR?f?RY3!+{xxA)s+)!RT$R48GED%NI4!R)PZoa-L~+`hjaVOBsGLrM zZje@I)bhKDueqvX#)6PW@%?5OAiZY=WeC^tVUD$o<1ROhBkYHqM1hiOFEw??@GLT9 z=v0=TB#r8DZ2+VLS)(vl(8G#nLr53$&HJ9sO7JpA5)> zNdHXZ|8W$Usf7fdo4;iD!JPkr%O%QaNqZO58k0lhY+?U}xaF=6FpWStcNrTvbSi42 zwa7<@r8FNgmnWe{(<^orO&uh!>oD;k_lN@q(SK6pXNG^Ewj8wKZk28CMOML9Gr!(S z6^-GCi&>sVYlW?cT8XZ`w6wBKt-3#A!j7p@)e^$tuj_;ptseJvD^kwSshk9Q1ZnAhVL(aO(A?zQkeUCqxM6` zdZi|<^*d8OO3AbvHUfOZWam#XoES37mv%g&pnrcugV1w4u9r=1sU4RyM}Xuz>?%{7 zBDKON1S`D%J;Dh*ql8#NFOmRC098+lP;P}91vr5c z{B9DK&NT!?(u2e)g;oiSez)HtomoC*Jb|s|*&s>23k|)rRE2cx29S8!lo`}&a>YF} z&?@p>yyuklZT$$b43a6kIz%{=t<4Ys0K9)(N@>Z4zU|LUn*!ec7>00OXL;U7+*4VW zR1{Y$!EkcwY@br(p7#3@{EmRf-8$tx8lZ*%*V`0b*LlbMsQu9|{n;})J3AubH;q(l zG3^!_1P|~^pxXt!kfb0AJ&4uQGnuXHdqL_4%F+b(k~;^m@p?=QXtDIaw1dm6J&s@=MIT-KC_n?21lV-%13Kt?24QIHb#n^|2zIF9xBnv1fZ+v|piZAx6b^`H z;qv(0rgbp62|)BvAiakAdv2&zLSkC9HPGc#9!}mhTRwpqdgH*%FlfnzeOrGNo!G>K zNTn5Mp%>rY-|~Qm^y7!e4N04C>*-q@etkw=fCyEMr@{R79{m})>mNpk#-_mWD z2Zc-&fx7+~PiBNG;JxdRH1E-TDiRS1S@ee+pSjhK6I4t3N}*60jAl5@w^mFPf(v+g zMM9$P;d+o1(KZwDHAy7XWF&tR!nM=?My^!o5!EyZzJ<=J$z18B{v$W#whu|#VbGvC&C3I-JzOGW@&5vaHI{Ue9Nu<4@9wIN5dUhBu))Y8bfBg&20x7!(W#VZ-@48)QglMMv4_t)8+8>eKSmZV28g`gU&fj8cFXd(BX2A|8-K zYX7@2V)pV+6JAs;%=>>O+XBVr-K`A?+2?XVR~A zoL}iORHg`4cZ{_soMPS`jrIFYcf7}&g<*AVzV}F=|B%u5Zc-o?T6QrfNlT; zGNqj)wN7z?h={6SL-=6T)kv|7fY5$9ijP{q%u0#pAIG0DZPr`3AiX1>6hs{=OwjWPwZR^QDe&~3 z@D&K1+MvJgWj|!*Y$p73mY7XU_oOyEf^MiHg%vWbQc;k~XdNxsHBKXRsxKTv_T8u* zn(UC1`xJp{=NvDrZwkIlSyX}Uj+;x_-q@BC3#5pTjn{uyAK1BZDmsbOI03_}s}>cO z>_JRlhMoYt-c2W8l@kT5!lKq4J$t0tcs81ruGx9|DEow7O3wuaMS@k(XZe*ANxARX z!LTxP5xBMg14F&Tf0rrecR3&tbzaDg6)}MwuHiin^2RRj-OEj^2wk{mI>)oP@XaqB zaHCv;8-0Ia1*5&VF=!0a@qv@UVt7wSu|FDFJ^ff9-Xit}tMNgexSdjR#C62M3QjNu z?DNsxaUb}KOIO@EuNe0VEXvRk)=)RUk%F>1Z;Vy*&rLQ)&9}g?7&(I~g(QsZzz2aR zGf}bO9ub<7zi*^6#~H(3DoWOHO2?(q-1`N{-(i1{ABm(f6(IF{;x7#Rg69cjf7au- zk_66>*^b=BPE)vZvy2OsMisOBl?p``SFGl)M^hz!_=wO$z=6k=Q4d}`J2xk{ z7f*k_*Is@huxNZO4H}}(&s{=q&bW4oM(Mt)$Nb;YBZ(>ye8+9CjxgByUG-Wy;ntuP zmHHTUBQRVM;$hLq`NtNmi#Ga?k@#hFnTW-9V>og7ybVAygwRW+SbY)m#!p5sRFRpy zT3$Btc2=|^>XjW4eAXftLf`v~oS$nJRwREB_pB{q8ZWl}@X=G9AJ2&c-hd~9v7pwy zu|LfM4Q>aA>w`Ezo%##Xo4PKGYa-y(4$?@!ZFp6>s-wS}N79t*Ni|_N%_*;Q-V|Zl>NBNPwGZrA+=NyKA0@`te_nrZ z-KdxTI5n$cTFZL=JA%@%#p#4`R#pe{z6_T={NvVhYYbgO*oM02a1|2Z-F!}^2@~eH z@SZ*PSktHX2P~APsP|V>>EY=qI2w^P`IRLq0Po!x=Xhfs0wloTr0&s7M;DZgH?Y~O z%8IhejX46eXhH?4LAE5U5~*0u7RP`3yeJQcq#;w)p8j=vJuU;{p8~c>QNLs#=)3H_ zk0-$$ZvR%r4?%Igl+xRN)ZlzDpr0O4y9{@^Q38zZ;xg4d{OBM`ZrGJF4^sGXzw=N( zcmx)R0vG$AllbfU*Lafu1T|%nEuqXRC3sym=|7*V49em!vz zz3QxD5Cf!jQNi{{Hw0D<#K&hDo&N>LVSZK$aZS7|8|qL}npTfyr*kFTaty9vtl^_} zm?8@&ulxQtRSPRnZf*OTSicx0Uw=xn8$V#ymNagtKazs@1g>hY+XHM^I(ca0m>9o$ z{WumTa@UV|bsYmLYA(=(Bzk`fakNg$)vf{7zD z`zfHub^LEh7!y+YlQk+74kK*>?9wR7D{f+^Et6tI91QOl6#4Yoi8p_p^}7Qk)=map zRHqOAjxet`nsJ%KqKD$T%mZ-2d~!@KLxDP@DKv1lQS=& zd%>7FHf$#$xHS4McG7=SpAYJujA;t^q1ibK`UFe{B+DEsDIvqNA7sM+`h zYVq5KdU)nWzQ=Tj>Xi4N4uAd^3)dKFRO17Waf&Q+5j zd{Q4fE4P=)(t<6}P9eQJ!czM%=`MHl=v%70bS=2>V&Q(~8d85Y-6jqc^Dtt>iE%^& zJ2Vre24;j}EaH-aK?tUSgt z5t}jTGI%u|N4PO zM}Lj6d5xZpRVaT*Ooh$D+%9(3=h{gX(o^gAjv?U3h$^c=r7P43MSIwrdv z0SC0mlqJu?-0H6Z3coxRyK+?A6Y50op#i)_B~5k^zFL1O2rOi&-_#&QRlQ6_Sq#^F zCkrG|>Unqk>A`2_k;IZ#4T*Sz$Z&Vg!y()2YhcLZU@cHFAWZQxjGRIz%GUs&g?z?U zu58q)HUDM+r)z>~F-O1dtiCOAA5OCA(e38WbWOggBgD?AXvNj-&(V9Tc0(gz>5uLXtCJnwaiPf{6RdzYO0z>#?-~Jg}v@{&!U9l zPRn|sV3hv88$p;48wBr9k0A~fu@;MUdj{PdNEB44CNzpl(zVWABZAlapwkT5crr2= zC8B<-K0?M3^Bw6%6ir$SMsI=Ah=RZpD}j)?C&_$CP%_nSZNUhz+}nnI zTBLT9L38Bij{_W-CVQAoPw`XBW^exraWl-`)PJ4vj@s3?17Txi3^QgTErfmq@Wb?s zF?+0v@uIjZXhuSpbdF%lQzoIiYnlMNhUO>B+ymPc%z4!*&q<g4AgKk#_b9x%Zg&D(xi(1FaYlRRksK^F z1Nj6j>RLu8_#f)5XQ&sO8H-_pegXW~Tj+nnk|U8kyG<+gS*+WQn@E8T`O%)WP+^_!7ctsK~5-S`-H33qreBzB?`E3)qAZ^;<#eVsLFK zhdg-Ch4}qaP~*|GRHG&j1Rh!~{KIepm8pQiB^n`BER3#6KDDZ-rX3fRlF|OrR=I!u z{keJRakboq>!T|~i#tOT2IM8}SQCxz*HbkLrgDe(Bj#8<8%Teg))L-;xColIGhV5&KjjFP69@U!Rw5$9Zn$F8 z&fd?}v#T~%nDTKhmWXQ889^An8Ewm&AM{ieARfldNt?%geSys2mG@)4Tu5KiYYB~C z#MegFx>9K1WMp%~M>%2^F#@y9V@jv@Ic$s13pkqzTj%^M`F{Oo8>v&^UA%wRDZ5RN z1ZW<^MHJQ~9H)*Fp1F3z@98Q;vd-`QNHB%VkAr&aq{@mPbU9XQrFuBSHSF78oi&@C zC=XY!khLBuH37uuaKaLUYez0&*@W`rji4Q^9Sx+?Emn|`WY~E~HpA@|>phnxr_n>f zxgMkF6#W)BR6?T=p6V~a$Mb)?Y6{Tsrf_k3F$RxnTv(i9a-aoO&Yb=Bo0M#C>oO_r z&1DaK)eDV>rUlcQDpZ)Bx6U8zoiRbIM*gE)B_nVR8!hQEDX#A5u_)*Cy-bA77d9_G zJ!vL2oiH~Mkd@ach8v+||NBD#yWuZBmc^;Mt7$1dywOYLYxA*5C*gm))ifA%{#>N( zX%Yf1*QEt8eeR(|;G+WkSmliE2qwzn^$zuh{sC2tp(~J^g;bis7KtI&j_C#`F*=lO zZgBy`!eMOvh+Da8~XR zkM*|fDL<4j0JU~tX!ST0{qn4`4J#ugRRJpdoD&LCVBk?D%KfFp^|`QXY(kWT^`&wj z4va+Eft=1?6SIHDsK>q1z;$X0$@9OMu-EdM+)D%zJ?b1FCO^2?@YLP_=~K>p?4R)w%7F$-cZ zcX8SYu9U<1%62>uZ6vYAi3;11o2nL<81HsO`^f~wJ6?Zj7QgngL4EJL{zjp#juAS8 zrgh!MR;s@;x*PyHOm1-ofYq&F=vCHbMRj3bN_SC1rT#P#NxwP&r{|Z&Da@47toP5Qzpdd=SB<} z-DiY-X8Wy?H&S;>jJvl2q=#SAUlmRW%rhdyg#P!|BN$t-l8LPRl6#0=ALy!_#&*It zL{753i5c9&v{4^Qm2*%!Q*8uytM3d2&?r4sz!ZEJ#JKaVCQ?nT`KF=wZGPWNfeG4ligkI#dK0z%Q9rk7Hsm=! zh4J1YXcVYsCo%{$p2ZT6&a1^eliq}6U3CGj*F6xCWbhWNv|Nzp(Z}cP6T_<_`#(PR z;Tfy1$SURG({bXDlp2zl^vCvjyvxIBlemBG(e+&2lbeKO{Fp9;*SM-ntw(5=17PXN zkdyhz_3|TXuTcc>Lpi+r(JAh9h=pFB?fekPk-jPCLQ`GJ@Lj$Z7}IXsU5`NfT2y}A zMgKoUHd#-#x#ZB6^xl(!PW5B1RWt59fBp+1zhb1b;8CT0Rx($UWww*qTjVS)|ni z@L|&<57xfYda((5sTe+T0_smE?yZ3Wd%{l{cG; zeOYjd9obFbDoXx7Ti_ONu51ri27!Nj9J1Zl=88v|zTqUuujt*!*RS`}l)tx2la%#T zhhA@k>5)rKh!)XqVyX;TJpl)KjJh1am9%5*h>lYaT_r!1Qe;;MeEIBlkPV*^V}43Q z&;9cUJ`|5_*eVLuq?RtSH*IvZKG#QeE1%{spW+JEDlX<2m(^t*ZDG@E4asZ9-C4^#k*M(0_uO9cXX3$B3k7kZ%9Ctt~R6I1CM?`0SF``e` zVfkx5q;elu*?`1C`3=b{m=bdmeFhytJ~#S}2+;I%2ml4UkOvu=kTgERQ{IC>%Ys2b z3-ms|gZwXj1E0P5j5_77mU@4A&M_SxrLi3W$JKut&TS?bT}NhY6GEM(i%3&(C2^ZS zxgW##DSBJkB*RRH{tWj;cYA6#-7by<55Y=RiCmQZQMHqJz(SjF}4 z7Nx^11D=s4o!r(JCRX2+`QxC7(qGVVZ8{K8AqQm_yEO`-W8=lh?^=JL^H6V0YyLFD zV;C)He-??OsG(e1n4g1luYh?1f4AtT1GA3|7C(vpXC}G;jbwQ&G2r!1l^B2}wlhAj^WCVznUo^Q%Y>Iw*l8^g;gJfm z&}3ecsj?9H5g5_E3Wq2?X|giq`L}5&f=gz7Ur$08pUuX4A9v6*Zad=WXEEp9E+V9b zj(%iI09xjwO@)S+h8B3!aBUIKz2Ia^UGO&cE4+A$X5feekcp!F{ zG}KpI60d4Esdq8cwg!K6bml~)?>e&|z5guMe%z^l1E~9eFsTX1$bc!E)IhtMgBWU>HxN-8 zkwoc@=#&dg*Fr(X8TJ}wiCWQlf+4NQaR~HmO8`^#Ba)(0wPIixy>n(*GR+7P|J1^? z+{s)~Ak_Li4*%c|6f7o@C6VrHF4|b9k;>sbzma3H^&Ec-ovbdl5zNsV&}^0M%HwAh zb-UVWZGfj|1#5!mG8I-(y55AR&=GY@P0( z_V}%qrs#fOrMdesG+S{FS{x+ON7D~6ruPB$ry|KKuSsmWZN9MjIXZ{4X+$|e6|xl; zufc$L;rxHc@3ubjHsFlAZgg89AX+fEET#7F{`~54GEH3BMvR6x8pmWN%bb%{QU-67 zvl`^>=p)(##z-$p#|h)DgJm1&Hs}%MfSyvUe&TF_g5IecGA94C5;&z=`RlsI24RKO zY0>b{i#Wvi|Ei8)`ts#sSTb#tR&lAE6A}zqqpyG2ZZK-Q9Bz^jO&-h+sd*};1BEcn z&_fv=7wK@r!DeC#P-@8WM!0UnhFy!UKVXb*6k@3UDv6;Dms%0_=XYf(4*L5XgW&vI z*icH9E@@({w3?1n=2lhHH4t97jDf)KEY`@Id+LR;*pajY1-#Tg96+G77Qd0JltGFD z15$s6RhuoX1Eaqr+`4qNU*l+&v?$u~27 zf{Mv|=K!gQXC$eJ-;IyoOTx@adqSO3PcJ!@U z`7Nc|yoGD&n86-*(Accc<_Pk8)1^meER@?9y3V7XQZRkmcmJ1uH$g8O*3?Mp;kAD% zq*s&dhylhLP%WrwAsu}LBf1mGX`x_2K`#0EFxX|c17{;JVW^-HO_b~4ac z78qfDN`7;v>dNnL&I7EzoX-aXFhsy7Fg!Snz!FZkpP2f*6=!9{ss<~=r$m4I-&P+V z7(|8wnhet~ck4?<1EdjvfSKh?3iot>m8vjF#KsRjyJ7>An(64rO2bFul5 zV(h~F-2<7nCm^0;l+~fff9!aL1%TCcRJV!}8GT)OCx;U*d-0W!R0L>=S=9?ZdxN<$ zB^xxXR$QZJ!60GuCJLDXCRk-PnHKvg}~pkv^x#wUb)j4HM~dVhqwlU zV>_jSaeZ%(wWlv?q~?F|7(mT?bv=o> z5txp(N!edeLW4Y?bQo+>eb2|&eDfB!e{=7U0D3FgA2awoAxjPz1p`JREe%n;>fv1L zf<+rQ#V!3?xw!n=9Nl-M>j0t_LGrZB$|4*(s(<4e5h*n&-DjBytXZc-?DpV_2>Yf6 z&+$q2L+oIFX@_?4SdxD`B5gekOOZVa8tN2iG?Go@av z%Nd?JGb@Xt*M(NxL_}k{Fl7d1V+cLq%A&~TxO3*YE&cKr1qFFED7n`$L?7}9XV!nW z*?-!JRHBr(yVZXstgRw#2&KV9I~xLY{nXAY!gnhMP$YckDg8cB?|Oq zDVnu-O*8vDrE8=#Y+e~V5xVh%bl&NWWa)iVa>`AEi(iLFMU`PzphnWDq7iJ(;7~;X zhwQ{mE^-WeR_ z9XHJdwFn+#$AhIbJ?JU@F(|0e-?t=J#QLZ#8{@WO@&x8q_F2Z1jK0H-g}tNQJC*9Q zutB^(QsQOK(5d;k+jIyf#t3~B6eLG#lCl$}pX&FxF^p}Mv}z6}wyKoUf4 zZ&BTH6Sn&^XXPyuUd*1jm=*Uzb1hi*~ z2f%-sw8&wziRvwIEb5j>3{!6PMOja2a8t=vAN-uvV^Y7h&CD1v;K3Knwc)OL|@O^t9Jq0&h!S@#-^8sd2?H zx1C*^U}%|#P50piIpEfss3aP#1=j^Ti>`Y02uc;gFTw^xdJWd0UnefCUC{_{~VMO2CNE1<#NA*3qmv=Yr#C#dh z=4{0x*q#rD7nbowG$w$|6vmRf;Y17`f{nJ{vq?rk837oMjs#w%^@7pgvxo2LqC=5K~Y+hS?eZ)JIDJhu0N*dBTF{s?iuXk(Dr zVS(wH3w*k&`X!m=b{Nh2z2kr497FSX-_#cBAU>~MCxL@tqpx1`Ocmod=S|dVS9ns! zV=;sVRLkR#9ie3ywyz!_EzD?o|Mlcp`iOqdRt5i)sai(FN}KeDVA?2%2mlcW=DOV} z`IRlagRRn~jJtIJ+2Y=Uz27p&TvoO8(RoO#Bw-KA_&)I=^1{f>AKHIIp!+;B3iLq& zuPjvb8rRAcL9iJM8jP#1hOwde{d#OUo3GO+4iR0*qe$eqeVd|oEtTZX>4<>syt;jE z;U&|^g0?JiL8-%lQ^<2g6M$Km4HG#H4p&tFr7`qcf?rh5+x|$_VS}9%8YPpi&N`#@ zvVhEz8NotFbNMHs1WSLs&9z*jtDT$mJ{O$xY%{x|8hKK*l1c%$KxRje=3~d>z|xKA zs54%>Etob6(n~XsI(YLSPGhkyMrHz2hrmJGqBK9ikZ6VTy^Pu(M|2e!42r^k6fhB* z$R+S{%_|u)+G2pBPc0V)Uid<9x#*zI+xWFMVAkpdcOEf>-o$^ev)j1y%4PRBS|QH~ zZm{C7Heej}Xn(;n4$K3sk}cN)Y(&%EmZitUocKRGq;V5n9jy(Ou-RExi!#k^8(rwx zXsn9iL#7}Z1O`Ap$`otU^q*RQ9aTJ%B9Yxi&dU~qAmeD9+QJ-8{AO$E<1cx;B)$qL6t+Kdpd<|3HnO+ zthGpVX~dYSmzg6;6FyVD3!9%P9qJv7qEj#d&gFFwBgMV436DGvIb|4gWYcq{TSXqr z@x;%CPe1B@hOH=QJCqZ>o_y2kx{m)B)7}MVf#1Pn)-ZoM9ts8(%va!D@t*3tiY?Wt&w?UXxbLIZ_Y`e%efmZ;_; z0-YsgHt@HJb!gdlxlwUXERvE6;n~4bA!feEx8y*OA+Y{c*!!AssyY=0yO(wqtBy6% zOs%TQc^iL=uOp#`vPrYtIa*#_OmY#=%WqsRVz^1%a#L-driUQw`>ENBwSYc=R_0ts zu`iO&g!=bCDhTmK+a{HFtx)b!qLRL3X2wtRrASB^4PwY7f2IQ{H1i-Gmj>MtkaC}! zC&zAmqIOYAFrR*--N&iOPvA-T2|T@Cy+f{E1}uNFk24WAWc|C|N@{VnwII_rTvzR^ z71hp5%7-Q+P@gQLaltf`umYnlXwTLbT`k(5JZ4ev_O?*vnTnp%0TQb3C<+sC0hB;p z8jcY}IdBl%fy64mnmo^s_-iWXcYfw;ma#}t`LobY4NR;~DDN}&7%vM%+x`!eF@crw z2k(Cp6Zv1#7uTQUY+dbW!!60@z1Vgl8>erb!6nV9QgeLGgA4-%A?(I}{01Fa_|#)M?uc!=x?o zZJl{MRE-+P9s4d@_KbabX=cutX6#!=6QYnM5!q6VWe8!ilPJa*qU;TlEFopfk`kh_ zMV6vu389eXj^2Ck+toe)oO7Pf^L(G*@Ao<9-*X;L8cVNl$0~o-vP;{{%TKXNH@R;O z*vE|<63FR?)`dH)OdCeizG?Dh#te&ETx8PBvu+J8slEg=(bHZs4WbO&qp^?6L**Na zG3C`ps1}U;*B&)<&hWMd;U>M7wx9Yc?{AeTTC@cyA$cs7f~GxRCTgJP{_ul_;q#2@~Q? zwyGR%WBYuSgd1`}7eC|)uuB;FJ)kHg+e(ByQ(2sTWTIe2{qs~f+T(zcs1P7g;TPdM zK_(;+Nsp9(!?bgV-s`KQS6mdUD^|14E-C~t6LGFdH<%_|k0p@9qt`zu}Ta5s6qbbuG0 zlwcv-Y2fVGtKsa~dx>WZ>_gpxC~TX$vBMhc5>ax{lEvZjS&joX4_wy!`1 z1E3+W)UI9(aiK6=H_a7YKHT#eJt+9isQTXJs)e&P`KOU{@$&#LO`Zi0jrh1k6b8>maN1^0e^iZ5awi8femAEEh(nbii@;_ne zg%P8HF9X&;ZW-*@uoDB7MByEMUv8UO_=+f@@>xwDwy(x-lq7969KNj)G_#=Bh-I52 z<{hVuM@repmhs&uHc)1p%%$@#jNOrcKhq=97R4T0vwUm~3DYW~lrcKWZniTDb&|So zkY^7Nr+)_Ai_z&D;kxl;wuoTaRq!a6McE%apOP^P%C2{o6lPbKNq#-sV(H^jW$&<| zaeYfFUf6%C^Gq>c`?&h$6w4FPvn4XpgVrs++O6dyyfe%u`8mp{EB50O@8T2_D?Y2e zYjI|}+0lf#clTgcuAx^bWjYV1YIXNq)43v@_qOeM%FLsX3LojL43Qjy|FgztVM=ut z7fKFHb7uxmz1Il4ueCw4TZgl#S~&R7_S9of6`$NSlVrhwn|TvKoeD!un|A4c!9x>@ zxwKy=g={CydfM~UUH7{w#Cz#-|OPv+dv16Mv$};>01bMNE#{ z)U8t1u!C@2hs#fKAD=Av5bTddd_M4D6gl&{2PoZoqAn%cvJ+(Ka@qApEU-$ajXcQM z1dBVFD5&+3o1fR7wkyl+R?qcP2Rpk1OD?N5+2T@9pHwr5swoHAT4m*-6_D>t_1=a{ zy6Mb$QfL_5OHbp@tEMMgVe>oBlsgw(jr_T=kX^nrRh2HCI~t!c5=&^qaPyv~UI;i#KZh(&t}%9U$Flq3y0Fymb;#^j5?r zkJLz)Y$x|9^1Eo=zmP#Nvl;g`=OhX#DfiUP-&2;cp>gT3@z@&lkJ2a$P;to>;GX`KtCYg-uy@AU@TP8C=(w zo%AU@TS};tyB8yXyG?hHXr+Our;eq`K2VEUKcDs$;a0FBFPqL-@|b>Daju=#ml?K-`7~zXO42~%vG|5v$a5z9)|hOs zT<^>ETyJ!`%h)AL4WNtHm=0EG_Dh%x9&a%*aPv;<+8th9JX5vP8QZ~yjW_{t-$`06 z{d6@$p{;&uYM5?LJ>Sw()2tPU4Gpht?=(1s$Yeuu~^y;MCqQt8cO9cGpX1{C|Z0&^w=295NH0AWV)wR#&D zv7hQSW4@ea&~M+&Ox=##Ni-V)lxwf&irnjF$Es{g4MO+?k~>p}haTOSvv0>==5rQp zra_+^V{&NTA-qT{9X;Vzh~@fZh2_{t#fBYvV_=FKzWF{bs1hA6k$tyyfIVDQ9zWft z^c4TX+v(0&vrT1_xcPmFR+*^$Nf}|~snq2m;v>%Reg{X9i?>7NZ2Nh6?3~;hb+)n1 zIjVdz$dfv}6&;e#_+ZZ2e3n})rIqJ2vRq?n2lOsCd?5&B-(BGhJ^W?RRFc7Xhs!)~ z1L-qZm-LP$`A7WIA0bAofmb!{R{WnxA!Ub!QcBi2F{BJ(%#^K&}+GMMI(Dv|uEex$MGDLc6}1)g~T*H4gL6$Jr%`660b*?-uo%Xd3zH=wN+^+>bpy z)tE`+Wfxt)#891yYxcO*+oZOZxMtcC-TdSp$+?~&JX7x5*Oc{+IBMEZqqye!aXEr{ zxYt+wU8@I2bi)kInD@yUZ?&Z5+yalA!i|j98TU3)-9E~9;k67wpNoQgJNE-eogO;( zeV8p3V=-$nZg>s3w&-(Q;4$rC_gYo&Gn7r8z&2OE@<^7)uQD`!CI)T&x6`L~tV~;L ztd>y&gp{aZ?CD&cV-c0kzIHqd@%Sn!hEb!?ui$#+!!zIP-gj4579yIi^oGKA& z>WFtUlHZY5p53_5++m4x+p?>*Q!ko+B6k`$#~95XK(~yFcvQ`UmW&TpRTI$e%V)=@ z*DCS zcZGBEUAct_yGj*90z|bf(wA)iG-Bnv>=+C9SO+l zBYxTVO&)?exMx#CgK7ldFC2$*d3t)+skh6fOtVz_s(zxsEVW__l|#R zi2u`vfFGgexvvk7ri$KI17YN_1VD|dB{>1$C3htNd<0NM<8 delta 3152 zcmZ{jc{tQ-8^@b*7zR_#G+Hc88YMZ*Z#FX~lQ7mgmXtOXk|n~>*kV#7OU9NdOHI}= zZ=sM-c178S$k+=xwki8Xc}KnPxjN2up1>nI&lpEc;_##>S zL7(Xki?d22mH$H0F$_Yi0bot7p@`PyT$LS%L z5=HMABw!H=-i3^LE#@~z7@TR|U809H;~1^_?k_Kh7Q-I-5Np(JjQB3^vESr8F}k_5 zqY}|u02gP-WGjzB1mP?nq{LhbG-Ifj{y8j=rlZHLOQc^V*Uz4*qzv;d_pR}#Ho*68In~;NPwc_-V zR&mNPh4B!9ssf*K;!?Dtj+20!mgSdu%jnoM(T2sv4fiJ(Shd=!*GD)Xa%fIsWtHuS zQHRYN_1>z9qJi4Gpb}_x2z2ir%g$wccW2|(XS~JnvIsutC`8_oH3nYhTa=xdrCP~O(5d)vg$6>=4}!+D@qXivQi!4n>@Lqn(k<}V_#7K z1}mCvxg#<_(?<7_G>@fTaT0H0Vy&uFb7nl&B6K;ddVDb6j7^|_Z9FTjdz~Fx4C2FX z@NWe6cS6b(WgL4nbM6J|O`JO``_R-ivW$A`fm+tNtc1}7roejTQG0RbC6Kx- z`>sE*Cn6%)!u)jCA&_Wf+|7G!i^PHu353yOjmAE-U}B4F^B?#bW@CmGFL#OcuroOJ zl!{1z4jOK!;5v2crgMH9m9rtuB5LeZ=F*ou#JlFuC7g#7Zb>tVD?5Dx@sCy0Xvhz5 zM?5SU+9JEUImeHdhE+2n-`}%v(rl@)2x3_n4r?v7T;l?PkGe_v4+~ds1C#qc)QYFMKTB44|pC7TYtBu=E<= zc5cu%rKLOaaE#xr4~~UgKV^fVUAe8EiPq^U^15e}JJceYbWh$6^`gk=v=pqCeWpxx z_hz_A>x0xJv}x)twS6(grOdv9PHdT~gsgHQOz=s4|F)ds=^*}%k$w;O$OlIWK5%n4 z^}JoZ0nu;ezG;%azt1yIt4jHW4j6u`o;s`t^9cZQrmv+{d{JaAIiKHsDt%aIZbR?# zL3%s4muwsczeq1YWcc>)-*y-M$Cz|ErYBr~kF!eCP;0t3G%A9X2?=Bx{y4@Pp zgm(LvL!7-0@up0@HBYOMQp%rH%rcMX0a1OEs3uR#&G1ITut*f+>ZeEb(@S}mp6!id zOsPEo__Dk3h~JBGXUb*eueT~Kzk&>of%jM)juDX!vQ4MiM}F>GjY$hv#4?SD^X4HU zkFF&eP`K7@p45r7#8%UmpoR)_j{ipM0|mutR(b5HjwyoofRgzn`wQTA9I6mBYvj=g0_dnbRz}Y2Ao3CUuh@#02`5VWh4RfLC2y7q$bQQR#~V+#`w3lq|X303|rVVm=TM|=;5 zST$yTS3GBk7lwXUPZ?vTz!bHNBoMUg;dJe`C>DBdQETVCufMbWyPo6jp#k{_V&U6K z_1*)7YpI(#*>Yy_17~US02&D&2eCrN7NiMZL5M=fB= z;hHt*U53xIlv~O})hKWHt&&|P3$_c16oNiQGp4Mg;LnG0VWB2m{bnad41*!%>v6nP zbxLbi;@rNdQp7V`@UG9V>xSEnyn#Z*Y{%@mSzCK?sQu*&o?di&8`rhU&DI$TU{M&< zT7ku4;#$NNK%}SB;&3ns%YGFDQj5b#21Brb&|(2)yT^j$_uQPYhDzU`hPwaTO2L}+5!gQ*Uakdrfnd95P^Lu7Nh~j65?3# zGVB+zAnD@F?VzA5AlzR>f-r%(u#p|g0wLh|^;BXpmj!}@!C2B-M>ZHs{()q7ih-r$ z65}O{g-``Akj5`O;!hq40*0@(9fpAAc3@Ee`3G8EA3$Mo|0{7QtOolw1U$1@{Rj>P z5P!-?@F+l|*k2qhAp!=7>_~C&w>aq-I6FlgJkAc51Z#mX>?%ocAQVFe$N&*XB4F6Y zQeasSz`iI2R&&H)NE*WFtS!`K6!b?yCSyrL{5lPX|ECX65?<6g9|0pwyk5hySERt2 zU?Pr9k_PW#5&_}t)_I7+1nV^tVVw`aU;qs1?>?c1Bd*IK)QD>%`r95DERjH1)0n>; z7Nk|@5MErUn znGCF-0)|KiG}Z?slJS^zy}~Jwaclt$EKdMX>K1NpUZ@|x2BmIv`Jx*N_t&@Z^s=FQ V`OxhzLV>mUYCx5gXr}tm{{jrhNo@cC diff --git a/Examples/Cosmo/plot_rho_mean.pdf b/Examples/Cosmo/plot_rho_mean.pdf index 950b3de94c119667bf7c6a11f1581f18b0d40dab..a43c239bc56c304cc949ca74615bde023852accf 100644 GIT binary patch delta 30 lcmdn^oMYp2j)oS-Ellc$T!sdQhK5E)2FAwQ4GozZm;j~N2 m_center_tag; + double m_rad; + double m_rho_mean; + + public: + CosmoHamTaggingCriterion(double dx, + std::array center_tag, + double rad, double rho_mean) + : m_dx(dx), m_center_tag(center_tag), m_rad(rad), + m_rho_mean(rho_mean){}; + + template void compute(Cell current_cell) const + { + auto Ham_abs_sum = current_cell.load_vars(c_Ham_abs_sum); + auto sqrt_gamma = current_cell.load_vars(c_sqrt_gamma); + + const Coordinates coords(current_cell, m_dx, m_center_tag); + const data_t r = coords.get_radius(); + + // Divide Ham_abs_sum by rho_mean + data_t criterion = Ham_abs_sum / m_rho_mean * sqrt_gamma * m_dx; + auto regrid = simd_compare_gt(r, m_rad); + + // data_t criterion = 0.0; + // auto regrid = simd_compare_lt(r, m_rad); + + criterion = simd_conditional(regrid, 0.0, criterion); + + // Write back into the flattened Chombo box + current_cell.store_vars(criterion, 0); + } +}; + +#endif /* COSMOHAMTAGGINGCRITERION_HPP_ */