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