[petsc-users] [petsc-maint] Create SIMPLE pc for Schur explicitly

[Ju Liu]

2018-03-06 14:49 GMT-08:00 Smith, Barry F. <bsmith at mcs.anl.gov>:

>
>
> > On Mar 6, 2018, at 4:37 PM, Ju Liu <liujuy at gmail.com> wrote:
> >
> >
> >
> > 2018-03-06 14:13 GMT-08:00 Smith, Barry F. <bsmith at mcs.anl.gov>:
> >
> >
> > > On Mar 6, 2018, at 4:09 PM, John <johnlucassaturday at gmail.com> wrote:
> > >
> > > Hi PETSc team:
> > >
> > > I am trying to create a SIMPLE type approximation for the Schur
> complement. It is a nonlinear code, so I have to update the Schur matrix
> once the tangent jacobian is updated. I am trying to figure out the best
> way of doing that.
> > >
> > > In terms of generating an algebraic form of Schur S, I am trying to do
> the following.
> > >
> > > // S = A10 x A01, I have performed a diagonal scaling on all matrices,
> so the diag(A00) = I.
> > > MatMatMult(A10, A01, MAT_INITIAL_MATRIX, S);
> > > MatScale(S, -1.0); // S = -1 x S
> > > MatAXPY(S, 1.0, A11, S, DIFFERENT_NONZERO_PATTERN); // S = A11 + S
> > >
> > > My questions are:
> > > (1) When I update the Schur complement, do I need to call
> > > MatDestroy(&S);
> > > MatMatMult(A10, A01, MAT_INITIAL_MATRIX, S);
> > > or shall I just call
> > > MatMatMult(A10, A01, MAT_REUSE_MATRIX, S);
> >
> >    The second form. No reason to delete the matrix and rebuild it each
> time (delete and rebuild will be a bit slower)
> >
> > >
> > > Basically, I am not sure how to properly use the flag MAT_REUSE_MATRIX.
> > >
> > > (2) When I update the Schur complement, the initial call of MatAXPY
> already incorporated the contribution of A00 in terms of nonzero patterns.
> So in my second call of MatAXPY, shall I use the flag
> SUBSET_NONZERO_PATTERN?
> >
> >    If it is a subset then you should use the subset flag, if they have
> same nonzero pattern you should use SAME_NONZERO_PATTERN (same nonzero is
> fastest)
> >
> > Thanks Barry! I guess I still need some clarifications on the second
> point.
> >
> > In the first assembly of S, I called MatAXPY to incorporate the entries
> of A11 into S (with differnet nonzero pattern for safety).
> >
> > Now, with new blocks of A, I update S by calling
> > MatMatMult(A10, A01, MAT_REUSE_MATRIX, S);
>
>    Ahh, you cannot do this if the MatAXPY introduces new nonzeros into the
> structure of S (the second MatMatMult() will likely crash or generate
> incorrect results since MAT_REUSE_MATRIX requires that each time the
> routine is called S has the exact same nonzero pattern.)
>
>    Now for most problems I would expect that A11 has a subset of the
> nonzeros of S, in that case you can use MatAXPY with subset each time and
> you can use MAT_REUSE_MATRIX. But you would need to verify for your problem
> that A11 has a subset of the nonzeros of S. This is the "ideal" case in
> terms of performance.
>
>    If A11 does not have a subset of the nonzeros of S then you (obviously)
> cannot use subset and you must destroy the S each time and create a new one
> with each MatMatMult().
>
> I see. Thanks a lot!

John


>    Barry
>
>
> > At this stage, what is the nonzero pattern of S? If S's pattern is
> determined by the most recent call of MatMatMult, I guess I will still use
> the different nonzero pattern. If the pattern of S is inherited from the
> first assembly, the A11's nonzero pattern is a subset of S's, so subset
> flag should be safe.
> >
> > Thanks!
> >
> > John
> >
> >
> >
> >
> >    Barry
> >
> > >
> > > Best regards,
> > >
> > > John
> >
> >
> >
> >
> > --
> > Ju Liu, Ph.D.
> > Postdoctoral Research Fellow
> > Stanford School of Medicine
> > http://ju-liu.github.io/
>
>


-- 
Ju Liu, Ph.D.
Postdoctoral Research Fellow
Stanford School of Medicine
http://ju-liu.github.io/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180306/a1e1a8c6/attachment.html>