From 4fd6596dbb49db087ecf6172f119eeaec173fd56 Mon Sep 17 00:00:00 2001 From: APaganini Date: Fri, 1 Apr 2022 13:48:40 +0100 Subject: [PATCH 1/9] compute eigvals using reduced matrices --- .../eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index 1e21917932..324c2b04d7 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -154,13 +154,22 @@ We define the Test Function :math:`\phi` and the Trial Function To build the weak formulation of our equation we need to build two PETSc matrices in the form of a generalized eigenvalue problem, :math:`A\psi = \lambda M\psi`. We impose the boundary conditions on the -mass matrix :math:`M`, since that is where we used integration by parts. :: +mass matrix :math:`M`, since that is where we used integration by parts. +For a more stable computation of the eigenvalues, we also create +submatrices by removing rows and columns associated with the boundary +conditions. :: 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 + bcnodes = Vcg.dof_dset.lgmap.apply(bc.nodes).astype(np.int32) + ISbd = PETSc.IS().createGeneral(bcnodes, comm=PETSc.COMM_SELF) + ISin= ISbd.complement(0, Vcg.dim()) + A = petsc_a.createSubMatrix(ISin, ISin) + M = petsc_m.createSubMatrix(ISin, ISin) + We can declare how many eigenpairs, eigenfunctions and eigenvalues, we want to find :: @@ -187,7 +196,7 @@ the PETSc parameters we previously declared. :: es = SLEPc.EPS().create(comm=COMM_WORLD) es.setDimensions(num_eigenvalues) - es.setOperators(petsc_a, petsc_m) + es.setOperators(A, M) es.setFromOptions() es.solve() @@ -214,7 +223,7 @@ solver:: and we gather the final eigenfunctions :: - eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi + eigenmodes_real.vector()[ISin.array], eigenmodes_imag.vector()[ISin.array] = vr, vi We can now list and show plots for the eigenvalues and eigenfunctions that were found. :: From d60c86089faf1260ee3145e75cd4045b1e1d33e8 Mon Sep 17 00:00:00 2001 From: APaganini Date: Fri, 1 Apr 2022 18:14:34 +0100 Subject: [PATCH 2/9] new motivation to remove DirRows2 --- demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index 324c2b04d7..7fed985222 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -155,9 +155,9 @@ To build the weak formulation of our equation we need to build two PETSc matrices in the form of a generalized eigenvalue problem, :math:`A\psi = \lambda M\psi`. We impose the boundary conditions on the mass matrix :math:`M`, since that is where we used integration by parts. -For a more stable computation of the eigenvalues, we also create -submatrices by removing rows and columns associated with the boundary -conditions. :: +To avoid spurious eigenvalues associated with the rows of :math:`M` +that are modified by the Dirichlet boundary conditions, we create +submatrices that omit these rows and the correspoing columns. :: a = beta*phi*psi.dx(0)*dx m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx From 251e8a619f962d96df3b7efb2fceabfc7bca869c Mon Sep 17 00:00:00 2001 From: APaganini Date: Fri, 1 Apr 2022 18:19:03 +0100 Subject: [PATCH 3/9] pyop2 and petsc_a.comm --- demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index 7fed985222..48e40af4b2 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -164,8 +164,8 @@ submatrices that omit these rows and the correspoing columns. :: petsc_a = assemble(a).M.handle petsc_m = assemble(m, bcs=bc).M.handle - bcnodes = Vcg.dof_dset.lgmap.apply(bc.nodes).astype(np.int32) - ISbd = PETSc.IS().createGeneral(bcnodes, comm=PETSc.COMM_SELF) + bcnodes = Vcg.dof_dset.lgmap.apply(bc.nodes).astype(pyop2.datatypes.IntTpye) + ISbd = PETSc.IS().createGeneral(bcnodes, comm=petsc_a.comm) ISin= ISbd.complement(0, Vcg.dim()) A = petsc_a.createSubMatrix(ISin, ISin) M = petsc_m.createSubMatrix(ISin, ISin) From 3d06769595cf13cecb97625550d6e84cf2c47946 Mon Sep 17 00:00:00 2001 From: APaganini Date: Fri, 1 Apr 2022 18:23:45 +0100 Subject: [PATCH 4/9] new ISbd.complement --- demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index 48e40af4b2..7ff2156d3d 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -166,7 +166,8 @@ submatrices that omit these rows and the correspoing columns. :: bcnodes = Vcg.dof_dset.lgmap.apply(bc.nodes).astype(pyop2.datatypes.IntTpye) ISbd = PETSc.IS().createGeneral(bcnodes, comm=petsc_a.comm) - ISin= ISbd.complement(0, Vcg.dim()) + (s, e), _ = petsc_a.getSizes() + ISin= ISbd.complement(s, e) A = petsc_a.createSubMatrix(ISin, ISin) M = petsc_m.createSubMatrix(ISin, ISin) From 92d65c14b9b9b8784a06ad68c83208e2d96a20eb Mon Sep 17 00:00:00 2001 From: APaganini Date: Fri, 1 Apr 2022 18:27:31 +0100 Subject: [PATCH 5/9] import pyop2 --- demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index 7ff2156d3d..6a68c70938 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -109,6 +109,7 @@ Firedrake. We import the Firedrake, PETSc, and SLEPc libraries. :: from firedrake import * from firedrake.petsc import PETSc + import pyop2 try: from slepc4py import SLEPc except ImportError: From 688141662c996e8174306ec4eaf292c12645297f Mon Sep 17 00:00:00 2001 From: APaganini Date: Fri, 1 Apr 2022 18:30:13 +0100 Subject: [PATCH 6/9] typo --- demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index 6a68c70938..6df811ebc7 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -165,7 +165,7 @@ submatrices that omit these rows and the correspoing columns. :: petsc_a = assemble(a).M.handle petsc_m = assemble(m, bcs=bc).M.handle - bcnodes = Vcg.dof_dset.lgmap.apply(bc.nodes).astype(pyop2.datatypes.IntTpye) + bcnodes = Vcg.dof_dset.lgmap.apply(bc.nodes).astype(pyop2.datatypes.IntType) ISbd = PETSc.IS().createGeneral(bcnodes, comm=petsc_a.comm) (s, e), _ = petsc_a.getSizes() ISin= ISbd.complement(s, e) From ee296bba5f619e7588623b0d20b9d38524d9fbee Mon Sep 17 00:00:00 2001 From: APaganini Date: Fri, 1 Apr 2022 18:35:09 +0100 Subject: [PATCH 7/9] added blank space --- demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index 6df811ebc7..d2e1b6b8b2 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -168,7 +168,7 @@ submatrices that omit these rows and the correspoing columns. :: bcnodes = Vcg.dof_dset.lgmap.apply(bc.nodes).astype(pyop2.datatypes.IntType) ISbd = PETSc.IS().createGeneral(bcnodes, comm=petsc_a.comm) (s, e), _ = petsc_a.getSizes() - ISin= ISbd.complement(s, e) + ISin = ISbd.complement(s, e) A = petsc_a.createSubMatrix(ISin, ISin) M = petsc_m.createSubMatrix(ISin, ISin) From 80ec0911eb0ed9547036e400afd0d832197e4b19 Mon Sep 17 00:00:00 2001 From: APaganini Date: Wed, 6 Apr 2022 17:00:17 +0100 Subject: [PATCH 8/9] parallel IScomplement --- demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index d2e1b6b8b2..e279f71805 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -167,7 +167,7 @@ submatrices that omit these rows and the correspoing columns. :: bcnodes = Vcg.dof_dset.lgmap.apply(bc.nodes).astype(pyop2.datatypes.IntType) ISbd = PETSc.IS().createGeneral(bcnodes, comm=petsc_a.comm) - (s, e), _ = petsc_a.getSizes() + s, e = petsc_a.getOwnershipRange()#(s, e), _ = petsc_a.getSizes() ISin = ISbd.complement(s, e) A = petsc_a.createSubMatrix(ISin, ISin) M = petsc_m.createSubMatrix(ISin, ISin) @@ -219,7 +219,7 @@ converge any eigenvalues at all. :: If we did, we go ahead and extract them from the SLEPc eigenvalue solver:: - vr, vi = petsc_a.getVecs() + vr, vi = A.getVecs() lam = es.getEigenpair(0, vr, vi) From c31b7f9bfee723324827a77b9c4ef3cb1842372f Mon Sep 17 00:00:00 2001 From: APaganini Date: Wed, 6 Apr 2022 17:05:27 +0100 Subject: [PATCH 9/9] removed debug comment --- demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index e279f71805..ced59f63cb 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -167,7 +167,7 @@ submatrices that omit these rows and the correspoing columns. :: bcnodes = Vcg.dof_dset.lgmap.apply(bc.nodes).astype(pyop2.datatypes.IntType) ISbd = PETSc.IS().createGeneral(bcnodes, comm=petsc_a.comm) - s, e = petsc_a.getOwnershipRange()#(s, e), _ = petsc_a.getSizes() + s, e = petsc_a.getOwnershipRange() ISin = ISbd.complement(s, e) A = petsc_a.createSubMatrix(ISin, ISin) M = petsc_m.createSubMatrix(ISin, ISin)