Eigenvalue solver in a Mixed Function Space

edit: I'm using the following parameters, and it seems to work now for some reason. From what I've read the Cayley preconditioner is similar to shift-and-invert but has an extra term or two in it. I think it's just a stronger preconditioner overall when you know the target value you wish to find. I hope this helps out someone else in the future!
PETScOptions.set("eps_target_magnitude")
PETScOptions.set("eps_target", 3.14)
PETScOptions.set("st_type", "cayley")

Hi all,

I'm having bit of an issue with this eigenvalue problem that I'm working on. I've been using FEniCS for a while now but I'm really stuck. Below, I've posted a reduced version of the problem I'm working on.

from dolfin import *
from mshr import *

# OPERATORS
zcross = lambda i: as_vector((-i[1],i[0]))

# BASIN & MESH
length = 1.0
width  = 1.0
resolution = 25
geometry  = Rectangle(Point(0.0, 0.0), Point(width, length))
mesh = generate_mesh(geometry, resolution)

# FUNCTION & VECTOR SPACES
DG = VectorElement("DG", mesh.ufl_cell(), 1)
CG = FiniteElement("CG", mesh.ufl_cell(), 2)
G = FunctionSpace(mesh, MixedElement((DG, CG)))

# TRIAL/TEST FUNCTIONS
(u, eta) = TrialFunctions(G)
(v, lmbda) = TestFunctions(G)

# EIGENVALUE SOLs
eigenvalues_real = []
eigenvalues_imag = []

# VARIABLES
F    = Constant("1.0")
Ro   = Constant("1.0")
Fcor = Constant("0.0")
#Fcor = Constant("1.0")

# SET-UP EQUATION
a = inner(v, Fcor*zcross(u))*dx \
    + inner(v, grad(eta))*dx \
    - inner(u, grad(lmbda))*dx

m = Ro*inner(v, u)*dx + Ro*F*lmbda*eta*dx

A = PETScMatrix()
assemble(a, tensor = A)
M = PETScMatrix()
assemble(m, tensor = M)


# SOLVE
n = 5
eigensolver = SLEPcEigenSolver(A, M)
eigensolver.parameters["spectrum"] = "target imaginary"
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
eigensolver.parameters["spectral_shift"] = 1.0
eigensolver.solve(n)

print ("Number of iterations: "), eigensolver.get_iteration_number()
nconv = eigensolver.get_number_converged()

for i in range(nconv):
    lambda_real, lambda_imag, x_real, x_imag = eigensolver.get_eigenpair(i)

    if lambda_imag == 0: # REMOVE ALL ZERO IMAGINARY EIGENVALUES
        pass

    if lambda_imag != 0:
        eigenvalues_real.append(lambda_real)
        eigenvalues_imag.append(lambda_imag)

        # DEFN FUNCTION SPACES
        eigenmodes_real = Function(G, x_real)
        eigenmodes_imag = Function(G, x_imag)
        ur, etar = eigenmodes_real.split()
        ui, etai = eigenmodes_imag.split()

# LIST EIGENVALUES
print ("Real Eigenfrequencies:"), eigenvalues_real
print ("Imaginary Eigenfrequencies:"), eigenvalues_imag

# SAVE LAST REAL
File("eigenmodes_real.pvd") << eigenmodes_real

# SAVE LAST IMAGINARY
File("eigenmodes_imag.pvd") << eigenmodes_imag

Essentially, it works for our test case when Fcor is set to zero. We obtain near zero real eigenvalues and imaginary values in the form of some constant*pi. This is expected. When we add in a curl term (in this code you can just make Fcor non-zero; I've commented out when Fcor=1) however, it fails to solve for anything. This even includes when Fcor is very small (1E-5), when should expect to see similar values to when Fcor is zero.

We believe this may have something to do with the PETSc matrix becoming anti-symmetric with this additional term, however, we're at a loss as to how to get around this. We've tried using many different SLEPc and PETSc parameters. The only thing that seems to work is using an LAPACK solver but it only solves for very coarse grids (so solutions can be found) but it fails for larger grid sizes.

If anyone could give me any direction on how to fix this, or have any ideas you'd like me to try, please let me know.

Thanks,

cekaufho

Hello,

I am having a similar issue, but I discover that in my case is not related to PETSc matrix, but with dolfin implementation of 'SLEPcEigenSolver' and with boundary conditions that it apply to the PETSc matrix.

One workaround that I am doing is convert PETSc matrix to Scipy Sparse matrix and then solving the problem using sparse.linalg.eigsn, that use the arpack algorithm.

I compared 'SLEPcEigenSolver' and Scipy 'eigsn' both with the same configuration based on the arpack algorithm and most of the tests Scipy implementation was more efficient and delivering better results and converging with few number of eigenmodes selected to be extracted.

Here is some code sample that I hope that can help you in some way:

petsc_Kt, petsc_phi = FC.PETScMatrix(), FC.PETScVector()

# assembly the matrix
FC.assemble_system(fc_jac, fc_a, A_tensor=petsc_Kt, b_tensor=petsc_phi)
mat_Kt = FC.as_backend_type(petsc_Kt).mat()

# scipy sparse matrix
sc_Kt  = sc.sparse.csr_matrix(mat_Kt.getValuesCSR()[::-1], shape=mat_Kt.size)

# get all dofs ids
np_Kt_ids = fc_sp_fun.sub.dofmap().dofs()

# implementation that I did to get the dofs where dolfin apply the bcs with bc.apply()
np_bcs_ids = get_bcs_ids(fc_bcs)

# select the free dofs
np_free_ids = np.delete(np_Kt_ids, np_bcs_ids)

# slice the stiffness matrix
sc_Kt_sld = sc_Kt[np_idx, :][:, np_idx]

# Number of eigenvalues to be extracted
eign_num = 20

# Scipy Eigen package
eig_vals, eig_vecs = sc.sparse.linalg.eigsh(sc_Kt_sld, k=eign_num, which='SA', mode='buckling')