[petsc-users] Solving eigenproblems in the form of Stokes equations

[Julian Andrej]

On Sun, Jul 12, 2015 at 1:07 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Fri, Jul 10, 2015 at 4:28 PM, Julian Andrej <juan at tf.uni-kiel.de> wrote:
>>
>> Hi,
>>
>> i'm trying to solve a generalized eigenvalue problem which is the
>> stokes equations discretized by finite elements (with fenics) and
>> producing a banded matrix (reordered structure) of the well known
>> block form
>>
>> [N    Q]
>> [QT  0]
>>
>> The domain is the unit square with dirichlet boundary condition
>> evaluating to zero at all boundary nodes.
>>
>> A is the block matrix in banded reordered structure and M is the mass
>> matrix.
>
>
> You have not written an equation here. Are you solving
>
>   Q^T N^{-1} Q x = \lambda M x

I don't have / want access to the underlying blocks of the matrix. So
in fact i give the solver

A = [N  Q]      M = [M 0]
      [QT 0]             [0  0]

With a little more reading of the documentation i was able to figure
out that the EPS object also takes regular KSP objects. So i got the
following settings working.

E = SLEPc.EPS().create()
st = E.getST()
st.setType('sinvert')
kspE = st.getKSP()
kspE.setType('gmres')
pc = kspE.getPC()
pc.setType('lu')
pc.setFactorSolverPackage('mumps')
E.setOperators(A, M)

These options give me an eigenvalue bundle at 1.0+j0 (which represent
the DirichletBC) and the correct eigenvalues following.
I'd like to understand the background a little more, how this produces
a solvable eigenproblem.

>
> If so, the Schur complement can be rank deficient by 1 if you have Neumann
> conditions on the pressure. I don't know how SLEPc handles this.
>
>    Matt
>
>>
>> I'm using slepc4py for the eigenvalue calculation.
>>
>> E = SLEPc.EPS().create()
>> E.setOperators(A, M)
>> E.setDimensions(NEV, PETSc.DECIDE)
>> E.setFromOptions()
>>
>> I always get the error
>> [0] Zero pivot row 1 value 0 tolerance 2.22045e-14
>>
>> I cannot find a combination which solves the eigenvalues for this problem.
>>
>> The system itself solves fine with a KSP solver object from PETSc
>> using tfmqr with the icc PC.
>>
>> I can assemble the matrix with a penalty term for the pressure and
>> calculate the eigenvalues with a direct solver, but i try to avoid
>> that.
>
>
>
>
> --
> What most experimenters take for granted before they begin their experiments
> is infinitely more interesting than any results to which their experiments
> lead.
> -- Norbert Wiener