From 141028548b480ea61d8da1952b227b061fb03fcd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 1 Nov 2023 10:40:38 +0100 Subject: [PATCH] docs: add theory part from HP report --- docs_sphinx/images/adi_scheme.svg | 1413 +++++++++++++++++ docs_sphinx/images/conservation_law.svg | 367 +++++ docs_sphinx/images/discretization_1D.svg | 563 +++++++ docs_sphinx/images/explicit_scheme.svg | 607 +++++++ .../images/grid_with_boundaries_example.png | Bin 0 -> 21250 bytes docs_sphinx/images/implicit_scheme.svg | 607 +++++++ docs_sphinx/images/numeric_overview.svg | 475 ++++++ docs_sphinx/theory.rst | 699 +++++++- 8 files changed, 4723 insertions(+), 8 deletions(-) create mode 100644 docs_sphinx/images/adi_scheme.svg create mode 100644 docs_sphinx/images/conservation_law.svg create mode 100644 docs_sphinx/images/discretization_1D.svg create mode 100644 docs_sphinx/images/explicit_scheme.svg create mode 100644 docs_sphinx/images/grid_with_boundaries_example.png create mode 100644 docs_sphinx/images/implicit_scheme.svg create mode 100644 docs_sphinx/images/numeric_overview.svg diff --git a/docs_sphinx/images/adi_scheme.svg b/docs_sphinx/images/adi_scheme.svg new file mode 100644 index 0000000..df5e2d8 --- /dev/null +++ b/docs_sphinx/images/adi_scheme.svg @@ -0,0 +1,1413 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs_sphinx/images/conservation_law.svg b/docs_sphinx/images/conservation_law.svg new file mode 100644 index 0000000..a8dec19 --- /dev/null +++ b/docs_sphinx/images/conservation_law.svg @@ -0,0 +1,367 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs_sphinx/images/discretization_1D.svg b/docs_sphinx/images/discretization_1D.svg new file mode 100644 index 0000000..8cccaa4 --- /dev/null +++ b/docs_sphinx/images/discretization_1D.svg @@ -0,0 +1,563 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs_sphinx/images/explicit_scheme.svg b/docs_sphinx/images/explicit_scheme.svg new file mode 100644 index 0000000..06cf984 --- /dev/null +++ b/docs_sphinx/images/explicit_scheme.svg @@ -0,0 +1,607 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs_sphinx/images/grid_with_boundaries_example.png b/docs_sphinx/images/grid_with_boundaries_example.png new file mode 100644 index 0000000000000000000000000000000000000000..60e4092000eeaee27d2389fcb6b9ff6ffa3281a6 GIT binary patch literal 21250 zcmd_S2{e`c+cvzL6|qwxvl5jMWr$D|il}7BSmr5n=Al9eNv1L+DJmiJOl2lhGLte( zl6mGkw*URS>v`9=-sgLsXRU92?|MJ$zPsJKz4!IIe#3bl=W!h8<#$CvX7?`oT?7JQ z_eI(B*9e5o1^Bm$Y%^Z@o0FQ4|7@|9ym+0AjI8UY;tv9Wm2mO=x$BOdzk57fLk5f1 zCugL4@~V5h*(MR{dL~vHd1>0a)2B|IaxD5DFB2pvDA?N4a^c_}Hm?>_fqWj#YwXl- z-@dJ`so{Sn92HQO*UYxn^IT3>%(G`Ec5J!SD=f4+ttnaze&-`2BfBkfx}SuG>Zhp% zq>yvtOLDqa`}}#s!@{M-akLV`cGNLl_TF88$Vd|u!erlqLN+v`)P zZKX=4-#dezP4;*I}hhx_Vc|G<>KCK(t zv}HS;79~%8|72Cgi*udQ|LG|((EIl0s&w0ZE*V{48Fl88iQ;8km)2os{#XE*SKK1zIr*{B6bFewaK7B z(#3zb^DS!3wQ5)n3MsApzK&VkAD~?voQ44dTmFwY=WMo&txG z(S)b~o*f%AmihC|&Ye4PS4on(zM`6}K^~p^ir>FKj%kpv`csxD<}z2c-OTik9s7fY zsd8E!#?v}1LKZhj>E zxwB7#H1jNLe|I=Vv_(fVVoWUu`O0X=vnSn_X87`HG&OOt(dI;NVs&+Oo0db&KdZHO z`myLF$ePs|!#BPUkq!3$+wuu*8$|Tn{WR#5Y+9yS^>iz1>rqWSQ$58q#Z#+nv?*J=DoPproY4JN;odgW$}c9+%(k zX|!uC&59B|2tX?dFNfD~VPda`FB_mN%?Rfaw>axxZ-=~`= zBqU1G)31@;5{%+?{fHgJu*-X4kwH{DRaKBj;;$&R-}@UG6}o8yLqpjMT(J!~x4W`s z+4)`O%;SR+La-S%EQ=Qt%r)HypPXg=+ts@^_h|=9q1AIaJ-z<9@va@a=|f{;ZqGEG z^|ke-a+&`V7Avc+9@%X(_q*e@T!j1D;_y^ZwIlJL@@$9en<|OCzsAN4q`rSo zPDn_2@q+cxAt{~+4lgTV#{K*EM|70d$B23p4;(lUusoqR;9u|DykBZ-AaJ?L^`9BlGLZ%Cn=*Mtx;|{Mb24sI#0iJ z_D%Fk@rMbiDk~rMS{POIFrA6BQ}0~A+j*gbwC)qVIt)C+(~FusZ7Es>L-jE(O=4B| zZsM@>MBq+G&ydcRS!xjn!_0_pf9UCv^8LKMaG+*g-&?&Z+5qvEeN zFfcHQWeq6D%w+aTH5s4c=Ds~x6V{p9PBRd=SEGB2Npr$wi=5?f4U1dyHouK&4aArg zZ`|nU%rZL^T!rgrD%OnJiNq*a`*98SO{{Nu#2{2ZA~kv)g9zyw58v&SGdBF(g-1-x z4GZ46B~xhS+OWTTyIoZ|6}42MBTnzw*jS9haDR}5`yx*+R%s_&@!D*vQ11~@L(9BD z8lTczr2BU-(Gfhu+>mL~jMab3iPe2s;yxgxmu^)Dh=m-l*bU+l`h z9m13`-|JD*)N;oVH$L8$8rwV?I?cq>@6TIIlyseKQSA&zK-tlIn1jQnQZC-=PiA5u zgHY?YZ?6&)xkN1oBGu~}$V+-l3~FiEwmQW7av?*?Jo=tL+p00hlJxE!H8u6YVBgxo zukKq9(R?wSe*M5Wi&Y|iX#J{1`-gMdtp*h(_kV4^doq^15xn) z{rh}+6<9X$sj5Q}^fZs`k@78Zh-_9XLqkI&Bk6|H~5eYy;i3{pe?7k$*{lX7nj*zEr^jclKKV) z2ejzxOO^c|4afN~^FNcnapML;R94q`OG`_ztp3zMwRD@e$73a)SeHL}fP7T+r;Z#t z6#wpBjk#d+^6leD-Run!D=RDbYooQ5#U|s-v@||3F)f}xzlC$UBaRh;3=Y4(y=4-$ zuTm(a)p2by7gUT9X&K4?OJ&3F-7p*_?zZ&ptx}P8YTH;{O^xmzUA=$J9DwB4uh$}K z%RgVX_D`Z_dLl@FgxU{50&A{dpL0;*%v`%(ld;{Gr+mOWlD<<@Hf&`;BY`R8$W6vL zw^4$(ooAY42JQSVM+rWw8zP@r-c_pA=4#3>68-MoyUG5_iQbY}@-*u)QqUK7bxhSR zS{wFQs~l1POGV0WiXKbNA8wUCckJ`)>(h-lJ13_ZHq!naxBB_BTKJM zQNiyd6p{1Yd90SVB-B=oXempBy|z-QBgvKT&Q8rV|N3cu`hSsVNvDV>`pFc<$UMAk1gIEUKByvlkIUU%h&jnHd^X zeKuOuI;ni|nyM-SUpTLh1#qXZw75N;)B}pQSee>d(xd_~jQ#p0X!-ri*KBKp@`w9v zhrXP4ntEm59P{$!OI$WH^HjWC0SQJYVk%tDJO^6c*u|%9hu&)!Mc!fVr*Xqw+}ij# zMejBkrka0eEL*RV68A0y)b~g_zCmnW^A!DVpVDpHwi#wPP4<;bx7|tZM=r9+DR}$# zV*TrxW;fDkC2XFhUS~M{*_-@8F+sB?Qot#n6b|SGJYfkaC=j4D8 zDYkBwR*g$EJII%-nS#P0KIx-VWUOpXw+%VVmBPtISM@H#jLA*qClPLT?DP3zn9+22 zBp`{JKsXy%hM|0urT7nieG!B-A1ED|vYBvKjg3z;40YF&Cr=6s3z@~79gDQ%Yc~;C zRb(kt%O~bbmj4tWH5z7|Ce61Ef#AEJx>UKkt}ce(;34rKH+S8LI$mU>E^T?2l(bhe zfR69;h~z)7h&a!@As!@E2LJZTlbnTu>3Y#Xo~?wtbO-q|ZwNE;96R=_Bf}_yN9#}F zbp1vkGs)=4+_8_A)&U6GsNwnU-Mc$0uy3&%DL`W*Q=166)a(zk5q?8BF4Ld5`4NF8 z$D-$#O`c@jl}nc%&fWd>>sL)p&CJZqS=Ps~vAjYw`}ZfSXC0#t9s^Y-nYVa4R*$B&Yb^zb4o#_Vt&7+y>$(`YV848-t zSu!Ej`dR8xk&yvBj06JHV>G(h1bC{9R zN0JRi>AISZjtA=CF3S?FmDOKTm06@~XW zvxjskJd&!&iD6-3xQKm$8`2|<{{?O9Tl+!hwY6D(FYGWE zOiWCSl26!1PJ?W5^MsnZy6KEHfneo5qG!B6sxJKbbBm0Q1jX25vo@Z*_xC%@>+9>k zBda6fc;m{A`uqX{0=K(!ucmCqelx6(HW7^CrjfgF!Q1}o7Q$UhA8|}-k&A;rH!^J9 zx^YgI#Nt2*P0)>E*Tt#9+Q1Yh!ri0`EdNh@ilJ)urCYb+H1j8aF-f?0kHl{xyr`3< zNFJPCOmXL|vvk`@IJ@Js@=j9d!nOcrqSpd%^MH*nv*zJ`t^gj8H;1g1@zR@MIs)M? zy&qn{uZ@P-yDjgOkz%B2HxrB{*tchXdU9CqPu}ozRPwa@_F+j#CpnXabT;3_Dz&5 z`!#KCZB5O>Tcr<#zB@{7l-T2zr9c0^Q%=HuqaoWGDPe#q!#0y{qa-x$z_;-nBK^dd zOD8rTWuWmdAK5(H)Yy2!to3+@PC9BFBr5)AnmiG6cMA*8Jfjf2-I(_&WgGf8CU%wWZc!&!a z^IfC}mVH%R{C#+uN!&Fzj`VYmEcnS2-I-6>K7cwq*>J>oOLRp~VcMN%Mpl>S%_^xW zeLO$?@Y=aIvvOcypz7jYVC=fuTHAsiiyTn_fjg*1NvF9{Nxryc#eJQew-8k*0|SGb zcx<;Lw|pas!ya{ZDk!CFg&7fcPvsMtwt9X-xE`#H z5b9m{OwQ{oXGhy%p6%u3h42aVa6n6xaMmL&h6^=%);3~PP}o@LLw|oU=7!c}4bjS{ zm4I``r#cIb5Q5*v#eFi(;^5%Wo;k`B!A0$t`Si?Y!U&^o6x~4jR!QH2Klxlok5-R( z2F9hQrz7V9nzoPI$D1`(V#j(BI}>RFrTzW=$L*J2r|@`}{zP$XCH!8?XEVW$BkeK; z0+wL$Vx1+u+bh27H*TCpAnNQ?M`32iosX<+U!WDWH_f=|-Qa-kqFVoG0%M^P{t3x;kAASZs$8yH;Lz9g2X?zw~ z59OIDpTLi))>NJKszX6rH{Xks^3itSvw2cwVBVc;T_1gV9(RV5Ccg5w^eJ=yBPLPg z7WdULo%PvNkJX9wHOs+%CQ2D$#^=J;y2#(I8sD9Unb^*2Y4yYhi62w7d7qkEf?N{o zH1uqCtgTjnfkXZ&kz%_KCDT%;St_???ui7MeWENLE0ao!V#@03ojI1hL-QQl?wWPx zJ#A8PB#~iy^T#5tX=o%d>4B&rhFw#1T^c(<)OfYjk0< zKcZJ^Iwc~?BL65SXX53`H)0N-?r)25YRA2CC_e9QZCGtj)02zf5d)^kHdtRXpuXnKRFWd5hzPr{!?DaE#e8@SeH~+U2U@V6buXMwD?84OHp-m-;@+z zCnkUI4>q62a*ju{^gYdc`T6-eD}A&e ze->6JJbGQ1Lyszb2mk!}^Jja}ibuu+Rl&>i6YF)B?$g*tT3u=y8c{JZkR|eNm2QtG z-%8k0c7m_!p3TdSX!lawd?1tP++3HXrM0y)AUo>ASR5#OklLBH7g17+JO2HFuZkKH zS`xoiO9+)2ESbbU@Q{%@0j4uC4a2_YB&AyprcBcw1dn<+wN)(f0mEbB0la zy?6zsDYiPuWP9YX>A=*&shU?8L(B!cz(Ee|dP=Yg^lj;R0t!{0qR3?)I81PkRmWvK zYieo=FjHH&cj8aM#6L?Hm~4gbLYzh^q65SN1eH0ifyA0_P}^2KClFD42W)Z={y>H( z`8(>d($B);JQ4?n16zIKb4d~YLpl#Et0OWM%z83**dl9d^ChtYX00fhyH;xI!ThjG z4dbNyL-3dc`Ds7@rzgwzw|#-nZiZG}Rlw8f;Rd0ly97IG%*f|$TejT0TdwQZlV92v zV~= zw<{!0m-+2Yj#lWXFGh96Tf88dU4XjnwH$3c*R=%%1^Phg)CF5iudz72(AuJ_`0ZPM zP+$PI8MOAO=Ob|Hmk!!6opEd&v1z;pL^SYK07ryd#4T0r6g~IRqc<{*ud))`3!O3q zv%mCD4h#>swX~2<=F67{rDgUdi9bk4;OH=3r+k2rhD<*;P#qHC{IM$dP@?>^sgT%} zS9qiS#Bb%AFJH_}O?PeGtXk*;76@nG-uLJuY*ZYu_f5)~HFZM0zwJZ@Gj9=$Q8kP5 z@+xYms%l(p7@MWUCUMuHU=%)Kzpd;TAEUUdnqa2~u@uBz$y`<{1f1z>akqx+W557j z6jg0e0i0FM^TDQXZCF0b1ri{UHjrGaz<%{tQtE*3-lP`ipKfiLc@cz06;LN21pXK>pu0d%(vX6m(YDP`l|*y;&q&94BxpM5Ab5s9>;UM66%D@*uT%>5|}f zD&{jdCAs-;hrml^Tl5_3LefTX*C{Bn(?~7Dl2{=UCMw2Sr1cnP>}C;n#Yrpmp_=H* znZPdt4|{++^tA30hz*~A{zU*qrp7mc-2>KPo9JLKJ$LnLaM!pXsR%nt?RP4cI<-XY z_!NH2ftvSUlFF){T z<;y?gAIxggF8G#_DI=)!$pk|yS!zkmcbNFDYe_mjD~r>CVqb>(oVIb+0X125-OkZk zupB~;we~N6U!Qh=Z)eQhrnLJ?u3yczCE}3(>CBQ_j%WOfYNBCvHOd0)7-ni=Wx632 zM3YjYLirzG-6wRJVRRKiEGz|Ic_*%!iT}aPtK!0eLf_l_L|FmYutHo`EJ{()=QwV znzOMWmWZF+9KDsD66@0vD8__XZd=uI=d4tEtoXEzA-l*}=8q;MbvRwTSPk~DEm>n8 zVOXcj5ufv6b~XmtKkoFDNx{%?=)@ul26`mwL1iQucK2 zUuTng@Fhx!7fAhbj#U07z+m2zm2~rNWZrRm_S;{RHF8jFEUd19YZKJTdi(ypIau~_ z73pckrKPb{kD+Hq-`(bW?%&_cq|-MnvA!4uT!UHnSnld=6?Df~$L#}{nzm4#X@2WV07v@O)ap(9Vk!T-I{rWoDz%G zZ?497Cf|(_budkJfy$+`ZuYud_TRh!OKp74aR|f3t8-nL=}-?>oY1auI?;pUo>gDI zbOUi9gke~Q@&20$lS5IKH5K3&Xm~W+SF**&&#ssKm~|dhKL{Sy!riihTtHCpyKRFQ zkl-`H+v>uMt=T@3aUf$w|K#@nGNQ|t#;KH(lXLjf{VXzaDRq4*b^Pzc$Q{g%Kcae9 z44s@7@#PjJt8-6_RluF7q-ZXZwry=Vw?nU+gA22vj2(*q-qk-21@_~nt1e+p0ZoHj zMxqClcZO7xq~_lIaq0DruunJtA|33N)D&JQ$wO3IX_i=Rj$L0dVk-lww=^?4&AUFu zE9r}jB|iI231Q$J@uiqYyV&24Tez2Ac2%}$vNlRPWUq1k`g}D!_edUxL)3GAk2TlU zd)tXhpaUd*p*C6llngnlGImJw;hGzT&82+f>__KM&}n9wD*O8&!S3I;4_N8aBiS~U z#!w-Oc=wjE?{0$=!KJGuQlXAqJ;6C8_e}ZJxUz zb5b5z0Eba9iX{3-lBKpWQ2NMSx`gBluxVp!J3rCeTd!&iVoQU?_3u4%b`B0xhy5ms zdD~wCHG`?Gl9rUb|2q%`iPKMSTK&6V70CTl*o59jG8DD2Oe&fm&RvB`D<%HQBkFO;`KbvG&Al@Oi zkcd}AERB<=^cr8Bll1-l`*&?aGCU%k`*wNjdz=NP=CY(F#Y7ZeAZwsJvXxnSENkQy z<>uz1{uL5xGZ!>lT}_VjV3GP4H5~Q);1-3bC20a7*FH2Xth(YUpPu^SKL7G`_f=;$ zuhO4%2>)H7GF89*cRPAVG`nBQ8ff_tBm|!~Ou};gutyPoL85@^%%EZ5kvH zh&m_`q@K}PAnvqsJ{TG}gY&S0kf1%ul(t)e;-^*an>WTl3RB1Yv%|sMJXV0Py0Ty$ zRVQdOP}OBoFurS!wURtvGX=<)i>LH%8F82QfYzq?CSz^w@5p({DoGTmvTr{LdYh6m zbZ{GiAW4?k?1YGh`BlDtU9*FE;IcyrU`f)^_@cI_`y@e5VXbj6pk54u`Yoz?l~nDXwX979q9sV&-+m>!$bB>6ia2kP!ltbX?1eG<|FqL$G*Z8dN2leFsiok*op5=xaVo;z|&pg;m zAdsj43k*tdmVa?}X*p%6aEC4B&d9_5{)~z~&n4I!mtNr(Mq} zL>%GbDu9OssTw-okv0cYvHms^RWCB-4ABR+k|z9#lw_lM)80P2#WYYl9*k}$>YSmv zsBwFjD0?{)!4`8#SGNxZ2tWp_B{Q3l`tfiS?%-MMG4Z4k-|u+>ZfxW__1?XzDVk0{ z=Q1~bxW03%9{zWn+5cI;W`537Fa;GUU6O@GN^CZnQGF>sLYm%xuY_>Kllz82*f1W) ztb!S#*VNEmWtU>)sp^r3R?lyq5H_l3>`=IL$$4e6GBAaV8`XKWbP%%` zYMNeHk<=qk8Y!Hh+b1c^ND0P7Ch;$`6!sr2qRG7E16>( z5)W^OE(wVuv#b%ZR{VwVbxqCVd!-N`o3chc|Fen(?&&|4H2!ye1NQ?zHu}H^rMKd# zh}fX49DbiB=4okZLHlt^Nc4_uffONfVl*orbr<3aN|l_haisYQI1>)I-9<#)kw44# zR6YSK_$ncx(_HYj9t4i29b{kwUudv+)(+a_X>NpaSXycWO-xLrSy%~##(Qi1bcl-P z;2o1%OlLNMKV}+7ZqdNmrUO)U?hYy?7iX$i~oNXw3(bJW` zn2eKqmh**??wgd9@IY|?bZ(k>MbVjKI|t;p#v7y0=pbBq5ml1HvZ9cfy+Ht9yY?hD zc4B@x$`ZO>Z9@Yu+(|uhWA-^?C?aah9*2g01L4rFOhNdV^)@}-2COST7s-gwgobK{onu1Gn z;|v8hP_q)eB;zjl(~P&LNlAM_ZjNo>w4A4-^qHQTlILiKVu3OTg*zjq%6}3(Jj010*U>zke!83jHbmX;g1 zwckkh;Hc_Y;I)7!>3CNNry`A{;~7$*1n-OGFex@M^m*8YE(ezE??V8{H^gi&yo`7L{`i#adI`=1$SGb# z2>x(gsVBNSIXQu@g2fkp7TDOh<`af{QZyR+C$Foj`Ws4dY$AL#sIvTml}Zbnx9rTm zt$`zg0|u>)ofuY;fl3E{S~Fo{Ai+ho%bpdV}NNw6|tg>`jx*{Pda zZ$6#?Tr)oRz3^h@pnEHqYD%cmUI=0keT5lO@8LcgB-UMF&m$fb6)m_dkg24rmSK4L z>eZ`9_8mKR%%-5nrbDLVKsxw+L>etEG4~b6B?gR+qz#*Ait`zRbBKCLkb9Y#ye*_jADDhGIf|+#5yBC93kz=Wx8mE z;rWvFRj&5Z$81tw#9-M0KhEd5ef*uHJi)z-^$vhJ85tR1=u?w>Pgr#4qT+;Lr0|U2 z0Pd%AeGomBfbDD4SKEu-$tVCNtir05?>{k|6${2&@4*cMi?m}Hxz(B z*=4P!s_crYs)K+ipP@ZsurZ5`>FYW}7t%Q-HaA*v%X$4Tp4;{RLg@A*D&@qS6Pu$$d&uru>VQq?M$gJ8_om>a}Zq z^UF|k!0Nh*>qObt&X98A{?s-QVGhlgg*wLPps^&i9P_HGuje^(glxMfCevP+QPlpo zar6D4`{JY_9mwd2hJ1GGvlrwAvouyVwp)Nk;G!`jJE)k=%`}NKaJ@HHlgw@+b|at1 zQG;XK4DFu?Cw;+x$v5&}@tFMY`(8u-{JV1}tozPm``PUmZZL{+mK8q->6u8H~`G+EkvgwB#Hw*SKQ?i>s4hV}ZRuK34fWh=%~4b1X4)@6MOFuV&g? zSfrHxJOPVLB4|PQa!D2$WF&ymVFY(@bC1W%%E}-p(1FKQ|8Vo3dHeX>fs=+LbgG`F zdja|)$r;X|6(NCl_`QlM^p!!Q>Y*R*p{en~Xm2&Zj$AZ&Q$zgIUA zeyFPgA0S$TcQ;2D1apu~m-hARc9l#H1A~JtZEa&~FP%9aK#5QmbYCcCsaPo5j1Pec zBj2^_0VUJB-CG(#QCe0%Ic)Q3Q0=%CkgKrk!UuD!i@@GJVCyKE&$N%tn)=hj69Kn@ zS358X9P7zq{4_dcOYbXLi(pnw4@MV2!8|$44p%n1OV*Iuk?c`kGdXm#d6Bv>p2@?w z4yAx(%H2b@6@pOptKs>AS6$3H9;;;~t^b8h`G>dP9@w@4Vzs8jLbKiyaR$My0UrQWPT4B}|8sJ3 z(rgBYrl` zvgj=qLuz&V^W_Dw8@hCOl2GnwhJhkY7;s3bGsNnMLK|N$d?W4`_x9~5kn~1rFiT1f zmv#_-&<{*H4z!<3Z38KZoX{zrn$)1YbNv3U>0_G+zDIAr^r|ylUTpCAZjqW~tGfZ- z`DVY$@m77%}Ccsa$N^sFpY~Os3JO`E-C)$gL#rAqu z>mNkAAc}(C{PD?i1L|Y_a>E)TLAIhMcAoov+VtBIr~aLYQ|Tl@$9X2#_^I(ZpRy}D zI=rebU+U|_&8}{O!;Z{ba!ND@ZxQOYk}`Xq4tomxs!q#(a6RCcydu;g;KJjsRGDUV zj@u8q+sRU>TY^ITJjgOSzbs_lennZi<;)KUYD2{74fC@J+% zI=&<^V%jiE6&|(cGU6U-to-0IZ^2eJ+orSb$+Gg}S`fM|z~$7SZKbM;C~4l6oeZmr z1&e@bXxDf)?Ee`c9{EXUd~w_fvVrEBs7jjnpC0aF^Kgkmeva34;J{ zFX_U)Tn7VbW1K{mSqS0F#7VesWnmOwyqUG-W9NsQ_HJ>aTgU7 zvLPPn`f-Sgy1d7xoIFs~D3aY~mye;rT@OaXP4KisZvgw#`lf`58TdJQwQ3{X?Q z5Ri;(gcijbnOxA$p&CV(o{<4&xPcE?)HO8mHXu*}0@}X}ptE zV@*#_gXTNBMFGIuKE$$P47V~fwYEHg^BYoq9o6aXEl)-6yfg0_p937saX0atlcm6J ze3z0k=UfD5Owsz<^4`6B!CmxN1f4(VPs5FtEqI)O^Z7?sw-R31suwxWVh53Uf7fp> zT1>^DTSC2or@A|m;fhJ;AwDp>SL3C7tSg!RWqGjCQ&4zB)I|aSzHx&e0Y0dnlnGRL zZFl_*H)nLXZcAdwtP}<#cEImE%xnC%L-M!#!N?D~HiDdkgDOYr08YijEE!!j99?gQ z!NI-pR4fES=EF@HT$>1IeWU=U|0jCyKWB9Q-xq}bMzcRp zHNrXqJ=@4`f9@+_|G&VcAbm*i6d?EO2=&aTf0-~hW|)LXQQl_RtEhzcSdakYFzVsy zjmAz&A2u6SLx3W#FrX5 z79ny=3!N@#&pBHfaBtfCwR-26wmxdrkMX1+T7 z01*JgjLXj#c}DjH%v{O^7E;dq|II%5Z|-r}h`=qvdiO6Lb%}F3`h$>s@z3l;?>b6v zfhD+C?^4P;kH1MzpL%BW9NP;Pzb;`|l>bH8p3M<|0D=bz1%d-yOj2D>!@~{#p;KEc zZ)T}O5yAqi9fwU*l5E4p9Uy)dVR*}GC7?kBd|XASaq*p(;I4Rh$I_fqQ6|6+Y;)QZsR>PNrx_`^dCGjB zyHSv8(s_7M$k(V_>q~y=uuMd~> zrn=@^`27S3i?|cOc>0<~4ygE%KR~eQe8^ua;;*4=g7GEUNBbX%Mh?TrlmLtZ4!EGO zaAnG+k(=*>Jv6^@sqddJlJOZgzGyx3y?^pcq=4La_s?UqJnEUUUE;cBWO*r`D_ymC zJ88h_*4QQ4^P(j5a$X@67$RA-D-BXBkPb#pD%6QD-@Us;w#V;02kb8GdO=KZqr)aR z9zmCB!Yr}yi7LQp9Te1d;Ka+?Dx}U|aAI{zB*AyfClqQHf>EBRdg!7Ny3y!GLNq?5 zX?CO622*A(QoyK|#yHmwy!7YXHq7Ght{g?d0d6i0DY-3Z+<_f{EX}Yk63i{gx7Lyr zaX-?TSpDHiYs+gaN}+<)BL(djH|?nYg@(Wr`&tFuCu{cGP#S$Fp}u%zg* zR_TcNThO~|0o@wDEf@h{{0A5X{LJrdU?AMBA3)`$I&If+GLE~s+j)uo0Ww9$V^2D0 z!qCj46^y1SJuyC~sR7iMY+j63fLa)5w)qXnsxbBs*t>mKz&-LcoS6>0G?K!&1;wnqz=F&P77JN$r!#cLC>H3b_>>)f>GG;W3$dAvBN`V=k0K> zdlQEsGe@A zk|S5T-8|dB{)F+E@Mrr=6!*@aP00c%PDKCQgd;SPZwR+<~R zR&m`6e@m&z*<8hYQU25!+M1XgC5N;?jO`p{aK^Rpva>ay5KW!sIK>bJPPN;51kQe@M=J@|lKDMwv-eWBvI?M8LGPO&< zax6Qgnue#w?40O04ccFX{TVN7(uwGGMsH$J(B!jfp^H*r5`1m^>9YwJZM`6O^80hL zZ$6wO?|;a3>zZc2oq1#K(COHzn#4czl{sPuIH*?Hi{}|9+H|NIdNSM|O^DCC{#kZ< zK7k|hTe;d+XEq>+!B zXB-Zr6XqAuN-kKo>}k5KNDQTq)SdHB&F$>2TbwwXR!zrR(DVH9;};@~WhG9KzV2Y4 ztgTHXtUC0D)Vy^17^wD}j~-&>pT?7b_J8^1y9XMfoWQAdQ55+AVjC;*3KJ9I3*p%I zAKB>lLp!*{+Co1{IeL0}{<$_Zh`)Lj8Xi73GxO%r9)clZzo!5m1n}U&1H2hEcW-V< zVxkcolsF74k7L{-6X4=Gu2YOo3`J)XcnG1hHs8@ACE5g$uX(IKkr~~ z|6Iaj4WJNLlp@_59UaNzn5O5$CxUU<#J_n{l}q9KNBj8A^i%eWuEnHB7fd&nQb`FtMYV7GN2-J9wkUFeCX!hV-2e zwq!^3PJ`kBGlMD&8WYH_JtWY0qhY)wW4K~D2vJ;l2YEg z=^)(Z&-6ZQcj|hMhGxezXh?_Ak7U1A-T~xFM@L6Q_g8&uYb(@v^zydsIY)M_Vi)*s z{0z=5EaY}q)zo~wQ-tO_F0Q2bcq#4fn}xfw-q5{2=f3K6=6OYhbZT3zgU;{A`{A3B zw6&dvqwTtnr$jJ+tL}RP+|@Lmq#(-pC_7uQUtPApX8W59G)_>0eV+8JrsZX2 z87~P;PIVlbW&z3Yh=zBE?h{`HY}-vqS+g5-v`>JGjvjITlG?|{#>a)8h0y5HA!Q?K zPe!R?8!)%4y|tpW~SUVnop zz$8V0%l_4{k4pA{8Q|SiXT{{?BwG1A+5B%@r=YoMcZ=k8Z%O-shll@R;Qa-K6J%`2o_NgDFDom#>lA7A6Z52lCg~f- zrK(`NI|O(O(5AcFbz1Jz%@2TJ#l_mQl@E9Pv~ombBTKXT`FU}ASaH=*s|_Ur_RRSF zvcTnnMT;~vFKpYqNyRna?Nz{_E1u?JqX#?MJH1cW-_YsyEG#&lqTerh22R&zr;^@@ zhA&^Rb=6QvTRwpnP*0LWFGY&~Ztm|$J} z&64=~z$D4Yy}Q|;WE220O&F<1G&GJJvJ!R_=^a00@cu9qo3VSa#ji#HLE_V7Sg*66Z9-p%L`j^u7u{~)xH$S*CGuHOqy3%3$J zZbY77D|5I<4BK8;9VGM_?$-$k;tM}Kl}d1LH+b@#lQI-sTndQy%O?!aTUw^BWM8^` zxd`3qiB6$~W0p-XF9@kn(DG=KP|7w5eWT5Y22kOPyHBVD6z|V>sc8!HH7P180vGN> z)DYE@Lt} zYr?o-osWPJlgh02cPzDJ7n}R~s)Jro3aPV?RucZ;ixO6RQ3Eel(41899Hdsw@Cv#Y zuHUhnyo=dIv*#)TgrqMvIn3dDkups^O5P9I5>`ESU%zp;K$AzW0V)Vrs`Q$w(#TWm zkOv$@GBnyzf#_nSJUJ_N2*=S{m=ShBECmdkOaFWUWi{6ndTDqT6BA65latXCZkQ1) zXR1}`*whYy54SIvV8}K0i+MpFWF@?(Xb7mg}J?aIHe!qNRtVZqKhNU8k^?64WJ(wJ2EN? z9m*FAdK7q$M!+&=Cj!q~$wD7FoZq{=kzC-*Mgq9*BbbPb+t@l-qEQiD<72a5dzKZi z@|DFQl(WHu-kokB2mXW=9VVM$ZOWVvwawEc5rn0skPykNaS^RHJV51ii*?(e7F)=8clg)XWS$8a!T&r*4DwLXD^Y&2aaADa`{+ zOz6VM>>97Hs(LZ9T|Stf2Ah=xb?*ukqS*tbRb=z%m)cq`9v(F9sW3HlOK!k+>p>Kq x|0@Xpzd-i?@4Vbycvz14DibOS-#hDDaxcE-?PvdH2z*YsctPQOisa4v{}(I5B|`uJ literal 0 HcmV?d00001 diff --git a/docs_sphinx/images/implicit_scheme.svg b/docs_sphinx/images/implicit_scheme.svg new file mode 100644 index 0000000..c8ff96c --- /dev/null +++ b/docs_sphinx/images/implicit_scheme.svg @@ -0,0 +1,607 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs_sphinx/images/numeric_overview.svg b/docs_sphinx/images/numeric_overview.svg new file mode 100644 index 0000000..da3d6ab --- /dev/null +++ b/docs_sphinx/images/numeric_overview.svg @@ -0,0 +1,475 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs_sphinx/theory.rst b/docs_sphinx/theory.rst index 4d3468f..fbc2891 100644 --- a/docs_sphinx/theory.rst +++ b/docs_sphinx/theory.rst @@ -1,15 +1,698 @@ +.. Converted content of report from Hannes Signer and Philipp Ungrund + (https://www.cs.uni-potsdam.de/bs/teaching/docs/lectures/2023/tug-project.pdf) + Theoretical Foundations ======================= -===================== -The Diffusion Problem -===================== +This section describes the theoretical foundation underlying this research +project by answering questions about the background with the mathematical and +scientific equations, the goal of numerically solving these, and their use +cases. It also discusses discretization approaches to the equations that are +needed to calculate results digitally. [BAEHR19]_ -================ -Numerical Solver -================ +Diffusion Equation +------------------ -**Backward Time-Centered Space (BTCS) Method** +The diffusion equation builds the theoretical bedrock to this research project. +It is a parabolic partial differential equation that is applied in various +scientific fields. It can describe high-level movement of particles suspended in +a liquid or gaseous medium over time and space based on Brownian molecular +motion. + +The derivation of the equation simply follows the conservation of mass theorem. +In the following, the derivation always refers to the one-dimensional case +first, but an extension to further dimensions is straightforward. First, let’s +imagine a long tube filled with a liquid or gaseous medium into which we add a +marker substance within a small volume element at the location x. + +.. figure:: images/conservation_law.svg + + Visualization of the law of mass conservation + +The boundaries of the volume element are still permeable, so that mass transport +is possible. We refer to the difference of the substance entering and leaving +the volume element as the flux :math:`\Phi` and the concentration of the marker +substance at location :math:`x` at time :math:`t` as :math:`C`. According to the +conservation of mass, it is not possible for this to be created or destroyed. +Therefore, a change in concentration over time must be due to a change in flux +over the volume element. [LOGAN15]_ [SALSA22]_ Formally this leads to + +.. _`mass conservation`: + +.. math:: + + \frac{\partial C}{\partial t} + \frac{\partial \Phi}{\partial x} = 0. + +Diffusion relies on random particle collisions based on thermal energy. +The higher the concentration of a substance, the higher the probability +of a particle collision. Fick’s first law describes this relationship as +follows + +.. _`ficks law`: + +.. math:: + + \Phi = -\alpha \cdot \frac{\partial C}{\partial x}. + +The law states that the flux moves in the direction of decreasing concentration +[1]_, proportional to the gradient of the concentration in the spatial direction +and a diffusion coefficient :math:`\alpha`. In this context, the diffusion +coefficient represents a measure of the mobility of particles and is given in +the SI unit :math:`\frac{m^2}{s}`. Substituting `ficks law`_ into `mass +conservation`_ now yields the diffusion equation for the one-dimensional case, +which gives the temporal change of the concentration as a function of the +spatial concentration gradient + +.. _`one dimensional PDE`: +.. math:: + + \begin{aligned} + \frac{\partial C}{\partial t} - \frac{\partial}{\partial x}\left(\alpha \cdot \frac{\partial C}{\partial x}\right) = 0 \\ + \frac{\partial C}{\partial t} = \alpha \cdot \frac{\partial^2 C}{\partial x^2}. + \label{eq:} + \end{aligned} + +An extension to two dimensions is easily possible and gives + +.. _`two dimensional PDE`: +.. math:: + + \begin{aligned} + \frac{\partial C}{\partial t} = \alpha_x \cdot \frac{\partial^2 C}{\partial x^2} + \alpha_y \cdot \frac{\partial^2 C}{\partial y^2}. + \label{eq:pde_two_dimensional} + \end{aligned} + +Since this equation now contains both a first order temporal and a +second order spatial derivative, it belongs to the group of partial +differential equations (PDE). In particular, the diffusion problem +belongs to the group of parabolic PDEs. With these, we describe the +evolution of the problem via a first-order time derivative. +[CHAPRA15]_ Furthermore, it describes diffusion under +the assumption of a constant diffusion coefficient in each direction. +Therefore, we distinguish once into the case of homogeneous but +anisotropic diffusion, where the diffusion coefficients are constant in +one spatial direction but may vary between them, and the case of +heterogeneous diffusion. In this case, the diffusion can also change +within one spatial direction. This distinction is important, because the +subsequent discretization of the problem for numerical solution differs +with respect to these two variants. -**Forward Time-Centered Space (FTCS) Method** +Numeric and Discretization +-------------------------- + +The solution of the diffusion equation described above can be solved +analytically for simple cases. For more complicated geometries or +heterogeneous diffusion behavior, this is no longer feasible. In order +to obtain sufficiently accurate solutions, various numerical methods are +available. Basically, we can distinguish between methods of finite +differences and finite elements. The former are mathematically easier to +handle, but limited to a simpler geometry, whereas the latter are more +difficult to implement, but can simulate arbitrarily complicated shapes. +[BAEHR19]_ + +Since this work is limited to the simulation of two-dimensional grids +with heterogeneous diffusion as the most complicated form, we restrict +ourselves to the application of finite differences and focus on a +performant computation. + +There are in principle two methods for solving finite differences with +opposite approaches, as well as one method that uses both methods +equally – the so-called Crank-Nicolson (CRNI) method. + +.. figure:: images/numeric_overview.svg + + Methods of finite differences + +The idea of finite differences is to replace the partial derivatives of +the PDE to be solved by the corresponding difference quotients for a +sufficiently finely discretized problem. Taylor’s theorem forms the +basis for this possibility. In essence, Taylor states that for a smooth +function [2]_ :math:`f: \mathbb{R} \rightarrow \mathbb{R}`, it is +possible to predict the function value :math:`f(x)` at one point using +the function value and its derivative at another point +:math:`a \in \mathbb{R}`. Formally we can write this as + +.. math:: + + \begin{aligned} + f(x) =& f(a) + f'(a) (x-a) + \frac{f''(a)}{2!}(x-a)^2 + ... + \frac{f^{(n)}(a)}{n!}(x-a)^n + R_n\\ + R_n =& \int_{a}^x \frac{(x-t)^n}{n!} f^{(n+1)}(t) dt. + \end{aligned} + +:math:`R_n` denotes the residual term, which also indicates the error in +the case of premature termination of the series. An approximation of the +first derivative is now done simply by truncating the Taylor series +after the first derivative and transforming accordingly +[CHAPRA15]_. This leads to + +.. _`first order approximation`: + +.. math:: + + \begin{aligned} + f(x_{i+1}) =& f(x_i) + f'(x_i) (x_{i+1}-x_i) + R_1 \\ + f'(x_{i}) =& \underbrace{\frac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i}}_{\text{first order approximation}} - \underbrace{\frac{R_1}{x_{i+1} - x_i}}_{\text{truncation error}}. + \end{aligned} + +The `first order approximation`_ shows a so-called forward finite difference +approximation in space, since we use the points at :math:`i` and :math:`i+1` for +the approximation of :math:`f'(x_i)`. Correspondingly, a backward and centered +finite difference approximation is possible applying the appropriate values. + +For the approximation of higher derivatives, the combination of several +Taylor series expansions is possible. For example, we set up the +following two Taylor series to derive the centered approximation of the +second derivative + +.. _`second order approximations`: +.. math:: + + \begin{aligned} + f(x_{i+1}) =& f(x_i) + f'(x_i) \underbrace{(x_{i+1} - x_i)}_{h} + \frac{f^{''}(x_i)}{2} \underbrace{(x_{i+1} - x_{i})^2}_{h^2} + R_2 \\ + f(x_{i-1}) =& f(x_i) - f'(x_i) \underbrace{(x_{i-1} - x_i)}_{h} + \frac{f^{''}(x_i)}{2} \underbrace{(x_{i-1} - x_{i})^2}_{h^2} +R_2. + \end{aligned} + +Adding [3]_ both `second order approximations`_ and rearranging them according to +the second derivative yields the corresponding approximation + + +.. _`second order approximation`: +.. math:: + + \begin{aligned} + f^{''}(x_i) = \frac{f(x_{i+1}) - 2\cdot f(x_i) + f(x_{i-1})}{h^2}. + \end{aligned} + +Another possibility of the derivation for the second order approximation results +from the following consideration. The second derivation is just another +derivation of the first one. In this respect we can represent the second +derivative also as the difference of the approximation of the first derivative +[CHAPRA15]_. This results in the following representation, which after +simplification agrees with equation `second order approximation`_ + +.. _`first order derivative`: +.. math:: + + \begin{aligned} + f^{''}(x_i) = \frac{\frac{f(x_{i+1}) - f(x_i)}{h} - \frac{f(x_i) - f(x_{i-1})}{h}}{h}. + \label{eq:first_order_derivative} + \end{aligned} + +The above equations are all related to a step size :math:`h`. Let us imagine a +bar, which we want to discretize by dividing it into individual cells. + +.. figure:: images/discretization_1D.svg + + Discretization of a one-dimensional surface + + +In the discretization, the step size :math:`h` would correspond to the distance +between the cell centers :math:`\Delta x`. Now, the question arises how to deal +with the boundaries? This question is a topic of its own, whereby numerous +boundary conditions exist. At this point, only the cases relevant to this work +will be presented. First, let us look at the problem intuitively and consider +what possible boundary conditions can occur. A common case is an environment +that is isolated from the outside world. That is, the change in concentration at +the boundary over time has to be zero. We can generalize this condition to the +*Neumann boundary condition*, where the derivative at the boundary is given. The +case of a closed system then arises with a derivative of zero, so that there can +be no flow across the boundary of the system. Thus, the Neumann boundary +condition for the left side yields the following formulation + +.. math:: + + \begin{aligned} + \text{Neumann condition for the left boundary:~} \frac{\partial C(x_{0-\frac{\Delta x}{2}}, t)}{\partial t} = g(t) \text{,~} t > 0. + \end{aligned} + +The second common case is a constant boundary. Again, the case can be +generalized, this time to the so-called *Dirichlet boundary condition*. Here, +the function value of the boundary is given instead of its derivative. For the +example we can write this boundary condition as + +.. math:: + + \begin{aligned} + \text{Dirichlet condition for the left boundary:~} C(x_{0 - \frac{\Delta x}{2}}, t) = h(t) \text{, ~} t > 0. + \end{aligned} + +In practise, constants are often used instead of time-dependent +functions. [LOGAN15]_ [CHAPRA15]_ [SALSA22]_ + + +With this knowledge, we can now turn to the concrete implementation of +finite differences for the explicit and implicit scheme. Centered +differences are always used for the spatial component. Only in the time +domain we distinguish into forward (FTCS method) and backward (BTCS +method) methods, which influence the corresponding solution +possibilities and stability properties. + +Explicit Scheme (FTCS) +~~~~~~~~~~~~~~~~~~~~~~ + +FTCS stands for *Forward Time, Centered Space* and is an explicit +procedure that can be solved iteratively. Explicit methods are generally +more accurate, but are limited in their time step, so that incorrectly +chosen time steps lead to instability of the method. Hoewever, there are +estimates, such as the Courant-Friedrich-Lewy stability conditions, +which gives the corresponding maximum possible time steps. +[CHAPRA15]_ + +The goal now is to approximate both sides of the `one dimensional PDE`_. As the +name of the method FTCS indicates, we use a forward approximation for the left +temporal component. For the right-hand side, we use a central approximation +using the `first order derivative`_. This results in the following approximation +for the inner cells (in the example the cells :math:`x_1`, :math:`x_2` and +:math:`x_3`) and a constant diffusion coefficient :math:`\alpha` + +.. _`approximate first order diffusion`: +.. math:: + + \begin{aligned} + \frac{C_i^{t+1}-C_i^t}{\Delta t}=\alpha \frac{\frac{C_{i+1}^t-C_i^t}{\Delta x}-\frac{C_i^t-C_{i-1}^t}{\Delta x}}{\Delta x}. + \label{eq:approx_first_order_diffusion} + \end{aligned} + +By rearranging this equation according to the concentration value for the next +time step :math:``C_i^{t+1}`, we get + +.. _`explicit scheme`: +.. math:: + + \begin{aligned} + C_i^{t+1} = C_i^t + \frac{\alpha \cdot \Delta t}{\Delta x^2} \left[C_{i+1}^t - 2C_i^t + C_{i-1}^t \right]. + \label{eq:explicit_scheme} + \end{aligned} + +At this point, it should be clear why this method is an explicit procedure. On +the right side of the `explicit scheme`_, there are only known quantities, so +that an explicit calculation rule is given. The basic procedure is shown in +Figure 4. In the case of an explicit procedure, we use the values of the +neighboring cells of the current time step to approximate the value of a cell +for the next time step. + +.. figure:: images/explicit_scheme.svg + + Illustration for using the existing values in the time domain (green) for approximation of the next time step (red) + + +As described above, the `explicit scheme`_ applies only to the inner cells. For +the edges we now have to differentiate between the already presented boundary +conditions of Dirichlet and Neumann. To do this, we look again at the +`approximate first order diffusion`_ and consider what would have to be changed +in the case of constant boundary conditions (Dirichlet). It quickly becomes +apparent that in the case of the left cell the value :math:`C_{i-1}` corresponds +to the constant value of the left boundary :math:`l` and in the case of the +right cell :math:`C_{i+1}` to the value of the right boundary :math:`r`. In +addition, the length difference between the cells is now no longer :math:`\Delta +x`, but only :math:`\frac{\Delta x}{2}`, as we go from the center of the cell to +the edge. For a given flux, the first derivative has to take this value. In our +case, this is approximated with the help of the difference quotient, so that +this is omitted in the case of a closed boundary (the flux is equal to zero) or +accordingly assumes a constant value. For the left-hand boundary, this results +in the following modifications, whereas the values for the right-hand boundary +can be derived equivalently + +.. _`boundary conditions`: +.. math:: + + \begin{aligned} + \text{Dirichlet for the left boundary:~}& \frac{C_0^{t+1}-C_0^t}{\Delta t}=\alpha \frac{\frac{C_{1}^t-C_0^t}{\Delta x}-\frac{C_0^t-l}{{\frac{\Delta x}{2}}}}{\Delta x}\\ + &C_{0}^{t+1} = C_{0}^t + \frac{\alpha \Delta t}{\Delta x^2} \left [C_{1}^t - 3 C_0^t + 2l \right]\\ + \text{Closed Neumann for the left boundary:~}& \frac{C_0^{t+1}-C_0^t}{\Delta t}=\alpha \frac{\frac{C_{1}^t-C_0^t}{\Delta x}-\cancelto{0}{l}}{\Delta x}\\ + &C_{0}^{t+1} = C_{0}^t + \frac{\alpha \Delta t}{\Delta x^2} \left [C_1^t - C_0^t \right]. + \label{eq:boundary_conditions} + \end{aligned} + +Again, it is straightforward to extend the `explicit scheme`_ to the second +dimension, so that the following formula is simply obtained for the inner cells + +.. math:: + + \begin{aligned} + C_{i,j}^{t+1} = C_{i,j}^t +& \frac{\alpha_x \cdot \Delta t}{\Delta x^2} \left[C_{i, j+1}^t - 2C_{i,j}^t + C_{i, j-1} \right] \\ + +& \frac{\alpha_y \cdot \Delta t}{\Delta y^2} \left[C_{i+1, j}^t - 2C_{i,j}^t + C_{i-1, j} \right]. + \end{aligned} + +Here, it is important to note that the indices :math:`i` and :math:`j` +for the two-dimensional cases are now used as usual in the matrix +notation. That means :math:`i` marks the elements in x-direction and +:math:`j` in y-direction. + +The previous equations referred to homogeneous diffusion coefficients, i.e. the +diffusion properties along a spatial direction are identical. However, we are +often interested in applications where different diffusion coefficients act in +each discretized cell. This extension of the problem also leads to a slightly +modified variant of our `one dimensional PDE`_, where the diffusion coefficient +now is also a function of the spatial component + +.. math:: + + \begin{aligned} + \frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left[\alpha(x) \frac{\partial C}{\partial x} \right]. + \end{aligned} + +For this case, the literature recommends discretizing the problem +directly at the boundaries of the grid cells, i.e. halfway between the +grid points [FROLKOVIC1990]_. If we +follow this scheme and approximate the first derivative for C at the +appropriate cell boundaries, we obtain the following expressions + +.. math:: + + \begin{aligned} + \begin{cases} + \alpha(x_{i+\frac{1}{2}}) \frac{\partial C}{\partial x}(x_{i+\frac{1}{2}}) = \alpha_{i+\frac{1}{2}} \left(\frac{C_{i+1} - C_i}{\Delta x}\right)\\ + \alpha(x_{i-\frac{1}{2}}) \frac{\partial C}{\partial x}(x_{i-\frac{1}{2}}) = \alpha_{i-\frac{1}{2}} \left(\frac{C_{i} - C_{i-1}}{\Delta x}\right). + \end{cases} + \end{aligned} + +With a further derivation we now obtain a spatially centered +approximation with + +.. math:: + + \begin{aligned} + \frac{\partial }{\partial x}\left(\alpha(x) \frac{\partial C}{\partial x} \right) \simeq& \frac{1}{\Delta x}\left[\alpha_{i+\frac{1}{2}}\left(\frac{C^t_{i+1}-C^t_i}{\Delta x}\right)-\alpha_{i-\frac{1}{2}}\left(\frac{C^t_i-C^t_{i-1}}{\Delta x}\right)\right] \\ + \frac{C^{t+1}_{i}-C^t_i}{\Delta t}=& \frac{1}{\Delta x^2}\left[\alpha_{i+\frac{1}{2}} C^t_{i+1}-\left(\alpha_{i+\frac{1}{2}}+\alpha_{i-\frac{1}{2}}\right) C^t_i+\alpha_{i-\frac{1}{2}} C^t_{i-1}\right]\\ + C^{t+1}_{i} =& C^t_i + \frac{\Delta t}{\Delta x^2} \left[\alpha_{i+\frac{1}{2}} C^t_{i+1}-\left(\alpha_{i+\frac{1}{2}}+\alpha_{i-\frac{1}{2}}\right) C^t_i+\alpha_{i-\frac{1}{2}} C^t_{i-1} \right]. + \end{aligned} + +The value for the inter-cell diffusion coefficients can be determined by +either the arithmetic or harmonic mean of both cells, with the harmonic +mean being preferred for the default case. Again, we can extend this +equation to two dimensions, resulting in the following iteration rule + +.. math:: + + \begin{aligned} + C_{i, j}^{t+1}= & C_{i, j}^t+ \frac{\Delta t}{\Delta x^2}\left[\alpha_{i, j+\frac{1}{2}}^x C_{i, j+1}^t-\left(\alpha_{i, j+\frac{1}{2}}^x+\alpha_{i, j-\frac{1}{2}}^x\right) C_{i, j}^t+\alpha_{i, j-\frac{1}{2}}^x C_{i, j-1}^t \right]+ \\ & \frac{\Delta t}{\Delta y^2}\left[\alpha_{i+\frac{1}{2}, j}^y C_{i+1, j}^t-\left(\alpha_{i+\frac{1}{2}, j}^y+\alpha_{i-\frac{1}{2}, j}^y\right) C_{i, j}^t+\alpha_{i-\frac{1}{2}, j}^y C_{i-1, j}^t \right]. + \end{aligned} + +The corresponding equations for the boundary conditions can be derived +analogously to the homogeneous case (cf. `boundary conditions`_) by adjusting +the relevant difference quotients to the respective boundary condition. As +described initially, the FTCS procedure is accurate but not stable for each time +step. The Courant-Friedrichs-Lewy stability condition states that the time step +must always be less than or equal the following value + +.. math:: + + \begin{aligned} + \Delta t \leq \frac{\min (\Delta x^2, \Delta y^2)}{4 \cdot \max (\alpha^x, \alpha^y)}. + \end{aligned} + +That is, the maximum time step is quadratically related to the +discretization and linearly to the maximum diffusion coefficient. +Especially for very fine-resolution grids, this condition has a +particularly strong effect on that required computing time. This is in +contrast to the implicit methods, which we will now look at in more +detail. Unlike the explicit methods, the implicit methods are +unconditionally stable. + +Implicit Scheme (BTCS) +~~~~~~~~~~~~~~~~~~~~~~ + +The main difference to the explicit method is that the discretization is not +done at the time step :math:`t` as in the `approximate first order diffusion`_, +but now at the time step :math:`t+1`. Hence, the neighboring cells in the next +time step are used for the approximation of the middle cell. + +.. figure:: images/implicit_scheme.svg + + Illustration of the implicit method, where the values of the neighboring cells in the next time step are used for the approximation + + +The `approximate first order diffusion`_ thus changes to + +.. math:: + + \begin{aligned} + \frac{C_i^{t+1}-C_i^t}{\Delta t} & =\alpha \frac{\frac{C_{i+1}^{t+1}-C_i^{t+1}}{\Delta x}-\frac{C_i^{t+1}-C_{i-1}^{t+1}}{\Delta x}}{\Delta x} \\ + & =\alpha \frac{C_{i-1}^{t+1}-2 C_i^{t+1}+C_{i+1}^{t+1}}{\Delta x^2} + \end{aligned} + +Now there is no possibility to change the equation to :math:`C^{t+1}_i` +so that all values are given. That means :math:`C^{t+1}_i` is only +implicitly indicated. If we know put all known and unkown values on one +side each and define :math:`s_x = \frac{\alpha \Delta t}{\Delta x^2}` +for simplicity, we get the following expression + +.. math:: + + \begin{aligned} + s_x \cdot C_{i+1}^{t+1} + (-1-s_x) \cdot C_i^{t+1} + s_x \cdot C_{i-1}^{t+1} &= -C_i^t. + \label{eq:implicit_inner_grid} + \end{aligned} + +This applies only to the inner grid without boundaries. To derive the equations +for the boundaries, the reader can equivalently use the procedure from the FTCS +method with given `boundary conditions`_. Thus, for constant boundaries with the +values :math:`l` and :math:`r` for the left and right boundaries, respectively, +the following two equations are obtained + +.. math:: + + \begin{aligned} + \left(-1-3 s_x\right) \cdot C_0^{t+1}+s_x \cdot C_1^{t+1} &= -C_0^t - 2s_x \cdot l \\ + s_x \cdot C_{n-1}^{t+1}+\left(-1-3 s_x\right) \cdot C_n^{t+1} &= -C_n^t - 2s_x \cdot r. + \end{aligned} + +The question now arises how values from the next time step can be used for the +approximation. Here, the answer lies in the simultaneous solution of the +equations. Let’s look again at the example in Figure for conservation law and +establish the implicit equations for the same. We obtain the following system of +five equations and five unknowns + +.. math:: + + \begin{aligned} + \begin{cases} + \left(-1-3 s_x\right) \cdot C_0^{t+1}+s_x \cdot C_1^{t+1} &= -C_0^t - 2s_x \cdot l \\ + s_x \cdot C_{2}^{t+1} + (-1-s_x) \cdot C_1^{t+1} + s_x \cdot C_{0}^{t+1} &= -C_1^t \\ + s_x \cdot C_{3}^{t+1} + (-1-s_x) \cdot C_2^{t+1} + s_x \cdot C_{1}^{t+1} &= -C_2^t \\ + s_x \cdot C_{4}^{t+1} + (-1-s_x) \cdot C_3^{t+1} + s_x \cdot C_{2}^{t+1} &= -C_3^t \\ + s_x \cdot C_{3}^{t+1}+\left(-1-3 s_x\right) \cdot C_4^{t+1} &= -C_4^t - 2s_x \cdot r. + \end{cases} + \end{aligned} + +If we transfer this system of equations to a matrix-vector system, we +get + +.. math:: + + \begin{aligned} + \begin{bmatrix} + -1-3s_x & s_x & 0 & 0 & 0 \\ + s_x & -1-3s_x & s_x & 0 & 0 \\ + 0 & s_x & -1-3s_x & s_x & 0 \\ + 0 & 0 & s_x & -1-3s_x & s_x \\ + 0 & 0 & 0 & s_x & -1-3s_x + \end{bmatrix} + \cdot + \begin{bmatrix} + -C_0^{t+1}\\ + -C_1^{t+1}\\ + -C_2^{t+1}\\ + -C_3^{t+1}\\ + -C_4^{t+1} + \end{bmatrix} + = + \begin{bmatrix} + -C_0^t - 2s_x \cdot l \\ + -C_1^t \\ + -C_2^t \\ + -C_3^t \\ + -C_4^t - 2s_x \cdot r + \end{bmatrix}. + \end{aligned} + +This system can now be solved very efficiently, since it is a so-called +tridiagonal system. Very fast solution algorithms exist for this, such +as the Thomas algorithm. + +In principle, this procedure can be extended again to the +two-dimensional case, but here no tridiagonal system arises any more, so +that the efficient solution algorithms are no longer applicable. +Therefore, one uses a trick and solves the equations in two half time +steps. In the first time step from :math:`t\rightarrow t+\frac{1}{2}` +one space direction is calculated implicitly and the other one +explicitly, whereas in the second time step from +:math:`t+\frac{1}{2} \rightarrow t+1` the whole process is reversed and +the other space direction is solved by the implicit method. This +approach is called the alternative-direction implicit method, which we +will now examine in more detail. + +First, we consider the `two dimensional PDE`_ again, which represents the PDE +for diffusion in the two-dimensional case with homogeneous and anisotropic +diffusion coefficients. As explained in the upper sections, we can approximate +the second derivative for each spatial direction as follows + +.. _`Taylor series`: +.. math:: + + \begin{aligned} + \frac{\partial^2 C^t_{i, j}}{\partial x^2} =& \frac{C^t_{i, j-1} - 2C^t_{i,j} + C^t_{i, j+1}}{\Delta x^2} \label{eq:taylor_2_x}\\ + \frac{\partial^2 C^t_{i, j}}{\partial y^2} =& \frac{C^t_{i-1, j} - 2C^t_{i,j} + C^t_{i+1, j}}{\Delta y^2}. \label{eq:taylor_2_y} + \end{aligned} + +Explicit solution methods, as can be used in the one-dimensional case, +require smaller time steps to remain stable when applied to solve +two-dimensional problems. This leads to a significant increase in the +computational effort. Implicit techniques also require more +computational capacity, since the matrices can no longer be represented +tridiagonally in the two- or three-dimensional case and thus cannot be +solved efficiently. A solution to this problem is provided by the ADI +(alternating-direction implicit) scheme, which transforms +two-dimensional problems back into tridiagonal matrices. Here, each time +step is performed in two half-steps, each time solving implicitly in one +axis direction. [CHAPRA15]_ For the above equations, +the ADI scheme is defined as follows + +.. math:: + + \begin{aligned} + \begin{cases} + \frac{C_{i,j}^{t+\frac{1}{2}}-C_{i,j}^t}{\frac{\Delta t}{2}} = \alpha_x \cdot \frac{\partial^2 C_{i,j}^{t+\frac{1}{2}}}{\partial x^2} + \alpha_y \frac{\partial^2 C_{i,j}^t}{\partial y^2} \\ + \frac{C_{i,j}^{t+1}-C_{i,j}^{t+\frac{1}{2}}}{\frac{\Delta t}{2}} = \alpha_x \cdot \frac{\partial^2 C_{i,j}^{t+\frac{1}{2}}}{\partial x^2} + \alpha_y \frac{\partial^2 C_{i,j}^{t+1}}{\partial y^2} + \end{cases}. + \label{eq:adi_scheme} + \end{aligned} + +Now the `Taylor series`_ can be substituted into the equation and the term +:math:`\frac{\Delta t}{2}` can be placed on the right side of the equation. The +following expression is generated when introducing :math:`s_x = \frac{\alpha_x +\cdot \Delta t}{2 \cdot \Delta x^2}` and :math:`s_y = \frac{\alpha_y \cdot +\Delta t}{2 \cdot \Delta y^2}` as new variables + +.. math:: + + \begin{aligned} + C_{i,j}^{t+\frac{1}{2}} - C_{i,j}^t &= s_x \cdot \left[~C_{i, j-1}^{t+\frac{1}{2}}-2 C_{i,j}^{t+\frac{1}{2}}+C_{i, j+1}^{t+\frac{1}{2}}\right] + s_y~ \cdot \left[~ C_{i-1, j}^t - 2\cdot C_{i,j}^t + C_{i+1,j}^t\right] \\ + C_{i,j}^{t+1} - C_{i,j}^{t+\frac{1}{2}} &= s_x \cdot \left[C_{i, j-1}^{t+\frac{1}{2}}-2 C_{i,j}^{t+\frac{1}{2}}+C_{i, j+1}^{t+\frac{1}{2}}\right] + s_y \cdot ~\left[ C_{i-1, j}^{t+1} - 2\cdot C_{i,j}^{t+1} + C_{i+1,j}^{t+1}\right]. + \end{aligned} + +As it can be obtained by the following Figure the above equations are only valid +for the inner grid. + +.. figure:: images/grid_with_boundaries_example.png + + Example grid with boundary conditions + +At the edge of +the grid, starting from the cell center, it is not possible to take a +full step :math:`dx`, but only a half step :math:`\frac{dx}{2}`. +Therefore, the approximation for the inner cells (shown as an example +for the x-direction) + +.. math:: + + \begin{aligned} + \frac{C_{i,j}^{t+1} - C_{i,j}^t}{\Delta t} = \frac{\frac{C_{i, j+1}^{t+1} - C_{i,j}^{t+1}}{\Delta x} - \frac{C_{i, j}^{t+1} - C_{i, j-1}^{t+1}}{\Delta x}}{\Delta x} = \frac{C_{i, j+1}^{t+1} - 2 \cdot C_{i, j}^{t+1} + C_{i, j-1}^{t+1}}{\Delta x^2} + \end{aligned} + +is no longer valid and we have to perform the same changes to the +boundary conditions as for the other methods by replacing the +corresponding concentration values or difference quotiens with the +respective boundary concentrations or fluxes depending on the boundary +condition (Neumann or Dirichlet). + +We can now quickly introduce the ADI scheme for heterogeneous diffusion +by using the same discretization logic as for the FTCS procedure and +discretizing the points between cells. + +This results in the following ADI scheme + +.. math:: + + \begin{aligned} + \begin{cases} + \frac{C_{i, j}^{t+\frac{1}{2}}-C_{i, j}^t}{\Delta \frac{t}{2}}= & \frac{1}{\Delta x^2}\left[\alpha_{i, j+ \frac{1}{2}} C_{i, j+1}^{t+ \frac{1}{2}}-\left(\alpha_{i, j+ \frac{1}{2}}+\alpha_{i, j- \frac{1}{2}}\right) C_{i, j}^{t+ \frac{1}{2}}+\alpha_{i, j-\frac{1}{2}} C_{i, j-1}^{t+\frac{1}{2}}\right]+ \\ & \frac{1}{\Delta y^2}\left[\alpha_{i+\frac{1}{2}, j} C_{i+1, j}^t-\left(\alpha_{i+\frac{1}{2}, j}+\alpha_{i-\frac{1}{2}, j}\right) C_{i, j}^t+\alpha_{i-\frac{1}{2}, j} C_{i-1, j}^t\right] \\ + \frac{C_{i, j}^{t+1}-C_{i, j}^{t+\frac{1}{2}}}{\Delta \frac{t}{2}}= & \frac{1}{\Delta x^2}\left[\alpha_{i, j+\frac{1}{2}} C_{i, j+1}^{t+\frac{1}{2}}-\left(\alpha_{i, j+\frac{1}{2}}+\alpha_{i, j-\frac{1}{2}}\right) C_{i, j}^{t+\frac{1}{2}}+\alpha_{i, j-\frac{1}{2}} C_{i, j-1}^{t+\frac{1}{2}}\right]+ \\ & \frac{1}{\Delta y^2}\left[\alpha_{i+\frac{1}{2}, j} C_{i+1, j}^{t+1}-\left(\alpha_{i+\frac{1}{2}, j}+\alpha_{i-\frac{1}{2}, j}\right) C_{i, j}^{t+1}+\alpha_{i-\frac{1}{2}, j} C_{i-1, j}^{t+1}\right]. + \end{cases} + \end{aligned} + + +.. math:: + + \begin{aligned} + \begin{cases} + -s_x \alpha_{i, j+\frac{1}{2}}^x C_{i, j+1}^{t+\frac{1}{2}}+\left(1+s_x\left(\alpha_{i, j+\frac{1}{2}}^x+\alpha_{i, j-\frac{1}{2}}^x\right)\right) C_{i, j}^{t+\frac{1}{2}}-s_x \alpha_{i, j-\frac{1}{2}}^x C_{i, j-1}^{t+\frac{1}{2}}= \\ s_y \alpha_{i+\frac{1}{2}, j}^y C_{i+1, j}^t+\left(1-s_y\left(\alpha_{i+\frac{1}{2}, j}^y+\alpha_{i-\frac{1}{2}, j}^y\right)\right) C_{i, j}^t+s_y \alpha_{i-\frac{1}{2}, j}^y C_{i-1, j}^t\\ + -s_y \alpha_{i+\frac{1}{2}, j}^y C_{i+1, j}^{t+1}+\left(1+s_y\left(\alpha_{i+\frac{1}{2}, j}^y+\alpha_{i-\frac{1}{2}, j}^y\right)\right) C_{i, j}^{t+1}-s_y \alpha_{i-\frac{1}{2}, j}^y C_{i-1, j}^{t+1}= \\ s_x \alpha_{i, j+\frac{1}{2}}^x C_{i, j+1}^{t+\frac{1}{2}}+\left(1-s_x\left(\alpha_{i, j+\frac{1}{2}}^x+\alpha_{i, j-\frac{1}{2}}^x\right)\right) C_{i, j}^{t+\frac{1}{2}}+s_x \alpha_{i, j-\frac{1}{2}}^x C_{i, j-1}^{t+\frac{1}{2}}. + \end{cases} + \end{aligned} + + +Rearranging the terms according to known and unknown quantities gives us again +the scheme as in the one-dimensional case and allows us to set up a tridiagonal +system for each row of the grid. We also have to take into account that the ADI +method always requires two calculation steps for each complete time step. +However, since the implicit method is unconditionally stable, it shows its +adavantage especially with larger time steps. The basic procedure of the ADI +scheme with the two time steps is shown here: + +.. figure:: images/adi_scheme.svg + + Illustration of the iterative procedure in the ADI process + +Crank-Nicolson method +--------------------- + +Another implicit method is the Crank-Nicolson (CRNI) method, which uses both +forward and backward approximations. As this method does not represent a main +goal of this work, it is only briefly mentioned here. Essentially, the +Crank-Nicolson method averages the two finite differences variants at the +corresponding point. For the one-dimensional case this results in + +.. math:: + + \begin{aligned} + \frac{C^{t+1}_i - C_i^{t}}{\Delta t} = \frac{1}{2} \cdot \left[\frac{C^{t}_{i+1} - 2C^t_i + C_{i+1}^t}{\Delta x^2} + \frac{C^{t+1}_{i+1} - 2C_i^{t+1} + C_i^{t+1}}{\Delta x^2} \right]. + \end{aligned} + +By transforming the terms accordingly, it is also possible to generate a +tridiagonal matrix, as in the BTCS method, which can be solved very +efficiently with the Thomas algorithm. [CHAPRA15]_ +Like the BTCS method, CRNI is unconditionally stable, but becomes +inaccurate if the time steps are too large. Compared to the FTCS and +BTCS methods, which have a temporal truncation error of +:math:`\mathcal{O}(\Delta t)`, the CRNI method with an error of +:math:`\mathcal{O}(\Delta t^2)` has an advantage and should be preferred +for time-accurate solutions. + +.. [LOGAN15] Logan, J. David. Applied Partial Differential Equations. + Undergraduate Texts in Mathematics. Cham: Springer International + Publishing, 2015. https://doi.org/10.1007/978-3-319-12493-3. + +.. [SALSA22] Salsa, Sandro, and Gianmaria Verzini. Partial Differential + Equations in Action. Vol. 147. UNITEXT. Cham: Springer + International Publishing, 2022. + https://doi.org/10.1007/978-3-031-21853-8. + +.. [CHAPRA15] Chapra, Steven C., and Raymond P. Canale. Numerical Methods for + Engineers, 2015. + + +.. [BAEHR19] Baehr, Hans Dieter, and Karl Stephan. Wärme- Und Stoffübertragung. + Berlin, Heidelberg: Springer Berlin Heidelberg, 2019. + https://doi.org/10.1007/978-3-662-58441-5. + +.. [FROLKOVIC1990] Frolkovič, Peter. “Numerical Recipes: The Art of Scientific + Computing.” Acta Applicandae Mathematicae 19, no. 3 (June + 1990): 297–99. https://doi.org/10.1007/BF01321860. + +.. [1] + The minus sign ensures that the direction of the flux follows the + decrease in concentration. + +.. [2] + A function is called smooth if it is arbitrarily often + differentiable. + +.. [3] + The remainder term is not considered further at this point, but is + :math:`\mathcal{O}(h^2)` in the centered variant. Thus, the centered + variant is more accurate than the forward or backward oriented + method, whose errors can be estimated with :math:`\mathcal{O}(h)`. + [CHAPRA15]_