9. Solving Stokes, in parallel#
Starting an MPI-cluster on the local machine:
from ipyparallel import Cluster
c = await Cluster(engines="mpi").start_and_connect(n=4, activate=True)
Starting 4 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>
Generate and distribute mesh, refine in parallel:
%%px
from ngsolve import *
from netgen.occ import *
from ngsolve.krylovspace import BramblePasciakCG
from mpi4py.MPI import COMM_WORLD as comm
box = Box((0,0,0), (2,0.41,0.41))
box.faces.name="wall"
box.faces.Min(X).name="inlet"
box.faces.Max(X).name="outlet"
cyl = Cylinder((0.2,0,0.2), Y, h=0.41,r=0.05)
cyl.faces.name="cyl"
shape = box-cyl
ngmesh = OCCGeometry(shape).GenerateMesh(maxh=0.05, comm=comm)
for r in range(1): ngmesh.Refine()
mesh = Mesh(ngmesh)
print (mesh.GetNE(VOL))
[stdout:3] 36792
[stdout:0] 0
[stdout:2] 36880
[stdout:1] 35600
Using parallel preconditioners from PETSc:
%%px
import ngsolve.ngs2petsc as n2p
import petsc4py.PETSc as psc
[0:execute]
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 1
----> 1 import ngsolve.ngs2petsc as n2p
2 import petsc4py.PETSc as psc
File /Applications/Netgen.app/Contents/Resources/lib/python3.13/site-packages/ngsolve/ngs2petsc.py:3
1 import ngsolve as ngs
2 import netgen.meshing as ngm
----> 3 import petsc4py.PETSc as psc
4 from mpi4py import MPI
5 import numpy as np
ModuleNotFoundError: No module named 'petsc4py'
[3:execute]
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 1
----> 1 import ngsolve.ngs2petsc as n2p
2 import petsc4py.PETSc as psc
File /Applications/Netgen.app/Contents/Resources/lib/python3.13/site-packages/ngsolve/ngs2petsc.py:3
1 import ngsolve as ngs
2 import netgen.meshing as ngm
----> 3 import petsc4py.PETSc as psc
4 from mpi4py import MPI
5 import numpy as np
ModuleNotFoundError: No module named 'petsc4py'
[1:execute]
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 1
----> 1 import ngsolve.ngs2petsc as n2p
2 import petsc4py.PETSc as psc
File /Applications/Netgen.app/Contents/Resources/lib/python3.13/site-packages/ngsolve/ngs2petsc.py:3
1 import ngsolve as ngs
2 import netgen.meshing as ngm
----> 3 import petsc4py.PETSc as psc
4 from mpi4py import MPI
5 import numpy as np
ModuleNotFoundError: No module named 'petsc4py'
[2:execute]
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 1
----> 1 import ngsolve.ngs2petsc as n2p
2 import petsc4py.PETSc as psc
File /Applications/Netgen.app/Contents/Resources/lib/python3.13/site-packages/ngsolve/ngs2petsc.py:3
1 import ngsolve as ngs
2 import netgen.meshing as ngm
----> 3 import petsc4py.PETSc as psc
4 from mpi4py import MPI
5 import numpy as np
ModuleNotFoundError: No module named 'petsc4py'
4 errors
A stabilized \(p=1\) method for Stokes:
\[\begin{split}
K = \left( \begin{array}{cc}
A & B^T \\
B & -C
\end{array} \right)
\qquad
A = \left( \begin{array}{ccc}
A_1 & \cdot & \cdot \\
\cdot & A_1 & \cdot \\
\cdot & \cdot & A_1 \\
\end{array} \right)
\end{split}\]
%%px
V = VectorH1(mesh, order=1, dirichlet="wall|inlet|cyl")
V1 = H1(mesh, order=1, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=1)
printmaster ("ndof = ", V.ndofglobal,'+',Q.ndofglobal,'=',
V.ndofglobal+Q.ndofglobal)
u,v = V.TnT()
u1,v1 = V1.TnT()
p,q = Q.TnT()
h = specialcf.mesh_size
bfa1 = BilinearForm(InnerProduct(grad(u1),grad(v1))*dx)
bfb = BilinearForm(div(u)*q*dx).Assemble()
bfc = BilinearForm(h*h*grad(p)*grad(q)*dx).Assemble()
prea1 = Preconditioner(bfa1, "gamg") # AMG precond from PETSc
bfa1.Assemble()
# make block-diagonal A matrix:
mata = sum( [Ri.T@bfa1.mat@Ri for Ri in V.restrictions] )
prea = sum( [Ei@prea1@Ei.T for Ei in V.embeddings])
bfschur = BilinearForm(p*q*dx, diagonal=True).Assemble()
preschur = bfschur.mat.Inverse()
[0:execute]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 4
2 V1 = H1(mesh, order=1, dirichlet="wall|inlet|cyl")
3 Q = H1(mesh, order=1)
----> 4 printmaster ("ndof = ", V.ndofglobal,'+',Q.ndofglobal,'=',
5 V.ndofglobal+Q.ndofglobal)
7 u,v = V.TnT()
8 u1,v1 = V1.TnT()
NameError: name 'printmaster' is not defined
[3:execute]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 4
2 V1 = H1(mesh, order=1, dirichlet="wall|inlet|cyl")
3 Q = H1(mesh, order=1)
----> 4 printmaster ("ndof = ", V.ndofglobal,'+',Q.ndofglobal,'=',
5 V.ndofglobal+Q.ndofglobal)
7 u,v = V.TnT()
8 u1,v1 = V1.TnT()
NameError: name 'printmaster' is not defined
[1:execute]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 4
2 V1 = H1(mesh, order=1, dirichlet="wall|inlet|cyl")
3 Q = H1(mesh, order=1)
----> 4 printmaster ("ndof = ", V.ndofglobal,'+',Q.ndofglobal,'=',
5 V.ndofglobal+Q.ndofglobal)
7 u,v = V.TnT()
8 u1,v1 = V1.TnT()
NameError: name 'printmaster' is not defined
[2:execute]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 4
2 V1 = H1(mesh, order=1, dirichlet="wall|inlet|cyl")
3 Q = H1(mesh, order=1)
----> 4 printmaster ("ndof = ", V.ndofglobal,'+',Q.ndofglobal,'=',
5 V.ndofglobal+Q.ndofglobal)
7 u,v = V.TnT()
8 u1,v1 = V1.TnT()
NameError: name 'printmaster' is not defined
4 errors
%%px
gfu = GridFunction(V)
gfp = GridFunction(Q)
uin = (1.5*4*y*(0.41-y)/(0.41*0.41)*z*(0.41-z)/0.41**2,0, 0)
gfu.Set(uin, definedon=mesh.Boundaries("inlet"))
resf = (-mata * gfu.vec).Evaluate()
resg = (-bfb.mat * gfu.vec).Evaluate()
sol = BramblePasciakCG (A=mata, B=bfb.mat, C=bfc.mat, f=resf, g=resg, \
preA=prea, preS=preschur, maxit=500,
printrates='\r' if comm.rank==0 else False)
gfu.vec.data += sol[0]
gfp.vec.data += sol[1]
gfu = c[:]["gfu"]
from ngsolve import *
from ngsolve.webgui import Draw
Draw (Norm(gfu[0]), gfu[0].space.mesh, clipping={"y":1, "z":0}, order=1)
Draw (Norm(gfu[2]), gfu[2].space.mesh, clipping={"y":1, "z":0}, order=1);
gfp = c[:]["gfp"]
Draw (gfp[0], order=1);