From 69b63727bc7108217e0ea1b93beda692f126e068 Mon Sep 17 00:00:00 2001 From: Jie Cheng Date: Fri, 27 Apr 2018 00:19:02 -0400 Subject: [PATCH] Parallelize the code --- time_dependent_navier_stokes/CMakeLists.txt | 17 +- time_dependent_navier_stokes/Readme.md | 9 +- time_dependent_navier_stokes/doc/Re100.png | Bin 63099 -> 0 bytes time_dependent_navier_stokes/doc/builds-on | 2 +- time_dependent_navier_stokes/doc/dependencies | 1 + time_dependent_navier_stokes/doc/scaling.png | Bin 0 -> 12247 bytes time_dependent_navier_stokes/doc/tooltip | 2 +- .../time_dependent_navier_stokes.cc | 1346 +++++++++-------- 8 files changed, 774 insertions(+), 603 deletions(-) delete mode 100644 time_dependent_navier_stokes/doc/Re100.png create mode 100644 time_dependent_navier_stokes/doc/scaling.png diff --git a/time_dependent_navier_stokes/CMakeLists.txt b/time_dependent_navier_stokes/CMakeLists.txt index 8eac2e3..7c58eee 100644 --- a/time_dependent_navier_stokes/CMakeLists.txt +++ b/time_dependent_navier_stokes/CMakeLists.txt @@ -37,13 +37,22 @@ ENDIF() # # Are all dependencies fulfilled? # -IF(NOT DEAL_II_WITH_UMFPACK) # keep in one line +IF(NOT DEAL_II_WITH_PETSC OR NOT DEAL_II_WITH_P4EST OR DEAL_II_PETSC_WITH_COMPLEX) # keep in one line MESSAGE(FATAL_ERROR " Error! This tutorial requires a deal.II library that was configured with the following options: - DEAL_II_WITH_UMFPACK = ON + DEAL_II_WITH_PETSC = ON + DEAL_II_PETSC_WITH_COMPLEX = OFF + DEAL_II_WITH_P4EST = ON +or + DEAL_II_WITH_TRILINOS = ON + DEAL_II_WITH_P4EST = ON However, the deal.II library found at ${DEAL_II_PATH} was configured with these options - DEAL_II_WITH_UMFPACK = ${DEAL_II_WITH_UMFPACK} -which conflict with the requirements." + DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC} + DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX} + DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST} + DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS} +which conflict with the requirements. +One or both of the aforementioned combinations of prerequisites are not met by your installation, but at least one is required for this tutorial step." ) ENDIF() diff --git a/time_dependent_navier_stokes/Readme.md b/time_dependent_navier_stokes/Readme.md index d592c5e..2e6e5c3 100644 --- a/time_dependent_navier_stokes/Readme.md +++ b/time_dependent_navier_stokes/Readme.md @@ -127,9 +127,12 @@ For more information about preconditioning incompressible Navier-Stokes equation to [1] and [2]. #### Test case #### -We test the code with a classical benchmark case, flow past a cylinder, with Reynold's -number 100. The geometry setup of the case can be found on -[this webpage](http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html). +We test the code with a classical benchmark case, flow past a cylinder, in both 2D and 3D. +The geometry setup of the case can be found on +[this webpage](http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow.html). +The video shows the 2D flow when $Re = 100$, where mesh refinement is periodically performed. +To test the parallel scaling, a 3D case with 1009804 degrees of freedom was ran for 10 time steps on different +number of (Xeon E5-2560) processors, results are shown in the graph. ### Acknowledgements ### Thanks go to Wolfgang Bangerth, Timo Heister and Martin Kronbichler for their helpful discussions diff --git a/time_dependent_navier_stokes/doc/Re100.png b/time_dependent_navier_stokes/doc/Re100.png deleted file mode 100644 index 77027e38fc34a9e7f54bf9f56c461ac517b4485c..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 63099 zcmeFYWn5Hm)HQsDkd%;A8VQvSsi7o98Uc~+5DCelVWg21lr9nJZWy{7rBfONhK8Zx zJ^b(Yem}k6p0AHT)R{BGnd|Jc_F8N2eNC8}irhooC%6CrJd}SUqYeOQZ~%ZTVLt$m zy!HgZSM&_e4hoXaqRt)vdEBZM$yU|r!db8>Xl zKjl;N>#N`J}Uz8>jP=)x%WmsWE#0DvfPxN2g@Dh4$D%VrscjOc>a{yjrX_z zyy9RffY5ZU_iOwVi)y;)n$ zwXxUR#EG0wN1JAg_FMT#X10ID{Xwln&$F| zoKruP!pQq;ZPV2r!ttn&b(QY!9IMgLr(m$Rx3z;%B7*MbpB$xxulHR#>8leyeiEY2 zw`Eo%wEB)+WW9ZTuReZ!F)Vtl)OdTjbM2!%7L$+=0-mAgx>|uDx%M*^7Ik;GS3I`U z+}SxfS;fW0sy7eDFRUuOn`SDl%#POjrfck|HU^Ud#>X}BsRZy_uGhI}HXW8)A%iJ` zgg?}&}1$E2i$g8<4J5q1ibL5mhcQRXNm z`5SYGUSE5s`rTffgQx(2O3R`4=eenjH^@52Rc^OE{dKp2Zku~-pRa@;{`K1bk(Gtk zW=!<+b!0nuxpF<5;LXzTJ}oqoRZzfHP*gm;zpr63Ql3SEi{2!BR>LySNTobIMeZ@; zeYuqzFEC~3#W6ZMTIO|TukXGb^hQ;+$F2D|?Qd{z+h7tO1VFg$-%8z3V)qXWEIZYE z+SQtN;gXV*uTU+XKk(y_kf3{S=)ru8s%<&Ea+s^JyN`)|^MxhnEwB`?Xa70JTD#n| zt6T!dpVGJ&!h{_#IjP0Z&#!GbFfgD|Q(9QqF=LexIXtW~?|Ze|J2ccqcYA=P@4590 z{54EY1d_K-pHnGAoW^FdXgOo9!U9iMPj4B9I4t$M!J{y6AqFxlDjwyg!F1KJ%tq41 z+(FPebdkL`1`{d|Jh_IdSddPor3tvQXV0IPSNZ+>5dxyv@od5E>66Iz{j;Y_AJUBOWi5t?gXR!~Fbwc?rLchl0jn_SP-=<2kHG3Nn(1 z{{R=hGg;iRm75WDalE;jlVaaKS**vgd>2V2xNfpDcCD!Bl0& zZD^ub$jW;vI{aShjTy(Hj0TKPPMQv<3Ws)?)}8(JBE}6mI2luY%FNt#=%bjx+2xOk zcMqNWy!?kmcHS?|tQ+6^vBW|5H@rQrp=uEc9^9d*#LF}nY-0+LvT|3`Liq7YR7!YDaV`gT4 z#=)_^Ptxfa`I<;Nbu-;d1)cupsk)(I35RwW{LZ_F#_a=8)!(mhtLx>35O7+1V+Ued zSCn=Nlt*6>z1FnRayB+LWA{(5<`zG$udnOM$_AmvGHx}&gW)}^x>{QCu+20#mR%dB zN6U)?c?AWdXrJLIFl~;ew2F0qL4LoEeDcc`2GgBAt6sF#FDNKr&#+o-`Z)38cD6{T zk{_78&zM_X{wEV(zI>_r@)Hy%elkVM>|=4kerAK@sE@fZL4C5!hU;JSg^zU2nSL5gclLql`@5}0id!n+u$ zIDEmyH3fXQM{p^PxVJVQ2HLEF+d`L-*_8Rs#YBUEPm^oKVtsqT zD5&2K&)l8BRHt!*q~aGfuDqs<{F+uz!G8_VGn-B)^^$Y8S}K96o2xVEtV)W2GW2+m zFFH>Z)F;o$8YOcf>x3JS_AD#8G47aGVoI5 zsL~MY*yfK(K#^a1Bd!i|&%wgJZ^9k&?CR~?x1YtlP})%j`uc;%r2zfQmt-Js@%bLb zWt@28AL@I0dj6iA3L%jU8fMms=>F4RRaN6(L%ye_MRMrXXh(@z4pPb?99O$fw$?d0hmi_PIQ0A32_>ZqB8xw*D4%Cm`ld6Hr}K+VR6E!_vD z1EyFvLts=D9n8NiQ~Cn|#`EXTN4>}1FtsmI+-X&R58A+4aSy>*It`78%pC*B)H(m-Jc*W*ujL(<$uqb@Hm#ZOsTS+RUu?%t@pl6Q10uTsMIid_Fq zPC;=zHPGKrZ6})>PapOE4e+RwBD=AX8VW-?G&cKO!EVLxGs%7JtK0!F9)kiKXv-*wsYq-1UX*-dU!5Y&;&=0h3FfqR9KJZ)LWh<1gM!#5<#J zM@LTgd6%4{gOO7t!K+oI{TUOFe5vK`s^!Aw>zmF$a*uz?MRzaW-MaC=e97zz5u+t6 zX}#Eqa7)=&?>V?=q%*qr38}u`w0#krCziSj>|j$jr<%mZatfGeI_-+yoRG z_v^zR50jIhyK=6s)#GXL~AINq*!am~wSPPXI?Ht5OX~^4nZkAyU?yM1Kpl>{|6Ts40Qpu*t`RiN@Yvoi}!KLVqlK-IY{_>?jyHp z$D1Q^&dz+zXHEM*e*Wx;otvh+cSYK@_A~;v)2&8@rKLT1T#d%of}ms?hL_~#_8f}& zxK2a`f}ZdvD2VsVBm7QlC!tMXY;;RKK;1gP)k5v+>gs;6o^WCF@Zm%4W^ZB6`n5;2 z+v)t)kGiAjqH{BRXPrKA7&cOJCJe-OC=m6qt=*$!4K#RN8Mh~_`k*nb0|n=6TwH*} z^)ewy^X>PP7vAA8sbkd!mM1+?G~{PXVA?Jta8G0bLl7_SV5l7zT&I?il_TA(8Q%UB zhVBs1D2i;{+<)-BTb>84`TR@08rzw-c0sOV%7V*>J#;-Ur}P7qQv`n0)lq`>Vb|Fs zHk_;Z{-La;B@2s6YA7g#be~#oYWVN|fJkuW_Ph1I1`S}>a2nL{_G0tt*Ejd3%z0z5 z&8~K#doCR|Yn>?S!5hzI6oo^=VhoxGqx00Oc2(aa5_W?I(#;G%ao>}lbcej$pw7or z_eVtFvxC6^??umezQ(<+VgSkd>TJJfr?e$%Yb(dR3sMwr(?4qY|)9ME#Z~v>p*Xj;J6VnXh-Q z9voESpv3lgPPKf@P4_Gy$Hu4{9aY~2 zIbx~-fxuddz|OjJw4UV|{LO;xFZ9YCS>EFS8bj0lzfGgx8^Kmf4<8c%miD$El7lr= zgfT(QY?`ck3ci>22WWWcu>%~0BJRs{P~Y3{^P~k-{eT5#fBI10&;Rp=Ca<&m3ca_# zzn!CY-v91;Yo|LNs~zC!!|c~l6G_6KXx+8zxJ*e3lvkBR(4B=%AW5%zo_kIJsp!*I%KMVxcOIq)PSNG)= z=yUGX7dcZI7YzUT=h^J+EV~~VW^)jkxqlZT2~h&fc2hKpOFJ#T6wZE?55x{-Mr*F9 zpohb#2D8x?ue%CNtGh)vhdQzPGQCDm!do}zhg!A;DzkeFBRh`$BX^(gbSkg%uCykI zx&C;bw@Svmi~Wj|GZT>@TR`M;sv^TP~5R=e6`2Pi?U+}k&}4t=_eR3-Ucm;OG}%tDRV`V7&`!u@#9K_wfGXpv>kbK%btypt9N^=;HG8 zSEb&}UX#JY>&1zoZJxhVl-*#=qhvpM=a9OE%>C3p_7I!gv+tbCpwaTBaCAvkLp!LM zv+n@g(82EtwN7=X4?knkwq9Qx-^-h6|F`FL$cll3gAPKoKUuPGnAQZ>ea;p_z)sC= zZo}ojjCO4Qe}0p@)}4udek6Y0exF8R3*{KTl}FA+e}AOAOzvBEj)!QiY#FAOh;^Gc zx8#$&$5B^f>T>U!YjaN3-;>C@p_uT~{&o0F$voZ1dgMUQS$sAY{ zDq7jB!kzUF? z?-#P0$?=SCzTJ66OGQjn$q;S*89B@Ko5>>O3DJG3ug&cIIi-CekEea7ucXDD&|`jk zl$pzJJ~|pzyHv8C;yC4W+9CJ$b&l{(5JOE0V((X59aw^>l^;d(AGrE008+Tip zH);v*xveKrEuINbc*Xx^ye27QU-4hkxoIj9n%WruJZ#SHe8=eSi`!ohtfz8O$Wy+w zw6v&QBT-tyJMqV#CE-)oKWF_gc+MLe>K_6GBUfW4baHU`PJySw$@meohx^JmkBn9I23RwBChZ9wn|Fz8hOwZdB+@2E}CI z@fAuS+3Y%BigfPEO>$YqPHg(iR1NitgrG|=jVzCO^SUG(nD4J=_arNYyHg0O7@veI z&iUw?;l4lgMm919`e@wHS|FsQN;V3Z2zg!o%`QiFw_#bpd_EBzEIguI)ec)Lf7Z~5 zIUXk$9w~H9f!H-RdSs&H?Fed#J17Dx7n>P9pr;H4r3~y{W!2P>8yIXTyAQpbT57)9 zTa_?!>SNY!bmsxBMsPvFiz5B{NB4E#y=&&g!0|tKI^K+&m^Ps6K2q*;&J~|64N*E{ zPBJb`5h-^`+Y@lTq=>GQo-Q$8xqXc6#4d< ztIH+ewV9(W`O&GS2@J97R;SDW*p(KhzU9H6rBQr##fm$X<526w^V>#G>BkIh*z;p# za!~oFE#%vc?_UTZ`_I~ol5uYp|BS@6xZVnBH#wN-$ipL+vj{ONR{RBChsZ=s-6~vDb3vtCT84EI99ge&2-w zEz(DBmCy9pp9_MyZ+~GEhLy$4^ju8{JAzDv#`z!4UD5**CHBD|;%U@{SWlyHv^YDF4A$Q@K4cuYy*L7P!>K=#GPH{0 z0VYabD^IqWOp934%j#@jU^(oxX(@&_H0W@uB-{FE7`(sAN_rhMEoKajq6{W4Khs-> z&6Ii+XZt1R2=9IFCa}88DzB3*+s(6)`?+SatY(9;Y<~=qV+N-7YHmvkO8lCNxwb*a zZx|ovFT8rhV`iIPMFXn=@EAwEkT0!6feIxM6K|^G7AFKApZD532)<=`h zl1qXC$9MJ=2!eVB3BUuJ`|Iy>$~DSeb2R!m?-mMAP>G#e!(0_NcxqRvCviit3=7gY9dVL*_r!w)MNT3@!a9k%Km z1U+@Ec3pWF#Iw8H`Y2{uOiFdAwjXm6466d<752tUNVX;vPhe&yy#on zH7~<)&HAIl3=lpc7+Nz)x;$f(!B06*6l|r?8yWB+uKebdbTrO9A&au=ugbV&=+?xGj)XMf#eX6}@=XIb1Bn_sl0lNZ{T- zn@W0|?2hoPOWc;|FIZDNz55jDU7hEXRfQ+kEhn&P(Zbl-9R+1q%n-8R)Z-`L{IML? z&BTINUfw4m_ODAmc=~J5bj{Q@2{P%YQyoc%LW=ffH_NoCkH6i&P=)t@C+AqEgQM~^ zk!94y&pW2ighv=aS?F;9xDxg8CHH_D%GG681a}GdL%2Rrc08hV_@k53ENP&-yXhI< zhc6^yy{1zwqJm_u@ZHkK9+(UUw$*k{7N-C z*B8i+ccDs^)RFdfMv8f(QHAr*ey*+wW)=38fr%MAN$yHLnPm)sR0TU*7>-KQ8jF8y zzOJ&+cgITob$t&tAw1)*EN4E|MfY!PY1Al~DfIBa0G@+Y2D_Y1vxulD9oQeQbm zo{Hfp5>N!-z@PuR3|yF>*>$iWtd%2*WL^RT=cAOq?p)Dq@$0hQJ-*s3&c`sB4)Zl7 zUlO|ZH)5rCSG)!PyB9zNW}EtYEjZ?Tfqf6gN0ZlZXmN#zs}jKDRsl_q;p{cedO3qG z>Bio+zZy3~sd(8fW1>|+_tYs)|4xMkK%Xaz6M?r=^0l#$AoKyt93~pG4N`&;`SqK) zUwCCBQEUn!J-bruKjce`7a&ES+`iZo1C}C}%yiG*Lt(I=k8rTr(JWiBk5O&EhI6q%qz+=nS z|2vMP&kFmH0IPNqdMECVW$q=UXEcL?cWq_<=D9{t-|7LI(1Z!s=I|S*lFtMxeL+z8Q?WE-~W0sGoC0?78199%pFa6%CSRJVuLGDr28%qJP>dZS^oV zEDk4P6sc*>>X@5c2OK<=2uLffr19cl9m0q7V@L*cL1OXX-422cW;k{F4V7`U2S-z~ zbwcu|3zmKep^0Si?BzubtFGb+BF(JCfo#S(B9-I~YJ}x71a81} zW77G8^x~9)S~6dnBcgo!5`H|Gd_B0Sv72e^<~1&Bd>Ur<$f$QIl?cM8 zDd1X{4xJzm6IN*CenoW9;4wUTo1-FEH|8?k*|#XnFiB-Hpe2%mPu`XrJUJReJMdaiNKnRC^%S$u#%k)`aUzZ*;gxQE zNb@#?%~Aw@n~fHOladtyxTTrzeS5BUBZ|T9oSTfvpcg~Cc0&@z$}herobOBe41?bJ zPRf^13c}@wUB;DZyl50d3m;kCV4VBgV)@FroLX%mvbw>v;)<^xJ#qer<(Z@-7lZKt zY|XtQxABXLr+Y38il$hD z2ZYWn-_Xh;gIHbLCc%Xu3S=>V5m?u>dUgjqoOVp!qy@DV>MY{7-q%MrMuM5_tijPXXAO3@nCaCf6l zz>WrObCxJ{ZIi8bH(ZkdZ;Fh&r=*PswsC5@h$BR!a7gVhPNbX#FF;OA~Hkd-^_7VI(?OZ4>))x zT6?BCJ;L)en)A3l>D`b}- zlfG=a9#Fyhji2^pMaPkL1swn^IDfkPZ3JF7qJC5!xa~g&Oq;D>1sO(uSBYii&h_vJ zkR~kOivih}OEAQFeU(d}!P`kLfnZca@)|%uec$y{t3dofUpT}7d0{;m1v;{6hHqFi zbCmIu7iJ`q^)YTN&XAozGWd(H%p#1mwTh3sAJIs*rpD|_{g5YOoyiwJe@HPPF3bZ_ z10l~aY1s$U8)WWNIApt;q~S{&LNW%IYPTt zd-d~^6CD(ypo)5mg$$C?iP|SPA(QP0admJLz8m;Wa_yn9`PS+{v+nM9t_40S-$?i< zF*h{9@ta?4acBp|Fd7%bc|e3IX%Z^~;}QM!!oq~GeEk`XWTDX85FX{T!F8u30yTsI zKuMtX8qmr2BL9(&egrEaS*6sR2RlaM9`))QRewa9jz4sQrnZFK!t5LQLHoXI3Z9)O zu{Kit4Y}#W6=hlDFcF1L=gVlPl7xT!hcak&F3tRiq{xhK>W*u^bL%*Qk_hH10q!z+PN+T)j`@~I~R9I_B%3a}3PnNyO_WDTIq>{667@R~OJT&yuQ@|mfOn}VoiHtT0j zlV8F+0Z%)h#D`OcXxLKNcB5mGAUtqGW?gjb{FM*ku@-JX2P|ZCJH7mv{Hn)E(tz{k zbMKDhaR}idZ5gz6WhG8);AW<_e3cXG1|;?FrJ#QpsjRxQC%M78m*;Px1WS|20@lOn zl<7}QNcl=)*C4vol3h}Gj#On4nF&xhidr()c{C-od{aAL+7?d#EC?1B*zH#nqSJL- z9>%rg`m}5R&->3$*gnXHD3FO%Cg2e(!b5#*%D>}*zmmkTw02=Fs+mVR(dKAsj{UKDV(XyaFkHji z=c84NMAs^+lf?U-L(~u^`CORPyv-VR|0sa$!FLxMRhQ&|omC~2JLdq?G6Hz2R6hOh zM5|_kSpRcxp$rPytehOBjEv;WQLOydw3+`@$`;THnC4_(V&iPHx6P%aKMVFV;7q8O zY}+W~4=G-ed6_hg3IF&qhp%&E(ClM$XC0pvg=&`cj49q-*diNVSNXpBcf7Y@aoSRq zL9)dph!=jN&rR1(vAQNIR}w>g#`uS}00qz(@4=^rC{Vpe3&_1}_rDe8Wr<&CDkxg$GRkA({sx-+sH(ny{p^Hr+5+e<9_~2YsNzN1>M5 zVXuZESXFLgi($#9CXJk5xpEE2S4awF&}t>sMwt$v4zl|pnh28Zck!lPGl7ZQXY@;_ z2hNZpyvN9PmHk}gRjL$^Vl4WkIZ64ej1aEO4`1u7Y#A+08Ruzgi!ANmv7>E?VIWCa z4a6K}B2&eze3R+R0!4yhR;r|{;K`>|T%PUjWlyUXKC1|3eRqG=OJLN#(VKdt%+z2J z=2lwm8maVpEXsPWgq~1dL^1483d++073!y$s;1C-ck-Uj$Lcj9xn0Nckvw`BYcIFY zW}@fJ5z#Z#%=|Y$`a}s~ek^X^`XJOo_j7xqwT_#if7@YZ((g!`Z%#zy>PaL^DJZ~d zzjyH})%95aUPLT1Zgwa%ZhobURbrJrfCKLd)E zEA{1{$zP`$-gcPnS=o@lQ4U!J(-l698w8mg5i^~b@L3b<9bO7MCOr%55$M<)lw?)U z7HwcuLO$x9O7C-1^Y&pzfXTE+yscq@Q@n3rxP!Jkf|A$H!4_s56Ft7Vl2@iM)$epZlInJhS!1GH*L zSE{hvAaI8Y{ww4hf8X|DIVSznV=2qGCSibW)%fo?JOQ3N)w?+tm7-s)MR6FV9=ZB` z1u#I?ic80yZa`T9&)wt0RBBloN#X>t$>uVCmu^s&=L(V|;yDE{4BNn93#HKuv%?pW=IAQO0U2=>k168YD^2Pn)!JsjDiHOaS?HYzqyE9w15tQG{Bwn!(YuGV$` zicc$l#bE{u!&fTak|vwI!>k6Bpl)XpvxG8W&*6PzDC}{|)EQZ&N{=hhe~X)l_CkeM5rlsf7N=}H ztj?13d^oI0o(1MfEW0A!+L*R0+W9-lw%yV`iV|j?r7G%nC6mhhMSj%_>>;+l&L=~( z)}#i=D6PFo71E*9aPc!o`<>qNEZuFtB=L0i3!aYCXeUn%D%vMuac5)M#Bq~E0}zRR zuD8tL%ww^aeFBkGA5lxDc;5H63fc<4hA0Shqpqf#m4A5X?@$)*L9cw2UIJL;pg!|p zm@q7j1zSNJCK2dLtjt`>^Y6!JODdt=Bs0#V*WOb=t4Pwc;-7U%KK`VwI4> z{z4{fW=a3#4?T!ZbdJYlSXpK;*iwX?nMksPOllQaU>6EnI-$ys(vdOi>KSsQ7Of|n)Ncd(dZ5T~r19r*apxGHR@e9@bp!E*|~tNCKO*-ojG^h8HI$F`7H zK-|eq!LL;8$?9-=mxbcjXoDIvWYJy{d*?O@xf=hBzE+{y_}#(;$tHtEdd74(Iy;G~ zC1+tBwoKYwf-M>M5~(fCpG1oksVk4N+pzHR&hko$%--&(e-iifg>yfH>*4Z=0H;%0 z@gDba=s%fgLW!wD9NqBgN*5+L*F+%0N3`~};dh%H8>=(nt!PK1Bq;7eZQCTDj70O1 z*~YDL*!e6e~u$l0elpum{Eio;-~ zx?)~&Hhq?JIfR8db2dXUl6Ye_rNw~+tPEc3^a@DFQm{JbB_Y_!A0~)=iS@~wnYO^) zb+RP{l5?1cGR*PJo&=REo~8^N#uuwqD@@POsnvO@5G=tZ!Maa{+wg9U$nf?by_T&A z{E581ksN<_N6;YTo85w?iYUZ9+x_Cfi7Bn|Px^G^x|qz{&e<@jS1e!5IbW9)s~Cir z$v3fU_R1sXtW<#xC_N$!~KZ|bkTu?+N z`@Ho-)-Q|9HD+k=`Z%2$j4SGt(6>P=NLgQx+gu0z+PU@Ve5q}g$9I!~H0tKeE}cZ& zUYTEmbd)JHBtpP}f`t{n?|jUJYixa@>D_d>`45dlXV5@-K@CSk&qlUj+YA;Ll$+Qs z(4&}&%g()CG}MQ&4@JF3`jd)&^Kq>W4%Jh6^5hZfj&*>GuV_C zMUA9gZ~DY?%n&@@GkX0y{OyJ3vZW5Tb}A)kUARAoMWro2$`NLzDHp2nVPlh|CI4Oz z0|lOAewVjnXvHeoc3D0ow0xa4iKJ*|_^HmTId&W~OI)cxv{7WyYYoo(=FHhJ{H0x+L7(;1b@iDx+e*}sJ0TJ-6xAr2vL+Of;&%9dg*q4Jy8{a zRdmUj{yBg>)n^08=WcTI1Ic!MtTwqB=sPJl)RzJ<6c*`>kJJerSjnwmtn6n@nogNG z-xYJpDI#Ofp5I^h%haPj^wszC1|Cn>8G@aYNuhpRxoMQxvfACGH zKv1Fr0&N>}EZSa%-`IN;?tHzp$AyJG#O&L4D?)TIv3~!Pb`$+vQduAjl~$*@-WfVVVovL{zAli$zL2z%t`2oaQSqL!PwcM; zm8OU^TnN`0l+|Tzlu%suJ87lPZK>$UI>yszJTsI*GI zCW&6AN5hX^!kgPurq=LM%r>1MTd#5fBskh4j>Ku1noJd&obqh5AnSpj4T;0P>?}@U zk!Jbb7jmpbbnk1megCpK0yYt3t;{;ld*x3%HXUD$6)pF@QKe1`)`;m)tkq`amP`b) z+lL>c)!<4N?sVz{Ej@OQG>P0cWy&HjVvSb8gsG9}Fw!dw8#kIKWTrjMqv|!1ZNAMo ziVzFuM|hWmvRtP$Z~xXdR>8Zvs!A*c~t29`tyUho?uV z9w4IV`FCO7*HJXIY`;}t_dpdLTW*NDhBTv#%))L3=OH*+PzV{~teBWGVn*&xEh4Gt zxvJPv0b;#*d#W(n)7`fNw`1GfHiWZJd6;Znsdjgh+ZRMJCH0Fir)6MuDep0j2}2{= zT`J*-9P6CQPjm=PL|I(8_LX#e(K-MBczzOW6rG!n3NDaY@MG*xPB0Yys&0XkG#h%h zZgIzOCAz_hSrNIu-ah=~5lIWz)Cy_8GYZM3d2Ck=So9s6LRnNaNs5~ALtkKgP*G_J zwT>0?f4fbU)XReL;^3I4v|>2ENFS>LEPq$||GN60M%k<^2BJ0@BV! zCjIk@R5sC=nObx6vK_iD`Ao4VeA2B2d-tu3=St&Rl(abxI8~VMauPj)9qknOwLGu*7?eJN7wD*EBYe(ML|FF5kSpm6J#S^q9DSQnFnyrrbzPnk zp(o8*+;oe1z4Ht?r&(4(36pJPrIgKHM>3JYj9xEkwk!{=Aa{+kQj#c;A;NIWT?9F)(v+w+dx`X zhMYznRS14=LvNj*>=9C>O^}YDz(Z48R|&r%S)ip(ao5!J9o;t(@6EWYqKb zXGaxsnph5M_q5W3117Hixi&MchYjKc{HE4<+?91lnG9~#vqwL;f9L9ZhZ1~;IK zsLXkS)2%E}mi}QZQD%L4BZn1Nmz7SW^`eiYR*ftYJ;heD?Y0t|Q^x-6Iq4eO7Wiu!uy7a3E^%%)0N zOH5RyvG|oKSW){VAh3hb0JV>y8R$c4IZ!eum2J}!V_wL!$yu{}TwbggF7*%Z!$~An z{Q!K(bWeF)vD)O7<(IVC8f>FpP4k_Mdfe$`^cNdEXgRH;r7faZntqXDY9|vAl_Qj22S3`XJ_}iuh!+ zc%p^mywf{~d?|!SjR-u1ltEJK!J0z~N2eu;4$}E$v4N-(%}G5+knQPh(sG|UNl_yD zvy!r!IPeTrA=;2v05*eU5;Z-d?WR3R^>mu=HG$VQH9MDNr zw)I6Q{4w7dWBTuNCN>dF(IPYCL8fkRtFL8s-#1}+cjwFmiJX;N1Rj9sH6_vDkC*&j zeyKxlYL)@+6?hep6%lr>YoN!pg~7I=O`{oj2NVkBt3EPi>vsr-Tc9}vI{KUm1b+|8 zglLr67FtV7jZx@KcdWo4XI&#>7y#xtKd_9e4VW{Wv-BSJP-5YJB9hg(Xa0%7xW+VX zp*FrbEL`&PA5cL48Sn>4sU>x^zv0Z#RPWA1s{Aqr)>qD2Cj*M(_bKEF6Ft1=2&>U! z!ZlF9Z*a};e6(k3i9I|QQG%{43)1gh@(P1)ns!5O)(@!?L$5d@19Z33e+}M20S% zkrIZ75-!2dy9aT;RNI7i{+^(@Cz!*)2MGfrkWi>-l}kY!l4Gvc)rHL&IOS-FFnBXN z#BCGCfyMssn0^5p%0s~>!Lfv`N^Xgm3y5O;hilpoppB|KI<#p@F zo{W(ysDJr;;kL_$3TXE(U<#GBV^03$schBD@sM`5@R?w~&T0xVWD@$?3b;C*8F% zy(@UjMVVem0GYzOh*QVl77>7TD*PYaeP!k- zDaRj7)13_gc2NAnvt<(@N_k;Y{%uE_PEx+jp_ytGxZ3381K6G(i3e&n6mIy=(hx=a z3~&ynk|abtz{2A=MKq)G050yxD`>yD43mjGa5KsAbh`Y8BN^6ughxP`l4?$3YwYj=@@?v=TO3d4tS;!CsaqjBX?A4AsFY&we|x&P=gb00AAFBp8q z!DX5LyI3s2-+R>3X~zABA?3PD+^?HP(h3R&m=(5x>>XZaG846#W7MzWGhu5lQFYZ| zUF_zbVT=(wHF2Q|8j1^o&zzol&A-8Or`df?1zj(iH6#Qm^-1aPr5m+=2xd>wi7IHMd`^Nfx+;YQ=Ln0+Ea52CfZd&WN?C z`{}WGpTem;4Riyjz;oKbkyjGg)1BUIT1rtKp*6F9`8 zvcamh23EKWcB2A2Pnk2pWCO!YNM>L-4E`3TGNDq9bm8z0yBPO#snO+lRzh? z-E^76{;k#*?7o(N?18NZGD zM=P7JBS-3k8-W(Au;PUJ0X6f`PZRAVl6OAT_S6a0Df+Qn9|W79(L5I?ln9xRp=0Qq zF5l@z5sp-oJ+~sAy=1upHz^2z0l5DVdBRx8cSWnL)~*Tt#)q=Zu?7TB&3*z5B2Dx3Y`{5tmt`h! zxmUElmnBNHL!o|{Q!te&%9NYm#!3KrLoPh}s+eu#f zk61BbGY{W+hPo=5G`Ku6BVmPFmiI#j|IK08P$W5bQ@oiALp#WBcd;fQ0R|~;MBY?j z(nsXC%dq2_pf%7_*o+!DPhzj1g`q8(2Lbq)jpUmO8;>BZ^YXHZF9NOK5yS}PcvwSz z2@%s|@?m_owely@Qmp}&zKkhew2TEDl^@p`#D`FQciwQScW)90guEge{UD`1jER~I zODCCxxzXIE#R`Ddv`ox4M3(w$CnvgWZTd?nAyyn20k<&){HXcYWmBEW#+%zpTZ9l0 zH6Bjqcf(|j7$mJE&m388x3+MCB!d0IoxQV$?}TpuuaChBz3@@G9-qb)O$nDg>3d3g za8#y~X7~Tl_LgB)ZEO7YTy%GLNl2$OBB3a0&>$TW(k(0rDFNx0QbM{@Kw285LrOrT zB_;o3dd|7t>y5APhkb4KX0hg&BcJia{oB*_wNW+19{dx3&p#$_TP0RgQ9_TmSty6d zsp>I1X=LegRy>)?$%#Pzi#Vt4(qbG0F6mM53ZCsWbp=j4cDKtZX2lyK%fA6{d&?y8 z3#jhx{CsS!PM4>x+aSC*ZrbkwLaJE|nvx#KT{_{X1lgcAafCl}FhIM}-puCo6j$9d zVMZ{xB!h_V{_&H1m+!^9QEt}PFhLlx@nAAQ(= z6Qfz6K}nT>Fs#2aeGg)B$IGMnyI!`s<2?3b_!)fJFUUt`-1f<=;ZI5~R1^AGiD&A0 zED8kG4^emCWac!KLHozyArB6`_c!)3KMsjWlf2-e^v9<-M7O46~ER6j>+3u7$#tn_BV+o6^x4Spc8MY^Wg1?pz=0zSU* zZWkRzw>NyUjk>cKkn&o3pkf_LieNqyeh4*UyjCe*k_u$nUH=KHuFUy_&5>M$zEo=+ z(;ipr>pU{=94%v_q}zA&lWke`8`q@YC{}^@WD1QiB58 zoKW(UtR;dQ`}Q-QUPHP4H|9H!uNd4sdgIF=_eMjn-C5hSp}8y$a*z5!h9Cp07_HV3 zuV0Jo%gcA_%s27LRpz4j9(2mpGisMXR)zwNSZ!Au>QCPes`&8oNtKd4r`HK<$X9uS z4w+0&`^{fZ2%hza(S3Y^C$7+byXx^(%G0>l8g{ch%$`y?NXmh>0E+iMo`hvtvX5tW z16nW0-4x9q9*l4fseh?ca7Yomw@V_Nk2la&oNmWMXsCDtXWn3>_C#}1l>dl&q`@28 zbPps2*`)10{rTq7@G{JxrPkwk-%o5TmU|G(%=6V57F|0h0Tl=bbm8l|LNgu2GsF0u zRgafkQj0log;5^g8D8};_0z(T6?36`CZUZe$K24k-SI~a6Do|~?xn~#AChq$C5p}! z;a1vI592ceR^#s}+!s2c!Gc_0(Y9O~1w0=UU!PseTJkY;ZaJ^ocjuz^$hh_{$31z- z!>mWGql4Ew+(a<%H(dB($0~;Vu6kdLD@UKg!b#m%r)HyuqU2S?`Sj8=W0^uliiQUwO-vPi>wwD1(mVD0?4a6e<-L;Q^(Hu@vZhS$CGgEg1lt zKW7&y4Mm0wi)}M^L1GTb2u#k*_h=3`Cj?$yDvRqHQ#t<|V{-VP_h~oegP9g;56Vd? z41lZ3^;X#uR3!#-P#%dF_i5wg=MYLDG1DgIlmtB^L)`V+w$WZOkmX}~a~>A2;wF;r zbdfxP*j*t$iH}`mNEO4n<;J3K_v%^DH|6I8#?a2B9U~)IR@=(N(X%FNR*u;CD4@&j z4(tY}hl0bJKfCC$bC7q6?0Pu+uvWtNRCj8POqGe!cjfAJDqTyo60a)`=Ze*pM3OWm z&90Q5_n4ojNZViiotHgUOLQ514%&P9*Vh{H>&j61X9z}1k8yV&on>_41FHx^If;Nv z=qQW`j!Om&TG#GeE_(7@P;K2@*(C|k&z(nvi)f>Kl%Y`owI`eMS>^~V1u~Dc{37~x zH`iG7a!U^98KK85Noi(CpKl;(jussD$t4>jo4Y||WW9VPfWAi23wRbPZ?pyB8CI4+AdG4qlro;q_l{#u{zrTC` zJ~}WE1qQ-^dhZA4sh~`_bInEK8m_2B9k&Y{oq6(dcVZ4ibmrm(3Fas~Kb(J-zU(;Z z7PW))(U~p=m+LrZSp?U1u~fw;>y+c(AiLuPk=H-V(&+RUvKO~{+`4Z+x=ZQdliu6>#K#wdRc){hgcqSv^IqK3mrW&95&QdOF=h&fC;2ij45lJ;YXT&#P*JJ8Vt?%DLXya>43;Rvs0>i6LY-VVurjCZT zW#dSO`T`UPe|e^36r%YclZ`e?nK&x#;!|C1?W+TI5;UaIvCs8lVGB|DYimykL`cAu z!)(-wAbS>ZIX$heo|qCa`OEO~>boR#cXUJBG3KPyw)3FW_V=vC+m;&Wdhd?LZY@2i zYZTniAYwL&1tp2nnhhO+fW=!$A*ie4PQnHV6p||IGrG4CZ5& zW>}qnu9PW1m?mgBP3}Vd^D^2PArAA#Q@d{1h=rX2?7w5AQxzS$@^<{KgG7(QL9xL@ zJ`Pao`Wnq@lbbNBl^^{Hm8S$*0FFCT@#L`S_4f0o-E-{*Cj|3%`;s@e6XiK2l z`whSe;e$fjdjR*e#|RM-gjU+m2+ccn;(xBgM?>lB>#MkA4$Cv;AVs3va(-2nAKg&a zFjT0mB&jc!iWW0J#dY@O)E=>$S|(s4No&MXtm2&$F>p3YM8$&m7yQ`4by!;D}QTlVYt}bMkHd-|FrpPS|RY=6U$)V`x^yc zc05cijA>-0#^1;s(^I8VuDBlobECDC0&3-*SU%xA5eX*NuyJb4STll(t99b273)~) z387)C#7IKMvvWa-JO9dd2e6l+A3j^9UZCFnU?q~D1(em+*47r5mMZBaH$Iz$Yv@6r zjvUNbzxT!fNGS!!-Q;>Cx3gdMmQ$$Z2A;Awjlv^qLMU70{>Ub+L-1=)+K;BwWppR6 zyuSy@wx*qa&OCRikmLV+ zMDns4;ECuY-TMAC0yxDIsLI`Xa+b{0iUsN(O?UI`oUY=cqEDxdTQ?E>i+_Xo(F7wD z_Uy*dO}wM4IS;4recthHrh)o8+mNj>yIeijSzkNa;3Ye{9E^=P?E(BW`f|V~D$~fi z6~-S3$WBK+9^ed;EpwEA;g`oEB`_7;B*mg#n?t#kiFj8ue7%--> zE!5GF9H7LnQsbx*zxJvwvR8l^jSk1lJF0Er*pv^lBv_*AM4%xLIVf)#Ej)t9f(2ap zav8(ANvB^1zbwzt#z=nj9N%!!(eak9IS&dk{avVq=_vn|+OM$RqR^x2OcVjKeY1>Z z6#sOW3()Q`RwVFF9HC;SH4{Cu0kNAD&4kwnO6R}l2v5@ZmYBZ+pA=PMfs}7?-PW2> zgis*Xx%9VHp6lJGmj~xK8KP5%vPQxlS4MFg_6$0KTMY)RmL`Y<#$*QXYB<6X-MDzW0 zj<>%?0{PcADa#V~E@&v$iVf_&zqr+3j|o?6 z({UFPz-@U(fGH5w#t_vVLW3M$83N(Mvzd_917~;KswY=xlQmU{-=gp?`1C(to2b>7 zTKbE(AB#oimy$2r8&JH)Ps}>Ej0!o*S|7E$C4k0K#c5)tFhvU*EjSzh1fK@Bt~8<3 zeFrHqoq7mT=Nh8R02{PCkwO~66MFvYctW=RGr&GvRC1t=$!$A~FSG%Z`lI-vO33DP zO*ysHmSja^6FKDc&~`HWU*R~_EM5bYap&dxHS6R7bWoxpC%kfShCGgBiVg{M{sJz7#AvcZ8?!c^)s1=+?I_*obBCh!yl71(5t6P|^}MK~8(nCN z$N3)nyFkPcbidgW_4v$2wA5}p3?DG*c}sm+3jj7O{{)WJZFQxWb>iSeD=LlRh1EM0xGIcb{C}f=~=>L zow6HApW&N&&GV^J=Z88e$=ggqHIs*~fSPkSAAdFiCOeq3U(7Vomw;`v4;Xz^fNXME z+tefpct|n;LP*v1o;WzIw%M)H6I@(|Gc~c=`Jgps6G+_7jjvU$d#q=MI_?4b=SeqYMLTe=F!MZvK=2GI{lk zR74HlH`AniA&B))hUm2+A0~-)DFv`p0Hp{s&kij%q3!|Ui@76~ zk@$x~L`U9Ue2eFw*jrJtb{65y!Y`GY3~}}=FJ=Xx-tH+YzghhZU>SaUzKflJX168( zi&~z|ddl?G>(JKPbVA7KHMRvW0-A5}79N(n@msq94PyWx z1QDbROIT}Vu{x!fIqnd7YjVb86v*rJ(5K~cSstqCgcs()nw2J+iIEY#pKR8jSLIWY zfLMgsCYC>n$_$-DfEr~~!`b<_qqMxS;EgB!~*fR?BA2#Z$}MN zWQjfnFsTZ_zcd7FnxIkPzz4(X3Q#ofgZdW#uUxkt5cx0kR7ff`=3??_4<$hbfx?M@Zl{6&+SrCuU(~m0Mf;tSLK0 zyb{D7O;>LIN5DboakhldBo_Zl+69n8o>bi_F&$i9c9CT?d`n6M4M+{Jk(C4KT5`XKQgm zP4D(}W#8~HLkS>*6ckGlWc8w#_YTILZEMm(*4K>xcE^`^=WozvW@b{ru$g!UF^M-5 zcrK+!PMO)u(-@!s(A#e?f{=u#U_ax({I^(pxr3^(f|*)Nx>m5~k2hjqxCP;f=xX=z z#hEYTt{n!X%FB_Q23|+J6}|Ser^B!+SU*bx`6^dLzRr?(6B5&)P>-cI#GGD{tA2Ln8)iDnU+M1(z#= z`zr7jHR;rX@s8;n$0!HOgW(B{3?HpdU|b7yeOPXAFr{&uIN||>-28bBpMACvc$=`Q zLtWt}RuYVjGcB!z@Tx%Q4Y}3IuOfBn=4CAC(+#}gd_#SRAXQAz$WvqZexIy;?XhzE zyx*4x&|Od6j(3J~d|7u7x!{R=4er8d#Fmeq0;5ms!w-+nzpR?y=R@QrWDNs;VEwfA z0do;v6ZUXm0eWXCv6+$#9H+izX-iLaL4(l7xHWN7n4qU(ZrM4XX=P(*M zBp|@XRC0C}ENfh+*KFYlL=;tzL}geg&htTh=bDCw zZb>;t^aIx1PLjK(;SFf+Yt`B3Pz=rSza9}cy zUczWH5iv2d&Nuk`K8a%tFc4u!GhDW>wMon^gze{zI z0M6Mzo(b4N@%9Y#^&?;$9fOAb`rO&Us+NI4=ZlCPQhWRH84Q_O~_eci_TQ-3M`K^#)^K z02E5W)|Sgp_W3u!Wxd+JxT@<>&oJ&F#^!%15fSeo@8_$GZhU^k2`z0|CwH#pp3A{T zhRHh$a9wWkh7o1P0axxn*Vkb<1N1H{D+U}@P&npFN*O>dAjMt-;8}dhd9Qw^p--Xi;l%gIt|6SbMc@l7XTs@G11K${#)TL~hp=3TI zoPTVtRmXsn8Efamp*yG$fEzw7AcSt;z6~GWm%r}IbtG^-(KI%WhFNyj04%hLlP*5Q zW8_Noe7uH9?t5fapsJtN_rN%5x!tr^6g_x32mJq5+m~AdRoiYvSBYS#0YWgZv60%_ z+gm+UfohaIV(Ic$@z^+VM|FKhNh* zD!$I$bXnSU0UA5GMcE988Ts<7nJsA>Fe)9v4)U&i{tcg{30$5LLmn)W@xGU)01oxwJm$;X=E6i`O%X7Q!+~ldII<`5huWIm+kASlJ=5` zz|~!t-gFbcdLMH#g$jXqQblL%zcsdqL^EMzTRJu|)02fv>?XKE^EPk%YHrw|rCi}q zaXM`Hn~v}wKqp$mV;^2)A3haz%dW~nrFNtoW)$w*bo#f{H148Q#0WVwRPpadA=)BInb1c7!x?~)b$R}H7ZonuED zlRHx3nuA{vee*@- zX;eRcsMs!9`1utg#=B9_<0T=!-}^-ACg`lXP83BsXF?+oKY<>$cyDa2%)pn$os5v z-?y9;!RNqL42+vbgP$38l^+8*BTXKyae>|;ZMjjRC&;qDFP}tt=ecm#%@N4jkcyn(1|K1J$Z0zRY(gh_KeYu987es7p z5*0iHMkiVJ*|$ZxmTlGlng`d|*zNv(r^Ae|LAMi3)~_<6_MzJozB$OFtYhZkYWB~N zs_NT(_ivn^!{Z!ByZ@>7C1GnbX66EMDC_T69t?IpQw+~z6gcJtMfi~y<*7g29{lT# zm&{;1{$3|5=d>`9iVhj+)}3U3Up}HCR_XaBig^$vcTI}9HK&!a)%%ThSYAG9WOMZ2 za~8pKvPkxQDt;BOLeQfrx-QV!p{!L|t7@5N^>VM?VlFNOffY+m8fvy?L}@;|f3H9R zHpIgG@xj{+&VT)wBGXE>K#Yz1d^&p9>HKZTjKN$)x!i$)8G3!=Qbmp$rzE;oBP@T9iOnhiUjWa_~m>Yl^p>+rOC z+75i}^fd4y|9}oAh(vcILfU0##;GwuCh(@^y5c*t-b%aze#7;q#wa)oX2{XKmKUu} zyRz;7yDn?Mz|-unRd(yTeSA!KmR{f5e&Zm$Vk2#(SAB}+=x0Sb8)!O6BZD=v=)v%( zy5+xX#}+JSk>!&}}CTGNNmz5PC7f-{w^>3ok7X2*JUTF*Ggy5*`ydgqydSgQMTF=$^ z&ZudpeYWe$D{*!AUfx;A00kH3x077$UXsB;sz*Uki zQKycsdx6`7i7M&c>PUOrb7`8zn{S*JQrm^1HRuQZlX+*oCp6pQlXR!e-!374K~HIr zrSQ>;5NxcuI-A8?r<=iO+34qgXn(bbi(R%Q*O~jelR_p=(S&jhUU=8;=K^Gz>sq&) zW9paqX`@ZL>&%_=UJUaLfBBn-xjeA?W+|>ZV%a%5xtXy!BlRfLA#-G0}b1@ahVmO-1O(tX5Fy;#?#vK0O4L z_l5+dY6dpO4+ywEC4Tes+#z7zN%5YVE2S92i}nG7XtVFj%EzAZ1BcE#jew@>Y*Q{R z_HTaVmX|LU1l1HG3}K-4PzKs^xPY#ml`ybm%QZT$Hndt5#Z8IOs;4YQMe{Yo^|; z*Y<6g64!VZe}$Jk5LI}=tU~H_n9jMkz~to{ugB1uBIydBEc+D@WMe~k&sPrwpq#ot znB_=B#2Eq4HZUINbO$t#eT6$0cSPCUG#SC1Q< z7Y3K%-Ev#@N1#qUL-#o(!S#Q-EkF*|N&xrn%hQ0H zJ0Nfs3b4_pK%qnBhMfwJg$nP}jiartypmEkvWyu%9Kxm*tMEIRy(--m;wr{*VXYsDC2HpYzJ$v@dso9)*4``Dt1M2jB2M1ojnH79>?q)w-fgoj9 z6*n@4EmmM=-kTMDouQtVp01^-iTwThH^M;iA^hc6_+s68MMc;^@wYBc;ICotjUM3ZiXhyl{`T6;VouTQ(?z>JHnEkHwoFk-E;v zuc$b|pB9K>d^+rLp8zm!~VuVKE2+n^7j!aXSrvudSHO$?J!j@28_5@$)A=B+JY{BKVN%#D3doL_FZ3z^>qG1i@U~5s zJ3?U8UY+AYMA4H1GXQiJcU}pG|JK)63OL>(INp18Y-;)>*4j;(fxF?^*=h_7`wqUx z469NAc&(?Thyk^jJr%6fU)WmzLRqVmF&X`cV0DYt7`bTv)-x1B*)t{(>|0zW3^)t=p7W7X$dnuCog{1wr59XF zk%adFrtY7(TaW18e)2>xmFre^YznC3a;Rs#L-Yafh&yOhK}Ch8KM+O9B*@dTad1ck zRAW`|iBHe|+)6fN;n%Hb-+X6peX!E~2T-ZOXQuwnjQFz$Y(9GKB8hijx5DJLG+^%Q z#jzDAR3-xFe1H89AG>*n#mn_OYXD>SBPZ+|sSD)rp}CnBc9;MR?8Mx>;~g0A2cl|{ zRMe)F;DCmZxw$N%teY|{u~^ptYMuA_Gd+-jagoe0hTT2DFX=dMElNKFcESGsegU5& z8`#qUh~O@gY;)d#=P3Y-dw`)(2Hv84we&Y~SDPB@p%YuMiw6cTf+#ReJ>@G&aJW*M z0OIO&{Yys31b7T!-hNrjpNFBWkN9n>a-%Y&y-C1k1>`)O>yr}_On|%!fs)d6$p78r zq8@541@C8oU_Cl95t5NXSE5^q0yYpQNuIgC7mGL1aL5@%K-ELP6=wiQrWZJ@1i*<7 zxa~_|6C{VFXtw7YNx|w+54ENZC7|yB(kif0(ma@})4;H8B{YQTPW04G1dvMa#*`~4xBwvLWD(529=eL)VCM1To=Tfg1{z)|4^4j+6i zKywGijPIFPoUo{n00=#xsCq$=0Pw)WF90)a?GR@`wt*Vf!~vs|LD201Q0z*V$K?Ut zf7S*BK@IMP(M=K{c^of-pN*53&KcaQ6hnQt;>{$8!p?q*CMpWPkV_ zs+Abj699*rE0_XbKsUS<9)=)T6DR%7%bvY}8(-DQ>y=oD z;lGQKyA}ryFKLFDAmS4S0A3c~DyfDwir`5Kv^N+ix!@xZsNSq-Dt^`g*nXHv-g5rC z7>pP9%U zG2lNwmY$&^fg2xIsQ|!o^-wczz;!l_8B+!GvonmCl`f+2)Ejt;w1Bk(mY5{{!}8Iw zF+{zT@p~Zg^dRo1rs5{H1X8RTo<2axAVts|!)iDqyqCcj56EOStR)#_5DaMpQ)VS2 z7|2?f!@@wfZ+8G~F8a5F#<1D9HjoUyko;2k-gw?Y}#bEu=ub+;tw>n zRsr?A7qk_gNPeIjV9^LC5KNS$*Y{a=qjiz|6DoHp&mim9uzmRVXR5GuKV!gU<#*

