Hi Francis, We are always spectacularly happy when people volunteer to write documentation, so yes please! Cheers, David On Thu, 1 Dec 2016 at 15:20 Francis Poulin <fpoulin@uwaterloo.ca> wrote:
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 <(519)%20888-4567>
________________________________________ 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_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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake