dolfin team mailing list archive

[Neilen Marais <nmarais@xxxxxxxxx>]

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]