Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Use SubMatrices to compute eigvals in demo qgbasinmodes #2402

Closed
wants to merge 9 commits into from
Closed

Use SubMatrices to compute eigvals in demo qgbasinmodes #2402

Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
4fd6596
compute eigvals using reduced matrices
APaganini Apr 1, 2022
d60c860
new motivation to remove DirRows2
APaganini Apr 1, 2022
251e8a6
pyop2 and petsc_a.comm
APaganini Apr 1, 2022
3d06769
new ISbd.complement
APaganini Apr 1, 2022
92d65c1
import pyop2
APaganini Apr 1, 2022
6881416
typo
APaganini Apr 1, 2022
ee296bb
added blank space
APaganini Apr 1, 2022
80ec091
parallel IScomplement
APaganini Apr 6, 2022
c31b7f9
removed debug comment
APaganini Apr 6, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 15 additions & 4 deletions demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -154,13 +155,23 @@ 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.
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
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.IntType)
ISbd = PETSc.IS().createGeneral(bcnodes, comm=petsc_a.comm)
s, e = petsc_a.getOwnershipRange()
ISin = ISbd.complement(s, e)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In parallel, ISbd.complement(s, e) may rise exception:

File PETSc/IS.pyx:234, in petsc4py.PETSc.IS.complement()

Error: error code 62
[0] ISComplement() at /home/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc-complex/src/petsc/src/vec/is/is/utils/iscoloring.c:815
[0] Invalid argument
[0] Index set must be sorted
[1:execute]

Maybe need add

s, e = petsc_a.getOwnershipRange()
bcnodes = bcnodes[np.logical_and(bcnodes < e, bcnodes >= s)]
bcnodes.sort()

to delete the ghost nodes and sort the array before creating the ISbd.

A = petsc_a.createSubMatrix(ISin, ISin)
M = petsc_m.createSubMatrix(ISin, ISin)
APaganini marked this conversation as resolved.
Show resolved Hide resolved

We can declare how many eigenpairs, eigenfunctions and eigenvalues, we
want to find ::

Expand All @@ -187,7 +198,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()

Expand All @@ -208,13 +219,13 @@ 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)

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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You should also insert the Dirichlet conditions into these eigenvectors (bc.apply)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this necessary, given that bc = DirichletBC(Vcg, 0.0, "on_boundary")?


We can now list and show plots for the eigenvalues and eigenfunctions
that were found. ::
Expand Down