eigenvalue problems with SLEPc
Hello, A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well. I have a good start, which can be found on github, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector. Are there any other examples that someone might be able to share that might to able to do this? 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
Dear Francis, On 30/11/16 03:23, Francis Poulin wrote:
Hello,
A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well.
I have a good start, which can be found on github,
https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py
Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector.
Doing this is strictly in the realm of SLEPc, rather than Firedrake itself. We're "just" providing the operators. The SLEPc documentation is therefore probably helpful: http://slepc.upv.es/documentation/manual.htm It describes everything in terms of the C API, but fortunately, slepc4py has somewhat useful documentation strings (type help(some_slepc_function) in ipython). The translation from C functions to slepc4py functions is quite mechanical: EPSFooBar(eps, ...) translates to: eps.fooBar(...) and so forth. Thanks, Lawrence
Hi Francis, I’m working on a project involving a SLEPc eigensolver in Firedrake so I’ve done some of what you mentioned. For instance: creating an eigensolver: es = SLEPc.EPS().create(comm=COMM_WORLD) setting type of eigenvalues: es.setWhichEigenpairs(self.eigensolver.Which.SMALLEST_REAL) setting number of eigenvalues to compute: es.setDimensions(self.num_eigenvalues) setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) getting eigenvectors (in my case the eigenvalues are all real, petsc_a is the matrix) vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lmbda = es.getEigenpair(i, vr, vi) eigenvalue[i] = lmbda.real eigenfunction[i].vector()[:] = vr Here eigenfunction is an array of Firedrake functions in the right space that has already been allocated. Hope that helps, Dave On Nov 30, 2016, at 2:02 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>> wrote: Dear Francis, On 30/11/16 03:23, Francis Poulin wrote: Hello, A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well. I have a good start, which can be found on github, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector. Doing this is strictly in the realm of SLEPc, rather than Firedrake itself. We're "just" providing the operators. The SLEPc documentation is therefore probably helpful: http://slepc.upv.es/documentation/manual.htm It describes everything in terms of the C API, but fortunately, slepc4py has somewhat useful documentation strings (type help(some_slepc_function) in ipython). The translation from C functions to slepc4py functions is quite mechanical: EPSFooBar(eps, ...) translates to: eps.fooBar(...) and so forth. Thanks, Lawrence ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited.
Hello David, Yes, that is very helpful. Thanks for getting back to me. That is very helpful! I am tried to copy what you've done, see below. But I do have a few questions. 1) Do you have a topetsc function defined like what I have at the top? 2) When you define the petsc matrix, petsc_a, how did you define that? 3) I guess you are solving a regular eigenvalue problem, or is it a generalized eigenvalue problem? If it's in Finite elements i presume it's generalized? In the case of a generalized eigenvalue problem I wonder if I should have es.setOperators(A,M) where A and M are petsc matrices. 4) If yes, then when you have vr, vi = petsc_A.getVec I guess this is solving the eigenvalue problem but should you give it an A and M? Sorry to bother you but I feel like I'm close, just missing an important detail. Cheers, Francis def topetsc(A): A.force_evaluation() return A.petscmat # Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx # Assemble Weak Form into a PETSc Matrix # Impose Boundary Conditions on Mass Matrix A = topetsc(assemble(a)) M = topetsc(assemble(m, bcs=[bc])) # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) # Set type of eigenvalues es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # Set number of eigenvalues to compute num_eigenvalues = 5 es.setDimensions(num_eigenvalues) #setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) # Fin eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lam = es.getEigenpair(i, vr, vi) eigenvalue[i] = lam.real eigenfunction[i].vector()[:] = vr ------------------ 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 David Bernstein [David.Bernstein@meliorinnovations.com] Sent: Wednesday, November 30, 2016 1:26 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, I’m working on a project involving a SLEPc eigensolver in Firedrake so I’ve done some of what you mentioned. For instance: creating an eigensolver: es = SLEPc.EPS().create(comm=COMM_WORLD) setting type of eigenvalues: es.setWhichEigenpairs(self.eigensolver.Which.SMALLEST_REAL) setting number of eigenvalues to compute: es.setDimensions(self.num_eigenvalues) setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) getting eigenvectors (in my case the eigenvalues are all real, petsc_a is the matrix) vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lmbda = es.getEigenpair(i, vr, vi) eigenvalue[i] = lmbda.real eigenfunction[i].vector()[:] = vr Here eigenfunction is an array of Firedrake functions in the right space that has already been allocated. Hope that helps, Dave On Nov 30, 2016, at 2:02 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>> wrote: Dear Francis, On 30/11/16 03:23, Francis Poulin wrote: Hello, A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well. I have a good start, which can be found on github, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector. Doing this is strictly in the realm of SLEPc, rather than Firedrake itself. We're "just" providing the operators. The SLEPc documentation is therefore probably helpful: http://slepc.upv.es/documentation/manual.htm It describes everything in terms of the C API, but fortunately, slepc4py has somewhat useful documentation strings (type help(some_slepc_function) in ipython). The translation from C functions to slepc4py functions is quite mechanical: EPSFooBar(eps, ...) translates to: eps.fooBar(...) and so forth. Thanks, Lawrence ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited.
The way I do it in my code is petsc_a = assemble(a, bcs=bcs).M.handle I have both standard and generalized problems, for the generalized problem do es.ProblemType.GHEP) and for standard: es.ProblemType.HEP) (HEP = hermitian eigenvalue problem, the actual problem type depends on your form and bcs of course) Yes, for the generalized problem the matrices are set by es.setOperators(petsc_a, petsc_m) where petsc_a and petsc_m are the stiffness and mass matrices respectively. To solve: es.solve() I am not entirely sure what the line vr, vi = petsc_a.getVecs() does other than return some initialized vectors in the right format. I copied and pasted it from somewhere, it seems to be the right thing but I have no idea why! Dave where a is a Firedrake form and bcs is a list of Dirichlet conditions On Nov 30, 2016, at 10:50 AM, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello David, Yes, that is very helpful. Thanks for getting back to me. That is very helpful! I am tried to copy what you've done, see below. But I do have a few questions. 1) Do you have a topetsc function defined like what I have at the top? 2) When you define the petsc matrix, petsc_a, how did you define that? 3) I guess you are solving a regular eigenvalue problem, or is it a generalized eigenvalue problem? If it's in Finite elements i presume it's generalized? In the case of a generalized eigenvalue problem I wonder if I should have es.setOperators(A,M) where A and M are petsc matrices. 4) If yes, then when you have vr, vi = petsc_A.getVec I guess this is solving the eigenvalue problem but should you give it an A and M? Sorry to bother you but I feel like I'm close, just missing an important detail. Cheers, Francis def topetsc(A): A.force_evaluation() return A.petscmat # Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx # Assemble Weak Form into a PETSc Matrix # Impose Boundary Conditions on Mass Matrix A = topetsc(assemble(a)) M = topetsc(assemble(m, bcs=[bc])) # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) # Set type of eigenvalues es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # Set number of eigenvalues to compute num_eigenvalues = 5 es.setDimensions(num_eigenvalues) #setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) # Fin eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lam = es.getEigenpair(i, vr, vi) eigenvalue[i] = lam.real eigenfunction[i].vector()[:] = vr ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of David Bernstein [David.Bernstein@meliorinnovations.com<mailto:David.Bernstein@meliorinnovations.com>] Sent: Wednesday, November 30, 2016 1:26 PM To: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, I’m working on a project involving a SLEPc eigensolver in Firedrake so I’ve done some of what you mentioned. For instance: creating an eigensolver: es = SLEPc.EPS().create(comm=COMM_WORLD) setting type of eigenvalues: es.setWhichEigenpairs(self.eigensolver.Which.SMALLEST_REAL) setting number of eigenvalues to compute: es.setDimensions(self.num_eigenvalues) setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) getting eigenvectors (in my case the eigenvalues are all real, petsc_a is the matrix) vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lmbda = es.getEigenpair(i, vr, vi) eigenvalue[i] = lmbda.real eigenfunction[i].vector()[:] = vr Here eigenfunction is an array of Firedrake functions in the right space that has already been allocated. Hope that helps, Dave On Nov 30, 2016, at 2:02 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>> wrote: Dear Francis, On 30/11/16 03:23, Francis Poulin wrote: Hello, A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well. I have a good start, which can be found on github, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector. Doing this is strictly in the realm of SLEPc, rather than Firedrake itself. We're "just" providing the operators. The SLEPc documentation is therefore probably helpful: http://slepc.upv.es/documentation/manual.htm It describes everything in terms of the C API, but fortunately, slepc4py has somewhat useful documentation strings (type help(some_slepc_function) in ipython). The translation from C functions to slepc4py functions is quite mechanical: EPSFooBar(eps, ...) translates to: eps.fooBar(...) and so forth. Thanks, Lawrence ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited.
Sorry, insert last line “where a is …” below second line. On Nov 30, 2016, at 11:01 AM, David Bernstein <david.bernstein@meliorinnovations.com<mailto:david.bernstein@meliorinnovations.com>> wrote: The way I do it in my code is petsc_a = assemble(a, bcs=bcs).M.handle I have both standard and generalized problems, for the generalized problem do es.ProblemType.GHEP) and for standard: es.ProblemType.HEP) (HEP = hermitian eigenvalue problem, the actual problem type depends on your form and bcs of course) Yes, for the generalized problem the matrices are set by es.setOperators(petsc_a, petsc_m) where petsc_a and petsc_m are the stiffness and mass matrices respectively. To solve: es.solve() I am not entirely sure what the line vr, vi = petsc_a.getVecs() does other than return some initialized vectors in the right format. I copied and pasted it from somewhere, it seems to be the right thing but I have no idea why! Dave where a is a Firedrake form and bcs is a list of Dirichlet conditions On Nov 30, 2016, at 10:50 AM, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello David, Yes, that is very helpful. Thanks for getting back to me. That is very helpful! I am tried to copy what you've done, see below. But I do have a few questions. 1) Do you have a topetsc function defined like what I have at the top? 2) When you define the petsc matrix, petsc_a, how did you define that? 3) I guess you are solving a regular eigenvalue problem, or is it a generalized eigenvalue problem? If it's in Finite elements i presume it's generalized? In the case of a generalized eigenvalue problem I wonder if I should have es.setOperators(A,M) where A and M are petsc matrices. 4) If yes, then when you have vr, vi = petsc_A.getVec I guess this is solving the eigenvalue problem but should you give it an A and M? Sorry to bother you but I feel like I'm close, just missing an important detail. Cheers, Francis def topetsc(A): A.force_evaluation() return A.petscmat # Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx # Assemble Weak Form into a PETSc Matrix # Impose Boundary Conditions on Mass Matrix A = topetsc(assemble(a)) M = topetsc(assemble(m, bcs=[bc])) # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) # Set type of eigenvalues es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # Set number of eigenvalues to compute num_eigenvalues = 5 es.setDimensions(num_eigenvalues) #setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) # Fin eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lam = es.getEigenpair(i, vr, vi) eigenvalue[i] = lam.real eigenfunction[i].vector()[:] = vr ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of David Bernstein [David.Bernstein@meliorinnovations.com<mailto:David.Bernstein@meliorinnovations.com>] Sent: Wednesday, November 30, 2016 1:26 PM To: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, I’m working on a project involving a SLEPc eigensolver in Firedrake so I’ve done some of what you mentioned. For instance: creating an eigensolver: es = SLEPc.EPS().create(comm=COMM_WORLD) setting type of eigenvalues: es.setWhichEigenpairs(self.eigensolver.Which.SMALLEST_REAL) setting number of eigenvalues to compute: es.setDimensions(self.num_eigenvalues) setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) getting eigenvectors (in my case the eigenvalues are all real, petsc_a is the matrix) vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lmbda = es.getEigenpair(i, vr, vi) eigenvalue[i] = lmbda.real eigenfunction[i].vector()[:] = vr Here eigenfunction is an array of Firedrake functions in the right space that has already been allocated. Hope that helps, Dave On Nov 30, 2016, at 2:02 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>> wrote: Dear Francis, On 30/11/16 03:23, Francis Poulin wrote: Hello, A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well. I have a good start, which can be found on github, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector. Doing this is strictly in the realm of SLEPc, rather than Firedrake itself. We're "just" providing the operators. The SLEPc documentation is therefore probably helpful: http://slepc.upv.es/documentation/manual.htm It describes everything in terms of the C API, but fortunately, slepc4py has somewhat useful documentation strings (type help(some_slepc_function) in ipython). The translation from C functions to slepc4py functions is quite mechanical: EPSFooBar(eps, ...) translates to: eps.fooBar(...) and so forth. Thanks, Lawrence ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited.
Hello David, Thanks again! I modified your code an think I have set it up correctly. One issue is that you define your petsc_a/m using M.handle. I don't know what that means but I presume I don't change it. Also, since I only integrate by parts on the mass matrix I impose those Dirichlet BCs on that matrix. If I impose it on both I get an eigenvalue of 1. Not sure why but there is probably a good reason to explain that. The good news is that it runs, the bad news is that it does not converge to any eigenvalues. Besides that I like the setup. Seems generic enough that one can apply it to a variety of problems. If you have any thoughts I'd be curios to hear them. Francis # Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx petsc_a = assemble(a).M.handle petsc_m = assemble(m, bcs=bc).M.handle # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) es.ProblemType.GHEP # Generalized eigenvalue problem es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # type of eigenvalues to find num_eigenvalues = 1 # Number of eigenvalues es.setDimensions(num_eigenvalues) es.setOperators(petsc_a, petsc_m) # Build Operators es.solve() # solve system nconv = es.getConverged() print nconv ------------------ 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 David Bernstein [David.Bernstein@meliorinnovations.com] Sent: Wednesday, November 30, 2016 2:03 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc Sorry, insert last line “where a is …” below second line. On Nov 30, 2016, at 11:01 AM, David Bernstein <david.bernstein@meliorinnovations.com<mailto:david.bernstein@meliorinnovations.com>> wrote: The way I do it in my code is petsc_a = assemble(a, bcs=bcs).M.handle I have both standard and generalized problems, for the generalized problem do es.ProblemType.GHEP) and for standard: es.ProblemType.HEP) (HEP = hermitian eigenvalue problem, the actual problem type depends on your form and bcs of course) Yes, for the generalized problem the matrices are set by es.setOperators(petsc_a, petsc_m) where petsc_a and petsc_m are the stiffness and mass matrices respectively. To solve: es.solve() I am not entirely sure what the line vr, vi = petsc_a.getVecs() does other than return some initialized vectors in the right format. I copied and pasted it from somewhere, it seems to be the right thing but I have no idea why! Dave where a is a Firedrake form and bcs is a list of Dirichlet conditions On Nov 30, 2016, at 10:50 AM, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello David, Yes, that is very helpful. Thanks for getting back to me. That is very helpful! I am tried to copy what you've done, see below. But I do have a few questions. 1) Do you have a topetsc function defined like what I have at the top? 2) When you define the petsc matrix, petsc_a, how did you define that? 3) I guess you are solving a regular eigenvalue problem, or is it a generalized eigenvalue problem? If it's in Finite elements i presume it's generalized? In the case of a generalized eigenvalue problem I wonder if I should have es.setOperators(A,M) where A and M are petsc matrices. 4) If yes, then when you have vr, vi = petsc_A.getVec I guess this is solving the eigenvalue problem but should you give it an A and M? Sorry to bother you but I feel like I'm close, just missing an important detail. Cheers, Francis def topetsc(A): A.force_evaluation() return A.petscmat # Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx # Assemble Weak Form into a PETSc Matrix # Impose Boundary Conditions on Mass Matrix A = topetsc(assemble(a)) M = topetsc(assemble(m, bcs=[bc])) # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) # Set type of eigenvalues es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # Set number of eigenvalues to compute num_eigenvalues = 5 es.setDimensions(num_eigenvalues) #setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) # Fin eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lam = es.getEigenpair(i, vr, vi) eigenvalue[i] = lam.real eigenfunction[i].vector()[:] = vr ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of David Bernstein [David.Bernstein@meliorinnovations.com<mailto:David.Bernstein@meliorinnovations.com>] Sent: Wednesday, November 30, 2016 1:26 PM To: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, I’m working on a project involving a SLEPc eigensolver in Firedrake so I’ve done some of what you mentioned. For instance: creating an eigensolver: es = SLEPc.EPS().create(comm=COMM_WORLD) setting type of eigenvalues: es.setWhichEigenpairs(self.eigensolver.Which.SMALLEST_REAL) setting number of eigenvalues to compute: es.setDimensions(self.num_eigenvalues) setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) getting eigenvectors (in my case the eigenvalues are all real, petsc_a is the matrix) vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lmbda = es.getEigenpair(i, vr, vi) eigenvalue[i] = lmbda.real eigenfunction[i].vector()[:] = vr Here eigenfunction is an array of Firedrake functions in the right space that has already been allocated. Hope that helps, Dave On Nov 30, 2016, at 2:02 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>> wrote: Dear Francis, On 30/11/16 03:23, Francis Poulin wrote: Hello, A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well. I have a good start, which can be found on github, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector. Doing this is strictly in the realm of SLEPc, rather than Firedrake itself. We're "just" providing the operators. The SLEPc documentation is therefore probably helpful: http://slepc.upv.es/documentation/manual.htm It describes everything in terms of the C API, but fortunately, slepc4py has somewhat useful documentation strings (type help(some_slepc_function) in ipython). The translation from C functions to slepc4py functions is quite mechanical: EPSFooBar(eps, ...) translates to: eps.fooBar(...) and so forth. Thanks, Lawrence ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited.
On 30/11/16 19:30, Francis Poulin wrote:
Hello David,
Thanks again!
I modified your code an think I have set it up correctly.
One issue is that you define your petsc_a/m using M.handle. I don't know what that means but I presume I don't change it.
Also, since I only integrate by parts on the mass matrix I impose those Dirichlet BCs on that matrix. If I impose it on both I get an eigenvalue of 1. Not sure why but there is probably a good reason to explain that.
The good news is that it runs, the bad news is that it does not converge to any eigenvalues.
Besides that I like the setup. Seems generic enough that one can apply it to a variety of problems.
If you have any thoughts I'd be curios to hear them.
Francis
# Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx
petsc_a = assemble(a).M.handle petsc_m = assemble(m, bcs=bc).M.handle
# Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) es.ProblemType.GHEP # Generalized eigenvalue problem
You need to write: es.setProblemType(es.ProblemType.GHEP) Lawrence
Hello, Sorry that I misunderstood. Now I have a problem with "Wrong value of eps" Traceback (most recent call last): File "basin_modes_qg.py", line 45, in <module> es.solve() # solve system File "SLEPc/EPS.pyx", line 1107, in slepc4py.SLEPc.EPS.solve (src/slepc4py.SLEPc.c:28270) petsc4py.PETSc.Error: error code 1 [0] EPSSolve() line 80 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssolve.c [0] EPSSetUp() line 175 in /home/fpoulin/software/firedrake/src/slepc/src/eps/interface/epssetup.c [0] EPSSetUp_KrylovSchur() line 99 in /home/fpoulin/software/firedrake/src/slepc/src/eps/impls/krylov/krylovschur/krylovschur.c [0] Wrong value of eps->which I have part of my code copied below. # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) es.setProblemType(es.ProblemType.GHEP) # Generalized eigenvalue problem es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # type of eigenvalues to find es.setDimensions(num_eigenvalues) es.setOperators(petsc_a, petsc_m) # Build Operators es.solve() # solve system I presume that there is something else that i did not set up correctly? Sorry for the bother. 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, November 30, 2016 2:33 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc On 30/11/16 19:30, Francis Poulin wrote:
Hello David,
Thanks again!
I modified your code an think I have set it up correctly.
One issue is that you define your petsc_a/m using M.handle. I don't know what that means but I presume I don't change it.
Also, since I only integrate by parts on the mass matrix I impose those Dirichlet BCs on that matrix. If I impose it on both I get an eigenvalue of 1. Not sure why but there is probably a good reason to explain that.
The good news is that it runs, the bad news is that it does not converge to any eigenvalues.
Besides that I like the setup. Seems generic enough that one can apply it to a variety of problems.
If you have any thoughts I'd be curios to hear them.
Francis
# Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx
petsc_a = assemble(a).M.handle petsc_m = assemble(m, bcs=bc).M.handle
# Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) es.ProblemType.GHEP # Generalized eigenvalue problem
You need to write: es.setProblemType(es.ProblemType.GHEP) Lawrence
Ok, I don’t know your problem but is it possible that you have the mass and stiffness matrices switched? It looks like you have defined “a” to be the mass matrix and “m” the stiffness. In that case they should be switched in the setOperators line. Also in my code I apply On Nov 30, 2016, at 11:30 AM, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello David, Thanks again! I modified your code an think I have set it up correctly. One issue is that you define your petsc_a/m using M.handle. I don't know what that means but I presume I don't change it. Also, since I only integrate by parts on the mass matrix I impose those Dirichlet BCs on that matrix. If I impose it on both I get an eigenvalue of 1. Not sure why but there is probably a good reason to explain that. The good news is that it runs, the bad news is that it does not converge to any eigenvalues. Besides that I like the setup. Seems generic enough that one can apply it to a variety of problems. If you have any thoughts I'd be curios to hear them. Francis # Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx petsc_a = assemble(a).M.handle petsc_m = assemble(m, bcs=bc).M.handle # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) es.ProblemType.GHEP # Generalized eigenvalue problem es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # type of eigenvalues to find num_eigenvalues = 1 # Number of eigenvalues es.setDimensions(num_eigenvalues) es.setOperators(petsc_a, petsc_m) # Build Operators es.solve() # solve system nconv = es.getConverged() print nconv ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of David Bernstein [David.Bernstein@meliorinnovations.com<mailto:David.Bernstein@meliorinnovations.com>] Sent: Wednesday, November 30, 2016 2:03 PM To: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Subject: Re: [firedrake] eigenvalue problems with SLEPc Sorry, insert last line “where a is …” below second line. On Nov 30, 2016, at 11:01 AM, David Bernstein <david.bernstein@meliorinnovations.com<mailto:david.bernstein@meliorinnovations.com>> wrote: The way I do it in my code is petsc_a = assemble(a, bcs=bcs).M.handle I have both standard and generalized problems, for the generalized problem do es.ProblemType.GHEP) and for standard: es.ProblemType.HEP) (HEP = hermitian eigenvalue problem, the actual problem type depends on your form and bcs of course) Yes, for the generalized problem the matrices are set by es.setOperators(petsc_a, petsc_m) where petsc_a and petsc_m are the stiffness and mass matrices respectively. To solve: es.solve() I am not entirely sure what the line vr, vi = petsc_a.getVecs() does other than return some initialized vectors in the right format. I copied and pasted it from somewhere, it seems to be the right thing but I have no idea why! Dave where a is a Firedrake form and bcs is a list of Dirichlet conditions On Nov 30, 2016, at 10:50 AM, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello David, Yes, that is very helpful. Thanks for getting back to me. That is very helpful! I am tried to copy what you've done, see below. But I do have a few questions. 1) Do you have a topetsc function defined like what I have at the top? 2) When you define the petsc matrix, petsc_a, how did you define that? 3) I guess you are solving a regular eigenvalue problem, or is it a generalized eigenvalue problem? If it's in Finite elements i presume it's generalized? In the case of a generalized eigenvalue problem I wonder if I should have es.setOperators(A,M) where A and M are petsc matrices. 4) If yes, then when you have vr, vi = petsc_A.getVec I guess this is solving the eigenvalue problem but should you give it an A and M? Sorry to bother you but I feel like I'm close, just missing an important detail. Cheers, Francis def topetsc(A): A.force_evaluation() return A.petscmat # Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx # Assemble Weak Form into a PETSc Matrix # Impose Boundary Conditions on Mass Matrix A = topetsc(assemble(a)) M = topetsc(assemble(m, bcs=[bc])) # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) # Set type of eigenvalues es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # Set number of eigenvalues to compute num_eigenvalues = 5 es.setDimensions(num_eigenvalues) #setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) # Fin eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lam = es.getEigenpair(i, vr, vi) eigenvalue[i] = lam.real eigenfunction[i].vector()[:] = vr ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of David Bernstein [David.Bernstein@meliorinnovations.com<mailto:David.Bernstein@meliorinnovations.com>] Sent: Wednesday, November 30, 2016 1:26 PM To: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, I’m working on a project involving a SLEPc eigensolver in Firedrake so I’ve done some of what you mentioned. For instance: creating an eigensolver: es = SLEPc.EPS().create(comm=COMM_WORLD) setting type of eigenvalues: es.setWhichEigenpairs(self.eigensolver.Which.SMALLEST_REAL) setting number of eigenvalues to compute: es.setDimensions(self.num_eigenvalues) setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) getting eigenvectors (in my case the eigenvalues are all real, petsc_a is the matrix) vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lmbda = es.getEigenpair(i, vr, vi) eigenvalue[i] = lmbda.real eigenfunction[i].vector()[:] = vr Here eigenfunction is an array of Firedrake functions in the right space that has already been allocated. Hope that helps, Dave On Nov 30, 2016, at 2:02 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>> wrote: Dear Francis, On 30/11/16 03:23, Francis Poulin wrote: Hello, A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well. I have a good start, which can be found on github, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector. Doing this is strictly in the realm of SLEPc, rather than Firedrake itself. We're "just" providing the operators. The SLEPc documentation is therefore probably helpful: http://slepc.upv.es/documentation/manual.htm It describes everything in terms of the C API, but fortunately, slepc4py has somewhat useful documentation strings (type help(some_slepc_function) in ipython). The translation from C functions to slepc4py functions is quite mechanical: EPSFooBar(eps, ...) translates to: eps.fooBar(...) and so forth. Thanks, Lawrence ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited._______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited.
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. Any suggestions would be greatly appreciated. If I can get this working I would be happy to put together an eigenvalues demo, if people might find this interesting, so that others can benefit from my pitfalls. 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 David Bernstein [David.Bernstein@meliorinnovations.com] Sent: Wednesday, November 30, 2016 3:50 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] eigenvalue problems with SLEPc Ok, I don’t know your problem but is it possible that you have the mass and stiffness matrices switched? It looks like you have defined “a” to be the mass matrix and “m” the stiffness. In that case they should be switched in the setOperators line. Also in my code I apply On Nov 30, 2016, at 11:30 AM, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello David, Thanks again! I modified your code an think I have set it up correctly. One issue is that you define your petsc_a/m using M.handle. I don't know what that means but I presume I don't change it. Also, since I only integrate by parts on the mass matrix I impose those Dirichlet BCs on that matrix. If I impose it on both I get an eigenvalue of 1. Not sure why but there is probably a good reason to explain that. The good news is that it runs, the bad news is that it does not converge to any eigenvalues. Besides that I like the setup. Seems generic enough that one can apply it to a variety of problems. If you have any thoughts I'd be curios to hear them. Francis # Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx petsc_a = assemble(a).M.handle petsc_m = assemble(m, bcs=bc).M.handle # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) es.ProblemType.GHEP # Generalized eigenvalue problem es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # type of eigenvalues to find num_eigenvalues = 1 # Number of eigenvalues es.setDimensions(num_eigenvalues) es.setOperators(petsc_a, petsc_m) # Build Operators es.solve() # solve system nconv = es.getConverged() print nconv ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of David Bernstein [David.Bernstein@meliorinnovations.com<mailto:David.Bernstein@meliorinnovations.com>] Sent: Wednesday, November 30, 2016 2:03 PM To: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Subject: Re: [firedrake] eigenvalue problems with SLEPc Sorry, insert last line “where a is …” below second line. On Nov 30, 2016, at 11:01 AM, David Bernstein <david.bernstein@meliorinnovations.com<mailto:david.bernstein@meliorinnovations.com>> wrote: The way I do it in my code is petsc_a = assemble(a, bcs=bcs).M.handle I have both standard and generalized problems, for the generalized problem do es.ProblemType.GHEP) and for standard: es.ProblemType.HEP) (HEP = hermitian eigenvalue problem, the actual problem type depends on your form and bcs of course) Yes, for the generalized problem the matrices are set by es.setOperators(petsc_a, petsc_m) where petsc_a and petsc_m are the stiffness and mass matrices respectively. To solve: es.solve() I am not entirely sure what the line vr, vi = petsc_a.getVecs() does other than return some initialized vectors in the right format. I copied and pasted it from somewhere, it seems to be the right thing but I have no idea why! Dave where a is a Firedrake form and bcs is a list of Dirichlet conditions On Nov 30, 2016, at 10:50 AM, Francis Poulin <fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca>> wrote: Hello David, Yes, that is very helpful. Thanks for getting back to me. That is very helpful! I am tried to copy what you've done, see below. But I do have a few questions. 1) Do you have a topetsc function defined like what I have at the top? 2) When you define the petsc matrix, petsc_a, how did you define that? 3) I guess you are solving a regular eigenvalue problem, or is it a generalized eigenvalue problem? If it's in Finite elements i presume it's generalized? In the case of a generalized eigenvalue problem I wonder if I should have es.setOperators(A,M) where A and M are petsc matrices. 4) If yes, then when you have vr, vi = petsc_A.getVec I guess this is solving the eigenvalue problem but should you give it an A and M? Sorry to bother you but I feel like I'm close, just missing an important detail. Cheers, Francis def topetsc(A): A.force_evaluation() return A.petscmat # Set up Eigenvalue Problem: QG Basin modes a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx # Assemble Weak Form into a PETSc Matrix # Impose Boundary Conditions on Mass Matrix A = topetsc(assemble(a)) M = topetsc(assemble(m, bcs=[bc])) # Create solver es = SLEPc.EPS().create(comm=COMM_WORLD) # Set type of eigenvalues es.setWhichEigenpairs(es.Which.SMALLEST_IMAGINARY) # Set number of eigenvalues to compute num_eigenvalues = 5 es.setDimensions(num_eigenvalues) #setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) # Fin eigenvalues an eigenvectors vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lam = es.getEigenpair(i, vr, vi) eigenvalue[i] = lam.real eigenfunction[i].vector()[:] = vr ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca<mailto:fpoulin@uwaterloo.ca> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> [firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>] on behalf of David Bernstein [David.Bernstein@meliorinnovations.com<mailto:David.Bernstein@meliorinnovations.com>] Sent: Wednesday, November 30, 2016 1:26 PM To: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Subject: Re: [firedrake] eigenvalue problems with SLEPc Hi Francis, I’m working on a project involving a SLEPc eigensolver in Firedrake so I’ve done some of what you mentioned. For instance: creating an eigensolver: es = SLEPc.EPS().create(comm=COMM_WORLD) setting type of eigenvalues: es.setWhichEigenpairs(self.eigensolver.Which.SMALLEST_REAL) setting number of eigenvalues to compute: es.setDimensions(self.num_eigenvalues) setting the matrix (petsc_a is a PETSc matrix) es.setOperators(petsc_a) getting eigenvectors (in my case the eigenvalues are all real, petsc_a is the matrix) vr, vi = petsc_a.getVecs() for i in range(num_eigenvalues): lmbda = es.getEigenpair(i, vr, vi) eigenvalue[i] = lmbda.real eigenfunction[i].vector()[:] = vr Here eigenfunction is an array of Firedrake functions in the right space that has already been allocated. Hope that helps, Dave On Nov 30, 2016, at 2:02 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>> wrote: Dear Francis, On 30/11/16 03:23, Francis Poulin wrote: Hello, A while ago someone pointed me towards test_slepc.py for an example on how to solve an eigenvalue problem. I have been using this as a model to compute Quasi-Geostrophic basin modes in a closed basin. We have a version working in FENicS and I am very keen to get it working in Firedrake as well. I have a good start, which can be found on github, https://github.com/francispoulin/firedrakeQG/blob/master/basin_modes_qg.py Unfortunately, I don't know how to do a lot of things. For example specify number of eigenvalues to determine, specify a target eigenvalue, what method to use, or even just look at the eigenvector. Doing this is strictly in the realm of SLEPc, rather than Firedrake itself. We're "just" providing the operators. The SLEPc documentation is therefore probably helpful: http://slepc.upv.es/documentation/manual.htm It describes everything in terms of the C API, but fortunately, slepc4py has somewhat useful documentation strings (type help(some_slepc_function) in ipython). The translation from C functions to slepc4py functions is quite mechanical: EPSFooBar(eps, ...) translates to: eps.fooBar(...) and so forth. Thanks, Lawrence ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited._______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake ©Melior Innovations Inc, 2016 all rights reserved. This e-mail is covered by the Electronic Communication Privacy Act, 18 U.S.C. Section 2510-2521 and may be legally privileged. If you have received this transmission in error, please notify the original sender immediately by return e-mail or by telephone at the telephone number provided above and delete/trash the original message from your system. Thank you for your assistance. The information contained herein is confidential information provided only for the use of the individual or entity for whom it was intended. If the reader of this message received it in error, is not the intended recipient, or if an attachment is made in error, the reader is hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited.
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
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_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
participants (3)
-
David Bernstein
-
Francis Poulin
-
Lawrence Mitchell