Eigenvalue Problem (NGS)¶
In this notebook we study how to solve an eigenvalue problem using NGSolve and SLEPC, in parallel.
Note to run the code in parallel using a notebook you first need to initialize ipcluster, using the command: ipcluster start –engines=MPI -n 4
[1]:
from ipyparallel import Client
c = Client()
c.ids
[1]:
[0, 1]
[2]:
%%px
#STDLIB
from math import pi
#NGSOLVE/NETGEN
from netgen.geom2d import unit_square
from ngsolve import *
from numsa.NGSlepc import *
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
npro = comm.size
H = []
E = []
for k in range(2,4):
print (rank, npro)
print("Mesh N. {}".format(k))
h = 2**(-k);
E = E + [h];
if comm.rank == 0:
ngmesh = unit_square.GenerateMesh(maxh=h).Distribute(comm)
else:
ngmesh = netgen.meshing.Mesh.Receive(comm)
mesh = Mesh(ngmesh)
fes = H1(mesh, order=1, dirichlet=".*")
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes)
a += grad(u)*grad(v)*dx
m = BilinearForm(fes)
m += u*v*dx
a.Assemble()
m.Assemble()
EP = SLEPcEigenProblem("GHEP","krylovschur")
EP.SpectralTransformation("sinvert")
PC = EP.KSP.getPC();
PC.setType("lu");
PC.setFactorSolverType("mumps");
EP.setOperators([a.mat,m.mat],fes.FreeDofs())
EP.setWhich(1);
EP.Solve()
lam, gfur, gfui = EP.getPairs(fes)
print([l/pi**2 for l in lam])
E = E + [abs(lam[0]/pi**2-2)];
[stdout:1] 1 2
Mesh N. 2
[(2.240498967068547+0j)]
1 2
Mesh N. 3
[(2.0458028065043323+0j)]
[stdout:0] 0 2
Mesh N. 2
[(2.240498967068547+0j)]
0 2
Mesh N. 3
[(2.0458028065043323+0j)]
[3]:
%%px
E
Out[0:13]: [0.25, 0.24049896706854712, 0.125, 0.04580280650433233]
Out[1:13]: [0.25, 0.24049896706854712, 0.125, 0.04580280650433233]
[4]:
%%px
#STDLIB
from math import pi
#NGSOLVE/NETGEN
from netgen.geom2d import unit_square
from ngsolve import *
import ngsolve.ngs2petsc as N2P
#MPI
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
npro = comm.size
if rank == 0:
print (rank, npro)
if comm.rank == 0:
ngmesh = unit_square.GenerateMesh(maxh=0.05).Distribute(comm)
else:
ngmesh = netgen.meshing.Mesh.Receive(comm)
mesh = Mesh(ngmesh)
fes = H1(mesh, order=1, dirichlet=".*")
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes)
a += grad(u)*grad(v)*dx
pre = Preconditioner(a, type="direct", inverse="masterinverse")
m = BilinearForm(fes)
m += u*v*dx
a.Assemble()
m.Assemble()
EP = SLEPcEigenProblem("GHEP","krylovschur")
EP.SpectralTransformation("sinvert")
PC = EP.KSP.getPC();
PC.setType("lu");
PC.setFactorSolverType("mumps");
EP.setOperators([a.mat,m.mat],fes.FreeDofs())
EP.setWhich(5);
EP.Solve()
lam, gfur, gfui = EP.getPairs(fes)
if rank == 0:
print([l/pi**2 for l in lam])
eig1 = GridFunction(fes); eig1.vec.data = gfur.vecs[0].data;
eig2 = GridFunction(fes); eig2.vec.data = gfur.vecs[1].data;
eig3 = GridFunction(fes); eig3.vec.data = gfur.vecs[2].data;
eig4 = GridFunction(fes); eig4.vec.data = gfur.vecs[3].data;
[stdout:0] 0 2
[(2.0063458676628505+0j), (5.039046175256581+0j), (5.0397702505874244+0j), (8.102182495028751+0j), (10.152804821204073+0j)]
[5]:
from ngsolve.webgui import Draw
eig1 = c[:]["eig1"]
eig2 = c[:]["eig2"]
eig3 = c[:]["eig3"]
eig4 = c[:]["eig4"]
Draw (eig1[0])
Draw (eig2[0])
Draw (eig3[0])
Draw (eig4[0])
[5]:
BaseWebGuiScene