VecGhost and PCFieldSplit

[Jed Brown]

Matthew Knepley wrote:
> On Wed, May 27, 2009 at 2:50 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
>> On May 27, 2009, at 2:23 PM, Jed Brown wrote:
>>
>>
>>> Gee, having vectors with a little extra space at the end is hard!
>>>
>>> Notation: J is the Jacobian, P is the associated preconditioning matrix
>>>
>>> PCFieldSplit gets the full row blocks (Afield) or the off-diagonal
>>> blocks (B,C) from J.  I implement these matrix-free and use VecGhost
>>> internally.  It also queries P for diagonal blocks (A,D).  I provide
>>> these as assembled AIJ and BAIJ matrices.  It then obtains work vectors
>>> by getting vecs from A and D which returns ordinary MPI vectors (no
>>> ghosts).
>>>
>>   You could provide the getvecs() method for the A and D matrices that you
>> provide to give you the ghosted vectors that you need?
> 
> 
> I agree with Barry. If you want to use ghosts, make everything use VecGhost.

Thanks guys, this occurred to me while riding home.  Hooray for dynamic
polymorphism.


I'm now using PCFieldSplit_Schur for non-Newtonian Stokes, with
everything completely matrix-free except the preconditioners for A and
S, so I'm finally phasing out my homegrown stuff.  I would like to add
an option -pc_fieldsplit_real_diagonal, to use diagonal blocks of J
preconditioned by diagonal blocks of Jp, instead of only using blocks of
Jp.  I implemented this using an array 'mat' obatained in the same way
as 'pmat', but you might have something else in mind.

I am running into an issue with order of events.  I want UserJacobian
to be able to determine what assembly is needed, but I can't do things
like call PCFieldSplitGetSubKSP because PCSetUp has not been called the
first time around.  The issue here is that I want to be able to
optionally replace the PC in kspschur with my own.  This is tricky
because we can't create the Schur complement matrix before PCSetUp, so
there won't be a subksp[0] and PCFieldSplitGetSubKSP fails.  Creating
the KSP before the Schur matrix exists creates an ownership issue that
would add complexity so I'm hesitant to do it unless it proves
necessary.


For the interface in these cases, I've been thinking about what
information to put in the preconditioning matrix (Jp and submatrices)
versus communicating directly to the PC.  My inclination is to put all
real data in Jp because it's messy for UserJacobian to talk directly to
the PC, and just move metadata (setting options, usually prior to
PCSetFromOptions) directly to the PC.  That is, I like to think of Jp,
less as a "matrix", and more as "data" that will be used by the
preconditioner.  It just so happens that an assembled matrix is a really
useful format (understood by almost all PCs) to store that "data" so
basing it on the Mat interface is appropriate.

To illustrate, consider the usual problem J = [A B; C D], S =
D-C*inv(A)*B.  I don't want PCFieldSplit to need to do anything special
based on how S gets preconditioned, and I don't want UserJacobian to
need to unwrap sub-KSPs from PCFieldSplit (creation issues aside, this
is a dependency that I think is avoidable).  But we have this Jp hanging
around and it's usual purpose is to route information from UserJacobian
through the KSP to the PC so maybe we can also use it to route
information through PCFieldSplit (via submatrices).  The query for Ap is
obvious, but what does a query for Dp mean?  The matrix that actually
sits there (often 0) is not actually useful.  One choice is to return
some other user-assembled matrix Sp (e.g. mass matrix for Stokes).  I
think this is a slight abuse of the interface since it specializes the
definition of Jp to a particular preconditioner.  In particular, it
won't work correctly with FieldSplit relaxation (but that would fail for
the class of problems we're discussing anyway).  Let's assume that
PCFieldSplit gets Sp somehow, perhaps by just by querying Dp.  Now what
if our preconditioner for S needs more information?  For example,
unscaled BFBt amounts to

  bfbt(S,Sp) = inv(Bp Cp) B A C inv(Bp Cp)

where the product Bp*Cp is actually formed.  Now "Sp" is no longer
needed as an assembled matrix, but assembled variants of Bp and Cp are
needed (B, A, and C are part of S).  An alternative is to recognize that
B*C corresponds to a Laplacian L, and actually assemble this as Lp.
Let's call this variation the least squares commutator.  It's normally
accompanied by a couple scaling vectors, but the unscaled version looks
like

  lsc(S,Sp) = ksp(L,Lp) B A C ksp(L,Lp)

How should we get L, Lp, and scaling vectors through the layers of
PCFieldSplit?  Another preconditioner for S is from Kay, Loghin, and
Wathen,

  klw(S,Sp) = inv(Mp) Ac ksp(L,Lp)

where Mp is a (lumped) pressure mass matrix and Ac is a representation
of A in the pressure space.  Now we need to get Mp, Ac, L, Lp through.
At this point, my preference is to compose this stuff with Sp.  This
avoids any need to unwrap internals of PCFieldSplit and PCFieldSplit
does not need to concern itself with what auxiliary information the
preconditioner for S might need.  It would also make it easy to write
alternate preconditioners for S, for example PCBFBt could be greatly
simplified, and lots of other variants could be created easily.


If we do indeed go this route, I see two remaining tricky bits.  First,
how to identify the various "Sp" showing up in a factorization?  When
pivoting is used, they are not uniquely identified by where they sit
(coordinates) in the factored matrix.

Second, how to avoid computing way more than needed in UserJacobian.
For example, if I want to support default, bfbt, lsc, and klw, I should
be able to provide Sp, Bp, Cp, L, Lp, Mp, Ac, and scaling, but not all
of these will be used in any particular run.  Before PCSetUp_FieldSplit
has been called, the submatrices don't even exist and none of the wiring
for the preconditioners has been set up.  One option is to assemble
everything the first time around, but this could could blow your memory.
Another is to assemble lazily (e.g. in MatGetSubMatrix) at least the
first time around (it's usually much faster to assemble several matrices
at once).  This goofs with profiling, but the first time is often
nonsense anyway.


Jed

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 260 bytes
Desc: OpenPGP digital signature
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090529/cdae0181/attachment.sig>