From 11dab4e0f3b9f98b71ac9559320cecbe32ff3f84 Mon Sep 17 00:00:00 2001 From: ratnania Date: Tue, 27 Feb 2024 15:21:46 +0100 Subject: [PATCH] update --- _toc.yml | 1 + chapter1/boundary-conditions.md | 64 ++++++++++ .../boundary-conditions/line-boundaries.png | Bin 0 -> 4086 bytes .../boundary-conditions/square-boundaries.png | Bin 0 -> 9725 bytes chapter2/poisson-nitsche.ipynb | 93 ++++++++------ chapter2/poisson.ipynb | 119 +++++++----------- chapter4/subdomains.md | 1 + 7 files changed, 169 insertions(+), 109 deletions(-) create mode 100644 chapter1/boundary-conditions.md create mode 100644 chapter1/images/boundary-conditions/line-boundaries.png create mode 100644 chapter1/images/boundary-conditions/square-boundaries.png diff --git a/_toc.yml b/_toc.yml index b98b174..0804ebd 100644 --- a/_toc.yml +++ b/_toc.yml @@ -33,6 +33,7 @@ parts: - file: chapter1/space - file: chapter1/rules - file: chapter1/poisson + - file: chapter1/boundary-conditions - caption: Linear Problems chapters: - file: chapter2/poisson diff --git a/chapter1/boundary-conditions.md b/chapter1/boundary-conditions.md new file mode 100644 index 0000000..6764b3d --- /dev/null +++ b/chapter1/boundary-conditions.md @@ -0,0 +1,64 @@ +# Boundary Conditions +*Author: Ahmed Ratnani* + +SymPDE & Psydac allows you to use both strong and weak boundary conditions. +We start first by explaining how to identify a boundary in **NCube** domains, such as **Line**, **Square** and a **Cube**. + +## Identifying a boundary in **NCube** domains + +Since we're dealing with **NCube** domains, the best way to identify a boundary is to use the couple **(axis, extremity)**. The axis is usually defined as perpendicular to the boundary, while the **extremity** gives the lower or upper bound of the interval associated to the **axis**. + +Let's see it through the following examples; + +For a **Line**, the following figure gives the value to access a boundary. + +![png](images/boundary-conditions/line-boundaries.png) + +For a **Square**, we have, + +![png](images/boundary-conditions/square-boundaries.png) + +For a **Cube**, it is similar and would not be helpful to visualize it here. + +## Essential boundary conditions + +Essential boundary conditions are usually treated in two ways, either in the strong or weak form. The latter can be achieved through the use of Nitsche's method. You can check the associated examples in the sequel. + +For the essential boundary conditions, usually, one needs one of the following expressions; + +$$ +\begin{align} +u &= 0, \quad \text{(homogeneous case)}, \\ +u &= f, \quad \text{(inhomogeneous case)}, \\ +\mathbf{v} &= 0, \quad \text{(homogeneous case)}, \\ +\mathbf{v} &= \mathbf{g}, \quad \text{(homogeneous case)}, +\end{align} +$$ + +where $u$ and $f$ are scalar functions while $\mathbf{v}$ and $\mathbf{g}$ denote vector functions. + + +Let `Gamma` denotes the boundary $\Gamma$. +The normal derivative is an available operator in SymPDE, it is provided as `Dn` operator. Alternatively, you can also define it manually, + +```python + nn = NormalVector('nn') + dn = lambda a: dot(grad(a), nn) +``` + + +we have the following equivalences, + +| Mathematical Expression | SymPDE Expression | Example | +| --------------------------------------------------------- | ------------------------------ | ---------------------------------- | +| $u = 0$ on $\Gamma$ | `EssentialBC(u, 0, Gamma)` | | +| $u = f$ on $\Gamma$ | `EssentialBC(u, f, Gamma)` | | +| $\partial_n u = 0$ on $\Gamma$ | `EssentialBC(dn(u), 0, Gamma)` | | +| $\partial_n u = f$ on $\Gamma$ | `EssentialBC(dn(u), f, Gamma)` | | +| $\mathbf{v} = 0$ on $\Gamma$ | `EssentialBC(v, 0, Gamma)` | | +| $\mathbf{v} = \mathbf{g}$ on $\Gamma$ | `EssentialBC(v, g, Gamma)` | | +| $\mathbf{v} \cdot \mathbf{n} = 0$ on $\Gamma$ | ? | | +| $\mathbf{v} \times \mathbf{n} = 0$ on $\Gamma$ | ? | | +| $\mathbf{v} \cdot \mathbf{n} = f,\mathbf{g}$ on $\Gamma$ | ? | | +| $\mathbf{v} \times \mathbf{n} = f,\mathbf{g}$ on $\Gamma$ | ? | | + diff --git a/chapter1/images/boundary-conditions/line-boundaries.png b/chapter1/images/boundary-conditions/line-boundaries.png new file mode 100644 index 0000000000000000000000000000000000000000..b5e38ba77d8339fb53ccf6d2b1f0708f6d4c99e0 GIT binary patch literal 4086 zcmd^CdsGu=7LP=tytGQISwT>d^|3tiAQABu5J6Ef$Wv>bXUbEvewn2BB<20D$)W9b}K=uAghH60$R}hCNrZ!ckAxDf9yGX&YYZYe)oR& ze)rzrz278T7#i$gJIE7`U898xjH5F2fn%SM}Tbb~KvpKJ|hawOkjE$Y7c~ zcRsOl=dcn`mD-&f=+5ClN~MBC=1=tELiK^k0&rXVsIs{aa`zf ze<3$yWyP{3lcs-#5-!%*#DwRluUMn~!>LEwX|<3RWuDaFxf zM50b0bbWc`@o;=0j7SqG&!GH7d~%3VCIM9fBtf;reluiQ6(WyRHZc?h_AH8KDW;U< zy7{mSPd1fD(hnfqzrXioHn5ki@h;{73vnVQLoyAST>?<4*HB!A6%v2EJ|H0qQ~|UC z64Y`Tc*FY-Lp2JCxuF4J15|1ZO?sJrp1J>U^RVQN+PylyIig?B;j2sr!>AlqV{3p> z?W-zNMNY6NFoj{T3{oS>gSi8dWRt!$H@NrSETbiYPdp)tfYAgd4Z#yGF2_S9yNAms z-ziTjRD+3Oau=@n_UgaICGTx6)et6y)&GG@?tE}KFes?x{CaUnZ8OK;rujxIaQ^}S z_%hH_0TLwHox_*nOGt;{^#Um)MBsE5QVO+zU>)A&_?RR%h%m|&g}iP8uZRWb)n*Hp zM2BRE6u%9`pzZLH9RMx}H3G(uKPiwS5((Zi0onva4TmcsF-|!x7`6Z#+$LnqzyM5* zUjkA>@7@-qMGnP(ijY1eAuZa5;;3cYa9`uCWrzGLq!~ZiGcvMyY z2gzV2%s;tDsX_GWi^X!ln5EP=5+NWr;2);M)yj{mDyPw`?H2|5N9bZ78G=9cSU}5c z3TLl+-)&Q~ky-J{D966dZIv5z!cng6EfYn41zpo8i2}nSuGx&79wgjsZ1b-izDXpw z_F#gjqJOQ>;o_fu`l-98C$zWU+V#Yl+S=yBJ2R5X&LocU4$0OtKh-l6E>+wOy7@)k zC4oJ7JQBIWjElW~or{-toOxWfYqg`!In?Ubj?r}Xdvv$0-f~9!+`6ZWJGD7waa(Vt z*&eu1+etgi*j^tK9Sd}7=HW^jx+bo@Wn6=O`NmokyTLYgoxC~t^oZx{xF%_gYkHPb z^YP^edi=Z#_oVlFfB*e4k5f~&Xe;sulE%*5%Veg~4cE7q=@ZePN1q;>yRK_-!R^sr z8#Rj)=3e;2#J0=*(}g4Nx2m9dG4wl`ADJ#9JHO@Zo@JC>t-AS12~=N-ri^-zKESPM zOgrdQ?6f>9E>Cdcz#4a%uwBD*X|hh;FRgUE*AZzPw_AH*xz76ejs>Z{+hBgyJ-9JX z*Y6lL?A9(k5Le<9DIOYavdzhny1$Ai2#-dk9+^_{HC?~Ly1k_Tx=Z`VtCw}n^&LNA zWA?P$m#pEL^BB%6cK2$UD!=JOGj50FPG*Q`xE()z;<9dWgUuc9JrAt+9?B{kGw=M6 zxcuYY&Ik9|MiiaWK6TJl*V`?LU6rETx&OxpLsBbe<-^IXKPwt0UH4{)jL{4~gSW;? z<8-J4?p$BDOv!V0NZOuw&JaxRo!EJ(d~9*K>+p!wW5X`6O-D!EsBkb<9m<}zGjmSF zjjHshKezaGWc5~bot_4YIKZxXkap0qILZ*+v1PrtDaChG?9h_#CPrp)tFVODG3QCk z+WPq)>5B#m&JVT!HmmvlzNAb(w52bKb-IYYp;7)vK$GsP0%_Mig|X-GGt37VW!J3J z4^5RE@xj_hc)5oUnOM-EK5$w{J@qT&OHZ^!zz2_YV z)1#Q&?+7|xVYL6aMEuov@kYNr*)9Es58@l{3Tw|bhc*1h*^p;u;Wh{RZQPvJZ&q2= zL;LQGZ6DYZs>DW|9g(~G*-)>{oQklBP}$vr=T9n9BhO`N3b2uO9kkpdkH2o>TWQum zTIF*sG`^9^bfcG_U%Ta+cE^oEWOnrZ(ZwbCO7+6FiKc`p7Zo4q?0VK;t=zhF?u=_0 zxi1b+mKLo1LSX3N7e2XMv6*q(wvLbOFI@!)wrlS%OO@%hi8s?ybE-dkHraXPlnsu1 z9x|$@<_Qd?Iq`aC+h#rA>g3VF6`g`7Y4^iT1p}3s%lW=1ldHqq`_@j$FaOGe#VpS1 zItA*L+tOojtKV_NqdJ7=i-d@W&HXI?Z_n?!z3>Qxx0d?U1R1NP)w++*OwYL&@`KBv zR*pk;W#jJIm!FREjc-~Vk@*5T+aB-DpHo+w)^yL9nrO@s$bTmtGh1jdU9td!f|j|on=*2*$V*G<9f;cvs7g>||79MRaFlny zcjgaP=~daPqG#4A=|y!XOx-D3eQS*2T$4?42Cby;JiWemk;8{ibH?`_cdC92s%3a0K*6VBxdlL9RPg}GgH1Jr!>h*sI8BIF~ literal 0 HcmV?d00001 diff --git a/chapter1/images/boundary-conditions/square-boundaries.png b/chapter1/images/boundary-conditions/square-boundaries.png new file mode 100644 index 0000000000000000000000000000000000000000..818c1173a76460e5553b460c3f18755f3056fb29 GIT binary patch literal 9725 zcmeHN2~?Bk(oTqw)=I3k7LY`2Ek#S&f?){)vPOS5Z~=t^MhFlHNz6h52!?tEQI>kG zVg*Erz0xRKHJ}h6E=UytX$=UJ5Q{(zLV^Y(o80#!!D`#P{O3RS-20!~bAa=GGxN^O zJM+Br&Nqh)zs&^0`Ag@+U@*h=>%6zZVDM1z`_WqlfKt&^eFo6QET5gMIBq1FLV}sw zdugA{v9{4PCd=I3+Z>A}(&^R_WMUM97{|1xl30KQ_;IudawM4)p_RefVpm()Vy&>a z?RK{2_8zv-AFQ3Vt(`-VR-PC|qUscGq9u|k6rwrS7iVn?RIS`;XJ>Bj37&TBVN+Rb z?ej)5mCXe-u9HKU!)h&!5f9YbL3&pK4uKRNO9NW4wsyAQ#}1Hu$qZ5?4J`AWmC~Ay zCo!008dc}1y|taSgO)*J5u>!!aHjWqKR4YMIawD- z$Gsb)xBFqexB7W(NeGYM8evO=a`1+>Ih}D_0)t48-b9NaQGntIuFj#=_N%o%MI>k` zj?OqOA&N0)M`soyY$GSmZDrPRtZfvV96@5fu{{=zMq!cZZxX|4R4OT4msOocB7;HW zyeSq*qv&?0lc9r+zo8u>Z2L_}A!;Bwf)%Zmb;9Y~gXrr>qVkK!Jdnw&JFiNEY(c4L1(l)nU&K%&KwSd0YVx=y1##DvV~;ESD0 zFXtC9U_qIgm5F{02(O4l9RQ4wz(Gf4t6x*Y#~|Dj6xD=64&9eN}F4itbX&j@n7 zIo2TxdL(kmkZ-^jEeqItjXOtujYQ%?D2K{~cEKLToT@p>uT?->eJ$r8QIO%+(@Y^p zLDeuEzz!-|o}fUo$l&;S=ooS2hzJNWp1+3|28~UPn63Muu7RQrnS&Jk>+pb`JNR>W zaMI=eH9R={4jk-tw0{d4bS3y-f`bu=fJT`ERtCv%EK&`PX}o zS?=8DgO;IN0DZG*yUu?6nsFG+=*4<(kL}#h;hvE1+3bQ_qpUsmk7Z!Zza2PZ{LaOj z-x?&PUi_ru0P5Uf`UqpME{mc+oJgaq?=RFhwIB>YVH>G&GugUN=)OY_Dq<*0;9KgbvkjNj!-} zI*H+^y9)98OC!xK&@){|6kE>gek3+ymgi)@q>B%vgecsNrWL$y|CZDf2yyp!o7(m= zWaFJmM}ED7W=3`2)n>v|ohff@9$wiePDR6{h|cpbU?|aYp)%*RanIt~539`+8LVLS zeL1HM$wwAh7WEQxGG4x){Bh#r@|m20amAuWM)<9K$0vQO#4gwEyxM<;sY)9HOfe6| zHrN47+xeng0u964`HeEO#aNceJYD#KM15UPWYf}0%MLYm{b<5NeOGR=;(4%F3MK?5 zTG4~0^rK+PJ66JLM^S4W`eL%j?ib@bG~Q8$!t?#(EyegAOTLX^eHgE5`DBLtcw1Ew zk}$%wTfntg8u=my1CvT7$gzDGE;}uB$#jBzO=0akY1jX-#koE?VR}P)BF!L2G8($n z6+6LF-;)03-z&vgv5v5a)z=vUqAA zLEQuw;L?M0^tnSy{U^QNIE$%Hc}Ux)6|Nnl*{koAk2%3aL!*)a4d;y8gKAg7!uF9P zQlr|B@X0G{+*g)OlihJQ;wS&fH z1?K)6J~#tNpe1nG5)_KgL)g2(ME0gs4q667@&XNJvGW9`YN1@|*~>$qVWB1ue331; zJvu^S*@4{;(wE%6g(v;g1;Qx|rXXYKGUv|7l0m-<$1y`hpYOt*A1N;ZjhBL{=}At0 zRr+D3k^+iQ)|^+Q^exOxSytP_tfD9ds&cS%^= zyKH&TJN=@M1vU2+BbD2Vl!VqSDMk`P8LL`_QC|0R-1H>kvcoeUgfjAUZU|mr@u(ri zM!=3AdfuCOU$0;2rYvcTE$Z+U+O~}Qmv|e9YQFk%u}Y|?e4D$q>WR6;(Q}`(KxlHp zjCOnfK05(vMQlM@&|A~tH+Upf-DmY6;IWl%H%{mib5TZI&SMH~Ahtx(Fuy#?1{0dR z;pP^8P~KTs5WTeJN@-+q>CVQQuLpnd&tk#2gr%ehUFc$cngSs_Ra4Ti*ZOb_Jbyx5 zaw2CckwsOf%1hxl^`#o$*264OtP9RTWR2u;LRz{;-bbGP>Y7=O!ZIbsP})HltBnmi z;8jk?s84W(?V9wD}T8_Hp1 zJfJ%KGy+j!WeY%t(I-3MvVu-}O&<&iV2Fh(!tO|5IP%7DHUg|s-2n`7>enErUmI4r zK&TK~LX3Z`Y+6R|2#OiT7o}w=KqnrPh;`J88w|C0)LuEZrHI+AdcLOhHj`f9| zSrmwqa~)q`M(#QM00+9hy|(Xrk)CR^@ZO6M4p7KZB*D1n@sVqg+4R}hJU#6zkF(X_ z?&zt$g7=LNINRQWW6H8Q8#kY5lxx{6(J}0qKFR3j`I*UtBujLuB^AtfsLZ zVm3z*ev2taE6**EbJ?O}df~RnTTYP4?Asa6vv%N2(eNyDSheyI&siF3z%$||-S%&Q zV`wf2@jIefdUjPVwDC~J_PuKxyIWu?J(|m!`OathhDl(+xfcL}Zge6E8&=0dwRo>u#iNOD%{*U1NJCI{*!JRqD-^Ws?54649rf^$2ovdz< z_wu_21ZL?uKgk~6SU4^bXvj%ur?Y6M$LCzzuokrJfk9t-^yU?(%`ECnriDk#>(@8i z6m<~7GTcm)EfNF6gBw4EOQA~#?UaTvT_DYp`5B394k`bME`>fKI8r-dJUCjIKJ6N3 ztjS+96Qw?f6CrpU!@qPrqB`3y=2o5X3Q+}U#T6=srZ>TxJdT0CZ zv@4U|&%&Nn85vFrh&h=ln1bpwCxmk8M&a!q#_ms4RIBuASu$mh`T~5P<_Rxo?;Hg?YKPIB9=MAV&F(KfAOrTWTDF#&oUthGu$8 zG7~V>q|yi`d*w9#j{Q&dMJ^fhZ4!^8?#`2z^!p)2lFpZ(B_^u*@X|4lj-tNg?J0FD zCi-ELQ#VQ)^+aPDMZo^Wl6-oyUx3-4doH(Tl#w$Of|&n~XnkHse)Cz@31i(+jh zpKzK|O~f@4w1THum2q#_Ke;?UOkF2m()QTva%n3wb)4xR(%@%yxu<5HRPiU9%+Mv0 zZ!lq5ZfHq5HU0UaD7FX}IOTpLB_^*hT0=7^Xl(is?B7Ua!T6G=@ph{T?fwjk$*#~Qe7v4$A$tk8Y%h)Ek_Eqe%3c=ir%Rjc_LXAX0tsLQ(rW7~b zSCn35yJ~O!7N>H33@x2t+1HvcD>)b!gGmpQSq^h)OS9gyI zvJxiR#=XpViJm87yU;@dgF8CfJsoL>>r5pWl6#JHnX9#+rYj~jFX+<4yTzqrC9%Bo zYp8c+&x)#lEV0fvkgggrX#%&r6G|3hkrU8Tcz~zrMtBbeNOQf}nq3~FU&$xin0yX_ zeSX(;+?{$`!jbhQ>*C?6mMQD+Tu>WXGAAAfwX7u5Jw)a9;o!q0|B|wr7atr_c*cfE zgZUlaWmz#E`Oo-4si5Uf>!7137ZCObIM&NB998?ipi77&|)s~WD!9(`< z_2g!z2A>_!kj8t?x_a&y!rot6D}1`kHD<}B{sv~i^Xr{jwGvyGp#qrY?pZOL()6Og_MdrAaz0(MWQglH5chDcQ+X&d(m?t7 z8vF+u;`9M9{!@N_n(*CZxWd)|m%d2Byds8*F{#K8#gf-So+1x$)cxUH3y%zC! z=~$gbTrv`U(nQ1@jV4V!ju$N{n5esyLl~2|37=Dw7m<>CyE=-3Do4|b2Y5r>D`U>8 zdwtmlZ!F|DXi7F^$?k+YfgVN5MP+h;`z;Opt1>Pt)ai!;{3KP`J=`8ql5aR z|NEmke#H9^{Bn>kk`3>=a2SDhx(*kaEFjKP6(NuV3>b~4DcCUE|6I>KcD6g?Z& z$b&p-lj}TPbw?EZkV7t@1;G;GWj1WSce+;!i zI+fdjamlUx`N->V)RV&qagAkE*U%#v(%}|d)^ljJv^lwZmK+D9pz};=i z{1p8bSZWEDmgbYojv>$wjKPH;zWE@+Qkb@Fd^DSL01P zM9nyj?)ktKHxZ(~CQp?gmTTI=i`;$8GD7iBX7~wpi8m{LI%zj(ILV!_%BJ7_QV1N3 zt&sQW<^65J$mIc#UHJ9}SETua=Y4Rsm7pSXQ! zy8X*Af)Px)u4+N>RS*ZphB*2;L{9USQ{k>SLCctH+b*x+ztw~uo&0zkA5hN-oycEVq7{7__7 ze_`Q|JkHd(9U}`wb!oHkF(2<&Fx*;np5VSzHc?=%a{JW1chNM17t{Q(Jsd8bNPvhZ)y?D%fzhZfUZUU||nxxFQx)=1u(yY~KHT_<7V{N&~KZTz{ZJusghkV0_9ngFPnqnM00A?iOyTG%%Y26rT>MA=9Mk~J= zsa4)iWirayRY5(??a4d8zBrU$ecTx)vb^*|S|l_+Bb4;{UHEDfU{S!VpsX1#3&|5Y ztht_Ajd{^kJhJ=3L8|oqWXFjJ^)vZ!TVGN2g<&86l1_T$_Kt~moOs)^fKT#%K2Z@1 zTN?Z7Di9RFrJxNa#nAxpr)9?pN`mKSIW-7z(ob8i zwJZ}aetM_gNHS5)?xwlM5T6F+YTD(vwpf19yMt530m~86f~=;DvFoXtoA|^kTaD^gF5-2z_ZK`vP);8Yl^)XpGX~ zn@96p17xlqTb_I%x)`+MA}s0o4;AfIUQ??DIFtEj^amB0LosQJ!)#G7eGlg63NVcd zPl~vih%WwB?*`4&vVJ{Jdp^X??U1@Ym9s@d&V`luYhdKH`@~sGj%JpZg0#C)18tpq3nXXQ}!p2(tk8V d6VJfc=e(1gd#rXB__r?DdY{eSrJlRL_$%v!hzI}x literal 0 HcmV?d00001 diff --git a/chapter2/poisson-nitsche.ipynb b/chapter2/poisson-nitsche.ipynb index 052e688..68d5626 100644 --- a/chapter2/poisson-nitsche.ipynb +++ b/chapter2/poisson-nitsche.ipynb @@ -18,8 +18,8 @@ "\n", "$$\n", "\\begin{align}\n", - " - \\nabla^2 u = f \\quad \\text{in~$\\Omega$}, \\quad \\quad \n", - " u = g \\quad \\text{on~$\\partial \\Omega$}.\n", + " - \\nabla^2 u = f \\quad \\text{in $\\Omega$}, \\quad \\quad \n", + " u = g \\quad \\text{on $\\partial \\Omega$}.\n", "\\end{align}\n", "$$" ] @@ -54,12 +54,12 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 18, "id": "ecd1330a", "metadata": {}, "outputs": [], "source": [ - "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.expr import BilinearForm, LinearForm, integral, Norm\n", "from sympde.expr import find, EssentialBC\n", "from sympde.topology import ScalarFunctionSpace, Square, element_of\n", "from sympde.calculus import grad, dot, Dn\n", @@ -80,16 +80,19 @@ "u,v = [element_of(V, name=i) for i in ['u', 'v']]\n", "\n", "# bilinear form\n", - "a = BilinearForm((u,v), integral(domain, dot(grad(v), grad(u)) + kappa*u*v - u*Dn(v) - v*Dn(u)))\n", + "a = BilinearForm((u,v), integral(domain, dot(grad(v), grad(u))))\n", + "an = BilinearForm((u,v), integral(domain.boundary, kappa*u*v - u*Dn(v) - v*Dn(u)))\n", "\n", "# linear form\n", - "g = sin(pi*x)*sin(pi*y)\n", - "f = 2*pi**2*sin(pi*x)*sin(pi*y)\n", + "ue = sin(pi*x)*sin(pi*y)\n", + "g = ue\n", + "f = 2*pi**2*sin(pi*x)*sin(pi*y)\n", "\n", - "l = LinearForm(v, integral(domain, f*v + kappa*g*v - g*Dn(v)))\n", + "l = LinearForm(v, integral(domain, f*v))\n", + "ln = LinearForm(v, integral(domain.boundary, kappa*g*v - g*Dn(v)))\n", "\n", "# Variational problem\n", - "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v))" + "equation = find(u, forall=v, lhs=a(u,v) + an(u,v), rhs=l(v) + ln(v))" ] }, { @@ -102,7 +105,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 19, "id": "3d86674b", "metadata": {}, "outputs": [], @@ -113,7 +116,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 20, "id": "94c83267", "metadata": {}, "outputs": [], @@ -125,7 +128,7 @@ "Vh = discretize(V, domain_h, degree=degree)\n", "\n", "# Discretize equation\n", - "#equation_h = discretize(equation, domain_h, [Vh, Vh])" + "equation_h = discretize(equation, domain_h, [Vh, Vh])" ] }, { @@ -138,53 +141,69 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 21, "id": "8f1b2d68", "metadata": {}, "outputs": [], "source": [ - "#uh = equation_h.solve()" + "equation_h.set_solver('gmres', info=False, tol=1e-8)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 22, "id": "e96f321c", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "uh = equation_h.solve(kappa=1e3)" + ] }, { - "cell_type": "code", - "execution_count": null, - "id": "c7ae0b9e", + "cell_type": "markdown", + "id": "d87fd79f", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "## Computing the error norm" + ] }, { - "cell_type": "code", - "execution_count": null, - "id": "9ad386a6", + "cell_type": "markdown", + "id": "08ee5cac", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "### Computing the $L^2$ norm" + ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 23, "id": "068bdb95", "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "aa88a918", - "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.00021402394834116976\n" + ] + } + ], + "source": [ + "u = element_of(V, name='u')\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(u - ue, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "l2norm_h = discretize(l2norm, domain_h, Vh)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(u=uh)\n", + "\n", + "# print the result\n", + "print(l2_error)" + ] } ], "metadata": { diff --git a/chapter2/poisson.ipynb b/chapter2/poisson.ipynb index de5c678..0e26eb1 100644 --- a/chapter2/poisson.ipynb +++ b/chapter2/poisson.ipynb @@ -12,9 +12,10 @@ "\n", "$$\n", "\\begin{align}\n", - " - \\nabla^2 u = f \\quad \\text{in $\\Omega$}, \\quad \\quad \n", - " u = 0 \\quad \\text{on $\\Gamma_D$}, \\quad \\quad \n", - " \\partial_n u = g \\quad \\text{on $\\Gamma_N$}.\n", + " - \\nabla^2 u = f \\quad &\\text{in $\\Omega$}, \\\\ \n", + " u = 0 \\quad &\\text{on $\\Gamma_0$}, \\\\\n", + " u = g_i \\quad &\\text{on $\\Gamma_I$}, \\\\\n", + " \\partial_n u = g_n \\quad &\\text{on $\\Gamma_N := \\partial \\Omega \\setminus \\left( \\Gamma_0 \\cup \\Gamma_I \\right)$}.\n", "\\end{align}\n", "$$\n", "\n", @@ -32,13 +33,7 @@ "\n", "- $V \\subset H^1(\\Omega)$, \n", "- $a(u,v) := \\int_{\\Omega} \\nabla u \\cdot \\nabla v ~ d\\Omega$,\n", - "- $l(v) := \\int_{\\Omega} f v ~ d\\Omega + \\int_{\\Gamma_N} g v ~ d\\Gamma$.\n", - "\n", - "## Examples\n", - "\n", - "- [Poisson on a square domain with Homogeneous Boundary Conditions](https://github.com/pyccel/sympde/blob/devel-documentation/docs/examples/2d_poisson_dir0.py)\n", - "- [Poisson on a square domain with Homogeneous and Neumann Boundary Conditions](https://github.com/pyccel/sympde/blob/devel-documentation/docs/examples/2d_poisson_dir0_neu.py)\n", - "- [Poisson on a square domain with Homogeneous, Inhomogeneous and Neumann Boundary Conditions](https://github.com/pyccel/sympde/blob/devel-documentation/docs/examples/2d_poisson_dir_neu.py)" + "- $l(v) := \\int_{\\Omega} f v ~ d\\Omega + \\int_{\\Gamma_N} g_n v ~ d\\Gamma$." ] }, { @@ -51,19 +46,26 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 8, "id": "d742586c", "metadata": {}, "outputs": [], "source": [ - "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.expr import BilinearForm, LinearForm, integral\n", "from sympde.expr import find, EssentialBC, Norm, SemiNorm\n", "from sympde.topology import ScalarFunctionSpace, Square, element_of\n", - "from sympde.calculus import grad, dot\n", + "from sympde.calculus import grad, dot, laplace\n", + "from sympde.topology import NormalVector, Union\n", "\n", - "from sympy import pi, sin\n", + "from sympy import pi, sin\n", + "\n", + "from psydac.api.discretization import discretize\n", "\n", "domain = Square()\n", + "Gamma_0 = domain.get_boundary(axis=0, ext=-1)\n", + "Gamma_i = domain.get_boundary(axis=0, ext=1)\n", + "Gamma_n = domain.boundary.complement(Union(Gamma_0, Gamma_i))\n", + "nn = NormalVector('nn')\n", "\n", "V = ScalarFunctionSpace('V', domain)\n", "\n", @@ -74,34 +76,25 @@ "# bilinear form\n", "a = BilinearForm((u,v), integral(domain , dot(grad(v), grad(u))))\n", "\n", + "# exact solution\n", + "ue = sin(pi*x) * (1+y*sin(pi*y/3))**2\n", + "L = lambda w: - laplace(w)\n", + "f = L(ue)\n", + "gi = ue\n", + "gn = ue\n", + "\n", "# linear form\n", - "f = 2*pi**2*sin(pi*x)*sin(pi*y)\n", "l = LinearForm(v, integral(domain, f*v))\n", "\n", + "# Boundary term for the Neumann BC\n", + "ln = LinearForm(v, integral(Gamma_n, v * dot(grad(gn), nn)))\n", + "\n", "# Dirichlet boundary conditions\n", - "bc = [EssentialBC(u, 0, domain.boundary)]\n", + "bc = [EssentialBC(u, 0, Gamma_0)]\n", + "bc += [EssentialBC(u, gi, Gamma_i)]\n", "\n", "# Variational problem\n", - "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)" - ] - }, - { - "cell_type": "markdown", - "id": "62ac1fd4", - "metadata": {}, - "source": [ - "\n", - "This very simple Python code reflects well the abstract mathematical framework needed for variational formulations.\n", - "The structure of the code is as follows,\n", - "\n", - "1. Create a domain.\n", - "2. Create a space of *scalar* functions over the domain.\n", - "3. Create elements from this function space. These elements will denote the test and trial functions.\n", - "4. Create the Bilinear and Linear forms, $a$ and $l$ respectively.\n", - "5. Create Essential Boundary Conditions.\n", - "6. Create the variational problem.\n", - "\n", - "Most of the time, you will need to follow the same steps, with some minor variants depending on the problem you're considering." + "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v)+ln(v), bc=bc)" ] }, { @@ -112,27 +105,9 @@ "## Discretization" ] }, - { - "cell_type": "markdown", - "id": "51095918", - "metadata": {}, - "source": [ - "We shall need the **discretize** function from **PsyDAC**." - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "id": "a2a0a2a1", - "metadata": {}, - "outputs": [], - "source": [ - "from psydac.api.discretization import discretize" - ] - }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 23, "id": "00e54163", "metadata": {}, "outputs": [], @@ -143,7 +118,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 24, "id": "5999c62b", "metadata": {}, "outputs": [], @@ -168,7 +143,17 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 25, + "id": "004dfdb3", + "metadata": {}, + "outputs": [], + "source": [ + "equation_h.set_solver('gmres', info=False, tol=1e-8)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, "id": "541192ee", "metadata": {}, "outputs": [], @@ -202,7 +187,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 27, "id": "5925c6cd", "metadata": {}, "outputs": [ @@ -210,13 +195,11 @@ "name": "stdout", "output_type": "stream", "text": [ - "9.49756716724664e-07\n" + "0.00042499840633644024\n" ] } ], "source": [ - "ue = sin(pi*x)*sin(pi*y)\n", - "\n", "u = element_of(V, name='u')\n", "\n", "# create the formal Norm object\n", @@ -242,7 +225,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 28, "id": "e5c1a8b8", "metadata": {}, "outputs": [ @@ -250,7 +233,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "9.768702410721585e-05\n" + "0.025283770887649267\n" ] } ], @@ -286,7 +269,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "9.769164089397408e-05\n" + "0.025287342563909892\n" ] } ], @@ -303,14 +286,6 @@ "# print the result\n", "print(h1_error)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9250897b", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/chapter4/subdomains.md b/chapter4/subdomains.md index ef83d5e..51cb34d 100644 --- a/chapter4/subdomains.md +++ b/chapter4/subdomains.md @@ -1,4 +1,5 @@ # Subdomains +*Author: Ahmed Ratnani* In this section, we consider a domain $\Omega$ which a union of multiple subdomain, *i.e.*