Hi,
I'm trying to solve the eigen values of a PEC cavity using the vector
wave eqn form of Maxwell's equations. I seem to have the problem set
up right using Nedelec elements, but I can't quite get the
shift-invert working for me. The spectrum of this particular problem
has a bunch of (unphysical) numerically zero eigenvalues that have to
be discarded. I found previously that using the Arnoldi method with
shift-invert and a pretty small shift gets rid of those modes quite
reliably.
I did confirm that the correct eigenvalues (as compared to an
analytical solution) are calculated if you manually discard the
numerically close to zero eigenvalues. I also managed to jimmy the
dolfin LU solver so that I could use the ARPACK wrapper from scipy.
Using ARPACK and shift-invert I indeed get the answers I expect.
In the code below I sort the ARPACK result, since it doesn't always
provide the eigenvalues in the correct order. Also note the commented
line:
#lu = dol.LUSolver(S - sigma*M)
Using it in that form results in dolfin crashing!
Running the code below results in the following output:
Eigenvalue solver (arnoldi) converged in 1 iterations.
[0.030560806979768547, 0.039273074500013903, 0.046222796597261383,
0.057890000263477112]
[-2.7842311789427754e-15, -2.7842311789427754e-15,
-2.1423834928313568e-15, -2.1423834928313568e-15]
[0.030560806979769109, 0.039273074500013806, 0.046222796597261161,
0.057890000263476911]
The first set of eigenvalues are the ARPACK results. The next are the
unfiltered SLEPc results, and the last the filtered SLEPc results.
Thanks
Neilen
Code as below:
import dolfin as dol
from dolfin import dot, curl, inner, dx
from scipy.sparse.linalg.eigen.arpack import speigs
# Test for PETSc and SLEPc
if not dol.has_la_backend("PETSc"):
print "DOLFIN has not been configured with PETSc. Exiting."
exit()
if not dol.has_slepc():
print "DOLFIN has not been configured with SLEPc. Exiting."
exit()
class CavityDims(object): pass
cdims = CavityDims()
cdims.a, cdims.b, cdims.c = 29,23,19
# Define mesh, function space
mesh = dol.UnitCube(1,1,1)
mesh.coordinates()[:] *= [cdims.a,cdims.b,cdims.c]
V = dol.FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 4)
# Define basis and bilinear form
u = dol.TrialFunction(V)
v = dol.TestFunction(V)
m = inner(v, u)*dx # Mass form
s = dot(curl(v), curl(u))*dx # Stiffness form
# Assemble smass form
M = dol.PETScMatrix()
S = dol.PETScMatrix()
dol.assemble(m, tensor=M, mesh=mesh)
dol.assemble(s, tensor=S, mesh=mesh)
sigma = 0.01
smat = S - sigma*M
#lu = dol.LUSolver(S - sigma*M)
lu = dol.LUSolver(smat)
lu.parameters["reuse_factorization"] = True
lu.parameters["report"] = False
bb = dol.Vector(M.size(0))
xx = dol.Vector(M.size(0))
def sigma_solve(b):
bb[:] = b
lu.solve(xx, bb)
return xx[:]
M_matvec = lambda x: M*x
arpack_eigs,v = speigs.ARPACK_gen_eigs(M_matvec, sigma_solve,
M.size(0), sigma, 51, ncv=91)
# Create eigensolver
esolver = dol.SLEPcEigenSolver(S,M)
esolver.parameters["spectrum"] = "smallest real"
esolver.parameters["solver"] = "arnoldi"
esolver.parameters["spectral_shift"] = sigma
esolver.parameters["spectral_transform"] = "shift-and-invert"
esolver.solve()
eigs = [esolver.get_eigenvalue(i)[0] for i in range(4)]
filtered_eigs = []
for i in range(S.size(0)):
ev_r, ev_i = esolver.get_eigenvalue(i)
if ev_r> 0.001:
filtered_eigs.append(ev_r)
print sorted(arpack_eigs)[0:4]
print eigs[0:4]
print filtered_eigs[0:4]
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@xxxxxxxxxxxxxxxxxxx
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
