On 01/12/16 00:27, Francis Poulin wrote:
Hello David,
Thanks for the idea. In your example you had your a matrix and then m matrix. I presume the second is the mass matrix.
In my case I do have a and then m. I suppose that even if I had it mixed up, I would get the inverse of the eigenvalue.
In case anyone has any ideas I am including two lines.
The first is my firedrake code which does not work.
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py
The second is a fenics version of the code which does work.
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg_feni...
To my non-expert eyes they seem to be doing the same thing but clearly they are not.
In the firedrake code you do: num_eigenvalues = 5 # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) es.setProblemType(es.ProblemType.GHEP) es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) es.setDimensions(num_eigenvalues) es.setOperators(petsc_a, petsc_m) es.solve() In the fenics code you do: PETScOptions.set("eps_gen_non_hermitian") PETScOptions.set("st_pc_factor_shift_type", "NONZERO") PETScOptions.set("eps_type", "krylov") PETScOptions.set("eps_largest_imaginary") PETScOptions.set("eps_tol", 1e-10) n = 1 eigensolver = SLEPcEigenSolver(A, M) eigensolver.solve(int(n)) So there are a bunch of differences: 1. You're looking for 5 eigenvalues, not 1, in the first snippet 2. You're looking for the smallest imaginary in the first, but the largest imaginary in the second 3. You've set the problem type to generalized hermitian in the first, but generalized non-hermitian in the second. 4. You've not set eps_type in the first, but set it to "krylov" in the second (I don't think that is a valid type, did you mean krylovschur?) If I do instead: num_eigenvalues = 5 # Number of eigenvalues # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) es.setDimensions(num_eigenvalues) es.setProblemType(es.ProblemType.GNHEP) es.setWhichEigenpairs(es.Which.LARGEST_IMAGINARY) es.setOperators(petsc_a, petsc_m) # Build Operators es.solve() # solve system I get some reasonable looking answers. Note that you can control things using PETSc options in the usual way as well: opts = PETSc.Options() opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") opts.setValue("eps_type", "krylovschur") opts.setValue("eps_largest_imaginary", None) opts.setValue("eps_tol", 1e-10) es = SLEPc.EPS().create(comm=COMM_WORLD) es.setDimensions(num_eigenvalues) es.setOperators(petsc_a, petsc_m) es.setFromOptions() es.solve() Cheers, Lawrence