Re: [firedrake] eigenvalue problems with SLEPc
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
Hello David, I have an undergraduate student working with me for another few weeks. I will get her to put together some demos and we will send them along as they progress. Sorry to ask another question but since the recent success, we were inspired to try and do something similar for the Shallow Water model. This is more complicated because it not only has the Rossby basin modes but it also has gravity wave basin modes. We have a rough version of the code below, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py At the moment we have set the Coriolis parameter to zero, so we are essentially solving a 2D wave equation but with a momentum equation and a continuity equation. I realize that I need to learn SLEPc and understand all of the options that you need for convergence but the problem that I'm having now is more fundamental. Traceback (most recent call last): File "basin_modes_sw.py", line 64, in <module> es.solve() File "SLEPc/EPS.pyx", line 1107, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:28270) petsc4py.PETSc.Error: error code 56 [0] EPSSolve() line 80 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssolve.c [0] EPSSetUp() line 268 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssetup.c [0] STSetUp() line 310 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/interface/stsolve.c [0] STSetUp_Shift() line 139 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/impls/shift/shift.c [0] KSPSetUp() line 393 in /tmp/pip-SPrpSR-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 968 in /tmp/pip-SPrpSR-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 92 in /tmp/pip-SPrpSR-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 260 in /tmp/pip-SPrpSR-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 19 in /tmp/pip-SPrpSR-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest If anyone has any ideas as to what I am doing wrong I would be very happy to try something and fix this up. In fenics we were able to get this example to work, in case someone wanted to see that code. 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
On 02/12/16 04:23, Francis Poulin wrote:
Hello David,
I have an undergraduate student working with me for another few weeks. I will get her to put together some demos and we will send them along as they progress.
Sorry to ask another question but since the recent success, we were inspired to try and do something similar for the Shallow Water model. This is more complicated because it not only has the Rossby basin modes but it also has gravity wave basin modes. We have a rough version of the code below,
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py
At the moment we have set the Coriolis parameter to zero, so we are essentially solving a 2D wave equation but with a momentum equation and a continuity equation.
I realize that I need to learn SLEPc and understand all of the options that you need for convergence but the problem that I'm having now is more fundamental.
Traceback (most recent call last): File "basin_modes_sw.py", line 64, in <module> es.solve() File "SLEPc/EPS.pyx", line 1107, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:28270) petsc4py.PETSc.Error: error code 56 [0] EPSSolve() line 80 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssolve.c [0] EPSSetUp() line 268 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssetup.c [0] STSetUp() line 310 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/interface/stsolve.c [0] STSetUp_Shift() line 139 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/impls/shift/shift.c [0] KSPSetUp() line 393 in /tmp/pip-SPrpSR-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 968 in /tmp/pip-SPrpSR-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 92 in /tmp/pip-SPrpSR-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 260 in /tmp/pip-SPrpSR-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 19 in /tmp/pip-SPrpSR-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest
I think you need mat_type='aij' as an extra argument to your assemble() call Cheers Stephan
If anyone has any ideas as to what I am doing wrong I would be very happy to try something and fix this up. In fenics we were able to get this example to work, in case someone wanted to see that code.
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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 2 Dec 2016, at 09:56, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
On 02/12/16 04:23, Francis Poulin wrote:
Hello David,
I have an undergraduate student working with me for another few weeks. I will get her to put together some demos and we will send them along as they progress.
Sorry to ask another question but since the recent success, we were inspired to try and do something similar for the Shallow Water model. This is more complicated because it not only has the Rossby basin modes but it also has gravity wave basin modes. We have a rough version of the code below,
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py
At the moment we have set the Coriolis parameter to zero, so we are essentially solving a 2D wave equation but with a momentum equation and a continuity equation.
I realize that I need to learn SLEPc and understand all of the options that you need for convergence but the problem that I'm having now is more fundamental.
Traceback (most recent call last): File "basin_modes_sw.py", line 64, in <module> es.solve() File "SLEPc/EPS.pyx", line 1107, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:28270) petsc4py.PETSc.Error: error code 56 [0] EPSSolve() line 80 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssolve.c [0] EPSSetUp() line 268 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssetup.c [0] STSetUp() line 310 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/interface/stsolve.c [0] STSetUp_Shift() line 139 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/impls/shift/shift.c [0] KSPSetUp() line 393 in /tmp/pip-SPrpSR-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 968 in /tmp/pip-SPrpSR-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 92 in /tmp/pip-SPrpSR-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 260 in /tmp/pip-SPrpSR-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 19 in /tmp/pip-SPrpSR-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest
I think you need mat_type='aij' as an extra argument to your assemble() call
Yes. Selecting matrix types finally has its own section in the manual: http://firedrakeproject.org/solving-interface.html#specifying-solution-metho... Lawrence
Thank you for both suggestions! I am happy to say that I now get a solution, not a correct one mind you, but a solution. After I solve the problem I need to split the solution to get the velocity or the free-surface vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eigenvaluer.append(lam.real) eigenvaluei.append(lam.imag) eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi The last line is not correct and yields the error Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/vector.py", line 182, in __setitem__ self.dat.data[idx] = value TypeError: 'tuple' object does not support item assignment Can someone tell me the correct syntax for this case? 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: Friday, December 02, 2016 5:05 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc
On 2 Dec 2016, at 09:56, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
On 02/12/16 04:23, Francis Poulin wrote:
Hello David,
I have an undergraduate student working with me for another few weeks. I will get her to put together some demos and we will send them along as they progress.
Sorry to ask another question but since the recent success, we were inspired to try and do something similar for the Shallow Water model. This is more complicated because it not only has the Rossby basin modes but it also has gravity wave basin modes. We have a rough version of the code below,
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py
At the moment we have set the Coriolis parameter to zero, so we are essentially solving a 2D wave equation but with a momentum equation and a continuity equation.
I realize that I need to learn SLEPc and understand all of the options that you need for convergence but the problem that I'm having now is more fundamental.
Traceback (most recent call last): File "basin_modes_sw.py", line 64, in <module> es.solve() File "SLEPc/EPS.pyx", line 1107, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:28270) petsc4py.PETSc.Error: error code 56 [0] EPSSolve() line 80 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssolve.c [0] EPSSetUp() line 268 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssetup.c [0] STSetUp() line 310 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/interface/stsolve.c [0] STSetUp_Shift() line 139 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/impls/shift/shift.c [0] KSPSetUp() line 393 in /tmp/pip-SPrpSR-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 968 in /tmp/pip-SPrpSR-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 92 in /tmp/pip-SPrpSR-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 260 in /tmp/pip-SPrpSR-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 19 in /tmp/pip-SPrpSR-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest
I think you need mat_type='aij' as an extra argument to your assemble() call
Yes. Selecting matrix types finally has its own section in the manual: http://firedrakeproject.org/solving-interface.html#specifying-solution-metho... Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hello again, I played around a bit more with my code, and after reading the documented I have restructured it so that it is more similar to what is discussed in the docs. The updated code is on the same github repo. Now I seem to be having problems plotting the function. I beileve that I have managed (very inefficiently) to store the real and imaginary parts of the free-surface heights in a proper function. See below. em_real, em_imag = Function(Z), Function(Z) u_real, u_imag = Function(Vdg), Function(Vdg) eta_real, eta_imag = Function(Vcg), Function(Vcg) eval_real, eval_imag = [], [] # Find eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() em_real, em_imag = MixedFunctionSpace(Z, vr), MixedFunctionSpace(Z, vr) for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eval_real.append(lam.real) eval_imag.append(lam.imag) u_real, eta_real = em_real.split() #u_imag, eta_imag = em_imag.split() p = plot(eta_real) p.show() When I try plotting as above I get the following error. Sorry if this is a silly mistake but I am hoping that if i can view the eigenfunction I can see whether it looks reasonable. Given the eigenvalues I don't believe it is since it differs from my calculated values in fenics by almost 200. Cheers, Francis Traceback (most recent call last): File "basin_modes_sw.py", line 85, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type ------------------ 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 Francis Poulin [fpoulin@uwaterloo.ca] Sent: Friday, December 02, 2016 7:37 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc Thank you for both suggestions! I am happy to say that I now get a solution, not a correct one mind you, but a solution. After I solve the problem I need to split the solution to get the velocity or the free-surface vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eigenvaluer.append(lam.real) eigenvaluei.append(lam.imag) eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi The last line is not correct and yields the error Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/vector.py", line 182, in __setitem__ self.dat.data[idx] = value TypeError: 'tuple' object does not support item assignment Can someone tell me the correct syntax for this case? 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: Friday, December 02, 2016 5:05 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc
On 2 Dec 2016, at 09:56, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
On 02/12/16 04:23, Francis Poulin wrote:
Hello David,
I have an undergraduate student working with me for another few weeks. I will get her to put together some demos and we will send them along as they progress.
Sorry to ask another question but since the recent success, we were inspired to try and do something similar for the Shallow Water model. This is more complicated because it not only has the Rossby basin modes but it also has gravity wave basin modes. We have a rough version of the code below,
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py
At the moment we have set the Coriolis parameter to zero, so we are essentially solving a 2D wave equation but with a momentum equation and a continuity equation.
I realize that I need to learn SLEPc and understand all of the options that you need for convergence but the problem that I'm having now is more fundamental.
Traceback (most recent call last): File "basin_modes_sw.py", line 64, in <module> es.solve() File "SLEPc/EPS.pyx", line 1107, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:28270) petsc4py.PETSc.Error: error code 56 [0] EPSSolve() line 80 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssolve.c [0] EPSSetUp() line 268 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssetup.c [0] STSetUp() line 310 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/interface/stsolve.c [0] STSetUp_Shift() line 139 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/impls/shift/shift.c [0] KSPSetUp() line 393 in /tmp/pip-SPrpSR-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 968 in /tmp/pip-SPrpSR-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 92 in /tmp/pip-SPrpSR-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 260 in /tmp/pip-SPrpSR-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 19 in /tmp/pip-SPrpSR-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest
I think you need mat_type='aij' as an extra argument to your assemble() call
Yes. Selecting matrix types finally has its own section in the manual: http://firedrakeproject.org/solving-interface.html#specifying-solution-metho... Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hello, Sorry for the bother but I am wondering if anyone had any advice on my problem. The code can be found at the following: https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py There are two problems: 1) I can't seem to plot the eigenfunctions. Maybe I didn't split them correctly? I'm not sure. 2) The eigenvalues are not at all the ones I'm looking for. I know from the fenics version that there is an eigenvalue with an imaginary part of 3.14. These might be valid but I suspect that they are very noisy and not physically interesting. Any thoughts or suggestions as to what i can try would be greatly appreciated. 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 Francis Poulin [fpoulin@uwaterloo.ca] Sent: Friday, December 02, 2016 3:41 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc Hello again, I played around a bit more with my code, and after reading the documented I have restructured it so that it is more similar to what is discussed in the docs. The updated code is on the same github repo. Now I seem to be having problems plotting the function. I beileve that I have managed (very inefficiently) to store the real and imaginary parts of the free-surface heights in a proper function. See below. em_real, em_imag = Function(Z), Function(Z) u_real, u_imag = Function(Vdg), Function(Vdg) eta_real, eta_imag = Function(Vcg), Function(Vcg) eval_real, eval_imag = [], [] # Find eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() em_real, em_imag = MixedFunctionSpace(Z, vr), MixedFunctionSpace(Z, vr) for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eval_real.append(lam.real) eval_imag.append(lam.imag) u_real, eta_real = em_real.split() #u_imag, eta_imag = em_imag.split() p = plot(eta_real) p.show() When I try plotting as above I get the following error. Sorry if this is a silly mistake but I am hoping that if i can view the eigenfunction I can see whether it looks reasonable. Given the eigenvalues I don't believe it is since it differs from my calculated values in fenics by almost 200. Cheers, Francis Traceback (most recent call last): File "basin_modes_sw.py", line 85, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type ------------------ 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 Francis Poulin [fpoulin@uwaterloo.ca] Sent: Friday, December 02, 2016 7:37 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc Thank you for both suggestions! I am happy to say that I now get a solution, not a correct one mind you, but a solution. After I solve the problem I need to split the solution to get the velocity or the free-surface vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eigenvaluer.append(lam.real) eigenvaluei.append(lam.imag) eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi The last line is not correct and yields the error Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/vector.py", line 182, in __setitem__ self.dat.data[idx] = value TypeError: 'tuple' object does not support item assignment Can someone tell me the correct syntax for this case? 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: Friday, December 02, 2016 5:05 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc
On 2 Dec 2016, at 09:56, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
On 02/12/16 04:23, Francis Poulin wrote:
Hello David,
I have an undergraduate student working with me for another few weeks. I will get her to put together some demos and we will send them along as they progress.
Sorry to ask another question but since the recent success, we were inspired to try and do something similar for the Shallow Water model. This is more complicated because it not only has the Rossby basin modes but it also has gravity wave basin modes. We have a rough version of the code below,
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py
At the moment we have set the Coriolis parameter to zero, so we are essentially solving a 2D wave equation but with a momentum equation and a continuity equation.
I realize that I need to learn SLEPc and understand all of the options that you need for convergence but the problem that I'm having now is more fundamental.
Traceback (most recent call last): File "basin_modes_sw.py", line 64, in <module> es.solve() File "SLEPc/EPS.pyx", line 1107, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:28270) petsc4py.PETSc.Error: error code 56 [0] EPSSolve() line 80 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssolve.c [0] EPSSetUp() line 268 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssetup.c [0] STSetUp() line 310 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/interface/stsolve.c [0] STSetUp_Shift() line 139 in /home/fpoulin/software/firedrake/src/slepc/src/sys/classes/st/impls/shift/shift.c [0] KSPSetUp() line 393 in /tmp/pip-SPrpSR-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 968 in /tmp/pip-SPrpSR-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 92 in /tmp/pip-SPrpSR-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 260 in /tmp/pip-SPrpSR-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 19 in /tmp/pip-SPrpSR-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest
I think you need mat_type='aij' as an extra argument to your assemble() call
Yes. Selecting matrix types finally has its own section in the manual: http://firedrakeproject.org/solving-interface.html#specifying-solution-metho... Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Francis, On 06/12/16 04:24, Francis Poulin wrote:
Hello,
Sorry for the bother but I am wondering if anyone had any advice on my problem. The code can be found at the following:
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py
There are two problems:
1) I can't seem to plot the eigenfunctions. Maybe I didn't split them correctly? I'm not sure.
From your submitted demo, it seems that you can now plot these, is that right?
2) The eigenvalues are not at all the ones I'm looking for. I know from the fenics version that there is an eigenvalue with an imaginary part of 3.14. These might be valid but I suspect that they are very noisy and not physically interesting.
I'm not sure about that. If I run a small problem and compute all the eigenvalues using scipy.linalg.eig, the maximum value of the computed real and imaginary parts are 0.108... which agrees, I think with your firedrake and fenics code. Thanks, Lawrence
Hello Lawrence, Thanks for the reply. 1) Sadly no. I commented out the lines because they weren't working. If you uncomment them you will see that plotting gives an error. Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type 2) The version that I have posted now uses lapack and when I run that i get numbers I believe. Some zero eigenavalues but the first non-zero ones have an imaginary part of 3.14, and then they increase. The fact that lapack gives the correct answers tells me that the weak form is set up correctly. The problem is that when I use krylovschur, which is fast, specify we want a target imaginary, with value 3.14, it gives the eigenvalue with the largest magnitude. I know that these are the easiest to converge to but SLEPc should allow us to pick a target. I am not sure if the options are being set correctly. opts.setValue("eps_type", "krylovschur") opts.setValue("eps_spectrum", "target imaginary") opts.setValue("eps_tol", 1e-10) opts.setValue("spectral_shift", 3.14) 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: Wednesday, December 07, 2016 9:36 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, On 06/12/16 04:24, Francis Poulin wrote:
Hello,
Sorry for the bother but I am wondering if anyone had any advice on my problem. The code can be found at the following:
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py
There are two problems:
1) I can't seem to plot the eigenfunctions. Maybe I didn't split them correctly? I'm not sure.
From your submitted demo, it seems that you can now plot these, is that right?
2) The eigenvalues are not at all the ones I'm looking for. I know from the fenics version that there is an eigenvalue with an imaginary part of 3.14. These might be valid but I suspect that they are very noisy and not physically interesting.
I'm not sure about that. If I run a small problem and compute all the eigenvalues using scipy.linalg.eig, the maximum value of the computed real and imaginary parts are 0.108... which agrees, I think with your firedrake and fenics code. Thanks, Lawrence
Hello again, Jemma was able to point out that if I defined my eigenvectors as Functions, not mixed functions, then that works much better in terms of generating plots. em_real, em_imag = Function(Z), Function(Z) #em_real, em_imag = MixedFunctionSpace(Z, vr), MixedFunctionSpace(Z, vr) I guess the outstanding problem is how to get slepc to converge to an imaginary target. 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 Francis Poulin [fpoulin@uwaterloo.ca] Sent: Wednesday, December 07, 2016 9:46 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc Hello Lawrence, Thanks for the reply. 1) Sadly no. I commented out the lines because they weren't working. If you uncomment them you will see that plotting gives an error. Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type 2) The version that I have posted now uses lapack and when I run that i get numbers I believe. Some zero eigenavalues but the first non-zero ones have an imaginary part of 3.14, and then they increase. The fact that lapack gives the correct answers tells me that the weak form is set up correctly. The problem is that when I use krylovschur, which is fast, specify we want a target imaginary, with value 3.14, it gives the eigenvalue with the largest magnitude. I know that these are the easiest to converge to but SLEPc should allow us to pick a target. I am not sure if the options are being set correctly. opts.setValue("eps_type", "krylovschur") opts.setValue("eps_spectrum", "target imaginary") opts.setValue("eps_tol", 1e-10) opts.setValue("spectral_shift", 3.14) 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: Wednesday, December 07, 2016 9:36 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, On 06/12/16 04:24, Francis Poulin wrote:
Hello,
Sorry for the bother but I am wondering if anyone had any advice on my problem. The code can be found at the following:
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_sw.py
There are two problems:
1) I can't seem to plot the eigenfunctions. Maybe I didn't split them correctly? I'm not sure.
From your submitted demo, it seems that you can now plot these, is that right?
2) The eigenvalues are not at all the ones I'm looking for. I know from the fenics version that there is an eigenvalue with an imaginary part of 3.14. These might be valid but I suspect that they are very noisy and not physically interesting.
I'm not sure about that. If I run a small problem and compute all the eigenvalues using scipy.linalg.eig, the maximum value of the computed real and imaginary parts are 0.108... which agrees, I think with your firedrake and fenics code. Thanks, Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 07/12/16 14:46, Francis Poulin wrote:
Hello Lawrence,
Thanks for the reply.
1) Sadly no. I commented out the lines because they weren't working. If you uncomment them you will see that plotting gives an error.
Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type
Aha, I see, I was confusing which code I was looking at.
2) The version that I have posted now uses lapack and when I run that i get numbers I believe. Some zero eigenavalues but the first non-zero ones have an imaginary part of 3.14, and then they increase.
The fact that lapack gives the correct answers tells me that the weak form is set up correctly. The problem is that when I use krylovschur, which is fast, specify we want a target imaginary, with value 3.14, it gives the eigenvalue with the largest magnitude. I know that these are the easiest to converge to but SLEPc should allow us to pick a target. I am not sure if the options are being set correctly.
opts.setValue("eps_type", "krylovschur") opts.setValue("eps_spectrum", "target imaginary") opts.setValue("eps_tol", 1e-10) opts.setValue("spectral_shift", 3.14)
No, I don't think this is right. If I use: opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") # opts.setValue("eps_type", "lapack") opts.setValue("eps_type", "krylovschur") opts.setValue("eps_target_magnitude", None) opts.setValue("eps_target", 3.14) opts.setValue("eps_tol", 1e-10) opts.setValue("st_type", "sinvert") Then I can converge a few eigenvalues around 3.14i Does that work for you? Cheers, Lawrence
Hello Lawrence, Thanks! I see that one major difference is you don't try and ask for the imaginary eigenavlue of magnitude 3.14 but for the eigenvalue of that magnitude. Since the eigenvalues are purely imaginary that is a clever way to get around that. Yes, that's working brilliantly. However, I now see that my plots are just plotting zero. I guess that means that I did something silly. Maybe I can get your advice? Using getEigenpair below we get petsc matrices vr an vi. If we can convert say vr to em_real, then we can split it and plot the different fields. Unfortunately, what worked in the QG code didn't seem to work here. # Find eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eval_real.append(lam.real) eval_imag.append(lam.imag) u_real, eta_real = em_real.split() #u_imag, eta_imag = em_imag.split() p = plot(eta_real) p.show() What do you suggest we try? 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: Wednesday, December 07, 2016 10:20 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc On 07/12/16 14:46, Francis Poulin wrote:
Hello Lawrence,
Thanks for the reply.
1) Sadly no. I commented out the lines because they weren't working. If you uncomment them you will see that plotting gives an error.
Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type
Aha, I see, I was confusing which code I was looking at.
2) The version that I have posted now uses lapack and when I run that i get numbers I believe. Some zero eigenavalues but the first non-zero ones have an imaginary part of 3.14, and then they increase.
The fact that lapack gives the correct answers tells me that the weak form is set up correctly. The problem is that when I use krylovschur, which is fast, specify we want a target imaginary, with value 3.14, it gives the eigenvalue with the largest magnitude. I know that these are the easiest to converge to but SLEPc should allow us to pick a target. I am not sure if the options are being set correctly.
opts.setValue("eps_type", "krylovschur") opts.setValue("eps_spectrum", "target imaginary") opts.setValue("eps_tol", 1e-10) opts.setValue("spectral_shift", 3.14)
No, I don't think this is right. If I use: opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") # opts.setValue("eps_type", "lapack") opts.setValue("eps_type", "krylovschur") opts.setValue("eps_target_magnitude", None) opts.setValue("eps_target", 3.14) opts.setValue("eps_tol", 1e-10) opts.setValue("st_type", "sinvert") Then I can converge a few eigenvalues around 3.14i Does that work for you? Cheers, Lawrence
Hi Francis, Your code doesn't assign anything to em_real and em_imag. Jemma ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 07 December 2016 15:34:50 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hello Lawrence, Thanks! I see that one major difference is you don't try and ask for the imaginary eigenavlue of magnitude 3.14 but for the eigenvalue of that magnitude. Since the eigenvalues are purely imaginary that is a clever way to get around that. Yes, that's working brilliantly. However, I now see that my plots are just plotting zero. I guess that means that I did something silly. Maybe I can get your advice? Using getEigenpair below we get petsc matrices vr an vi. If we can convert say vr to em_real, then we can split it and plot the different fields. Unfortunately, what worked in the QG code didn't seem to work here. # Find eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eval_real.append(lam.real) eval_imag.append(lam.imag) u_real, eta_real = em_real.split() #u_imag, eta_imag = em_imag.split() p = plot(eta_real) p.show() What do you suggest we try? 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: Wednesday, December 07, 2016 10:20 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc On 07/12/16 14:46, Francis Poulin wrote:
Hello Lawrence,
Thanks for the reply.
1) Sadly no. I commented out the lines because they weren't working. If you uncomment them you will see that plotting gives an error.
Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type
Aha, I see, I was confusing which code I was looking at.
2) The version that I have posted now uses lapack and when I run that i get numbers I believe. Some zero eigenavalues but the first non-zero ones have an imaginary part of 3.14, and then they increase.
The fact that lapack gives the correct answers tells me that the weak form is set up correctly. The problem is that when I use krylovschur, which is fast, specify we want a target imaginary, with value 3.14, it gives the eigenvalue with the largest magnitude. I know that these are the easiest to converge to but SLEPc should allow us to pick a target. I am not sure if the options are being set correctly.
opts.setValue("eps_type", "krylovschur") opts.setValue("eps_spectrum", "target imaginary") opts.setValue("eps_tol", 1e-10) opts.setValue("spectral_shift", 3.14)
No, I don't think this is right. If I use: opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") # opts.setValue("eps_type", "lapack") opts.setValue("eps_type", "krylovschur") opts.setValue("eps_target_magnitude", None) opts.setValue("eps_target", 3.14) opts.setValue("eps_tol", 1e-10) opts.setValue("st_type", "sinvert") Then I can converge a few eigenvalues around 3.14i Does that work for you? Cheers, Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
More helpfully, comparing it to your qg code, it looks like you need an equivalent to the line: eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi when the eigenmodes are functions on a mixed space. J ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Shipton, Jemma <j.shipton@imperial.ac.uk> Sent: 07 December 2016 15:57:33 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, Your code doesn't assign anything to em_real and em_imag. Jemma ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 07 December 2016 15:34:50 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hello Lawrence, Thanks! I see that one major difference is you don't try and ask for the imaginary eigenavlue of magnitude 3.14 but for the eigenvalue of that magnitude. Since the eigenvalues are purely imaginary that is a clever way to get around that. Yes, that's working brilliantly. However, I now see that my plots are just plotting zero. I guess that means that I did something silly. Maybe I can get your advice? Using getEigenpair below we get petsc matrices vr an vi. If we can convert say vr to em_real, then we can split it and plot the different fields. Unfortunately, what worked in the QG code didn't seem to work here. # Find eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eval_real.append(lam.real) eval_imag.append(lam.imag) u_real, eta_real = em_real.split() #u_imag, eta_imag = em_imag.split() p = plot(eta_real) p.show() What do you suggest we try? 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: Wednesday, December 07, 2016 10:20 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc On 07/12/16 14:46, Francis Poulin wrote:
Hello Lawrence,
Thanks for the reply.
1) Sadly no. I commented out the lines because they weren't working. If you uncomment them you will see that plotting gives an error.
Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type
Aha, I see, I was confusing which code I was looking at.
2) The version that I have posted now uses lapack and when I run that i get numbers I believe. Some zero eigenavalues but the first non-zero ones have an imaginary part of 3.14, and then they increase.
The fact that lapack gives the correct answers tells me that the weak form is set up correctly. The problem is that when I use krylovschur, which is fast, specify we want a target imaginary, with value 3.14, it gives the eigenvalue with the largest magnitude. I know that these are the easiest to converge to but SLEPc should allow us to pick a target. I am not sure if the options are being set correctly.
opts.setValue("eps_type", "krylovschur") opts.setValue("eps_spectrum", "target imaginary") opts.setValue("eps_tol", 1e-10) opts.setValue("spectral_shift", 3.14)
No, I don't think this is right. If I use: opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") # opts.setValue("eps_type", "lapack") opts.setValue("eps_type", "krylovschur") opts.setValue("eps_target_magnitude", None) opts.setValue("eps_target", 3.14) opts.setValue("eps_tol", 1e-10) opts.setValue("st_type", "sinvert") Then I can converge a few eigenvalues around 3.14i Does that work for you? Cheers, Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Agreed. Unfortunately, the things that we have tried haven't seemed to work. We need to assign the vr say into a vector and then split it, or maybe vice versa. Not sure what is best to do in firedrake. 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 Shipton, Jemma [j.shipton@imperial.ac.uk] Sent: Wednesday, December 07, 2016 11:01 AM To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc More helpfully, comparing it to your qg code, it looks like you need an equivalent to the line: eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi when the eigenmodes are functions on a mixed space. J ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Shipton, Jemma <j.shipton@imperial.ac.uk> Sent: 07 December 2016 15:57:33 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, Your code doesn't assign anything to em_real and em_imag. Jemma ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 07 December 2016 15:34:50 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hello Lawrence, Thanks! I see that one major difference is you don't try and ask for the imaginary eigenavlue of magnitude 3.14 but for the eigenvalue of that magnitude. Since the eigenvalues are purely imaginary that is a clever way to get around that. Yes, that's working brilliantly. However, I now see that my plots are just plotting zero. I guess that means that I did something silly. Maybe I can get your advice? Using getEigenpair below we get petsc matrices vr an vi. If we can convert say vr to em_real, then we can split it and plot the different fields. Unfortunately, what worked in the QG code didn't seem to work here. # Find eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eval_real.append(lam.real) eval_imag.append(lam.imag) u_real, eta_real = em_real.split() #u_imag, eta_imag = em_imag.split() p = plot(eta_real) p.show() What do you suggest we try? 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: Wednesday, December 07, 2016 10:20 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc On 07/12/16 14:46, Francis Poulin wrote:
Hello Lawrence,
Thanks for the reply.
1) Sadly no. I commented out the lines because they weren't working. If you uncomment them you will see that plotting gives an error.
Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type
Aha, I see, I was confusing which code I was looking at.
2) The version that I have posted now uses lapack and when I run that i get numbers I believe. Some zero eigenavalues but the first non-zero ones have an imaginary part of 3.14, and then they increase.
The fact that lapack gives the correct answers tells me that the weak form is set up correctly. The problem is that when I use krylovschur, which is fast, specify we want a target imaginary, with value 3.14, it gives the eigenvalue with the largest magnitude. I know that these are the easiest to converge to but SLEPc should allow us to pick a target. I am not sure if the options are being set correctly.
opts.setValue("eps_type", "krylovschur") opts.setValue("eps_spectrum", "target imaginary") opts.setValue("eps_tol", 1e-10) opts.setValue("spectral_shift", 3.14)
No, I don't think this is right. If I use: opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") # opts.setValue("eps_type", "lapack") opts.setValue("eps_type", "krylovschur") opts.setValue("eps_target_magnitude", None) opts.setValue("eps_target", 3.14) opts.setValue("eps_tol", 1e-10) opts.setValue("st_type", "sinvert") Then I can converge a few eigenvalues around 3.14i Does that work for you? Cheers, Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Francis, I've just pushed a branch doing as Lawrence suggests. I can plot the eigenmodes in paraview now. Not sure what they're meant to look like! Hope it works for you now... J ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 07 December 2016 16:02:18 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Agreed. Unfortunately, the things that we have tried haven't seemed to work. We need to assign the vr say into a vector and then split it, or maybe vice versa. Not sure what is best to do in firedrake. 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 Shipton, Jemma [j.shipton@imperial.ac.uk] Sent: Wednesday, December 07, 2016 11:01 AM To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc More helpfully, comparing it to your qg code, it looks like you need an equivalent to the line: eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi when the eigenmodes are functions on a mixed space. J ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Shipton, Jemma <j.shipton@imperial.ac.uk> Sent: 07 December 2016 15:57:33 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, Your code doesn't assign anything to em_real and em_imag. Jemma ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 07 December 2016 15:34:50 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hello Lawrence, Thanks! I see that one major difference is you don't try and ask for the imaginary eigenavlue of magnitude 3.14 but for the eigenvalue of that magnitude. Since the eigenvalues are purely imaginary that is a clever way to get around that. Yes, that's working brilliantly. However, I now see that my plots are just plotting zero. I guess that means that I did something silly. Maybe I can get your advice? Using getEigenpair below we get petsc matrices vr an vi. If we can convert say vr to em_real, then we can split it and plot the different fields. Unfortunately, what worked in the QG code didn't seem to work here. # Find eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eval_real.append(lam.real) eval_imag.append(lam.imag) u_real, eta_real = em_real.split() #u_imag, eta_imag = em_imag.split() p = plot(eta_real) p.show() What do you suggest we try? 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: Wednesday, December 07, 2016 10:20 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc On 07/12/16 14:46, Francis Poulin wrote:
Hello Lawrence,
Thanks for the reply.
1) Sadly no. I commented out the lines because they weren't working. If you uncomment them you will see that plotting gives an error.
Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type
Aha, I see, I was confusing which code I was looking at.
2) The version that I have posted now uses lapack and when I run that i get numbers I believe. Some zero eigenavalues but the first non-zero ones have an imaginary part of 3.14, and then they increase.
The fact that lapack gives the correct answers tells me that the weak form is set up correctly. The problem is that when I use krylovschur, which is fast, specify we want a target imaginary, with value 3.14, it gives the eigenvalue with the largest magnitude. I know that these are the easiest to converge to but SLEPc should allow us to pick a target. I am not sure if the options are being set correctly.
opts.setValue("eps_type", "krylovschur") opts.setValue("eps_spectrum", "target imaginary") opts.setValue("eps_tol", 1e-10) opts.setValue("spectral_shift", 3.14)
No, I don't think this is right. If I use: opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") # opts.setValue("eps_type", "lapack") opts.setValue("eps_type", "krylovschur") opts.setValue("eps_target_magnitude", None) opts.setValue("eps_target", 3.14) opts.setValue("eps_tol", 1e-10) opts.setValue("st_type", "sinvert") Then I can converge a few eigenvalues around 3.14i Does that work for you? Cheers, Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hello Lawrence and Jemma, Thanks a lot for the fix. That does work very well. 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 Shipton, Jemma [j.shipton@imperial.ac.uk] Sent: Wednesday, December 07, 2016 11:25 AM To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, I've just pushed a branch doing as Lawrence suggests. I can plot the eigenmodes in paraview now. Not sure what they're meant to look like! Hope it works for you now... J ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 07 December 2016 16:02:18 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Agreed. Unfortunately, the things that we have tried haven't seemed to work. We need to assign the vr say into a vector and then split it, or maybe vice versa. Not sure what is best to do in firedrake. 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 Shipton, Jemma [j.shipton@imperial.ac.uk] Sent: Wednesday, December 07, 2016 11:01 AM To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc More helpfully, comparing it to your qg code, it looks like you need an equivalent to the line: eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi when the eigenmodes are functions on a mixed space. J ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Shipton, Jemma <j.shipton@imperial.ac.uk> Sent: 07 December 2016 15:57:33 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, Your code doesn't assign anything to em_real and em_imag. Jemma ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Francis Poulin <fpoulin@uwaterloo.ca> Sent: 07 December 2016 15:34:50 To: firedrake Subject: Re: [firedrake] eigenvalue problems with SLEPc Hello Lawrence, Thanks! I see that one major difference is you don't try and ask for the imaginary eigenavlue of magnitude 3.14 but for the eigenvalue of that magnitude. Since the eigenvalues are purely imaginary that is a clever way to get around that. Yes, that's working brilliantly. However, I now see that my plots are just plotting zero. I guess that means that I did something silly. Maybe I can get your advice? Using getEigenpair below we get petsc matrices vr an vi. If we can convert say vr to em_real, then we can split it and plot the different fields. Unfortunately, what worked in the QG code didn't seem to work here. # Find eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eval_real.append(lam.real) eval_imag.append(lam.imag) u_real, eta_real = em_real.split() #u_imag, eta_imag = em_imag.split() p = plot(eta_real) p.show() What do you suggest we try? 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: Wednesday, December 07, 2016 10:20 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc On 07/12/16 14:46, Francis Poulin wrote:
Hello Lawrence,
Thanks for the reply.
1) Sadly no. I commented out the lines because they weren't working. If you uncomment them you will see that plotting gives an error.
Traceback (most recent call last): File "basin_modes_sw.py", line 86, in <module> p = plot(eta_real) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.py", line 125, in plot raise TypeError("Unexpected type") TypeError: Unexpected type
Aha, I see, I was confusing which code I was looking at.
2) The version that I have posted now uses lapack and when I run that i get numbers I believe. Some zero eigenavalues but the first non-zero ones have an imaginary part of 3.14, and then they increase.
The fact that lapack gives the correct answers tells me that the weak form is set up correctly. The problem is that when I use krylovschur, which is fast, specify we want a target imaginary, with value 3.14, it gives the eigenvalue with the largest magnitude. I know that these are the easiest to converge to but SLEPc should allow us to pick a target. I am not sure if the options are being set correctly.
opts.setValue("eps_type", "krylovschur") opts.setValue("eps_spectrum", "target imaginary") opts.setValue("eps_tol", 1e-10) opts.setValue("spectral_shift", 3.14)
No, I don't think this is right. If I use: opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") # opts.setValue("eps_type", "lapack") opts.setValue("eps_type", "krylovschur") opts.setValue("eps_target_magnitude", None) opts.setValue("eps_target", 3.14) opts.setValue("eps_tol", 1e-10) opts.setValue("st_type", "sinvert") Then I can converge a few eigenvalues around 3.14i Does that work for you? Cheers, Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 07/12/16 15:57, Shipton, Jemma wrote:
Hi Francis,
Your code doesn't assign anything to em_real and em_imag.
# Find eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(nconv): lam = es.getEigenpair(i, vr, vi) print lam eval_real.append(lam.real) eval_imag.append(lam.imag)
u_real, eta_real = em_real.split() #u_imag, eta_imag = em_imag.split()
p = plot(eta_real) p.show()
Also, em_real is a MixedFunctionSpace, not a Function. I think you want: em_real, em_imag = Function(Z), Function(Z) then you can do: for i in range(nconv): with em_real.dat.vec as vr: with em_imag.dat.vec as vi: lam = es.getEigenpair(i, vr, vi) print lam u_real, eta_real = em_real.split() plot(eta_real) Unfortunately we don't currently support plotting of vector-valued Functions (u_real) in this example. You could output them for visualisation in paraview instead: output = File("eigenmodes.pvd") u_real, eta_real = em_real.split() u_real.rename("R(u)") eta_real.rename("R(eta)") for i in range(nconv): ... output.write(u_real, eta_real, time=i) LAwrence
participants (5)
-
David Ham
-
Francis Poulin
-
Lawrence Mitchell
-
Shipton, Jemma
-
Stephan Kramer