From 687448e97df1a932c5541e48a55b2d3a3f95a3f2 Mon Sep 17 00:00:00 2001 From: nmannall Date: Thu, 4 Jul 2024 11:43:10 +0100 Subject: [PATCH] Fix all MPI ranks should use the same sigma max The CFS parameter sigma max values for PMLs should be the same for all PMLs. Therefore it needs to be broadcast to all MPI ranks. --- box_half_model.h5 | Bin 49648 -> 0 bytes box_half_model.in | 11 ---- gprMax/grid/cuda_grid.py | 10 ++- gprMax/grid/fdtd_grid.py | 129 +++++++++++++++++++++++++++++++++++-- gprMax/grid/mpi_grid.py | 77 ++++++++++++++++++++-- gprMax/grid/opencl_grid.py | 4 ++ gprMax/pml.py | 39 ++++++++++- 7 files changed, 244 insertions(+), 26 deletions(-) delete mode 100644 box_half_model.h5 delete mode 100644 box_half_model.in diff --git a/box_half_model.h5 b/box_half_model.h5 deleted file mode 100644 index 2d732fe5b301a91821baead99cd86e742e8f2a74..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 49648 zcmeF&30zIzzc~ETppip3d2Ot%Xo1RLD@u6wydg6iO&GCnQryN~Y#> z)>hXzjE1T5GTU*=z4zd%e#FZ+A~cISo0n zUyF>4n53A}U#opTe)cXs=YHq(eeT;H)Vm_yo2K-p`}+M}5f@Vs`?YOoZ@yA*`??ap zGlur2laPmp7+12d{l2ufnwYM}?=_jej6doBs0es_xPe~>Jhk_z`^Hma?>OrVMpWdS zs9)=mEB{LV8hgJ!kNfrGzpyCxj@Q4BYlYqx#l`xIsfx*n&5ei@%vvcXrYG_lWJgKKlRaL+>E|mG{r>NdMK&r1?u1hH&}+UOU78MLWMjpde=VD~H>sJhv_* z=Cd>;a@NxM;fsHVK}?C;=kMK)>n-C-YbP-1+&d`G3`~|8%>*M?_!07q1Ll938lLRbLPO zqwNRxR``3WAQtr3mS63wa&ZY>85kTL82nG|E9<)cQ(Ud&#jJkz21NWH9LKk^*KLB? z!1eFxKc7E|zs5oE(%&8WhwbqT;d_&-ojoI8`h+8Fv*X)&?B>!dFusCQib{kr{iUHn?&dbj`nB=&2(js3OW z*DvqB4gbYdTukc!)xW;y&dR?oK(XGlDkk>#{`Fp(V!FNM`gY&fuivA%FZaKm_VtHo z|866$x2b==A4vVaA4Es|x*hzx_k#hy3iWpW*AV-b!@#S zrtdiW4&Tbu+|+!wxme$Awzry?cSz*YRr6;p)(xH?5xy`)H+*?Wr0y);d9xOVnnXp+ zn)6quxn2C~%dh3H{&8tvJ@>z=iE+Pg{9Y7#Q?%01*4+57f1`V|{z`%d{>x8)*8fD{ zPXzu%;7dcD?WS`pW(Nz2N`(k@?RZ!SCx#>>ti= z{pVWyQ~Lk52>iaze$9^jUatS``kANMyP@yJ?Eh^6_-{I%zIlkh$LsI&I3PS~{;%n= zzTN)Ad6s`LkF#s%ucgfO@9DqyJl&`NdyjqRxo`h}(*GF|`1^cQ-*YviZ(^l443R7U znhyJKoDY-y|C|r&8TosMjlKK*_r|aOxZefMd&~bx{%1wt@8kEkIuH{Bz4x`gXn6dd zMi%pk{+;8|J7ezMoA2G5@7Ojp zdp`a>UVoo=W@at=56wIOHE-MdF#9!^```VXeGX4)uX}dH_3vrl{C?jJQoQ#Mj6QeK z=UVzw32ytZv7)T2(>v#^E2btRrWV^fBBZ$Aynb!(*RS7SLrSK%JqhtfgWpA__fd!5 z8Y#W$pXC3q2zYL~>=k;W*6VR%XG6lm6Ah1g+Pz#gTfF+;h;R6Gc3Q)yjTK&7mKS*~ z);M4PQ)W-SMYfk$>>FpV%`!=zjoXVn4Lpx}mTs@}vj_ujVilXpOlqz^Av&S(0aK|cM@ z>G%X$8u?Uwweq>@=;rhB4d187*WagdTbR$#eye;Dt{;$FsTB z%KKo^!CQtN@p9DO@oEL)zQe?2ecx&;`_>;(_YJ!|!gtqvV_$<6Hoj(tZoZ3G@qMSy zoa(##S+MUD<4E5OOH{G7d z51kyyUov$8-|xo?zEVj7fB)VM{Dp5a_)f(;_+wA!@uypt@uh|xiG zQ-Z-(=L8n>ngz2ZuLvRsb_qICo(Kw7z7>4E_e0<`M~YdqQ=ZWptioKnJd}CYsl)8E zHegnDnljFpZI~^(t_-i)i}}(xCPo#E=xNp!y`rAks*9P`kFPO{`rl=|;+`^M6W=j&*8gB+ z(51WKk-s`!j>p-ngBmp1JHv%K~=wq8Js>mpBPeY60WK;?z_%i4=(SG{ zg3JMM=5`g>ePjsKd87fqndrdz>3VRFtRW;jjo~==vCvlD1}^D#fHEDfaPc#DxJQcz zHwCaz?DQm9VCN4zg@I6SbqIXxwE))IEP>0cSHN*-4V>K)4}X42f=kR(V8W_Yn0R$F z%y!%giz>3P0ZBs}zbIIsoO?RY19=$KZmcCt*xXHQbVY7Al>; z0M~zPfHiKdFl8r!TfSU@DiIy9-^W`pW!*ivOz$z=d+izIXZJue|MxJ^<_jD(<|oWm z?1vOnCDAdp{>bEjEDCf{K(Hdh_^s{3KQqNFFmsu^8I7kN_y*V7omgpnL z7z6Zjx)E}+H9@lno1@}>mgvYkYjo#<9fH>#k=->Hw4lQcg|)jQ{aaorpofPHz6sEM zX^7~sNvOkd3i=HFQE~7LWU?U$Jv}l9Eollx-uLIDw@Qo9ESCs$YDyGxUcMaB;%M}= zH3r>(9*bPo#iQ-&i74n`5*li^9x)BcXsX`^H2eNWWV9#^MSR$V=ESC>ILR$YGifVY zD4U56q-CL2$!zp|^$yhWaVL5hk%NA8eKr%N9(FTWNB(c2& z)!Z*d(@e_H!I=GM*M$S9=Ep&_()JK)46i_;`G?V{)+6Zd%cICn{y1_SRf#^^oj|o- zCy}$@6k6b2h1QNcjjC*_kCMdTaMSF3meX&{Y)Jampq3~pF4-XFF%hG z)Gr|SLl@8{coBVmc@dS)zJz+(E}`4A>e1)R^=Q$^2DBro0ma^FK=%zAQA0>0QrX{# z)?aT#b9)-mb=4+hZPSEg9hy+CZ4+{^YeMnHO=yu)6MFl!5t(0VL}uF>(V|I>2&pup zb+rwsXnF&RkZ(YK2kOz+QS~UZ;1cRGzJwG@FQSLC7g2o31!Por9!aa7M`xy;LzcVh zP_F1Ka{g9}lnrW8muC&inR*7D46Q~NmYqgL;Z{IBS-$``Z{RBe#l}PW?F{Ijd z6zLWmK|ALjM#-ZpkoEO)bYsgwba}!7bm(pw@<`c-78;cz!_&p+s%H^e-By5vA^9lz z{$BJbU=Qloo{L6!<)D*icA}joJ5b~1?P$xROcXnF8(MgA3o;bTKo_TPM!cpp` zcxVGk+noGE$%=?sD2G1H z?T^yDBoR^;M+Pq6pi<-q_+e2GG!A_VXJtQt%1^uC3Y!kNbn9g(E-!)uq8s6`Pv_yO z+cnT=^C@^(`WQ@-JOrcieJ)FURGuei5 zCTuZy$&`3+Wz-ja60FspC-irkCsI0lK@=}5gKs@F!7pzke8GGV#azL=!tUZPTVLP@`JeEt+u}q4^(QYR6v)!mN~C1qV3IslomeXjBWu%i z$)xf6BsAWDgnu7RB+HG-0ihXj>}NqrD=djqgEbj{(3a>8w6}A6LGrW zOdeLb5MxzW@(H?H20!hS+iaj2TH@*Xj}3E5k{llrA|rPMa+E(IRvB8bm5{D7n062-y*(N|xuT5ZI2mwY^}bPvwonu8M;?ZB_Jv+(B;TXC^^Iu`GjhVA7y zU}Uf!PaU3!Bd5mU%3-UqiO(`zaX1_wNDaf@av|7$eITxz=8rX}OvZ-qSZwoyho7r> z;yf8Q+^pn;XUE&(tYQm%?3xLF<~ACSxTKGJn`nA&v#3Y+f@l+4Ejl*im?+}r0a44; zB9UEUuE@+POSD9Gvnb4Qy{IB%jp%1?xTxT6kVr`ki4ry3MXw^QMOxcOiTKxtilSxY zMcwM3gwGG&6iRQe6IPTL3jK|fh3NKFp;4}mFm*+9YxUW2t2yqGYAv=h_a5zFVq5nzD>jxeDG}vNz{g{ZcGqd<-G+0_opFs! zYK@3_yX`9DE78g1jT|d}LXn+a zrNq`AAH;SY8p8fuq0VLpYO&)jby%ZU!&#*qefGsP19nc=Xtq7nnAQ4h%G%8y%LbQN zvh@-+tgM?IduXBqE4RXl&53tm%O;LvJB-J(XZyRe{JkD*p{W-eea4%;HI&DGVBF8a3c zg6QPz4pH6dS0ZsM37lLu01pq=z*&Z)aH5_C?iufbZFcZ**C#(*wrviUaEQQlYHRRk z@pafCD;4MIZ^cn~2fn{%FIM&{!OaTg`0CSR_{gQxxbeg}EGllq(PbE$m9*o9$8KR# z{Q&R0`V5CWe1q-9zu-d~#Yn1#B)L`CpU7Cqletv`h;V`m5j6}Z-F$U&r%8+K@zf>h zr}fB2Qv+hK!-#DDZbIfx8%wSnu_8)8Y>AA!BMFLiA>K#aNcmNF;`z>t?33^%^->IR z{0hk8hZ9M8jUU;YJdO1Cm`D0&;NaBC>vGIH6@r zNM>0S2})f?3Z|?eHwQ!$QTZyOVH!h*=B^>TpRFa@N^ztsH;#mA$CKK$c#_Z)PueFW zkiZQI1k@xD%ZCX>^f`fOh$WK4?-R(f&IEGpPy)$~N+7`_6Nr05JlQ-so*cgtNBX(O zk!_Z-1&CPh$hBJ62{}6UKzZYjb z&%i!6<8i%O2=3qMfqflD;ZikeEP1(E^ddGzl%OduDz~y|ahC7T$nrKYes?gl>$E&; zcgBppWj~2k(Trx@#-y`n9_F)aZXRc|UN*2=qi?cHQeUzSa^hg(6*;i#geoAV+Tcc^ zA?Pl)0Gel=z>IcJQ1BT5IywN%pA`bCGr~bQjRvn}6M%zJGU)b619^+LfZVk0VE+Cb zP~MaWEFKkuU7z=Z&jSvF#oCpCZ*UrvjyVh5O)i4|!<)eU0U|K$OB>kOaScQt=mcde z?*dEvM*uv33chTA1xDGv13%7x0__vNgEyUGP<)^SOtX`Mg?0 zs3&g#*%^k=qhK_2x^4vD{1^jeG)-XYC{t*tX9gRT%;D(Q=Fsx&Sg0Lu0l$v5giG5j z;jYP8)bJ!3HYKuz`%e4P1BA z8a{}%hH-`Kv^UQ|wew@p9dQ?2 z*l-gRr(FeOqcAw=-UJlBod-30&w%VPCqQ#y1+bJX15Z%_h}x6`Qch%ofi0WBV!9sC z^|4?Fe+e-BG8;TgnFP+g_5jTSTfl5I1T{z1LDhLVpq=uO9XYjw{Ul$_UMR_B+x%l# ze`!~il#8=svr?GOXKnnMrPD=@HQl1KRTlVqas-w%$i~S(D)IJPSFn2VTfEjxmcWFe z#HH1Mh$L*u9BWT<6HX>01A@u>Nl|2{W<0S@NFk>?(}~b1o7|bVhy2)8OnPb#5`d49 z7<`&ssXIqncQ=wj^DxQNX(u<&brQhbCr>YRlfY3uM108yQk(aU1kMnrSEozT1kj&G zSj*9Gii&jCwSlx|n+h#)8%zuDsL{*e8r1N$7S))oL-$=6PB*IRQ%S#3blW;Z8dqXO zXH^>0fMceV*=tT6S6k4{6RhZiuh!Hc-_dhbp(PI*v+Fi)D<-;3Vd;ziYEyy>ik-n8+IH!b<>O`nhOp?BxAR@7wv02?JnBeC@*SvdyB*a6wsc9AH7!xKqN#xv)GODVYBZbD zfNo>@`m+&T@!gPad^?IN+|s9_6CqP@ zi5`t4i%BRM*6vT`LnmJG!G(wjwvnTtf`voeI%*(U^##ABGS$wDUJpd!m?TCsy~ zPGd*CU(0?uvW2~uQOGV;tYj@6>sd>QYi#oT$E?EQkL<)W2@ql=4<==)fO9SCz~Qnk z*tFLGjIcEY6+5lK(Rv4^jep;q_CmLRyIkLdO3e6wyPR(C=Cld_b~(e|&AIOXcDcT}>w_NucDeuV zt^e!htn|Hq_x*qGIjjG5@rVoh6OUhC5B`0{Nr-p$&cgp0AHT+hgt*<7KjY(Hj*rZ; z{}bZzyN`c<{)xbUe*}Knf{$xG*e&KAUQwjf2P#7k}FO|M6sEtt&N0!hvmL~^!Z9cK&bI9t%4 zvjq*DEzsa>foP1EcVXjX?>5dBe3uUK9?#i=>!|NL4GrLa5C)e7@C;pq2kFA57&!g9TA5*@+ z&$+BH9|`eQK2n@5*uFK*hvjU+P0ki5aJFE%ZKcl!&K7*(Y{4PU791{m;uFc)g3CiC zcvEWSd7%layziVX5bZJGz2|Jf2PJ2oKW7Uvc2DM+bGD#=bOev?h~-HaB=ZdKrSt4L zTQK))KJRkV0UnRD1)%}wcy67oysM@iJnzYmc(oeuc(eK9zRQ2g_`2yT`#K%tY(d)y z-zoEqeYK-)e0LhS`9530_syF%)pyMEVBbZfBYi(}w%~C@vhOm-bYDS4wl7~P&v&X( zneW?{qrNkqpYg5V*5G?f*yd|h)a83z^O^5^$4|b;q$Kz|d2)Pz&KBJJq0UcPIf8G` z*@6#et@u{CF8q5Rz4_){6ZuCc2J)o?7Vy_{wji@8fq$E`1?;;F{=1?beBV=fd|jI| z{shhzlt)$bW85zBCv3!g>xnn`{zVV@hVifYRXyMNWVMvw1!oJ+tQsU(z}W&r&K4{k zI7X0LZYj8Q+)40pl&9dj3=oVB@)xvmwm@@JgdluYjNsI)B*B2&8wJuOTLe{{EtsvH zC%AC4R51T_h2Sk`3ra1{32Hc7@S3v)HHuw=QyZQL#;p{h2#GbC}^ei+-Rd+k%E zn6m{S`3Ey|vm`s8vjz9h4rIm42D7D2n(RT&78F+-vTB?ym{?@Pj{M@x4%q6!W-Sr0 zF}o(S9r82T-JC6mqT%cn&K5*2iD$hyTQDa-l}+PpL2~ML_LWyIdxf(FMra>fTY89{ zYgNh0y{~51JU`DK8s5zMXHj<2gbsF-@m=;qKsUQX^qO6f_=PpC5eH#`(m-S&59E}U zz^oC2ffHv7#-G;#&HRzzh{PD6`p6t?d}jmNES$lK-4j5*H4l942Z6HWR8VZm*@Dba zz?d%v8zq;63fZ+_8)plS%uWUgHK{<#HUr!eWdg-@J3+Sp9^g8;5X=wX2L_gu1K)nf z0JH8C7&xXD%;Ic8*3BmH;S~ncjN8Ga4V}P5={{i3b%P~4d%*dekAQda2RNxH0S^#q zc&1biUfnhT+H$sF=+Pmt^PvV*Fw=pv$6=R{(J>SE$x6UA8x@N@%Ny-{$n`d`ZJiD)dP0~y@!~y1tznw9p7;9W?vea5Vj( zKB|i|Krf&XI%8vkZfKjM-C~w#_$zB9|JV){b~>Vf%Pwf_6*m;#>5l5|d7(YsJhbbJ z07)rAbWnQ|+HN%kePsR7p+z%LWO5J^md!y~!cY{}Js&j;ScI-vN1*M~qEO}96w8v(4Ey!yZeuc%oU#E8xw{cnEKEaXA2%VEvjsW*wxE?8 zwxTthEx47Kh16xT(SUV3P-o9h)D@hA67J`sMz9CTAKr@!hU6p5_yW{oy`9)`(h48qu+< zT>8Bco#qUInqd=)n9zg*$2XxvMonn+$R^~X*@PavZ$zt3H=-lSjYuP?5zSI=M6aqF z(4zSb$WyifJEaM`&@$&*2A9Vr6HJ(Qs70#pG6VIXRoGoy< zeiq&MQHz8cwdnZx8l*M%44NNSjr?YvMhhaU(ASlxkTPcru9}@dQ#o5ORN**MB}Wm+ zK7y1(4kKOt3iR@BIigz+qWca9(3i($X#2)}XyV9HwCYqba`7rcE;kC$ve0}~erqr4 zLfrZ3%0=BE2i-ch69pRVKua^XqxAQgC{C~qEk3mcjUJkT`pwvk@+#BNz#*yVvi>OyI`aS~6?MUb`W^62*=0DFvjw$l8sSXk3sAGF2409h z1!X*rLG=NLVCAiSaMF)L*fwS_ykN8oz7NQPmriBC7^O70cS$l#=4^qH?HYJ*(Nb7{ zeF2=FJ_ou51;A#5iO|Tx7qV+6z|3ck&~l~~d~(AW^3INg8y;&zi`QyU>6H>Roy*<( zhD*W*zpucg<|SzBe;@Rpc@3;jYyky^bzoFfC8(<25BB=xfs0Dp!PDodz$iBXOxYU= zW-9~%fj0wcA2@-zuEt=LrWWX$rvT)gzpy3FH(53Rv)mkF5nDeeg^jQXVD|;TaBQc=LXMS^uRt!8r~09M8dlO8are!&5l7 ztpV?DxPpfy-N6>6FYwrs&-g*II7xabO@1aS5O7(EbY`iN=WmCSQxU^R#!p=m@I;TK z-5W*JBaMhBFd?dbW@LEMSQ7Zml2pmtkoWs+$yyV8GIWIl@i^m1-hXu>3Fn+i*&-KW zq3ue7_PdgE*IY@9+&B{EIF5K|jU(SbyAq97S2EYyl}Mg(A=WM~EGx?rhIfH zFF0FZ_tBo1J+vdQ%Wa8Af(>yqvnH1hTatszEXX^Rv1IH@GcvQmgsihOCb-&&j1?PA zw1yiHtGJP*;Ful}JQ_}_2I!KQ;o3x!(IUAkG{}R9p(JVL5c2JWD!Fx7g+!)vwxDbP zk-4Toe5>S$1@e}h$ z*nP}Bd~$LZ&b)d9Px;=C!#=g)K`h2OZ7tYEqY;k|zKCmubr`$U;Nk;SIMw$AK6dda zUf@}QZ*?EQFn1q*YFmtpzU5=bmOVH#DF=_T-GQM;7Ot7F6Ub|1iA%Suh@36o|_f`(t_C$yikb;C3}%{N$Af7JM3q zv-nQ9+20nk+br;vcP4m)&1k$bOCLXS*TvHhYvQFkYWUkZ73`5U03WuH!~1%qaPnDk ztRVeW^!36Uk)Omf(d62DqKC{4QN`vq(TbPNB0bp)qIpKuqLJ&4i7HADh&I?2iIQIA zibjNFi7Xau5{2>Ci>%$&iuSDz7wujVB$^zHL=R?95KXkQ7G3l+5V`8AiwdRWMAtn( z3Ehr%3e6|i2|HF53dP2z2<>+H3757F6aH)=t+&f2v}Oiq`+7c&=PQqS!5>s$Civko zS0G)mS#U}HfZ%v$gTSlwu3+!Z?}DUJicDp@1~Us9GJY{OjE9~(({unb$2$USjL9f5nXU`@mdV@Qn%NiL>QvC0V^A{aN8&Io7mMk)3);iQQN_ zh&{by2)jH_o#oHbVmn6YuMe_|M3KNa-IhZt-V-{I&b!#B9EOu-i+5)uyL0aS+Vw zE0_~&DhQt~&aZoRzcutzqHt}uNT|DJfGA7PMik~bQ)H5zC~7;pLv->&xoCXI1<`}0 z9U|%0S0V>(39Nf^0Okj2U@6WPBxzb;{YfrZCYgtaeDTA<8FTQ72@!aofd57Brp4(+;1*()$~6N)g67``Yn|!?&FBxl1ka%94EV*Y$4X-}R-qE&*)`ovHY!r1~>`vs&aY7se?9Zv3-FCpEX zGSWPC1(~7}O&a&DB06RT0S)ru!eI9u@IRXhptOd!jX6Ua)= z7U25{#Og}|$^4l>R((t$_iiL`_s9fd$k_r(!vr#~KAzl}6i>uBTX2@M1$JX&NzlbL zWSLP6aag#D$noH($wqW}GAd;3agRJ>B zjcj!HBhvE_S5oa?F=3@{t&jR*^BSL&%m`^@%X@( zV0`_B2R`jM3h$Da#;z^RBHr;7(TOrSQPzmDEuUXWGe5^~U~FXw<6b1s7Brc$2gXfe zr>aM@sl(G*tvC5>(v9Qnv}XFF<76)ABr_=X&KSfvU!7Hb37y@p`UK?@Ld z#t97LY{9*60GJsCfX)6P;NA9c@b*eHxY$1dD2z@9Mf^0dY|$35o3jNU59WYzoGs|M zR}9#1`$3Y@VQ@vO610vw4KA3R1wV~00*12%uT(@}+^04$yZsutSl$UjR^0`=Y##yH z7f*rbj#uD|kcUyzK6g1>)!G@)>&|$qi%;s#t+q8jjUV<`= zTs#PxGK1kpeKi>TdMK1RrU7ZF7W7iqhVwXEQ0F!rhH<92^O6$Qi&xoGl3F zY(e!kBREOS7>a3_KuyjTICHk(i?TU9`N|xQsv8S;#aTc-3rqNpTEc^KtYFatE9k}9 zf)dUaq)XetHg_A?!`Xu488&eF9vdjZ*@Egs8z`J^164R%(2uhPv*WCx0cQ)mvaDc( zlod?kY{4YX7PN4-z=^X36`U<#I9nh-#uPeow!oRQ1v@xf@QAYo_n#ZUcFq>OeWwq5 zX6iw?Q^TPtXA5p}wxEi$1q^2k>N#7G$Jv6m$*OR@u?j4GJrG{qrwG$HTcE(%0&~t5 zSU;11zML&E;cUST$xqH|GFeWhR)}x(SHzdf=843-*I0z~R$uus>-M2>jpyv>98l zGSv|5KdcUxa<<@n;zxF?e+PR|shXWpmdmEQ#jtaSxUydli?g;^8Y-YtCsPwBjZLkr*H2b?X4OB_lTX3Za<;&Tvjydr z_Vf{F3yL{gFpaYXX`C%sCN_>{>$%Z^*5j$3rS<*Jm}H+o-~=W1+|&k;LF*9yJenK zfwKk2Ia`p$*#a@n7OdfH0X^!&K3;n?@Zf{IMOG+4)ozwJIZHmsfe=$ zb%U&^4`&NpI9u?t$&|W3F{Z~jTad`vf)dUae7K=cb54$+sq1y=H*am4@l%tYt&+UsY{3f77Km}SKyLGHQkjxXmT|Tq zowEf`BG!?{j#A>81>6tQyH1E(M!e@gB;A}xj_Gf&i z>>mE|ss(3-RN{e(d+{^Q77URI#gj4}@Q^bDu}fmJ=tRXL(FdIjVQYJ(AaP&}li(v{ zQjaRKQSw%7*p+GQ-B)YbE6295*Ha4F%t4hb&!V0+|Jlx(akk)H#77nrH2sHLyq88oQdC{`*kdL!>k6ydi_Rua+K*20`mS)V;nJPFhA|pW7eh6lT{2VFuV1iDz2Wotn+*%hk{jsAAzovvTkGfB zYBWZ%jScc4lAhgHG#d(SdKzp*v|+kQUPDw}PQ&mq>l#i6UTsjcHfq$FBGFj2D5=3{ z!RH237p8Gru4UuT56Mrcpt zp?O;x;q85m`(n;F^7jcFx8ROOBeAZ=bb7Z@%SXJ)d9ORH-xRcmRENp5d8 z6y9%6NqpVB(&=aOrbp5(vLypr7Th1w!iaQQChCrENnU2rf)_ird|%se6QdVXPTtGifZ zYgoUw)=A4cThl`xw<@&1X?@l5qxHppDd7!U1>r;aK|+DKhA?vfaA9loXrXcbSfOx) zgAhqg5Y`I#LU}Pip}lsHP-kD5@Ji`2q3*Ug;p&X_!hB(x@F2TQI7WS!aK`9-;qCSN zgl@Kng%>PO3C}0i30Ipn2@7r7gax~A3Rlm1D2z{jCH$!UML2M4KT({Qj40B6favCY zRZ-+M4N*WrndrT~t?0``7txaso+5>5tmwAvRM8UoS)zWyVWOStOGN_> z*NUe5B#DZSZV=fE(nUwqvqXLxIigQ~`J$7eQjy`xL!xbgl_IZ=XG98bFNg|vwTL9v zUly5E-4w-IJP;+meJ(mC@qzPGKSicVk~lnG7RMYIfIs~hgb&83V`l>$>?=7EuU8&} z@A{0zBTm@jr&C;ToUsQ!Y|Y2>S5CzH-}>X&qS?4kAMoji-|@?)egvJCCc!OoWJ}KgqBC+3kzK7uYF=rQFUxev zXO)p;N!w`Raln-LWLuJx`F6w>JCl>D6UfUjZ&LD5Kvu*}Bwm)&$V~A-qWUbDOnVQwTmT6>HQ93zhvW=j|?c~^^UF5tL=NLQkN%Y}j z63so2PR5iImxv?8Gp3T1=Two5Gc{y-_c=1ysGgXuY$n-6M6B#D6HtDQd^G4JiO24c zr_4hl`Lvs?+y0VN_`N0aYM)5wwQuB0hZuEcB&c4M6n(&XybDL=sE4^CJ#}~>onWm( zC)KFZuU=}@?3y|iglbWPx7yTf-EeBHp-)%sA4M6b(NvC%p_cVe@afm<3_WQ9>5Yhq^lJ5F`a)|e zo$2XM2Tz|)J9o~cs^^1f=c(B=^+YgjJsnEF70sh|(P7lkWg)GXSWG`24X2a6me3nl zBI!1?l%6=bjAjpBL4U9-Y0Khhx^??1I=W~zO-_iRic{CnYVEc3O4V9Anu(>KZ^Y7l zgX1XAD2^(7#L-f}IQorCbM@k=u}mDbZHuK(#>7$uw3a?NyN23Jt)b3st7!t)rpNMV zI$+>R`fkN?IaZn+GXSEJ2w~FcGRAZ8*ELZCYw{Kh%wY}@hH05VmPgD z(WFa)htU01O0@rNIXeHjBvs%2oxC-BO{&u$5a|Wi$nMe2#O2T#GN|(qnQ|$Q`+n3` za&1o{ar`u&bS?5F%Z3>fSqnvyb+;2g*}e-OzhH^q-EtCbxw4-rD4oJeIUHkGqxM+!zPWoeKqJ=aRr&>Fr?km3_cf@eJrbAOs%|-v*ik-U79U{h+Fc0`$Bw7`8g< zK&MDU=%8!?CBHbr;Da8p<~jp)J*Prn{n>EWy9Llfei@XSz7~3DC&8sZQsB+(O|VjO zEA+U%9j-r_18=}QxOYbp+-+V4@OyQdw6O{pi~CB14`xa%w&QF{Sy>S=)MEL&lG zD1nlJSKvo@9qyXb2~#HBf!;m$;qq;d;qkT4;K-4$;N74%(B;5;c*XxSG&%JRo=W)% zeZ0ldgFF3D*d$4GvRn%B2KGmuGh~qS7FqO-$f4;H3MjOnBAO~T0KL%}h+MRl(3?j} zsAq>VI<-y(otroam6JheNQ^2vt}z&WS~(c$JR6KAL=8cb#0N;ouFApxDEEaInE`xM%)N z=>7N_jDB|+zP*Z}+xHe&W7G&^?JvTL4`*R|MKyG-JqgQ%$KdWc6>zEC0hrlQ0;LWW z!0@0w@J!=Qcrq~y{)pTH$91Q{!Z|5$s$mjrNsooA_bi9mR~Nw-gF~QMz;t*pZX$GS z;X$79c(?=E!=a_-@Vl=e9G9g7T}}*v!zT}bn_Q&f#i*~~g!ohNPS^oP4{HW0`c>fk z@_j%)Y$sTAegoK090|@$oCx;Y*?{VDHPCRNhrKlCB#W|R*wR~)tga|eFlm0d=%J=J zcGt+o%knU;u@NUK-J0amVQVrcNkGQ^2qF5(Yf0GkjYLanJ87StPs*wflBJ_hk`E)! z5$6{zQwoe z7CjQ8ODneN(T6FcXkzwg`lHyGeswaV?_0;xDG63|_hlO@@!gKz@^Pd;Bb{m11Xp@P z!j1ag7*F$K-094(9(2nbFWPXNv)f%h^s}5Vy+4#sCqCfQ>?#2b-N?{J1D2W=v2?gD zprUnvj%eV%mLUu2a}!8U8bCT@IHXI}xb+u+I^+SG;ST8acPyQ?iKXrmEZsMsp*Gb5 ziUtX&&tg6um+nhv@8QvcW*>UsjyGL%*^9agJ*n`d2YtQHok|BxpzZeK>D4de=wOZ3^Z_c-$vY^xT#!~xiGa4ajO5=i!>EJpe zYNj)qN`)KHvc{3L%s`*oM2w&h%5~{cDIK~bU>M!LTa&8YQm3_gLutyQAvF4oD#hPb z=u2B=x;9}T)xWGrpC~C%sVQ>wZIKLp`%Idiu#}<`*GbUFZQ?XrM~qh8`9{8O{!FqQ zJ`mOJH)L<_E3$U-bJ8@Zn|$kfNE8p;BVtQ#le2D}B+#^joKk8hDo@(T^}CqV-E1XG zTAN62dp!wzdVy4Q*O9wowWQOmntVW~NY~OzVtV)}ne(E8^ca;B@0k5W^>Qh3)+#2q z;tI&JcY8_AjNL>i+C@IY9c1t2EYcjbjg0x6L5^+ROw^6jNMy$b(vp!(K8#;SvcD&g zL5*=FSa%KK)vqLU-BL2wDS~|axqvKwHkTyRIV8Oz*_67F>eS>ekk z;v2k)I(8=_b2l<6%bAQ4b0A5}ZAgfeCE1;BM)J&!iONkwlD}~zIps2fM5ySH1>E@0 zs2xfY@&}VadsRqF)<7bjqd=}6l_m3~Vs!pB7RElu%Z5M3 z1K!@n`QJLRsN)*;z21fgD~j;)aZPyO!i#v-t+V)^XEh#w-~@i=eiY}+9Ktnw%kZI) zV!UleKDTW*cDCM$<(6jQRlBxeuj);B?dXko)c*B&y+k5zbBM(j7gu60`ADo#vH+)5 zhhV8~fq34wY1qbW5_Ubp-~k3cc&v##PMGeBW9K>GPo37-rfe*BK_-|TKN>5SjKms> zBe2>kZ7e@m6Nf|%#Q~Xvv8K`>Y*9wu zU7n8lrE_pzPZ(}ah{QAXSK(JG@i^Xe9X_yl11{~@gpY@9#i8G};}3N?cynnUE<09? z$9C+;yLAs^gN2ni`t@lXwW|(0``6>MI<0uY&^BD5c^zj9x^UBu2iWf2Q~a6#8ozq_ z5s%IMfm@6E5jk0D;ua@MnqDiA9`1Y7P6Jg*Yk?a1Y^zBU9%+;0j1ff6Z4{ZzeQ)~1 z851%!YAmTAXHBjRv?tTAI+4J(o-|EO zA&*{fATG(NWQp7+g4c7`%8zt%Yu*+Tduc0~qMb>Gg=UeRdE3eO%h}}1*BwN9z%C*? zG>0r#%_V2Wb`#56yNS=?Jw!EZFS$M>kJy~dBX>RW$?EESvT9TT@rx-SN@okm$5#cU zbzmWBQ7I&bGKJ*&(*lxpynu*@7LeC61!P1rTh-cAW^0;LW(Ri?%6yMJ! z@)vVRzl>cZ%4;X#-O487Gq#iMjhV!6&^EGX{1)xMNHvOR+AShx%5zEg(I8@M=}(-)CldcKJ~5l%LDt`OCcABI$hP^W z#Av}NA~8~%v>hBwsEQ&PJyw!wdlN&e?4*I|q8ZI<5EK?sTHBch0K?h#a+Eu z+Th$`PjA^jx4i=#rumdCboKd=AM9f-N%CP<7y0mTtB<|jQ=c&O-6wp#65nLoOg@pd zfZzAdn*WRA%)ce%@=bTH;h&qaji1xAk6-M2f*$~5_yL@=d}s1JfAjhreq_*PetzT) ze)gpbzG`C~-|ty7zf(rUkJY=&|LoDjpRuEtpBnv=ADTDBAFKbucf>#WmmdA&Pu?d- zvcD*jD|e}okprq^#fs_Vm)kSRr}}e9y_4Ex<^x^wl#xD3?J**2ZkUo;&n(DWlWoW) zX7=PG*Tv)u&Qda)wTujRa3!nfd64Tyy~x~V78xGFB^R3l@=%92S+$ZRkKS8BzH(hn z9?kM4|Gi#EzR=x3#xggNLB5+wYu_!T31b^sth=4;d%lCba&{M4X1bePT(XB$*4j%Z zhK7)d6`^GGqkUxR(0(%U^#PLTI7psPJVXw=9wxW;93~aXF!E4s7`a9!oV2$JCl@io z$@}i%WUo#*`DGxCgo$CKyjB>wD*iCJY48yFOXm=2?0Arjbv;0CFxpR^`VdN9$POU` ztoM@h^Y)M_KZD3cOLvhQSMDGcRtJ)CZd=K{lLN@>)tksj&IU5)&N|X@rY{+`Viidr zUQRAJ;zOEj1LS%KE*brfNmfQq?62hA$fY}$ku%C1NTU}E$p3c?JE8J73tml-f;;( zYiToo7*z9Da z|IPE+6}ic$#9+)@;hqQ@Y)!xGlVa1k$U45$@@Z4&^2Se{T<_^!8u!h+6t*+EHcZ;s zRe9xb*Zpg;UG(AeU9F6Qu7Fu(T}lIuU1bI0u6rAMy4uuVbRDl9?)tj&ch}5sa?;Os zDpIHI>e8BYP3fyf9jRrHzI5M9Q>oHxYiZ1rMbhjxCuw=EyVQR>OFDM}klvG!Q` zqF>cp>3jEU(v*ZE=~>}*Y3-*Q(%VWmrT(hK*hqd_LZMZ`@hxxH9QdEAihRlY8L)O<~v zzV)hf$B!%0gW;E@>xc5CPTV}{&-9DZPdzzOmRz=Uu3n~ety?Be`&n*L2rySabL7F|5v0m>EH>ep74lNeEN{|qFty|)D|S2?Gz}r$nlq^ zy;~<8*Igw&;o&3w?9G)zj+eB?+C^%o>A-Ph9ZXl7oFV*G}Tg zMda|SFaG6Uo})mvA{BBlNS)+;n?<&4(I)>k>X9G+7?J-OT960qY{??S#pL{PN3yiu zmCW?>B9A<0lTBwBm);kkp;G=NN??EQsH|vIj8Y7xjtl~higVG zdFx6n`NklQoR%C%CVz+{H!hARson8pZgM<1voM~lxEfCmosTCI_s5f|Zt>*9u{iQs zK^)oZ7)MTOi6uFUVoAr07*aBHnzXh&J)xD3CYz6)B7+m7$o~=}Nrhu4$+;U&kXP0o zCo`TNA)8#o$+wXcx$VY&QoA~YbjaUB{tVqkrkQLfyIQu8FWmjf%)<5L_-9}8|34G( zzlqTN@6YyoF%R0=@1`(H|2*(QTPF*4HN9?EYrf-Us}k>J@oJzw`(ADP1+JTy=QkNI zP37I~4$AEIXML(|r#=yF?Qg8wR0fjUdaDn&+duSccQ(G^p=oZ}ws-Vu`xe!dcICXF zHa(r*wztj+?Q1e)+wUj_x8HJ&Xz$Nk(Qe>!vV9}n)n4&g-2QR?iFSTOLHkhH-}d~( zKkX?-kJ}dwe`?>Op(qF*QWtiaC7Htg?+b(_cWw$@?$!ur z?`;+8+(bgNBaei?j9&?VlzbKTc*=;598(hMtyU8?H%KuNFNFyCWKC z6^mpF?u&M-^^1P1ycZqM`ysO0B`c2dRTdXrQ5QG!=ZK}N^~656O~nD{?ZmsDJBjnM z8RE5_-r_&U){1wOZxw$J3lUEs=HmEAr^F!hjCgN(idai0OFZ14FJ2*gLtJ*MQXIFl zQM~@PK&*T6t~jRhk=WnwrTBQ%uz1w!m)ONtPSSo>S)#KNiNnVNy545l5YK6 z$s6Vs$>xe;$?2df$)$)!Nz$`+Nn|RL3|#Axcx(4dqK4l}It{)_Cbx}C-t3U4J{hY} zHahCmBdAH`wCGUH`wggB;t&>5cJ4bkfvP+Nq(6&Wx+4^^P^u2A4YMq#>FPS$B_4 z`tgYVC+Meh>tEBQ4?oiV<|DLT$sbzNLKZd6QAFGIRZxk?G=vXppj#46biqLz72cVT z_V^j1&{<|k@ZJhJ_uHdY9~{tMbr&>kwI@ocXQM-o5b1O;M<~V@jc?qDa@K7{bAomu z?eifh#^v#egZmsJqe{POGR?M>BuxC3;kJl5fxe$pnIxU z(M;JJXxZda#M7=sv`Y(V};sVBc6{Uf6&^T=@erlBJ39EfV2Ia@Byw1 z{%-1yWwg97F_VSw%5ibb5P{D>@W$Wnli20Wh8Fe9fkF)Phpj@ zX#CA81}AQa#oA}%uv%?AHhgpj8+=c|uawT>pK3`sXx2GgtC@@!tEAx4Zz(wJQ7S%m zCk@*rp2tsF>G-+&1>Dwk0Uz6+fsd{BW)7B5 z%)x)C931*J2d`4k#m*+VSYVrr6LfR&PUT!&@HPi;EX~0M5ji-(HU~S@Pu$PX#siV?zdZxOf4dzMYN_(dTjN%QSrHMJk@&l7cUsPsY0g&SB*RN!Wr) z#6Q<1;E?z6_{HuxYZ|L_FPesc^z_#TcIzd4Ar?uX(zA$#ya z-JSSC%2v#Mxe04;^TWC8SL4sxBwoH9;6cd7W@R4u>IP?gaC`|SH{0Qg7z@1Qju8$S znUA-A*TT2o%)p(6s`wsT33tz#gdY`+pq})%NX6kPYLCIF_h1W>gjAu`2MUq(csg2s zJq~>j4MUG3{Lt4D7c_ii0a`6kL|gAk=?z8a=`5M0bk&k1YKE^X&&SxDFnf|f^wkTA z{%Be7%1sZ9>~RM7hu4CdrC}iEM?8?O$^@=0*FbP%4d}bv0VG%MgRJZUFp%{bJoEnt zB6ldl=?1Eh%$*5?pJ+oHcYRp+&=iJ!v4NgY?G1 z21t}O!1$?+Fx#XNE|}2>cQYEHa!@0ToznY$TW9h3{Mg=cbVV6)^F1U9!|2CEuI&8~tUl`7#fn+kZ+u^gW4Ermz-mB6QdH=!!K z82%AlhZ|W%Fmm83{BgJt9v`_3xx51SxFQcmSYCqbBXVJPTQ)qYm<973GvJa#=`gY| z4bBjzz)kPY!4tAc&}?QRtT#CWg_d!!$0i2eag2t3wox#8-bvVHc^rPSh=7lO9)yJz zp-{1C5B#0H6XqFigMXj;L+a~#_^rto-tb%rZ%^SvD^&o0wzA>$jb1PUxWP?!%iyNI zCD6ao93c=@m&v}}3=+_p=>-18Jr*wGHSo0{uZd5S_Y2B6oHiQ`C#+9 zY~XEv9%_}d>2o>9RdVEq>0SGoqc7(rlB;0ZPycL3{lTLEIg0JIfq zfyw*TK+6h6@Ym@#QGM$zk^A2R;x(s(NFA#tTB@!P2P)HuU+5I^ZdwrW@v%2?>99R< zbZ8dw+F*oNG}Ov_6PUysy6D6!E3W72uGJQvSE-XM?6svsyocXlcI~#gYr5iobR!z&6Y@lDC3Zge@M9_f&QFNqb9KD`MqL(P1ry~rq=<{nY z(dFh>=)2pC=&`h$w8Nwdn*Fkx_8zFEH@Y{{VUjyEe{Czh;%PhmeV>r_EEm%jl@xs; z4$+yWcWJ|Oopjeo7d>bCJ^I1uJ^B}QpMDtEP2cG3rUNW{Xz}G9S}(PSZm{d2U#;w+ z_b%+AYY%tR4c+(YTU+kaEA#HrR&7#xQFj;p`)VhBD&;O6bPm&aCZZ=5P_)bj2_61e zL_>Wc%}f^1TZh`{J#$*=5SL~;{K##(C#jK+PN}Cy&eqb45^vGx!mH^19ToIWMj5R! zxrFwiZqUh1MReZrtMpNoEA+_2e7fh_B^o*8(w|>u(YnPMbn?=4`llq7es(6APF$Qs zAFoNE>wM#Bwn8lZJUN;+caEYphfmPksiU-+d>FkeV?SNg7EHf>v6C+0Zlkjj{OO3t ze)RcOYiNfX%jxe=A+77fr3;1`^cJo=y}M)?{mR{eKGSGV=drBmdvfOV(HJB8g4_c7 z%>y0UJ9`fO)^#Q={i8-KA=v1 zx=StCCZa6EnkkM_9o5iKP8D#jQ_(T`l>Xr?YRA!3%G@P^3YZ>E4U`|H#CwA&U!gxm zO<75i*ICrO#--GfWtP;I7xO7mvj&xSNRi3|Ba)Y*r;@hRcFA!@sYD|_U2@vvn51OJ za!FL8i6p)1yLkJV%i=CqBk`<}A57N{Q_io4S)BZ;cO15f4)^SSSMK`c&0GWBDDI`J zDcph`1>D_Js<~brZQLtbUEJC1e(t4~kKCUzW85zhay);5GEY1>jhE^-lgHRRmzQ3u z%X>SZ&&!$cSWhdo;MLYy^ZX9l^U5L@^VY>X@M>N;^6U>g^F9>2@Pe{ThM z^(F)Ey&bPO(OM@s`B!Dx{}>Kpqd{7nxYR(x8Q(0i^gJio)KnqaN8gv|uN#$=Ca6-C zNqW?(4GSj>;aJq#&3;tof)FaMDUwRBPo(AwGN_SRSE%&QWmMv~M#^EPm`dNjrvs;(ru}?4hUE9;DAbKSn>F7e#Y?V`%TpGj!aeB-+_JjgAh>papkxXq_eb z^p*TVT0C~0{=Tb}KKiVZUd^keZJHYCd&^qrN?M4xG8bhn)InWk^HA%8`RG)F9-1?- z0GXK?pwg{|sOzi|QoU-7*55Kg1vgDmUa}c7-(imWZ7k5cUJImg)DoRmv_kert&r+7 zD^y}_jkW|>qkmD>$RNcUdBjywvKv5jhoa_{yXET1GEDrzp+B%d2x^}pV!4*6GQL{pI_kpVN<+sW#c@VXLE=O z#f5}>A|#6U>>yrMoFp1oCJ`FOxkQdqF){SBj`%o7M5q%z#ILxQ#Ok3h#6~+A@Oifq zsJS%_c&KOs+s*U9rv^hHuV)FyBNqadAtz9^#uK~{aKSG#5}Zo+1=A+^gKvS`L5V0B z6q+0YXCjV)@W4pmJM}cEyBP=mG7|x(_Z(QTI}KQlUI2Ar*`V*+MR0a&0eIF~2v*x) z2P@)lg7gPvz)GtM$wz!L#xa0ER}vf~h~jCap0rS@{>Bp8o=M`M<#()*mp0{(xO><6t;_94Ox% z2h%=|1Hrd(@aEw-@VGn@Xa8Ca{IH&Lq%#geWr-KH;w$T9)w!I(tRQG}%uO5RZ z{~iF7vHPI-K^M56gMpg+G23};5AUD zc^L%eUIgkIS>RUid7x5$4zQn{0Y0Orfx+lWaJMZS1cmPfd;V+%`fJyLs~3F0qY?)2 z%5((moz~#^GXy5*q6 zH3f*ftq7F`mZB+Zt5MRj1{5;A1)b>WK;yAAT3{hXryC!jT8Cb=CvO0KA9;=J7k@w+ zo4%kMp`)lN_&2gLlfhX@viQSWd2G8_3A-Oq#-giJ@DEBAUwW>FFTb6RWnRp{-|1O+ z#kJX(*sq0`*l6RKPj&FC!@BsGoE}a%z5xGwtB<=p4Y4@h2v2J>#-cYS*zUI}zW&Dy z`+PLVDzpXOo^FW``dH!d5i9&Y&Kh5vVuQ!SZLo5W4KA5yizlzN#hVV<;@i=-_}fWa z%m}i@avrw0Pu>=Ds%`LLz73xF#2SC*TjMj;R(O(v6~5(ZiI?xOz%bezx12G<9mh@a zg^ecI#L5^u4jE$0Gy{CoQXlWF(!;r?^Ko#(Jbd)IHjdDri%%|}gDdxIV!^4I*f?4P zJ04QUoHf&MrRh{$@?IIw$XCMuu@vx)XR`S0b{RZf{2M7<9Yx=FeL;?zAJCtxuMwd? zfW~8NnJzB{|DFR0nq;H9059uOPp*otZB!{jg4$)Bxo%HLF zDthhn3-s+9VYK&o585|fh5i|IjapOoMly5^bGJL^@eH6M(a&=x%vJUgAq}a-Azmr5 z_@Q4*VfzubBeI_2+=OcMZWMUpw$Y-~zs?@Iaj1YEWbw0G2BT1Cz>da3m)Z z7_(!6QePriy*?ETv}J(3>KDON#$_<*Uj!C!E&=__O7L907Np;91XsLUfJ&PH6wH=@ zF^GZZL`^~Z_kggaj{*A93q0eVgIh9#;QHz}AT@3XWM_Q@H)B788>_y7JL4l@DDo#z z8vF%f^~S+0r+8v%0g}edwdXBmJfh~ zj6M+G_ykaz4}tTl``}ktC-``h26IA0ppw}R>L)dW{|$2g|Np5I+=ABxx8TKkH$`ys z&x7^}Zh;Zl^xEs+ojdK0D)H?auLiuz?$vtj=DM}-`6ko;W%6#XaAmgF<33f-=wZV1 z@>?rU$@65-ht-F@k{|kcG0kqYt~9gstQ)=RwOTc0!i+S?b5N((^NVwWmuqIMm$p){ z*I$OLm_b{eAL^J9O5*bO+8H`_kd5mi} z3m6XDDi}Y?#f&xUqzpfqZbpLtE5=&y-waz)f%)Z{Dzlr=WOnYJ#~g!Z%*dprOrIm3 zOf_2`)9`^WGoWfGbM@Q9%(jhDOqaMb%*NfR%vt+#ndQ^2Gr!T5%;nRXm=_Yo%q_bg zFo(&P%vocfnPG+hn5S1Ou`E8Rv5Z2rSUb1svqT3hSkoMruv(vbu;%+g*5}lwY-1SC_Kc^A%RGDd9$T)tn%9t8*AT zP(6x0&nKQe{mX<%Yb28`>tDd0LKL$NE34Vj;dj`DEn>E3!F~3ST0i@f@_V*&{tq_Z zCCmA~Mwv5qS$zUdpTkLAqsLj>WXfTlx8o!YIB^bUGdMK%<}5kBmZMy`m9sJ|gmd%m zQO?&#r#P20&Tu@-Q#gk6vN-Jid`{-18=M)@m7G(%8ad9l1swi~yPV+4N1W|`FFDgr z4Rczoe{mk!$Z`KBDsz37sB>SK&*tt3*5yVh7;_W9T5*>+EaCP)cIDoD#^#nWNbbg2 ze%zQ9+qhH4_i~@iJHlPLDT@1ZR~&amZ4%d*m(ESp&*eHYuW(OS7IO!8RdFTZjogC1 zcCKj};uc@+;p%Jma}7Sdt9p60?Ss6%hd%Pe z`$u_Uf`2^aN_ir%P?`8oFpVJ9XA*m(wTPca^NG?QhQxv~bHc{hj`;7i17V`$N}PVd zAYA$hqEurAasKE!qF6tGV9M+wytVcbKEV-0*~yc{#*fj&sqlEh)F6p){hUfzzsw-o ze&!PMY%dc%aYe-TX{AI^Llv<-wx0NSq?x#Ose>5!KoiC5?h*dq9ufW>{X~27Yr^rt zN5al%gvc%VLp-pO1*|!WAWBaKO!Jxsbi*{jE0HE3muds8migfNT0?MQrWv?9Yz6l9 z*#q`F2aurQ0(fgZfw+zhyp})^+`SyM#QB234I6>N+O42Bco)z}3jqZ$4*}Bm7?}F) zB>2OM2B(W+!H8Z0n0X@!Om|KNf1joUnUpNB!S5nyvn&8r)2@PKxf>H}U@0)vssvp3 z8eo600mPKw0XyW|!N)x!5b=lxcuf~rGTaURBuwm07xn?|Z!f@z=neRg{sHJ6{Q^vm zj)ImGzkwZ727bRU3sH~)T%|M_PQR`K^Omc^pa0b0&vFe|v1=AwZ#f5^`ZO2HRnCJo z|LMWRl?L#au`y&!GK1V_7VucTHB2e6gAFl@VBmHKSh>atF5$YsQcHJOGsg?goxy^O z6}V92H38rCc*D3SBs_U%CH$7L2G%C8gWFGUfSWh^!$rUIuPK9mE8^>UzJdp5kZH5&>CvS1~V1>ZJj!cMJB z=;WILJCiO<%t1(p%b%Qw8wb+h>Yh}%TaW^aE+)fO{^#JJRuasm5}`OC0WKbnhYs80 zAdO=n+wwFV3pxeQltjXd6W`-o|1miHA{@Scdk{W-845iQ?tyaCcS3ITR_H3@4+qx! zK@*qNP}7-&E&BlUc4Wg@$sW)<$aw4ZuBV&bFkS)SZ|8_HDCj5gkwMAs*Y=={)0=yCacq&mULee_v} zJckaW_x)JNjgrWGq8xAjhpx8 z;xA^Gu$fUl)^fOvkB?lzZK+qW&$c4mGUWz-wY?alZ#S{*ky0E{REA$aD92xIDzT?e z6@Ia>8dvMz!r$#{@UI!Qcx`?y4tJ}=>XYlSLrOhf{JI_&1UBH*iUwSAzX21*jaXw% zBR=5Vi0up-akph7{%>X@-d5Kz!BsWjspbthU09DbE$Z>=t9AICMjfu+Sc{{s*WiX% zx9}^2TR41PHTKl1!q+@1aj$I!Hd83aX7@_*o9q&-zvCuWGb_f?E!VM#KOu#B0=&664~H->Ve=!oxF$avH%`sMIV&=7Tu?gh%1*<1ohevx_#Ebpp2hcO zCSs2L8SG~ohx<)p@McCd-n}#m@0@lLhiM(huFMGRIr9+CxfF{3Mef1+`8#oo=Qb?v z_Qy;6*W75+8SStGRJNf z##nB-J}%g)i~kwS#lJN)aTnCUC7fwE#6blcgeqcjjVwNSVjMAQN6=X2FtR=T1|2W$ zN9%}3=t_nZ(XkYYc5g@42bxe?*Dds3wG5Hli;&TaeDv=~HYy~}qepX-kpBBv)H4`~ za^^-LgMwhRci9%Sc6JgGIQ=~F**c5xbGbx(w7f#ZpDZG}FWw{;_LdWO zURD#&73&By`$nSV${j-UXe*K2-%h0Y3yFsnVj{YnA}+=vqQ&bj@%D7*1XIyPI9uH# zycXOie2VT9=R&%Pp-zPV3iFP^wQ z?H&>ESV|o4=^|1qI*Hk5?-F=FCjQewgkBRxJoc6lwoXlUCv@ zznNIL7y&`7jhswcG1))F?cw}|^uRRm*g1;KYLBl2`gh~E!x5bF6w1UPn;(AT~~ z9DSHi{LZ^XO!CYnp7&)DNrf4Nl5;vi-AyGv#U>Ny7A6s!DiVnC)$v54N-Uv~8%>N_ zMG;?LpCF{KjuM))!U%=q`-yEhn0WVMClRr98}Tc{pSb+Vk9cOkhB$wDIbrn`5{i5- zVb#kZo^Eg_K3-i$)Z04{53xO=xW<}re`7{8#TXHDXD=Y`cIgn0&decliJ8RjpK8R_ z_9?{L8%jicq8t%h_m|gSJ<7Am80LYKuX&>KKAvd!0nh%;UEY$VBHqWcX5JHxI-W&j zId9zIIxoI0pQj#}#rwM}mG|8*fu}1Q&1+wCl;<@P%ws5S=Dj|#f_LsLix)H9k;l@r z`iu++%5E%G`xwFEHz7&{PJDA*}2rlLS3rur9P!tZbCi2Y(Y(W zZ%t_~u%}x57E|Xa2TE+`L^)VGQ@@g2sFKsJ)R{^*D%;MTI;ZbWJqpGo|hCXJ0US%wo)>?*Fdtr=cPC%D?(iHMMor) zwvc0c9CIRP8giFk-psvy?i}}=Qw4YCru$r#*Q4CPAF8}u8$I4kgN3|F|5&`Kr~P>M z4MTWw4bJXP74w_M3N`kSD5U zO(AMDHHf~YT0}QfkC=VIgqZf;nm8A*gt#Z;O5}-|geVXa(-^CWOS}!l8{cij%B_2d z(XxZY)3?Wn#oAGXrf&@K>f#wf?m-ffXOu>K{4ayp*quXY+U679mkNodf7c0SU@38^ zuafBZsU>)IjYQVM7UFTSfVekHLd-ma2|3w&gz3%)1pEFIVwqPzu{-A_kv#H-_(Xmn zlJY(i8hzi1Z41VTw-YM0>5Vc#i;@M?AIO8_QbllK0*#MuQ2}1XQ^BJ6X<&4_I!IZl z0p|Uj0e)A{0!ce&gQ=6X0Po6NAhFc}XUgY+QvLZLDp3zK^)CS0W(Giat0DN5WCRXg zH3k=NnSf2jrr=<*8Thxu9JJY5fZnGTpy!AsSfgkKrX95c3!hm5XB%tKve_CqL|FsZ zRBJFj-Wnuqu?7(n=U4Yz0mlPYpkmY#JlSjsmS493V}0h}rHnb4|IZZIzcc}UHOAm+ zpb_{sW&my;)dwwodf=J%d_Z&Ofd}igK{jh*Ua!t8U9uew75mvfw1bw@n$kwbN z=JSh)7{(>SaO^x`dn$p@8H^+vUWO9V@(qNNA(v>Mw}@!IqDw3rolLZZzT%Zx)bYAD z#PcLC9C!=UD!7`nUW-cAq>`>Z1oiap1xl4iQSlAq)TGCA=-ty6(t-p?9|+h%KdCrL zv)3ollu<5iDOXHejMdSbv_$l5ZV#<|?j^nP{TG^LB{PBPDxv@CC-@*0P4slrJhZ;i z5LxM3A`rC@rG0cly1t(1Vlx+AGbd4QrY~AM$sbj2+m5ou!3Y>0Lh?tBp`AM-k#y>5 z^n9Y{B#V`ZW<5KH_U}$Z4x<;4UU)WI_2nYkvaJ9ONej^gw}2F4Z=xv=%g|=6D)gAF zL5}~`BkE!k8Y^u^`)k`!Rz(LY&l02l12mHJx{F-KyO3`2eRPfU04d*_Sj*W@(eH|9 z$Wrk+LaY~Pf~ZvKGg8x14ry-#Q;_zNok^c8LJ`;N%+5#&2G zfjO)EL|^8Nq5r1*LdOPvp<9=Jqe12$^d0{}s_x@R@4`4be0Lnhz8^<@U&qn=$Kz=K zm2p(HY8*{|_Xq9U@CUgy{ziGfexc=aejy{BF|_-|^y=dy|$7t)n2dMhjeYEmn7rL8` zQH%$L7GWVWaBD|O7n&#HX+##EY7z5i6%xEELtVASsQ$n;^ls*5G=A|S5^7{26nq}d zDLaRq-$-;3JEwxYLd*CF2vJ}9%4f&OMVq6<5$(SA(>w4z`( zDwsAE|!N&_kc3Idvyvlr&QtR=W4L_?Q}Th#SHkIo&`@_n+@ysX~9J{+A#X5 z4s;CDg@0uA;Njy7;G}o@u;0@VK94towr$2x{*4Ko{@WDlj+?=~AI;%?$^!m7ZwX_4 ztl*?kD>x_K8XlTr13SZQ;PY-9_+Xwb{Jp{!vJTn8Ez!2{)k$0E6=Vwu4_jC-ZwtGt zY~Z+$4NQMx4HZdixTo3*hRR#PPESi{zQ+P)M4Q7?XUyRCc}mcqr2qq-$-?^WGEhbQ8<<=j1;M+%0By|=VDh!sK-ORY_{BU03%Yy2MU^ga z%$foP%LHJz!yRC(R|md~RDjoY#lR}05L};m3Cz8b0R;9bU{dKB;H(o37XCX5W>oA2 z44(j?@^~e6 z8EDS*iztkF8FBoJP~WB!6wazdit4qf^L8UT>fM5twF%ItSrRnM2P3QPQgm}u4_dnP zF$(ax(@gIu4 zAcKwnn}pY_l*LQu$l++Q9R9OS9`F1jkL`CU;BTD@__MwuKC?v;|2n6Lx7<|3xLOgv zyRL{06BV)SMn#;et%#FE3b@H%0pEQik58|kIGW{fv$`Bk@{z?l$0lJSQwC>q{-LwW z#?k#LzfoY*PgKYoMFq9r&`YCc3p ztM8-6uufDENu%Mt6Me1BcGNzp8U6of0sc1;wy2aZ8d4F~9#lVDTj2Mr_V;IbN6n*! zOG#JxhFH6Gj)!WYljjQEhL8>HhUkNfoED!lbh@v2q9KSqyKee{!%l;YT@Cxne>xT! z#W*gGuWk6zmgf}OV()a5C)b#>_OYWcG3NBD|51bde5ZOHa-WklQQUBSqt3Dw@3ZQ) z@5VNsTgG43azt@iMoD4AiKw`SfXt+Z5HiiF@8AWe!i@QieJT}BwwI^X=X8E+xW~V? z?9R*TRKI3+@$)~3H8@kK03?tUHXY`DnK<*ami<8ICT+drmeHmN%6xE}rz>~c0g z?9R#lD(C#y(_E(>b#ZP2%6IghD>uE3QgLNGOLtOP?{vpz(|eZ(^VYfkCv*AEzWKJx z!Xi~%Lrp*2K6$3@_JYnCZj0u@J7y2On^g^_xp}4Pxt*IU-|TlZcO>gI;}>x zNxDm#(Feok(BE&}RAfz?&)Usxscmj|^Jvj_zt=mZ<@7z&ygW_K{rfJ@JIy0R^P*va z8-2T|`O4&<&2VI{d%jS?y=cB(%f$4G7M{!z_l^6LTvcW#wk&18ah>+E+3 zS!%r zQ=V(>ScY4({jL^`^hvFY3-a85-}3jE+WECvW#en2{m<3_MYfg)yyev7IO}J z?5Y^_cr_`x`Bh6%YhJ^NR(238?wW6IxxcZy+TUAv?^l(2WSw2+b#3okkG^fI+vVn#wJ0248=A_xOybc#b@~?#{Sysx3Zq*yEJFs;BY4Q*A;Wm9`bVW?nPcO&;GfV0-Y3 z^cML{ORxH^%iAuLYqneJo%Rfn)OhkEO4^POUTKqkd(F#xK;Ogc*uD0W<_HfBQzgdn z?d-O$|EvWRqr|iOS*}-&Pk`3}<4et_Z@y~}EYxY=MJ`|{#kP9RpJpnEZKT>q?U{_^ zd{2gCoDfu0O(^iC`V5n7e@0_+fS^%qDA>UdXUzJ#->U}i5Q)BEbmomaieLG(6O=narw`6ilHwf17nF9CsRvprr62=tsMN9=-Qn2aa z2EkNjLp#e=>Ltt4X#aGxld&sGy<;rKs$>0zDCVN#1I%K^m5wdS3PJ_R64sf+kC~GO zxkA3sK;ZkRg1K@v!Dt*?+Hu=`YRCKDdZxQiREPCvj}Gx;6&4ld*CD~0Lfe?REL)Ac z%pk6b&@V~S@yK!p>wdT_bJ8s@;bh1A%)p;HESK%eg(iU=Lbq?JtW*;&OMa`FFl~t% z^MUAxmv6+K4*cR5)85vBS+m?h*u@H9#(-}A*R zbq^!KwZK8)+%u+v;cM4ejF&%|iWx_RZo#1)(oLqU_xY(z_SY-I)b%@U$H_2-4m z`)zB4M8R93x3?C1un=>4@0rsp2^^)t3YrxbNj<75w((CaO$iK-QzO zeVfkag=w?L^NmHp*ZqWzBiig&VY69J8IB_3B~e1}z5sTU_ zUMys-^$~&01+0#x1iN^vy3oW%U9?JPyYTznt?Z?ZzXg%U4Mnfj;#m%R@>qw){xa{p zc_lnl_^0FFl})U*fo4qE$ZNvP(ZlS7Z><77EvaZpttorQ++Oy&+M6Qdx>@4)im%uo zmGn8wSoI>wNqKS8UoCd=B3aJX@TUTUP&IMq2PU&_u@c9^ZL1JE2t`(Lz8!Hd)7VdM zZerOEM~VV1)kI$(yl2^GIplcI>)^Q^hO0-?P|#ikyRZNc2xIS9Iye8TOJH5_bHswrFU4k9f+tWcI79 zSk9v`b@3eymiSDjEN4tXo|FFrhz#yDi>?)$v&CFp(TLw!p(xCVqcnT3NQ+@3M$=!i zqtul+cQQ1^tNZIF2Ho6cw{<9TMC+=0#b&o3@YsnDK{lqlk{2hwoXMWQ} zK`SVsnZalF8=E@jW93JpS&IwVBQh;4+h?al8adm9*WHB79h=M9^0Bh4D@MmfS1vsh z8OOHYk!E8MK6b zeBY?3$nLc8M3pM1ruA7zz$RsJQ}R|eKBL5;+&V-7$=^gp2Wr`EF(23qnc1RY!AFrp z%SHC~rbFz6;ZV_bceXH2DTp1!JI4(8)+w@bXk&6hrR*sO4hzxFS)z7kACrg&!uT~2 z!t`~RU0rsQacOL!cy+@ZwuwlIV;kQrjB4}~%?L|nZ8JH}T-ItM@|$1A>hY~&>#~YP zCSEs14IN_kJLfv~)fY{o@3OU`ZGq3&Tm0s*yB}4HuADIy>d<_)xyxOa-~15KBXUUC zaK)D`u2yEJT$?NU^6rT6Y)uTy;M4^cl~5-1O+G2iJ0Q>A@pCPcn7vq3+}OZsm+G*4 z|I-rL7a55DUejXV7-O>+#j1(+f7&QweN<;>_?}@8^4TKErPHF!$$%YF>c&p~gN3d~ zDkdySC`R$1h49_LNftjXjind)UFd9lykp-THTIY33p&ao3`7}BK5MPzdzLqN*I`*? zDioYqAS@wp2N{o9Z@0)ZXW7IE3vV(y{;z>E4W~kD!>|k~38zTnQ%FK7V{_Siy=!M` z)L=-7hGU2Z>ZHjK8iWW%hRR7JTZT$h_IlTDl2Ftk6%u9Y6zN3z`uAMV@B6x*``5EH zXQ^o1jGf$U=@K44cPb2*l@rA|7IXX8W>V)yB8pQC7pdxtrS@V?src?xQQtHp(b|O( zywB|ee>X>$`!A9exn_J}MZvr29NiJpS@JU6$I4EWKYug-`aV?p&c<1EXwX{}lTjws zEq%)a0}4f#UmnKp3%jJ(wU&#vJXH~uHRVVL)Xwlg%PLXQ-6iG!!;4alRv-TJtfT1h zysy;7rGXzxSCo$4YQv?a_M#H|Ug|G(lD3ssb0sTF(djdLsgK9+>|u;5Cz>Fo$d9kQ z)Fylp&_fq|=G-XDlM@h+8jhBU7uOo}5t0&;g&rLjLk^)_qI!9#YUBTP;DM_EW z2GaAjMbzA1nZJ$MfiwFPd5Kvmd>L6!SIm@Q(xGNvx?n2R&2XkcE!VNb>@{ldOQfku zN!;oj<-TW*QX1t*wJtp4W3K-PUpG95F+_%s*z^#ln`F{TmFGNDvkcuHoux~UC{q{H zE4-@pB1~Ud$kbd?$StuSdt&Z{l4wOZnAB ziL~XbKG^lj@_#bcOYS995q@wn-;g?$uCMc@B@@H>q>VGU`;w!yPdy7BFB;23eNI7Y zXDOTP_6sYQE75y1M^mFRU!Ktrhm}u9(!BEp)!Ys$-?-`g|K#-*}kbz7mIb zmZyO2%uuqptDbdt9O46^CuzH)8)RNy@&WE?0w3UC-2<(tN-S+!= zu=16qaNO){$Z8LT8@(^^pWl7H^+&ABiNvhac~JhwjGmqr&2?9dgJE0(@3O?KKR+M8 zX>FkHbItf8yQ7%Yt40fb24Hhx1ou;EgrS~V(jUGR(N16T>$kXCss6=J5jk?OL6>5y zER;=l;$QY^QgSkr`9usu!&qJTx$X`miGgQb+zxpsD&Wm08~mX3ngJtO9$C|0MX;6H!9hmp}Iao?>uFyW#z zD+m&kbHy*Pa`ty{d+N&+!wLl&3u7@SyBWr*)}r6oTVQA~hLt#m2z7gU@rJ!Ro%T`- zv*RPn4el6Yv-55UG0w+f?Hu;3jUWkf+BnuL86U{>01x2o`!gGfzIZIN6Kjx7xAvj! zg-Xz_@e|m5&q5u&yOI@^IZ(cFCVB}A$j8b;@?lU7dnbjl^tpH7RLerF?+?QV2hKqL z8#iJy7>dTt?f{PJWUP-b+Ffraxv?>j)f9r-HEqmE>nRlN9ED#NtViMy56W+A${j|SeOGq3xeO*ShP5vXJ ztTSO&a5Gw27D|fDT8N<70zQ?5qgLW^FiCg7wP|&v@^C#f{`3hpCuU&k1_Kxrk~pYWZ5OGPZ45@m5|nkB;i{>72Y2-iPoi$ z1G)G=c;EJ|APTq&PIM5b4BIDoY%`sHJRrkQnusyXY%-navQH9RKa4Nl)By@JnuU5t zbhxKW4?Nkz;BA!+w_er34C|e#`La^H^d$C5?RlG7= z8_k0CY4$3H(|r<|S!^Oqj0$2;7Bynh88iBN)o#4#7Kx?Dh2@deD)h@jJ$|yX7={aE zu}xl6^4iZ6y-&;2f{tToYTwSTWeAyaO$9LBvlvWHfY|&32Gll?8YeZ}vA7*BuZu(j zGXwlj@F%M_bc5C_TToC{3>Wm=ups>!6YFh}dEWC0k!G>pcBmtqn`{K&y2{dsKlrzq&`x5J{~0f|RR9t32J#mqD@cIw{-IWoX1 zY&r-t*@q(YaHx}?in&wGC3Z_QNNpD*VIHBF=3@vS?rgv%K{B+ZC>cjuE8&t=yR+1KKzbFn-Dk&(7LU)(pNR71f);bMq)Zc~=0K z7>5ax^_AGJBRN=kJs2LVj>cjmIeZuxB>c1`7DkzdV;U=Exv`GIh!9nBW3egL1$uzi z{6xI09mCc+ZZDUwl}GsI1dUmq=o@5ZjWaIq7N$mLf>*sE?zEnYvqSu`WH1eO|D1}(JqaXZ-%&jIRT)ajF?9dq z0t8P`CFkc93;nHc@Qd1$Fjgn_z? zAZ$`FmIOY9kt=Rvaa}1>ayG@eRy}b4nmy}oUyeB=dRW7(WOAr=3yv(PfD~nU?0U&@ zw$5m{$rR{NOEHGV{SaoKnuZIeKejd9n>@23YB&V7g8_BUK3mY!QtfD@Vjq359eMS?m0d0_mjg#25X0KLW9+{k?_ zHE1wpm8qV<#7fu{yP#Z6HHJhgrh|CGVtkdR4{yznu+^Rg&@|^23+wV?0Z-!D`Fcwb zTwDPU_Krr2gM@{rjv)!%;V`g%IBbqpgm;hY&~S=5yhzZ7itvYIM@biZxO$r~SM@Vn zAgfBvjAVG~B_YUGc}T8FbHMm$B}Nq5KtS$%{2SB2DW!-_%)Z20W}Bk*)dNtV;f`gU zQgB}ogSDNf8F7pjnmeqA=<-9@bv5pOzGoZz5IO+)WEg+fK0r+Cd&x CUDAPML: + return super()._construct_pml(pml_ID, thickness, CUDAPML) + def set_blocks_per_grid(self): """Set the blocks per grid size used for updating the electric and magnetic field arrays on a GPU. """ - self.bpg = (int(np.ceil(((self.nx + 1) * (self.ny + 1) * (self.nz + 1)) / self.tpb[0])), 1, 1) + self.bpg = ( + int(np.ceil(((self.nx + 1) * (self.ny + 1) * (self.nz + 1)) / self.tpb[0])), + 1, + 1, + ) def htod_geometry_arrays(self): """Initialise an array for cell edge IDs (ID) on compute device.""" diff --git a/gprMax/grid/fdtd_grid.py b/gprMax/grid/fdtd_grid.py index bd5105a2..82297f94 100644 --- a/gprMax/grid/fdtd_grid.py +++ b/gprMax/grid/fdtd_grid.py @@ -21,17 +21,19 @@ import itertools import logging import sys from collections import OrderedDict -from typing import Any, Iterable, List, Union +from typing import Any, Iterable, List, Tuple, Union import numpy as np from terminaltables import SingleTable from tqdm import tqdm +from typing_extensions import TypeVar from gprMax import config +from gprMax.cython.pml_build import pml_average_er_mr from gprMax.cython.yee_cell_build import build_electric_components, build_magnetic_components from gprMax.fractals import FractalVolume from gprMax.materials import ListMaterial, Material, PeplinskiSoil, RangeMaterial, process_materials -from gprMax.pml import CFS, PML, build_pml, print_pml_info +from gprMax.pml import CFS, PML, print_pml_info from gprMax.receivers import Rx from gprMax.snapshots import Snapshot from gprMax.sources import HertzianDipole, MagneticDipole, Source, TransmissionLine, VoltageSource @@ -152,7 +154,7 @@ class FDTDGrid: self.initialise_dispersive_update_coeff_array() self._build_materials() - def _build_pmls(self, build_pml_func=build_pml) -> None: + def _build_pmls(self) -> None: pbar = tqdm( total=sum(1 for value in self.pmls["thickness"].values() if value > 0), desc=f"Building PML boundaries [{self.name}]", @@ -162,12 +164,127 @@ class FDTDGrid: ) for pml_id, thickness in self.pmls["thickness"].items(): if thickness > 0: - # TODO: Consider making this a method of FDTDGrid (that - # can be overriden in MPIGrid) - build_pml_func(self, pml_id, thickness) + pml = self._construct_pml(pml_id, thickness) + averageer, averagemr = self._calculate_average_pml_material_properties(pml) + pml.calculate_update_coeffs(averageer, averagemr) + self.pmls["slabs"].append(pml) pbar.update() pbar.close() + PmlType = TypeVar("PmlType", bound=PML) + + def _construct_pml(self, pml_ID: str, thickness: int, pml_type: type[PmlType] = PML) -> PmlType: + """Builds instances of the PML and calculates the initial parameters and + coefficients including setting profile (based on underlying material + er and mr from solid array). + + Args: + G: FDTDGrid class describing a grid in a model. + pml_ID: string identifier of PML slab. + thickness: int with thickness of PML slab in cells. + """ + if pml_ID == "x0": + pml = pml_type( + self, + ID=pml_ID, + direction="xminus", + xs=0, + xf=thickness, + ys=0, + yf=self.ny, + zs=0, + zf=self.nz, + ) + elif pml_ID == "xmax": + pml = pml_type( + self, + ID=pml_ID, + direction="xplus", + xs=self.nx - thickness, + xf=self.nx, + ys=0, + yf=self.ny, + zs=0, + zf=self.nz, + ) + elif pml_ID == "y0": + pml = pml_type( + self, + ID=pml_ID, + direction="yminus", + xs=0, + xf=self.nx, + ys=0, + yf=thickness, + zs=0, + zf=self.nz, + ) + elif pml_ID == "ymax": + pml = pml_type( + self, + ID=pml_ID, + direction="yplus", + xs=0, + xf=self.nx, + ys=self.ny - thickness, + yf=self.ny, + zs=0, + zf=self.nz, + ) + elif pml_ID == "z0": + pml = pml_type( + self, + ID=pml_ID, + direction="zminus", + xs=0, + xf=self.nx, + ys=0, + yf=self.ny, + zs=0, + zf=thickness, + ) + elif pml_ID == "zmax": + pml = pml_type( + self, + ID=pml_ID, + direction="zplus", + xs=0, + xf=self.nx, + ys=0, + yf=self.ny, + zs=self.nz - thickness, + zf=self.nz, + ) + else: + raise ValueError(f"Unknown PML ID '{pml_ID}'") + + return pml + + def _calculate_average_pml_material_properties(self, pml: PML) -> Tuple[float, float]: + # Arrays to hold values of permittivity and permeability (avoids accessing + # Material class in Cython.) + ers = np.zeros(len(self.materials)) + mrs = np.zeros(len(self.materials)) + + for i, m in enumerate(self.materials): + ers[i] = m.er + mrs[i] = m.mr + + if pml.ID[0] == "x": + n1 = self.ny + n2 = self.nz + solid = self.solid[pml.xs, :, :] + elif pml.ID[0] == "y": + n1 = self.nx + n2 = self.nz + solid = self.solid[:, pml.ys, :] + elif pml.ID[0] == "z": + n1 = self.nx + n2 = self.ny + solid = self.solid[:, :, pml.zs] + + return pml_average_er_mr(n1, n2, config.get_model_config().ompthreads, solid, ers, mrs) + def _build_components(self) -> None: # Build the model, i.e. set the material properties (ID) for every edge # of every Yee cell diff --git a/gprMax/grid/mpi_grid.py b/gprMax/grid/mpi_grid.py index 8ac54f4b..2b4dda20 100644 --- a/gprMax/grid/mpi_grid.py +++ b/gprMax/grid/mpi_grid.py @@ -19,15 +19,17 @@ import itertools import logging from enum import IntEnum, unique -from typing import List, Optional, TypeVar, Union +from typing import List, Optional, Tuple, TypeVar, Union import numpy as np import numpy.typing as npt from mpi4py import MPI from numpy import ndarray +from gprMax import config +from gprMax.cython.pml_build import pml_sum_er_mr from gprMax.grid.fdtd_grid import FDTDGrid -from gprMax.pml import build_pml_mpi +from gprMax.pml import MPIPML, PML from gprMax.receivers import Rx from gprMax.sources import Source @@ -62,12 +64,13 @@ class MPIGrid(FDTDGrid): self.x_comm = comm.Sub([False, True, True]) self.y_comm = comm.Sub([True, False, True]) self.z_comm = comm.Sub([True, True, False]) + self.pml_comm = MPI.COMM_NULL self.mpi_tasks = np.array(self.comm.dims) self.lower_extent: npt.NDArray[np.intc] = np.zeros(3, dtype=int) self.upper_extent: npt.NDArray[np.intc] = np.zeros(3, dtype=int) - self.negative_halo_offset: npt.NDArray[np.intc] = np.zeros(3, dtype=int) + self.negative_halo_offset: npt.NDArray[np.bool_] = np.zeros(3, dtype=int) self.global_size: npt.NDArray[np.intc] = np.zeros(3, dtype=int) self.neighbours = np.full((3, 2), -1, dtype=int) @@ -262,8 +265,58 @@ class MPIGrid(FDTDGrid): self._halo_swap_array(self.Hy) self._halo_swap_array(self.Hz) - def _build_pmls(self): - super()._build_pmls(build_pml_func=build_pml_mpi) + def _construct_pml(self, pml_ID: str, thickness: int) -> MPIPML: + pml = super()._construct_pml(pml_ID, thickness, MPIPML) + if pml.ID[0] == "x": + pml.comm = self.x_comm + elif pml.ID[0] == "y": + pml.comm = self.y_comm + elif pml.ID[0] == "z": + pml.comm = self.z_comm + pml.global_comm = self.pml_comm + + return pml + + def _calculate_average_pml_material_properties(self, pml: MPIPML) -> Tuple[float, float]: + # Arrays to hold values of permittivity and permeability (avoids + # accessing Material class in Cython.) + ers = np.zeros(len(self.materials)) + mrs = np.zeros(len(self.materials)) + + for i, m in enumerate(self.materials): + ers[i] = m.er + mrs[i] = m.mr + + # Need to account for the negative halo (remove it) to avoid + # double counting. The solid array does not have a positive halo + # so we don't need to consider that. + if pml.ID[0] == "x": + o1 = self.negative_halo_offset[1] + o2 = self.negative_halo_offset[2] + n1 = self.ny - o1 + n2 = self.nz - o2 + solid = self.solid[pml.xs, o1:, o2:] + elif pml.ID[0] == "y": + o1 = self.negative_halo_offset[0] + o2 = self.negative_halo_offset[2] + n1 = self.nx - o1 + n2 = self.nz - o2 + solid = self.solid[o1:, pml.ys, o2:] + elif pml.ID[0] == "z": + o1 = self.negative_halo_offset[0] + o2 = self.negative_halo_offset[1] + n1 = self.nx - o1 + n2 = self.ny - o2 + solid = self.solid[o1:, o2:, pml.zs] + else: + raise ValueError(f"Unknown PML ID '{pml.ID}'") + + sumer, summr = pml_sum_er_mr(n1, n2, config.get_model_config().ompthreads, solid, ers, mrs) + n = pml.comm.allreduce(n1 * n2, MPI.SUM) + averageer = pml.comm.allreduce(sumer, MPI.SUM) / n + averagemr = pml.comm.allreduce(summr, MPI.SUM) / n + + return averageer, averagemr def build(self): if any(self.global_size + 1 < self.mpi_tasks): @@ -276,6 +329,20 @@ class MPIGrid(FDTDGrid): self.set_halo_map() self.scatter_grid() + # TODO: Check PML is not thicker than the grid size + + # Get PMLs present in this grid + pmls = [ + PML.boundaryIDs.index(key) for key, value in self.pmls["thickness"].items() if value > 0 + ] + if len(pmls) > 0: + # Use PML ID as the key to ensure rank 0 is always the same + # PML. This is needed to ensure the CFS sigma.max parameter + # is calculated using the first PML present. + self.pml_comm = self.comm.Split(0, pmls[0]) + else: + self.pml_comm = self.comm.Split(MPI.UNDEFINED) + super().build() def has_neighbour(self, dim: Dim, dir: Dir) -> bool: diff --git a/gprMax/grid/opencl_grid.py b/gprMax/grid/opencl_grid.py index 93a80e88..8e68d6c8 100644 --- a/gprMax/grid/opencl_grid.py +++ b/gprMax/grid/opencl_grid.py @@ -19,6 +19,7 @@ from importlib import import_module from gprMax.grid.fdtd_grid import FDTDGrid +from gprMax.pml import PML, OpenCLPML class OpenCLGrid(FDTDGrid): @@ -29,6 +30,9 @@ class OpenCLGrid(FDTDGrid): self.clarray = import_module("pyopencl.array") + def _construct_pml(self, pml_ID: str, thickness: int) -> OpenCLPML: + return super()._construct_pml(pml_ID, thickness, OpenCLPML) + def htod_geometry_arrays(self, queue): """Initialise an array for cell edge IDs (ID) on compute device. diff --git a/gprMax/pml.py b/gprMax/pml.py index f6e641b5..4e8b52f8 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -91,7 +91,7 @@ class CFS: self.kappa = CFSParameter(ID="kappa", scalingprofile="constant", min=1, max=1) self.sigma = CFSParameter(ID="sigma", scalingprofile="quartic", min=0, max=None) - def calculate_sigmamax(self, d, er, mr, G): + def calculate_sigmamax(self, d, er, mr): """Calculates an optimum value for sigma max based on underlying material properties. @@ -205,7 +205,7 @@ class PML: # x-axis, y-axis, or z-axis directions = ["xminus", "yminus", "zminus", "xplus", "yplus", "zplus"] - def __init__(self, G, ID=None, direction=None, xs=0, xf=0, ys=0, yf=0, zs=0, zf=0): + def __init__(self, G, ID: str, direction: str, xs=0, xf=0, ys=0, yf=0, zs=0, zf=0): """ Args: G: FDTDGrid class describing a grid in a model. @@ -345,7 +345,7 @@ class PML: for x, cfs in enumerate(self.CFS): if not cfs.sigma.max: - cfs.calculate_sigmamax(self.d, er, mr, self.G) + cfs.calculate_sigmamax(self.d, er, mr) Ealpha, Halpha = cfs.calculate_values(self.thickness, cfs.alpha) Ekappa, Hkappa = cfs.calculate_values(self.thickness, cfs.kappa) Esigma, Hsigma = cfs.calculate_values(self.thickness, cfs.sigma) @@ -682,6 +682,39 @@ class OpenCLPML(PML): event.wait() +class MPIPML(PML): + comm: MPI.Cartcomm + global_comm: MPI.Comm + + COORDINATOR_RANK = 0 + + def calculate_update_coeffs(self, er: float, mr: float): + """Calculates electric and magnetic update coefficients for the PML. + + Args: + er: float of average permittivity of underlying material + mr: float of average permeability of underlying material + """ + for cfs in self.CFS: + if not cfs.sigma.max: + if self.global_comm.rank == self.COORDINATOR_RANK: + cfs.calculate_sigmamax(self.d, er, mr) + buffer = np.array([cfs.sigma.max]) + else: + buffer = np.empty(1) + + # Needs to be non-blocking because some ranks will + # contain multiple PMLs, but the material properties for + # a PML cannot be calculated until all ranks have + # completed that stage. Therefore a blocking broadcast + # would wait for ranks that are stuck calculating the + # material properties of the PML. + self.global_comm.Ibcast(buffer, self.COORDINATOR_RANK).Wait() + cfs.sigma.max = buffer[0] + + super().calculate_update_coeffs(er, mr) + + def print_pml_info(G): """Prints information about PMLs.