From 566aa5ee96ac41c411187bca748ef056da3767df Mon Sep 17 00:00:00 2001 From: Saveliy Date: Tue, 21 Nov 2023 22:07:22 +0300 Subject: [PATCH] remove Cell_number & add OpenMP --- PlotScript/convergence.py | 25 ++- PlotScript/src/Release/sample.exe | Bin 38912 -> 47616 bytes include/FDTD.h | 16 +- include/test_FDTD.h | 144 ++------------- samples/sample.cpp | 2 +- sln/FDTD/FDTD.vcxproj | 4 + sln/FDTD_test/FDTD_test.vcxproj | 4 + sln/gtest/gtest.vcxproj | 4 + sln/sample/sample.vcxproj | 4 + src/FDTD.cpp | 289 ++++++++++++++++++++++-------- src/test_FDTD.cpp | 133 +++++++++++++- 11 files changed, 392 insertions(+), 233 deletions(-) diff --git a/PlotScript/convergence.py b/PlotScript/convergence.py index 302bde8..983e714 100644 --- a/PlotScript/convergence.py +++ b/PlotScript/convergence.py @@ -54,12 +54,16 @@ def execute_cpp(field_1, field_2, field_to_plot): numbers = [float(line.strip()) for line in file] convergences = [] - for n in range(0, 10): + nums = [] + for n in range(0, 5): with open("Source.txt", "w") as file: - tmp_n_0 = numbers[0] * (1.5 ** n) + mult_1 = 2 + mult_2 = 4 + + tmp_n_0 = numbers[0] * (mult_1 ** n) file.write(str(tmp_n_0) + "\n") - tmp_n_1 = numbers[1] * (1.5 ** n) + tmp_n_1 = numbers[1] * (mult_1 ** n) file.write(str(tmp_n_1) + "\n") file.write(str(numbers[2]) + "\n") @@ -67,9 +71,20 @@ def execute_cpp(field_1, field_2, field_to_plot): file.write(str(numbers[4]) + "\n") file.write(str(numbers[5]) + "\n") - tmp_n_6 = numbers[6] / (1.5 ** n) + tmp_n_6 = numbers[6] / (mult_2 ** n) file.write(str(tmp_n_6) + "\n") file.write(str(numbers[7]) + "\n") - convergences.append(execute_cpp("Ex", "Bz", "Ex")) \ No newline at end of file + convergences.append(float(execute_cpp("Ex", "Bz", "Ex"))) + convers = [] + for n in range(0, 4): + convers.append(convergences[n] / convergences[n+1]) + nums.append(n) + print(convers) + plt.plot(nums, convers) + plt.xlabel('n') + plt.ylabel('E/mult') + plt.title("Plot") + plt.grid(True) + plt.show() diff --git a/PlotScript/src/Release/sample.exe b/PlotScript/src/Release/sample.exe index 5759498b28161f8ab3360a9bd1c5d51e7f570ad7..7fda1d43de17638ae4a2afdda89e3e5886eecab2 100644 GIT binary patch literal 47616 zcmeIb4SW>U)jxidY{CKwv(YRKin3^sC<1|?1_QE7c41aF7$IO#G>OT^q=Y0*HV_ms zbctpi7SpONeX4C#pj54*_2B`v?IuAIUW9-z4ZhS=Z98dfn_8;@Tle=p_s(o4A-ve< zd48Y&=lAcz%(<`Uo_p>&=bn4#&LsQRCKkgOGvNw{8QTL$KRdsl;}*r(mXXoX-<&`KHqIB?zm(zus((to0HlGLnePePkVFrf0=h*_r z4&x^WR;F0wfG5Sg{1!*9W4YrN+5btJM>!=ce@%ZEnN{xg{bBs13H8K2cXMKX~`RZ8-g16PuhJ`r3rRMq)n8SIJrnRNxO4rh8(tT zu0)~!uuqo5M|v)_CDKr$8YSGB`)jWt6ju#{!mXi*4=Hbitq$O6Qp@t7_~?;aIUKewrzE#xN+Ggu zR6G-9x#56sgsen`b1Y4=vR?+eQrW*R&5`f8!?7~|jy4)h{5%0?N-`hMrF=k_X+w&k z_)u0(dX$kE+w3VP-7Vo5S#g-%Q(EMf{uuWj8e{w>xHCi2rW~-)6D*CA6}QDb<$&CB z5`~Y+$NuP6Lh_VWIkQ80;0lzHgC-FCcX#G%vXs^5&g_#Oka*Quf2XdK43xSzenkeR z`bu|j@nm=A5$XOQHONb4z98JY}l!XwjZ^+7_?9A@kSkIK!-f`K=K0x(G*m^m* z*i1#Tr|fepZ>rZJ37nLOu`^`l7K=RPnC$e?$YE?3IRJ91{HEY0tdgIQ7XT;8i0a!{U3!*~B1 zV0-TK!PtKv(3z9m%qC zRQCVVQoTr4E|rxx)OTQu<>0cY1&wpAo4?vl{-hh^n$^}E1hz!A8C3q2ShD_{pHQV!`rp3d?++SHpz6PJ520SQ&$cl%P< zGw2*d`NWH(NS)-ELl~D;j7hsY(RA1v1$)Q?j0f5Xn%Rx_yJ#4K) zIfzbd+4dO??QEN zE~-)qFu+bF`l*DAx&&2!w>A?zpw?iAt_NlaU)VQSfz|k7VV`D z-`>R2+xrS=@>WexG0IUE(o5qO0QLf{pSMOjozTgSw?h5$^I5|&w-Ug5DI*5PB@%e^51Aiy}+X|A>T>|h)H#aZ3Nl|^hM2E zEF3~v86CE|gtv0xS?CgTmaZRBa_|ScNT`BJ0sITOD8~d+uYj1#flU|XP?F;306EU! zw<+jEjt5X(`z5l2N-4r&>#1VY8L%CGL^WO!Ft*2~d{BFLaMoXR0(r*K4H;qE9W)mKl2!=JZB`4V0h%5YSDI-BUH zL08(vxEshf>U(f6F@ew=RsZuV9gPJFM=3?aqZw#l4(_E+iTb~b)HiGOFF$wv*bhpm zt&E2SWh)tYCHU?WQm5;pRGh8WR}B)W+V5=Ubngx3-Vy8DUZuT0i2^BJ^?Dy_mO`D zCv9CRs6HWWB9W|L_bn!E+QB;%^)gr%CeOAmoRbVMN7@wN88Bi&C&nze)JG_69!!{d z9wcRit<85b1|3TEt5Lqm{{HCN75lTk#dcw~1sbg0&|vyu^*C|Lzpvo7Hp-bi+?i0H{G625`ty55w`^)w9yAs9N6Y_IkA2h*E9lwNMXf(ic# z&Zr^xAl2tA(DE;VB5XRFeeCm-Fjm%=`{?Aaq%K^FNjuH{t$v#4Z^!`!ODmXcaIZ&4_kjNBGwN}{;70FOV2x4@*F<=sY6i9SG00EsNVKB;z$QpW?}1((bG*K--U3ipySlZ$XrJjd{9pzoXvL-WBC4$FD04A z$L>%<9_2W91%=ZoHAF=*=2BXePB&UhkO?amD_FRzXaQV~66%G8R1WZ^e@=sqSs%MU zbW0`0#SJPf2XkrVax`qc32dT;JBnptS({~}2#V8aXW1!&Nk5Bn^CANmd%=@LZ59Ts z2DTv=MQHK%>g{|j)}b!_l=#rk_l>>}LcoVODvIqF-fxuMCCYkW2%w(h24>?ihS}H= zhVHcy=koZz4Mx|V-v*~~NQ!0vh<~OhwN}H{*+W|F#WL52)hfq0rEici&hseyHDmlCNoQtX?FCfO4k3$@XQNT2y*73a zpEMM~88XgwWbl8xfSB=E%1-Lg#(nhw`E7>zS%G0yPSBumZwwI#-)wWubPSQtEbVrL zGh2ySSvg_gfh>m2=^;OMKArvW9E_^O*(o<1_thZVS0VfFDP^^iUD^%)UJSA^kJZA` z14VEt#>(Pb_7Jhusl0~0iLiAMAN5+S7I`vhWLns| zOr(w_S0(WzX(VZzvJhH3x_6}99z!EX%Eb4yNDA7s%)m;Tt?ONe(Mua_pKNZ$Zbl#6 zSZ0!ow3jGC5H1r0)pwtI?-!ixHG`@83B2-#P;CLW|DZk#t7K)d1zUDLQlk_U6=U&z zF#=qT7G}r+sL?{MMmPD$B&Eq>D|6U5lk!4QK_cm`Od)7qE-5;hX|l#GKu8@#(ttHa$pLE&1VgRyCbDzzscgss zd!T2q%9seT(*y>H0CiAXA4Ep#gm0<`OM4Onhw%7(HwI16y^fmkj-WkE-tkQsxkNcQ zKML-)v^%=z>S$o_|EeBrWwoO?2|dY3(}hN}xm-##nSsrPjF701=@r5vw84)pmfG=Q zYc?EsSOdPrD{X3r372+v_52DDTnf0laH(>;t4~6sjjpb+b-bwHR4~DcDI;5zm1TW$ z=D`}wpbk)r8bQ>=UJW+gU0u#blX<@r%Z%P6>>hXEp@JirgXLd#m-rd2L;gTEg0dOO2gLcu=5Q|9}7!g zdkZf(8`Xibwoq?^FfZ%yOzHF}E;T!|HCvjG64? zdj5%36HOrbgmVC+^vmSln+Kw@aywX(-9h`vFhooP|kCv71 z*CU4KJX{gp-4I7LdVeeJzCF|%%ZGz=l{UO)(QcOl65-8hoB!kaKDFZRH@JIiMue?# ztjCmXfsc;u=^G{<`@m@ zNNG1V4`|a&=G(a>-D%;R45gpu5p;!CE_yyap2CK&ji!W1I=3X^2DieDHNtk@?c=iK zGA{m5gwKb|43CC_*+A|C2NDGwKcyZCe2#q#(Nk`}5dzZfb9! zOjKa5^aGGL;f=x2z!~r?Tet{>1M2l<~M-xFB9Q&wW5wo(z-V~94ng4uRg&=rjZP&742 zg;m?%!(k$?ArxtwTD0MLModTg@v!w*KnDLXrsHMkIfUiE5Qm)YWVO?Q=m3$RBwU3q zt8aZMyzHR$nLmJd>n~7pFL3~V8#(4K2E8e}ZxvoiEWQaM1cFrO+eWZHPO<`yKwT zVx=vAfj)D7`qqmHCYZey9#<_a{<1yLpIv;Cb*2<(78D*Q+}7YqGCnkb&;NyHbGPJt zABA^@trCQ*0dd#+`(4}Ee%BAgK!>ewW1CfzieAi%x>V4Ul!d)l3@QB6fx^Eu7S`pH z%Ls%*Skh?hH=sh3L4;x-z>q^Y&mhx8`8BHG5L}{f(-v@$lt-oz9MmcZ4jfccxq+vJ zn!TW>eZqBfkYfN#qI#51@jMVt!F?A(m$LXnnBS$e~xqLt_YQbt~ii~Ht%KLcgN>20VD_x*^)B<6tMaAX9(IXm47H~^-fB)+bsr3^iI)cLyZ zCeBwo4HfuG`v}V4Jc>CRY$U$s>3r>>x#{MH0~@Bwnc(aJS@QIHl-JbVe4(;1PoDC= zC+6>0rebon%SEpdXHmiC5gs+6OTs{K*6<#7lKf2RqHd+hQ~sur1{aZ&5Fv5Uia=K2 z6XluZKoxpfXIq9d9^yzqaHpB}QZ|0l83U)}%Bl|O%DO#H=}P-f+~lpe}JeM@FHm1SWQ9@a5XpHKECmGcz2YympqKX?KV(-nkH~zG$y-lNaL$dQ>U)T7vJ&J>hx75L>xA<`g9zs% z&8^hf;ef~vDYd2#w|lx>QfgKwZi_;=E!~CNZOyn9HQ~0Z9*&A`f^NrNV=gTYXe1J- zS*(*m+~>|{OB~!i<@S1M0#3xd=S#vog|)<|QeuqkZ-Z>zUf+9>SYBW*gXOMMYu^v! z5RkW#Y8lY(ru|a`THh`Ygt-p8U01spMl**(*=8I^%B6iU7dDm?vG4GyOR{7qf9jMH zJt+m0k!?eUC#94!mZefgZb}uv%8nAL0uuQ95F&SUP|7i+wCE|F9BUfK%D0qDiR?V0 zFh>psX%N+?NAO{~$&B>~xZtD}^G%kp^&&)gI5}(`L*XtEoCcYjO==lvpkqgfjs65f z7ij3$*(Ghh8k#Hmj#l$dhTAhwN}ESeA|_REC$X4h0DB)= z7l_Ca>KZt6y|HR|5?BR|5Q6mr`a?dUl=`Lr%jgXWe5|C@lhmD4s8?&EELVEip?okq zdBZx+`*NCWbe@t7ph3afG<1|p3$gYslh#pso`K?aPLt>@{21Y4c8=r>4S}KP=gR#K zDYq4k=eIef)UI$F?;t$MCXjImMr;P((tek_0c6EqfISOa+Yl!SlNxT}5)&(JT0t@! zv?joMquMDHP`CeDw|)2UT0tag{_;CF!e554YBEWORJkmFfWRQROOSz336i9d+=U@y zN&`uukh%9r=9-8%rFWz%3B_mCnmV-By9JAv zPGw?|@lekOxFy`8je$9;HN9BZJ8nQ)9z`aXmdyyfxK3BCEB*>`qlwsuoj=ze%wZN< zXIzAJ2Gk@wl{=rb+rl?$Tz&9DgSR)MShKnaQ!QU(e1L*5yU##U`5NPc{jTv?VRs= z-OB4~+BDkc>dD5QtoqJWq{B#T1A8%Q%}`#o^Bro=I_}Paq=zz5O#K-ZR)50pPw@K> z`TbG45B~A};ekW1#+?m3&dOGP3rF6g^l&HIskEwxp;=9O%pJ_OsT47$*wvjwQVP`X z4NeKx6Zsp5B=SC$QzHDMO(_(GgeNxiBj8ZRIs%V|5OOICOs+t~J_Pv^TD_iX$auV_ zapi7qLnvC0GU_Tlnq8pgN5Xb>X(Vh@7iwYfd=7%ZhBc^m2#8!{JhQ|ZUs->qdG39# zOELV1)vBvOqp~T5)^X^P=uL-+ou*64Ik$ZyY5<|4yE?&0}2@p71FgGHI<(FM|RrBxAR?A{P%c3osJ%bLh z=SWx5uD5a$kCFoplPCBsvn&1M?7$V-ff6$=i{daz&qqnOd=4A@Qp#3pqe#^kFgAzs z4|>SCgOd4XFl_Jy9IllfD8_)v!C8Mo4wmuS;bxpJ#<^M?Znh&=J#`hbs5(779dx+a ziU%rmp7x0PHkDDjdSl$mdvb8jrJRa-!U>jhYOt@LxhaK`Y9%}uOFynQW>EGQy_d)> zzl*_A5=pgM(kNO6VkwEokz})Np3FYq&2D83iNPGVvH*wV21~)&&N*k($*FIFK(Nal zxYQk3V|E8NSn%Wyt8Y8D17tcr9!&ZtXahlQL-tK?!@emtCPMx6n4GM_k@=7(Y+S0B z04EN%w|lV3lY<3IhFe)>mr^lM`bvb(va45o0K=s=<33p1BhM4)*|(eT@6W@TvSfM4 z+3n%Zl#Qo8J$|LwdDl+=(E|@G>*=`S6e5xE(<3wR5swK^V0Tv9w28VX9-H7W4WtMU zC4EICa&}oMHGAkeibuKDLwg7w1$zMtaAf{3{1EqBa?7b0u6miBYp{`fL?wn%VHi^7 z*I3l6yRM}KRCNbU(>xf$9ITcF2Vt2Q|F-Nzwx}lfE42aNoMET5;-{jEk1x7zvMQ3iX2=7AwIOmqs~sHaqXZd z)nhP_I1k?4EeGWOUMzm195k2^1c=-rRYFN@BZcI^2pkwU%Q5`fIrc>s^23$+ZR!fl zLXlzQ`_JMi2;YS#t@7ZMERpnv^uYaGo<>q%wIli9phr3AKGs8TDuC+I?%?=DcV?Fx z+f*dmnGi#*7uv!7;RgRI8vJIQbsr#<2;_dT*e-$?}35e<7E3!#jc$# z7qxnVu`RE{F;O0=ugZd9is=l0NXM1&NXAAtJKYL!+gOh-5NsyckS_HoORsrKCxD)Wvz+Xl~+WUv!fk z{M0a>R*7EmzEYY^%>YtQ(IYAKpuvd~0W?*(GT)XS`ZEGDc0z6pUM<8!`8nk&N{C9N zzZ-Of!G|hksndl`Ae{{krA>O?Pbx21(dsGcz?m57`BsPY{2MN%OYox=9C8Fr7ob}F zQSe4u8Q;!%0&b{OdA^hve+2i$_#>i4H0R6=ttZzYetH;e)7ajkH7EODGfU4ea%#=W zMV-7QykFzM)2U%vKm%`M%7}lPdb?#3Tw7PNJgoEgX$n>3oaveC!Y0;U^Ky3YubwQoH)0?9ZhNkhPZCPcUYPd zA(EuwLI>8+pzd`uA1~+Yq!gSCPTpzcfItH|`(2suNn3sa14Pq7j=Ca^1~VIXdfBC* ziROkp+K4WdA?93Z6rn_X`jk?|)4;01(}Q*%bC7i*b8?Ld&oKR`qNFYV1AHTk=!QmNcmjP|~DcL-zq&QjdnomNas>UYHSnjEXun z#D<_31_Sa17Ac-0Y#-7_q%-p^>7m~PiAT9D6|a|Z%OMmPM=XrpgFb7rA+V#{@i&Ox z4pLV#lIhenO%jEXOcyvEnFltkQXVGphmb?+Ay9&Q8YWK>gbl*q;)oF@hgxlQgs~d| z!$XM%%RnLy0W>s4V%{d2a%LV_|0tOf|9&&|a!Wc{jL)#3YEo6Y3ubbWpw$i+3I4Jf zmkZclGQkiDib6Crx!8-UoO@0PDdgf!>7j+(NN{5dJAqSRA>t=FP#&eG^v2+vx5;K; zY4(+8%2!RKa;AAuLfsB3M@c4}^B@xzNuF}V@e1%k(=q1=H5Xkb?bdk;(HY>H5T(5M zE6F%sDOz{Hwb=3~(3GH?WQVw!&vF7Qx9_m|S;age*xkm>#eUpZI`{Vzo5|c~3EmTn zxk6JqnHv~-*#EN=g!#Y&dw>Z6KoI53d`H^68l%8RrQJ--+wvz(|LvsivYLw(Sd%K_ z&fj>$!Va=c`p*R&`d^REkp9yFfdW!}JnvcVSO$-x6qH4Klt1E`jeRGcA7X1QgoW1h zko4_iailHuP^c73r@CZMpg{Ems=7Uax-NPWgqqiUd_DOdb$bw=8=|mgzxmP`&JYgW zhLwXJ|m}R8`j>Z}c{x3;?JW5HbvSlAa&frp- z1@7mw07hUzF$sN43&}i1ort!>{KvGAJ3_mY7R-2LFOEPqU@zvw*1m?|9a|y}2>6R6 zY&uEL-%ncQ6iq41Q!slhfHksgq2?ShP@ok@JVjkxyI|$ta47E}2`WXc!I&qyF$1H? z+{68RkWIg)rgQi8W0oAC4z&(ux#YC%9RyN48atHAaJbaFLN4!DV75hQkd^*lVF-ts7ge*^m2zrQ)6~Yf8 z-W)oeH;4K?Mc?DhSGI(JF~Sh?NHD!AEtvmx-Vw(@%HkB_bnbCNP9_Ij0}n=VHOy@# zenf*FqJ>nMkC2l_Nav4%Gep(CA8&7^3rr5>p!EC}5;a0WqRU$(DLTgG#`6|(s28X(s1sS3 z<;cP$#_xA=T#!}j01bXI#+f@kFg%?OZX8=+C8%%lZwO|Kb}*BIAk!PXyNCFYoNrIb zLp{np)T7Q!(a+5jV2yMIZyQkXcd%~KlRY?Y-lQphE7G~*b6Ik5k6PMt*_joe{KNlW z>VAYlr|CYe1T@`;;&I(y?8JEr*vsXmvht&y5D9LrAY5JUqGK>kc!Phj8L|Z9$aOBh zlN$iCj2)yb9W;A6-{Aruv5XKjBzcN{M7riw_<*BnVGhFwg>)X_k|->vrv1x>_7gVE zA195!P}hD4BcjQ=_8%ibbRM9lHR~ABHLPbW&Myg#P=3&Zm{ZCV3JB45F)=^+x8mS1 z_h-3h%iUq_{c=ZGR$diLN6dbW=nL%+)>EIPR4iXG4VF4?k9-Y5+b=4J_`l(Uw7%vZ za5R?5+}$~b+UNFwkERgc3?8ub&}LLLc)>n{7aUpdpx1-kVT-m}CX<$TU|_91)HTi( zOKlifa16Y=gIJrF#Qox=0l&Ce_{GweyUCrK1R9&vadbb6Z_ItM1}qGpCQUu3gZ%d8 z5f{%#E*{p_UasQzz(bC#tw~Pld9b2vLr6Gndfx%Hg&6Pkam-zbFl>T#WVyLJ`?ls! zaaW40y~Q*g{l-oRb{62sPLHhD9GPv>mV0@F7!uAoM{sFaZ%#OODPlD?s9`>Pbp+=} z!=H*QJ(tK4UD0fooERU zF4F9azTpH0^x8Sy%+11dG7FdlSJK#Wr2wOKX6Iq`1SZGd;MA5zm`yU{?+n0j``q8Url`8x8>tb_`L)M|*%ya!p=0gM^jN z*UwJp+nmm~jdb=)o9S4T&MQaI9VQ!`be%>a*Ll9cBbpCjpQo616gY(;N6aBionOIW z(3U_*hU^M1^6dj|I<-9yq7CYR9gV>}>;zSLC=DJaiLHd|^reJ4#zsOEjk{Py@U_Lc zmlAx7QoPSOoc(OY!kNqBaSX`0?dL-iPP3mu8p(k@$y$9+lI)1S-q5?r_gsn*JA?hq z!d_)$Utx&x#Ax%4N-^;a-lY7oZa~+O0sZ(=;#3#i&&Pn$TmI)ZppOEbwno+s=t{1D zuQY-I12(iqCjV%W^MAy6lHJUUhXGBDA7VJ`(dfW=dF;nIPcxo;1N6UXKgn_%e9S?% z^FfTGW;=OP=V?Cqf_&1N*Ln z4JS5hXA!C}`39^;>QM8!?7Yk;pXxOWY8u{xa$5g2^Jy_*GQ#FCu~L}O=@Aop z=pvG*U*Y~WO(^zbmh%0WJd6bG$Izw<9xdW@U5PZkUAjLno)6M@XkP*wGFa=E$vCZx zXHzfG+&5@{M)?*VlGNdB@;Wo^)3_Cya;}+;XPOazQ!euFL$!zvp_PJ#Z&!Jgr^w!b zeV8;?diG-`SwTDEFcg@lC9&`o=2$M`yJUQ0LZAILyPHSPLTCtR4o!Mfo2B+?%>otQ zHCZ5Ld}`wivz3&$45TD3(;moRw&z~mSZElT|44s`dcI{tdpcQIZy%DL=Pss8IfOh- zo`qlO!Sb8;-mv_B+bpabbr_BzH)?zQH;IDEg?GqHxKGIg58_u`S$qrb`Te2 z&==uiM%bVU7M$GZ!A>r1xIpkdhfzI`o{EC|!rr4}o_$CaKJN&=cii5G`0yStwtq+k zu`CCG)P>L++BzklcRW`Dl0*1M3aPVrB<954xTX$(gG66!o9fTfBDXC3TPBtq84!#R_-|lD#2%1`kNsdYaKIV75lYg+c}A>+l)zeuB_9nV z?b~S2&NUop#~lcN7(ma`4KTXw7TRQ!Q8t*aIlI00zysn5F>`zcPx3xzKeqarx|y_ z12U?l%i7l8(4G%5@si~mz!^e&Tw(?pc&{%bM^j$nQPSrKtZR5$YU zz&BEpZTW*MEjjpu&-BGXWIkYM~v1|DD51jo+Tz!Pg=;C0rVX5{~HW+Pt= zSBM1E3DQrS0pW|C(ZCa1_+zO7w*1RR-V|9kMeIE3$$!tt6H{R1wG{{tS{~-M9WqUZ zH18B&F+4EtGuj@YO0Is;we8q(@3ZO6zi^ADFKg&=}~d)lN7@#Il^=U$B{9X@cohB zKjHU3<8C;k*^4hE9a6tYpcD1@DiOkeMVQe4L12>oUxs~g>TP%+eOR6GpExJjQ%PYw z$HpPdZ2lF8gKp)RifU{R6ox3^AKXaJELMG-)Lz!pufA%6Vu^IM>7X*>vR4i&Hm1SYr595pC*(lX=dDrrrdu+4rv3rO2OP3m;Kq=9pijTHU{g}c;e5#}WV$!Z5SPw-8c8ufvZz#U9@6+kUA;r3??8Y)gf zVSHz036HB+0t1E`uc{{;M7}zL(0Zu&Lb~l1b>>m(5rly`R{aa!ky0O{(zls8t~enm zdx49oUm`=34fT11ahQ|Dfn-H}0Oi#?k%&{PNbLC&Fi}~&-wWy0zCV^GI{R}-(k7hF zz826v^K*M6zLcZhiO7RpO2WZgu&;SCO4@We5*iNF-q;>X-=RhXeYTn)`RLH^An0l? znm`Bftq|n(oDc=tkb?r-sK8VdXpgN$ZE%kDn(e#sc^0(w;32gcDe8qFm%eZ9&tN`$ zJ-I{u6s>>{7ih1Z<4m3s6Smq3*IzKadRcr!wL{%6%0i$=tDlWRiM_bt{WA3}4wKb( zK)fN1L#~Qh{j3f0*`+3ck7!_-eG@+GqrQL=>BnD)Qm+LyN?Y$Zl;{ObD3MtY)4@>Y zVU(fMxcoI}oJj$5)AHNIc&q|VO(Ix;AoM#7k!S*(hgwa1TrS60^8u-b>w zXIo?JWW?2PqHAgcTEkb=>5KA7YMqFTOh2HmMzHsf57G&0l6)Ok{oGHMp?T#*?Oho1 zwck_^aH;d9;;o((jN;qKRiEe7LgRu7zoGb#1@+U=MtX@<8GDjM!ru|bi_fv!Q2@5` zIRJ_hMD)J^_Q!6-JqjNagPU%6!zc=pJ~(a+IaOVkO(3F0zmI@gy*tQvP37w#u4_EM+pcM?1 z{e}4aB)(S-J&^Zgd&7bG$u4L6fe;IV544%XAe`dbvhTV>Y^S$SruC}~H|^T=tBZ9ob& z+a7x#LS%>9V*^@f54gqTRUtAt8>wcyDT-PiNLLaP7b893G1>G)uqENMY%LKL^Q@pS zOjCn@Uk3Q4ee>K7PK2$8F!jRkCPyTle@!Krixb53jfTCHgTBs|{lzI6oS~(%f?N+w zFxUeBWIu7aiXy3WzIY{lA?86;3%C2j(d%D8dv^LZIe5=A`6r$D=h|ae5`C}CeDZ^- zKbV40W?QXK_Is1jN;Y0W*usf&CF3nTyqR5xm=cm_P`Frx@nV;VdnoR*p&Dk)ZkwEr zaJtJPr?<%bD|qxbaOBGoNgrqmqKD&+u_rvpg58tlFFFxzr6@-Jc6{Kf?J%F@LHHhK z=U?Xp8*hf4hg6)xgm(DwuBbM=?iQPkG`tL;L7J!N8!ocnoTAdVRuG5=PIxrwq1mx@ zXaQdxrSjKk^(GPYPqG3456xxEV>4c?G-MCQ{-^;|zBrA`rD26u;yC_{bgjpwS{fKgaXO2Fq!ET<*3k{4c zhzzQfL1RleDPSYTd77hc=-_SfkB>D_GeihabnMwmt#aAtU;a^F1hujvYQX<2z~>=y zAd+^#ToU)N^*P}5_v58V?W==!eEy+NdITK7cL=Ca7c~l5vxluS5rh*G4^w3h$}#FM zxzF=hN++0wudWc!@b*{_Jsu^w!4A|Uf?g7z;X(!NzN%`0aFCY!gH7gg>38J6vtO5Pj1XCs#)D#JKld!vZ;2g zveXv*cK50B-{Z?yYpi%H%l-p=zL@)m+`qvm{?x~y({cd*ZoXOlo(QXH;MySE!At&v zCyjD&t_sk50ewsg9f}KE8QCd6KGcx;L3NLCJR!$}zKs~Xq)pU)bx7TPq-kfUioPpTiA=o75W1=CcSOx<&Q?zjB<@`LEWFCN_%XtEbOf_kf0(cz5ItH$^RWf6i9dj39rCc zVu4Ezi@b|#)}h4sni`;w>XchPw8^7kIC8NBgXOnf)^>ec?JY{|k5Jm*VQ!D5cj8G; z1F;H1j?dd{z6^hdtv%LHNW4_?;(D_mES%5 zegnTZ^Ls15mmAW}JS^ym=iwu~-97yNAWzBV_p$t*%q-~fz*MJ zA%D_~gD}X!2W=Rt9_N$dB{x0((i2)7<8Pt!*bst0;&t=_{R{3f2OYbjUWSo`d7B5} zY?3k1;wNSq(K^Xbq&~i>hZ%NR9o|$gP?!9KtUG+$+hEC>bXw8A9avy{Ugk~>;yH-- zZ2M7!zEp}A>+}U5%%k{ck4N%t6i?CmW7laHpI-2HN5Q5f`?uFa>%M~Bx5*n@cA`Sq z`b}EQ@cB`vg`ma{e9hbfn%1CiUnBEMdP`!tyLtVsW+2A|X5+ZycG;0-?$<*GE;thMXW zVE;Cee^`TCHF&oM^E7y^2FGi#RqF@S;sOgSrPj{N8vIy;rswn$0yb&!=QOxYgT)%$ zufg{;Xw%?K4c@51|Inb0{q7<0#n00v&bBul)H0aadat(@H{29?IW+%s3 zhqiShZXW{m@^6-le0&Ii*VEH=7_a5;)63}@=LzeA0&mkL0(MUka9PEwipsSWHg8># zx5`&mS)t+UJ45-ZEyLrvv#{7!SY2JX&Q{@FYb*Cw+~q4JR5iZpvWmNGzRF5ld1b|2 ztp64{H;qeVF#mn&4@wrd1ud;ez9Fm&^{p`}gFh`4rG=AE&c#am|s@3b(zkK6)tUR{#Bp{+c5Ugqcs`{f{a4sBq2GP{4k*$YOa6|avwM|Bua=>e zivTW%<@l4dg+iF^uW9u(NY(VeJ)D4WJe<{6zk|y!5N43l{xAH@p}}}b%u)@yJK%*EH!{qbe-XcwFTDL1DSaU z*Bjz`4&l;qlQ2-a9T&$>5$)H6a{A9$Zs#yKsT|=98S;jtTE*rDr5VdL8QM0MC!9MCdFM-0 z9onHW(tm^}2UouVm!3xu>c%X}KuQ|IdAN2N(gt!FGvYV)V<`STpo7MOe#Ub3XTVD| zZ#CqdiH2_ArJos>5pS2FE-#?nh=x96-jMOs@lo-!;WPG!$`W7gxQ6ye$4B(hPcJu= zKB9}}C(5^A$=3|JFo&J-M{>Hykgpk~R_(W&?QTaWF{0>8yDiIF-XQZ3(8h^%i>qD5wwm?3Hlrfek zxlBVD{U;+lv~9{GTsm%J`-I>4)5}m=P6#rOeE4u@{LxtKGURt6d?p$Q3E?D~>1U*a z+A;o&baaLovpobGF{F&KJmJhh8U2@o@Nj$}Tsm%J`-FeEAC+w~r6l5Po5#a3gMxBp}K@i$8GEn z;Wz&DGL%+;Iw9o4*E{2n+V3*t_aQvA?^Ky^>9~#U6Mo}QFGFcNssHfF;PajFN9FSj z`K5+1RU#bv&)6>E)qlwdMAt{xZ$lc{&u0NGfK}_D1Gp@R{1Jc6B8TF-kRy%CO!Il` zd@03c2=a0aA(2LZ(iv%}V+YjJQqs~m9u;^*S{=(+&)98|GV6*Fgks}4A(U`hS#hDS zFy$Ibf{(~j0jH&<6_wlt+?xcw)3TQ2rsQQgs7&&ZG)lA$PV-Q?)WK;u@c?@>D9z3B zWeiE9MEj7mETqYU(>T7IL1{T0U*6y}UanwJS}reFIyjA&t3q0n3H`4oo}f|qI-46b8SQr2%SbYbwiqtBVmQ-=2JYE&cAh?=HR*EMy&zXgsX;6{WZxyzS0GX}9vW zyNqeMqFnc&G*M1P+AgB0PSCg&C)w6V>G2fWQ)fv(hK>MgV>^@5OwC9iNK2=*WTY8n zU>2p({?la;Jgx4vJS_y-L|~f67P7@0KNp>9Mm}*9^9Q?yTW))`)b3A3y$;nTRPo!`WK4}vF}qnWuliD7?$ z?ZQ|TlVtsK8dj;$%TzweZi42N6ByHT~am?yVV%Av~GV7|F zq_ZLEQLa(URXl>(bAiVPJhKFzs3bE>0^TIxO`0`=CC#){34WNt5AX#1h#C>YM$C^X zh)x53-bWL|{sudOYdvoh@?waOvtnFcj7}roNXCF*KAJ zH`dI?&Prfomt4?i=^hh`uNqxIV`=$5=J-o+yB60%T=Ao97*G54&2&%8*kD7UX>^M= znc{9hf#}3~GS1(L9rGc}jouyCW%7)=BAx^?YE(4IR!-vh(QN#zX>5Gd70l8%raQiK zRLER4q96{6LTzpWmee3_+i{u4q+zV4RpEB^jp*dWxbNoRty6pa-E4kzy~!STt9kZF z;1=?kGp^(kwyf(KmfdkJbM3i~&E1*7GPcfU>5bPz?=#s|rE}O7K0mvn_}lD?S@rCS zsJmHm-`MV?&I>{ov`_M3-pbf~=wxal0#00(QK4jr*tLB`N!mRcB~>qWg9nc2>uuxE z{)p?VP$p`$0P)GQqoc1k`QkFnd8776KgOOLHPRISgILRqq%rGaEQ}?ZP0T#gR7E3A zeXc_@sT*%wmZ2;gQZ~^R&upb8c5&_n?BZF8?Bef!qtDuXVJNZ6!e#Ig;7ffRx2HHgeHeha zWzo_1njit1zEnl@1^Q)kKwNP>OLUBlV`D!vu`wTy;xtmajB5|&fySY6W8Fl2`u!24;#&w#~K-2UwY^52yu<6*qv12gqrNh;uP`Q)fMYChn2(#%CSDN7E&!ecUjl6O zg{&zHBI^W>AEG_7J-|=#u#2(rdK%iMc+LXIDW!q#p>boKs4qsbFMOlfr^SivkGU3> z1b&c>6Zl85zJ-iUy%F+7xGpfU3&5*b@QTvb;@YaUF*I(h7d6&|K52L_GP8?N--4#my4 zL9_@b@f7V)JO@{u*3Qtlu})NU6eE3z8aIlKgT9R`hTX^=!Ny^pg8k@@0{=(yd3Yq| z9n3vTMpxxPreQPKxYdkZO5X8$0hb_LY>1B}`*9g)L3zUS!%D^`RiP}(6Xb2@0FI1e zBku$>LvE^tor_AKxqc)|z#PEke%8y(8nuU!-52v9e1ciAOuAq+Fd>a@`00?sITzC% z@(SN==np2Fby50{pqsdc`{jg0s=nSxMwHHu8-Er=>O&&MSdI-TM8=`vF5-LR8EEbt2LG;M=+NPe|;n5_Q`>?Q%LJWT7e-g{(0EnZ^H*Oq#b+#dWE#DhO|$RR*$qQLt6Zc&K~_Eyxq{DofMUFDT`7V^iZC}J#`<;Zb4mJ@MhtE=oYX4+Bd@#uMl1BM1td={)S0$11#qtM9BWu{)Q6silS0 z`M&DHGGC3|DR?RFZjvpGRDdMwYt~(+)eXhWyUXhXMd-xWYW!x`e8)V;@1z_>aOB@eQj~uqaq~uKR8B=%#n&c#l$v ze@$8lMIP#B#m=j%0dK!{XR?8SI^z(XSwMKe)_)m~#*}-FQc_k4Svs3+>GH87dR|ev zx3D_Gg>%*)wK%msYM_U0AX1YnxnQ?tp|HfEQ7Wp*x!w zF>`rZWXf@{mtdX9N#L$ldPD=yCDSwP4*c=0rOUH&X3R>PR$N}r*flJFO;P3QDtI^5 z)n(pl#%|D4N-C={BjdzbQNFjLm^p??uYucC?#(a5biuL=`K#FStediVD^3sg$}VU5 zKvkUY^H#4etAKZe9W9X(@#^x6Drpn^veVLwy^zJqb=cF$&#J4-Ev>Fx>oOke%%|i9 z)rG6Qa$!Yrxwjfmmgcb4-ql5fZ!U#aSFZ7*UkJe}W5ta&vhr)HyhUXt=v=9&#@NM3 zFRHDs_Es3YO*VeGR8U(%&q~qHCCe7)y0cvojv7<)X9zCOX5bM&^O;{&NTh-ha4*X$ zN-6=Wy}sIN>=-R$-a0T6PrnL%mDtWkEzF~uSCkLoEJ3X*)I#kV=3^yA<&`yF?BLf?AlrEK~}Q1x(vg@ zwlZX|L0jx`R$IZ5VG}1GvbQR~1gAW)-_P=ED-7`qqH4-2h^M3+I3`=l#px6aGvQ;ILYu&OHGx30>|M$&kr|22?|LSHGC0oYvtuF8s%vb$=lz4;Z8 z%gVL+h1GYh!RwyH$1-0TPS4&Q$?{gLDXXrmSdC%Bu8rWCR!XmJEOM>DL88d68wlrn zNlN(uSd!g$IRI`T0aa5xqediGdz)rbvvL6(Lc3F3Yd>-UZWy9{Wm(>BPOl za_W9QL?1)US#@V&5j{E4`zneG9VfiRC2%vJ69L&L8lb3xrdfmdP+7r|`U>x)@tI?2 zU(3<$7mmf|7gdA%oC@$6GJzU*F=B?$%fhO%sjF+It}UyWijq^SYb$(ZtG!dpXH1{=WO|cvx7>!Wm@7Z z&nm2f_Ia~3eSVM9Je8HJYO5A#R^L_OtA-`{h;o*fRr_iS%a>JvLB-e`i(OVB%tvv= zCgFfI*vZ&^u}i(aGiKYO7kYhJ!r0_gR~C8k;QB|gx%xQhNo;eByGDytRxj`t60JGa z-Wn42C8nEez18b-ywxz(n1hPY3yc>UpLjZUcQKn}NQE7RWFnFMjZ&9#{^nw+9X05H zmqM&EhSmW>*_L`B(K1<(22kYm-dTIsUEXS)l8uBTyRx{pTr6K2D}ju_@=qoEZpS8jubf@zLB^^>PV{;4Vl0Ig zM|yM~$r$^XvqMlH=nzRw5-$NaF*QIX+QaU|DUV@9l1%Qv^qVtQ*fFqM}9`R9} zIRLE0MP=yl(PLUT1UP=2h6nIMJTvbnoPbVT?`=hSz#rp^2cBlY_i=rO^lm`ccu|+m zS^cLL#yfn>bcw)Y0bHboa{v`BjQ1AVaV=~-TSn*5KEw6u_XJFzfM>G^6O?g%jBqDl zA1)KAa-7zpx`6swLW0Lr zFrG*!i1!cKBM1{bK1tNYJBIA(EAYG-wBbEN)_^m2

+k-qD%7A{_c{L^{E8T#q4) z_Ym10Tr@6thmozCf_IUD-*{G$&Mf`{*EXaRd<)kegxS?NQ-Sk^1>eWG0P1HA2}aLA z-+-Us2wbTMefEC$)APo3|9qmIVLx7h$1Z_4z{hXBX z92A|8dIr~FlqdKyE)`*d`*8Im+y_`UU+}6PaC(*q8_xq7&-572`_MU{ojAke`JsS& z5iUTO;77Qs5H_CGF`nb0^F2RXgtJ{JL-3cl_90C0H@MzJ*m#DC&NiJ`g0l(W1HnJw z>POgkzQ=e*$9R5+&h`B57VHr`CZH(~=a&#BcqJ|wVdME7;~5;|`5Zd8bM5Uo|BUhk z7vgG0m|!`sT?iY`anbp%KVipYJ?JDzhXHpY9RE$|FD{a+V!&oCOlLn|!$tDn0eJZ; zj5Bc3S;nVubs!u9{0**dgz;(ICDG%|D=Wc0L!tX^D*iY z+<=SPst4Swga>QfGeSGiS$?Lo|ljq%#FltU=O^tW;o#pFnKu}RLcQQ6wudgX8^{y_gnYy~HsJgPI zvcxwP8*A4W)~uekX2xV2Hd4w;u;;wo*cwpTY;%0owKYCBTtyArW#_<_ei^Tanl)Z* zESLG#X>p`fd%snSR=mYI)n#k2#JkH|6Uj7YyXsJt9zJ-yYrN(71_ADKCl}VZE7nx5 z@>Wl_)s{JmX!~gHoauvTnlnApIaHrBUFRDjrw_;JDfQ56 zq^eK#J=Omd+h*RDyv?>PW1D?j-M0E|joWr^+qG@aw$QeYZJpcLcJp@2c6ocw_Pp%{ z+e^1sZLiy2zrAsL)Ap^~o44=WzH9rQ?V;@*+dH@0cBJlT+|jgS>yG9fJ9q5bv1dnU zN5_uN9bG%Rcc?r1c4Yj-{*$Vo)cvIYCzhSbJDYdz*%{i|v9oh02nmVzW^$VAn|C&M vHupCtKVg5O>51+q(w?k)a?g|9Pg!GZq+{CQaX delta 15658 zcmeHue_WJR_WzxS83Y7o5C#DS1{@U?aTGNa#TglNP*74#G=HFwl0Qk&@aNzH>hv&% zUM;t5!z3g1lWewq%(WVSQ(2MKWP2pAY}+=U_qoqAI9mJpzCM3_|JmpD zy7!!W?z!ijd(OE(o@b!uTdAg1Y8)qZvfViooy)RH!ZPNJp5reAQ=Xn2T1bIBIizsN`F#03@6vqUvqzTn zY&Os9F<82=kU!AlS?Q6nyl+U7w0UNr z2uP6&`2Mhg(i8JZX66=Uc*xHW*Zg$&r}M?|p<=ick;fIu+R-P>vIw)RGaa>GnjG?nu%%``f39bVriAwhw@U{{^TO~%>GipMMfe`+ z>rs0mdP|9Od3?k@(w0&Dfr#0Xc@%#=;%@1kkvuWdA?+B+mqywp+erRjk@n<&X3`8d z!Gui)lsNLHg?zE_oX%z#Vj{ooO~-mutd%b*22=2EJ@C_+60o== zhl~iiCL&YLiU2D-=)eEp_4B^HBBiqyp57}*vzssPHCHo(U+iU;Qit;|dkvSS7xU!k zWa%vnw?+?>-k-rAh@QA=%x-K6e)SNZe_OluONR)p_oi% zYX1$TzjlScf2Qb{T{SfC`q!#@z^#chcuGuck9}~q*O0+;V&Z%JYZma%p?pEi2I=V; z{KJ?&y{6Obe=Ekx7qP-lavp4)5wdqEEwta=%}b3bW6Ln?hJL}BYhl?i7RwXEdGnQA!6*hbcg(~a1SiA_>7|p`J#lq()0rUUBVmE zQ&ai*K4y)BU+Xhi`sGw^?3*c7HSj5YXG%S$@+bR_)$HIO_1z~8ox&eT947sc&!0*> zpqb1kB+ZaM8_s#sd})P)-$*Kmy*U}K^|th54C`4g>~dsY_MS5H;$$w>Lwxe&s5;ty zX%ld`3*m~V?c8QE^}-&5O_Hc(4lOb&JD;yLW#r}(*{!R?IY(pC>^oOBqx_Q6Da~j6 z68#0Pl`Ys_zu7@7c?ZyJ%ySo7^W1q7a_^50{<&%JB=Et{L2ak#obeg|8*NAJmKL^d z72kzi>_Xo6|43141NG=%d_JG?WpKVkN2AZUERPqb%rb)@_wa!O37vg)ogSlmNnW_1w^}e=l&{=BA~S_~d@+JwaBy>9Fyw{Z?tL{HuQBC2KlQH*b=D zKaoFc9xff3$UiYxXeRKy{u87x2J`LxM@y|X{z`vYvxQF?kT(7?a*8Vt{w}oojOEyb z8ZX(+7Mgrs=4*MG9*28Y6FM7x#vRn@kijdWHEMP)KR#ek=y)Ysg%kN721JLW<#5e6 zAvUjU+cVH8`6lw-0|!n3IWQDEh^=;uiRtS-&;;jPARoyQH__oe_=!Ip5B`C|;Tq4k z3>;~qF~Z6)V2vlqM2GCQBDj=qt!4^&EfN)Go=bg}tuO4EP3MTO9hiohvLrL?oIKC;LevPU4jOWm{ z)&8i_*UyHvwLj|gb#usUxA7~QlpOIr7`#t^;Ggy;G0F|1nuTM>^1AdiDJ|#3bLnZC z&=*66=C}3yv!Rou561C<8Q=H67wq+(e0|<$7xVhEi+NGxttOvw?I=pHEfkXK4Sdot zQ}``H!3Mr&SU=4u{>ZSUQjf9RH*8Q(^eS%VBYCD}=)(6#V{zlYlrO*Ot?9*BrsrDW zcRrtYFmPNYjCIiKaJ!oj4sZKMn1b7QH|}kI(yPg|H~cG#>hDDTji|4S`p;1(UbQUM zghWE?Wv|C@B4b3BMw+7M6`85h06lNW9AP}!mwfcy^+d;+KL-`a7c#GSpYF@w%`B15 zMDn2{^LxI6vl`d3H@qt&`L>Z;rLQ8mW>lQiX9VvvYMs=P!1s*Wm+nHeXpLh>jdav{ z20CgxIvg%X*LCE>Hx?35|8W(L`!~D`!}-k7&q|#=`H!OqNCpcZaMwy_zH6poHeFC| z=nay7_E$^v!!F>s#SyOGegoHRkP}DYssY&F2$=7>HzeP+JWhW0JRFnnT5BqFm30=n zmU|0bYdZ>D%Y?Px`H%rw9Bu0z?GMIGq2vK-}NMD1@h_;_cTp`JiVeoUg~h_ zldWvY5`A(V_=J`nT{ObIK!a23T;|6%`J(>#E^JKrNFZj?VzokDu-JX|u z#jTBZl(oCHk&d!9u+pIy{vI8A;c(pSGz@FIk*hao&)}$hxKn7r(U_}`(w@;D-%-b$ zw{(&|p&~E)GkLN-PhapZX|?H(pAEXGh=(P1q2e+q`-PhG4td8;As3pU!%O%?eMz_n zbyT=V*uSHh<|JRVLx2{Ze3<44;co;vF{d_pvYE7i0J-3K&tv~vF?n#==o7ori(>moH6y)qqeZy$bWiOp+S$UZs+vUGg zjDHn`?dC}uh?swjwurAYB+oq+$Ifvw5HU_e$Q?fov52eIDgRD?{84g`T_$&*vjv?) zl=JFe!+fI^uiNzb=Yr0X8w>?Bxp@T_x)5p#tv{%_aQGz(YxpI&*enZ);1|hsORRg6 zuihh{wd?bt^`g!9y3HAMPAEAm{7>>>fbe+11*R&Y{5jtpTfLe;=j~z!cKK0CU5acu zsw=|&Hd;WB?gX9G2ozG0KS+s+o*=cjn^wnQMmLOd#{7PuXwjht-*Al<| zY*#=~DL1-jhG>y15V3QJ7<^DdPJ$}{C=l`v5A2{A)>90x%Eu|DHg~QDfwVcZk#N(l z+OvPJZ^%Z!vUtxUkak(5gkApEl199A)!)?b+WXV>i7Xn%5VN;)jE;=7sQhKzpUC@8Y#BwqWzS#(W`#`g{`l(d73l+;C6HE4nzT~*F|+z0v;G3qf@~)7BHX+aL)u#H=A4XU9zJG$mKyfuIYS zUpe;Uz{6a$uTxGbn++^-@VoDgc4A`KW&@T+_E8U6(oja5$x(_gB9!UNzNnwpmHjV$ zLmz+6A->{}2*!obE}syil5rLAiLK<_!0=A3u)Aj?9m^cz+$renJ}sQxQO$d*WyNyN7ss4^0Z9S5oKHt z$bz2>^QSiVGR+IU*oKVxt8M6on8T66;r6C#Fb~7vEyvgl6o)Wqo+hJ3y3!P zsXLMK>M>VDduCc-sFCA2V*V6r$~AEzD7jJ~PpyDUY(YqSAI1Dr9Q^+jH--y%kFTXA zc3a3rDt4z~tq1+~>@IynvAW}&xD|Zd6}qy2)Hj6sgOAW03fOSY>vw$mPi4dT$$qn$ z_M0l&Z&rx!SJC*G=j&h65APfOr-*u7KOVLqwYX>ra@vJ&T_yM?q?<@BeC;CUuTzQH zn?m(t0zV+?U!lHyx%>Oubv@mJKZp=tAEF6=hu19VBmJSBKe6C1nnqr~Fh}}kA-}ND zEd6N!zrJv2yd9To@30;SOPjZX@ZN-9^)~R~rRkcUd_ieq|G)GEy+DgoB;J ziTf?!_%fpNzRjCT^EGeszm{grTvmoxrGIGjb^Wjg>R($o<&thDNNC{Q!q@OWp*#x&NW z*ZVNU)qB<%(B4XIr)XDYjPP8~T%upCb0Nx(+BVw4KK#c~{+qJhG1E#(E-@=9P%}Qal+P_6 zl9+@J9E=~X>*39Ty_6U4!nhESy#x7Eg8*% z!h1LuyOS&PHJ|ZN0(j_K`5quX;{nlZI~A&B9$^Q*|1}cLRbURmD_7P_y)g%$u?{qZ zxyMTqnr{w*C1GUAvuKIa8xHu!35vsw+c#eEDE+rlTh@IIpWM5&|L*ulS-W8(cpEy*-yntuL(B# zXu;Q>&%)&$!TJ<6@M-*z@A>!3JM6($JDUj^TvID7d@Oq_lsC<2RO!iE@>YeEMbKyLQZZO$i^j zG&$*;#i087ABm&Tv}(@L4+39-c1}RPVQIgmkbfYR?6Nwf%pyn4eslDuu=I zPnM<(7=)#auLYg>q=z6m*dlzc|soA;}jAkFQ~%PTCJJ$zq9inO3N zKUt9+@%LEb7C&kHF_wQ`F;dd@<_Y%~OABIo<^7YanO+S@{ zh!=8)JMqcEBxd}6fa|^JcN3&Zu+O+Np5MH`U;J97cf}pO1@U~uvNh7MS^V&_v9XN^ zu&uUDV|VH5fnM)nmvx&L@|(+s56egUlJ0SQil+RqyCU`j)t~xHjJ|faLZ}7dw@7#9 zO59w`?^&LnF%6X3wvbD@@jzfObLp~y*>0XQRb|z-nJ($fXwuqr=~9b%>+-(NAP}ss z#DB~r&&@NZAsevY(@zzKh4k4y(=$f*13LBlJl$X6ef@u7q}{wu)Z0b9Q`FCkxKsw0AqKww-#bO)w5VSZb%p4dE9&BA3i>?)cZhn8sQZfg zYth>(YQ>-z1ioI>4pEO4oTrUQMqmFQ1<51oiK0Gl50>4AR0ddLLVU5k0e z=Oe1WUH39ywsMH1kLHiA%#>2k@yjcdtPe$td~rWiL_?}hbq`6S_Q;pJ+emqK?fgsK zZ4eVGqlWs%1A;`3IeN|So=npFu!r|8=iX{{!f&u2Q3^hy-Bu1=W!y+NpuAF3gdW_$+~6II&y!N_#Y0Ggw_ zZN0{BYnNW@?6}px?Tirf`wsr=)k)IJ9sHfuy@9h4p!(%3${`yZvUG|Tg3p+aj#ID`y>ECYd`;>Sx<2D*Rl9KdR6HKGM$ys{=-=Fh_-c!7+EVy?3;?KcR3|tI$90ufI~T z11h~mg+)*LITZcuJ(bZYScM&`L6g+BS%qP$qB^xbQEf|VyW~kR ze;oX(qV;OWdA0q#3J0hIOV$1vDm-1;WVd;M5=I{3P-AN zunJRE7^lL{I~1E{F{=jd^TVR~O1tt7==ZSsZl#}_*2jg+m`bRVjFEcD0o9*dsC`g^}{F3R+QcA7kwh)CJyl zvU1O6UE)T+SSSWwfzpnWwwbY_&0G_dG?9kTPeME* z1Sq>*5sXQwz<8XE-4E$3Lg}LMd$#C|o4UK01$YTc>n8L~3>+D#xhp^)jctL3pYwJN zNkPc8Ir{^2DDrRV^keOiLJz$f1Zuk~q;zc2IRk_+_%ZPD4<_2}g9Ev4HLJhHSx-X+l;ufvFlJK+ujp88D(P04FiT?0*8)HErR` zoC#so=#JZpzkefH+*G3*kB0wo0KZcn_0gc8ONh9wLZ$=V`!Td~E6HrC(HR3|wgXq> zV!8>$>QCN4>-H;>+R$ynfzTx+QKn_7#{-#B1pMdBqftL+p2Eo+GmB5>zC^{aT zKQ{$yjSy6XZV!0(2MEf!3WERX1b(M+)JFqoyp;P01{02^#*VA^;~E zG5eoDPmuus69+`9k!qSjqqTMn&;ZypAAf3zVnAbT4%R2|R(~RZMm+diO8!P|gws+e z@aih_%2gmUD|oF+zo3~d7BSw7@#60<&DE@;l(Boec{Y>-6@P%~!OF|HpZ?m?Jc>v- z!5F>FWLPX^OP1o+fi~HhkvpTSUdF5Q*A#JFNgD#;R@`JeNozx6$%xC^Ota zBVybgv|P|kf`$}RtjQMw({5ud5&|u^(*$4EZM0&+X9Z1oR?Mzg@D+j9ihV0%jrupJ zIG!5d0;O$UBkqF@+ZeOJB9*v)^@^&M$xguXAa$@}fty;TgJxG}7!vKCMAUK=2Dt zNQN$6E85JX{4H@trRv^xycK{iW6jz%WoFB8psci=zqozLu!?z%nd{LUgi1QXqV&N` zUmni%_(-c?)3cKL?dVD?L{(hRr|uZ)?3#w%7$QiE)*%u0Kf9@ zc4WkFCox)x!S~9M!K`T*KCGg!?L24ay-v{!KMiE;EtJh@#~URU{Jdl#<4eeRVVItk zCNZ`T<(L>lhTTbna{_qcV?$W%$OslYqgSV)BhnLI*|U_KWY~c}&rU=&dIuZACIbl% zGxb5>$2L)I$=X^3t_-Rqtu{CZ3{q50(|n%m8k0-{>Oe_zqk0-8B(ih?xSLR=KQ@g& z?pid6JQA{#v55$NS`-jFiXqGck93YUEd+NdD&p?{&{QKytOi5cHdMP%v|&@G@%o1b z_NSf~K*&0S>O9K3)Ei>J)TPTnB~4AoeE!Zu(L;&zOHfL`A@XPzkKeq4M@SZ;YhAd1 zMxcsDks^!c@hI6c6gVY7O*CWAqr{1JFw}&>uc1hwkbpkEN=_1X?gd-g6Q~ZMNRg&e zeo#)+xcKXGpXH#C#L+e^3zU;;+XS3w$Iy;S3S&tlqgc|KgwD8**!CDtuga)WLlK#o z_7KW&R3#{p;T|haA}ObdywysThpKj}Iv4VyhD2w$#KK2}SW;3WS?afpC4H`83GYcP z_IwcQg&BmDhcIn%Fl)_c>|>M_5Ss$~0^1rzM&-}SD2EJqL5rMi%*6*0?3Y}E{GilqvC3it~F?NF_TX;a~T@N2+-eI{e;iW!>* z{&mHC?8D16o&13<(YT+7G#~*`N_ooW-jbu0fBNtMsk)U1Hx7`RT6soehNcBw(R5kR zc~r(|sjZdoYRr&A+&p<>tgaMv#&s!ysBf+`0Rfk~YvBYoYuW zBlCi01z&prEwTw4256oDS_Wt}pwaH_f1W!CtP)q1|9mt1&o{IGd^7w1_swj~6UvhO zzkf43absUw=m}HPq~JX*b{tRGaUSF|S~4VnU$hJ$a6|!D6!3jZ5)l)g6=36$1U~QC zkwhW%m>=!K|MhI9LhZu~j|@@J90I3~2&sQOB6QHG^rPJZHXltVo+C#^=SN3JD=qLG z6Ff(bCh*P2i08zSWA}*vw%YStC|zVpxuuHI`Gh?!}p^Y9DO(*33i~^ffM`! zr5Jc8;2@)NBeMYFU!hnzVdX~qQe`{B|? zV*&k_Xo6SF2r1|UJ5lxmCwQg5GPc7B$64OT&2xI>8kv z2Y^=soD+D-1wldH|DbigAJWyudpk?t)Hs z6ik;ffZ!<8!OhzI3cEEBJEpX=wAWbM@FyNYdAP_>U0RN0q z2iyy|cQ%p|_+`Mqppc;O62@lFK{A3~3ivz<(a!@W<1Jn62Y~ZboZx;HCpcyvTg{GP zFn-`*Gf+Ih@goOYkJ1L5zV&>JQi4RIZ$9*ok%ZH0)n1erfSVRFb~g&)mQv-OOgP<_ zm!MQZP6gmOl+HSQUvJ`(FJ^1o8QX#%5~5F(zIafRRf)+l=s(nUsQu8@Lmh`Y4|N?f zG{-fUG?zA4H1BL~Y~J78)ZEhC+T79XZSHJNd)o4J*3(u#{AA)Xb~xm4#o@}s)rV^i z?>yXic>m$1!!3tf51%~jIox)*{czkfre}(tDS77VGo80sk6y&&~mWmVEaMSQx#9OJk|abJ5+RNKdfTSMX=z+!Ah& #include #include +#include namespace FDTD_Const { const double C = 3e10; } -class Cell_number -{ -private: - int border; - int current; - int start; -public: - explicit Cell_number(int max_num, int start_num = 0) : border(max_num), current(start_num), start(start_num) {} - int operator+ (int) const; - int operator- (int) const; - int operator* (); - Cell_number& operator++ (); - bool operator< (int); -}; - class Field { private: diff --git a/include/test_FDTD.h b/include/test_FDTD.h index afb7a71..bc8ddd5 100644 --- a/include/test_FDTD.h +++ b/include/test_FDTD.h @@ -5,12 +5,8 @@ #include #include #include - #include "FDTD.h" -//sin(2.0 * M_PI * (x - size_wave[0]) / (size_wave[1] - size_wave[0])) -//sin(2.0 * M_PI * (x - size_wave[0] - FDTD_Const::C * _t) / (size_wave[1] - size_wave[0]))) - enum class Axis {X, Y}; class Test_FDTD @@ -21,139 +17,21 @@ class Test_FDTD bool shifted; double max_abs_error = 0.0; - void initial_filling(FDTD& _test, Component field_1, Component field_2, double size_d, double size_wave[2], - std::function& _init_function) - { - if (axis == Axis::X) - { - double x = 0.0; - double x_b = 0.0; - for (int i = 0; i < _test.get_Ni(); x += size_d, ++i) - { - for (int j = 0; j < _test.get_Nj(); ++j) - { - _test.get_field(field_1)(i, j) = sign * _init_function(x, size_wave); - if (shifted) - { - x_b = x + size_d / 2.0; - } - else x_b = x; - _test.get_field(field_2)(i, j) = _init_function(x_b, size_wave); - } - } - } - else - { - double y = 0.0; - double y_b = 0.0; - for (int j = 0; j < _test.get_Nj(); y += size_d, ++j) - { - for (int i = 0; i < _test.get_Ni(); ++i) - { - _test.get_field(field_1)(i, j) = sign * _init_function(y, size_wave); - if (shifted) - { - y_b = y + size_d / 2.0; - } - else y_b = y; - _test.get_field(field_2)(i, j) = _init_function(y_b, size_wave); - } - } - } - } + void initial_filling(FDTD& _test, Component field_1, Component field_2, double size_d, double size_wave[2], + std::function& _init_function); - void start(FDTD& _test, double _t) - { - if (shifted) - { - _test.shifted_update_field(_t); - } - else _test.update_field(_t); - } - double get_shift(Component _field, double step) - { - if (static_cast(_field) > static_cast(Component::EZ)) - { - sign = 1.0; - return step / 2; - } - } + void start(FDTD& _test, double _t); - void get_max_abs_error(Field& this_field, Component field, double _size_d[2], double size_wave[2], - std::function& _true_function, double _t) - { - double shift = 0.0; - if (axis == Axis::X) - { - shift = get_shift(field, _size_d[0]); - double extr_n = 0.0; - double x = (shifted) ? shift : 0.0; - int j = 0; - for (int i = 0; i < this_field.get_Ni(); ++i, x += _size_d[0]) - { - double this_n = fabs(sign * this_field(i, j) - _true_function(x, _t, size_wave)); - if (this_n > extr_n) - extr_n = this_n; - } - max_abs_error = extr_n; - } - else - { - shift = get_shift(field, _size_d[1]); - double extr_n = 0.0; - double y = (shifted) ? shift : 0.0; - int i = 0; - for (int j = 0; j < this_field.get_Nj(); ++j, y += _size_d[1]) - { - double this_n = fabs(sign * this_field(i, j) - _true_function(y, _t, size_wave)); - if (this_n > extr_n) - extr_n = this_n; - } - max_abs_error = extr_n; - } - } + double get_shift(Component _field, double step); + + void get_max_abs_error(Field& this_field, Component field, double _size_d[2], double size_wave[2], + std::function& _true_function, double _t); public: Test_FDTD(FDTD& test, Component field_1, Component field_2, Component field_3, - double size_x[2], double size_y[2], double size_d[2], double t, - std::function& init_function, - std::function& true_function, bool _shifted) : shifted(_shifted) - { - if (field_1 == Component::EY && field_2 == Component::BZ || field_1 == Component::EZ && field_2 == Component::BX) - { - sign = 1.0; - } - else if (field_1 == Component::EX && field_2 == Component::BZ || field_1 == Component::EZ && field_2 == Component::BY) - { - sign = -1.0; - } - - if (field_1 == Component::EY && field_2 == Component::BZ || field_1 == Component::EZ && field_2 == Component::BY) - { - axis = Axis::X; - initial_filling(test, field_1, field_2, size_d[0], size_x, init_function); - - start(test, t); - - get_max_abs_error(test.get_field(field_3), field_3, size_d, size_x, true_function, t); - } - else if (field_1 == Component::EX && field_2 == Component::BZ || field_1 == Component::EZ && field_2 == Component::BX) - { - axis = Axis::Y; - initial_filling(test, field_1, field_2, size_d[1], size_y, init_function); - - start(test, t); - - get_max_abs_error(test.get_field(field_3), field_3, size_d, size_y, true_function, t); - } - else - { - throw std::exception("ERROR: invalid selected fields"); - } - } + double size_x[2], double size_y[2], double size_d[2], double t, + std::function& init_function, + std::function& true_function, bool _shifted); - double get_max_abs_error() - { - return max_abs_error; - } + double get_max_abs_error() { return max_abs_error; } }; diff --git a/samples/sample.cpp b/samples/sample.cpp index 17c45a3..11f347f 100644 --- a/samples/sample.cpp +++ b/samples/sample.cpp @@ -1,5 +1,5 @@ #define _USE_MATH_DEFINES -//#define __DEBUG__ +#define __DEBUG__ #include #include diff --git a/sln/FDTD/FDTD.vcxproj b/sln/FDTD/FDTD.vcxproj index 6c76308..eaa3556 100644 --- a/sln/FDTD/FDTD.vcxproj +++ b/sln/FDTD/FDTD.vcxproj @@ -92,6 +92,7 @@ Level3 %(PreprocessorDefinitions);WIN32;_WINDOWS;CMAKE_INTDIR="Debug" $(IntDir) + true %(PreprocessorDefinitions);WIN32;_DEBUG;_WINDOWS;CMAKE_INTDIR=\"Debug\" @@ -125,6 +126,7 @@ $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"Release\" @@ -158,6 +160,7 @@ $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"MinSizeRel\" @@ -190,6 +193,7 @@ Level3 %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR="RelWithDebInfo" $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"RelWithDebInfo\" diff --git a/sln/FDTD_test/FDTD_test.vcxproj b/sln/FDTD_test/FDTD_test.vcxproj index f8b42e6..75c37ec 100644 --- a/sln/FDTD_test/FDTD_test.vcxproj +++ b/sln/FDTD_test/FDTD_test.vcxproj @@ -100,6 +100,7 @@ Level3 %(PreprocessorDefinitions);WIN32;_WINDOWS;CMAKE_INTDIR="Debug" $(IntDir) + true %(PreprocessorDefinitions);WIN32;_DEBUG;_WINDOWS;CMAKE_INTDIR=\"Debug\" @@ -143,6 +144,7 @@ $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"Release\" @@ -186,6 +188,7 @@ $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"MinSizeRel\" @@ -228,6 +231,7 @@ Level3 %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR="RelWithDebInfo" $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"RelWithDebInfo\" diff --git a/sln/gtest/gtest.vcxproj b/sln/gtest/gtest.vcxproj index fbc2137..7c48291 100644 --- a/sln/gtest/gtest.vcxproj +++ b/sln/gtest/gtest.vcxproj @@ -92,6 +92,7 @@ Level3 %(PreprocessorDefinitions);WIN32;_WINDOWS;CMAKE_INTDIR="Debug" $(IntDir) + true %(PreprocessorDefinitions);WIN32;_DEBUG;_WINDOWS;CMAKE_INTDIR=\"Debug\" @@ -125,6 +126,7 @@ $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"Release\" @@ -158,6 +160,7 @@ $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"MinSizeRel\" @@ -190,6 +193,7 @@ Level3 %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR="RelWithDebInfo" $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"RelWithDebInfo\" diff --git a/sln/sample/sample.vcxproj b/sln/sample/sample.vcxproj index 3864acf..c9357d9 100644 --- a/sln/sample/sample.vcxproj +++ b/sln/sample/sample.vcxproj @@ -100,6 +100,7 @@ Level3 %(PreprocessorDefinitions);WIN32;_WINDOWS;CMAKE_INTDIR="Debug" $(IntDir) + true %(PreprocessorDefinitions);WIN32;_DEBUG;_WINDOWS;CMAKE_INTDIR=\"Debug\" @@ -143,6 +144,7 @@ $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"Release\" @@ -186,6 +188,7 @@ $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"MinSizeRel\" @@ -228,6 +231,7 @@ Level3 %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR="RelWithDebInfo" $(IntDir) + true %(PreprocessorDefinitions);WIN32;_WINDOWS;NDEBUG;CMAKE_INTDIR=\"RelWithDebInfo\" diff --git a/src/FDTD.cpp b/src/FDTD.cpp index 327133c..35251d8 100644 --- a/src/FDTD.cpp +++ b/src/FDTD.cpp @@ -1,49 +1,5 @@ #include "FDTD.h" -int Cell_number::operator+ (int number) const -{ - if (current + number >= border) - { - return -1 + number; - } - else - { - return current + number; - } -} -int Cell_number::operator- (int number) const -{ - if (current - number < 0) - { - return border - 1; - } - else - { - return current - number; - } -} - -int Cell_number::operator* () -{ - return current; -} - -Cell_number& Cell_number::operator++ () -{ - ++current; - if (current > border) - { - current = start; - } - return *this; -} - -bool Cell_number::operator< (int other) -{ - return current < other; -} - - Field::Field(const int _Ni = 1, const int _Nj = 1) : Ni(_Ni), Nj(_Nj) { int size = Ni * Nj; @@ -103,66 +59,243 @@ void FDTD::update_field(const double& time) { for (double t = 0; t < time; t += dt) { - for (Cell_number j(Nj); j < Nj; ++j) - { - for (Cell_number i(Ni); i < Ni; ++i) - { - Ex(*i, *j) += FDTD_Const::C * dt * (Bz(*i, j + 1) - Bz(*i, j - 1)) / (2.0 * dy); + int j; + int i; - Ey(*i, *j) -= FDTD_Const::C * dt * (Bz(i + 1, *j) - Bz(i - 1, *j)) / (2.0 * dx); + // Left boundary conditions + Ex(0, 0) += FDTD_Const::C * dt * (Bz(0, 1) - Bz(0, Nj - 1)) / (2.0 * dy); + Ez(0, 0) += FDTD_Const::C * dt * ((By(1, 0) - By(Ni - 1, 0)) / (2.0 * dx) - (Bx(0, 1) - Bx(0, Nj - 1)) / (2.0 * dy)); + for (i = 1; i < Ni - 1; ++i) + { + Ex(i, 0) += FDTD_Const::C * dt * (Bz(i, 1) - Bz(i, Nj - 1)) / (2.0 * dy); + Ez(i, 0) += FDTD_Const::C * dt * ((By(i + 1, 0) - By(i - 1, 0)) / (2.0 * dx) - (Bx(i, 1) - Bx(i, Nj - 1)) / (2.0 * dy)); + } + Ex(Ni - 1, 0) += FDTD_Const::C * dt * (Bz(Ni - 1, 1) - Bz(Ni - 1, Nj - 1)) / (2.0 * dy); + Ez(Ni - 1, 0) += FDTD_Const::C * dt * ((By(0, 0) - By(Ni - 2, 0)) / (2.0 * dx) - (Bx(Ni - 1, 1) - Bx(Ni - 1, Nj - 1)) / (2.0 * dy)); - Ez(*i, *j) += FDTD_Const::C * dt * ((By(i + 1, *j) - By(i - 1, *j)) / (2.0 * dx) - (Bx(*i, j + 1) - Bx(*i, j - 1)) / (2.0 * dy)); + // Basic parallelized cycles + #pragma omp parallel private (i,j) default(shared) + { + #pragma omp for schedule(static) + for (j = 1; j < Nj - 1; ++j) + { + for (i = 0; i < Ni; ++i) + { + Ex(i, j) += FDTD_Const::C * dt * (Bz(i, j + 1) - Bz(i, j - 1)) / (2.0 * dy); + } + } + #pragma omp for schedule(static) + for (j = 0; j < Nj; ++j) + { + Ey(0, j) -= FDTD_Const::C * dt * (Bz(1, j) - Bz(Ni - 1, j)) / (2.0 * dx); + for (i = 1; i < Ni - 1; ++i) + { + Ey(i, j) -= FDTD_Const::C * dt * (Bz(i + 1, j) - Bz(i - 1, j)) / (2.0 * dx); + } + Ey(Ni - 1, j) -= FDTD_Const::C * dt * (Bz(0, j) - Bz(Ni - 2, j)) / (2.0 * dx); + } + #pragma omp for schedule(static) + for (j = 1; j < Nj - 1; ++j) + { + Ez(0, j) += FDTD_Const::C * dt * ((By(1, j) - By(Ni - 1, j)) / (2.0 * dx) - (Bx(0, j + 1) - Bx(0, j - 1)) / (2.0 * dy)); + for (i = 1; i < Ni - 1; ++i) + { + Ez(i, j) += FDTD_Const::C * dt * ((By(i + 1, j) - By(i - 1, j)) / (2.0 * dx) - (Bx(i, j + 1) - Bx(i, j - 1)) / (2.0 * dy)); + } + Ez(Ni - 1, j) += FDTD_Const::C * dt * ((By(0, j) - By(Ni - 2, j)) / (2.0 * dx) - (Bx(Ni - 1, j + 1) - Bx(Ni - 1, j - 1)) / (2.0 * dy)); } } - for (Cell_number j(Nj); j < Nj; ++j) + // Right boundary conditions + Ex(0, Nj - 1) += FDTD_Const::C * dt * (Bz(0, 0) - Bz(0, Nj - 2)) / (2.0 * dy); + Ez(0, Nj - 1) += FDTD_Const::C * dt * ((By(1, Nj - 1) - By(Ni - 1, Nj - 1)) / (2.0 * dx) - (Bx(0, 0) - Bx(0, Nj - 2)) / (2.0 * dy)); + for (i = 1; i < Ni - 1; ++i) { - for (Cell_number i(Ni); i < Ni; ++i) - { - Bx(*i, *j) -= FDTD_Const::C * dt * (Ez(*i, j + 1) - Ez(*i, j - 1)) / (2.0 * dy); + Ex(i, Nj - 1) += FDTD_Const::C * dt * (Bz(i, 0) - Bz(i, Nj - 2)) / (2.0 * dy); + Ez(i, Nj - 1) += FDTD_Const::C * dt * ((By(i + 1, Nj - 1) - By(i - 1, Nj - 1)) / (2.0 * dx) - (Bx(i, 0) - Bx(i, Nj - 2)) / (2.0 * dy)); + } + Ex(Ni - 1, Nj - 1) += FDTD_Const::C * dt * (Bz(Ni - 1, 0) - Bz(Ni - 1, Nj - 2)) / (2.0 * dy); + Ez(Ni - 1, Nj - 1) += FDTD_Const::C * dt * ((By(0, Nj - 1) - By(Ni - 2, Nj - 1)) / (2.0 * dx) - (Bx(Ni - 1, 0) - Bx(Ni - 1, Nj - 2)) / (2.0 * dy)); - By(*i, *j) += FDTD_Const::C * dt * (Ez(i + 1, *j) - Ez(i - 1, *j)) / (2.0 * dx); - Bz(*i, *j) -= FDTD_Const::C * dt * ((Ey(i + 1, *j) - Ey(i - 1, *j)) / (2.0 * dx) - (Ex(*i, j + 1) - Ex(*i, j - 1)) / (2.0 * dy)); + // Left boundary conditions + Bx(0, 0) -= FDTD_Const::C * dt * (Ez(0, 1) - Ez(0, Nj - 1)) / (2.0 * dy); + Bz(0, 0) -= FDTD_Const::C * dt * ((Ey(1, 0) - Ey(Ni - 1, 0)) / (2.0 * dx) - (Ex(0, 1) - Ex(0, Nj - 1)) / (2.0 * dy)); + for (i = 1; i < Ni - 1; ++i) + { + Bx(i, 0) -= FDTD_Const::C * dt * (Ez(i, 1) - Ez(i, Nj - 1)) / (2.0 * dy); + Bz(i, 0) -= FDTD_Const::C * dt * ((Ey(i + 1, 0) - Ey(i - 1, 0)) / (2.0 * dx) - (Ex(i, 1) - Ex(i, Nj - 1)) / (2.0 * dy)); + } + Bx(Ni - 1, 0) -= FDTD_Const::C * dt * (Ez(Ni - 1, 1) - Ez(Ni - 1, Nj - 1)) / (2.0 * dy); + Bz(Ni - 1, 0) -= FDTD_Const::C * dt * ((Ey(0, 0) - Ey(Ni - 2, 0)) / (2.0 * dx) - (Ex(Ni - 1, 1) - Ex(Ni - 1, Nj - 1)) / (2.0 * dy)); + + // Basic parallelized cycles + #pragma omp parallel private (i,j) default(shared) + { + #pragma omp for schedule(static) + for (j = 1; j < Nj - 1; ++j) + { + for (i = 0; i < Ni; ++i) + { + Bx(i, j) -= FDTD_Const::C * dt * (Ez(i, j + 1) - Ez(i, j - 1)) / (2.0 * dy); + } + } + #pragma omp for schedule(static) + for (j = 0; j < Nj; ++j) + { + By(0, j) += FDTD_Const::C * dt * (Ez(1, j) - Ez(Ni - 1, j)) / (2.0 * dx); + for (i = 1; i < Ni - 1; ++i) + { + By(i, j) += FDTD_Const::C * dt * (Ez(i + 1, j) - Ez(i - 1, j)) / (2.0 * dx); + } + By(Ni - 1, j) += FDTD_Const::C * dt * (Ez(0, j) - Ez(Ni - 2, j)) / (2.0 * dx); + } + #pragma omp for schedule(static) + for (j = 1; j < Nj - 1; ++j) + { + Bz(0, j) -= FDTD_Const::C * dt * ((Ey(1, j) - Ey(Ni - 1, j)) / (2.0 * dx) - (Ex(0, j + 1) - Ex(0, j - 1)) / (2.0 * dy)); + for (i = 1; i < Ni - 1; ++i) + { + Bz(i, j) -= FDTD_Const::C * dt * ((Ey(i + 1, j) - Ey(i - 1, j)) / (2.0 * dx) - (Ex(i, j + 1) - Ex(i, j - 1)) / (2.0 * dy)); + } + Bz(Ni - 1, j) -= FDTD_Const::C * dt * ((Ey(0, j) - Ey(Ni - 2, j)) / (2.0 * dx) - (Ex(Ni - 1, j + 1) - Ex(Ni - 1, j - 1)) / (2.0 * dy)); } } + // Right boundary conditions + Bx(0, Nj - 1) -= FDTD_Const::C * dt * (Ez(0, 0) - Ez(0, Nj - 2)) / (2.0 * dy); + Bz(0, Nj - 1) -= FDTD_Const::C * dt * ((Ey(1, Nj - 1) - Ey(Ni - 1, Nj - 1)) / (2.0 * dx) - (Ex(0, 0) - Ex(0, Nj - 2)) / (2.0 * dy)); + for (i = 1; i < Ni - 1; ++i) + { + Bx(i, Nj - 1) -= FDTD_Const::C * dt * (Ez(i, 0) - Ez(i, Nj - 2)) / (2.0 * dy); + Bz(i, Nj - 1) -= FDTD_Const::C * dt * ((Ey(i + 1, Nj - 1) - Ey(i - 1, Nj - 1)) / (2.0 * dx) - (Ex(i, 0) - Ex(i, Nj - 2)) / (2.0 * dy)); + } + Bx(Ni - 1, Nj - 1) -= FDTD_Const::C * dt * (Ez(Ni - 1, 0) - Ez(Ni - 1, Nj - 2)) / (2.0 * dy); + Bz(Ni - 1, Nj - 1) -= FDTD_Const::C * dt * ((Ey(0, Nj - 1) - Ey(Ni - 2, Nj - 1)) / (2.0 * dx) - (Ex(Ni - 1, 0) - Ex(Ni - 1, Nj - 2)) / (2.0 * dy)); } } + void FDTD::shifted_update_field(const double& time) { - for (double t = 0; t < time; t += dt) + for (double t = 0; t <= time; t += dt) { - for (Cell_number j(Nj); j < Nj; ++j) + int j; + int i; + + // Basic parallelized cycles + #pragma omp parallel private (i,j) default(shared) { - for (Cell_number i(Ni); i < Ni; ++i) + #pragma omp for schedule(static) + for (j = 0; j < Nj - 1; ++j) { - Bx(*i, *j) -= FDTD_Const::C * dt/2 * (Ez(*i, j + 1) - Ez(*i, *j)) / dy; - - By(*i, *j) += FDTD_Const::C * dt/2 * (Ez(i + 1, *j) - Ez(*i, *j)) / dx; - - Bz(*i, *j) -= FDTD_Const::C * dt/2 * ((Ey(i + 1, *j) - Ey(*i, *j)) / dx - (Ex(*i, j + 1) - Ex(*i, *j)) / dy); + for (i = 0; i < Ni; ++i) + { + Bx(i, j) -= FDTD_Const::C * dt / 2.0 * (Ez(i, j + 1) - Ez(i, j)) / dy; + } + } + #pragma omp for schedule(static) + for (j = 0; j < Nj; ++j) + { + for (i = 0; i < Ni - 1; ++i) + { + By(i, j) += FDTD_Const::C * dt / 2.0 * (Ez(i + 1, j) - Ez(i, j)) / dx; + } + By(Ni - 1, j) += FDTD_Const::C * dt / 2.0 * (Ez(0, j) - Ez(Ni - 1, j)) / dx; + } + #pragma omp for schedule(static) + for (j = 0; j < Nj - 1; ++j) + { + for (i = 0; i < Ni - 1; ++i) + { + Bz(i, j) -= FDTD_Const::C * dt / 2.0 * ((Ey(i + 1, j) - Ey(i, j)) / dx - (Ex(i, j + 1) - Ex(i, j)) / dy); + } + Bz(Ni - 1, j) -= FDTD_Const::C * dt / 2.0 * ((Ey(0, j) - Ey(Ni - 1, j)) / dx - (Ex(Ni - 1, j + 1) - Ex(Ni - 1, j)) / dy); } } - for (Cell_number j(Nj); j < Nj; ++j) + // Right boundary conditions + for (i = 0; i < Ni - 1; ++i) { - for (Cell_number i(Ni); i < Ni; ++i) - { - Ex(*i, *j) += FDTD_Const::C * dt * (Bz(*i, *j) - Bz(*i, j - 1)) / dy; + Bx(i, Nj - 1) -= FDTD_Const::C * dt / 2.0 * (Ez(i, 0) - Ez(i, Nj - 1)) / dy; + Bz(i, Nj - 1) -= FDTD_Const::C * dt / 2.0 * ((Ey(i + 1, Nj - 1) - Ey(i, Nj - 1)) / dx - (Ex(i, 0) - Ex(i, Nj - 1)) / dy); + } + Bx(Ni - 1, Nj - 1) -= FDTD_Const::C * dt / 2.0 * (Ez(Ni - 1, 0) - Ez(Ni - 1, Nj - 1)) / dy; + Bz(Ni - 1, Nj - 1) -= FDTD_Const::C * dt / 2.0 * ((Ey(0, Nj - 1) - Ey(Ni - 1, Nj - 1)) / dx - (Ex(Ni - 1, 0) - Ex(Ni - 1, Nj - 1)) / dy); - Ey(*i, *j) -= FDTD_Const::C * dt * (Bz(*i, *j) - Bz(i - 1, *j)) / dx; - Ez(*i, *j) += FDTD_Const::C * dt * ((By(*i, *j) - By(i - 1, *j)) / dx - (Bx(*i, *j) - Bx(*i, j - 1)) / dy); - } + // Left boundary conditions + Ex(0, 0) += FDTD_Const::C * dt * (Bz(0, 0) - Bz(0, Nj - 1)) / dy; + Ez(0, 0) += FDTD_Const::C * dt * ((By(0, 0) - By(Ni - 1, 0)) / dx - (Bx(0, 0) - Bx(0, Nj - 1)) / dy); + for (i = 1; i < Ni; ++i) + { + Ex(i, 0) += FDTD_Const::C * dt * (Bz(i, 0) - Bz(i, Nj - 1)) / dy; + Ez(i, 0) += FDTD_Const::C * dt * ((By(i, 0) - By(i - 1, 0)) / dx - (Bx(i, 0) - Bx(i, Nj - 1)) / dy); } - for (Cell_number j(Nj); j < Nj; ++j) + // Basic parallelized cycles + #pragma omp parallel private (i,j) default(shared) { - for (Cell_number i(Ni); i < Ni; ++i) + #pragma omp for schedule(static) + for (j = 1; j < Nj; ++j) { - Bx(*i, *j) -= FDTD_Const::C * dt/2 * (Ez(*i, j + 1) - Ez(*i, *j)) / dy; - - By(*i, *j) += FDTD_Const::C * dt/2 * (Ez(i + 1, *j) - Ez(*i, *j)) / dx; + for (i = 0; i < Ni; ++i) + { + Ex(i, j) += FDTD_Const::C * dt * (Bz(i, j) - Bz(i, j - 1)) / dy; + } + } + #pragma omp for schedule(static) + for (j = 0; j < Nj; ++j) + { + Ey(0, j) -= FDTD_Const::C * dt * (Bz(0, j) - Bz(Ni - 1, j)) / dx; + for (i = 1; i < Ni; ++i) + { + Ey(i, j) -= FDTD_Const::C * dt * (Bz(i, j) - Bz(i - 1, j)) / dx; + } + } + #pragma omp for schedule(static) + for (j = 1; j < Nj; ++j) + { + Ez(0, j) += FDTD_Const::C * dt * ((By(0, j) - By(Ni - 1, j)) / dx - (Bx(0, j) - Bx(0, j - 1)) / dy); + for (i = 1; i < Ni; ++i) + { + Ez(i, j) += FDTD_Const::C * dt * ((By(i, j) - By(i - 1, j)) / dx - (Bx(i, j) - Bx(i, j - 1)) / dy); + } + } + } - Bz(*i, *j) -= FDTD_Const::C * dt/2 * ((Ey(i + 1, *j) - Ey(*i, *j)) / dx - (Ex(*i, j + 1) - Ex(*i, *j)) / dy); + // Basic parallelized cycles + #pragma omp parallel private (i,j) default(shared) + { + #pragma omp for schedule(static) + for (j = 0; j < Nj - 1; ++j) + { + for (i = 0; i < Ni; ++i) + { + Bx(i, j) -= FDTD_Const::C * dt / 2.0 * (Ez(i, j + 1) - Ez(i, j)) / dy; + } + } + #pragma omp for schedule(static) + for (j = 0; j < Nj; ++j) + { + for (i = 0; i < Ni - 1; ++i) + { + By(i, j) += FDTD_Const::C * dt / 2.0 * (Ez(i + 1, j) - Ez(i, j)) / dx; + } + By(Ni - 1, j) += FDTD_Const::C * dt / 2.0 * (Ez(0, j) - Ez(Ni - 1, j)) / dx; } + #pragma omp for schedule(static) + for (j = 0; j < Nj - 1; ++j) + { + for (i = 0; i < Ni - 1; ++i) + { + Bz(i, j) -= FDTD_Const::C * dt / 2.0 * ((Ey(i + 1, j) - Ey(i, j)) / dx - (Ex(i, j + 1) - Ex(i, j)) / dy); + } + Bz(Ni - 1, j) -= FDTD_Const::C * dt / 2.0 * ((Ey(0, j) - Ey(Ni - 1, j)) / dx - (Ex(Ni - 1, j + 1) - Ex(Ni - 1, j)) / dy); + } + } + // Right boundary conditions + for (i = 0; i < Ni - 1; ++i) + { + Bx(i, Nj - 1) -= FDTD_Const::C * dt / 2.0 * (Ez(i, 0) - Ez(i, Nj - 1)) / dy; + Bz(i, Nj - 1) -= FDTD_Const::C * dt / 2.0 * ((Ey(i + 1, Nj - 1) - Ey(i, Nj - 1)) / dx - (Ex(i, 0) - Ex(i, Nj - 1)) / dy); } + Bx(Ni - 1, Nj - 1) -= FDTD_Const::C * dt / 2.0 * (Ez(Ni - 1, 0) - Ez(Ni - 1, Nj - 1)) / dy; + Bz(Ni - 1, Nj - 1) -= FDTD_Const::C * dt / 2.0 * ((Ey(0, Nj - 1) - Ey(Ni - 1, Nj - 1)) / dx - (Ex(Ni - 1, 0) - Ex(Ni - 1, Nj - 1)) / dy); } } diff --git a/src/test_FDTD.cpp b/src/test_FDTD.cpp index 8d1c8b6..4dc5d52 100644 --- a/src/test_FDTD.cpp +++ b/src/test_FDTD.cpp @@ -1 +1,132 @@ - +#include "test_FDTD.h" + +void Test_FDTD::initial_filling(FDTD& _test, Component field_1, Component field_2, double size_d, double size_wave[2], + std::function& _init_function) +{ + if (axis == Axis::X) + { + double x = 0.0; + double x_b = 0.0; + for (int i = 0; i < _test.get_Ni(); x += size_d, ++i) + { + for (int j = 0; j < _test.get_Nj(); ++j) + { + _test.get_field(field_1)(i, j) = sign * _init_function(x, size_wave); + if (shifted) + { + x_b = x + size_d / 2.0; + } + else x_b = x; + _test.get_field(field_2)(i, j) = _init_function(x_b, size_wave); + } + } + } + else + { + double y = 0.0; + double y_b = 0.0; + for (int j = 0; j < _test.get_Nj(); y += size_d, ++j) + { + for (int i = 0; i < _test.get_Ni(); ++i) + { + _test.get_field(field_1)(i, j) = sign * _init_function(y, size_wave); + if (shifted) + { + y_b = y + size_d / 2.0; + } + else y_b = y; + _test.get_field(field_2)(i, j) = _init_function(y_b, size_wave); + } + } + } +} + +void Test_FDTD::start(FDTD& _test, double _t) +{ + if (shifted) + { + _test.shifted_update_field(_t); + } + else _test.update_field(_t); +} +double Test_FDTD::get_shift(Component _field, double step) +{ + if (static_cast(_field) > static_cast(Component::EZ)) + { + sign = 1.0; + return step / 2; + } + return 0.0; +} + +void Test_FDTD::get_max_abs_error(Field& this_field, Component field, double _size_d[2], double size_wave[2], + std::function& _true_function, double _t) +{ + double shift = 0.0; + if (axis == Axis::X) + { + shift = get_shift(field, _size_d[0]); + double extr_n = 0.0; + double x = (shifted) ? shift : 0.0; + int j = 0; + for (int i = 0; i < this_field.get_Ni(); ++i, x += _size_d[0]) + { + double this_n = fabs(sign * this_field(i, j) - _true_function(x, _t, size_wave)); + if (this_n > extr_n) + extr_n = this_n; + } + max_abs_error = extr_n; + } + else + { + shift = get_shift(field, _size_d[1]); + double extr_n = 0.0; + double y = (shifted) ? shift : 0.0; + int i = 0; + for (int j = 0; j < this_field.get_Nj(); ++j, y += _size_d[1]) + { + double this_n = fabs(sign * this_field(i, j) - _true_function(y, _t, size_wave)); + if (this_n > extr_n) + extr_n = this_n; + } + max_abs_error = extr_n; + } +} + +Test_FDTD::Test_FDTD(FDTD& test, Component field_1, Component field_2, Component field_3, + double size_x[2], double size_y[2], double size_d[2], double t, + std::function& init_function, + std::function& true_function, bool _shifted) : shifted(_shifted) +{ + if (field_1 == Component::EY && field_2 == Component::BZ || field_1 == Component::EZ && field_2 == Component::BX) + { + sign = 1.0; + } + else if (field_1 == Component::EX && field_2 == Component::BZ || field_1 == Component::EZ && field_2 == Component::BY) + { + sign = -1.0; + } + + if (field_1 == Component::EY && field_2 == Component::BZ || field_1 == Component::EZ && field_2 == Component::BY) + { + axis = Axis::X; + initial_filling(test, field_1, field_2, size_d[0], size_x, init_function); + + start(test, t); + + get_max_abs_error(test.get_field(field_3), field_3, size_d, size_x, true_function, t); + } + else if (field_1 == Component::EX && field_2 == Component::BZ || field_1 == Component::EZ && field_2 == Component::BX) + { + axis = Axis::Y; + initial_filling(test, field_1, field_2, size_d[1], size_y, init_function); + + start(test, t); + + get_max_abs_error(test.get_field(field_3), field_3, size_d, size_y, true_function, t); + } + else + { + throw std::exception("ERROR: invalid selected fields"); + } +}