rt|i)mcZowF#r4aJ*g;V&X`*dfA-dibHuT&wi(Zbo2zB& z)~#o+tzMmE%GlWK_1sBztnlN3;2SiQ-Bd|R9655@`I)1lppKldbUyObV?saLp?Azm zp!!hrvpcbHz+Y|nrn?IK#+4txHD%8~`Tps$zHZ`Ar5L?q*U@50hiBF~)uMhAA}7pZ z(`z;|Y9vO#P1no%K0d|gCRC}I+NGfpJ08Yz;9iCHzb<18^u}3I%Uw3JT*}S{eEV~N zwFO6cS86IhR*g~U*grR-ePPg4@wo{*Jw`-d_J{R^T2f5x+ue#_nOfeYXyySSF@~oU zOS@sg;Ij3oSE@**?ds>+$k;q_m0OELnydAM%nre$Lfc&wNW?vO|FSnRDVj3kfw%Gx7&*s_1zkykH)t8h#JZeLbxQ zc_qePsIj?)-mD!`v!g?68 zwO%>3eovmwx7nqacHh_Ka6nrEas^Ic)66R=!Q)HDEeC28{cS)^0w2~L02i+e=N@X< z81a$Q45Z4wJSMe7ytCZG-57gP*}Vx|#LmNkoZtetLtpu{@MPLbd~6oo;igwR(snR_ zAxtBgzLFCnrTl>rsXlh{JH@H;H^h$9s`%wrKc(owZt;YO6Cq;SjJGu?E#dF*>JXN} z11=6pe-E6y(An9$S`nv}i`@n#j`RxkS*-Sx6AxJV2@YvbT@fKwja(kIi*B8II?4DouiUcgqQ+G7II?j)I zHZlF$oaX)kh=D$$ZxI#4<*D007fvYO>Q}boGcbDUXWc-a{_UZDO^ zyK+L0te}MJR&V(>Rr6+rDCRrtryeb7>@H9yI~asS9piA=1U>`edSBdcKz66;c%c@I z`Gb zon5H^52{l$(nPnara~CurL95$W!G}s7wIACU#=IjKzG+vd>^-y_9t8ol#WNdVkk9` zDrTd!{EL7N%91#@TQoSMx~)DME=Ag)_Jx>ZpjPmsOUu$^kz-bXe$5Gc%aer(;=?Ev zNf=fu37JyG)r5LQU(n4l5Wl=lS^JBk&(2=(N+txme^Z7cs+?Dqf$EXeZk>edx(T?S z5hy7Vy^6fzNa(H-CX=^Iu0`1&A;p<5&DTo}>QkKL!Y%t=|MfUeM2-~s&A~1y ztx%j^N<_oGM_rFdPc5^5@x7521Qptg1njb@Xg>+EXO^|iiG7V-7m>^JmP^|hm)q8a zy>R&U*N_492UVi4o9_HRGZHSPLq689)z*B&NYCY-GMo@@|5u7GiqT(TCY_`Oq_RS!h!I7pAlT$*Io)AyLkxW9fyr>JWQ z%L{PM#pHND?4OMAn#5@*4kAN%jR*Snqs;-g7f4e;q);Lq7EA@I%Q@#l=WTXUXPnu2 zqH}o|-zaCAaeqXq9XnM~GEMR^wXht2XkLR+!Py%c_sh;sJKH0K%Gbjd&Y;j3DXZC$q^(uQX!N;$;hasSW8DI6sWNR2Lbz4IuPy3X>8QHDk>|JrOzP= zCejZ5j?&Lbwli^k|6+@mKKY|wD*UcqH7dmUK~=9fxiw{rlj!OSw;rdgDJaFhB3K%2 zNr|)+Yo$sP+fa-hnU5EV)6RaFS*I#Zz&)sFnp7Kt zt@+p&xdIVHN+Oe9KzS$?B}~HWP@1%-YRS)Y;djf<#YHl8*%2L*3Zd)h0-V1qW)e|` zaTU?C(J8r$P^Vu!miD}j64~92maY9XvOV}$HXUOG_?A{)oAXzL!%vK@FmQMxo-my- zaQC0&uxW&k*yA)P9X2%?@{w8TJ<_Tk2{5hGaR>F&_kL~Xn~ti-iJ?@Kmg0k3NoSoy zzqT+e$qNnv&1WA68Q@%S$Cd--We+Y z?zjyOF1h!M5n*hrt@R!Tk`#|_qG>WpBY3L$l6ZQ%R=cYF1~Th-WDJZ`!gC9^EN(@8 zebE8@)0{Dq5mvb4&7E<#$=EhV{*z$x19AwbUDv)};>nZo zeERkj@gZC<^E{J2bhm^WVq}*8IT!4y0P;i#2v!0HKP)G=UEYF#7Lqb~5G3$ooU;X# z48T`&*d+xE0m90oKsImn!6lk0`ObCP;IY7?XEmr5xt!eWW~!5wwQeEA(G=!jr!9_* zpAiutwO&IP#H?r(w(tq#bN!A}RZ<|hrJuk&+KP7v!}3C)?-MkR#lysIrLPvjbuK1r!tcm6eDJVyv&db+2xEP8@;NpX{L_gqSXjB@U_@FH^53U}`=e`QuHr?;}tm z0EgXD2;+K*3R;{Zd-cB<`2KSP16-kX^`$k%`GpolMBhr9q}BRtuYS-)Px5CYG68v1QJ zg+{cy|8YA7`X@1TJWjqBk925(F)uDWJM)6Qb^flhA!%5uYp<-iBcZJ5{tGP@{QL*l-aTzQ;PpMFz6$k9KcOYU;OVCyK;y0cW&rV+d3+O(XWj77AYD7v; zaPf^$q^-7Cya0zjWp4(3wP-|%za}xVKX8gZR@Fg!k>c&gV+e30X2Q}6QS&5`@^xRd)6VJ*wq2edVJw(@yfQZ;XPK381Jx4 zx@hjkv*h#U)I~PxSq0xyjGNd<2is)QNM3%z(1+12ceJDlns&chjM03)CBtOFf#CMIt2Mjq$+cNIhwrO7RUG*HPGrRv(dlmJ+h0IK#(hNpX(m+XwSF!< zj%c&b($eZ6Rx#3CQBF;9he%s5&2?_pDdM}%*^=H1J;@p;k0^Av21hois``1Ro0)lb zkRoo0o8gWrkcMIxzRcC`(>PqLzYJn%Eu;8D4lZe*8Z8Tnil%Za3+-Sk!=#&URj!(>(oA=MLg<371!8g+`FE^ zYOmdK312PKw4GDMRf>SD3Mv%idY|2r4|=pDU0l95h`YY^f6rV<`KkuRiYi9*w6s?| zp$SM5xU>QvpsK%k(;oeHx6i5T18QF|8j`u{f4Kk}Nh7G`y`_9!R(9JYc2H?Fu8Jyx zQs0TwoW5PQcK?01^+T{%ULI_2YAT8)?!MR^L${Yjr;D|4m(ym4^zuSd3HG|CYi?fl zklX2W7qOHiAcp%%EFW69_k8<4dy_^S|HVm>x6&k+lVkRsxf1lE`MA^#g~J#`38nYF zb2o4cPEf0hN5pw5?u7YOW2x%1WZz#=+7wqK| z(mLEI&vql>mC8?{DB z(d3U1y>#uJ9CK

