From 153e3654756b838409159e247a547c7eae1746fd Mon Sep 17 00:00:00 2001
From: Hoffmann <hoffmann@mathematik.tu-darmstadt.de>
Date: Sat, 16 Nov 2024 12:47:32 +0100
Subject: [PATCH] last changes
---
2d/__pycache__/myfun.cpython-312.pyc | Bin 82170 -> 86678 bytes
2d/compute_errorestimates.py | 40 ++---
2d/main.py | 2 +-
2d/myfun.py | 254 ++++++++++++++++++---------
2d/plot_aposti_estimator.py | 10 +-
2d/plot_meshes.py | 5 +
2d/scaling_space_errorestimates.py | 45 +++--
2d/scaling_time_errorestimates.py | 50 ++++--
8 files changed, 271 insertions(+), 135 deletions(-)
diff --git a/2d/__pycache__/myfun.cpython-312.pyc b/2d/__pycache__/myfun.cpython-312.pyc
index 39b78f767f2996c22ecb00cf1bd39265a7a5aea7..0cda536201c5bb60894286bed115d358031edc8f 100644
GIT binary patch
delta 17792
zcmey>$U1E)E8l5eUM>a(1_l8e^Yqzb6Zs?<4{cNrX5?#Oh_XyEO|?p6N-=9;ncTq`
z!<FLJ!Vu+M$*AeR`32*9!O52&h;p4`Wnh@jkj_xcJXt|fjIV~lg&|g$fuWYQhB=FQ
z@<a<^QP$OrAoCd*7#Vs*SwKqJN?0fFTBfEJYsA1%%U;3`k!E11VXR@7WME)`X{hC>
zVFMX*S!8mb6gM;2KroxjtH*?ep#*9!149iPTni@?gF8bCV+%tH(`qJ&!L?jW3^iP<
zAsmKU?sSG)o^*y<-gJgqKCnwPC4Yk)Qp29YoUOsYP{f?hP|KgER>PCUJXz32n2~dG
z!6F&GdVZKCFq#qM8U7TO6xMX66lN&NlFpRQTq`iS!A?ZQg<%3?>;?vgT0xlCY6NNo
zL0$ymTA>=IES|{|Ekx=;Toftc6!sdS6lSovPzno(u3}(dsA6Scs1>OZfyh9dTPy+R
zvw)mmEKnl?(E*j`fb+rfMM*V$S<KMTfKZd?S+J?IK>{D$+H_`wcWbz_cqSXV3NwOZ
zR-L7W9ptDQVTf~Tgdx_V=->xMpoOqE5vKQSDdC5>hk>Dnw?+gWf1)*_Da_eS3=BmE
zHM|h>_(hOh1x^?i!rDk;NNIx?rZR=4ii?4vRutrYP|8$9casQ;%TfHQi|$vJ8cwj^
zML{HrKB$ign9c|B2a+mqB-aq<28ddMZqP;e5TXOc^a+eTn^_npF!uaoVL(mmV%XAp
z2{;=}chF(v608-g;j7`U5v$=zX9;G|<nWuEpr+3)pP5&ZT2!2xToONdx|)l@O2#5N
z1_p+ejJH^dQ*+XaLHSjoNR)wr;TA`Hd}dx|Nql^j29hQYbvemq1_p-DoD2*M4GcGA
z<vUz@{ffjH7#RFkPPSK9tM_4GU?{Q%5g<W`S?sx$$vKI|#YJum3=B;m?Vu8+h2b+3
zBc}pbb5RgT2dFqG;sCKg1zr(1hy}7plcPuq#1;S%pqy6Z3}W$t2vC9qd*T*vPG(+e
zePU63acW{wa)zdOkswHcJctkmd5H(?k@$lA%)FA~Tf7BDnYoEMFfLO;$t`AQ=i(v-
zka{POdfw#BqU6$Cs0y$XZ*i6sWyXWexy2Handg0rH8(TGwfGjhH$>Fp7IR{0#Vz*C
zyp+u3)Z!u~kU@;Dx7c0dGxJhXD{irX)Zb#vzQqG^2Si)3_+*3C;*&pUaMtsIO$WzO
zZfbExQ7QuiLp}ooLopvSD9Ar>Fo=jX_}mc@pTRcSyTSJXpU5W$MplImrW+h0*Eysw
za!Ai`y2>H{ftit2;R6Q)H*bUU9bTa+#uExBn{~&2Vqj)f=wQ6VD>*}PhT}|?%e?X+
zo!tFgom>-)diZYCqng6b%&PE(i$PRk2J>W(2H)>2OsonY89)RZgP_D@&IZp1+`J7g
zH+c9OTyO9RG`J(#^_f9|RlkGr1~<<H#xAoPLL!~c4|oLny*s_<FwN(f$+19qHt!W4
zwHrbr9nKF#L?<LpFrA?|`Q>B&`a6PR6AC9wd}ff~)bDV7ARstFaU#ov@EL&<qb>?4
zbU5A+lboJAF?WIEMKQGn<tqvoR$UY`>TtgyEH;C2hT=??8G$njXNJt=oM_+S@_=7>
zf@5d-48^XR86^uGXO^!}oK>@;WP{_%@*Rq+Y7Vd*2t3Gs!tfB!4FU1#d=vQ=FfXWg
zT*$VdaE0K)k`<0C5?8uhkT<*{V01%3^16WXMFHjO0$LXZv{v};P+aN1!~KBcPOk&y
zCkhW%UC<7?A`pB-KzKUWM6Mam3m9j*EihaVxX@$;;|j%<EElA8uL|gWVPg<gTA;c@
z<g$p~1jZY}A`=QHncfhQobEHxXMxBS5w#m43ia1T)GvytuVA?%qI*X|dWP$4t_6lO
zc_uhMfIAW4LLyvOdm%LGqI=SX<dlo<DHpU;uLz`Jb0x@yR|NDP2-L#@f>yR%6VUr)
zD8Z@!K~RBJ|APx711o2{ccb?drs*6LIc5keP3FDAB6mkvbh^bviy1zbg*6*|9`FcH
z@SKse!sG(C!37qB%@>7T8Ci?s7#J8PON*+o<}xraY})K1n#IUn6avb@AmZfa!(x$)
z+@Rc3R0t9*nyfCN##O?=z;F?iqn1uSA*Q;yM&bc4W5MKDwQBhy50ENP5a9(<fG7`(
zK;_!x|7z;qZ$PrXAOcjt`4#zrn5YFFxTq*9XJBA>4U!B75g^47ad0UID%PG)Uaf8^
z335540OLu^%u9(c&de>yNv$XfnfzBhO&C-z-r@iiIw_@zIYr?h|CVY@WbJ2QV7R%N
zTk|xRdL{z{gC<WAsE8{96?V5+3i6Ag<ybPvYET(elr#CaaTM!3P+_t;%47$VT2TV1
zXkso)%qa%7H53$}1wVUGW_m`6e_GK(1_p-No6XF*nHU#NcCxe+DT;@i2GR?0j6#&i
z<TaL}tjie~7?y56W|_;#STNbdI+P1k?G}NE1DhvXXEJJowJ3v3SO5xI6gxR`GxOq8
zGK))!HZd?TT-dB`>(0#hd~&{ng)zt#MW7ff0{OEDTm^s>X|flA4O;>-aVf|vY~TXc
zwJ2foEeC7X9SjT%*EUNz#xdD$1qBfEE!Lv^(!7)+5D#pYCSy?_$R&tKzQql$@}TDI
zW?*38X4<^fnTtti1IU(*AamHkVO*3j`J!tm+aU%91`(#sI&QJdj2k9*c`7l6PhRWk
z!*-5=fkA_5^MB7*OqLfwTH!=`C>T7o*JLjO*I(;F((6D3xF!SFZs5x9++-zRbrDcl
zfvcoj>?K8sY55hIMW-gm`Wms_Vqjo!WZFE{_YV_W2S|Iu<mmxMjG2>91^BZ;oI6=Q
zP|5%t?4Z<si#4$*HPN-`EJ*QH5OD=WfGq$glN%uJ&B>X8@~jUT7#NZ__XZX+Dc%K%
zf|9o;QxV9XTa2E?pcJgo0Iu7M?oE~tQDkcXwf*i-b_+2ucmmRZu+WpBh08)zSf7EK
z&YRDN$T2ZyOnw)p!U6K0CU?=r$x`7;0yjZMg0gxM*o~Jb`-JPVzF=Trn7X+xJdTMQ
z;$4uU%*oFqV_3g3Ffgp!>=6~kD7zLU3UdBRMsRvTYT^{#n|v@vj`8&5d(q*7KTuNw
zRPM~=)EIB<a+4p_u`>Rdtk^3$sh)E(b8HnO#JF>l=f>KYgA#xy7bK~G^6V{^qSUmc
zqH7@6fFs}(hy@ONCPq+)21KwhP8Nt$U}a}yU^uqfG;R;0G~D$?nC=ryQDVG1xxZI=
zvQMU<^ly-PAb(<4$#xGEK_KVKftCLwp!_~ex%%V@Ev)h@nf!_h@M~dwFu5(!hM<cd
zPD*ma5}<-9O6-q7b_y^~ceG@b-8?r*l#z*naq^aA8%5+0Vqho+wH*~2zyOv6L1FQD
zvPsI%SWv)2G8Q;H9S3DuKBRcyNH0oE!4em@m|T2|!0{mliVATMAps(gEuY+#s=+GH
z$iQ%8^Y+v|On5zeF-w*43sKQ?F-w*0E6DG}MGr{%H<)rz^sq_uf*gt?@-JqovVR9D
z=L3hV$Yl9C&dK(@{Dd6O{sW{55*=_&lWmf?u*4q7j-LqSvYVx{eb}@aL9PN76_6%9
z*nOgm*dq`WTE8}LFYMy1mjr2(1`$#qLWYrnVI^x3TD(pNB}Q=M@*_nqmq%)GacVra
z=w){CHN3^_;%h`Ok{Kr-tk&VuVPs%<&&0q`d~5QaRE^E@HEZzaMui4BCO+oL6MB`v
z`3qlmRA`W6<7Z}IC?YPt6&mEY1;ENdeion1H$!f6b%TK}C{Lz?gg^nS$qX*tikLw`
z#sUgtaQ)f<FQh^46K0;=R$r)T1X6>n37o`nX%c0g?9h;dLlek4Nb3Q`CLeK-dyPT%
zfSWz6j0_CW((o1+xOR=t%P(;SHz11GK+4bxMzD)Oh8Ka$L^%%EB6CIt22SS9sf~B>
zM}%XG5|avX5#iXP#HI>zDRB`2QmzJ5PB`N^wkWZygOrn!@j#k1K$?ilcp&AP5apA1
zweB_rN3Ad@v{@nP0a1=YLz97_n3aKn0bG8<N>5Ob=`c?|*Y*d~0B~uDC<wu+1gimh
z%#)QnzT+?e><oCMVs$~i0mvT?pu_-9GvHzXTmpd00#Gra$&6GCAZLCvkS<U)TO<o&
zp(Q3ckQkw&oni8!YU#<`UA%10j0_C&%$w!9_A=ouSnu?zvO9v@OiKKMG&zAZ5f{H8
z<<1c0lXtZ$ZI+$r!f3OS)vpMg@xg;6s0A#n8~_J6DB93-=bc_vMz_ru`X@1IdV<0Q
zxg3DUH5N_mo*-v<LF<tI&zzfOC%Q7?Fau;e7BkX)AZ9QyK*oH)Q4Ee<a6E&exJUz(
zczGdZIwD#%K~6=BRxOa2Hi$rv+sWsqXtVh;GBDUMZ~i|efmsREt6RwgX%>KM4$Rh6
z;|yiS?8!@K6i9;<7X^UK#jdJ!a>z_A`4&*g3u>emf%*VNpawI<d~j2K`sBGY^(??@
zz<nmL3&E{eNCOt68Qd!ax6YG5y%&zW{G!~%9B<d62a{Q5nXra4GB8AKww?8cnGMvH
zDcU+&f1W6$VSjXTz+x$Dc*mz`8Av~*_v2jzZpwqrJ^``?)Xu!cT$E9Ki#0hvH7yO?
zhM&A}u@o0`Mq)`(JR<``$L7QHl33ItK_Smk1hEXH6*Mvk@l!M?G{BAQqm#WC%j=<9
zT@(Y7#$s<#EJz}5@{Gk=tVxUv3|lrIU);o2e*|PWVw7(pWB`)2I3uwjRg(uY(g*Gu
zfvmm7lAo7aTm){+f%`@fXMlXnoRV4s?l>2L<CF!|)DtNJjZPMU#u$t0L5f&1OHzw&
zae&D9+{B8aCk&HMtWsx;o%~@{5o<Xk1H;43>8tfvg}||noEa`oUcNz|@!sU)8|0bz
zSth^Opa^j|BENyN9}mbgATNXzy#hIa6Dh>N-8itL!OjP}JOdOY;E-hlDd#Jy0ttY+
zOK2`HDg#Lvfs6<Dd=wZ&85oN4L4t*h3=H6L$9V6eqRDeNYO~dVDo2*h=Qaj1$$(-*
zlNBi*!7c#1C4RE?7I{&$*aN4{<jMJ4RM{FC85k5<HqY8(&8SodN)FhK0mWre1LI`h
zaEZzM+j!WT85tP#SvD(flVITn#aR*9#>UAmyRytd{jFO}86~%vl9O+-B<E!27S({_
zg#$dwS&*MogytM@r@mwI(_LzevnTWHwr5+;$iQIBvN>?KCbK9wUK>DR0_wR#LM(W4
z|Gs3_4U7y7(VM^T`^#tmPr;fTMe9J}wH`!(%|){dJk+~jvdTe4*3FCz427G$4*D`G
zZUXrnpVrAc4w-XpXJlZgV_{$@cAV^gP-(LC;kTOLGz2O+V5M*<$l_outuezRvW!zF
z`yWvYhL;|i9FRf}Tq1z<fXB9qK!wIF<}y&R0P;4-#9~n1QBX)wP*6~C0N3HS*wRx=
zKx%d~GBET_K5^ugJUFp|2P{}`F?k!_V)8Zu4R0HQ2&19}lP?~XWj(~mz%YCA=cB!h
z|0Xvc)8Kv!vSB{h&6BqrvtT{J$iT2{^2cMx*piqT7^bpLet1h^^4;T}jA@hAPn0tX
zOrC#2gQ@7v<kKhiGAd4HJY~trSi@w)z`)QmeKO;z>x|5kZ=RB86rTL^lrdxGWZl#D
zhDBRIiJk{*9sw7o;35;`iz1NYiq=hDc-q1dQD%ar!6hoVgae1wE|7yCfYKCGiYp`x
zkAZ~1z5x%41y7bb<IT1k6j=$Av(NM~UYgu{%}+V~lDO>w!2^Z|MNY6ElDiNbav?PC
zLVUu-(8P<uNf%OcE(PaI&OWQX`OMj1M#hxM3g_(@YbVE_uViM7{W<yed40xflesVK
z<&>!5>J+ZwuMwDReMxDv^2LKp(u_SHYPeFEL4#klTs1sU7S9C6o<EcSU+!jXnB0BE
zr5+R(MX8{~#skSQ;5>>g7{H|=xMDdC@)>vx>@0`{F7qUc&VdR9=A!)aTg=J%Ik#9+
z5=#=ngGNOcK&r5W-z}Ch@5CH$$n*qDSz?Y4csQ^K9+*WJL1wZS7bF%Hr+T3Ws|YBi
zK*AQ3x+NwXtQMWT?y4i(T~Hujp8WOdGY=wy_lCUkg5(vZ8-!L`EX=zoZ!&`gU1S!^
zgUSD|iEfU%zLKf_J;>G%AmSq^88H`~0r?Ib5TFqsP3~K)rSbVWDMjEkjBq1pP_zh~
zCcya$oG>ndwB7>wk*PGJ7&OzQppXCtuz}E`BOoOpVgtxib_|=FZq+c0wt__4Ktvaa
z=mrrrlZEev)c1jy{UBljC{hI=seBShU^0k66phs&KTH9MO$8CtK*V$qF#|-*1QD}9
z#B2~T2Sm&TnJ8K`55xwCJE&ZR`13f(6Kq8p;rXz_71Hzv4b~Nba$}J#$hyTKt>6L=
zoOQshWpJIn1SGK(M1ZWj#aWVBl9L)=ke6Px3>3Gb;EAce$Ly1%@AI)eW@KOxXWg8C
z-<=uKPy$t%ScaS~KGI^0nf&LGLM%ATf}2N&Kz4xAt|lL1u)63j$oVX8-mWf?6bn`l
z%Hp@!6H6<=;$R1Z^E#+Bgp2`#GJMg+$(4@{Szj<RFnDcV|M(HREGRcYD`H-xJO!?A
zt0u2|Ez1UODr8Q+_*z91WOLC=kfXrulq(?CRS<CvL?D7{<z%@xnrz_2mps|`jWpwx
z$(e7wxK@Hxfry;Vd*76@YHkJbLBRnl{lO;u0EHo1U)=t)Jmbm9)t_eQpw}3;*g(Tm
zg$707avqDe%ah|juM!7o(_{g6i6My$>?`@n3SXo|!8$Q|Xdo-EPtO11Sr4ypz(eVp
zLXfNta%RypP*kubrIsXmyMo91Z?S?nu8^Vu+zLU90MJ<PEhbOTTO5g{74aa&MSnoS
z0QNdK+Q7~cVqjoM0M}Ovlg+<MvobL;FznwP{#A{MF=}$%cPTq?h5)x3K;8Tz&<H4W
zya19?z|Q*zG8tTrgPJ(t)N+fd(4Z)8^0n`JtgK873{N(5{!n2uhyp1Dx6?rR8a(=>
z$x(!`1JuI<+ff7>%PTUU9Q#X336~W`44{E6kY6T)raIRC5@Tg&Vqo|&`NXd)x}b1|
zEr<ZM1VJ(OIe`(h;HHKpg)ut@G_+g8Qp=hrSi?A-k&#h>C55S$ZL;i7aqSlHiWasS
zrW#g|&mn6&O5i5eFu=rXm?y{o6lI5~VxC<1N6ZQ?$qZY9Q3CQR*eP&1_8R8ZAP0el
z)>#=CK($c{%jAQ9G&Ha&DFGELs0J`k7W}JR-^RjF!xqe-$?8{x7O<cOuqGE$z=HZ*
zkbngZpP{EuaIk|y9M;qU*;s@nRuMs7%)y||zyNCKfH0)3!{?J&kdv648lRhAl#^On
z#L2|Kz{$qIP^`!}`P?rf)?}s*#s`yc{C~j1%EX$?^nrtchp)kNdp#SY8!O|2>BqPk
zZ5idJv+yt`DuIWz>p`iW$v8!m@fKG>QGQlxa!F=>UJ;~ppWe>HsKnSheH9O*oMJmj
zC1+|%dTMfhPJU5Q2PmC&f=Xt_TTF=+MK#l3@-W)4KVW2FsFL-a{(*;4WV!<{Bd0%F
z%^(KJw4j>d7IS7^NfDxY05<_Z$)O08o{Qjh#4VP@iqzsFNF_0S1}~#oJ<e)EvItQ*
z$b*W0a5%%$_btYZqK6=tgWBOm{Gikg&QwIy7vQQwiiv^Y7Gr^65xmYQdIU0V+4ML*
zMl)7rCI$wn?KAlp4>LtA236B6kb)T8hC#GS*&yW(q@XWL%mH`KK)DQ3*nm1E#h@M%
zDBfW+lPoEjiMd5;ObiTu)3*sSPGwA(?jyvgASm80(Jt98)h^vG177IHKl$Mk5vESI
z$sf&y>l+vtx}-YTI+#18Iv6|H(^yj&TR6I;JJ>o{I;1+7z#>d799=RfYM5I%x+I|@
z>}f11EG-;e;vF0v3>~~3+#OsUJRO{A>?y1*99<F};vEtlk{zra(j78sY$<Fl9JK-|
z?6rb5e4SE~3}FnQ^>*S647EZ?d0VJP0G4rMm>C#qg=<A>MQa3WgsXT#YtecnSr`zD
z?mF2(1|coF6RY8$Jkdr}Nw7vNg`<iCvgi!9hKiA)r-_B3hCi4=lhbedI#osoUeHJ}
zD3z~dEMjGv{!*CHit*m`k0OjJjOo)oMHo}qG(m|#Yx-spMq8Ln#&j`J#uO!x%u2>0
zP>ENh4H^yy#|iS>3&fzT=?6s_jVAYh=A5kcPuQ2SlckfXhM|*r4dZHPi5ARI1TsvM
z5u9N)q2rU_^r-`~9$aOBQYyyiAlL@!98eYqwdsB|Ff57N;J6`qhsFWJ6OLC@{U!?}
ziA<j_#%SgZt{54MKpj#|4oEWy)MdQInw%V;YX}}Ff%G>)V*^E?Sz~Y?64Ijr`wBje
z4YrjrZ@Q~EqdiX-Xs(fMIs+3!C)?!4WC=FVI<-!g>5ANpN|S4o1t-hIa!;0P<CrcX
z!Kf|^?gFxdt-mFJT3W>$OwW>FG!w#{!vfpKRtWN?$Mm}rjF$Cf7s@Lxr&V5*sk~5A
zdr`W!0Xl33uGzq4FDoQ#feQ_AzX6mtAoaRzktHZpbwPv?hyVw&A&6xUBESI;3U0_)
z33${594JdbB?c&Dav;GK4eD937Nr))8{A?A)$#_QA^eim^rFO^c!OJvA>b}5q)~Pm
zq?ReeyT}G)g)NA%0}<f3bpWv(L1{=9+#CQmwxNW=<cX;w^}Ok+CGjO0sU?Z=MXAL_
z?o12}@u2Wn#lXPu|Nrj>h7Vi}y!;KWNK1@9F^ICNcQ8H|_U+)fFA~tf`CKBrgZG)V
z-vthd=OU#SMM^t3zq2t2)-hfXsO?~VAS%|ubwf~eLgGa6PYg_)${mh(R5Vv8u5exK
zchQMQu_Iza$_&#RJYq96F7haKdVOGJ;#B^^!N4apL2-iPM3wHS4yFejy!||#JQJc8
zFfDLg$a+y!^(u$j2W?SS^$&X6Uq~@>Gbw-)dXXu}zc})Yuq>l0TNfJxLy_5ZJ6XmM
zN4O5yvIcPcfa45YD1uwX;K&5U8I~!~=^!H+r%nGU%Q#!e8e}y%6|#U+VUgYRWpa$#
zjQ1d_W=;Pg$LLrOHWQ;c2(|=lG1v(p>u#~8WTvGVK$?W0&Y}UhkiNx|SXyClixV`G
z14<1BkVF7(nS-O<2jmQJ+#~$p3rYas_Tcm#@{H=NK}-w`_1mAwGtOd$1e+E}6OLfx
zS7TIV+&uk=GGhvoVgrzonB@{!@s{a%DvbS#-~a)&%~wLQWf3?JfSm?Omeb`_8C4kX
zO$Te*Ha%07u?@)vJkDV|2#Uqf=?<2Ro=7_IxQXo$NJltG2WaWQ^fzjZp$G>SMS^^T
zrFa0jjO_?WYs~asbw;UraE1qWr^JfDk%pShZ?P3+<bw*mC{Tb#gDP#-j1rIl*zfSj
z1veIMv6f`!rW!yfL&V4vBt}8m8ZvE)9Aau3j3%t9ObiSww#REQ7BhlJ*)<t&@q==4
zZem4zp#iu@u#)8#Yf(nNTUgP*>0dM%tyy!J7#PlPH_>7gVpIYTYGi=C3>wM+RrcUU
zUQs&8(IDd1^bBoAZ$EhJmVto*6gZ%E?dNsO3=Hio?X2x=?d<Iw?VRmG?ZWLM(-~^G
zQW&#A%f#;@Eiz{AWQA3^wUY(eMW;9DFp5r>)L{%}?3|vf!>Df5310^fD&66V!K={`
zwQ&h(<rcV|+{poIV8li;FmwrZaCESA2z9V^aHg@QfSMD+9UL939YP)KU=dJrLZm~u
zLnMs})Rd^@MKsqtIXZ<v4Za#)n2q2T0b?y+Eq@8~bSnc!S$;-X<zCBJ!;{57Iq;&W
z6?YW_Xgxk?y*W&lk)bD=g(2@*4MVL!4P&h!J$xUJ;d@Y{2-Wv&EgW4i-?O)HP{H>m
zNWKTH$_Ir4C_IotA`!G?x`v^KF`XH_{<(;0GGn&nbbnDsZ69!Hi>IIg^_Ubv$rw@@
z<>#iRC+gYc<R>TQ6x-?9<Rs=Mr6k&Y0cA!=&POs4d(J1$#IGO|VI_me<kAgN)3qfS
zx$BD%<rj3GA5v<8@(m~_D8O4S#M$$mje+4aDAE0BV3?V?f^l}v%7_h#tD`T<S?n;p
zC~E_4X(WSI#xXE3On<M;Xo=)S?8PB*R<pA+FceWW&^SOQ;tjN-5^@8I3uJaFs4N1_
zNo#TxfyypO?+rYq0}4aPgcZ0n)&#Ykz(um?t$LQ+#0tY(AkwG^UIrI|YgTYG3{hsI
zmBs9td1;`@KuD1b9(I9`;DXIF29qY>5s)H8N&qd}2bG|ZHWfFB3BnEz4%21z7{&SP
zm>3wIgW6M|e(V1q(;f8~CF+G#mIN)2SQxP)b#eTbpzRSGBX*=+Hg&zI>Ut#Uc*enu
z6BU=;BCbS6Uyn?@7@2q>Df@C{&PBJJ%R;#w><{<^I+*VWicKk=VK}*Lvg{`YPEMT;
z$NS=e9Uk{3LOML32-$QvKT%L!khp?zVaf{Q#rzj}WIIA8s7xpYE$%g+XkI_VWU_6C
z%L5+%ey>ii87woxXYwpaxhSr6l}G!dFej(Z2Pp;tkq(v{{K6e9clZ?-sI8E^%x~Pm
za)&w&uSdAtj5L?ya`+094J<1yFN&L7<uUyr%*ko?L5zV{U;@j8kPiC`EK=Z^?_^LF
z2ChR47#JAf14rjr7#P|a+8NuKz_prKJG^5H?$)xWGt_c`y0+8P^ch8^nPROO7`h}n
zSUQ+GBsv&7SfO28sZQ3(j_D#Cof4qbS<5+j!v|4DuE`Idh_bPP#U=}i2@6bMjJ?6Y
zP|KadQp=OVQp-E}0J~WI1jg7y3=Fk=H4udiHQXR0YPf3nKr{%~@*`?OG*$et`m0Mi
zg{_7kveT%BKZOaT4^%jV3h7#b8Uct*3V5qiu|$mkL`4cqwz>ubLoq@D6Q~X^=71{z
z^NTDYyIf?L7&;kh1YiYzr!?Glafkv1oC+lKQd8I<io76-QrJ7CL9PS25@sTr3>TvR
z4N`@$52UXPYG)prCZwHmXfnvgaWOC;jO&st;RUA*Q0>m&$pVfs2AJn+1#1LTII>+B
z7>epjI43tgm#t@psHkDC;Y4g)12<rxrFRM^Xk%8fa1B%?l+TT-10q_)oWjy2mgfdj
zi>d}>Dx&$51uEzmCf~ovF9<am7IY<?lP_MBWd*xvvf?F?`YLV)hAy!@CXnGBVjbcg
zQXSGA3~8(>TrC{6LP#Tz{53)_uYy}akZ}>78sQpIP&k7}M|e6UJEYQ>Qn-oI%UdH{
zlY`LAJAtw1D0oZNZ5D<a_AF4D2@7s0UC#sa2O~ob2UZpE4Qq@H6PSAJKql32q%#FG
zX!3xETpS7a1Zs>JRrzE=)r&E#dYRn#NJOv*)VwUJ1XU#9(bg@XR*O^)X#9b7IzuND
z6GJC+4SfG)73=g3;*84mtThZNjGZjtFhWXNoh-1>>|}(=r7(d;jDi_7nf)NvfJPF)
z>s>V2z@5cgd~Wfs@%d$`Me*PT<G0un%Zk$DlZ!xo_gm0qDj~tQn9Je~AVYKw+@PTo
zNW=Z!bdb-SZP^(Zib0c4KN=Von61>@;J8}nqO#eRqU~iH%XTIoP(0vqBI2S;=tZZ{
z%XVQG!s9QSCtOI(zi6I+Q91uYamhu+640p5q<T&+P@iKe({zST#>o>KM2$EY7=jrP
z4I=1d%PqF7(%b^5Bhi)~fn8MZ2yzi9x<M{dTPeB0aJ9@uCDR>@7Zoj`?Hx^Sr12C)
z?-)FC37TQK#g<tdoS)}f1n#wfr((eqTcC`fDO1!6ib!xL3|zm1x?GS!PEa4Zz6d;!
z!j+t#R~BEAT9jLSiwo>zFuMqB8+b63G3OSji){!djldJMMVCP1Dxi&H!I0%W;5Gq#
zDh9H!Ke(s~WPdY=00mi5JBS5}v!WIz1_t3`8OW%M0(2rqK><AIGX1<fBkT0_CXC!B
z6POqnL_ndmhlzpV|E~sy4?GM)A`M<Q1ce*Ckj7j-GiZUfkb#F=rmHD1szp5kw*gER
zFwQieTyt4a<uijKr+SCu4PlWE7cg;0RBDFd<cQA<a=hvj7;lJ3%n+RDKErXM*NpN7
zg)^%zifBz>d?Ky2!ej&EO7jhdE3GypUY53)!0}-6g;^r?H-yD!Fitd|kvPd_fyoNS
zh2|>~7ujqu*}=Hcd`IFYn*$~%7!R7CNIYb7M^=7**vznn&MR0Jd#=d1ENec&4O48Q
z+XGpJB}&WH7OJfjUJ<fEV~52>Bj<~{&X+Y^4hEkPIgxQ8JpQ6z{6)|B%dQC*HR}^D
zBxYV#%DO0<jbhx2hz%JRb?q){*j-k!KgfK*<%GnAps0(UQ5RjJE;~eDREfS2lX6)$
zb%NUy5vdu5kUeRW{U$Iz5D}fhGRfnPxa9P*iDfg=76i`BT*0=Ya);o`nu{7X7v*d&
zi`yPxza;K{!t6p&aQ#Wk3*m8>yyGq;Wn9QC0_Bd((u+xD9~c;_7#*2DF)$=EIx^kR
z)SAinK+9l*;7ZL6ZaXXw1n#uCXyS2E%j1mEd9#yd7lI-#L`Gc<ioWa}b0s$ZdTjE=
z*yIZ-Sr@WPE|ivC%r3u}QgJ!9@}hU;Ma`<2d=FIBCis14V-PQ8te?PqM_p${;ELkK
zH5Ua`CrHdtnNfN}Kx78T0@0a#7X{QNusjj7B&28(%O_hoUiA;`ik#{n+!z@+g|Bl+
zUgVIRp>l;o{sxB-jCYMg{zE9E7OTRC2u228@d;HISh<xhuqaKQt1rzK!33IGo1A|_
zl$8m*qy6A{F|$}328LRe8pax?8t@)kP<x}6HHEc?HH8_(u3=4K1ks>+8&nF`vVnKL
z3)%=zzIU8+^38P7>Af<HoP3C4A5^4IK6qY?BlalB=n|I6AJb(fH@*~8fttt+u?bWm
z!gQsufvZ^^gyPA5?rd^3xRo;YL_sWx1I+{?^ua8F?zJaui3x1)Iu@mjJyA6*$xMt4
zARNrV!q9V^g`tvNlVftjBhh+DZ3AA20-6icWP=odpaKo4mxGwK;Q$v!kWvj)S{22B
zDk|_Y0`UAiXeA@K$byWhgSsE^<rE#D@(4tDgUT2cMh1q4dWH{d41z)v6gx^MB+g)*
z5IVzgM&b<93w(+V?himaO7HLrO)%_^0Bu@T?O?pY!`qQQA+a;7)3<~12EWh*#|eoO
zU1kW*Fq|ne(X+F@gXIQfY<5b}^oWTOGg2nUUFA^uz|PF7`h|;uPZ(kp3#)7g;~j3^
ze#_2!%YM5~JIF@Z1)>WI7fQ_3zs#+*BK#7!$sVTtoI5%9^X=q2P<TS{V9ANp3n8H=
zGY-~Ywg|rvm3Y}A5n8ByVqoHxoxlj{&{$5ioNhPKZh_(g$Au~@1XmcY6j`W!Sy*p_
z_$6VBJx=>QcY5yk+v#`0@I>HAlM6w?suw~+F9wC3w7zTuHq$2YhP2ESkB_WOys{ry
zSy*L1@G$W3O^-2SY*u&;o@WHFSI}fG0<AkJ0xf(f0?nmsauh{P|8K@9!400}1I^<X
z{hO|2&Zx;cnTdfRaeII{<9tS>MhW&F{B(XxMs2n#P@_Z-G(1VnfF?*sHAsgsk`8=B
z-(WjxK|0LPbzmPl0lSK^ZaSwm<DDSTK!_$YByA#+C8Xg}1R8X{#RpnJ5f7OgjxR7o
z^pA@`^)s}=g4|vKIf$`gda?~;XhZ@N1B0d@Bz=QwZLl6}qpL+upw#RPB3wWOyn_yQ
zMA0me7`PUObkC;?*fN^3&Sqj@D4gzP%UI0lGku>eqbuW-=`41P=JlX0-bLVjprBo=
zkll{`pqvHTt_j|USp?c%3EoW!-XRHD-3i*D2w8p%-W>xTf&e#2K)dORK>Op0K<lC)
zGxDH`_ael!Ja{tRmyv;?D4emLfdRY_xCp$k6SOF{2(+NI2(%m&JeO7kT7y{xT4GrQ
zS|?cqTKiapSfvMACkR<B1zH{iNtWRCOW;%rN}!P0Y<5tnfTrt<K$FDaiCpj`k5*9-
zNCq^P3>jbn4g3~?hAWCdZIU8Tp-}`%$VIcK$JsN=axG+HU;y=Qic6=r+A|ulJ1{aD
zesW`)zTKYjmlq!+-zOzqMuATmjC_p3pOj=71wT13sxb<F3Sv}Yl>Fq!$ib-aNsNn8
j;gggWqx7dJMiWNmPkD@jjC`MTco`)=D~K`*g1hMey$wG1
delta 14120
zcmbQXl=W95E8l5eUM>a(28P(J8tH<(6Zs?<_iR)TX5>q?Ok+whZDEPBN-=9;nB2h_
z!&1qp>9+YD<9or)q4M^OthbmGODiVNR8V2{VPIg$-+V|Ri&69zZ%$@jYGP4*acW{w
zaz;@Q0|P_bW^<)TMs7a_28N;l5D_?ewz3*m2m=E{9|HqJ@zTjQN~)VbC_mt3^q+ju
zu$p@%qn{>YksJd9gXUyUDcQ+6Mrx8PK`J>J7#JEDZpg}axb*rJi8C-T_-RgFA*IF`
zJ^8$pjAl3k1H%fCI(7yIh8Bj;OpKfgU{yu&AZwXG1jvZ=$#;z`B|)wRTXKshEi*4A
zzBn_tASbn=C}^_2ahecF$1RSMqRjY|(!`uwOa&#AR~S!ZEn{F{n6}x@<TMwf`(!SA
zDXv?L*+ohqgS{sk*hjIpF)%RbZ0@(;!Nh!v(RH$qlMHJY0|P_E=3FOkCdRJGl`eLo
zx7aiDQZkcMiy<02K)SOgpK#G(oyfqzkhJ-yOD-d$|KxPHP_6`!9Ehmhe9SGAksGX5
z86?;_*~&wjbq)gqL+|D&4|it9^vR37EsQ~y6@fyw$Q`6g07QTkX|flAwe^AoK%st%
zEi*4AwZgTiV=}LgHS0nK28PL-eSPAXZ0CW(o%t4PQGRJ&N)d<$HcOMSs1y{oj0_A6
z3Pqv}3=FrpGxIV_;-Tg&VPIf5!MOQ_9~YC*ERZd;LFTX*$%EJ(lUV{o*;X?!FkEEZ
z93L3V%s6ZE-Vi0mtjP~TeAu=#Ffcr2+-wy3ib-`RNI9HH4+Vp#_L}TPd?1;bAY#U3
z#|UA@?UO?z)I~Z$0tz74u$L4irsY>;7Hyt9A;O65Fara_KgP``BK|P3<%1M;Og<HD
z#Hcp;Z?r!f#DSCh<E0G1LAncM0Bd4VYNBh=Hjv_dAYw0w09()iVjTn#hbGU7k!L*x
z3WUuEVhWi&j)6q2L4IW_0@-tm(X$wo(i9rNfmC#yfq|h)rm!?Ir6{qav?x{2CMQ2R
zF{jv0&n71^Hz_62&X$RRq3Fcq%6J2VvmnC|j`1YIF%p7Y2*ZocGcYiiZ)Qr6V`6lh
zteB+20SX39?xJ0jeUp>~4uOpK0uf-p?3rAjq|17dfq@}t^Twn&CT>WufaKjLOQpuJ
z-eX{3$lhF?8pS9(9V80!=1N9zN-9RAgQDY;H)hB&c1#vb4;Orhnns{<os(yzdt;Yl
z>jHWI>Ey%`(a9gWI47HBR53yf?4Eof!zKolI5fE+$tVjHO)N#JX-P%<LGA%Z)n*V2
z980f3tT!OyEyxw1a!mmO92^|Lv2lw%sWd4mC)K;?0|Ns?^X81qJ&e+DuN6H6X+!l`
zwi4SckOLo0PAt)vegZNa6bd+$&ITz3`9)v)IRT|}U`pjCC**RfFs3sFGpuCtEAq#$
zk8$qgojEpy+{`#{vT?2(wgAsoVw(@L`6npAr_Wbo<lKBASCo<ICBx+Bc{Yk@K?o{H
z6&k<*9Dbna%vNGtFgZQ{XDldYAi0ts<mL^aJivz(Q5@+-i78kj>lTxXZxJ}c{)6-}
zFfuR{F@gwW%O~$G&|u|YWMG)O`BlLlCcGX#Sggvn5#%Ytkpxn@38a*8B!QG}hA9O{
zk|H>gioStdjUy3&G;INC`aU_a1f+?dpv$*{l!KxJq+Ak9oPo^T22&&oi8J3aA2#h*
zAcug;nj#Jm3+$wSpd^eX+71@0GH&1es-}yxo(UAzEFgjzM6fb4Fsx)PLW|8hkcr?Z
z<wuHAE|1jW;?#I-k;?4iYj}&<#n*^X)G{(m{?ewyCCbRau!f0&q4?J1_XQf8gWA_H
zD&Wb9pi)$Fa^rl-$pYytOhw?BR=|@T;p*AWfPxMj-3oZ}BwRiBS+IJL_r)hWE|#0T
zuG_!}lrddFLZC3#WCoXjMTnAy8B`rKAQ$*xr!rk&n!K~CP)`~Z%*Z;y2_0NU;nI1D
zX>xIQ4i23li;HpStiJ+sz6{7lQ2VIp9Vk~qiqBhI;Ho-4FTcdK2vn37y$5*-tt17z
z4P-+R*nxNoQP&~`Mh1pspg2rtU|{&w!0?faVRBg}-{#FdhwvAx6Z(|c9)cVSjz_%3
zDoE)gkWz3I;w@G|N*}|Nf^ril5jhs5=?O>^aXA*G{3%E|z8pKDPl^2*SP{eIw-a`o
zf)l_WP)x8wk`AKGLxeOdBoV=r(c}qzN=z@ACNoa@gJ}S`<U|yj;HbiCz-y+-VUxe(
zFaYcfc!U>YGob!0$RC>U^aL(Ez~u(G1OXKun#@Rr2Xc`i57GrrXlx)BT0&z7i4m#>
z7$z%E)nL<MWMH_%v^j9<UM6XLr7!<XRknX1_YxPEAm#r-%JIb||4dbO24)6^qMx9+
zRN5Rc$Ayt~C97YN>*Qs#<R!pC2#N@_3^|#9rYa-zX6D(Gm^6(*K1I&|@YutmiH!|p
znek-(GSSWd=K3+>&<(N|i|%wzux^lr;Oq~ME^thN;|vsCMS_eB44S-<QXLURLLfb8
zQ6vl!69EzE@idunfi{~NBLl-vrp-nR5}1`D85kH=GC`UH;JO2|EwyE_GNZ=i8;c91
zL5hnkLFQstWi+{ciI#jesB{9gUyDE;haym$8)E(pkVM_&3rqAY!D=Re_+S@;8@Q0R
zEl3R9F#|R5Z?WWNrnuhX$jdLvP0aCjEjl^bY^e#WBO?QY_~!hjZ<yIY9iXCllM`2p
zLOKWQCfBT#vVeD&iuyrKPe`B1y9nGJ02{s$WDTeVdW*RzqxcqUa(-%BT2cPwb1Rit
zJsBAoTsMDRnZ%;z3<`3NB8Vj*EuhXV#51m-U;sC~*G;ZnC$EQUWsw_58jF2J?jQ+|
z$*0z7vHCDFFcfb7xvq&#U@gdSaBBcma2Hif?%paH3bNuBOMYHzaS^zsCJ9mmu@dAf
z=9JVDa7VZZ9C0k5a!RBq6Qn2uM5KZgv1FE{7Tw|ik@2~S6-8$!|J<s<=ssC#TM=s_
zBLl;{%`>;@u?m4>7CEW!ntXG&Jmc}nKX=PBonf9VwMUV0?PQxh67>}f3=APfmqC_s
zB87w$$VRXO!A=FcG6>{VaEP&il=Bs(fCNC@Bs5nRC4eN1Kt_OjIiP-FQ6xw#nvsFw
z7He*1O1yVb3`k(|<vr?bX`uRxdDGrNCK*sTYqBE65ZDc1*LY6O-7hbS7BArB)-ZYD
zepR+?Mh1p^%$v{bw`Nqz07W!*V?eP~l*Kr?AxT1v3)D1@Pfjb%OD@V~WMFv1%)n5b
zG5NxN)6G2x?=y3QVyp=4g6zp<N3+ahLC$B&D7nRyoP3KVIVUr>C<zoR9GQ6~sYM0(
zIYnsB19$fGCyO6bV{Dylcg&t`A|nICFXqk7$26Hm!LgbK3MEjF9ukH&lMkLuW}U^z
zz#y?%^VDBP19)21<S3c}3bC0W0&FgtUEomeoIc-;QFL<l86MWTj0_BVn`_SaGAhmi
z1qVL;lV6=R=UT|fz+lb7z)<Wsx$%tBWZ!ddjlt;%R6@WSTlOH^Z5SCC&<Z_JyG<fJ
zwWQ2@vS7b3W6<Qr^J)h0@<NjXQp|zMfFh9YB5=9@6%)6Z%RuD-C=~1`Uq3I$x|ETD
z!F%$r^RFDii4Ht`!Fr3y+wc~Xx6v(TXJ<nY0cqYss!CAb65K7i#gU#`l9E|m;##zZ
zk%1v%vg*ZN#^;kaU)12f2C|?N?D)ygE?Tf|Vq{=Qo~(T77~3pn28OAulOG;bn9P6K
zld)-X<mGZkoyiw3Yp_gUEMlMh`|@5!o5{bgI9Q}GR<SZLq%c)+GceRJ)G&b(Mh!y^
zTMY+@1;SOJU|?ln0MRv!HB2@P3=BP+C;z^3T~Qvy2dU}gOJN3?Q^m!=01^YKnZVei
zKl$@jW5(9W>euW|pZC7qy5Re<4R)G5ShFO!Fa{U7puhtAwP?oV>DMeA5oIn|8eAfS
z3rKLtF9G@gBq+f$rMN<3WCKVD95}+D5VDypaNV1&0+a}*O-{VtH(BPPzjE3oahn5-
z2NVyo94I@)dm%9BLU7E5*tm<q@fQOVE+l7O3e23Gctd;g-Wyyx35gfO5-%hrUkppW
zkeYTOJ>z0(=7p^63pu41%E~Y1R9wueyqH>fF|6vMO4a6`8^Mf>Qzk#SX~#Hqviz+|
zW~Nxd$+K_iGrpUA`qo}XgULH@%Wxq=ee%USDtrhoD3~{IynUESnyH7QhAV{`6za8H
zH9Sxj&jh9(k;(n{x*2ComcH*&zZK+Fe^5@~f#hj$se&ym!9!f&+HD(14|w!$2Z#kO
zlO>9Ff{Fp=qWtn(%*pvVw^&jVOA?E~*>5*U6_!Z2#Zu;(nBxr|X}rZ!mYCyHgc>77
zdq8Hg7Z)TJ6{mWk$BhUmZXglF$iOh!@0j>xy9bVJvq16mb#m*2XD(@k<LHLG@`B_Q
zrW=G-S}e@FC~vZ%c!%l6$^(KMs}C4%ti34jGI_%z32r3*<_!-IGS%Mz*?kj4+ya&Q
z%thNlegp>xXk1E@`xa|ye11+!5jb}tTnZj=1Sc<Wh6ksyy&$beK|W<F%_!yuwX_uy
zz(7F(GMHJk4x|J`<bphB$FTXslNx5xe2{1Xh$sOOr63|{a^s7T`U()U5=2yiqE-Ns
zVQN4EwIBjfgeQXhPzMsL2N4Y*q7g(ifrw@h(E=h`K|~vfXa|`nTGRnzgF_!w{zLq^
z5#$NBqKxo-SS0}IQGkXXi$Gbr$QEQ>FGwpmWq@-bxDgI+Mf8Cr`auN9x?7wjnI$=?
z@dbJ5MH4`UtLWr_V-l0)WY{L3e96alj*)@k3d`oZFWs3TZ8lKNj%Db}{;d|H<m9Zk
z3JKuc4{oV}3c?~#3fAN+0{a0RFDF3GXL0j(b%CT?uzFBwaf>~%v;r&+b}%@Rf{N)P
zaH0m49pIEYIgVFyGV?o8*2|0x4D6d_-#ucN1!XO0OM(|EgMk|i36o{M$+CglHSUw`
zzo}?~98h!x<bQV%u^&Vn01*d41R}&HPwx4q$p%h?4U>0&lV;pI`T936uE`)(AVPDq
z_V-d&&G{fcC{$KL2L-?;JO%|NY6pGt+Q0IQZIfUAnW2MTx7}g`wFU|eia>Q7TJ;9j
zHevGFzpKPS+B8|feQrov1N(|&a^F8GQLs*o?i<L;DU<L1^Q?zgh>#c&g5-0MGm9>Q
zVuURzwItEowFn#&tRRjnq?~|Pmq_sd8l=9(<mq{fBeAq19;CSF1t=K6UI)h>*f~NB
z3=9bgpgK}v^3wm(tnU~Z7-}{j{jbKv=ra8cBcqfZIAegSCWxm%L$a`uUT}&5yYCgq
zXmGOw)I<ZPnp;eT21Q=eU6>g4SU)l{Ff80&&cvv~WZ(i)3LXyw<#6zzm?lRN!WK|B
z5<Ku*1R7W@QkZ_4g;7cxmn}ta7^kr^8nS+7WMJ4j-Hes-3M1R}08T~`6;N}iouQqv
zovEF(oohNn4QmQxb_@eU5pxY|7gwHOCu0`l<cYz;lizRPtncDVVX9>-0i_1eJR3s_
z^J>rtGXn#I3q$NA28LSp4z3QSG^P}m7LG2?4u%fS4z3QyG?o<B7LHnu8nzndPA*Uw
z)v$xe5>O<A&FN%-iPx~KhMHai8u&w%4P;pj3K3}gpo9}F%fP?@*AAC20R<XZ8ZO6<
zFdE!)U|?WKVVf*}Sk9zG5Udi}ASA^lpy&fDfE&fKnrU)go|qQc99{;7TDG1kEDTHx
z^&B-UumGuH1J$I#44UkIMJ|k!<^D>7X9NoW=x4(7S<y4hyan={f=^;WPGWLud~SYG
zPHH9QEXyTU28Lo#68+J@FhhAp;sVB*$ul(;7%wPXA-J%3h2g@|g;p0N^(Ig37Mc8S
zF6Z<+5{$y0Or4CKEHw<B%xf4|Lqk58L6gZ(lkpZScrZ<q@s?<MYDs)?YI1&FN_<Lc
zQD#|UNoHB9CUm4n9#k58F-(`2WR$5tPLxY8gIo%lVEfU)a7EE<L-7v79l-|#PY7O7
z^lo5)OlX37cbZ&CX&clRhoo)Lz$zp|ftwTHQUaXvIY8BqCO0JSgJP%%Thf06%0!@C
ze2Y6Jvn(D|G{u9m8^|A^-~-p53DAtK-~gRh;X@7nqF<n%704eTcl~K#_`t!y!8hHN
zl~G-gg*BO}gYgNUUI+6%G4BrU8`EPs7$3-Tg2X@aFz|>^DD1YIp>&a3<^qch*tf~7
zlQ$%YNP$M_rZd#Aq%dZ?fHNCQEo+`&4HG!CiA)Y~<4~(*O=qZOOJ}HMN6K(au=M4^
z5W9$hp_XHE{zCD3R#0|=w4iG_5gDh3vxW(lV`{llm}|LHm}_}zxLg=wtr-|<c~e+w
z7;AVzdQ+G|B%B8}50rl4B5WXekO(Nv*YeeHGl5#~d=S^9F)%<3r~+jM?iwCsNd|@#
zh;bk}o*r8kh8k`})~(@!7{R~*>JZd&*YMz!<C&cQS&XTMdvg0{h58yEm_9~^o+T`x
zP5{^zh&v!;4Py<m7=+E5&I}iah`_uK;h~U7YQX^jH&+^>3?c|&v!*krGuQH$@Jv=b
zEHgP_vnW4*4I4c41rXxO{59-walsk^sG)+;U}9iks1=&rd04K#MhNCdSTVxDPy{O2
zia@1c5f3O8!-q#uGLRfd3N+PT1a6*#S}u@k7Tg<L2a*MiNN5V*Vu$4kaNY)IByj!o
zosl7e0a`&NRor4qssuN?izb8g-T+mz+-0dnC8-thWr;Zzw~#oM;1MP8AXRn|Xfz4Y
z<G#h5l37+PJUQT)SUtF}S@Z*xK|#b#Q2qjSE*cm<a4~T3c9b;L+~5)HWPHHGf1O9>
zB9F|Bz{@=H9gGh+xf`5rfQUQ1LK6(TBR(<kv&wfc-VhR<t~*h80pmp>r3Ho;g;YD(
zo`|b1NL!J7S=^|D=K~vqxI_o{4KeW!?i+%#pXwQyIdwW5ZwQOcV4R^ilVwKWO!ftg
z3ltZ!EC^i4K9he1;|j%<EGq(6vR{<YzbI_b;c|mVe5UyVli9W_1TRW!U*ys5@Wv48
z^!~uk%&GH%gF#$riSmM+%Sy%<#f(ALDxz>Xxa#{kJ2@w4U*V9u!NCP#U*(YdAkWV#
z|3QU;o2SFEQ>CM@Q+<MBr`CkP3563wW++Y!pAk4QY60T~PNfBo3lbN)tPou2zQS>#
z*NVaof-6flIIb+;P`I*chvEfgn;nS<7<Z-|Fx;7bAaQ5b3BePFCq+&qo)o_j7<5wV
zLP%)91>3L-9N`U)Pk8zJBReDe<2&PL1kNbB%&XMkb_a{t11^CHf}NTZ+-4X~^t{L=
z+u;0!o44Ph(_(^ix9tTMiS3@8jBc!qozo}qG1@Y+PCv)Tn5YDv|IGxY7AE5qO~zYX
z1x5K;smUdo`FTZ<KFIVyenuq+P|xHROJ06nswUel_M+mFqRf=kTkOdYMo}Ti1kTiy
z^wi}1ocyAqB2cOUP1<WR-eO9uC`y{XnV->y@$~c${EW)}Xq_1`NXiFwW^OTO=9LuH
zf^>s>G~mnuD%grZEw3VYx8@d0Vnu3k5u{%;JxPGktR832MzRRepF!jycx7~pF{9`V
z$U@K%R}lwjfL<Atg^1|lfO|K>ObiUS7z@OT;N6^}vmoQZ9URHY0mq~#%gV4%_ZMX3
zXO&@MU|`&yEXa76$+`#B&18Y}Aiz~8IOl`Mf!H9GCb)UbUY3{x9s&hbl8{z%=k!Do
zMk!W#CI$x6=?x-`QyIOci-|HS@ba~bwu`lkw@b83wo6TS6lIiVWS@L7O1Qp+fuT#H
zgS~^LL!yJJgCmVKg|UUBOR|H#gSA7VgBdKs)WXpvg`$SJg`-OxMFiAn6zu>t8~HkT
zI=DM{JGj!=K~2Xlu@2D=#tyL#@eZ~Q$quPB_7t`jj#>do1JVxEIO}BZlmNwKjR34@
zjbUbBs1>Xgsuf05h%i;Ka<x{ZhJW$~Em0-G8j%$CDh>t)Q1u3rWn}0nV_~S_4`$Hh
z@S9%f&FH`jnmYxh=aq~_9~q~g6l1hvJU;!RIHL-q-*ioJ#uPRsa12kMEY4^PlkuPa
zO`I`>7bLTiu?W-xFH)Jl(U(yKrXX<oLJ3Br>GRbXO~9%&8Np2sR8?}pYzz#=pm~rV
z4Gc@-HaKob-l2Mc<AmUapx`U2A)uz-^nOW3GjDL$iV-{#tH}ZBPlLutZ?PsP$LAV?
z$NnJ$te~OrBGB?F@Stfis3w6fkAO}7LVOlF-BgOvo(EKGOl6<Wz{JqWJ~<{&f&(=6
z*2!AMJ$XW&@^k@TM#;%9Ch<&Gt`MBOuAF;vLKh2T)$|5qMp=GPpmBi`J$MY1W4fR;
zqXaJxtdUg~pKCZ>Lz>Y{7?%s!B0yemojyyN(Xu}8LVm&Jw8D!;bsrcQ9GD!LJ~J>B
zGdeMSWMD{#P_kf-D^n=b2L^^Prf`r7SEdM%3fTydSTs`%h$9ojgwg_Kh2$1+p#UDx
z0%ZY6A5*%>2o$eoObqr6MQR|XHi%FK5n3R^5=0=fKPbS#p;`nEWTY|X>5(#wvJQ|A
z=pOJOVup8-DM+_Dh_C<=pq6Qo6^LaGnva$S_dUS_JTPMV0U1U`0S6`q1~*X9fFkk#
z|3A||$}oy~39_npFg_Rd?clgC641f<Tq3-K_qlY)1rCYlBBd8aN;^2evoQ#QLqITw
zse|>7k~*}q3k6kI9P%Af6I^D9KH(Ldo-50!AfznFs{TQZfk&{x2hxCLna(K6D6)OC
zEaPiNZcwT&(gV2_o~HN7Gpe%XvN13e=}&(t&luteR|8wP2M#E3;DL)B@VE&$R6*^J
zVo<*inh+t4gyYjeMlx1TU!uS`TgU`t88{`eKvK8G^gu;MZAPf7+UbiG869Q7W@3yK
zfGts-ro<@911|np$`W%7iaMq%Dlv+xAO>f^vET}_2OJd$N4iaqQere<^=4vVNZQ`7
z#5ju?5&+5|weSFVr_QL#IBj~q8e<BQv?eUigQcfWmsDr$R|I<oJmdmh(+EyTU|S&x
zVERdQMis{6)4`f%P8ZQ&Y(p{=#Y(m%pfLBFzR-lx6G;t<y=+TCY63uNz$?_IPuFA&
zMYy6U2$bNkL;%S8`sE<Kp&+k<ilHJ<mR<?1Vi1`I)YjA#!<I>Iu@z<HgNmSFkP#uE
zij6g+1SA0VG&oU!5&*bDy~SFRnVV_=p$x$@9STK=F*H!tfvn9z4y%h=j3%scObiU&
z+qtwEiy0Xsr`PE)Dlop9K3|8?nl+V)fnn43`#Ow5j7s2v=|qrkLBr9YfgMO%O8~hI
zMC_d|tjFkW3knNR{Q|-uHfV6-^A|=2hUpA7Y$=S{8Vn3Y#h|7PNC5*w4RZ}EtOBZK
znC!SvWcq&tMv>{$4H!cvzuUkyeY+l`{N%IToDQJU6Rwi6hM|TL)Vrzy&yLix=hf6O
z)N<4?)^d_zc0E3`YdC5+L6scHwp#WYc96~K%)tyrp!G+Z%tfF;g0z^xJp`^>EV+pl
zhG5dD2p)4q;M@x?`XJ38P*Da_X?%+%S17T<q!^SsK*0*Z(|gSr#q0B!7#L203M|kd
zcmu-+9tIB14#%dlC;W;F)K*Ad<~QzOxx=rxL}|IkLX8!si}klCZP(bSvBT`Lk<&#b
zrz1wkEe=|o@VM+8dL=C4dRW}Wu(%8H8JEK{FFI#l=FjS2dB86`!KJfuhRQ`gxlZN}
zYzzV-9V|Eag~7%ysHeIeGfEb?%&fe~r_sUufQNqq%Y={)`wJ{mka7|<zY9qf3=E*5
zDrQhp{d|v!fdM++#?sE(F4!*AF5E6ModML>pMFu6QFMB*1*3=|W2`3wLzhqoO9xYj
zPzPfNE3{Q5oWhjCoX(WO0wtN!nbMhS**jS$3&x44bqaynK(!o5Eglw#7_3cG%bCJ5
z`RqbT=Z_2wwcMRh*&6m5Zjeq8u0`sB!{vEkNwrG^F|;;$gM^4s4Ku`o8s-`fh#}M0
zxiX5U3&SK**s?*{j2oc<WD7*JNGOFJrVdpelsrMnv`aA08HWlk28I%DP-h$N2AC%y
z!b~;Hlh5CkumTM%bP09{b_jQfbTFi`q;Rxw)bb(*{c3n>cwyd)0W})=YWZsgYPd?c
zCwCu~5vqX>axgGJG7dyorbYl}8_0a{zz0{)XBLJUwk*-f4;P8ZrZ9urlQ2b$3=<f8
zd_dB{44Rza)}$l0Lb}VAQI+F0sG!whVwf!OS4y!6)GICmC7@*JSSvVDfkyB^rRC?S
z$pRBKK<zJXmQLo?Op`xWh{$y^f{Fl8_@yv{ToKGr1g>eo)h6rodOJo@g(6Vn0z6EA
zi@7Y`05S^<S{(~12#-$(Ii>y|D7$JgFfjaRV3?CMKW%2(!r&Exi^EoUT$Te*d~L2d
z;CInH=%QTkg|P68vf<Dsl_ob*YYJQ$gIiPJ8X4Rl0T;mFnL+RjizcXj2(CLsi$F`x
zi@*&mt$a{%45}(1eHM_1Z*e8(=at2mq!#5CgU9%bz{cByRDu^WfU9fpoC~%Af6&+o
zcw7VF|6<ViBxn>cAz|{u|Dv{hpg~Y1OG=p-7#@MbVJag7!~Y-Pf=fuG!RrpM&=lhd
zg_F&?V?Q%Uu_|>iJ{Jj?Ja4IVy@>d9--*5p6fcUXbZ|c56`Y_r!EvHWcT@+{13vx^
z<~xF76AC9wd}83^H1BY{As{$GaU#nEaPR)2fI^4k6LB^0NN&jr$IIe+9UfoU7zBj+
z%R0*@q|FGNn7M#$LFER)g*6uyOfHI<T;@03!FGw??tt0}!GoG7(k?lL)L)1I6^s`$
zvo6MGUyR7PkY9Q!qV$HM%7U~NfeSM?ux$X%Sk_!LaJi`Fa#7KBBHsfgjTM3mHCMQ8
zu-p;2(dMF#!$l>BBSy!~4w{`%1C1GIUUm+-5*mIzH2Pv_^o5wD3&}Yba`P@G=U<E|
zxExw|(Yf%VV$lS?XF}l}&Ua#zHC70&FkWnPkw>v3VnWIc(+52KQ;epYPc)xlaamYp
zLC{5E^)*cE*;le}U|r33QP}Jvk6EYl2U$K&^A8FPoWj>RBrkGE&QP9FI8*H!hvElQ
zDORNq77RRM7g$86zjkMAR=5V5Rb&KDerYlnfhMtvK(lB?potkxj-ue{eIASw+~96F
zXi})?)%1-XjGC<FObiSh+wXfY&WAQ#QCn{Fyco6F6xkUViqt_(B&^LSuo@+h8XY7x
zs7+w78WoTleRMVO<`_r~qw4e-K8$y4L77vN8B*;cchtc{1)$dGEk5u#9b`TrzQ7RC
zJt!)MwJX5}s!#vq%NQCF$;7~*DTv54kURx${(>_WqWZQ0<s0y>4LcAE+z1A@u0Yln
z)q%u7^%m5^>GS;<O<C)i7#Jj`-}GZFX7rx!>(A)QST}vLKcjg)cpDOEn-pZ*4rqf>
z5oi}t5qK*PXh%>HcvKR!iwCki3AED&vLp_?W*gk^0v8CN%^yXe?Hom*3GyP)d>Uk?
z3o#`Co(KSq$rpjvz7=^g)-x~^fmbGh7I77U7DN?+)-!?o!bPA(Lq(u9K}Dc7Jw>2(
zIz@<;AE2cykZE(!JTW9KfhV%TNeh&!AY(+}?mK8g05Wt58geO;2L%$SGv-qS>KPP)
zTA+}26sRFt1ZwpZfy$5~Q0glJ1zr)T$q#NNfD7AS95%W6DWy57c15kAF&qX4hT`1m
zl7Wmy(}MyTe|ZWp%6?K)Wt99Bz^Kb8@F|0lk5TxOk|d+xCkr)3!B0VqDvXk!{1`bH
i6+Ve^F)DnL5@nSBWT4C__Q{8lmr>%gf+(XP*fsz%sJ-F<
diff --git a/2d/compute_errorestimates.py b/2d/compute_errorestimates.py
index c4e7ade..eebd88b 100644
--- a/2d/compute_errorestimates.py
+++ b/2d/compute_errorestimates.py
@@ -14,11 +14,10 @@ method = 'wasserstein' # using log-mean convective flux
boundary = 'Neumann'
interpol = 'morley'
-test = 'manufactured1' # 'diff', 'blow-up', 'manufactured' or 'manufactured1'
-fineness = 0 # number of refinements of mesh before calculating numerical solution
+test = 'manufactured1' # 'diff', 'blow-up', 'manufactured' or 'manufactured1'
+fineness = 1 # number of refinements of mesh before calculating numerical solution
Nt = 8 # number of time steps
-
# initialize different test cases
if test == 'blow-up' :
# CONVECTION DEOMINATED REGIME i.e. BLOW-UP
@@ -94,6 +93,8 @@ for i in range(fineness) :
K_dual = np.array(K_dual)
[tri_dual,K_dual] = my.orientation_dualmesh(tri_dual,K_dual) # enforce positive orientation of all triangles in dual mesh
+[K_inter,tri_inter,points_inter,TinK] = my.intersect_mesh(points,K,tri)
+
# compute mesh topology for faster find_simplex routine for dual mesh
[oriented_edges,structure,edge_to_index] = my.get_edges_structure(K_dual)
my.init_indexE(numtri_dual)
@@ -109,6 +110,9 @@ print('Nt',Nt)
max1 = 0 # need this for more efficiency calculating theta_omega
max2 = 0 # need this for more efficiency calculating theta_omega
+
+toc = time.time()
+
for n in range(Nt) : # range(Nt) :
progress = n/(Nt)*100
@@ -122,9 +126,6 @@ for n in range(Nt) : # range(Nt) :
# compute c related objects
v_0 = lambda x: my.grad_cfunc(cc_0,tri_dual,K_dual,points_dual,x,oriented_edges,structure,edge_to_index) # piecewiese constant
- [lin_vx_0,lin_vy_0] = my.get_linv(K_dual,nppoints_dual,numtri_dual,tri_dual,v_0,nodal_avrg_choice) # get linear coefficients
- Laplace_c_0 = lambda x : my.approx_Laplace_cn(lin_vx_0,lin_vy_0,x,oriented_edges,structure,edge_to_index) # approximation of Laplace(c)=div(v)
-
v_m1 = lambda x : v_0(x) # as morley0 uses c0 instaed of c^(n-1) in n = 0 case
# load data
@@ -142,15 +143,15 @@ for n in range(Nt) : # range(Nt) :
ht_m1 = 0
# COMPUTE A POSTERIORI ERROR ESTIMATES
- theta_early = my.get_early_time_error(K,tri,v_0,qq0_0,qq0_p1,betaKE_0,betaKE_p1)
+ theta_early = my.get_early_time_error(K,tri,v_0,qq0_0,qq0_p1,betaKE_0,betaKE_p1)
print('extra term for early times computed.')
- theta_res0 = my.get_theta_res(cc_0,Laplace_c_0,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,K_dual,tri_dual,points_dual,oriented_edges,structure,edge_to_index)
+ theta_res0 = my.get_theta_res(cc_0,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,K_dual,tri_dual,points_dual,oriented_edges,structure,edge_to_index,TinK)
print('elementwise residual computed.')
- theta_diff0 = my.get_theta_diff(qq0_0,betaKE_0,K,tri)
+ theta_diff0 = my.get_theta_diff(qq0_0,betaKE_0,K,tri)
print('grad morley jump terms computed.')
[theta_time0,theta_time1] = my.get_theta_time(qq0_0,qq0_p1,betaKE_0,betaKE_p1,rho_m1,rho_0,rho_p1,ht_m1,ht_0,K,tri,n)
print('time contributions computed.')
- [theta_Omega,max1,max2] = my.get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,n)
+ [theta_Omega,max1,max2] = my.get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,qq0_0,qq0_p1,K,tri,n)
print('global convective term computed.')
# SAVE A POSTERIORI ERROR ESTIMATES
@@ -171,7 +172,6 @@ for n in range(Nt) : # range(Nt) :
# time-update objects
cc_m1 = cc_0
v_m1 = v_0
- Laplace_c_m1 = Laplace_c_0
rho_m1 = rho_0
rho_0 = rho_p1
qq0_0 = qq0_p1
@@ -187,8 +187,6 @@ for n in range(Nt) : # range(Nt) :
# compute c related objects
v_0 = lambda x: my.grad_cfunc(cc_0,tri_dual,K_dual,points_dual,x,oriented_edges,structure,edge_to_index) # piecewiese constant
- [lin_vx_0,lin_vy_0] = my.get_linv(K_dual,nppoints_dual,numtri_dual,tri_dual,v_0,nodal_avrg_choice) # get linear coefficients
- Laplace_c_0 = lambda x : my.approx_Laplace_cn(lin_vx_0,lin_vy_0,x,oriented_edges,structure,edge_to_index)
# dummy variables
betaKE_m1 = 0*np.array(betaKE_0)
@@ -200,7 +198,6 @@ for n in range(Nt) : # range(Nt) :
# time-update objects
v_m2 = v_m1
v_m1 = v_0
- Laplace_c_m1 = Laplace_c_0
rho_m1 = rho_0
rho_0 = rho_p1
qq0_0 = qq0_p1
@@ -217,20 +214,17 @@ for n in range(Nt) : # range(Nt) :
# compute c related objects
v_0 = lambda x: my.grad_cfunc(cc_0,tri_dual,K_dual,points_dual,x,oriented_edges,structure,edge_to_index) # piecewiese constant
- [lin_vx_0,lin_vy_0] = my.get_linv(K_dual,nppoints_dual,numtri_dual,tri_dual,v_0,nodal_avrg_choice) # get linear coefficients
- Laplace_c_0 = lambda x : my.approx_Laplace_cn(lin_vx_0,lin_vy_0,x,oriented_edges,structure,edge_to_index)
-
# COMPUTE A POSTERIORI ERROR ESTIMATES
- theta_res0 = my.get_theta_res(cc_m1,Laplace_c_m1,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,K_dual,tri_dual,points_dual,oriented_edges,structure,edge_to_index)
+ theta_res0 = my.get_theta_res(cc_m1,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,K_dual,tri_dual,points_dual,oriented_edges,structure,edge_to_index,TinK)
print('elementwise residual computed.')
theta_diff0 = my.get_theta_diff(qq0_0,betaKE_0,K,tri)
print('grad morley jump terms computed.')
[theta_time0,theta_time1] = my.get_theta_time(qq0_0,qq0_p1,betaKE_0,betaKE_p1,rho_m1,rho_0,rho_p1,ht_m1,ht_0,K,tri,n)
print('time contributions computed.')
- conv_terms = my.get_conv_terms(v_0,qq0_p1,betaKE_p1,rho_p1,K,tri)
+ [conv_terms,jump_cs] = my.get_conv_terms(v_0,qq0_p1,betaKE_p1,rho_p1,K,tri,TinK)
print('convective term computed.')
- [theta_Omega,max1,max2] = my.get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,n)
+ [theta_Omega,max1,max2] = my.get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,qq0_0,qq0_p1,K,tri,n)
print('global convective term computed.')
# SAVE A POSTERIORI ERROR ESTIMATES
@@ -240,10 +234,16 @@ for n in range(Nt) : # range(Nt) :
data.append(theta_time0)
data.append(theta_time1)
data.append(conv_terms)
+ data.append(jump_cs)
data.append(theta_Omega)
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
pickle.dump(data,open('pickle_files/'+pickle_name,'wb')) # store data
+ elapsed = time.time() - toc
+ print('This step took',"%.2f" % round(elapsed/60, 2), 'minutes.')
+ toc = time.time()
+
+
progress = 100
print("%.2f" % progress, 'procent of progress made.')
elapsed = time.time() - tic
diff --git a/2d/main.py b/2d/main.py
index 929dc8a..645d31a 100644
--- a/2d/main.py
+++ b/2d/main.py
@@ -16,7 +16,7 @@ boundary = 'Neumann'
interpol = 'morley'
test = 'manufactured1' #'diff', 'blow-up', 'manufactured' or 'manufactured1'
-fineness = 3 # number of refinements of mesh before calculating numerical solution
+fineness = 1 # number of refinements of mesh before calculating numerical solution
Nt = 8 # number of time steps
print('FV-FE Algorithm')
diff --git a/2d/myfun.py b/2d/myfun.py
index d16221f..b3eedfe 100644
--- a/2d/myfun.py
+++ b/2d/myfun.py
@@ -108,6 +108,73 @@ def circumcenter(K) :
return CC
+def intersect_mesh(points,K,tri) :
+# input: primal mesh
+# output: intersected mesh of primal and dual mesh
+
+ inter_points = []
+ primal_points = len(points)
+
+ for pt in points :
+ inter_points.append(pt)
+
+
+ CCs = []
+ for i in range(len(K)) :
+ CC = circumcenter(K[i])
+ CCs.append(CC)
+ inter_points.append(CC)
+
+ circum_points = len(CCs)
+
+ class myclass :
+ def __init__(self) :
+ self.simplices = []
+ tri_inter = myclass()
+
+ TinK = [] # get triangles of intersected mesh that are in primal element K
+ midEs = []
+ K_inter = []
+ for i in range(len(K)) : # go through triangles
+ if i == 0 :
+ aux = []
+ else :
+ TinK.append(np.array(aux))
+ aux = []
+
+ indices = [[1,2],[0,2],[0,1]]
+ for j in range(3) : # go through edges
+
+ E = K[i][indices[j]] # get edge
+ E_index = tri.simplices[i][indices[j]] # get indices of edge
+ midE = 1/2*(E[0]+E[1]) # compute midpoint of edge E
+
+ if midE[0] == 0 or midE[0] == 1 or midE[1] == 0 or midE[1] == 1: # E is boundary edge
+ K_inter.append([CCs[i],E[0],E[1]]) # only one triangle for boundary edges
+ tri_inter.simplices.append([primal_points+i,E_index[0],E_index[1]])
+ aux.append([CCs[i],E[0],E[1]])
+
+ else : # interior edge
+ k = linear_search(midEs, midE) # check if midE was already used
+ if k < 0 :
+ midEs.append(midE) # update list of midEs
+ inter_points.append(midE)
+ k = len(midEs)-1 # get index for midE
+
+ K_inter.append([CCs[i],E[0],midE]) # add triangle one
+ tri_inter.simplices.append([primal_points+i,E_index[0],primal_points+circum_points+k])
+ aux.append([CCs[i],E[0],midE])
+
+ K_inter.append([CCs[i],E[1],midE]) # add triangle two
+ tri_inter.simplices.append([primal_points+i,E_index[1],primal_points+circum_points+k])
+ aux.append([CCs[i],E[1],midE])
+
+ TinK.append(np.array(aux))
+ points_inter = np.array(inter_points)
+
+ return [K_inter,tri_inter,points_inter,TinK]
+
+
def are_neighbors(Ki,Li) :
# need this for function getdualmesh(...)
@@ -1161,7 +1228,7 @@ def finitevolumescheme_rho(u_old,ht,numtri,tri,K,eps,v,f) :
[vKE,areaE] = getvKE(K[i],E,v) # avrg over v*nKE ( = 0 if E on boundary due to neumann boundary w.r.t. c )
dE = getdistE(K[i],K[neighbors[j]]) # get distance between circumcenters of triangles
- if np.abs(u[i]-u[neighbors[j]]) < 10**-6 or u[i] < 10**-6 or u[neighbors[j]] < 10**-6 : # in this case flux=0 to avoid divide by 0 or log(0)
+ if np.abs(u[i]-u[neighbors[j]]) < 10**-10 or u[i] < 10**-10 or u[neighbors[j]] < 10**-10 : # in this case flux=0 to avoid divide by 0 or log(0)
sum_EK[i] = sum_EK[i] - eps*areaE/dE*(u[neighbors[j]]-u[i])
else :
sum_EK[i] = sum_EK[i] - eps*areaE/dE*(u[neighbors[j]]-u[i]) + vKE*areaE*((u[i]-u[neighbors[j]])/(np.log(u[i])-np.log(u[neighbors[j]])))
@@ -1199,7 +1266,7 @@ def finitevolumescheme_rho(u_old,ht,numtri,tri,K,eps,v,f) :
[vKE,areaE] = getvKE(K[i],E,v) # avrg over v*nKE
dE = getdistE(K[i],K[neighbors[j]]) # get distance between circumcenters of triangles
- if np.abs(u[i]-u[neighbors[j]]) < 10**-6 or u[i] < 10**-6 or u[neighbors[j]] < 10**-6 : # in this case flux=0 to avoid divide by 0 or log(0)
+ if np.abs(u[i]-u[neighbors[j]]) < 10**-10 or u[i] < 10**-10 or u[neighbors[j]] < 10**-10 : # in this case flux=0 to avoid divide by 0 or log(0)
valK = valK + eps*(areaE/dE)
valL = valL - eps*(areaE/dE)
else :
@@ -1475,7 +1542,7 @@ def grad_morley_interpol(K,tri,qq0,betaKE,x,indexK) :
return np.array([qq0[3*i],qq0[3*i+1]])+II # piecewise grad_morley = grad_q0 + beta*(...) = [a,b]+beta*(...)
-def Laplace_morley(K,tri,qq0,betaKE,x,indexK) :
+def Laplace_morley(K,tri,betaKE,x,indexK) :
# input : triangulation K,tri / coefficients qq0, betaKE / point x and index of the triangles in which is lives if known
# output : fucntion value of the Laplacian of morley interpolant at point x
@@ -1487,51 +1554,9 @@ def Laplace_morley(K,tri,qq0,betaKE,x,indexK) :
II = 0
for j in range(3) : # go through edges
- aux = lambda x : [qq0[3*i],qq0[3*i+1]] # grad_q0
- div_grad_q0 = lambda x : get_second_derivative(aux,K,tri,x,index=i) # divergence of the component-wise linear interpolation of the piecewise constant gradient grad_q0
II += betaKE[i][j]*(LaplacebubbleE(K[i],j,x)*bubbleK(K[i],x)+LaplacebubbleK(K[i],x)*bubbleE(K[i],j,x)+2*np.dot(gradbubbleE(K[i],j,x),gradbubbleK(K[i],x))) # product rule on beta contributions
- return div_grad_q0(x)+II # piecewise Laplace_morley = "Laplace" q0 + alpha*(...)+beta*(...)
-
-
-def get_second_derivative(grad_q0,K,tri,x,index) :
-# need this for Laplacian of Morley, which is in turn needed to compute element-wise reisduals for the a posteriori error estimates
-# input : piecewise constant grad_q0 / triangulation K,tri / point x and the index of the triangle in which it lives if known
-# output: divergence of linear interpolation of piecewise constant grad_q0
-
- if index > -0.5 : # if primal index
- i = index
- else :
- i = tri.find_simplex(x)
-
- bx = []
- by = []
- for j in range(3) :
- pindex = tri.simplices[i][j] # get index of vertex
-
- neighborstri = find_neighbors(pindex, tri) # indices of triangles that touch vertex associated with pindex
-
- # get nodal values
- vertex_valx = 0
- vertex_valy = 0
- for k in neighborstri : # go thorugh neighboring points
- midk = 1/3*(K[k][0]+K[k][1]+K[k][2]) # middle of triangle, as function is piecewise constant it doesent matter where we evaluate
- vertex_valx += grad_q0(midk)[0] # x-direction
- vertex_valy += grad_q0(midk)[1] # y-direction
- vertex_valx = vertex_valx/len(neighborstri) # arithmetic mean
- vertex_valy = vertex_valy/len(neighborstri) # arithmetic mean
- bx.append(vertex_valx)
- by.append(vertex_valy)
-
- # get linear system of equation to solve for linear inteprolation a*x+b*y+c = val
- M = [[K[i][0][0],K[i][0][1],1],[K[i][1][0],K[i][1][1],1],[K[i][2][0],K[i][2][1],1]]
-
- lin_vx = np.linalg.solve(M,bx)
- lin_vy = np.linalg.solve(M,by)
-
- div = lin_vx[0]+lin_vy[1] # get divergence of linear inteprolation
-
- return np.array(div)
+ return II # piecewise Laplace_morley = Laplace q0 + alpha*(...)+beta*(...), where Laplace q0 = 0 on each intersected element
def getmorleyplot(n,ht,K,tri,qq0,betaKE,h,rholim,fineness) :
@@ -1603,7 +1628,7 @@ def diam(K) :
return val_max
-def get_theta_res(cc_m1,Laplace_v_m1,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,K_dual,tri_dual,points_dual,oriented_edges,structure,edge_to_index) :
+def get_theta_res(cc_m1,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,K_dual,tri_dual,points_dual,oriented_edges,structure,edge_to_index,TinK) :
# get local residual for each triangle K
# input: c^n-1, Laplace(c^n-1), morley^n, morley^n+1, time step size ht_0, primal and dual triangulation, objects to make find_simplex for dual mesh more efficient
# output: the local residual part of the a posteriori error estimator
@@ -1616,13 +1641,17 @@ def get_theta_res(cc_m1,Laplace_v_m1,qq0_0,qq0_p1,betaKE_0,betaKE_p1,ht_0,K,tri,
grad_morley_0 = lambda x: grad_morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i) # get grad_morley^n
- Laplace_morley_0 = lambda x: Laplace_morley(K,tri,qq0_0,betaKE_0,x,indexK=i) # get Laplace_morley^n
+ Laplace_morley_0 = lambda x: Laplace_morley(K,tri,betaKE_0,x,indexK=i) # get Laplace_morley^n
grad_c_m1 = lambda x : grad_cfunc(cc_m1,tri_dual,K_dual,points_dual,x,oriented_edges,structure,edge_to_index) #get grad_c^n-1
# GET L2 NORM OF LOCAL RESIDUAL
- aux_0 = lambda x : ((morley_p1(x)-morley_0(x))/ht_0 + np.dot(grad_morley_0(x),grad_c_m1(x)) + morley_0(x)*Laplace_v_m1(x) - Laplace_morley_0(x))**2
- [integral_0,avrg,areaK] = avrgK(K[i],aux_0)
+ aux_0 = lambda x : ((morley_p1(x)-morley_0(x))/ht_0 + np.dot(grad_morley_0(x),grad_c_m1(x)) - Laplace_morley_0(x))**2
+
+ integral_0 = 0
+ for T in TinK[i] :
+ [integral,avrg,areaK] = avrgK(T,aux_0)
+ integral_0 += integral
hK = diam(K[i]) # get diameter of K
@@ -1735,7 +1764,7 @@ def get_max_q0_beta(rhoFV,betaKE) :
return [nodal_max,beta_max]
-def get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,n) :
+def get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,qq0_0,qq0_p1,K,tri,n) :
# get theta_omega term from masterthesis
# input: max1 and max2 that were already calculated in the previous time step, morley^n, morley^n+1, FV sol rho^n, FV sol rho^n+1
@@ -1744,57 +1773,122 @@ def get_theta_omega(max1,max2,rho_0,rho_p1,betaKE_0,betaKE_p1,n) :
if n == 0 : # slightly different situation at early times
max1 = get_maximum_morley(rho_p1,betaKE_p1) # no difference
- max2 = get_maximum_morley(np.array(rho_0)-np.array(rho_p1),np.array(betaKE_0)-np.array(betaKE_p1)) # difference
+ infnorm = get_maximum_morley(np.array(rho_0)-np.array(rho_p1),np.array(betaKE_0)-np.array(betaKE_p1)) # difference
+
+ val = []
+ for i in range(len(K)) :
+ morley_0 = lambda x: morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i)
+ morley_p1 = lambda x: morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i)
+
+ aux = lambda x : (morley_p1(x) - morley_0(x))**2
- val = (max1+max2)*max2 # = theta_omega
+ [max2,avrg,areaK] = avrgK(K[i],aux)
+ max2 = np.sqrt(max2)
+
+ val.append((max1+infnorm)*max2) # = theta_omega
else :
+
max3 = max1 # update terms from previous time step
max4 = max2
- max1 = get_maximum_morley(rho_p1,betaKE_p1) # no difference
- max2 = get_maximum_morley(np.array(rho_0)-np.array(rho_p1),np.array(betaKE_0)-np.array(betaKE_p1)) # difference
+ max1 = get_maximum_morley(rho_p1,betaKE_p1) # no difference
+ infnorm = get_maximum_morley(np.array(rho_0)-np.array(rho_p1),np.array(betaKE_0)-np.array(betaKE_p1)) # difference
+
+ val = []
+ for i in range(len(K)) :
+ morley_0 = lambda x: morley_interpol(K,tri,qq0_0,betaKE_0,x,indexK=i)
+ morley_p1 = lambda x: morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i)
- val = (max1+max2)*max2+max3*max4 # = theta_omega
+ aux = lambda x : (morley_p1(x) - morley_0(x))**2
+
+ [max2,avrg,areaK] = avrgK(K[i],aux)
+ max2 = np.sqrt(max2)
+
+ val.append((max1+infnorm)*max2+max3*max4) # = theta_omega
return [val,max1,max2] # max1 and max2 will be used to calculate theta_omega for the next time step
-def get_conv_terms(v_0,qq0_p1,betaKE_p1,rho_p1,K,tri) :
-# get convection contributions for each triangle K
+def get_conv_terms(v_0,qq0_p1,betaKE_p1,rho_p1,K,tri,TinK) :
+# get convection contributions for each triangle K and the normal gradient jump term of the FE solution
-# input: v^n=grad_c^n, morley^n+1, FV sol rho^n+1, primal triangulation
-# output: convection contributions to the a posteriori error estimator
+# input: v^n=grad_c^n, morley^n+1, FV sol rho^n+1, primal triangulation, intersected triangles T that are in primal element K TinK
+# output: convection contributions to the a posteriori error estimator and the normal gradient jump term of the FE solution
conv_terms = []
- for i in range(len(K)) :
+ jump_terms = []
+ for i in range(len(K)) : # go through primal mesh
+ morley_p1 = lambda x: morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i)
+
+ neighbors = tri.neighbors[i] # get neighbors in primal mesh
+ indices = [[1,2],[0,2],[0,1]]
+
+ val1 = 0
+ val2 = 0
+ for l in range(len(TinK[i])): # go through intersected elements in primal element K[i]
+ for j in range(3) : # go through intersected edges of intersected element TinK[i][l]
+ S = TinK[i][l][indices[j]] # get edge in intersected mesh
+
+ normalTS = unitouternormal(TinK[i][l],S) # get normal vectors of intersected mesh
+
+ midS = 1/2*(S[0]+S[1])
+ if midS[0] == 0 or midS[1] == 0 or midS[0] == 1 or midS[1] == 1 : # S on boundary
+ jump_c = np.dot(v_0(midS),normalTS) # get boundary jump term, works like this as grad_c is piecewise constant
+
+ # jump terms are relevant on these edges
+ aux = lambda x : (jump_c*morley_p1(x))**2
+ [integral,avrg,areaS] =avrgE(S,aux)
+ val2 += integral*areaS
+
+ else : # S in interior
+ k = isSonE(S,K[i]) # check if S is on primal edge E and if so, which index k corresponds to edge E, else its on an dual edge
+ if k > -0.5: # S is part of an interior primal edge E
+ if np.abs(rho_p1[i]-rho_p1[neighbors[k]]) < 10**-10 or rho_p1[i] < 10**-10 or rho_p1[neighbors[k]] < 10**-10 : # in this case set flux to 0 to avoid divide by 0 or log(0)
+ F_E_over_areaE = 0
+ else :
+ F_E_over_areaE = (rho_p1[i]-rho_p1[neighbors[k]])/(np.log(rho_p1[i])-np.log(rho_p1[neighbors[k]]))
+
+ aux = lambda x : (np.dot(v_0(x),normalTS)*(morley_p1(x)-F_E_over_areaE))**2 # get L2 norm
+ [integral,avrg,areaS] =avrgE(S,aux)
+ val1 += integral*areaS
- morley_p1 = lambda x: morley_interpol(K,tri,qq0_p1,betaKE_p1,x,indexK=i) # get morley^n+1
+
+ else : # S is part of an interior dual edge
+ jump_c = np.dot(v_0(midS+10**(-8)*normalTS)-v_0(midS-10**(-8)*normalTS),normalTS) # get interior jump term, works like this as grad_c is piecewise constant
- indices = [[1,2],[0,2],[0,1]]
- neighbors = tri.neighbors[i] # get nieghbors of triangle subject to index i
+ # jump terms are relevant on these edges
+ aux = lambda x : (jump_c*morley_p1(x))**2
+ [integral,avrg,areaS] =avrgE(S,aux)
+ val2 += integral*areaS
+
+ conv_terms.append(val1)
+ jump_terms.append(val2)
- hK = diam(K[i]) # get approximated diameter of the triangle
+ return [conv_terms,jump_terms]
+
- val = 0
- for j in range(3) : # go through all edges
- E = K[i][indices[j]]
+def isSonE(S,K) :
+# input: intersected edge S, primal element K
+# output: index k of primal edge if S lies in primal edge E \subset \partial K, and k=-1 if not
+# check if S lies in some edge of primal mesh
- normalKE = unitouternormal(K[i],E) # get normal vector subject to edge E
+ indices = [[1,2],[0,2],[0,1]]
- # COMPUTE convective term F_E/areaE
- if np.abs(rho_p1[i]-rho_p1[neighbors[j]]) < 10**-6 or rho_p1[i] < 10**-6 or rho_p1[neighbors[j]] < 10**-6 : # in this case centered flux to avoid divide by 0 or log(0)
- F_E_over_areaE = 1/2*(rho_p1[i]+rho_p1[neighbors[j]])
- else : # log-mean flux
- F_E_over_areaE = (rho_p1[i]-rho_p1[neighbors[j]])/(np.log(rho_p1[i])-np.log(rho_p1[neighbors[j]]))
+ midEs = []
+ for k in range(3) :
+ E = K[indices[k]]
+ midE = 1/2*(E[0]+E[1])
+ midEs.append(np.array(midE))
- aux = lambda x : (np.dot(v_0(x),normalKE)*(morley_p1(x)-F_E_over_areaE))**2 # edge L2 norm
- [integral,avrg,areaE] =avrgE(E,aux)
+ for k in range(3) : # go through primal edges
+ if np.linalg.norm(S[0] - midEs[k]) < 10**-10 or np.linalg.norm(S[1] - midEs[k]) < 10**-10 : # check if one of the endpoints of S is a midpoint of E
+ for j in range(3) : # go through vertices of K
+ if np.linalg.norm(S[0] - K[j]) < 10**-10 or np.linalg.norm(S[1] - K[j]) < 10**-10 : # check if one of the endpoints of S is an endpoint of E
+ return k
+
+ return -1
- val += np.sqrt(integral)*hK/np.sqrt(areaE) # sum over edges
-
- conv_terms.append(val**2)
- return conv_terms
def get_morley0_coeff(rho_0,nppoints,numtri,tri,K) :
diff --git a/2d/plot_aposti_estimator.py b/2d/plot_aposti_estimator.py
index 2be7a9a..3de8463 100644
--- a/2d/plot_aposti_estimator.py
+++ b/2d/plot_aposti_estimator.py
@@ -63,11 +63,11 @@ if n == 0 : # slightly different situation at early times
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
[theta_early,theta_diff0,theta_time0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
- [theta_res1,theta_diff1,theta_time1,theta_conv1,theta_conv1,theta_Omega1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res1,theta_diff1,theta_time1,theta_conv1,theta_conv1,jump_cs1,theta_Omega1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
# compute \theta_K for each element K
for i in range(len(K)) :
- theta_K.append(theta_early[i]+theta_res1[i]+theta_diff0[i]+theta_diff1[i]+theta_time0[i]+theta_conv1[i])
+ theta_K.append(theta_early[i]+theta_res1[i]+theta_diff0[i]+theta_diff1[i]+theta_time0[i]+theta_conv1[i]+jump_cs1[i])
@@ -75,13 +75,13 @@ else : # n=/=0
# load a posteriori error estimates
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
- [theta_res0,theta_diff0,theta_time0,theta_conv0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res0,theta_diff0,theta_time0,theta_conv0,theta_conv0,jump_cs0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
- [theta_res1,theta_diff1,theta_time1,theta_conv1,theta_conv1,theta_Omega] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res1,theta_diff1,theta_time1,theta_conv1,theta_conv1,jump_cs1,theta_Omega] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
# compute \theta_K for each element K
for i in range(len(theta_res1)) :
- theta_K.append(theta_res0[i]+theta_res1[i]+theta_diff0[i]+theta_diff1[i]+theta_time0[i]+theta_time1[i]+theta_conv0[i]+theta_conv1[i])
+ theta_K.append(theta_res0[i]+theta_res1[i]+theta_diff0[i]+theta_diff1[i]+theta_time0[i]+theta_time1[i]+theta_conv0[i]+theta_conv1[i]+theta_Omega0[i]+jump_cs0[i]+jump_cs1[i]) # get theta_K as sum over all K
theta_K = np.array(theta_K)
diff --git a/2d/plot_meshes.py b/2d/plot_meshes.py
index dcceb29..42755ae 100644
--- a/2d/plot_meshes.py
+++ b/2d/plot_meshes.py
@@ -30,6 +30,7 @@ vor = Voronoi(nppoints)
K_dual = np.array(K_dual)
[tri_dual,K_dual] = my.orientation_dualmesh(tri_dual,K_dual) # enforce positive orientation of all triangles in dual mesh
+[K_inter,tri_inter,points_inter,TinK] = my.intersect_mesh(points,K,tri)
# PLOT THE PRIMAL MESH
plt.triplot(nppoints[:,0], nppoints[:,1], tri.simplices,color='steelblue')
@@ -41,6 +42,10 @@ plt.show()
# PLOT THE DUAL MESH
plt.triplot(nppoints_dual[:,0], nppoints_dual[:,1], tri_dual.simplices,color='green')
+plt.show()
+
+# PLOT THE INTERSECTED MESH
+plt.triplot(points_inter[:,0], points_inter[:,1], tri_inter.simplices,color='orange')
# plt.savefig('test.png', bbox_inches='tight')
diff --git a/2d/scaling_space_errorestimates.py b/2d/scaling_space_errorestimates.py
index 25b6a43..dd07b79 100644
--- a/2d/scaling_space_errorestimates.py
+++ b/2d/scaling_space_errorestimates.py
@@ -17,10 +17,10 @@ method = 'wasserstein' # using log-mean convective flux
boundary = 'Neumann'
interpol = 'morley'
-test = 'manufactured' #'diff', 'blow-up', 'manufactured' or 'manufactured1'
+test = 'manufactured1' #'diff', 'blow-up', 'manufactured' or 'manufactured1'
Nt = 8 # fixed number of time steps
-steps = 5 # number of space refinements used
+steps = 2 # number of space refinements used
# initialize different test cases
if test == 'blow-up' :
@@ -93,6 +93,7 @@ theta_diff_max = []
theta_time_max = []
theta_time_val_max = []
theta_conv_max = []
+theta_jumpcs_max = []
theta_Omega_max = []
@@ -113,6 +114,7 @@ for fineness in range(steps) : # varibale in space
theta_time_sum = 0
theta_time_val_sum = 0
theta_conv_sum = 0
+ theta_jumpcs_sum = 0
theta_Omega_sum = 0
for n in range(Nt) : # go through intervals inbetween time steps
@@ -123,7 +125,7 @@ for fineness in range(steps) : # varibale in space
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
[theta_early,theta_res0,theta_diff0,theta_time0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
- [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_Omega1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,jump_cs1,theta_Omega1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
# initialize objects
theta_early_aux = 0
@@ -132,7 +134,8 @@ for fineness in range(steps) : # varibale in space
theta_time_aux = 0
theta_time_val_aux = 0
theta_conv_aux = 0
- theta_conv_aux = 0
+ theta_jumpcs_aux = 0
+ theta_omega_aux = 0
# sum over all elements
for i in range(len(theta_early)) :
@@ -141,6 +144,8 @@ for fineness in range(steps) : # varibale in space
theta_diff_aux += theta_diff0[i]+theta_diff1[i]
theta_time_aux += theta_time0[i]
theta_conv_aux += theta_conv1[i]
+ theta_jumpcs_aux += jump_cs1[i]
+ theta_omega_aux += theta_Omega0[i]**2
# sum over time steps, corresponds to L^2 norm in time
theta_early_sum += ht*theta_early_aux
@@ -148,13 +153,14 @@ for fineness in range(steps) : # varibale in space
theta_diff_sum += ht*theta_diff_aux
theta_time_sum += ht*theta_time_aux
theta_conv_sum += ht*theta_conv_aux
- theta_Omega_sum += ht*theta_Omega0**2
+ theta_jumpcs_sum += ht*theta_jumpcs_aux
+ theta_Omega_sum += ht*theta_omega_aux
elif n == Nt-1 :
# load a posteriori error estimates
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
- [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,jump_cs0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
# initialize objects
theta_res_aux = 0
@@ -162,7 +168,9 @@ for fineness in range(steps) : # varibale in space
theta_time_aux = 0
theta_time_val_aux = 0
theta_conv_aux = 0
+ theta_jumpcs_aux = 0
theta_delta_aux = 0
+ theta_omega_aux = 0
# sum over all elements
for i in range(len(theta_res1)) :
@@ -171,6 +179,8 @@ for fineness in range(steps) : # varibale in space
theta_time_aux += theta_time0[i]
theta_time_val_aux += theta_time_val0[i]
theta_conv_aux += theta_conv0[i]
+ theta_jumpcs_aux += jump_cs0[i]
+ theta_omega_aux += theta_Omega0[i]**2
# sum over time steps, corresponds to L^2 norm in time
theta_res_sum += ht*theta_res_aux
@@ -178,15 +188,16 @@ for fineness in range(steps) : # varibale in space
theta_time_sum += ht*theta_time_aux
theta_time_val_sum += ht*theta_time_val_aux
theta_conv_sum += ht*theta_conv_aux
- theta_Omega_sum += ht*theta_Omega0**2
+ theta_jumpcs_sum += ht*theta_jumpcs_aux
+ theta_Omega_sum += ht*theta_omega_aux
else : # n=/=0
# load a posteriori error estimates
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
- [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,jump_cs0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
- [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_Omega] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,jump_cs1,theta_Omega] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
# initialize objects
theta_res_aux = 0
@@ -194,7 +205,8 @@ for fineness in range(steps) : # varibale in space
theta_time_aux = 0
theta_time_val_aux = 0
theta_conv_aux = 0
- theta_conv_aux = 0
+ theta_jumpcs_aux = 0
+ theta_omega_aux = 0
# sum over all elements
for i in range(len(theta_res1)) :
@@ -203,6 +215,8 @@ for fineness in range(steps) : # varibale in space
theta_time_aux += theta_time0[i]+theta_time1[i]
theta_time_val_aux += theta_time_val0[i]+theta_time_val1[i]
theta_conv_aux += theta_conv0[i] + theta_conv1[i]
+ theta_jumpcs_aux += jump_cs0[i] + jump_cs1[i]
+ theta_omega_aux += theta_Omega0[i]**2
# sum over time steps, corresponds to L^2 norm in time
theta_res_sum += ht*theta_res_aux
@@ -210,7 +224,8 @@ for fineness in range(steps) : # varibale in space
theta_time_sum += ht*theta_time_aux
theta_time_val_sum += ht*theta_time_val_aux
theta_conv_sum += ht*theta_conv_aux
- theta_Omega_sum += ht*theta_Omega0**2
+ theta_jumpcs_sum += ht*theta_jumpcs_aux
+ theta_Omega_sum += ht*theta_omega_aux
# save results in order to compute EOC later
theta_early_max.append(theta_early_sum)
@@ -219,6 +234,7 @@ for fineness in range(steps) : # varibale in space
theta_time_max.append(theta_time_sum)
theta_time_val_max.append(theta_time_val_sum)
theta_conv_max.append(theta_conv_sum)
+ theta_jumpcs_max.append(theta_jumpcs_sum)
theta_Omega_max.append(theta_Omega_sum)
@@ -231,7 +247,8 @@ print('theta_diff_max ',np.round(theta_diff_max,40))
print('theta_time_max ',np.round(theta_time_max,7))
print('theta_time_val_max ',np.round(theta_time_val_max,5))
print('theta_conv_max ',np.round(theta_conv_max,5))
-print('theta_Omega_max ',np.round(theta_Omega_max,5))
+print('theta_jumpcs_max ',np.round(theta_jumpcs_max,5))
+print('theta_Omega_max ',np.round(theta_Omega_max,20))
print('mesh sizes ',hh)
# initialize objects
@@ -243,6 +260,7 @@ eoc_diff = []
eoc_time = []
eoc_time_val = []
eoc_conv = []
+eoc_jumpc = []
eoc_Omega = []
# go through all mesh sizes we used
@@ -257,6 +275,7 @@ for i in range(steps-1) :
aux_time = (np.log((theta_time_max[i])/(theta_time_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
aux_time_val = (np.log((theta_time_val_max[i])/(theta_time_val_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
aux_delta = (np.log((theta_conv_max[i])/(theta_conv_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_jump = (np.log((theta_jumpcs_max[i])/(theta_jumpcs_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
aux_Omega = (np.log((theta_Omega_max[i])/(theta_Omega_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
# save the EOCs
@@ -268,6 +287,7 @@ for i in range(steps-1) :
eoc_time.append(aux_time)
eoc_time_val.append(aux_time_val)
eoc_conv.append(aux_delta)
+ eoc_jumpc.append(aux_jump)
eoc_Omega.append(aux_Omega)
@@ -283,6 +303,7 @@ print('eoc_diff ',np.round(eoc_diff,5))
print('eoc_time ',np.round(eoc_time,5))
print('eoc_time_val ',np.round(eoc_time_val,5))
print('eoc_conv ',np.round(eoc_conv,5))
+print('eoc_jumpc ',np.round(eoc_jumpc,5))
print('eoc_Omega ',np.round(eoc_Omega,5))
diff --git a/2d/scaling_time_errorestimates.py b/2d/scaling_time_errorestimates.py
index c5718be..2c247b0 100644
--- a/2d/scaling_time_errorestimates.py
+++ b/2d/scaling_time_errorestimates.py
@@ -19,7 +19,7 @@ method = 'wasserstein' # using log-mean convective flux
boundary = 'Neumann'
interpol = 'morley'
-test = 'manufactured1' #'diff', 'blow-up', 'manufactured' or 'manufactured1'
+test = 'manufactured1' # 'diff' or 'blow-up' or 'manufactured'
# initialize different test cases
@@ -89,10 +89,10 @@ theta_diff_max = []
theta_time_max = []
theta_time_val_max = []
theta_conv_max = []
-theta_delta_max = []
+theta_jumpc_max = []
theta_Omega_max = []
-fineness = 4 # fix a number of refinements to the zerost mesh, i.e. fix a spatial mesh size
+fineness = 3 # fix a number of refinements to the zerost mesh, i.e. fix a spatial mesh size
steps = 8 # number of different time step sizes used
for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4, 6 ,8 , 16, 32, 64]
@@ -108,7 +108,7 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_time_val_sum = 0
theta_conv_sum = 0
- theta_delta_sum = 0
+ theta_jumpc_sum = 0
theta_Omega_sum = 0
# load 'exact' error
@@ -125,7 +125,7 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
[theta_early,theta_res0,theta_diff0,theta_time0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
- [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_Omega1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_jumpc1,theta_Omega1] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
# initialize objects
theta_early_aux = 0
@@ -133,7 +133,8 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_diff_aux = 0
theta_time_aux = 0
theta_conv_aux = 0
- theta_delta_aux = 0
+ theta_jumpc_aux = 0
+ theta_omega_aux = 0
# sum over all elements
for i in range(len(theta_early)) :
@@ -142,20 +143,23 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_diff_aux += theta_diff0[i]+theta_diff1[i]
theta_time_aux += theta_time0[i]
theta_conv_aux += theta_conv1[i]
+ theta_jumpc_aux += theta_jumpc1[i]
+ theta_omega_aux += theta_Omega0[i]**2
+
# sum over time steps, corresponds to L^2 norm in time
theta_early_sum += ht*theta_early_aux
theta_res_sum += ht*theta_res_aux
theta_diff_sum += ht*theta_diff_aux
theta_time_sum += ht*theta_time_aux
- theta_delta_sum += ht*theta_delta_aux
- theta_Omega_sum += ht*theta_Omega0**2
+ theta_jumpc_sum += ht*theta_jumpc_aux
+ theta_Omega_sum += ht*theta_omega_aux
elif n == Nt-1 :
# load a posteriori error estimates
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
- [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_jumpc0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
# initialize objects
theta_res_aux = 0
@@ -163,7 +167,8 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_time_aux = 0
theta_time_val_aux = 0
theta_conv_aux = 0
- theta_delta_aux = 0
+ theta_jumpc_aux = 0
+ theta_omega_aux = 0
# sum over all elements
for i in range(len(theta_res1)) :
@@ -172,6 +177,8 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_time_aux += theta_time0[i]
theta_time_val_aux += theta_time_val0[i]
theta_conv_aux += theta_conv0[i]
+ theta_jumpc_aux += theta_jumpc0[i]
+ theta_omega_aux += theta_Omega0[i]**2
# sum over time steps, corresponds to L^2 norm in time
theta_res_sum += ht*theta_res_aux
@@ -179,15 +186,15 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_time_sum += ht*theta_time_aux
theta_time_val_sum += ht*theta_time_val_aux
theta_conv_sum += ht*theta_conv_aux
- theta_Omega_sum += ht*theta_Omega0**2
+ theta_Omega_sum += ht*theta_omega_aux
else : # n=/=0
# load a posteriori error estimates
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n)+'.p'
- [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res0,theta_diff0,theta_time0,theta_time_val0,theta_conv0,theta_jumpc0,theta_Omega0] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
pickle_name = test+'_fineness'+str(fineness)+'_Nt'+str(Nt)+'_'+interpol+'_a posti error at time step'+str(n+1)+'.p'
- [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_Omega] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
+ [theta_res1,theta_diff1,theta_time1,theta_time_val1,theta_conv1,theta_jumpc1,theta_Omega] = pickle.load(open('pickle_files/'+pickle_name,'rb'))
# initialize objects
theta_res_aux = 0
@@ -195,7 +202,8 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_time_aux = 0
theta_time_val_aux = 0
theta_conv_aux = 0
- theta_delta_aux = 0
+ theta_jumpc_aux = 0
+ theta_omega_aux = 0
# sum over all elements
for i in range(len(theta_res1)) :
@@ -204,6 +212,8 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_time_aux += theta_time0[i]+theta_time1[i]
theta_time_val_aux += theta_time_val0[i]+theta_time_val1[i]
theta_conv_aux += theta_conv0[i]+theta_conv1[i]
+ theta_jumpc_aux += theta_jumpc0[i]+theta_jumpc1[i]
+ theta_omega_aux += theta_Omega0[i]**2
# sum over time steps, corresponds to L^2 norm in time
theta_res_sum += ht*theta_res_aux
@@ -211,7 +221,7 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_time_sum += ht*theta_time_aux
theta_time_val_sum += ht*theta_time_val_aux
theta_conv_sum += ht*theta_conv_aux
- theta_Omega_sum += ht*theta_Omega0**2
+ theta_Omega_sum += ht*theta_omega_aux
# save results in order to compute EOC later
theta_early_max.append(theta_early_sum)
@@ -220,6 +230,7 @@ for Nt in [2,3,4,6,8,16,32,64]: # go through numbers of time steps, [1, 2, 3, 4,
theta_time_max.append(theta_time_sum)
theta_time_val_max.append(theta_time_val_sum)
theta_conv_max.append(theta_conv_sum)
+ theta_jumpc_max.append(theta_jumpc_sum)
theta_Omega_max.append(theta_Omega_sum)
# print results so far
@@ -231,7 +242,8 @@ print('theta_diff_max ',np.round(theta_diff_max,40))
print('theta_time_max ',np.round(theta_time_max,6))
print('theta_time_val_max ',np.round(theta_time_val_max,5))
print('theta_conv_max ',np.round(theta_conv_max,5))
-print('theta_Omega_max ',np.round(theta_Omega_max,5))
+print('theta_jumpc_max ',np.round(theta_jumpc_max,5))
+print('theta_Omega_max ',np.round(theta_Omega_max,15))
print('mesh sizes ',hh)
# initialize objects
@@ -243,6 +255,7 @@ eoc_diff = []
eoc_time = []
eoc_time_val = []
eoc_conv = []
+eoc_jumpc = []
eoc_Omega = []
# go through all numbers of number of step sizes we used
@@ -257,6 +270,7 @@ for i in range(steps-1) :
aux_time = (np.log((theta_time_max[i])/(theta_time_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
aux_time_val = (np.log((theta_time_val_max[i])/(theta_time_val_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
aux_conv = (np.log((theta_conv_max[i])/(theta_conv_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
+ aux_jumpc = (np.log((theta_jumpc_max[i])/(theta_jumpc_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
aux_Omega = (np.log((theta_Omega_max[i])/(theta_Omega_max[i+1])))/(np.log((hh[i])/(hh[i+1])))
# save the EOCs
@@ -268,6 +282,7 @@ for i in range(steps-1) :
eoc_time.append(aux_time)
eoc_time_val.append(aux_time_val)
eoc_conv.append(aux_conv)
+ eoc_jumpc.append(aux_jumpc)
eoc_Omega.append(aux_Omega)
@@ -283,7 +298,8 @@ print('eoc_diff ',np.round(eoc_diff,5))
print('eoc_time ',np.round(eoc_time,5))
print('eoc_time_val ',np.round(eoc_time_val,5))
print('eoc_conv ',np.round(eoc_conv,5))
-print('eoc_Omega ',np.round(eoc_Omega,5))
+print('eoc_jumpc ',np.round(eoc_jumpc,5))
+print('eoc_Omega ',np.round(eoc_Omega,20))
# plot the EOCs
--
GitLab