Hello David,

I am very sorry that I mis-spoke.  Clearly there were many differences and those differences make all the difference.

Thank you very much for presenting the two methods, that helps a great deal on how to set the PETSc options while keeping the basic structure.  This is very, very helpful.

One last question I hope. Would you be interested in me putting together a demo for other firedrakers?  If yes I should be able to do it within a week or so.

Cheers, Francis


------------------
Francis Poulin
Associate Professor
Department of Applied Mathematics
University of Waterloo

email:           fpoulin@uwaterloo.ca
Web:            https://uwaterloo.ca/poulin-research-group/
Telephone:  +1 519 888 4567 x32637


________________________________________
From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Lawrence Mitchell [lawrence.mitchell@imperial.ac.uk]
Sent: Thursday, December 01, 2016 4:54 AM
To: firedrake@imperial.ac.uk
Subject: Re: [firedrake] eigenvalue problems with SLEPc

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_fenics.py
>
> 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


_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake