Use SubMatrices to compute eigvals in demo qgbasinmodes

[APaganini]

To compute the eigenvalues more robustly, create submatrices by removing rows and columns associated with degrees of freedom. This trick was suggested (on a different example) by @jandrej , @wence- , @francesco-ballarin and Matt Knepley, among others.

Disclaimer I haven't run the demo to check that it runs because I don't know how to convert an *.py.rst into a *.py I can run.

[jandrej]

Sweet. I like the complement approach of getting the dofs.

[APaganini]

Sweet. I like the complement approach of getting the dofs.

It was Lawrence's idea :)

[wence-]

To compute the eigenvalues more robustly, create submatrices by removing rows and columns associated with degrees of freedom. This trick was suggested (on a different example) by @jandrej , @wence- , @francesco-ballarin and Matt Knepley, among others.

Disclaimer I haven't run the demo to check that it runs because I don't know how to convert an *.py.rst into a *.py I can run.

Run make demos in the demos directory

[wence-]

I think this needs a bit of work to be correct in parallel.

[wence-]

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

[APaganini]

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

[APaganini]

Applied most of @wence- suggestions. Now it doesn't run because I don't know how to properly implement ISin = ISbd.complement(s, e)

[wence-]

I probably got the syntax wrong for getting the size of the matrix

[APaganini]

I corrected (s, e), _ = petsc_a.getSizes() into petsc_a.getOwnershipRange(). Now it runs. The result is slightly different from the original version. More specifially, the original output is

Leading eigenvalue is: (1.560678669632054e-13+0.10979244518390797j)

whereas after the modification from this pull request the output becomes

Leading eigenvalue is: (1.1395615354126143e-14+0.10979281297538487j)

The real part is different (but still very close to 0), whereas the imaginary parts differ by -3.67791477e-7. I assume this is fine, but my opinion on this is uninformed.

TODO:

• Resolve conversation about correct application of DirichletBC (see above)
• Fix the naming convention, which @wence- said they are rather off from normal Python. Suggestions are welcome.

[lrtfm]

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.

[dham]

If this is to be done it should be done in LinearEigenSolver. The right thing to do, though, is RestrictedFunctionSpace, which is hopefully coming.