From 40e4f4e05b0127909325c8dba0891cced7adb02e Mon Sep 17 00:00:00 2001 From: RadostW <51670923+RadostW@users.noreply.github.com> Date: Wed, 3 Jul 2024 11:48:36 +0200 Subject: [PATCH] Add advection-diffusion as another usage example (#1138) Co-authored-by: Tom Gustafsson --- docs/examples/ex50.py | 113 +++++++++++++++++++++++ docs/examples/meshes/cylinder_stokes.msh | Bin 0 -> 32931 bytes docs/listofexamples.rst | 18 ++++ tests/test_examples.py | 8 ++ 4 files changed, 139 insertions(+) create mode 100644 docs/examples/ex50.py create mode 100644 docs/examples/meshes/cylinder_stokes.msh diff --git a/docs/examples/ex50.py b/docs/examples/ex50.py new file mode 100644 index 00000000..b56e89cd --- /dev/null +++ b/docs/examples/ex50.py @@ -0,0 +1,113 @@ +# Copyright (C) 2024 Radost Waszkiewicz and Jan Turczynowicz +# This software is published under BSD-3-clause license + +# Folowing code solves advection-diffusion problem for +# temperature distribution around a cold sphere in warm +# liquid. Liquid flow is modelled using Stokes flow field. +# Thermal difusivity to advection ratio is controlled by +# Peclet number. + +import numpy as np + +from pathlib import Path +from skfem import MeshTri, Basis, ElementTriP1, BilinearForm +from skfem import asm, solve, condense +from skfem.helpers import grad, dot + +# Define the Peclet number +peclet = 30 + +# +# Code for creating mesh which we load +# + +# floor_depth = 5.0 +# floor_width = 5.0 +# ball_size = 1.0 +# ball_segments = 100 +# mesh_size = 0.01 +# far_mesh = 0.5 + +# box_points = [ +# ([0, -floor_depth], far_mesh), +# ([floor_width, -floor_depth], far_mesh), +# ([floor_width, floor_depth], far_mesh), +# ([0, floor_depth], mesh_size), +# ] + +# phi_values = np.linspace(0, np.pi, ball_segments) +# ball_points = ball_size * np.column_stack((np.sin(phi_values), np.cos(phi_values))) +# mesh_boundary = np.vstack(( +# np.array([p for p,s in box_points]) +# , ball_points)) + +# # Create the geometry and mesh using pygmsh +# with pygmsh.geo.Geometry() as geom: +# poly = geom.add_polygon( +# mesh_boundary, +# mesh_size=([s for p,s in box_points]) + ([mesh_size] * len(ball_points)), +# ) + +# raw_mesh = geom.generate_mesh() + +# # Convert the mesh to a skfem MeshTri object and define boundaries +# mesh = MeshTri( +# raw_mesh.points[:, :2].T, raw_mesh.cells_dict["triangle"].T +# ).with_boundaries( +# { +# "left": lambda x: np.isclose(x[0], 0), # Left boundary condition +# "bottom": lambda x: np.isclose(x[1], -floor_depth), # Bottom boundary condition +# "ball": lambda x: x[0] ** 2 + x[1] ** 2 < 1.1 * ball_size**2, +# } +# ) + +mesh = MeshTri.load(Path(__file__).parent / 'meshes' / 'cylinder_stokes.msh') + +# Define the basis for the finite element method +basis = Basis(mesh, ElementTriP1()) + + +@BilinearForm +def advection(k, l, m): + """Advection bilinear form.""" + + # Coordinate fields + r, z = m.x + + u = 1 # velocity scale + a = 1 # ball size + + # Stokes flow around a sphere of size `a` + w = r ** 2 + z ** 2 + v_r = ((3 * a * r * z * u) / (4 * w**0.5)) * ((a / w) ** 2 - (1 / w)) + v_z = u + ((3 * a * u) / (4 * w**0.5)) * ( + (2 * a**2 + 3 * r**2) / (3 * w) - ((a * r) / w) ** 2 - 2 + ) + + return (l * v_r * grad(k)[0] + l * v_z * grad(k)[1]) * 2 * np.pi * r + + +@BilinearForm +def claplace(u, v, w): + """Laplace operator in cylindrical coordinates.""" + r = abs(w.x[1]) + return dot(grad(u), grad(v)) * 2 * np.pi * r + + +# Identify the interior degrees of freedom +interior = basis.complement_dofs(basis.get_dofs({"bottom", "ball"})) + +# Assemble the system matrix +A = asm(claplace, basis) + peclet * asm(advection, basis) + +# Boundary condition +u = basis.zeros() +u[basis.get_dofs("bottom")] = 1.0 +u[basis.get_dofs("ball")] = 0.0 + +u = solve(*condense(A, x=u, I=interior)) + +if __name__ == "__main__": + + mesh.draw(boundaries=True).show() + basis.plot(u, shading='gouraud', cmap='viridis').show() diff --git a/docs/examples/meshes/cylinder_stokes.msh b/docs/examples/meshes/cylinder_stokes.msh new file mode 100644 index 0000000000000000000000000000000000000000..3d9481df301b3ec75cd4a28d735ee30baf02564b GIT binary patch literal 32931 zcmeI5b$k_PyRIil0wICm9#Y&bXpkA4;O^ew4#8am6nA&GmKG@#30expi&MNnDYCHQ zuruVoXWifG{(j%@oU_l_{q|w~v99ZR-dVHWd2QA^K%22@>z?h(ckR)+Ww>vF-1&m? z1r_o2G7Q5vW4SJ2QkJn=*Ra+-KNUrl8td%0sF!`H{bC#=I{b$mk6~F%91{zS4SIt< zU>qc>|hQsCzuP&4dwy!g89JwU@%w!EC?0?3xh?#qF@ME3@i?o084@< zXo01`(qI{|ELaXK4^{vxf}vm~urgQ$tO`~GtAjPbnqV!kHdqI&3)TbcgAKrjU?Z?G z*aU0}HUpc3Ex?vwEAUG&3~UXy0o#J@!1iDVup`(B>u@lfcQ~6mTjy z4V(_n0B3@;z}es&a4t9xoDVJl7lMnx#o!WfDYy(=4z2)Kf~&yQ;2Ll(xDNbBB#HEh zt?2ll=28BSpE8Uj|0M6?v^Mrk=>7h#&uO!Mjbrn&Z+_(D<2Cr<)x#<`%pI<FJ_&2$;m%1T_o1R^^eU2b_Xp^1)_?LW*-|Vyx+Sh&jmZx?;!C&$?ua)`5r}D_?nD0rf**><% ziDV|4NvxkUe6qjfTMV)D?)rZ^iv8^~Z&!58_as(OZTow9eeQh4J$9bg=gtqlZ>M>E z?tJm+`diFY-qG8d9}E71Hb38IaooS;kN@7|-t=xU|B|1(#jF~4{ady($_-<^ zf91|?Mm_mUzMId}Q^QtW|4ZJ$tAL*^i7bQ;+NH3S^>xxlkxw@!Kg*qowqIM=jAPmF z`iee2mZdq{Pv-rUu|B=aTXg;v>T4J&t}HvZ>Gi_mIlpY5uH(6vdwA%I^dDbNO*glA zSd$I!%GER?nO}RWO@;a!rWZG=_@6CaDz=k9G|(%~&Y%g!SBF}m?~nRA`66S!DfjyL z;Nm5-xZ;!?9^mAg4a*<1#l7FnG$H#KYgyXF`W8z@bjuR|EG`U#DSwiYadMLoIm@Zw%-(c5@2=m3QC`| z;w$I;t&5-8Q)X;{)i1cv;r?%4JLlh4_vQQ^KC!Inj|<+&zc+z%{%OnV-|N5ivKgy{ zUGezN`DYJ%7I;4WS@Xaz4cp8v^3FL=WaL}d&OrfI+;}_BT$vQlIo_{wy!DB<-80V) zaiu63;U5`wvQMw?m&fn(tXqJ!`Iv!*68?I z`qw*s*U5V?8khRnsJK?DJ~t2FOMBTl|KK|9-xPix+j>8L#S!o1Tn`LmbYQ@WZTq~e zA@Kqa9{lNrbDqq1&%b++!rQtNqr%!_F}bc7My=f$e%-Oz%j!2c@59_r-Rs`+`+;@h z9eiu%FOZ;jk*4nBBilVYTVcx^bM>sTTbGuIVX@tM^k12|a5pcj@Vbrf%LTaSMU47e zei-%pYjee&1Qn)Maj&06A3xi+CY6`frSX)EiBtPJ=Q$Kw>3P3?FU|bHIqC(Laj!pn zJM7s4&4XkeQLz7O5!3BR^*el_2CR^@BvmIZfjIoCU*>!W4~XT-MF#PJR~liuGs zf9U!=y(TY?W7SU`J2h}REl%xWU(AC%?BxeOucE2`*%o_Y|U$jTz_izogdL4G}L{+3k&}+_PiPO z`{elEu5rsfch1vj;<{L=k3TghPueng$8`62@EvjUc<`>9=F(I9-o|L~)Vbb(OZ)bZ zd-2HhYH+6NjLOfQ{Ju>WYb2d<+Z>c<)zQ<<-0RhxCsR)JnexgU+^C~}%3$~Ju99zS zEvq!bWu8vl^X;^xZ=LhREtVns;5C1m0UxKmzA^cQllRS0GDeQp-j>%F!5x;*bN^1b zQa0c856@htPmxEfJ9g&%)-VDWUE7!IPZk-w>K7-UZPl?CaX0;G4yZFG zZHm+ho&FjbJA3*H`vR!xZjeyW`o3YcI90Y!!qz%7MeNd z1z&$_l~2yry>q+oBkA|9jAQJ)VRm1aEdH)>?r~LMW%UHtQ+ivuX5SwFVZQtL4b{CD z?u}YUE?lkH;?_%dzW*P0Ue4I>GTQ_wp+ZQJ{{AkYJIcDqSaqiy( z2VSJT{xJ1(v-Q4lm;FN??cw#_zS?t7c%7%_>QhO7_4(61@6uH0)umvT>t@FnM>b~O z?;a0Z%14erTls}KE#K8kF&6*n)U$q-ny&_T^s@SnS+zR;uDFpre}|cf$nww=k(c=mvH%6sZ#hQeVqTJ(>^`MPCfNL!)-I`>r8%o7rK9U#W^>o_5t71 zX1$>w*QGz|p1%dRHVzBR8`l~-;z04>6FZ&uIhuXep){4ftow1dKh99?H|P8VcZ>>~ zmFk%J^n8Zz+KznY96x4p$FJv2eq?rf`D+fXn7Oa}x{lndd}Mm3gjVQ*BSB#e z-R;@o$M*Mq&)zaG?;dz<;b!-Iys*6A=<;38n&-A=ix|@Qopb(4iCYvpJ>?HGQ`Ne^ z7aZ=s?^JvAd%2d05?V{68TqM_QExp9tVy<}%vrPRq&#=>9rm)ApZlz1sqo9bR-4?3 zrv7os{rhL$v9J#@oA_HPUeEPw65`&clpgU*wZn~~=8w7V+luqt`^`hm4*#5P>Pyo< zSJ7qV?z`)`?#pg{SKc+OT}|>PzFq8zQ~z(D{1C5f$#>?MrrR$K4|A{Qmx3?HuNZ#b zELb4T^vH?s-(5eRtvxVNR#!jz z-pmTT66Xk5zHOMCf3xRylHy-h)^*G84 zBlhgy#K_kKJ?oJK`BY#Ja50z^;Q%zp3Vaa=~kyl&6y%!KkUK|Yusj5sq` z3-tV+EXaHI&m82lg27;Ekn2b`#Jqmb{^fk;JmY+tjL+o&*}h#6vwd?SW;?PU2O^&f z>g68X(uR5X22Y zw##C~#XzJlGU$0dkyIK+Jk`{;~cQ zQO@&J4E;&Cs z;CRjtwp%jfIbM1r=6FtnI4#I=$#IhqdDfHjgX4nZj^mv3zcxO{{%MZ51;}>eIAecv zeXon-c>TQYxyZA<*-wp;p9$9XtjiXa~j4h2~s zwlmw4?bZax@w$5<=Jm7PQX|iP;`cN|-U8W<{5;z`2<2aZg}@eIey}gdaoQ4b0g&U2 z=Wm6)=XnQu9G@O>2G9hiM7Kk_<$aSPG;t+5UI1l8y$9c32 zdCs5qh}pkG5l;X)-;yHcI4Oyk_odv3^MED5PGA>sB)A&v46>cs{%kjXj_VAsgZ1Y; z!@&myBV0>?e`tibG05K$yuWe0 z&Otf*$#dRu-jBm^@xWlv0>1)%!2+NKvb{NPd4E`ja;`Hi5OZ9x9VQ~511NKjKl=FT*7ctwn31YVMJj6a=UXb&6 zKH@kYV1XBh*)Ky;&g)!-nElZVadU7r$o{N_nB#2? zVqQP{3AEz;571Fc-*mlI!Q!$a{f# zK#uEg5POa*_TMIyd#;-t=ij0{2ACA&xc<(gJc&p7X2cmmj&rW7Taah{6C!5)w<6|x zP#rhyB3z-;VNda45*}xdU+M zJo2Z&QeYXd9C#Y^26-Lb5uX8BZ=OE_`CmX@H|Ies#f0m%mfkzE@D50Bj3R27d+Fp3@Mo z1+Rjtp?SeQ{CV7@r#kE(Q6o-#puIAyD9?wO_2zi7m=NXd@%dRG z*WW}a&w;$>`kol&p7*^KI6eu=2O=H>vOFouDdWIqS!H$$BP7Ij^7V8m~VE z$~kX2PB}hPqMYp?f$wGer}Fqb=RH548s#1Exn&;pPlIyLddiQpxm=vIDayt+;e?ii|^y# zL&}br>lVvD{Y59|Z#ePqI+bl1-qJUpZ~C4c+qCW+(ki4=>o(!(qY`;@=k*Qt%^PL; ziw1M*emZjuoZm)F@mS(qL+zBeHF#N^2p4q7ck`EA1D4s|> zu{iH*c8d29+mqREcqOMqF@<t z&nnJ)iJjv3u|0?VhF4BX6myB^7U$Y+r}#V4cCOjB=eOVR3Z@t?vw-7be}ilxQ57TPJ!TieUoZ+MlbM6rT6 z*I7Fi>f~K!CC3e;vN+d0J5|-myUc2i8%A~U8sas@Yl+tuuOnVpyq?ht|e1JHAH`pou zcCdZ0{f5^NN)(5R4-+3QK0xJ&$r_*3y`;?Koji2o`6Qv8+pYw$ zbFXhZ_w}}O4{tm7>$Y=mZaeqkwsX&IJNMVNb1!W>_szC*k8C^l!?ttpYdiP3wsTKw zJNK`)bFXSU_ocRT4{AI2o3?XrX*>6kwsX&DJNJjSb1!JS?)yIN_Ws${uT|IZ@UETz zj+?)z&bO}wFMmoDpC5i>DAzDV|C^wRjrwwBmu{>BNJ? zzYtF^oE{cp>q^;zh)Z ziie076E7}aLcFB7DQ=0E5-%-YM!c+eIq~x16~rrwhl*DcuPk20ahF-ual@!4UR}I~ zcumJ$W-Z4JqqcY*@w(#m#OsSU5N{~nNW8Il6Y-|v&BU9Fw-9eB-b(yS@i6h$;%&s+ zinkMQFWy1Cqj)Fr&f;CfyNY)c?=Id$yr+1$crWqZ;(f&XiuV)mFFrtgp!gv1!Qw;2 zhl&pqA1*$^ahEyLal;rTK3Y6Ne2n;5@p0l`iH{ebAU;uilK5otDdJPbr-@G&pCLX| ze3tlZ@j2pi#pj977hfR0P<)a2V(}&7OU0LoFBe}SzEXUZ_-gSr;%mj%iLV#mAih!j zYw>TyH;I2M{+;74bF<@yvBhzhxz%yQ*e3qH_;&Ff;y;M*6yGJjTYQgrr1)O(ed7DY z4~QQWKP3L6_+jy%#E*y{75`cMnD}w=6XGYuPl=xvKO_E&_*wCD;^)OLh+h=HBz{@^ ziukYMSH-W1Ul+e2epCFG_-*mu#P5jT6~8C`yZC+aKg1u1KNNo?{#e{4{zUw#_%re6 z;xEMi6n`oHO8m9>8}YZ|@5JAWe-QsD{=dKZGmOuhzvzbix7z&mXy3No-|`RsKJWcF zEz?t?m_aE{cp>q^;zh)Z ziie07bKGSXcib>aIPNk_ikpuAcm2ViVN{YjR2HuyURAuBcy;j_j=Rj7jz|5M*c^A6 zwH-H%I^uQ3>xtJFZy?^#ahKW1al>dV-bB2q<1VwAa2?-t)99_hHt{Ji&zZjAFzn;6Cg@r&Y@ z#4n3q5&u>Es`xeW>*AmHe$kEbQrhH|_-pYu;%~*@iN6>BApTMO^WHDIG5){#gTL^u z-9G=RUs9O?$;6Y3rw~smo=QBmcpCAv;(_An#Dm1Y5dWXQ^K+R6<$4Q=7ZxufUQ|3p zyqI`#@e<-C#Z7Tbyp(up@iO9N#mkA87q1{*Q9M-qzkcUu82|as@9+GlelD|xv`I_x zR^nee?lQw1H;mTeZ5(%*Z5=m^cH-^DJBa`L?fhJ3f2scf@qyxl#0QHH5g#f(OnkWb z2=UK5zv#x8D{V4Qe7^Vs@rB}x#21S%5nn34%yE~w+;PKLA-+<4mH2A$HR5Z<*NLwe z-yptG{A=-V#5ak5EB^1X^D~USQs;f*`^6849~3_%{-gL|@t+)bnMWKqjHBW|iysp| zE`CD%r1&ZE)8c2ue-S?`eop-J&M&$#9!i@$5`Qf25`QB8RQ#FvbMY7Ae~N$J`9(Lz z=fCn9-46e?f911fr%s=@exJ8~_Rw;fIig3{r}v!Vxx{me=Mm2f$xTYl_zruPt6jysmgX@%rKo#2bn?5^pTtM7*hZGx6r)EyVxz z*3U5h_15qI8{#*`Z;9U)|4sak_+9aP;=haE7ym>2f%rr5N8*pgUE)u~pNc;de=hz) w{7>