[petsc-dev] Note about CPR preconditioner in PCFIELDSPLIT man page?
Matthew Knepley
knepley at gmail.com
Sun Mar 3 10:34:36 CST 2019
On Sun, Mar 3, 2019 at 12:58 AM Jed Brown via petsc-dev <
petsc-dev at mcs.anl.gov> wrote:
> My take is that someone looking for CPR is more likely to end up on the
> PCFIELDSPLIT manual page than the PCGALERKIN page. The solver you have
> configured in your mail is not the CPR we have been asked about here.
> See Barry's message and example code below.
>
Thanks for retrieving this Jed. I am sure Richard and I both have the same
question. Perhaps I am being an idiot.
I am supposing that R is just a restriction to some subset of dofs, so its
just binary, so that R A P just selects that
submatrix. Thus my question is:
Why is this not equivalent to FS with the first block doing BAMG, and the
second identity, with additive combination?
Then you use COMPOSITE for the rest.
Thanks,
Matt
> ---------- Forwarded message ----------
> From: Barry Smith <bsmith at mcs.anl.gov>
> To: "Sébastien Loisel" <sloisel at gmail.com>
> Cc: Robert Annewandter <robert.annewandter at opengosim.com>, "Hammond,
> Glenn E" <gehammo at sandia.gov>, Jed Brown <jed at jedbrown.org>, Paolo Orsini
> <paolo.orsini at opengosim.com>, Matthew Knepley <knepley at gmail.com>
> Bcc:
> Date: Fri, 27 Jan 2017 17:50:39 -0600
> Subject: Re: [EXTERNAL] CPR preconditioning
>
> Ok, I have implemented the algorithm using PETSc PCCOMPOSITE and
> PCGALERKIN and get identical iterations as the code you sent. PCFIELDSPLIT
> is not intended for this type of solver composition. Here is the algorithm
> written in "two-step" form
>
> x_1/2 = P KSPSolve( R A P, using BoomerAMG) R b
> x = x_1/2 + PCApply( A, using Hypre PILUT preconditioner) ( b
> - A x_1/2)
>
> PCCOMPOSITE with a type of multiplicative handles the two steps and
> PCGALERKIN handles the P KSPSolve(R A P) R business.
>
> You will need to use the master version of PETSc because I had to add a
> feature to PCGALERKIN to allow the solver to be easily used for and A that
> changes values for later solvers.
>
>
> Here is the output from -ksp_view
>
> KSP Object: 1 MPI processes
> type: fgmres
> GMRES: restart=100, using Classical (unmodified) Gram-Schmidt
> Orthogonalization with no iterative refinement
> GMRES: happy breakdown tolerance 1e-30
> maximum iterations=10000, initial guess is zero
> tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
> right preconditioning
> using UNPRECONDITIONED norm type for convergence test
> PC Object: 1 MPI processes
> type: composite
> Composite PC type - MULTIPLICATIVE
> PCs on composite preconditioner follow
> ---------------------------------
> PC Object: (sub_0_) 1 MPI processes
> type: galerkin
> Galerkin PC
> KSP on Galerkin follow
> ---------------------------------
> KSP Object: (sub_0_galerkin_) 1 MPI processes
> type: richardson
> Richardson: damping factor=1.
> maximum iterations=10000, initial guess is zero
> tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
> left preconditioning
> using PRECONDITIONED norm type for convergence test
> PC Object: (sub_0_galerkin_) 1 MPI processes
> type: hypre
> HYPRE BoomerAMG preconditioning
> HYPRE BoomerAMG: Cycle type V
> HYPRE BoomerAMG: Maximum number of levels 25
> HYPRE BoomerAMG: Maximum number of iterations PER hypre call 1
> HYPRE BoomerAMG: Convergence tolerance PER hypre call 0.
> HYPRE BoomerAMG: Threshold for strong coupling 0.25
> HYPRE BoomerAMG: Interpolation truncation factor 0.
> HYPRE BoomerAMG: Interpolation: max elements per row 0
> HYPRE BoomerAMG: Number of levels of aggressive coarsening 0
> HYPRE BoomerAMG: Number of paths for aggressive coarsening 1
> HYPRE BoomerAMG: Maximum row sums 0.9
> HYPRE BoomerAMG: Sweeps down 1
> HYPRE BoomerAMG: Sweeps up 1
> HYPRE BoomerAMG: Sweeps on coarse 1
> HYPRE BoomerAMG: Relax down symmetric-SOR/Jacobi
> HYPRE BoomerAMG: Relax up symmetric-SOR/Jacobi
> HYPRE BoomerAMG: Relax on coarse Gaussian-elimination
> HYPRE BoomerAMG: Relax weight (all) 1.
> HYPRE BoomerAMG: Outer relax weight (all) 1.
> HYPRE BoomerAMG: Using CF-relaxation
> HYPRE BoomerAMG: Not using more complex smoothers.
> HYPRE BoomerAMG: Measure type local
> HYPRE BoomerAMG: Coarsen type Falgout
> HYPRE BoomerAMG: Interpolation type classical
> linear system matrix = precond matrix:
> Mat Object: 1 MPI processes
> type: seqaij
> rows=50, cols=50
> total: nonzeros=244, allocated nonzeros=244
> total number of mallocs used during MatSetValues calls =0
> not using I-node routines
> linear system matrix = precond matrix:
> Mat Object: 1 MPI processes
> type: seqaij
> rows=100, cols=100, bs=2
> total: nonzeros=976, allocated nonzeros=100000
> total number of mallocs used during MatSetValues calls =0
> using I-node routines: found 50 nodes, limit used is 5
> PC Object: (sub_1_) 1 MPI processes
> type: hypre
> HYPRE Pilut preconditioning
> HYPRE Pilut: maximum number of iterations 1
> HYPRE Pilut: drop tolerance 0.1
> HYPRE Pilut: default factor row size
> linear system matrix = precond matrix:
> Mat Object: 1 MPI processes
> type: seqaij
> rows=100, cols=100, bs=2
> total: nonzeros=976, allocated nonzeros=100000
> total number of mallocs used during MatSetValues calls =0
> using I-node routines: found 50 nodes, limit used is 5
>
> Note that one can change the solvers on the two "stages" and all their
> options such as tolerances etc using the options database a proper prefixes
> which you can find in the output above. For example to use PETSc's ILU
> instead of hypre's just run with -sub_1_pc_type ilu or to use a direct
> solver instead of boomer -sub_0_galerkin_pc_type lu
>
> I've attached a new version of testmain2.c that runs your solver and also
> my version.
>
> Given the "unique" nature of your R = [ I I ... I] and P = [0 I 0 0 ... 0]
> I am not sure that it makes sense to include this preconditioner directly
> in PETSc as a new PC type; so if you are serious about using it you can
> take what I am sending back and modify it for your needs. As a numerical
> analyst who works on linear solvers I am also not convinced that this is
> likely to be a particular good preconditioner
>
> Let me know if you have any questions,
>
> Barry
>
>
>
>
> > On Jan 23, 2017, at 8:25 AM, Sébastien Loisel <sloisel at gmail.com> wrote:
> >
> > Hi Barry,
> >
> > Thanks for your email. Thanks for pointing out my PetSC sillyness. I
> think what happened is I played with matrices as well as with
> preconditioners so I initially implemented it as MATSHELL and at the end
> wrapped it in a PCSHELL. :)
> >
> > On Mon, Jan 23, 2017 at 2:26 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> > I've started to look at this; it is a little weird using MATSHELL within
> your PCSHELL, seems unnecessary. Not necessarily wrong, just strange. I'll
> continue to try to understand it.
> >
> > Barry
> >
> >
> >
> > > On Jan 20, 2017, at 5:48 PM, Sébastien Loisel <sloisel at gmail.com>
> wrote:
> > >
> > > OK I'm attaching the prototype.
> > >
> > > It won't be 100% plug-and-play for you because in principle T1 and T2
> are built on top of "sub-preconditioners" (AMG for T1 and BILU for T2) and
> judging from the PetSC architecture elsewhere I must assume you're going to
> want to expose those in a rational way. At present, I've hard-coded some
> sub-preconditioners, and in particular for T2 I had to resort to PILUT
> because I didn't have a BILU(0) handy.
> > >
> > > Also, I broke the PetSC law and my functions don't return integers, so
> I also assume you're going to want to tweak that... Sorry!
> > >
> > > On Fri, Jan 20, 2017 at 11:29 PM, Barry Smith <bsmith at mcs.anl.gov>
> wrote:
> > >
> > > Sure, email the PCSHELL version, best case scenario I only need
> change out the PCSHELL and it takes me 5 minutes :-)
> > >
> > >
> > > > On Jan 20, 2017, at 5:07 PM, Sébastien Loisel <sloisel at gmail.com>
> wrote:
> > > >
> > > > Hi all,
> > > >
> > > > Thanks for your emails. I'm willing to help in whatever way. We have
> a "PCSHELL" prototype we can provide on request, although PetSC experts can
> no doubt do a better job than I did.
> > > >
> > > > Thanks,
> > > >
> > > > On Fri, Jan 20, 2017 at 9:50 PM, Robert Annewandter <
> robert.annewandter at opengosim.com> wrote:
> > > > Indeed that would be very helpful! We're very happy to support to
> test things out, provide feedback etc
> > > >
> > > > Many thanks!
> > > > Robert
> > > >
> > > >
> > > >
> > > > On 20/01/17 21:22, Hammond, Glenn E wrote:
> > > >> That sounds great. I do know that there is also much interest
> state-side in CPR preconditioning within PFLOTRAN. I have a Sandia
> colleague in Carlsbad, NM who has been asking about it. I am sure that
> Sebastien and/or Robert will help out in any way possible.
> > > >>
> > > >> Thanks,
> > > >>
> > > >> Glenn
> > > >>
> > > >>
> > > >>> -----Original Message-----
> > > >>> From: Barry Smith [
> > > >>> mailto:bsmith at mcs.anl.gov
> > > >>> ]
> > > >>> Sent: Friday, January 20, 2017 12:58 PM
> > > >>> To: Hammond, Glenn E
> > > >>> <gehammo at sandia.gov>
> > > >>>
> > > >>> Cc: Sébastien Loisel
> > > >>> <sloisel at gmail.com>
> > > >>> ; Robert Annewandter
> > > >>>
> > > >>> <robert.annewandter at opengosim.com>; Jed Brown <jed at jedbrown.org>
> > > >>> ;
> > > >>> Paolo Orsini
> > > >>> <paolo.orsini at opengosim.com>
> > > >>> ; Matthew Knepley
> > > >>>
> > > >>> <knepley at gmail.com>
> > > >>>
> > > >>> Subject: Re: [EXTERNAL] CPR preconditioning
> > > >>>
> > > >>>
> > > >>> Glenn,
> > > >>>
> > > >>> Sorry about the delay in processing this, too much going on ...
> > > >>>
> > > >>> I think the best thing is for us (the PETSc developers) to
> implement a CPR
> > > >>> preconditioner "directly" as its own PC and then have you guys try
> it out. I am
> > > >>> planning to do this.
> > > >>>
> > > >>> Barry
> > > >>>
> > > >>>
> > > >>>> On Jan 20, 2017, at 2:50 PM, Hammond, Glenn E <gehammo at sandia.gov
> >
> > > >>> wrote:
> > > >>>
> > > >>>> Barry, Jed or Matt,
> > > >>>>
> > > >>>> Do you have any suggestions for how to address the limitations of
> > > >>>>
> > > >>> PetscFieldSplit() discussed below. Will they need to manipulate
> the matrices
> > > >>> manually?
> > > >>>
> > > >>>> Thanks,
> > > >>>>
> > > >>>> Glenn
> > > >>>>
> > > >>>> From: Sébastien Loisel [
> > > >>>> mailto:sloisel at gmail.com
> > > >>>> ]
> > > >>>> Sent: Wednesday, January 11, 2017 3:33 AM
> > > >>>> To: Barry Smith
> > > >>>> <bsmith at mcs.anl.gov>
> > > >>>> ; Robert Annewandter
> > > >>>>
> > > >>>> <robert.annewandter at opengosim.com>
> > > >>>> ; Hammond, Glenn E
> > > >>>>
> > > >>>> <gehammo at sandia.gov>; Jed Brown <jed at jedbrown.org>
> > > >>>> ; Paolo Orsini
> > > >>>>
> > > >>>> <paolo.orsini at opengosim.com>
> > > >>>>
> > > >>>> Subject: [EXTERNAL] CPR preconditioning
> > > >>>>
> > > >>>> Dear Friends,
> > > >>>>
> > > >>>> Paolo has asked me to write this email to clarify issues
> surrounding
> > > >>>> the CPR preconditioner that is widely used in multiphase flow. I
> know
> > > >>>> Barry from a long time ago but we only met once when I was a PhD
> > > >>>> student so I would be shocked if he remembered me. :)
> > > >>>>
> > > >>>> I'm a math assistant professor and one of my areas of
> specialization is linear
> > > >>>>
> > > >>> algebra and preconditioning.
> > > >>>
> > > >>>> The main issue that is useful to clarify is the following. There
> was a proposal
> > > >>>>
> > > >>> to use PetSC's PETSCFIELDSPLIT in conjunction with PCCOMPOSITE in
> order to
> > > >>> implement CPR preconditioning. Although this is morally the right
> idea, this
> > > >>> seems to be currently impossible because PETSCFIELDSPLIT lacks the
> > > >>> capabilities it would need to implement the T1 preconditioner.
> This is due to
> > > >>> limitations in the API exposed by PETSCFIELDSPLIT (and no doubt
> limitations
> > > >>> in the underlying implementation).
> > > >>>
> > > >>>> In order to be as clear as possible, please allow me to describe
> > > >>>>
> > > >>> unambiguously the first of the two parts of the CPR preconditioner
> using
> > > >>> some MATLAB snippets. Let I denote the N by N identity, and Z the
> N by N
> > > >>> zero matrix. Put WT = [I I I] and C = [Z;Z;I]. The pressure matrix
> is Ap =
> > > >>> WT*A*C, and the T1 preconditioner is C*inv(Ap)*WT, where inv(Ap)
> is to be
> > > >>> implemented with AMG.
> > > >>>
> > > >>>> This T1 preconditioner is the one that would have to be
> implemented by
> > > >>>>
> > > >>> PETSCFIELDSPLIT. The following limitations in PETSCFIELDSPLIT
> prevents one
> > > >>> to implement T1:
> > > >>>
> > > >>>> • One must select the "pressure" by specifying an IS object to
> > > >>>>
> > > >>> PCFieldSplitSetIS(). However, since WT = [I I I], the pressure is
> obtained by
> > > >>> summing over the three fields. As far as I can tell, an IS object
> does not allow
> > > >>> one to sum over several entries to obtain the pressure field.
> > > >>>
> > > >>>> • The pressure matrix is Ap = WT*A*C; note that the matrix WT
> on
> > > >>>>
> > > >>> the left is different from the matrix C on the right. However,
> PCFIELDSPLIT
> > > >>> has no notion of a "left-IS" vs "right-IS"; morally, only diagonal
> blocks of A can
> > > >>> be used by PCFIELDSPLIT.
> > > >>>
> > > >>>> • PCFIELDSPLIT offers a range of hard-coded block structures
> for the
> > > >>>>
> > > >>> final preconditioner, but the choice T1 = C*inv(Ap)*WT is not one
> of these
> > > >>> choices. Indeed, this first stage CPR preconditioner T1 is
> *singular*, but there
> > > >>> is no obvious way for PCFIELDSPLIT to produce a singular
> preconditioner.
> > > >>>
> > > >>>> Note that the documentation for PETSCFIELDSPLIT says that "The
> > > >>>>
> > > >>> Constrained Pressure Preconditioner (CPR) does not appear to be
> currently
> > > >>> implementable directly with PCFIELDSPLIT". Unless there are very
> significant
> > > >>> capabilities that are not documented, I don't see how CPR can be
> > > >>> implemented with PETSCFIELDSPLIT.
> > > >>>
> > > >>>> Elsewhere, someone proposed putting the two preconditioners T1
> and T2
> > > >>>> on each side of A, e.g. T1*A*T2. That is a very bad idea because
> T1 is
> > > >>>> singular and hence T1*A*T2 is also singular. The correct CPR
> > > >>>> preconditioner is nonsingular despite the fact that T1 is
> singular,
> > > >>>> and MCPR is given by the formula MCPR = T2*(I-A*T1)+T1, where T2 =
> > > >>>> BILU(0) of A. (There is a proof, due to Felix Kwok, that BILU(0)
> works
> > > >>>> even though ILU(0) craps out on vanishing diagonal entries.)
> > > >>>>
> > > >>>> I'm also attaching a sample MATLAB code that runs the CPR
> preconditioner
> > > >>>>
> > > >>> on some fabricated random matrix A. I emphasize that this is not a
> realistic
> > > >>> matrix, but it still illustrates that the algorithm works, and
> that MCPR is better
> > > >>> than T2 alone. Note that T1 cannot be used alone since it is
> singular. Further
> > > >>> gains are expected when the Robert's realistic code with correct
> physics will
> > > >>> come online.
> > > >>>
> > > >>>> <image003.jpg>
> > > >>>> I hope this clarifies some things.
> > > >>>>
> > > >>>>
> > > >>>> Sébastien Loisel
> > > >>>> Assistant Professor
> > > >>>> Department of Mathematics, Heriot-Watt University Riccarton, EH14
> 4AS,
> > > >>>> United Kingdom
> > > >>>> web:
> > > >>>> http://www.ma.hw.ac.uk/~loisel/
> > > >>>>
> > > >>>> email: S.Loisel at
> > > >>>> hw.ac.uk
> > > >>>>
> > > >>>> phone:
> > > >>>> +44 131 451 3234
> > > >>>>
> > > >>>> fax:
> > > >>>> +44 131 451 3249
> > > >
> > > >
> > > >
> > > >
> > > > --
> > > > Sébastien Loisel
> > > > Assistant Professor
> > > > Department of Mathematics, Heriot-Watt University
> > > > Riccarton, EH14 4AS, United Kingdom
> > > > web: http://www.ma.hw.ac.uk/~loisel/
> > > > email: S.Loisel at hw.ac.uk
> > > > phone: +44 131 451 3234
> > > > fax: +44 131 451 3249
> > > >
> > >
> > >
> > >
> > >
> > > --
> > > Sébastien Loisel
> > > Assistant Professor
> > > Department of Mathematics, Heriot-Watt University
> > > Riccarton, EH14 4AS, United Kingdom
> > > web: http://www.ma.hw.ac.uk/~loisel/
> > > email: S.Loisel at hw.ac.uk
> > > phone: +44 131 451 3234
> > > fax: +44 131 451 3249
> > >
> > > <petsc.zip>
> >
> >
> >
> >
> > --
> > Sébastien Loisel
> > Assistant Professor
> > Department of Mathematics, Heriot-Watt University
> > Riccarton, EH14 4AS, United Kingdom
> > web: http://www.ma.hw.ac.uk/~loisel/
> > email: S.Loisel at hw.ac.uk
> > phone: +44 131 451 3234
> > fax: +44 131 451 3249
> >
>
>
--
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
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20190303/5dc913cb/attachment-0001.html>
More information about the petsc-dev
mailing list