1(CPpo*4sCG&AhZq>6tcUB>cPGWZ(K~lXJc_Bq)PJW)bVzN&} zkw!8~ZxcoA;-k1sDAK1^1LcCkQ*268@ORme4O#*n2TJB9Tqd_pB1Pw78?tHNSc7^JVMl7ePNeGDp71rPqQBf`HVwpl!D0oXx@Pq)WRLe zO1B8mZjquaZ$5rbQ^gwZl0XW+w!GE-iCnitue}8te6#|k5et4@8K`@A+kQ+Xs~)4X zL+w_0ZG8WYfeNCxd3siSD2yj*cI4>T&tV9ATsKrGJK+J;PavB+S1)RN5@}^gQORfd z3E@WauaWlcLgs%to=Kjns@Xp%WlWVNyd3Lr5q1o13VmVK0^uoHJn3K^v#gHm8OPxf z6JwpUVZJs@$3x3b-AL2T!h&XG>xd$GGpcyNFbsN~K24?M@MxK2BST8=+Zqp`wdMXe;;e`llBmgtk4f?DXOg4x30ZYsE1b&H~~ zI;4S@io$~|UeaD~{S?*lcPUwg9D9Ur34P$b99}Bz`(71O!c#B()lCiJT&;a^raP!o zPCket*H>`wvpq`@)S}t%;d`sN`iNb9ri@2Bjpt8Xn6`HuypZVw3$SPs z`UMU7l7FK_NG9wkZ$v+CkN!6`S2RDfKqJ(~;NB(L=qyFlrom<6Fcp-{FpujnoiB@It=QjE%ik&ZKyS7XxH0iD!zO6C1A{Y9R-^y_AZt)pdbBal6rss+o zVD`;_6?c9}V3GwF28pA&Le^2#Xm@AozpF;k6J?=blt>~EiiDZ8+cNNpf4JkKY31#k(fQc$T%#Kz1$O1wu{_q$HCEIaOPW$4nPl|D!1oVAH9d9Qi# zvYLfm4@WnNdAo<*Rb%HDN?o4Ht;*;`(^+kL;@|NYY8NT2VREawjKm zZ26_9T`6CceQj#4-P;t4WUatDg#BrlzkU&{dR52$$<$3M;6rmCTCp{$Q+O@cS<07E z2XYK3G^+DdQ-1aR@pd!hk;fD;Ex=!!PBS7z%gl{t>puL1YLn~H&l$z_TeHV3k`n37 zZe4h*7%_O{LnLY*vigXwrX))3X-0olMhot=wffsT9M03wVKo##QHkE-csw88WdBvx z{4WKK&09N2$^4Q3@)udw^~K)G=@NoZ?EhIH?XmIw@76nb_;IRVQKf-sTp2-yu`rsX z^Br7Q0+SIkrxIs{8n0IC9Ws;F?{&Q4t}6phy%ePKcjVon_y*G4Tvn#*JI^@Jn6io@91K-z3 zoV%n1OJBk9W}aE{mqHSU6!RBKnh(pDkTO=4bIh=hBZHHbxA@fU56E_E)YjF1rnyzt z;Oa$}7oc)VVzBG&B#6msor=t;)6T9&=AXoopW5+DEB;=o+2HfXdn=)pQ1TWr>cN#x z-h5z|@tGp`;1+>Ch8c=VO#(AjQoFjlAI`L4vP+dRcl)Ka>dn01foe6B2}k_I+rzb$ zekD>^pg?bw(M>Y3G2eZ()WZ}o(ZDxpkup{9Vk^Eu3{F)RaV-U}TBf>JHu*;Il%jI7 znTc>|*Bcl}S%Qy{>LR^~;|CGrJsdu}aS{pM>{2b}T-q8Pjlvsdo8b#Q@uDH1r*7Yf zdG5x!&)SWchqk}HzjG0>Q-jcbiSqT4v?s7Z+4!2pWPhW15P%(64b7d>0{LZsq z!VmLXIh&vn1|)QuY>U-2|^ zLr;p4us`{LRFgugynqm{q7hYeq+a@37ImEld-uNBUs+V8i8HMqc}%-8n=5zOY;4=| zS6xy&x;HdF@Y{CB@-hz~wDFu1SVYTJlJV@=TYV2Vjs!Uj`_0h`U*fL>02lvB-jHr> z@$iweS+hCk^QuZRC(Se6uZppt{~VEyabK@ihQpod)4fQ;oN;v7WL4tKk5y(57+Q7B zgCF?)e)l4itayf|xoeiIdhW|+1P|%1tpHC^^@@;adFjBKP~u~CF>AswQ6a1;4Q2Wo z)t4@mkGf1zkN7QIR6ihVH`ERG03&C7B)soJmwY8$)PsV){XLh7E2ZA~O6nVVsei>Q zCLqzi_QyJi9jG!E2uf~6M(I?7j!fD$%aZ_?^OFHx?y5l#bM~$6aR!5Vgz9?~u~jnn z7m7%agK>(}$S@um%+VUDB1!wN;GQBH;e+ywO?Jn>Lhr}fS(R#oeM#ek5d;~p%&0r0 z1}3l_3c_r-AuMbRXI0U;ys7&hLH3{NRO}JaUW=$qFV&GY(lh6U5gMpwK#OotqLt|c}uT5yLe<{mimSsDm?C6lsaF;SfPcxzf!Oo(lg zqp11>ZN;1~asb`E0RPSHs^D)aUnz%qzT;ZEU}}=vvu{nixGTO*b=??nV(lai*d1)p z^;uXEqydSaSE4Po7j+Yagjg3&^$If{3_O+pc>vl^2K zyp;0pZ?5XwbZGRHS$Ai(cB+moyM{PLXtRu?GCfi%`|g#975*Bpbe|P?r*rXAm`~fx z{;n!#SXY5(PW!Nt&Vv$R?Q(x1?OyiTJ`+PZE|=Y-d6iqT2R>RShDUUFRiZtcpY%{Z z{NeYZserw?+kK0J3w@i)_7;@owSxVYPJ-n1UFg5nZXbLXf$lc}UoyZL6Jcfa*o0sA zyg|jsQnGkU9T!!cpMWBeg%J&E)!%XT|M{}&A)nL@ru-*Ff<@(w*VpxU#JaV|COI;F z!BMRZyQ?KK(Vp3=mr}ym4kwoBE76PFrndcEyI)l&Ooi(&g8d`(4c*w9;`<#|_>z7U z5|!EzdAItoJSs(qa6zP5<-e=24LJ)phSOE0BVW=M`uC6wXL|J2r=F3ZTP9&A(}%T< zelx4X3$BO(ehyF#hXT7dx+5vgL4|*cr2B2)!3ZMdIHLRoTy-EJqEHTm*9!#Aa@({f z3I})7K~^3$D=q@M(5Ok3q#*U?LErXh7Fm@U{HRo$%bBgREQV|UqpUteJ9^VC?gtgs z0e^r>ocq>v6;Fe1cpz5TM;xW#Kx5J{j~vSIawCpf%3{)0Ic-iBW-4Bt(b6n3=gh$K zl^Io1EGf_Mz4!*^MRw!4L55A&1Dd;*qgCn=vt1#9@7o`*R4+m7A}w9`pTutM0^(=< zir&V(fo*N~;8!1_AU`+=4JvI`lMEcc1Y?s{fQo+$OUq?|Q5d_T{ZA(1CpnNt^6$l* zTU3fhxVifjZCi!fMh2JGhj+u8sH({XO#Dv?K3C#f{f5j5vSs*^{7J(arYCU7qzp^M zCQKhCMUXOj=j@_p5T(!U2-QYuhYoiXiRq^h*8J;*PK5pYDk*nU8Qg$0q>W~Hb zlpda*<{*zFMaH+z2RsMXYOVx6{)IL>J3Bc4KvPTWHC(JY4YHv&v66ok(7pVP27#&u z9#Mx1V72wqC)0?wsXcRvwX~uzDM1y;f#yF&eq7B+Dav$dPDr)YuTtXyg}PTFR{7x} zx+6PigZNClqtrQ_(YGcoi>1rtI_an8QFR|>aTFs&4)qJowc}bStlsc${U4^zIx4C! zeD`ODPU(;i32EsD5d>-J9vYEGN@_$phAv4#DUk+&0T~(;sMuTOI~XWy@lq^4!xnU z{{B0SIU8-1ZuJzrIukz%Zb`U1>7@RdYM5r5;?%{u1DrL97=k%C}XV~$q+TY7P=1g_EGG20`S7oYFh1-?f$;Q_`#*V&JPlQ+YIIWYSI?J z%gz766_o_ zevx-eGYhDmZ7wP+Dx{LgDRn^22q$-l$f3d}r9^PC5;)pkXut&$ppn^!uE6Bc*!UO(kDa zM<4ywPhP#h^+%X$mA=oszQ<5oiA(KCvG3N1k}GdHgaX6af;?&#=fWwAqG?pc>j zloGx%WkQsHsUHtTQZg!mBil6JR^Q78Fy81? z*MbAISUu>%txUdr@miTbp42Ef>ejl2%B*A4$idHE;1YZ0Ur z|Ci?uI<&Jsp9CNzaQu=3w$H=eE0FnobGP%DnUqN;5=h*ZqZvGufPBgWy)98*2p@}y zdmWDDdfUF@!$q(eVJoYq5Dx$QD>u6f7Q_g(s1Kw{E35If-Sl4qrDsSq`+yb+;tP)i zC26tFw1}XB@}SOfV1z|9y&`6;T;e(!<3raNCv^sn%ch*FChE&`Of!y<;918AOI?bVHpsAqR_=H>l1sXV$a{I3I<>*~m^fiytda_yV{@6$-nGcv$>^_VbD z>7e{kIFIWqeG{QkW;3vzDl0pPKFQHcA8kj@{Q)F^v4SRr|P~RuJZr z-_Wf05%E0zt-Hemx9Km;$;a$pHbNQ0cK}6%iPilf20Jwoi`!h}_{MKw05f-!nnVH- z%InKg&hdLz@4rS&`ag4{(sEB|qpAxdDLzSTaNtFVG9lFtb?SrMB22C{WT4EGzI=9$ zR4^m%Eo-)r<^9qVvz*j!FaLmgJuB$HllT>NnP_T2@UHdU{_KZ_>X{(dwLm6Zo%$pH zot1i-8BVxrus6pGCvEdju6DH4r8~Y$zRr5i3@k+mLa7hO8s}}LGD_Uw!2e2ocdH_a zLA1dmyt8tJ(H@O-VZ=OaRQguR>gM@Kr*8+ zwzb;rh8cYa!K_2EvLk6uTGf)#OalHDq8szaajA`Ol18(hbX5M+D4n>gw^jDECW1sk z?X&#;ek$|?PJz5PO@UIMG!G(#_MQWq+|0VNeXyiI6YWBG)Mtn>gVJ{!G&pJzk3h6k65+aW2%aOW4nB7!9AD^q@{e0vcvTB)HLSAPVn3m zf0m;u>(pu`-(oIGnWrBLYN)n&!ryyF3TVD33uy``eb`w)XiZ)E6|{2)r@wWM6j61| zl$YduKV(|aVywolRM1h@(tq22PX<{pVszvH!w&|KfmZ^-t`=kSOv*KiSd^dsh&&!f zK%l6ik@mf5yV$MkDH;S!+K%@>c3Y=Jrg6TpiFex*ef}MSpIb=oCp}v#f(@*d(~A>% zl_S5?BR|kmjtPaA=Gsl4ocx}Ydhx)`TuK#Bv%CtHC|Q4rVcMOpHw)8d3W3*k$`U*6+>$|BCtFLhT~b2E|2&to0Bwax0-FtVD5x36-fp8!$yYQow~Xa?qwN{%*7ht0G;2)(=X$ zX`M01hWAASRPi4XQ6kq-{BS5=Kc5IlKwM~LNHwoOXJ=72hi`yUNn-jIwSa#2+CghN zN{xK(MvipU2ekm<{V&uMzIL-uMKPf7`74fYEXcdQUqc9I*rt?6V|v8kVw@hsJ^e7r zJ&vT4+E>@lN~AU?P^XQIDAv2weG_>CJ^|x{(wh>PkU`(7_1g+*gVvxZa{wOkBKm!P zJ9fW(J=6ZS?M2=HIMK7(_U&}|tH0L>sGF)=ysvl!RKZNf5OmFNfb>GE3@XSYNjwC# z=w{%3j2MG_(&F2bcSO^`+hxL?bJD}Au&N8JnVbVDlFy{0)u47Q;O#uQ_qCcA?T<63 zvGJrg<8a%QKIjp(I;sQu+dR0{dGLF#(%gnd774-&;;cWt*9|#&FnT_ia`FHEdkrJ7d?g%+98ZF^_0Dhb>-6nch6*HzvUp)!gN%x zKMKDxx&3%hRcC=Sn{mK^uuQ;K^G{d!Ql38k-fJ@h$<>_Dg(+678dJboEcoeajkCXMUTfiv7f@id|Bfj?#`+`9p1*hs-FDXl|cp5RkA4JT&$y)K|0-}Po(j`TzP%F zke1wZI-`;4k=9|`rX>#)q0P$sHPU5)g^*2RUC&qFsSmgldL!9aEb z5e~!M+-|^Q7jRMCYQ3>C$fx4L{Nlrnl=zzD*Mt5gP9Ke(Af6tRyj@)#Q7bpE^Hi?a zWMVRR9a*`TA9Av!--vdQg9nFjBPd~AmB!DVSl;enzkJ!wUjJ|Q4f53f!p9gTIU+ya za<2f@zRIKKukKD1?R&yH{aALi%ptUQjtXSmMGdOnZZ?sYF_8#GUVwJQKWz7I;jsN_ zQV8>w7S+@5OrmU>wm*4x7Mc8=nYgq@mE64)Z@upV6tD!WY4%BQJCj9*RRqV!X! z-%PFiYIywXjtBgMJ`QX6XX}OnemZtvbWZU&s8J$5!p0W}?iEkUKdlCc2?eHlzJ@P; z`{VVM2b0z1=*&S1^?q*&_q`_p1K3$J#VckVVP7|99rMVJqXo0FyDVA#8dLsAk)y1v zP2J>}Ms;2x4U&l9ZVCa@PGwXU-|Ow$p5T7sOIG1UT+gq+1oN|DjD7;j07$W1Zj#uq zyPk*xb-cf5eKGPgsx2I-$c9GaY0b^v6BDCkf(Bp6w4o-*{Gfz}v~kGk_PY;0!*ik$ z1n9m{B)#-qUc?G6+0TH&jr|i^a8JJAnrJzQ#rwxZIH*)*PWI926y8*NFXpuSYWYdR zx_EN>X7l0$SweY}Lfo5YGG?kPn&K_04%yiwRLwa(_4{HhQ`_LiG?q&RSYI9&__My7 zcQ$?^$=YKL)~nD6lNM$i=zM0F96wP0+W49ViMnC3CpL&gH4z!LGC zKc-;ke)!-Wj|So7oaDc6|EmSaUNTl?(7h>@1()!EFL^^V55?SF!nr{{%qLirZGzC- zFIwU^qlEfDURkz1#ITp@yY1eV@Ta<<_h^-twQuXIQJOQW5@Rd&tH&N;%BFi~9lZIC zR9%;xB{|S|fm@E2*`MeEY8Ud=zV%_{pw+4>FMBhPooTeg-toT;n3RsAimU{SG)^|r z$v~6O;ZbcOA1+iqXCm}OjtiO^%DxwWkYTt(fSOafIZdSx3dX-B%&G9V@McTEEP(5Q z60su3%bx!Hd*yH7QWwiMhkmw4Y) zs@|*VDb753Fdg4&SBc8rWLG^z2j%*P%I%cEuV$UNgYzw z*4N$SLHB^Z+zl63Dg*on=x>n)j3v)aI9+z_@856ncO%j*f6KHywL?|cen60oHS`{9 zOFlsEvwkJ29ChdG={eRo_OjT&T4H5RHpbMwYn%Tc-73TELCGGby~`3YLl4$aUkeANQ?PD@gBpS zND(~jphHJ2%OM+x&3YEmY)e|~Bd36u;{1tE$?TPEI|pyxgYjsJ!YCBCygS;r2fOOf zkWqKP!>5B7>ISFYnkVnnFJ37z<*Cba>E8E`U+GVFufGh5`Z_0(L`0b5RY;PTEcKW$ zNNt`c8E}}Y^wp*ZUviN-v*M*Js--!cZ)GwCdi0V%xz9<=Hyrs1r{gq;gywx~(i;IxPu-0pE)NmvP62fkB{aU4r@wjv48<1DrFB4BuR0?p8>AM||5)K%A02;u%u z+~{VGU-{{|*GCi!p-- z0kt*I@XuP9*fdV0euWHj0G^VB^hW>vy*2BR%|+tYF}eQ=IaMDx&{83elp|$zo*r{y zgPqK!qx4+YIX&LUjP%EOx+NTTBEd>ajSg8KL9ILrAIvLhQ9-~$)^+~2-Aw9$*@+Sd zsgfsi;Kq*kzYbgQevQO=__1Ss#ATHPQ>38VoTt6KJA@;5@hTGvS&rit_E}ZOhQqFX zlao!4(w8o5p8mR($RJ4y`XTAux&r`082W;JTOO6~Hs~IHN|i0&LDOYgPpLUS1#or? zG?st9!~_O?(*HI%2iO`5wZ{HytyeEJa%5jqL0%usRJ=xwC?P2MeOGdD3$y)UAti*I zBG_SB*;>N+xH8e~6RJEe7QoK`VJ~Qw=K~RE=M8j0InksGEt%T0P3UGh((Iq)V3BRB zRAmw?FE*l@AhL8Ix*S({Xczi3p+j$yyEj6;wHW((Jx;)~gQ$Lggel|@HNCP=lGUOK z{rT>{FJ~gUj?$DPdY$oB^4Z?_)uZsDD$QK(P{FT8!=iuGv|&WsPwm#*M8g*L(e;%E z7y!7ozOm58u{D@=69eq<&(Gdp0|K@JhnPD60?V}3iw$5OqYp-8ZzRhf(OtJ!fYrLB zwY9tz@By%D&^t2M@G284Kk92j@Z*{`T#I)e{jW6ApWYgnU+b0U_F7+1P_Ymu2IDA~ zBSDyRB3Ug?{7`Rzug!I|fN5}LdnZvgKw=Bk{afzfJ&oLGCQ~i2 zOC&HncCv9=p%S)vkb7Utbw#rPt`oWh1g=;%4OHp#jZENhc9v~o)SQdfK1U_W%Onb z*v0JuaJrMi-V->ME;}xTnFG&=)b9xUI*>D;ta1K#)WwhH0NyMz=Fo7K!qauzCs7>N zqcr+Z)%^NLH^u41StZ<^E49(CPkXKV32iI;or<1-@o$uKl_@GLRlZ(BA9oG4x=T)9 zisx)UoRa_Q*c-Uu&GnJsAyoIaJx&F>eHh;RvLKQ zR3pv%6a3!S*_qbMaxNF_Xi3XW&w|fpO;sTu*-;4sj6a+5hhqhuL(6R~7+=?@64Bs) zBbTa$P%eBXsOkg@3XjGkaw>`&tdVCBq^&zozYWxY$~7x2 z3R9?<1%np*v)Z;5zi(AGGWNn{##F_!{7d1If4;-tcUe?kX}VpiuWk|EBI&qr@%?CD zrBa-5$8R@%$=hB}!iKx9{YD)VDJui~?Kh|g!PHav5VJ)z=><%aTy=>9A9;d<)vgAj zJCQ*Jl$TLbhjeuAAsq!RayXkEg8G{^R#`MZEn*F%xi=^Cs64oJr$rPL=KpE@wrT-} zX+Uq!M-?Ds6q*?xFKDS5Azn{D~z6bjo)-@A-{UYFX(fV4?pV1Z}}!4A!`kN%8Dl!VJL)YPxc61 z;_mR5@b|z&2H-Yg(BhQ@v&L!v%~XewR+GxBZ<<$kK)?kaoLW5ToM`U&jT~QuTR9=67t|jD@$RpOVo8i+Bt9ISe!^vgNB>AX-0e?%fJ1Fza zfTUsV7raFlC;v?uGAViR8+O{S!#@8R;gg8%FZ@r`#8W$Xf9jy7Ia$|y_6Z-L8CBmu zdk0p**8ZB}Vy&O~u{+p3c_`h?+f*8`@_5~*rCW!SS>smdkGtEtu!H_Vp+W5j87xJ%>x{56XnC+Q&EAf0aUoI}LDVYGPGkNIcP0TDu$DEAk=eHj` z@ZhaO? z$A$1n*~@ln?wJTWi?jslt035MiUV{qMZL1uTG{ibEO!HBgceHv0xQ*cZMtSjO1xNCj@ut z=+F-_Y&xm^0(u66gR=9_75hI*rj@7LuO=Psx2_Ww4%FY7uCe1Nk3F89lD``xiNiln z|LIskLdH6Jki_P%?AXvtL-bKlf%C{$EVL^$wg1p+XpitgPna7e5{#|Ar-ubeiHXGH z6A9bDJ|>K=_X3_~{nJC&LkW$vI`DP)g*e}@h3nUyOQ;S2KjzMgbeR7g8S-~tGHNkc zUqcW-&p!n$EukIjdaSF1ZW!3T0zv(A*)AVk>!NVB5&jR_p?V}9?}%=M1`|3ZF-Nu_ z=qBV3z4IwsS6oCrt`>C3i8W9DOtQ6l|4`qt4}T0;Rlw8*{F0gIwDNBpX^E08On8u` zh=z^V<6P6S*zE$1zT2hsofdvmYQ%SFs_SjPfEoT9@<-dT3(@nx=KC;El7uG8X(V{1 z?FUIS2oT#Dpa~3^6Lvqf7aVVPWizUlMx{H~S|P(bVzxn+W74yFTiayyh8J04K9HG- zchyQ<8Phgh5LV8w#Gi@D0nX$)%8ynDs_MPNE)5c_j@o_N%X#n5zq(9(6h10?@}iD* zrX_+8;D<*FWPyj(ne*aey(k>@_3J}4a}6F5tze(F{M|3{OPmOpiK$wf>yao;3d+xqA^rDM%-Z4?JAe z^S*Gfv5x#ICwJ+&PFIFa`QeY_;=sZE#c%~FB=-bUPiA`bQLD_wi5?cPIj7JLVX_D5 zC6U5hnd>XSCtdfuHEiW`s>{E3{ECRysGG&wu-~#sS3Iw*(LF}&4rU;Aky*jUM&uEOIo*I`)ej>&#no!DyN7J5(nd_T~JXrwIycYhy6WVbIK5dCH9 zvx#Evrs?$A?{2Hzp-?ss2?C_^q6~O85nyV=4&Sr_S<-cW3M*OcF4gwdRv!R9cz-uC znL6~6Zw1T*eH)u(8_-n-T8#cF=b|J47Bg}G3&<@IGfbj#LamE%ew+7!`dw-c&+yr2v1qm)KnXrufITNs@bWCOO~g(;Xdm z^*V(Db_TGz;fEOsOoPR)^CJdFft6~8vHrULYG_%%{zj0~f~GkdN0l4}z&xvbv)j#- znwDO3^tY`r2bJy!45O0@52=KATUp&h2A0kA>bwilR}1sG_XqXoPA@Iy-u`%R@x#s{ z&3!L|*3HdpIIT{~K9Cp&kH&U%yTRG9Iw`b8z^smRz!0o|Cf1~br_-l9Ao4qixri;` z*edFf(7ir4+(dVJp)_{^fy&ol?cYt_v@n{p`qcTZ>U+!Iv;0yl*MMvf6QN^^g@v^k z1Gme7od`xkShCD}0?aktU030h^P2KAw>g88O_Jki_Z6bXonA$|qVR zNRVoF8s~NO=u#;EwAVlMC|58>IK@9L)U&Lj6He8qfe{X;tK2n(L5|QzMZFQc_$S zCvT8UcUweCg$id3X+xQX?~@i=DA)=Y7Yhr;K3i&Cq4m!&(evT1NVs(@t2p71`S#O5 zM!Q0QbIgRaZE}5p?%X<6-o%=@;j+DH#ygEf<_kN?hcORZJQACU{5;w5)l2Zh3K zatLda{NqsWEr+v=PY+aC+j=L=FxuMM*USK~)gkh0M`9o#7JGy5<&I~4fu66bC_G+j zAlKR;LX*~obb+MHRn^7~+x>oN;MOPzEO=)vtjY6^Oqe#pbR7nqeU{sT(3^zM-;}DY zlnSLW#D&JfZZuw%Gn`~7GT~7=DK+xJ{lGWK=B{6YawYjDS0!(?0_SuTtNh!N2iv?~ z%=k$~S4gL?CmeQ&p4nGcxmL;j+#XyL^q-UPI|yZbgOu*~-E2D$DgR>YzUI3g6T3vB zi&YkcVXTd^*_GWPy0O$m67FM9aJBy#CT3t}7H6OvL^nyyR9(7hA_t2GRE%Nh1#Y4*Ut7rPSHOhxRyGxG_4vwYSo9(o-~Q> z6jj3@d*4W8ci`29JZi-|9{Whro`WEptU28NZ}c)jCsX?Gee9kb4>tX61@;p4#o1HV zfgtRx)7kwsdzqh9A1wp-qL9*paOwGg4(jm^(S{`Jyh}rB%FC5ESn6a>ABhh3FebFc z9l|%z$$v}vWD=!G99xh7oX+~jeg@nsf(cg<@>MTS;<1-f@c1(-vgkrQ7X$pS-)jjZ z6d5t8fp@PW{BFRw-yrReR%ZS^dDg!)88xSZk(0-vzZ>K|+h)46^ypxr#fss+Cc*rd zjMG3J>iCWl6?AfV$q-M2`sC59Si-G&OOf&kb`!S6@k8%c(|qYDIXjBI@)YH@HkoLo zv{1wF;0iBlIbiGLcHzZsY>snoU~yL#rV%$<8huxozuZ>YZ%5;_%QW88^Z9%K?BIlDSFAf<4oB$pw$ZI4`$f8FsHxD<&i~+n6oc|0WtIC^ z2}(B|10_0DIU@|5wAk8cYjK*yE1D_jE@o@tPXyr~%NNlHP0yBc9a6pJYV`8Y9(QvD z2BoimUb+JF>u0d?>TFyrT@6`XN+6ANKp2%89UB7u(ey|r5+46MAU<~Tr*uFF8^V8# zE9KMXit5aw7d$S~HTNUOlQA`Er1Wz}CHCJ`>oI!)MUQ(XMqh9VTte>bLw|1ybj-BZ z^8H-!4=HPovGsi(^X(RY?|4CUqq=ujK zhA>AoANypcr03D>P0w3gi{bdK>Z1Xxl@342Mj6>Wtz+&eQhgabpX#K!yN*)n!w19U z&sGVA-K_NRZ5lX=+pv>2a9|p{_Zw2Ntz_*?iu(SqA~V`|>A|5LZ6#$Bu3mlr0@MpbCybW3#KZgKwlN4sEznYIVB#TG)<@WmA zKzIF0TU#7*iAtLK*+wWNJ&Vn8JlrMBh3JB9;IP7%9u-=07dnI#;&!pKIV#Qs=YaP~ z@IX2y={_RPk|$6I4(pr#xF>P{$B#w-oqK9~hlA!)#GkT z_)fM|Kv)OOGh~S*fORo@GVshdlTQJ6n73v3KZ5s-wfA*HJd=q_C3ukJ#jAWeUUN(2nseNw{>p z#BWi;2j?;Yci0feD{akcUr~eliyW(ge+p!~3i(o*{=4I%uA2JhWOr@c?ghZiE*A(W zH6zs~>Ycu9(*Q+3aiOU>#;AVl2fj>?G+vJFUwG38*f-I$QRkz0UVFcoGTTPTJ4@;i&>AF4P(}e}>Gcsr( zJ6GA40&VVM_$ci6~5BW2s)~r zR&VR#I)jho^%Y&=(XlP|cGf3K++;ruhxGIvBHiAa`+zA6IuJjPRL}8YzS?l|0Lra< zF@UV>{(u0p0>yhi*^i9R(bkVs*LnF~+WRNf(*UYm%(6xy-Y{u#W8fOmlm;GO+3pf)>Y}F2{Tz;g_5&;8R4p;wz}dDB%brq>xBIx|{$~VU=AR z7Y@iw-RvGDqYKvvv-OQcn&nJ-aO}{k>qyp)YEyYfm23C(?F3#>Eb(^zgWB7e*c<&} z1_>-^cCxtlq0|4Wi$I$UP@krf4Pow|m-LXkJ-q^nzYf|cAJH22c`vaogMUJ6#U3JNv z=CMEjp8gV{w1C-?&`n+641aFJC6iyo=Xv|N#meKT z!FF#Z&MJJgKueRkkmyFuZN!M21+z%y+tJhJ4cb6~KPs@!ZEEU>M8W?qwHFWar|%}D zR3vcTJSy%(L@Jto8#oGz{77ES3Vy{(>Aq`7|^VH`@NU*oQB+ikrBO$WG>MsOO9CNiF4} zgtyevWJ?FsWR^2N7G)W|q2K4bb3QUS?rWwx{qcmeFu`g*;Az55)UJc7eg*v#6i^c~ zX3utnaM6fAVp{-j*GMz(HA=MdE+rDL-GE3rk<4o?5g~g!y%Umsi_N;uh*0CL#mdMwLB)5>@a?GK zSpx9--&j0f@nxupi&=$`-Ama2>^bNglgu{66{(bl93OzN&St(lCUjmm?Bfk}UtHC6 z4pF!$4F8&hm3)(MldyXhJ&O!ZzHhs~MN^XXi*MuWZEP}MD);_hb*7uoih|sb((7Yu z%_%M=?k25{3gzrkG&NDi`F0a;Pt;%8w)!is`~Kf2K|%?SRPfq8-9ZG~;>9DjMN)d$ z&|M78TZCzFPF0E5K$qD|U3fhiGA{CV2G$TtQcBL@7mpf_B@?D|5B!p%D_#kkQ5F*_ z8b=se5!8!1&5y$Kcqyol_r(Q!drjtah3Duu{Cn`v%9FMg*}{#eh)_!RZr|+nwZcCM z&PUeCNu0tIx(=2(^<6J{l$3{pPb|f-hl}v|*9OwCq}oBG8{GsQX)znLLO^~0jME27 z*z$#vxjf(SuYHDT8N)2_B@Uwd4RDBaqVBKl|HJO8p0ngH9W5-5t*6Woxpcjckt6Ag zxxFaI5&qTHxL#d9@-*y-U2O!`Bs5Jn8{4Kbe!f#V*Tz647vQnlfB_ zD&8jkI(+9d@EK)=$u`z3QDyazM*0*HQ6Njf!Nbnd-XAiH+8Vp0qqt25z{eEcc*62{ zdPYX_{FbDZFC_n03t%-IAXl!zI!rxGblb9X8CrkD+RnkXK{18aiI7G<70j}b`hYIY zdAQG?GUyZ_H#W-~5Tm^p=yp5dhAQjY;S{0aOY+FIS%ctpGdV4nPKhZZ`mWIemMEH~ z;-z@=sEYp40$(X>p@AiVV((mB&O7FJ978#w`_$km-x@8AyioD(3}1N3{Fe2V%oKMD zZ6~3W0G{w2gjmd40JQxE`B}R-=xIM}WCZxEtgKmZ&bJdJ4h@RR)jHqd+F5 zcZ8l(=Nc0$-iF(*ne|!@Q(e2~KHE`mkXG1v=gwK!i!MMQNH;;XilbjZhE|jS0$^E9D3X563lPoDkG4aMh3T5_k?#j{~`h5vg^Ya8px45Q#^-j2f(%%(nx4ureZxl>G?lszu#%7ZJ`XDnE zjl5+g@+Rq<<^J&Roj-1Uvo9O#*mT~$@cATz-1%M*VQZ{uD`XsqP$>LZD_Rt`p>KlM0xm9xhpxT2?o*EKd$3J3`B@bRgr zs^VNnKm9#ClvY^VY8;I&;o()h=BnL))ZbShIYy*Mi7a|DM4&_i-_at11srknA*Bzl zAbY#mDkmAUr}#Mrb+8em{->JsBS{|^3mdlhcr2TL8yoSP%>8PWfbDOs>I~@==^2?y z{G{2%8k%^gJIU+ABgeDMF9+ol2{>5sqf02Ja98~Wv;MaEZZy*qQ#f94? z`idA7DfvkvD9Qzc&>eOiUmQOYr7{gfLT&)e3;5MBVAn|wc*?|}q9tw#s#-P4X8;KI z3>aWG6|WiHJk?oJKU}JB!=aFTsxU36P`KW3L#i78l{12lqij}c4Kk+QmqnhI%{KnM zOLi1oiHsfdB@89ywP<&Pr%OPo=bXyi19rX?eq4|N-8bNh)T(=w+)P1eCz$ovMw#VJ zhN^1BX`pc>UQg{CO=Z*$L`Xw4bp*?h3ox#IT=dAjBM7Yn|4x_8!=(dTx8vIRpY_ltR7nFszzK%u?1^4?YAOiX}NOl*^+tR8?B0Ll5|q z9p7uO2BPx5rR3-d_k8=A}Z|Q&Vulior&^;PVXO6}s?1`W5 zPld(|u<|kMS`0p5?{0$EyMRa6nX>DLj? zPou&~p>fKF3ZGPOE7#wK$q$iT&?`-~O&2r3GG}`}=0=flWlIe6b9 z92k`#o)4o~!UzVq;1bn6B-o~9LVp(_87lKVkFD2b5wu1|1|zbRQ7{~yp%t{X-|gF} zzR)&KafN2zrGLPaNpB7sk^-1{riW|nlxeL?YP>;?-sSCiY4N3btX1IOb|7$0kT4#0xEzij36x+Ntzpj_A|0&8x{elg zqI&e|Dg+P~zYcbRy7P4|*mVFoSV>O$g~Fa`=SkR3)L6A7o4xa3diwO~-=fg7BaN%j zy~d%ONrh z#Az!{xV0oBEuC@aG)HF%YzyjGsjy`Qe%a;UTIPc<=MLVY@Yc;OuA?#=zO_t+Mm8F7 zGQxG3KUEdoqST0o@05@XHvb~cLn3@#4D7+iczwiLFFqPMb=J2jIwR!0w3LB+*CE;d z2dv{)S>APvXz=Vhcmd3`jg1ZUwonHL?1F-ViZ^doE`v3*_+e1UpAX!CN!|W4^v}UT zYLU>|GT9(rxtNjNH=Cp1U86PGpYWP*;!V(J`a&em{65Jd0TX!mq$>tkikRigRqFTctXNq@bOun$r;eHs1 zU>W#3BcDl7BpbIL);<3GZ$$o|9;J`FP_XB}prxf`EV@iDDBK9d}>yT~; z35;BAIFy6p3F$#dXs_*&QLEiPRpp8yC(fci3nPYtVeAKYs@{$zHalB`8zTm6K22}m z+*M=5=iwoU`NF`U44!P_WwV;O>#j|A_UB8b%I{*-;9?~`|-nRpvESR%z16>MAMwrhLSL0!A87 zn0?H2+JmM->jbY3SiYY2W2bhilf5kzyK>jW)8lG#ZhZeT5YF*0M`kyaFn3aj(#xMI zcuJ$yRQ2>|*{tbCaf`xZ6F1JnFDj{bK@~xsN3}cyS{sctAT4f1G=>OZC-5W3EiFWy zCakL-NV|vLX)DkXlueaKnc?Y04zh$ZUiS`y)(c{kSdpIu(Zo08_K=smy}i9)RVugs z$v_&j4-zYufRIKd{C^@(uKx;zJO_*D&qHCieZ0zPs zbE3ler{?O_Z^yg)DL|=5Ydewu!0jTK z*a+Y@1hJ|FGvS^d9y;`BF~^To*Rc80dNI3v5ou{-)Msi?@d35`YQ3$Mfx83zk;4ca zE+c@?CgTaMCzm&TYqak;KS&~rr&nesjOdd7C5@!mS;f@-xY^Ht@23QJ3_BS?s-vX) zbN9dIE@_HV25#v}LX-;5%TtngdP4B}QF&_0(I0HCEiE0zJ@TG2Asv4wj6e3!e=ni5T8PHs?pbv|CEpNv&=uvMFJ@f$EFz2;SVdKQsT--6 zzG>@QVqXZYMOr+iUrj@XT+&mjbrSQ0__gQ=d1f<9koD42p{EI)Dnl1s_7*gpnQbtL?v|IcYgjx<-ho7YhIVS>s8T zfJ%{(TOQU>cn3Tk%OSgn)8xYovnDnUyRAQ4_xO}cQHz3fj@dt6oOfqv2p56 zq7A`CLsVuizZ=6?jvp|No)Fd7dd@_Rsxm4W8^TT3 zc`2`}f66E7VLY+ZcRshBRXyXc?L=AJa7O;jyx5bd}ki_ z(W6g-4e8!C+{&Be=dA4;eB5V-PJsF=c)UjBD9i{w_k(f8jq|vD{Y(U zmkEj0Q1`X2g0Lkv1{?$}+<)aK{uR)?lPUD}V9wBZb-(Mtp&^9~$v;@jC-aJkB==k~ z+UQmt*i^3C!j1q3AQE$@_OoCUcN&(9warsLIxwpQ73(E=b*OJO3%y{_T08|Miw+w0 zh|m0XUk7m?Bt1{|Rolufj_#f0ClGW`+#WWA4{YrZDKK$RA)TXjrd6zB)fu|AL zpS1WyHLnqi=h9s7)}G5X%_C%h2Rm*#{S}mN%^?W!#WyAC8$Bz8x%VD$L88I1j~-9_ zZGLEpWl$3y$u)WF{FRgFKAeg+yX6CYW~i4I((x>0ASGz4_g~kUhc$)OG^~%*IKQe_ zNUf`2)A<|~s1+FrCp>;dg7K;sO_0cm{ZN2aaS=*WqC;Cju{=d(r!5_z!co>EuUH+YcRceP$aMw zgOx+CkeNZJdv{&Q!Ko zy}h?A!9aPPb3tZ7|Hqs_J^7;M|Vb;SDY-TQ5A~RP^mxfPt&%_diz=WUH!lLddmEyUzUvo@i*C5mWj; z5~(eib&r}45m~-cruG&20me^p@#mdjQ;Po-p(tNr%nHF=6vb35wd$Cp?7_du>3apF z6L)dDSEjTFQ>aIcX$jXJk+U!P`YNL`W63|uJ3jt2q?^GHR!Ud~5zHdW%yHC8SbCja zJ=8rH0J?B6;O_pKq@JdFw(&ci3OP>v^zStROyQfQRsJ8e)=9R0u7pKv<0u-7uAUf{ z5PC|xigP&o*25sv=iT4z@UK%JCiD9WD7SYXA3H}-2}^Lb9XhsCHZ12#ORgP4e?p;c z^ISTE2dmh{W^y-3x&L;Enq`0Hr*equOyzSeuwOT!`}PLQC5+|1HjYq36bkcjke6md@EpA-|{sgh575xeplbLIjP)}t)!%l?{mogS;sfXZz|pm zw#MA$ruq-h`2XnH=?mNzEjLAXwX}5Gy4zuw+OW%!Oqw`+S{t~N>ZlN))xjG-8RCX9 zm-`LZyz`m8@JB_vF1rwouyp?;J4-Op{CVQGuE9$U7%Q6>37WJzca>lxxzhUoq}=9f zcc4vGR)u%Z$Y*0=9K(6+P`XDk%c)`mn+1l9aIIM+;NcT@q+e>hRwIJ-(9!}dtP{?X zsV_(EbH1M|#J#vJt9qu}aCm6eTHa7&F!1kH9{r3_dF zn(=bz{Sr43!+-QX-RxxW`bAvg1X2#cr8~y-X8cH=Ob9zXTFx^)^9OQdMbRAHoGQ`$ zjBEMAM~8vkpqBO&{)6TlXO6SVcA=fe*(c4;1@#p2j9gX-4xTEqonU;YvHxO4N2^A+ z$f7dx#bd{1o3yh!T|GAvb?e6QSvJZA;j zVYwz$P6^ZdmH*j9&E0l-8-?6wZ7GbUMjE+QA~)1>{Zz~o3kyHAHGSa|-+6n;pz~~{ z<}{cn?Rpexkvg79oIYpeM%wT10x~*yR*Fi-m(NH_CBzNBq^gHcYKW%M=Nu+FPiM78 zIjrrplNjBUF0QO>FU8is>ddN&wihqGU`@%pc9w#pqax<^YEql%2PI1+T@0a&Al&yc z`*2uf1H&l=2Kz-)`<>I#$N>+HkC|){kv44Xis`3BUnRS-_NwEdY^PXo)lj4qa1uw= zE)IlY=!Ba@mAexokGy4+Jag=0SfuS!|LPY$vH76|?+C4t4zDfrXMP*gL>b_x%JV`c zr`3vZ(O;Q6ux&MZHVrc3A_Fgumfm^SV1ey^hiP(XUn++Z6pm674RINRyUTB-2N(Z7<~Pnm~5`8Ahpk*&00Cz~~q z)k<~hBxRTB)9tVj$wIMSMK&&C!5i^jMf`%KPnVA2qs zasQ?wvce=Ev9IkFiz#&366nJ>74)?AwLPuOt%FWoB%_ zU#m=wZMAvS@x;)=gL*f5n;s|qtG&{DZo&rS$wDjkoGPzx@C}KYFH5FX7dZ#8g>+Pz zIah4>x(mPgQjxOtFppHkjGT@mWaK>ADIvWzb=BI!w9Jp(Ei^VbHGDP7-N^abFl7o0 zRRlw}|A0R#jU%{s&93+s<8#cw>`{!K|IiI)~0n!1ZW zg*?hoQ&z(Vu4Bl~sN-+uiY_5DG{++!+iv;>jyZ@Tgrq(G|8gl#bScJ8RU;seJxx?ArxLxY0bzqLf;NzSY#?kG5a zaXzG;Z+ymTz~lx~R*!lgbG2oj)ANN?c{5IWnusk@OygA3t&%VvlGol{jGaX^mVLIQ zF(G+k705at<>Ol2RQY@nkCM|YHx*m9%9xsxe6mdqF@8tJO!ZUek zL$f>=^A0<{4_rymun%Fw8RlAbc&d+v>1ag(wCLLoUoa4JM;db%da#1`Ldj3z53y8LOe}x zy^!ew$gtM$!4skG$2-j~WV_+}Pqjh?l^i`+ajAgo;5ydC4z z%)uFp6&Z5Sy_myZ5Z1XNt=H(0U8YUxCNS;A1US1|dnUg_-)FB@{(9sVjv5~U8V`~m ztvvtgY%uPV(@{!Dku5E^DW8AgqGe@SbeFPJ-5U`kx4Ay1j@2pKzSu6E(=o?*agzF{ z95n^&tI{{!>fb)flH&@eC^`nF1t}RNWXyHg%+7Q^k@GC-kE&Kt6ZOBtLb5r6Nopx0 zJ>!Qt6E8?W%A+=*)It$KI1#XEcYW?fT%mk=jm=|fT~c-uH$5dXI)tU?+L7(3bs{zXf@ENbr^R^;XH1}3a5TNd{MweBuCN2 zsQ(M&$(%Gj=421Mww@@Jn)A4MsLWa%2_}CrGxVx(Vg6nJFp1N_n4&~cbZSTq&QQ56 zVV%*<(8>lA?zYtGhD=*UYe)2a-dovwyuZJ2ghTqGOnBq&5~j&7ThT=P(Kv?imr)bq zMJjmJF-*(G)yEb*6YZ&tICXV()3O`nbsKM7-IZsc&fiGNpn%KXUoSulm8fPd$gN@B z+$5q;3R1}>cfR_JR&7oXa9pI44W~CIKNXLqbX-Z9C_m)ug$VB~BSJE{S#zaz(i-u! z%V8lHL$)Kf`!ovH~6D6hsX?iUtLqNkpg-;c!tZxqMc zY-h80Ad!|+)(rt!G{a<*D~#0- zakAb4?<$$_;yFgaWKdEokm{L#|0E>^+mav%&~X`$sxu*9#Tx(Q zCIO86QY1y$Cn~W=ac0ebd8KGb6TUw zW^U!{zKuVGEvGGf@bkfmt2vO0RKNNj6}D;5B`t(dh7;x>KIbzYVsyLjCm^m!6Wqrr z;d$a@bi=hr3)J3C(`Jfy=NOS3P!ycr1jy>&v7*Sv8<7|IPbibQ*XZB&*b_=-zwD}w zbeS5FA=|yem%5--f+NY0D7&y=+TZwxN<>AC+>3jGJjaEqeDSOa#t#sb!-njH$;u;rVBw zctLY=XDSsVqkF;E#ph!NQ*JkGs&wKCgt}sA65D3_F{kJ#HwG(@pLrQ!3Mt25XHrs* z?1I?QR4a`3#MqL1I}P64wU#aFne0AV;D7!hA9RtQgDBrY+L=akH6_)lkR$KHRo1Q> zDYzyKyx4g)@E5^3NePO4WpRxrag+cFI%-ual8AJv^7D6~DjrW-x;PP7-wHITy zJ(mSHN|K+{P(L)8x?Jo7vV^vvMH;)6Chd~mSGKsYfFPjP__CyA&=vlgL&#~JdEVdo zdHx5i#GTdn%b^!|-G8{xrp%1H>J+LQ++&*2oq3)9{sFa0C3S_DRC!<06ak=|t6X#4 zfm@~BlEkWdKfLEtHvA8|VJP#P&xu5Z>_m|)Epn=ibibV4kUV@cOOL4Umt85%P@UO6 zw|o-+%%^*jh9% zs2cia=_QB<@O|dQKTx7pS*v>Y*D~GbPN7(`l9i<|MrmbE-OsM8kFx72vc$!yYI2uD zy@i^){El-JvSiPOpF(a%V?lkXuMc!jgibbR0YNC}E&l?U(!~&Q=7-}T-Re;!(9zj> zDNw2Q$V0-@hk4{xb`t{Du~H7#7E(V9ZRBOe%9NerS>+)cRi68a$uN8UW}2|6U(T`?+7^Fgf?E4@I z)PfV^7xzoPwplCe7`!tBb^$p)7O^3+K8fhO5>->RLdxRF0uQX@5tqWKQ_ja1Bm`eV zVuVQ4UZib&zp~hq$uh*Lf_HNt)eY~ohlcrQ+l3{b1_#Z-mMfS?)im2S?}f7L0oY@6 z?OIN)H!q3NqTcb-dTfo3Gs_W(odz%5-ivEslULaT`A8XFw&kddxYG0bb(TjL)#UJ7 zzGb2(UvCjM2W}KSX-W{VwnQ=n|Q0UP3x80)>wmYkEa)%75`Q3U_c6X@) zgUEYBsYy33FniU1R(@`EwHf8`gSJqH?@SshS$0@Cv|b>SOIL$e6Z3ZVa{cJzkt@4=$=O?Xu9ZHS ziFlE&byk#S%?VH%&o1fx7V&w zc|jYsX3dOYe>Rg8iJ?7qgywe$py<>)@e)e1UmJ)?q^YNyLP@v3W8+kf`eS&6_K`PR z*>VMX`mB57l*)K$<4VZNMuKtM$4(31XEMbxQQVXZo0;Szz5~`d zT3a0v>MOI~G5Dz$epQR5Gm%EDZ66bHe0s_&U0I8)46jIV?khuK@rN<86*_(lXPfUp z`c^2A)mWGS2@CcTdm?sbK@e*2^L)d;iBsykNU59FLiN@pLXy^5(yG2c)jr(AikTNYVE9@i@%RmL zGb6+J@=PP7kRx?^$5NsHqJ}3^or)dWu=m0K^zC`%`OSag%hA2IC&JH_*_)SkuIX0m z;3l$=(S~}JTRt?u>LmMZzZoLm$aPx`gW&FVn`!BV1xJFaH))&9p6(13Y|2dT6dAMF5MK!ctksF0^8`IUS? zX=hdC6OeuC@<%sf0x6Yan9a?a7r~4jj{&S7>TpYIEuJTVVP291L!Kg*5mORSOdXwG zsKDcT^Gy&5BNFIFTLk3^jkz^g2A0{RAHXjlL zHqs1SR$7|-px7#T&nfWh+U?tQrG=vOSKhzm`++tCzNF^mpR@LsW%p?>BgI7@)^4z_ z+tT>q_s4%P+f#Slzs5s75ki7*q@_zK z?o!{kK2Pq8VGXb4($!~ef9f0>UppmC@@iS%ID2$y-V|~BpovJ|X9Zj-k}ovo?H%A0 zJKOL^Mdee-i_dZ{Hy+NOdpzyA*kxfEnk@KedpZ5ABb^0%E~A7hBbHr5|JPirhHRH& zDl`42n!P)R$*Ii(%pO$qh=z-pPI`7iQj=EAUBY2TdI|-$m&kRI5)sY|`m!%ve_4HH zAScd^oo?vduaa*(lr$X_6nHdnbLY-p^MZjQ4=qx6@LqVLjp5)(+iah&(!f}jMVx~- zCq8WLsDe?{*?b+5i?A_oazDE76@?j(uQv~Up!i_BfiCF&a#oM&L*+n!vWi=^K}|X? z#SJs&9YkBB&|^`BMHe6Ve8$RI=PdL!&$zXdZtxq5M}d=6g0WZb*qMIqpmD=r;!;1eJtW_@nR zAgH-Jbkmg?mz5*KsLDvD${EzPe^m}-OyM&l0VV1G&nLS7>(eb#0>)UG<>3X78s>kp zy*{?Pn?A)41Tk~LqgF0O-AO0**U65f>cBVF=a!c>e5>D^^)#RCqi7<1c!7kHo>Rri z$w_N#?Xp0~34?r_vAO_kU5$T@!XWXtLNIa7V%|4av8l~7Z+4!`$f$c{Bsx5tf{cW4 z_E;>Z3Il+1K!AH_bF8hkwP`(~x829Lm@w^2LYS@YS>yDJNtxB^>^(;4|AURu9lC

