1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
|
'''
This module contains all the functions related to the SLEPc eigenvalue
solver (EPS/PEP) interface for NGSolve
'''
from petsc4py import PETSc
try:
from slepc4py import SLEPc
except ImportError:
import warnings
warnings.warn("Import Warning: it was not possible to import SLEPc")
SLEPc = None
from ngsolve import GridFunction
from ngsPETSc import Matrix, VectorMapping
class EigenSolver():
"""
This calss creates a SLEPc Eigen Problem Solver (EPS/PEP) from NGSolve
variational problem pencil, i.e.
a0(u,v)+lam*a1(u,v)+(lam^2)*a2(u,v)+ ... = 0
Inspired by Firedrake Eigensolver class.
:arg pencil: tuple containing the bilinear forms a: V x V -> K composing
the pencil, e.g. (m,a) with a = BilinearForm(grad(u),grad(v)*dx) and
m = BilinearForm(-1*u*v*dx)
:arg fes: finite element space V
:arg nev: number of requested eigenvalue
:arg ncv: dimension of the internal subspace used by SLEPc,
by Default by SLEPc.DECIDE
:arg solverParameters: parameters to be passed to the KSP solver
:arg optionsPrefix: special solver options prefix for this specific Krylov solver
"""
if SLEPc is not None:
def __init__(self, pencil, fes, nev, ncv=SLEPc.DECIDE, optionsPrefix=None,
solverParameters=None):
self.comm = PETSc.COMM_WORLD.tompi4py()
if not isinstance(pencil, tuple): pencil=tuple([pencil])
self.penLength = len(pencil)
self.fes = fes
self.nev = nev
self.ncv = ncv
self.solverParameters = solverParameters
self.optionsPrefix = optionsPrefix
options_object = PETSc.Options()
if solverParameters is not None:
for optName, optValue in self.solverParameters.items():
options_object[optName] = optValue
self.pencilMats = []
self.pencilFlags = []
for a in pencil:
self.pencilFlags += [a.flags.ToDict()]
self.pencilMats += [Matrix(a.Assemble().mat, fes).mat]
self.pencilMats[-1].setOptionsPrefix(self.optionsPrefix)
self.pencilMats[-1].setFromOptions()
self.eps = None
self.pep = None
if self.penLength > 2:
self.isEPS = False
else:
self.isEPS = True
self.setUpEPS()
def setUpEPS(self):
'''
This function setup a SLEPc EPS if the pencil has shape either (m,a)
or is simply a single matrix.
'''
self.eps = SLEPc.EPS().create()
self.eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
if self.penLength == 1:
flag0 = self.pencilFlags[0]
if "symmetric" in flag0.keys():
if flag0["symmetric"]:
self.eps.setProblemType(SLEPc.EPS.ProblemType.HEP)
else:
self.eps.setProblemType(SLEPc.EPS.ProblemType.NHEP)
else:
self.eps.setProblemType(SLEPc.EPS.ProblemType.NHEP)
self.eps.setOperators(self.pencilMats[0])
else:
flag0 = self.pencilFlags[0]
flag1 = self.pencilFlags[1]
if "symmetric" in flag0.keys() and "symmetric" in flag1.keys():
if flag0["symmetric"] and flag1["symmetric"]:
self.eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
else:
self.eps.setProblemType(SLEPc.EPS.ProblemType.GNHEP)
else:
self.eps.setProblemType(SLEPc.EPS.ProblemType.GNHEP)
self.pencilMats[0].scale(-1)
self.eps.setOperators(self.pencilMats[1], self.pencilMats[0])
self.eps.setDimensions(self.nev, self.ncv)
self.eps.setOptionsPrefix(self.optionsPrefix)
self.eps.setFromOptions()
def solve(self):
'''
This function solve the eigenprobelm
'''
self.eps.solve()
self.nconv = self.eps.getConverged()
if self.nconv == 0:
raise RuntimeError("Did not converge any eigenvalues.")
return self.nconv
def view(self):
'''
This function setup display the information about SLEPc EPS/PEP
'''
self.eps.view()
def eigenValue(self, i):
'''
This function return the eigenvalue of the eigenproblem
:arg i: index of the eigenvalue we are intrested in.
'''
lam = None
if self.isEPS:
lam = self.eps.getEigenvalue(i)
return lam
def eigenFunction(self, i):
'''
This function return the eigenfunction of the eigenproblem
:arg i: index of the eigenfunction we are intrested in.
'''
self.vecMap = VectorMapping(self.fes)
eigenModeReal = GridFunction(self.fes)
eigenModeImag = GridFunction(self.fes)
eignModePETScReal = self.pencilMats[0].createVecLeft()
eignModePETScImag = self.pencilMats[0].createVecLeft()
if self.isEPS:
self.eps.getEigenvector(i, eignModePETScReal, eignModePETScImag)
self.vecMap.ngsVec(eignModePETScReal,eigenModeReal.vec)
self.vecMap.ngsVec(eignModePETScImag,eigenModeImag.vec)
return eigenModeReal, eigenModeImag
def eigenValues(self, indeces):
'''
This function return the eigenvalues of the eigenproblem
:arg indeces: indeces of the eigenvalues we are intrested in.
'''
lams = []
if self.isEPS:
for i in indeces:
lams = lams + [self.eps.getEigenvalue(i)]
return lams
def eigenFunctions(self, indeces):
'''
This function return a multidim with
the eigenfunctions of the eigenproblem
:arg indeces: indeces of the eigenfunctions we are intrested in.
'''
self.vecMap = VectorMapping(self.fes)
eigenModesReal = GridFunction(self.fes, multidim=self.nev)
eigenModesImag = GridFunction(self.fes, multidim=self.nev)
k = 0
eignModePETScReal = self.pencilMats[0].createVecLeft()
eignModePETScImag = self.pencilMats[0].createVecLeft()
for i in indeces:
if self.isEPS:
self.eps.getEigenvector(i, eignModePETScReal, eignModePETScImag)
self.vecMap.ngsVec(eignModePETScReal,eigenModesReal.vecs[k])
self.vecMap.ngsVec(eignModePETScImag,eigenModesImag.vecs[k])
k += 1
return eigenModesReal, eigenModesImag
|