Yy5hUbFYieI{%?@Cj6ODjm!1`W{yvmoj_aDoa=*FzcjZXIV4cb6tVD zM30H;-#C1ceSrFme}nUX4$1!>CbRL}2yQh`&$Bc%G^KVe%+5pQf`eCnXM%q8p~IGU znPAJ#5uwStzwpXJ{KQEVaZ`|WUsUVL6;#gJ!lOkH8s84ZDbw?z&_ErUX^{K7anaN9{T5dDx4AVLApq>o$0$MYN=DM-&wT zml^YC5wF|k508i-$^@;+z@FW?(c^Fbd#Yu-J1#I3MJ|A{pDdgea|f|LLhaUQ+eS}$ z)4&^ZBCu=N0CR1%mDWJ&SqFN0ZiSZzHSng@Rtc+j@7+U1-&L%vxM3`OYYz5zw+n)H zgbj_2XYM=IRoq(~NE%!N;vX`k$h3^V`Q%y(&%t&xs)&}K#U^!+1dx$obbC=$Q0^aHR@|8K+Lx*RJ=hS8_};oi)(0K=29Xm6Lf z8o2dw0f^w!7VW1(7mGjYWCAiL{jz{ZEvt1k;OY8-vHSdVrFhVvwp+kH9&3C~8n`eo zg@=ci2JCvBGOyIGKg~#G19g!9&OjT)0kMQ2Mn*<*R+^iU%>*FKIKH#HTXuK8$3k5A zBup^{DjDDEaPQP_Vl&V4++0j={#8_o8${ndY6Zr8QQ(LV@dgYP35&4KmxeG}qU*t? zYqgc}z>%Xz8-n&{gC-YVmjJ@D+A6jLEzkJ;nNA^gU*f?aU_)H!v{a3U(0bIAmyy34 z*o;4ZDV7(4*Z;{f(Zo#dI@`oyFnM5FSL~1_&;gE4X0VNz$tCY)GvAZP!O3|+QLJlA0K--XhP1H(YEWwb+G+IFufH-#WcTDcO{^rm} zZ*eQ2WK#Gul?WzcTs56+?go{u$&=Q9Xb_vx8sG15#XyO*E=kuV5m?y`6B6vFuZN~s z2b$HNj=ww({1V(1)rMH`_hnwo*U^k8X!|pd#Jz<+3vs8u%L2?$v8c8RFZo9SRfeU! zZf?Gj4GI#^^YX|cSHc?aVN|TGdAyc~!*+IdEW{g4V7H8WBZ%$MNXe?geVfY10|<9inT!q_z2_1g98>GAbuPzZ!+luQ})g#)#O zXVGG61_mtjoU%Cil`7ItHjLJP7fd|C^;S=WA9Di(CS>`od}1@P2nKGBaj?urpn^$X z%(D->z#n|*n6wSN2US?p40~X=#^H&TkE<84*yAV}RQ=wx`1WkjzJo=;Mt`-{@OS<} z`>)t|N*;eb7iVX~&l9Ehhn)&uxz_!;S)4K@UsCfsiiVfq0Jp2QN-r^0yo4+OD>}Lf zRD)COhcsP^lm4}ol+*%jt_Ay0y6xO3zZh0FOfl1Xvczz4ak1Y`wn!)aoMtWLb^Yjx z^1-%Nnu9gvgXh5BVDhcA)5HNpY?XkiyxQv8!1mA8zycCI@0A~Ei+OgBE(<(E-RdnY z_W5-^y&RpV(pPJjD@XcFPbGs7Kf!90H!%F)yEAnwMgR>wLO2{y)vN!XPdBtneb;V6 zsbpnkC1dX|tlA?L4xJYjxV>vvX~k=Bl=fGiH^c@)o=gL>~V2%hB)Q@ZIPlFUMhYTG>V9!2HSu|kR#Ii z@?aU2sQBmwq;x5kO@CASUz-quDY?H+7Q%}y{R_U`8$kt_2!48Ue|I(FnTQ@Mx!>ke z3l&BgyynZ7FYy4-evo$7;TgcP;ko5t%?AQh@OgFh5ztuK-h5x5$z1uU@kZ63nWvIY zJp!zU52J-nSeq4ubQIJLrQ2u`5s@xeKM0%LKr=q)?ORrRdwWxm;RQ+F{DWN$kX>~d zsYJoF`&qm&;`D0OR){KNC>82vU`mgp@M-1kPYehGOXPcngKks6^KKDtpK*-yk1D3x z+S+a4f=1W=_;y+Rof5Pb(% zprb&%&?fbl7#H-#$i(uAhdFF=sJ;x!=>awaX8N_Qt@og|tc(fRkhY?;^to{zJmBYw2jFRT*-!(;+oJTl3fCqzZH+y`&nzTFM9P6-JKNh1(vA&9)!k_v24 zF-Oi`i)D#!)U&j7+4ZKNpt#mhkGO%angkH{`M^JI)@YXmnhYA@y(Xrnk&fp?|L(yH z+WjX*YvlAGM>O0@Ujt(8r=DDmW15+%H8rNjiaZY;AjS&|6iH0zUK>&%2H*?4oT(3eW zp0J5huQbIRsRC#93t~yOW75uD7mr8KZSUJt z8SLt4?RIbWKB`4=1Bk?S0!IDUu1?L)Vj%h>DJum6J3(RQ{Y6G6GiKriOT zbnBB;Ni0Qi)g+CD+b8lWi0#&!_XdOh9)_#+3=OA7hN(k%ee?~Dj67>$IwA>~=69;e zA@+z<_#L?z5XUJO0@is0;$aBG7etC;!4I%`An+ys@yd}q>frxPWbXsecN=bM!W;yJ v3%Y;*hMP4VnSzyrW+gFa(KYn{`DtG{^|cEjCo@3=fj>nQ=Md603L+BHFf;~|(k1=S4N?+ANJ>lh&^>hL zaQ1*c@B5zbJAa(tcaDGX6YhQQd+oK?y4G6PT0V+UNxU26H?Xj<@T8^0l(DcbYhYnr zvcb6w{_@Uias~^FBVJnUxr)QcYMh^X=wST*f_t>)!)%o6F%xiPi4xrovY>u>jpo5jN-mM@8SXEjrjmp4;hzJ4`WmfFMPuZe^;D3qM~ zsC_lQXg`()S?ZUmPQ|P>6yo*dV6PQ z>z9Ps*t;Yo_B)H+p$rOYDk@UR$GlzyR1tWWbv(ku!*40c%HA}%iUs+2#`*AJXh;Zc zud9fN$n)pVU%g7iQCOpY@IYQcfldSt&O-Z_a7s>&-272nB*(9gj!csBL~?OoZ$ijr zYLN{a?uQQz^Z3a9*ZimJyuH1XB?DXi4SoiEE(?y1&idZc))tD0i;FYrjJ=MBSLM8) zX95?96F%JuB?^1>>Q#4lcX@d^9UWbsy4tUhh>(z5{t|RG)2RHK2vV!F)i*_l(~W+? z!NFTMCS11VEiCdMKW1k)FF5()j0H)+@w}H`M<~e3W@l%=xO%!bYLR-{m-!+!HMOuKm`?IF z_NUrObk4znLqTUai+X2*Fh@z6*+T(=EZ(nCQNdn)EGHZF1jEC_%b{(Jvotp%CB;Ye z%d*w;TUii;h4@~g<1xkF9ke%2#i_VJ7qnDmAdp*oXe2EUd`lF%eDkxj%@I{qRX3lSzJ1$_jJtmQ`W8Y+uU3{gJSGO^WToyuG%}K9 z^3;*XF|Xr|%@gLH5I!aggIe|MtgI0AQ1?joI^vTW;6LvNJ;O}jy?e)9Bq=W5qQcG- ztD?2W!boyEQSNjtAee!kp5CeQ;PiOkPI_u)4as$mo~M1Lj|Nz_uMWlbGw$ z($e}8+ZyACCMG!0g7zyzu(PAVrluyzzmCCd(C1*E*LutIoQ8Tr%DMHb4k@oiX zDjWEy=p?wYiF?>1bhYO53(Mf^@x^U59N}I>Xu^7IMSK$0tc)fA{(?^uj}t z%uzjNg6D32efFWDp|ZhV4o7>dqd@xY2@Pjw#Mo1-X^%U|Q`>u+FN}>dkBAgaY8MyZ z?k{8|m$A^+&0Sa_YN|czAvZiSG7>R>i}U2k6aA;TYABO0U%u=dl$kz_ta*l39npNO>N$(xS&8u zD^C|eU)ta5`{Bc--LB?l@5*S@O#1};Ki_7$ar@S-KJ}crg#}VlQX@k{pNf<)V^h=Z zja0c!lc`#LT3dm*En8y+1qEK|lRR?@~DOw4HAj8HmIE`y)jHJD5$N&7~?aW zQeh_yw|(`o>lBU4Cts_G8{6pi`gIH>W;(vq}{%nbs9U#+dtU%x79vk!rgI!QtMh}>UEN$K?T zls6{#@4qq1Dv6AUhzJiaw_sxLd;18Ov|4t_w{M($eDDGZ5OWhCui)@~fEKzXabcdz z_HJ%&c6KGl5!l$+TDfY-d-v`IeEFi21Hr=jVI;N)b#yFOv0JuRQ&V$sanaDw(AQ6D zR@Tze*6!`;=}T6A>gnMzw7PeNV#yi{YvTz7A0NM|q2WtlU^b$oqeE9$SHA=SU%w~p z=3;NJmfjop)K%*VGrZ9e`1(y=r^>?aw>j8YSQC=-HZ{akLqlKHvR4-uH8zwWDQ~Ua zHW}u%zgcP4XobvKJfQkqO=r)2|E@43vEM{qMkXsiA71h}uCu90f&fB}d$l}KOOF19 zib~j+#yHQb5#Kgq+RsjO5s;b%}swyh9lDnni80^Ec3ewUf z8Z;FZgCZl-!}x3d{Y^7d)0>Apxw*N&^jy9HrIL`4M8w994GyOBMo?I*vysKd#_seq zd%yO&q|YE)w}}Vgru~2s@-MMw@4t=$hFS;-Y3Ob(3Syo<_xK|9_qul+?)7VbmfgSm zXI+rx|9Ih&pa?ZaJz&d}_4h+9=UawUWB>W~5&oZVFA4tNEEMzjaj_xDa_BDh6T|r; zb-o>hhLiSuw|`r7OW@Bcsl{g#mkL5@Xz$zpjgg2Da!a5K{%qy7H_mm+MT`()9v=8R zl!St0v?3zce(#uw?RpK|oOe*yhjZb@Fxx-`G5hQg_^#ou3l48WEoUayb4MLCBCU;I znpsC^481zv|45kb_Fk3G>a=iR0v-%KmRA}?Ac%|Ezpft%`7&6_uEY-}R?7eAY-Dv%caeI$(#YJT>b7h@MSPI&caMxjM?|ENrz+!Kzg{r9Dr>Q` z#q!-_xzcZgN;q?2k+JCd`QZkCOisR(lWRXXC>tACvQ1Y=P+Ki3DoQF{T3SlwyYRAO zs_&G!+w3BTnOS+*Snv4=qkf&Lc^DWNXytxSP8LPR!4`aGBgPl!NaM;f_RVW^rw9)H zz6fW30H)+?DPN$K3zW6HyPHSysg}AwvEm(fR8mq>1UT|CZxI}VMicrmew9DJhCCX4 zFS>^hi;|M+s=AkKWj<)UctIDy%E;K<)~2nN46h{md{&&V=PLKr-Z4nwM%3@Ou7DJW zrY5aq_s}s0z;R1Ui;3Oj>}<-2!;-(=xAw=RLU3CPedh@SQ{`$sv5bJ=fsz=quy*=$ zi#iP*l{bK;z~5c(f5)#w8jO6Q%s=fUHTa~2u0nj6>qsYmc$%91&Fyom`#@Y#&mZa9 ze=v}ax12p_J`i||=(j50->iZ+FCVG(OumRqQ=+iqiowE|9&TBe-})0J(o6?ucTu;- zZ$y3`O=co1&pf?3zllU@Z752VGUI)k9!XKY&-doG{xx9Qmq?FR(+d2(>`Yw6-}nUY z61USg-~-#oX16VBoYL|vV=F8#unED6lJ}-;CWKtoS)3Lw_`4_3;XjtfAfq((Rby@+?>z0Q}@^Z;v+$l(JTM+Mn&Jp zymby_%uwgijo*L#aObGws2jKnE4vG-*K;LXB4Yeu4547q`H4G7BzMa7RLR%8&AaUs zf;o$>7hnTd0xHTwR`Q%AYB?o@<(Ia9+le3c82W$N0E zFL`M&KzzrNuCUNiUv`lrl2@+$(T`NTe{EIV{#g7&-Y#+fyEjK4!4Uvlzxv_d9PzXb zDxWQ&OzxT&d46jkj(ozwHa{0t#%ug&=kUmP_gYM=iR2OyP~XC9E9ete+lt{kJi~M$9bA6vuDF^hrO2(#qy63 zk9>q=TNDlSjLZvX^I-|e3e~4 zS2gA^JKm49?~dF&&3ZIOGMs~oPx>99FzT+bFEVv)^zn0Q`x1H)t~R+K(Q;{9TU%ja zVW4VJH!w6b^rnC=Y#&yor!z1yF$oC?v9Orx>x;ZAe`)Z{UkW&ho|)o>ll1!+)niL3IMpt{9a*=D~w%+A)fj9gY$wym{wjPEw7$)jjb&#T$}S{fQ(V`2~*1}_7e+3Rf; zbz+$D2n=ll<%Wha5hcjdMbyK0cJ?(KGQY5(uC9)-;N;}|yWgaL6@-oE=H^CnY}N!wqHkHc4jU% z)`s2dr|uoJeI_N>@b3muaB*?L+zwO2Osx&V#TEljB5G@E4;vG{TGUFQAHDd>7O-ej z1tiNf-mTdsKGBrNoXeE0t|Nv|K~7GJzs$HdEt6EL&BVKF$(EOw*U7u3L*` zr@E?>5e!*^+Qm5j@xDogKelC*D_-V$ytlTt2Cn(jHG(YBfc2U)G_$J5f0Lty6sM&% znUs(a5zWrcF&TDRycTnMStKa_!Nm2mR6(KRN6M?`K{|WmbbmtjI^*0+SJ%___I6M~ zRUfY*RBt6}@UpS7v9ijI93LI2B##9lOjls7hF}p!%QKIG(a2u4PBYCC2^C(zA{ZvL)mao76 z>0)vqG)`c5v71h0-*|1gNf*V^!%i@=0=>wqrb51H<<8B`0p^R|7WR{SMy_=*8^Clu z6^Nq0k8+dq$rIPTm0>_mnU8*#Ujqb28SJr7mn)_2eqX>M3DD7b9Q~cmp?!yjhB7iT z=I7__Py;PUz^$x|mIij6UI$QfoTZeTKeRFTiw`CFdwX|-?fH%uFMieJEBkcm7nqHe z1$UhWj?n(z{{vedPHEvLp zRewCM1!7A`UICxJM&+_GY?R=GBPt>?`@NX3f{LzWGVzN!*tD@kI#mjjg}!>`w=wKj zp9qRWCT8ZWj11104v%{{W9HVrXdwAUY3he(j>FzmXjScZj zPZYFakzGx~@$&Kl#S5qh>&K{femp;6AW;z@;;7!+CElUL;7(2n&e=}7u(+fo3}K(! zhbKb7Xl8Hp#z%>*YQjMSjxhdGZ4iH%8jT68k7E*mIEY zlHKam^5cB;sQCN$h_2_sFs_q5^n&dA#C3DliidI&Tnd-Sj%sjXVq*BrhBgw;4igf} zo-cm*()Y18T|m{ijIY4ee^9*(zQTRsvPcA+;r@LghmF6R@I<^~c9SWpPc(Hf`OmPr zzv=g#_-cbl8qqpUQ;7L+;pbmCub#L6NAT~4x4(Gi>Bd+$J3AZ0XK_n_mHZ9D(eyyH zWalXp%p%IF=VHNwI$`u?AZjpwV$T$Za;#uD&N%W~|B)|SVg;hY>L|eRMst_D*ibp4 z2mZG4ptcT*MWo27iQZ^QaGr@;bb=%iY0hA)L4v<9c8CfDCM?4&xhQB)MlA{>eg&vE zK+z;rl^ylC;80e+{ErL)$bX4@r^VOO#YN{-kV?cW(mHo}>FAR@z=$v(mJ(nw7bg^i zFk9Q$AlC)%2ArI@h4?9u|{Wv zbFSxyMP)IK>rZx3PmdydT9t`3dRp5v>ypXINeVu5LIKvBZC_pxEo#A(kZ;n~hYXtu zU7mUP%E$0R1affLTz>=3kQ%&`wU z7wDPCm*)4P7SY2XO32iMfYF9AlY+J&nki}#iIfCIW&SOwWpdfZ>M0%wMi27-Yzykj z9``#@->6y_)ly$6>S#y4ioVc<<0r8)b1NK5wID>rroWfFCM;ObY@_<$6cxaMpIG0j z{~y@~zzTaZ?)$NW93moKJQ9OX9+0R=%JA^JzVWNy4KnVt&s{cFu142{qPAE(L2}d3 zc&?*m+Sx9~oG17TI=~SCFaSp=LLVO8;7X`S`izp~eXK4ZrgLtr*(@UX{&lD%cyk2< zI@o&LS<-gxhV2=@`ze>^`FV4aaMX!ITU7M);zL5F-J&Mq^H#6&!0apXOYT+cysh1ufnn*0pLMO!|%PqvQXQBq zfL+7C`#2~%xCs2a#Vpp8kup0SkG|adE^|n53YpdnVjp z!hrJNsXGoR0_zqcPyDlvm(7qRF)3Fj6PjJn(c&7~+6d!bdZxcXFphP50gy1UWq-5nslD?k?}4QPbq~0kxY8G5U1k>=X%~I6asNe4(t|+Scaf zsfP9w0y6IG?n?VU!^6_Fl79jentI@yQ(c1~=9fr&+poRyTZGw=a$l_$Ls}Xq#F_LUN z;V2tBoq_-bZjB{|;jv{aa6|76i<7pRleTrYeFSJw0VY0oGAikdl!NF|C-)2}HwI^S zGTSaA@MGPA$5ma{eIkER`%pG>DvAva4A@WAef}kU=7Ji?1-WlyVF*Xa#! zWp0@%#<`iEU%fQ;)S_pQM2&BmnM(&xlyBrXV5;3ktGa9pS0cQ&=Q??>d-Q`WnUq z;_KuEk zj^Rn>er8aF(9jV75>;f~xn*N^sfsOV83q}oXpxf;eyqgn{jmO{y$>${uAUQONbVr! zZ+XaAz-@s5yD$4{MR3>(ePmE&>^VeM@gC>=GYliQ4_D>Q;KG*6HiLw3FqVYb`c$o* z0tQ-gLinqDdnMCQ6%|Kw42;AO=&aMX+q~qn!d<*~f z`k;?tlOw`dF0{9Fsx4swChG9I4!VM?d(MO&u`6n9kk>u8dYA2l_X@+3oh!yF!xZOQ z!t%~z%d9+_%ED_%wZ$^vf@^$_MQI2%NdD;qZ z3J@AOElb)1JRX3;u#4eGWOS0&u*4?%;e3EdPt<0^M9VJ1%TKMQenBd)MA^sQL6GY^ z$gEz&ubtHQ_O;k$*9C(zD*D~6BEVsV_C~_E#(In#ZP??Q|EPSj z)S5AQ6mj{?!ZP{TF@4#r2U(EIz`hK)Do5QxU;RZu7`ztifq&y(tVkCUTG&jMd3LDY zTI;=XFHn{ti^X8Y$2}}8pCg;n@WE3c2!_Ajexl5d$BGJ5YWq@DV6%0y`%_uIQ@HR5 zkgkA>*=hhh8_O0l{m(_>tHoj5-$Z0n6L^nuyy^GHp!z7vC$EuOx32L_LdA@y=-vx| zciBjq_Mb0E(61a4*8@TzF&zcix9Tn8_o+9|NrlXtjcNDv!fq%pOJ&ldMmpSFxS7Fg z>^`lnt%Q8BoeKyt(9lv=R(5c3P*WRdXyBMwEq+Ny7eaZ}h1R~xI4#Hi5fiKi<9_O+ zWl#zK=C#e5P$jYoZ%gU0i@r#y_e~yhc~Q|%{U)F$1hutp&N=J+!iI$zSgTrHPZtDUGiTi@M2hUgX zLX|)^qz9F#Qrj(tH_Le9ataf@`nUFpJ{(?J4AOcK;ZF;_hRkSaXzUTLwHNQhR+eA) zjai3AhAO>a$2(|f$-vbWE4v>JNQ#Pz1{3)D_!GC2gO!yWXu6>5QL&TCn(ua6dOE4u zjssJ2x5hOrbfdyPs?!~o-jX$%oU46ySjpigsJMTopCEj+I#)5eFV*HAbSU&5bant? z4D!{JX@&OEmvrvYn~`*XX|2%RwzuOp*~siXWJor@-k%w&ls(#RQAdmC$xb}gW+f2k z@9%$)_R$?&YtfRga%mYEjO7Iv>c1=m)Ia~`^lfYEB@Tw<`G^=>vuph08|KAtlYWZD)=H{prr-)8Vze}YZxrS zBG+#;XR2 z-(3dbK%|+KN@>(DnEQ#i(_mA*00Bpy@DPzp*2AhHOLe}dy!?oj#gUHX%amuu>xTL;g+L0`u|Wau z28y84NrumWLY^+4tBxKjD6lC6Zdu&+2~eeQ6Riaf6Wd`cg%dh?=tOx>BB*9po&0`Y zmK4KT=?Cqk;}I)_vDrimKGPCivW(x9QFShF4;rp})l z8Rj}rS@_;O-`F{~`ORQG$PO=WP(w?p^dsrfGoa1$z83p+fJlXov^DyVnPBV406@2_ zDXc!WQcBP^ugIo=F$wyfq`a| za@g>zo-ZZrC{)i3RH7Gs%+pICKEsSna=jaqI+;U^b(F6>z_hshyqTpS+;I?M{G&RUO#A z;Jq80tUgv-Lb)R-nl>}Bhk_1AVWWIdrMyYU2zeG zl7&LmP&lBRS+c;w!eRYy#`mT|GIDZot?HL-LwpFj$Hrp8Mi9tzB2!V&9&FAhyLrLsEp#lRizQAq1y}fN?ZQb0^*w7#b;=Ym50-&LQ4p^9( z38%EEu?s3HqVKu{1EPNTRr8dc)lyF?s2_srS=x6c1MXUddxB{=>f59sbb*|v+I-*LLf8}vsg%Sg944DLfGnK3Kp$`lZYw1ijKij|;&iqpl=v)~Q>={hf}95fc%Yv98A!DtE8}{F3nnGi zjG4{iEw@9diHKbOxi}C4D74o5Q2cz)J^?R$%8*N}HUTHn&#BoZBqdW>gZVYlO-F4o z^hX>RHsC6VNUH08{CEb&@K6dXU0sS$0gZ=u@x4}9MMXu$HgA9d9pgvjQK&hKGQd4m zZha;?36GAB29t*^-jvG9KOC)fE&R{2V`D z0+6T{m2M{QEkSixS0^G?M$!zs3q~cl3Bk1A+^ho;iYE&(#WHqghDP)(+${t7Azg}!55~E(*+23!2 zvboRaHFB98yLz0qf1JGGN=39l2nYH#5+`Cs{ZCpTIEYcMMaR5B5*q16q#{|r;o8SF$G7C zR65!KUj(^kv5?DwVihZx)_WgLl(;Bad|wsW*Qa7{zgO*g4A}n^rTZ#MN&T4d{SGIilR=(>3PYSF?!|7mF1yK zvyIwoRAtPa3QrFkc5^IhPI(0cw9FXTXuH8=D9C(d5BnSApE#3~gX9j3DSr#0Wo2ar zqf~i&zyZp@ypZoVfMi$18Rtei;{;JCR6i1MDXs@oUfbe~0BBc7N>G4p3WwRA99RY7 z&6YN=EBhxKilFr6+WfPWAPLg&8gz2z}Pep8yP>3x**ypQzR~`3WX!Wh2=gNfRrp`#w(B z)u4CtduSDu7`&r9<~eFkwg8*5%mfBj*e_vYR7VDpsoDjn_a2>XDxGjUraIfCn!aCN zW!Oeym6E~8&#&pDc(T(-3Ui39FfJWUjgN_s-`j`{vsGro$HODT^~Ao0FD5DJ3`sP( zd-rY_oR)Y2V4LvKBAw&f;$lXa4;Tsr(^3+Q7Si`p*TH1n4;f!f-zDZ4Q1@k429}oJ zhNO7ZR8@Ui8sq>5xSz~WskT+XwkK<`jU$az)zq9Epp_47p($w;=;=$-0XUjS~k z0wS0jl#r6QpKp)SMc7MAlR_X|yu9lzbCA=x>=E1n9Jiv2?%*@Q} o{Uet_#{WOl>sVL=0cQlZG9+PGPP@TiC=W|o94eM4^77;V0t;eg+yDRo literal 0 HcmV?d00001 diff --git a/time_dependent_navier_stokes/doc/tooltip b/time_dependent_navier_stokes/doc/tooltip index ca55f62..ceecfc9 100644 --- a/time_dependent_navier_stokes/doc/tooltip +++ b/time_dependent_navier_stokes/doc/tooltip @@ -1 +1 @@ -Solving time-dependent incompressible Navier-Stokes problem with Grad-Div stablization using IMEX scheme. +Solving time-dependent incompressible Navier-Stokes problem in parallel with Grad-Div stablization using IMEX scheme. diff --git a/time_dependent_navier_stokes/time_dependent_navier_stokes.cc b/time_dependent_navier_stokes/time_dependent_navier_stokes.cc index 85f0976..215cc1f 100644 --- a/time_dependent_navier_stokes/time_dependent_navier_stokes.cc +++ b/time_dependent_navier_stokes/time_dependent_navier_stokes.cc @@ -1,9 +1,10 @@ #include #include #include +#include #include -#include #include +#include #include #include @@ -14,10 +15,15 @@ #include #include #include -#include +#include + +#include +#include +#include +#include +#include #include -#include #include #include #include @@ -38,9 +44,10 @@ #include #include #include -#include -#include +#include +#include +#include #include #include @@ -51,12 +58,14 @@ namespace fluid using namespace dealii; // @sect3{Create the triangulation} - - // The code to create triangulation is copied from Martin Kronbichler's code - // (https://github.com/kronbichler/adaflo/blob/master/tests/flow_past_cylinder.cc) + // The code to create triangulation is copied from + // [Martin Kronbichler's + // code](https://github.com/kronbichler/adaflo/blob/master/tests/flow_past_cylinder.cc) // with very few modifications. - // Helper function used in both 2d and 3d: - void create_triangulation_2d(Triangulation<2> &tria, bool compute_in_2d = true) + // + // @sect4{Helper function} + void create_triangulation_2d(Triangulation<2> &tria, + bool compute_in_2d = true) { SphericalManifold<2> boundary(Point<2>(0.5, 0.2)); Triangulation<2> left, middle, right, tmp, tmp2; @@ -75,11 +84,30 @@ namespace fluid // Create middle part first as a hyper shell. GridGenerator::hyper_shell(middle, Point<2>(0.5, 0.2), 0.05, 0.2, 4, true); - middle.set_manifold(0, boundary); + middle.reset_all_manifolds(); + for (Triangulation<2>::cell_iterator cell = middle.begin(); + cell != middle.end(); + ++cell) + for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) + { + bool is_inner_rim = true; + for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_face; ++v) + { + Point<2> &vertex = cell->face(f)->vertex(v); + if (std::abs(vertex.distance(Point<2>(0.5, 0.2)) - 0.05) > 1e-10) + { + is_inner_rim = false; + break; + } + } + if (is_inner_rim) + cell->face(f)->set_manifold_id(1); + } + middle.set_manifold(1, boundary); middle.refine_global(1); // Then move the vertices to the points where we want them to be to create a - // slightly asymmetric cube with a hole + // slightly asymmetric cube with a hole: for (Triangulation<2>::cell_iterator cell = middle.begin(); cell != middle.end(); ++cell) @@ -134,109 +162,158 @@ namespace fluid // Left domain is requred in 3d only. if (compute_in_2d) - { - GridGenerator::merge_triangulations(tmp2, right, tria); - } + { + GridGenerator::merge_triangulations(tmp2, right, tria); + } else - { - GridGenerator::merge_triangulations(left, tmp2, tmp); - GridGenerator::merge_triangulations(tmp, right, tria); - } + { + GridGenerator::merge_triangulations(left, tmp2, tmp); + GridGenerator::merge_triangulations(tmp, right, tria); + } } - // Create 2D triangulation: + // @sect4{2D flow around cylinder triangulation} void create_triangulation(Triangulation<2> &tria) { create_triangulation_2d(tria); - // Set the cylinder boundary to 1, the right boundary (outflow) to 2, the rest to 0. + // Set the left boundary (inflow) to 0, the right boundary (outflow) to 1, + // upper to 2, lower to 3 and the cylindrical surface to 4. for (Triangulation<2>::active_cell_iterator cell = tria.begin(); cell != tria.end(); ++cell) - { - for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) { - if (cell->face(f)->at_boundary()) + for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) { - if (std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12) - { - cell->face(f)->set_all_boundary_ids(2); - } - else if (Point<2>(0.5, 0.2).distance(cell->face(f)->center()) <= 0.05) + if (cell->face(f)->at_boundary()) { - cell->face(f)->set_all_manifold_ids(10); - cell->face(f)->set_all_boundary_ids(1); + if (std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(1); + } + else if (std::abs(cell->face(f)->center()[0] - 0.3) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(0); + } + else if (std::abs(cell->face(f)->center()[1] - 0.41) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(3); + } + else if (std::abs(cell->face(f)->center()[1]) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(2); + } + else + { + cell->face(f)->set_all_boundary_ids(4); + } } - else - { - cell->face(f)->set_all_boundary_ids(0); - } } } - } } - // Create 3D triangulation: + // @sect4{3D flow around cylinder triangulation} void create_triangulation(Triangulation<3> &tria) { Triangulation<2> tria_2d; create_triangulation_2d(tria_2d, false); GridGenerator::extrude_triangulation(tria_2d, 5, 0.41, tria); - // Set the cylinder boundary to 1, the right boundary (outflow) to 2, the rest to 0. + // Set the ids of the boundaries in x direction to 0 and 1; y direction to 2 and 3; + // z direction to 4 and 5; the cylindrical surface 6. for (Triangulation<3>::active_cell_iterator cell = tria.begin(); - cell != tria.end(); ++cell) - { - for (unsigned int f = 0; f::faces_per_cell; ++f) + cell != tria.end(); + ++cell) { - if (cell->face(f)->at_boundary()) - { - if (std::abs(cell->face(f)->center()[0]-2.5) < 1e-12) - { - cell->face(f)->set_all_boundary_ids(2); - } - else if (Point<3>(0.5, 0.2, cell->face(f)->center()[2]).distance - (cell->face(f)->center()) <= 0.05) + for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; ++f) { - cell->face(f)->set_all_manifold_ids(10); - cell->face(f)->set_all_boundary_ids(1); - } - else - { - cell->face(f)->set_all_boundary_ids(0); + if (cell->face(f)->at_boundary()) + { + if (std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(1); + } + else if (std::abs(cell->face(f)->center()[0]) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(0); + } + else if (std::abs(cell->face(f)->center()[1] - 0.41) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(3); + } + else if (std::abs(cell->face(f)->center()[1]) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(2); + } + else if (std::abs(cell->face(f)->center()[2] - 0.41) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(5); + } + else if (std::abs(cell->face(f)->center()[2]) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(4); + } + else + { + cell->face(f)->set_all_boundary_ids(6); + } + } } - } } - } } // @sect3{Time stepping} + // This class is pretty much self-explanatory. class Time { public: - Time(const double time_end, const double delta_t) - : timestep(0), time_current(0.0), time_end(time_end), delta_t(delta_t) + Time(const double time_end, + const double delta_t, + const double output_interval, + const double refinement_interval) + : timestep(0), + time_current(0.0), + time_end(time_end), + delta_t(delta_t), + output_interval(output_interval), + refinement_interval(refinement_interval) { } - virtual ~Time() {} double current() const { return time_current; } double end() const { return time_end; } double get_delta_t() const { return delta_t; } unsigned int get_timestep() const { return timestep; } - void increment() - { - time_current += delta_t; - ++timestep; - } + bool time_to_output() const; + bool time_to_refine() const; + void increment(); private: unsigned int timestep; double time_current; const double time_end; const double delta_t; + const double output_interval; + const double refinement_interval; }; - // @sect3{Boundary values} + bool Time::time_to_output() const + { + unsigned int delta = output_interval / delta_t; + return (timestep >= delta && timestep % delta == 0); + } + + bool Time::time_to_refine() const + { + unsigned int delta = refinement_interval / delta_t; + return (timestep >= delta && timestep % delta == 0); + } - // Dirichlet boundary conditions for the velocity inlet and walls + void Time::increment() + { + time_current += delta_t; + ++timestep; + } + + // @sect3{Boundary values} + // Dirichlet boundary conditions for the velocity inlet and walls. template class BoundaryValues : public Function { @@ -255,11 +332,21 @@ namespace fluid { Assert(component < this->n_components, ExcIndexRange(component, 0, this->n_components)); - if (component == 0 && std::abs(p[0] - 0.3) < 1e-10) + double left_boundary = (dim == 2 ? 0.3 : 0.0); + if (component == 0 && std::abs(p[0] - left_boundary) < 1e-10) { - double U = 1.5; - double y = p[1]; - return 4 * U * y * (0.41 - y) / (0.41 * 0.41); + // For a parabolic velocity profile, $U_\mathrm{avg} = 2/3 + // U_\mathrm{max}$ + // in 2D, and $U_\mathrm{avg} = 4/9 U_\mathrm{max}$ in 3D. + // If $\nu = 0.001$, $D = 0.1$, then $Re = 100 U_\mathrm{avg}$. + double Uavg = 1.0; + double Umax = (dim == 2 ? 3 * Uavg / 2 : 9 * Uavg / 4); + double value = 4 * Umax * p[1] * (0.41 - p[1]) / (0.41 * 0.41); + if (dim == 3) + { + value *= 4 * p[2] * (0.41 - p[2]) / (0.41 * 0.41); + } + return value; } return 0; } @@ -272,377 +359,429 @@ namespace fluid values(c) = BoundaryValues::value(p, c); } - // @sect3{Preconditioners} - - // The LHS of the system matrix is the same as Stokes equation for IMEX scheme. - // A block preconditioner as in step-22 is used here. - - // @sect4{Inner preconditioner} - - // Adapted from step-22, used to solve for ${\tilde{A}}^{-1}$ - template - struct InnerPreconditioner; - - template <> - struct InnerPreconditioner<2> - { - typedef SparseDirectUMFPACK type; - }; - - template <> - struct InnerPreconditioner<3> - { - typedef SparseILU type; - }; - - // @sect4{Inverse matrix} - - // This is used for ${\tilde{S}}^{-1}$ and ${\tilde{A}}^{-1}$, which are symmetric so we use CG - // solver inside - template - class InverseMatrix : public Subscriptor - { - public: - InverseMatrix(const MatrixType &m, - const PreconditionerType &preconditioner); - void vmult(Vector &dst, const Vector &src) const; - - private: - const SmartPointer matrix; - const SmartPointer preconditioner; - }; - - template - InverseMatrix::InverseMatrix( - const MatrixType &m, const PreconditionerType &preconditioner) - : matrix(&m), preconditioner(&preconditioner) - { - } - - template - void InverseMatrix::vmult( - Vector &dst, const Vector &src) const - { - SolverControl solver_control(src.size(), 1e-6 * src.l2_norm()); - SolverCG<> cg(solver_control); - dst = 0; - cg.solve(*matrix, dst, src, *preconditioner); - } - - // @sect4{Approximate Schur complement of mass matrix} - - // The Schur complement of mass matrix is written as $S_M = BM^{-1}B^T$ - // Similar to step-20, we use $B(diag(M))^{-1}B^T$ to approximate it. - class ApproximateMassSchur : public Subscriptor + // @sect3{Block preconditioner} + // + // The block Schur preconditioner can be written as the product of three + // matrices: + // $ + // P^{-1} = \begin{pmatrix} \tilde{A}^{-1} & 0\\ 0 & I\end{pmatrix} + // \begin{pmatrix} I & -B^T\\ 0 & I\end{pmatrix} + // \begin{pmatrix} I & 0\\ 0 & \tilde{S}^{-1}\end{pmatrix} + // $ + // $\tilde{A}$ is symmetric since the convection term is eliminated from the + // LHS. + // $\tilde{S}^{-1}$ is the inverse of the Schur complement of $\tilde{A}$, + // which consists of a reaction term, a diffusion term, a Grad-Div term + // and a convection term. + // In practice, the convection contribution is ignored, namely + // $\tilde{S}^{-1} = -(\nu + \gamma)M_p^{-1} - + // \frac{1}{\Delta{t}}{[B(diag(M_u))^{-1}B^T]}^{-1}$ + // where $M_p$ is the pressure mass, and + // ${[B(diag(M_u))^{-1}B^T]}$ is an approximation to the Schur complement of + // (velocity) mass matrix $BM_u^{-1}B^T$. + // + // Same as the tutorials, we define a vmult operation for the block + // preconditioner + // instead of write it as a matrix. It can be seen from the above definition, + // the result of the vmult operation of the block preconditioner can be + // obtained + // from the results of the vmult operations of $M_u^{-1}$, $M_p^{-1}$, + // $\tilde{A}^{-1}$, which can be transformed into solving three symmetric + // linear + // systems. + class BlockSchurPreconditioner : public Subscriptor { public: - ApproximateMassSchur(const BlockSparseMatrix &M); - void vmult(Vector &dst, const Vector &src) const; - - private: - const SmartPointer> mass_matrix; - mutable Vector tmp1, tmp2; - }; - - ApproximateMassSchur::ApproximateMassSchur( - const BlockSparseMatrix &M) - : mass_matrix(&M), tmp1(M.block(0, 0).m()), tmp2(M.block(0, 0).m()) - { - } - - void ApproximateMassSchur::vmult(Vector &dst, - const Vector &src) const - { - mass_matrix->block(0, 1).vmult(tmp1, src); - mass_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1); - mass_matrix->block(1, 0).vmult(dst, tmp2); - } - - // @sect4{The inverse matrix of the system Schur complement} + BlockSchurPreconditioner( + TimerOutput &timer, + double gamma, + double viscosity, + double dt, + const std::vector &owned_partitioning, + const PETScWrappers::MPI::BlockSparseMatrix &system, + const PETScWrappers::MPI::BlockSparseMatrix &mass, + PETScWrappers::MPI::BlockSparseMatrix &schur); + + void vmult(PETScWrappers::MPI::BlockVector &dst, + const PETScWrappers::MPI::BlockVector &src) const; - // The inverse of the total Schur complement is the sum of the inverse of - // diffusion, Grad-Div term, and mass Schur complements. Note that the first - // two components add up to $\Delta{t}(\nu + \gamma)M_p^{-1}$ as introduced in step-57, - // in which the additional $\Delta{t}$ comes from the time discretization, - // and the last component is obtained by wrapping a InverseMatrix<\code> - // around ApproximateMassSchur<\code>. - template - class SchurComplementInverse : public Subscriptor - { - public: - SchurComplementInverse( - double gamma, double viscosity, double dt, - const InverseMatrix &Sm_inv, - const InverseMatrix, PreconditionerMp> &Mp_inv); - void vmult(Vector &dst, const Vector &src) const; private: + TimerOutput &timer; const double gamma; const double viscosity; const double dt; - const SmartPointer> Sm_inverse; - const SmartPointer, - PreconditionerMp>> Mp_inverse; - }; - - template - SchurComplementInverse::SchurComplementInverse( - double gamma, double viscosity, double dt, - const InverseMatrix &Sm_inv, - const InverseMatrix, PreconditionerMp> &Mp_inv) : - gamma(gamma), viscosity(viscosity), dt(dt), Sm_inverse(&Sm_inv), Mp_inverse(&Mp_inv) - { - } - template - void SchurComplementInverse::vmult( - Vector &dst, const Vector &src) const - { - Vector tmp(src.size()); - Sm_inverse->vmult(dst, src); - Mp_inverse->vmult(tmp, src); - tmp *= (viscosity + gamma) * dt; - dst += tmp; - } - - // @sect4{The block Schur preconditioner} - - // The block Schur preconditioner has the same form as in step-22, which is written as - // $P^{-1} = [\tilde{A}}^{-1}, 0; {\tilde{S}}^{-1}B{\tilde{A}}^{-1}, -{\tilde{S}}^{-1}]$ - // Note that ${\tilde{A}}^{-1}$ has contributions from the diffusion, Grad-Div and mass terms. - // This class has three template arguments: PreconditionerA is needed for ${\tilde{A}}^{-1}$, - // PreconditionerSm and PreconditionerMp are used in the inverse of the Schur complement - // of $\tilde{A}$, namely ${\tilde{S}}^{-1}$. - template - class BlockSchurPreconditioner : public Subscriptor - { - public: - BlockSchurPreconditioner( - const BlockSparseMatrix &system_m, - const InverseMatrix, PreconditionerA> &A_inv, - const SchurComplementInverse &S_inv); - void vmult(BlockVector &dst, const BlockVector &src) const; - - private: - const SmartPointer> system_matrix; - const SmartPointer< - const InverseMatrix, PreconditionerA>> A_inverse; - const SmartPointer< - const SchurComplementInverse> S_inverse; - mutable Vector tmp; + const SmartPointer + system_matrix; + const SmartPointer mass_matrix; + // As discussed, ${[B(diag(M_u))^{-1}B^T]}$ and its inverse + // need to be computed. + // We can either explicitly compute it out as a matrix, or define + // it as a class with a vmult operation. + // The second approach saves some computation to construct the matrix, + // but leads to slow convergence in CG solver because it is impossible + // to apply a preconditioner. We go with the first route. + const SmartPointer mass_schur; }; - template - BlockSchurPreconditioner:: - BlockSchurPreconditioner( - const BlockSparseMatrix &system_m, - const InverseMatrix, PreconditionerA> &A_inv, - const SchurComplementInverse &S_inv) - : system_matrix(&system_m), A_inverse(&A_inv), S_inverse(&S_inv), - tmp(system_matrix->block(1, 1).m()) + // @sect4{BlockSchurPreconditioner::BlockSchurPreconditioner} + // + // Input parameters and system matrix, mass matrix as well as the mass schur + // matrix are needed in the preconditioner. In addition, we pass the + // partitioning information into this class because we need to create some + // temporary block vectors inside. + BlockSchurPreconditioner::BlockSchurPreconditioner( + TimerOutput &timer, + double gamma, + double viscosity, + double dt, + const std::vector &owned_partitioning, + const PETScWrappers::MPI::BlockSparseMatrix &system, + const PETScWrappers::MPI::BlockSparseMatrix &mass, + PETScWrappers::MPI::BlockSparseMatrix &schur) + : timer(timer), + gamma(gamma), + viscosity(viscosity), + dt(dt), + system_matrix(&system), + mass_matrix(&mass), + mass_schur(&schur) { + TimerOutput::Scope timer_section(timer, "CG for Sm"); + // The schur complemete of mass matrix is actually being computed here. + PETScWrappers::MPI::BlockVector tmp1, tmp2; + tmp1.reinit(owned_partitioning, mass_matrix->get_mpi_communicator()); + tmp2.reinit(owned_partitioning, mass_matrix->get_mpi_communicator()); + tmp1 = 1; + tmp2 = 0; + // Jacobi preconditioner of matrix A is by definition ${diag(A)}^{-1}$, + // this is exactly what we want to compute. + PETScWrappers::PreconditionJacobi jacobi(mass_matrix->block(0, 0)); + jacobi.vmult(tmp2.block(0), tmp1.block(0)); + system_matrix->block(1, 0).mmult( + mass_schur->block(1, 1), system_matrix->block(0, 1), tmp2.block(0)); } - template - void BlockSchurPreconditioner::vmult( - BlockVector &dst, const BlockVector &src) const + // @sect4{BlockSchurPreconditioner::vmult} + // + // The vmult operation strictly follows the definition of + // BlockSchurPreconditioner + // introduced above. Conceptually it computes $u = P^{-1}v$. + void BlockSchurPreconditioner::vmult( + PETScWrappers::MPI::BlockVector &dst, + const PETScWrappers::MPI::BlockVector &src) const { - A_inverse->vmult(dst.block(0), src.block(0)); - system_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1)); - tmp *= -1; - S_inverse->vmult(dst.block(1), tmp); + // Temporary vectors + PETScWrappers::MPI::Vector utmp(src.block(0)); + PETScWrappers::MPI::Vector tmp(src.block(1)); + tmp = 0; + // This block computes $u_1 = \tilde{S}^{-1} v_1$, + // where CG solvers are used for $M_p^{-1}$ and $S_m^{-1}$. + { + TimerOutput::Scope timer_section(timer, "CG for Mp"); + SolverControl mp_control(src.block(1).size(), + 1e-6 * src.block(1).l2_norm()); + PETScWrappers::SolverCG cg_mp(mp_control, + mass_schur->get_mpi_communicator()); + // $-(\nu + \gamma)M_p^{-1}v_1$ + PETScWrappers::PreconditionBlockJacobi Mp_preconditioner; + Mp_preconditioner.initialize(mass_matrix->block(1, 1)); + cg_mp.solve( + mass_matrix->block(1, 1), tmp, src.block(1), Mp_preconditioner); + tmp *= -(viscosity + gamma); + } + // $-\frac{1}{dt}S_m^{-1}v_1$ + { + TimerOutput::Scope timer_section(timer, "CG for Sm"); + SolverControl sm_control(src.block(1).size(), + 1e-6 * src.block(1).l2_norm()); + PETScWrappers::SolverCG cg_sm(sm_control, + mass_schur->get_mpi_communicator()); + // PreconditionBlockJacobi works find on Sm if we do not refine the mesh. + // Because after refine_mesh is called, zero entries will be created on + // the diagonal (not sure why), which prevents PreconditionBlockJacobi + // from being used. + PETScWrappers::PreconditionNone Sm_preconditioner; + Sm_preconditioner.initialize(mass_schur->block(1, 1)); + cg_sm.solve( + mass_schur->block(1, 1), dst.block(1), src.block(1), Sm_preconditioner); + dst.block(1) *= -1 / dt; + } + // Adding up these two, we get $\tilde{S}^{-1}v_1$. + dst.block(1) += tmp; + // Compute $v_0 - B^T\tilde{S}^{-1}v_1$ based on $u_1$. + system_matrix->block(0, 1).vmult(utmp, dst.block(1)); + utmp *= -1.0; + utmp += src.block(0); + // Finally, compute the product of $\tilde{A}^{-1}$ and utmp + // using another CG solver. + { + TimerOutput::Scope timer_section(timer, "CG for A"); + SolverControl a_control(src.block(0).size(), + 1e-6 * src.block(0).l2_norm()); + PETScWrappers::SolverCG cg_a(a_control, + mass_schur->get_mpi_communicator()); + // We do not use any preconditioner for this block, which is of course + // slow, + // only because the performance of the only two preconditioners available + // PreconditionBlockJacobi and PreconditionBoomerAMG are even worse than + // none. + PETScWrappers::PreconditionNone A_preconditioner; + A_preconditioner.initialize(system_matrix->block(0, 0)); + cg_a.solve( + system_matrix->block(0, 0), dst.block(0), utmp, A_preconditioner); + } } - // @sect3{The time-dependent Navier-Stokes class template} + // @sect3{The incompressible Navier-Stokes solver} + // + // Parallel incompressible Navier Stokes equation solver using + // implicit-explicit time scheme. + // This program is built upon dealii tutorials step-57, step-40, step-22, + // and step-20. + // The system equation is written in the incremental form, and we treat + // the convection term explicitly. Therefore the system equation is linear + // and symmetric, which does not need to be solved with Newton's iteration. + // The system is further stablized and preconditioned with Grad-Div method, + // where GMRES solver is used as the outer solver. template - class NavierStokes + class InsIMEX { public: - NavierStokes(const unsigned int degree); + InsIMEX(parallel::distributed::Triangulation &); void run(); + ~InsIMEX() { timer.print_summary(); } private: - void setup(); - void assemble(bool assemble_lhs); - - std::pair solve_linear_system(bool update_preconditioner); - void output_results(const unsigned int index) const; - void process_solution(std::ofstream& out) const; - const ConstraintMatrix &get_constraints() const; - + void setup_dofs(); + void make_constraints(); + void initialize_system(); + void assemble(bool use_nonzero_constraints, bool assemble_system); + std::pair solve(bool use_nonzero_constraints, + bool assemble_system); + void refine_mesh(const unsigned int, const unsigned int); + void output_results(const unsigned int) const; double viscosity; double gamma; const unsigned int degree; std::vector dofs_per_block; - Triangulation triangulation; + parallel::distributed::Triangulation &triangulation; FESystem fe; DoFHandler dof_handler; - QGauss quadrature_formula; - QGauss face_quadrature_formula; + QGauss volume_quad_formula; + QGauss face_quad_formula; ConstraintMatrix zero_constraints; ConstraintMatrix nonzero_constraints; BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix system_matrix; - // We need both velocity mass and pressure mass, so we use a block sparse matrix to store it. - BlockSparseMatrix mass_matrix; - - BlockVector solution; - BlockVector solution_increment; - BlockVector system_rhs; + // System matrix to be solved + PETScWrappers::MPI::BlockSparseMatrix system_matrix; + // Mass matrix is a block matrix which includes both velocity + // mass matrix and pressure mass matrix. + PETScWrappers::MPI::BlockSparseMatrix mass_matrix; + // The schur complement of mass matrix is not a block matrix. + // However, because we want to reuse the partition we created + // for the system matrix, it is defined as a block matrix + // where only one block is actually used. + PETScWrappers::MPI::BlockSparseMatrix mass_schur; + // The latest known solution. + PETScWrappers::MPI::BlockVector present_solution; + // The increment at a certain time step. + PETScWrappers::MPI::BlockVector solution_increment; + // System RHS + PETScWrappers::MPI::BlockVector system_rhs; + + MPI_Comm mpi_communicator; + + ConditionalOStream pcout; + + // The IndexSets of owned velocity and pressure respectively. + std::vector owned_partitioning; + + // The IndexSets of relevant velocity and pressure respectively. + std::vector relevant_partitioning; + + // The IndexSet of all relevant dofs. + IndexSet locally_relevant_dofs; + + // The BlockSchurPreconditioner for the entire system. + std::shared_ptr preconditioner; Time time; mutable TimerOutput timer; - - // We use shared pointers for all the preconditioning-related stuff - std::shared_ptr approximate_Sm; - std::shared_ptr preconditioner_Sm; - std::shared_ptr> Sm_inverse; - - std::shared_ptr> preconditioner_Mp; - std::shared_ptr, SparseILU>> Mp_inverse; - - std::shared_ptr>> S_inverse; - - std::shared_ptr::type> preconditioner_A; - std::shared_ptr, - typename InnerPreconditioner::type>> A_inverse; - - std::shared_ptr::type, PreconditionIdentity, SparseILU>> preconditioner; }; - // @sect4{NavierStokes::NavierStokes} + // @sect4{InsIMEX::InsIMEX} template - NavierStokes::NavierStokes(const unsigned int degree) + InsIMEX::InsIMEX(parallel::distributed::Triangulation &tria) : viscosity(0.001), gamma(1), - degree(degree), - triangulation(Triangulation::maximum_smoothing), + degree(1), + triangulation(tria), fe(FE_Q(degree + 1), dim, FE_Q(degree), 1), dof_handler(triangulation), - quadrature_formula(degree+2), - face_quadrature_formula(degree+2), - time(1e-2, 1e-3), - timer(std::cout, TimerOutput::summary, TimerOutput::wall_times) + volume_quad_formula(degree + 2), + face_quad_formula(degree + 2), + mpi_communicator(MPI_COMM_WORLD), + pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0), + time(1e0, 1e-3, 1e-2, 1e-2), + timer( + mpi_communicator, pcout, TimerOutput::never, TimerOutput::wall_times) { } - // @sect4{NavierStokes::setup} + // @sect4{InsIMEX::setup_dofs} template - void NavierStokes::setup() + void InsIMEX::setup_dofs() { - timer.enter_subsection("Setup system"); + // The first step is to associate DoFs with a given mesh. dof_handler.distribute_dofs(fe); - DoFRenumbering::Cuthill_McKee(dof_handler); - // We renumber the components to have all velocity DoFs come before // the pressure DoFs to be able to split the solution vector in two blocks - // which are separately accessed + // which are separately accessed in the block preconditioner. + DoFRenumbering::Cuthill_McKee(dof_handler); std::vector block_component(dim + 1, 0); block_component[dim] = 1; DoFRenumbering::component_wise(dof_handler, block_component); - dofs_per_block.resize(2); DoFTools::count_dofs_per_block( dof_handler, dofs_per_block, block_component); + // Partitioning. unsigned int dof_u = dofs_per_block[0]; unsigned int dof_p = dofs_per_block[1]; + owned_partitioning.resize(2); + owned_partitioning[0] = dof_handler.locally_owned_dofs().get_view(0, dof_u); + owned_partitioning[1] = + dof_handler.locally_owned_dofs().get_view(dof_u, dof_u + dof_p); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + relevant_partitioning.resize(2); + relevant_partitioning[0] = locally_relevant_dofs.get_view(0, dof_u); + relevant_partitioning[1] = + locally_relevant_dofs.get_view(dof_u, dof_u + dof_p); + pcout << " Number of active fluid cells: " + << triangulation.n_global_active_cells() << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() << " (" + << dof_u << '+' << dof_p << ')' << std::endl; + } + + // @sect4{InsIMEX::make_constraints} + template + void InsIMEX::make_constraints() + { + // Because the equation is written in incremental form, two constraints + // are needed: nonzero constraint and zero constraint. + nonzero_constraints.clear(); + zero_constraints.clear(); + nonzero_constraints.reinit(locally_relevant_dofs); + zero_constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints); + DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints); + + // Apply Dirichlet boundary conditions on all boundaries except for the + // outlet. + std::vector dirichlet_bc_ids; + if (dim == 2) + dirichlet_bc_ids = std::vector{0, 2, 3, 4}; + else + dirichlet_bc_ids = std::vector{0, 2, 3, 4, 5, 6}; - // The Dirichlet boundary condition is applied to boundaries 0 and 1. FEValuesExtractors::Vector velocities(0); - { - nonzero_constraints.clear(); - - DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints); - VectorTools::interpolate_boundary_values(dof_handler, - 0, - BoundaryValues(), - nonzero_constraints, - fe.component_mask(velocities)); - VectorTools::interpolate_boundary_values(dof_handler, - 1, - BoundaryValues(), - nonzero_constraints, - fe.component_mask(velocities)); - } + for (auto id : dirichlet_bc_ids) + { + VectorTools::interpolate_boundary_values(dof_handler, + id, + BoundaryValues(), + nonzero_constraints, + fe.component_mask(velocities)); + VectorTools::interpolate_boundary_values( + dof_handler, + id, + Functions::ZeroFunction(dim + 1), + zero_constraints, + fe.component_mask(velocities)); + } nonzero_constraints.close(); - - { - zero_constraints.clear(); - - DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints); - VectorTools::interpolate_boundary_values( - dof_handler, - 0, - Functions::ZeroFunction(dim + 1), - zero_constraints, - fe.component_mask(velocities)); - VectorTools::interpolate_boundary_values( - dof_handler, - 1, - Functions::ZeroFunction(dim + 1), - zero_constraints, - fe.component_mask(velocities)); - } zero_constraints.close(); + } - std::cout << " Number of active cells: " << triangulation.n_active_cells() - << std::endl - << " Number of vertices: " << triangulation.n_vertices() - << std::endl - << " Number of degrees of freedom: " << dof_handler.n_dofs() - << " (" << dof_u << '+' << dof_p << ')' << std::endl; + // @sect4{InsIMEX::initialize_system} + template + void InsIMEX::initialize_system() + { + preconditioner.reset(); + system_matrix.clear(); + mass_matrix.clear(); + mass_schur.clear(); BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block); DoFTools::make_sparsity_pattern(dof_handler, dsp, nonzero_constraints); sparsity_pattern.copy_from(dsp); - - system_matrix.reinit(sparsity_pattern); - mass_matrix.reinit(sparsity_pattern); - - solution.reinit(dofs_per_block); - solution_increment.reinit(dofs_per_block); - system_rhs.reinit(dofs_per_block); - - timer.leave_subsection(); + SparsityTools::distribute_sparsity_pattern( + dsp, + dof_handler.locally_owned_dofs_per_processor(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(owned_partitioning, dsp, mpi_communicator); + mass_matrix.reinit(owned_partitioning, dsp, mpi_communicator); + + // Only the $(1, 1)$ block in the mass schur matrix is used. + // Compute the sparsity pattern for mass schur in advance. + // The only nonzero block has the same sparsity pattern as $BB^T$. + BlockDynamicSparsityPattern schur_dsp(dofs_per_block, dofs_per_block); + schur_dsp.block(1, 1).compute_mmult_pattern(sparsity_pattern.block(1, 0), + sparsity_pattern.block(0, 1)); + mass_schur.reinit(owned_partitioning, schur_dsp, mpi_communicator); + + // present_solution is ghosted because it is used in the + // output and mesh refinement functions. + present_solution.reinit( + owned_partitioning, relevant_partitioning, mpi_communicator); + // solution_increment is non-ghosted because the linear solver needs + // a completely distributed vector. + solution_increment.reinit(owned_partitioning, mpi_communicator); + // system_rhs is non-ghosted because it is only used in the linear + // solver and residual evaluation. + system_rhs.reinit(owned_partitioning, mpi_communicator); } - // @sect4{NavierStokes::setup} - - // A helper function to determine which constrint to use based on the current timestep + // @sect4{InsIMEX::assemble} + // + // Assemble the system matrix, mass matrix, and the RHS. + // It can be used to assemble the entire system or only the RHS. + // An additional option is added to determine whether nonzero + // constraints or zero constraints should be used. + // Note that we only need to assemble the LHS for twice: once with the nonzero + // constraint + // and once for zero constraint. But we must assemble the RHS at every time + // step. template - const ConstraintMatrix &NavierStokes::get_constraints() const + void InsIMEX::assemble(bool use_nonzero_constraints, + bool assemble_system) { - return time.get_timestep() == 0 ? nonzero_constraints : zero_constraints; - } + TimerOutput::Scope timer_section(timer, "Assemble system"); - // @sect4{NavierStokes::assemble} - - // Note that we only need to assemble the LHS for twice: once with the nonzero constraint - // and once for zero constraint. But we must assemble the RHS at every time step. - template - void NavierStokes::assemble(bool assemble_lhs) - { - timer.enter_subsection("Assemble system"); - if (assemble_lhs) + if (assemble_system) { system_matrix = 0; mass_matrix = 0; } - system_rhs = 0; FEValues fe_values(fe, - quadrature_formula, + volume_quad_formula, update_values | update_quadrature_points | update_JxW_values | update_gradients); + FEFaceValues fe_face_values(fe, + face_quad_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_q_points = volume_quad_formula.size(); const FEValuesExtractors::Vector velocities(0); const FEValuesExtractors::Scalar pressure(dim); @@ -663,225 +802,212 @@ namespace fluid std::vector> grad_phi_u(dofs_per_cell); std::vector phi_p(dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), - endc = dof_handler.end(); - - for (; cell != endc; ++cell) + for (auto cell = dof_handler.begin_active(); cell != dof_handler.end(); + ++cell) { - fe_values.reinit(cell); + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); - local_matrix = 0; - local_rhs = 0; - local_mass_matrix = 0; + if (assemble_system) + { + local_matrix = 0; + local_mass_matrix = 0; + } + local_rhs = 0; - fe_values[velocities].get_function_values(solution, - current_velocity_values); + fe_values[velocities].get_function_values(present_solution, + current_velocity_values); - fe_values[velocities].get_function_gradients( - solution, current_velocity_gradients); + fe_values[velocities].get_function_gradients( + present_solution, current_velocity_gradients); - fe_values[velocities].get_function_divergences( - solution, current_velocity_divergences); + fe_values[velocities].get_function_divergences( + present_solution, current_velocity_divergences); - fe_values[pressure].get_function_values(solution, - current_pressure_values); + fe_values[pressure].get_function_values(present_solution, + current_pressure_values); - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int k = 0; k < dofs_per_cell; ++k) + // Assemble the system matrix and mass matrix simultaneouly. + // The mass matrix only uses the $(0, 0)$ and $(1, 1)$ blocks. + for (unsigned int q = 0; q < n_q_points; ++q) { - div_phi_u[k] = fe_values[velocities].divergence(k, q); - grad_phi_u[k] = fe_values[velocities].gradient(k, q); - phi_u[k] = fe_values[velocities].value(k, q); - phi_p[k] = fe_values[pressure].value(k, q); - } + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + div_phi_u[k] = fe_values[velocities].divergence(k, q); + grad_phi_u[k] = fe_values[velocities].gradient(k, q); + phi_u[k] = fe_values[velocities].value(k, q); + phi_p[k] = fe_values[pressure].value(k, q); + } - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - if (assemble_lhs) + for (unsigned int i = 0; i < dofs_per_cell; ++i) { - for (unsigned int j = 0; j < dofs_per_cell; ++j) + if (assemble_system) { - // $LHS = a((u, p), (v, q))*dt + m(u, v) - // = ((grad_v, nu*grad_u) - (div_v, p) - (q, div_u))*dt + - // m(u, v)$ plus Grad-Div term. - local_matrix(i, j) += - ((viscosity * - scalar_product(grad_phi_u[j], grad_phi_u[i]) - - div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] + - gamma*div_phi_u[j]*div_phi_u[i]) * - time.get_delta_t() + - phi_u[i] * phi_u[j]) * - fe_values.JxW(q); - // Besides the velocity and pressure mass matrices, we also - // assemble $B^T$ and $B$ into the block mass matrix for convenience - // because we need to use them to compute the Schur complement. - // As a result $M = [M_u, B^T; B, M_p]$. - local_mass_matrix(i, j) += - (phi_u[i] * phi_u[j] + phi_p[i] * phi_p[j] - - div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) * - fe_values.JxW(q); + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + local_matrix(i, j) += + (viscosity * + scalar_product(grad_phi_u[j], grad_phi_u[i]) - + div_phi_u[i] * phi_p[j] - + phi_p[i] * div_phi_u[j] + + gamma * div_phi_u[j] * div_phi_u[i] + + phi_u[i] * phi_u[j] / time.get_delta_t()) * + fe_values.JxW(q); + local_mass_matrix(i, j) += + (phi_u[i] * phi_u[j] + phi_p[i] * phi_p[j]) * + fe_values.JxW(q); + } } + local_rhs(i) -= + (viscosity * scalar_product(current_velocity_gradients[q], + grad_phi_u[i]) - + current_velocity_divergences[q] * phi_p[i] - + current_pressure_values[q] * div_phi_u[i] + + gamma * current_velocity_divergences[q] * div_phi_u[i] + + current_velocity_values[q] * + current_velocity_gradients[q] * phi_u[i]) * + fe_values.JxW(q); } - // $RHS = - dt*[ a((u_prev, p_prev), (v, q)) + c(u_prev; u_prev, v)]$ - // plus Grad-Div term. - local_rhs(i) -= - (viscosity * scalar_product(current_velocity_gradients[q], - grad_phi_u[i]) - - current_velocity_divergences[q] * phi_p[i] - - current_pressure_values[q] * div_phi_u[i] + - current_velocity_gradients[q] * current_velocity_values[q] * - phi_u[i] + - gamma * current_velocity_divergences[q] * div_phi_u[i]) * - fe_values.JxW(q) * time.get_delta_t(); } - } - cell->get_dof_indices(local_dof_indices); + cell->get_dof_indices(local_dof_indices); - const ConstraintMatrix &constraints_used = get_constraints(); - - if (assemble_lhs) - { - constraints_used.distribute_local_to_global(local_matrix, - local_rhs, - local_dof_indices, - system_matrix, - system_rhs); - constraints_used.distribute_local_to_global(local_mass_matrix, - local_dof_indices, - mass_matrix); - } - else - { - constraints_used.distribute_local_to_global( - local_rhs, local_dof_indices, system_rhs); + const ConstraintMatrix &constraints_used = + use_nonzero_constraints ? nonzero_constraints : zero_constraints; + if (assemble_system) + { + constraints_used.distribute_local_to_global(local_matrix, + local_rhs, + local_dof_indices, + system_matrix, + system_rhs); + constraints_used.distribute_local_to_global( + local_mass_matrix, local_dof_indices, mass_matrix); + } + else + { + constraints_used.distribute_local_to_global( + local_rhs, local_dof_indices, system_rhs); + } } } - timer.leave_subsection(); - } - // @sect4{NavierStokes::solve_linear_system} + if (assemble_system) + { + system_matrix.compress(VectorOperation::add); + mass_matrix.compress(VectorOperation::add); + } + system_rhs.compress(VectorOperation::add); + } - // Only updates the preconditioners when we assemble the LHS of the system. + // @sect4{InsIMEX::solve} + // Solve the linear system using FGMRES solver with block preconditioner. + // After solving the linear system, the same ConstraintMatrix as used + // in assembly must be used again, to set the constrained value. + // The second argument is used to determine whether the block + // preconditioner should be reset or not. template - std::pair NavierStokes::solve_linear_system(bool update_precondition) + std::pair + InsIMEX::solve(bool use_nonzero_constraints, bool assemble_system) { - const ConstraintMatrix &constraints_used = get_constraints(); + if (assemble_system) + { + preconditioner.reset(new BlockSchurPreconditioner(timer, + gamma, + viscosity, + time.get_delta_t(), + owned_partitioning, + system_matrix, + mass_matrix, + mass_schur)); + } - if (update_precondition) - { - timer.enter_subsection("Precondition linear system"); - - preconditioner.reset(); - A_inverse.reset(); - preconditioner_A.reset(); - S_inverse.reset(); - Mp_inverse.reset(); - preconditioner_Mp.reset(); - Sm_inverse.reset(); - preconditioner_Sm.reset(); - approximate_Sm.reset(); - - approximate_Sm.reset(new ApproximateMassSchur(mass_matrix)); - preconditioner_Sm.reset(new PreconditionIdentity()); - Sm_inverse.reset(new InverseMatrix - (*approximate_Sm, *preconditioner_Sm)); - preconditioner_Mp.reset(new SparseILU()); - preconditioner_Mp->initialize(mass_matrix.block(1,1)); - Mp_inverse.reset(new InverseMatrix, SparseILU> - (mass_matrix.block(1,1), *preconditioner_Mp)); - S_inverse.reset(new SchurComplementInverse>(gamma, viscosity, time.get_delta_t(), *Sm_inverse, *Mp_inverse)); - preconditioner_A.reset(new typename InnerPreconditioner::type()); - preconditioner_A->initialize(system_matrix.block(0,0), - typename InnerPreconditioner::type::AdditionalData()); - - A_inverse.reset(new InverseMatrix, - typename InnerPreconditioner::type>(system_matrix.block(0,0), *preconditioner_A)); - preconditioner.reset(new BlockSchurPreconditioner< - typename InnerPreconditioner::type, PreconditionIdentity, - SparseILU>(system_matrix, *A_inverse, *S_inverse)); - - timer.leave_subsection(); - } + SolverControl solver_control( + system_matrix.m(), 1e-8 * system_rhs.l2_norm(), true); + // Because PETScWrappers::SolverGMRES only accepts preconditioner + // derived from PETScWrappers::PreconditionBase, + // we use dealii SolverFGMRES. + GrowingVectorMemory vector_memory; + SolverFGMRES gmres(solver_control, + vector_memory); - // Solve with GMRES solver. - timer.enter_subsection("Solve linear system"); - SolverControl solver_control(system_matrix.m(), - 1e-8 * system_rhs.l2_norm()); - GrowingVectorMemory> vector_memory; - SolverGMRES>::AdditionalData gmres_data; - gmres_data.max_n_tmp_vectors = 100; - SolverGMRES> gmres( - solver_control, vector_memory, gmres_data); + // The solution vector must be non-ghosted gmres.solve(system_matrix, solution_increment, system_rhs, *preconditioner); + const ConstraintMatrix &constraints_used = + use_nonzero_constraints ? nonzero_constraints : zero_constraints; constraints_used.distribute(solution_increment); - timer.leave_subsection(); return {solver_control.last_step(), solver_control.last_value()}; } - // @sect4{NavierStokes::run} - + // @sect4{InsIMEX::run} template - void NavierStokes::run() + void InsIMEX::run() { - create_triangulation(triangulation); - triangulation.refine_global(2); - setup(); - - std::ofstream out("grid.eps"); - GridOut grid_out; - grid_out.write_eps(triangulation, out); - - std::ofstream out2("force.txt"); - out2 << std::setw(13) << std::left << "Time/s" - << std::setw(13) << std::left << " Drag" << std::setw(13) - << std::left << " Lift" << std::endl; - - // In IMEX scheme we do not need to implement the Newton's method, what we need - // to do at every time step is simple: - // 1. Solve for the solution increment; 2. Update the solution. - output_results(time.get_timestep()); - while (time.current() <= time.end()) + pcout << "Running with PETSc on " + << Utilities::MPI::n_mpi_processes(mpi_communicator) + << " MPI rank(s)..." << std::endl; + + triangulation.refine_global(0); + setup_dofs(); + make_constraints(); + initialize_system(); + + // Time loop. + bool refined = false; + while (time.end() - time.current() > 1e-12) { - std::cout << "*****************************************" << std::endl; - std::cout << "Time = " << time.current() << std::endl; - - assemble(time.get_timestep() < 2); - - auto state = solve_linear_system(time.get_timestep() < 2); - solution.add(1.0, solution_increment); - - // solution is distributed using nonzero_constraints all the time - nonzero_constraints.distribute(solution); - solution_increment = 0; - - std::cout << " FGMRES steps = " << state.first - << " residual = " << std::setw(6) << state.second << std::endl; - + if (time.get_timestep() == 0) + { + output_results(0); + } time.increment(); - - if (time.get_timestep() % 1 == 0) - { - output_results(time.get_timestep()); - process_solution(out2); - } + std::cout.precision(6); + std::cout.width(12); + pcout << std::string(96, '*') << std::endl + << "Time step = " << time.get_timestep() + << ", at t = " << std::scientific << time.current() << std::endl; + // Resetting + solution_increment = 0; + // Only use nonzero constraints at the very first time step + bool apply_nonzero_constraints = (time.get_timestep() == 1); + // We have to assemble the LHS for the initial two time steps: + // once using nonzero_constraints, once using zero_constraints, + // as well as the steps imediately after mesh refinement. + bool assemble_system = (time.get_timestep() < 3 || refined); + refined = false; + assemble(apply_nonzero_constraints, assemble_system); + auto state = solve(apply_nonzero_constraints, assemble_system); + // Note we have to use a non-ghosted vector to do the addition. + PETScWrappers::MPI::BlockVector tmp; + tmp.reinit(owned_partitioning, mpi_communicator); + tmp = present_solution; + tmp += solution_increment; + present_solution = tmp; + pcout << std::scientific << std::left << " GMRES_ITR = " << std::setw(3) + << state.first << " GMRES_RES = " << state.second << std::endl; + // Output + if (time.time_to_output()) + { + output_results(time.get_timestep()); + } + if (time.time_to_refine()) + { + refine_mesh(0, 4); + refined = true; + } } - - out2.close(); } - // @sect4{NavierStokes::output_result} - + // @sect4{InsIMEX::output_result} + // template - void NavierStokes::output_results(const unsigned int output_index) const + void InsIMEX::output_results(const unsigned int output_index) const { - timer.enter_subsection("Output"); - std::cout << " Writing results..." << std::endl; + TimerOutput::Scope timer_section(timer, "Output results"); + pcout << "Writing results..." << std::endl; std::vector solution_names(dim, "velocity"); solution_names.push_back("pressure"); @@ -892,91 +1018,123 @@ namespace fluid DataComponentInterpretation::component_is_scalar); DataOut data_out; data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, + // vector to be output must be ghosted + data_out.add_data_vector(present_solution, solution_names, DataOut::type_dof_data, data_component_interpretation); - data_out.build_patches(); - - std::ostringstream filename; - filename << "Re100-" - << Utilities::int_to_string(output_index, 6) << ".vtu"; - std::ofstream output(filename.str().c_str()); - data_out.write_vtu(output); - timer.leave_subsection(); - } - - // @sect4{NavierStokes::process_solution} + // Partition + Vector subdomain(triangulation.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + { + subdomain(i) = triangulation.locally_owned_subdomain(); + } + data_out.add_data_vector(subdomain, "subdomain"); - // This function is used to calculate the drag and lift coefficients on the cylinder. - // We first calculate the traction of the fluid, which is nothing but the product of the - // stress tensor and the normal of the cylindrical surface, and then integrate it along - // the cylindrical surface and negate it. - template - void NavierStokes::process_solution(std::ofstream& out) const - { - timer.enter_subsection("Process solution"); - - Tensor<1, dim> force; + data_out.build_patches(degree + 1); - FEFaceValues fe_face_values(fe, - face_quadrature_formula, - update_values | update_quadrature_points | - update_JxW_values | update_normal_vectors | - update_gradients); + std::string basename = + "navierstokes" + Utilities::int_to_string(output_index, 6) + "-"; - const unsigned int n_q_points = face_quadrature_formula.size(); + std::string filename = + basename + + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4) + + ".vtu"; - const FEValuesExtractors::Vector velocities(0); - const FEValuesExtractors::Scalar pressure(dim); + std::ofstream output(filename); + data_out.write_vtu(output); - std::vector p(n_q_points); - std::vector> grad_sym_v(n_q_points); + static std::vector> times_and_names; + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + { + times_and_names.push_back( + {time.current(), + basename + Utilities::int_to_string(i, 4) + ".vtu"}); + } + std::ofstream pvd_output("navierstokes.pvd"); + DataOutBase::write_pvd_record(pvd_output, times_and_names); + } + } - for (auto cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell) - { - for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) + // @sect4{InsIMEX::refine_mesh} + // + template + void InsIMEX::refine_mesh(const unsigned int min_grid_level, + const unsigned int max_grid_level) + { + TimerOutput::Scope timer_section(timer, "Refine mesh"); + pcout << "Refining mesh..." << std::endl; + + Vector estimated_error_per_cell(triangulation.n_active_cells()); + FEValuesExtractors::Vector velocity(0); + KellyErrorEstimator::estimate(dof_handler, + face_quad_formula, + typename FunctionMap::type(), + present_solution, + estimated_error_per_cell, + fe.component_mask(velocity)); + parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction( + triangulation, estimated_error_per_cell, 0.6, 0.4); + if (triangulation.n_levels() > max_grid_level) { - if (cell->face(f)->at_boundary() && cell->face(f)->boundary_id() == 1) - { - fe_face_values.reinit(cell, f); - fe_face_values[pressure].get_function_values(solution, p); - fe_face_values[velocities].get_function_symmetric_gradients(solution, grad_sym_v); - for (unsigned int q = 0; q < n_q_points; ++q) + for (auto cell = triangulation.begin_active(max_grid_level); + cell != triangulation.end(); + ++cell) { - const Tensor<1, dim> &N = fe_face_values.normal_vector(q); - SymmetricTensor<2, dim> stress = -p[q]*Physics::Elasticity::StandardTensors::I - + viscosity*grad_sym_v[q]; - force -= stress*N*fe_face_values.JxW(q); + cell->clear_refine_flag(); } - } } - } + for (auto cell = triangulation.begin_active(min_grid_level); + cell != triangulation.end_active(min_grid_level); + ++cell) + { + cell->clear_coarsen_flag(); + } + + // Prepare to transfer + parallel::distributed::SolutionTransfer + trans(dof_handler); - double drag_coef = 2*force[0]/(0.1); - double lift_coef = 2*force[dim-1]/(0.1); + triangulation.prepare_coarsening_and_refinement(); - out.precision(6); - out.width(12); - - out << std::scientific << std::left << - time.current() << " " << drag_coef << " " << lift_coef << std::endl; + trans.prepare_for_coarsening_and_refinement(present_solution); - timer.leave_subsection(); + // Refine the mesh + triangulation.execute_coarsening_and_refinement(); + + // Reinitialize the system + setup_dofs(); + make_constraints(); + initialize_system(); + + // Transfer solution + // Need a non-ghosted vector for interpolation + PETScWrappers::MPI::BlockVector tmp(solution_increment); + tmp = 0; + trans.interpolate(tmp); + present_solution = tmp; } } // @sect3{main function} - -int main() +// +int main(int argc, char *argv[]) { try { using namespace dealii; using namespace fluid; - NavierStokes<2> flow(/* degree = */ 1); + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + parallel::distributed::Triangulation<2> tria(MPI_COMM_WORLD); + create_triangulation(tria); + InsIMEX<2> flow(tria); flow.run(); } catch (std::exception &exc) -- 2.39.5