<div dir="ltr"><div dir="ltr">On Sun, Mar 3, 2019 at 12:58 AM Jed Brown via petsc-dev <<a href="mailto:petsc-dev@mcs.anl.gov">petsc-dev@mcs.anl.gov</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">My take is that someone looking for CPR is more likely to end up on the<br>
PCFIELDSPLIT manual page than the PCGALERKIN page. The solver you have<br>
configured in your mail is not the CPR we have been asked about here.<br>
See Barry's message and example code below.<br></blockquote><div><br></div><div>Thanks for retrieving this Jed. I am sure Richard and I both have the same question. Perhaps I am being an idiot.</div><div>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</div><div>submatrix. Thus my question is:</div><div><br></div><div> Why is this not equivalent to FS with the first block doing BAMG, and the second identity, with additive combination?</div><div><br></div><div>Then you use COMPOSITE for the rest.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">---------- Forwarded message ----------<br>From: Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>><br>To: "Sébastien Loisel" <<a href="mailto:sloisel@gmail.com" target="_blank">sloisel@gmail.com</a>><br>Cc: Robert Annewandter <<a href="mailto:robert.annewandter@opengosim.com" target="_blank">robert.annewandter@opengosim.com</a>>, "Hammond, Glenn E" <<a href="mailto:gehammo@sandia.gov" target="_blank">gehammo@sandia.gov</a>>, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>>, Paolo Orsini <<a href="mailto:paolo.orsini@opengosim.com" target="_blank">paolo.orsini@opengosim.com</a>>, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br>Bcc: <br>Date: Fri, 27 Jan 2017 17:50:39 -0600<br>Subject: Re: [EXTERNAL] CPR preconditioning<br><br>
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<br>
<br>
x_1/2 = P KSPSolve( R A P, using BoomerAMG) R b<br>
x = x_1/2 + PCApply( A, using Hypre PILUT preconditioner) ( b - A x_1/2)<br>
<br>
PCCOMPOSITE with a type of multiplicative handles the two steps and PCGALERKIN handles the P KSPSolve(R A P) R business.<br>
<br>
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.<br>
<br>
<br>
Here is the output from -ksp_view <br>
<br>
KSP Object: 1 MPI processes<br>
type: fgmres<br>
GMRES: restart=100, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement<br>
GMRES: happy breakdown tolerance 1e-30<br>
maximum iterations=10000, initial guess is zero<br>
tolerances: relative=1e-05, absolute=1e-50, divergence=10000.<br>
right preconditioning<br>
using UNPRECONDITIONED norm type for convergence test<br>
PC Object: 1 MPI processes<br>
type: composite<br>
Composite PC type - MULTIPLICATIVE<br>
PCs on composite preconditioner follow<br>
---------------------------------<br>
PC Object: (sub_0_) 1 MPI processes<br>
type: galerkin<br>
Galerkin PC<br>
KSP on Galerkin follow<br>
---------------------------------<br>
KSP Object: (sub_0_galerkin_) 1 MPI processes<br>
type: richardson<br>
Richardson: damping factor=1.<br>
maximum iterations=10000, initial guess is zero<br>
tolerances: relative=1e-05, absolute=1e-50, divergence=10000.<br>
left preconditioning<br>
using PRECONDITIONED norm type for convergence test<br>
PC Object: (sub_0_galerkin_) 1 MPI processes<br>
type: hypre<br>
HYPRE BoomerAMG preconditioning<br>
HYPRE BoomerAMG: Cycle type V<br>
HYPRE BoomerAMG: Maximum number of levels 25<br>
HYPRE BoomerAMG: Maximum number of iterations PER hypre call 1<br>
HYPRE BoomerAMG: Convergence tolerance PER hypre call 0.<br>
HYPRE BoomerAMG: Threshold for strong coupling 0.25<br>
HYPRE BoomerAMG: Interpolation truncation factor 0.<br>
HYPRE BoomerAMG: Interpolation: max elements per row 0<br>
HYPRE BoomerAMG: Number of levels of aggressive coarsening 0<br>
HYPRE BoomerAMG: Number of paths for aggressive coarsening 1<br>
HYPRE BoomerAMG: Maximum row sums 0.9<br>
HYPRE BoomerAMG: Sweeps down 1<br>
HYPRE BoomerAMG: Sweeps up 1<br>
HYPRE BoomerAMG: Sweeps on coarse 1<br>
HYPRE BoomerAMG: Relax down symmetric-SOR/Jacobi<br>
HYPRE BoomerAMG: Relax up symmetric-SOR/Jacobi<br>
HYPRE BoomerAMG: Relax on coarse Gaussian-elimination<br>
HYPRE BoomerAMG: Relax weight (all) 1.<br>
HYPRE BoomerAMG: Outer relax weight (all) 1.<br>
HYPRE BoomerAMG: Using CF-relaxation<br>
HYPRE BoomerAMG: Not using more complex smoothers.<br>
HYPRE BoomerAMG: Measure type local<br>
HYPRE BoomerAMG: Coarsen type Falgout<br>
HYPRE BoomerAMG: Interpolation type classical<br>
linear system matrix = precond matrix:<br>
Mat Object: 1 MPI processes<br>
type: seqaij<br>
rows=50, cols=50<br>
total: nonzeros=244, allocated nonzeros=244<br>
total number of mallocs used during MatSetValues calls =0<br>
not using I-node routines<br>
linear system matrix = precond matrix:<br>
Mat Object: 1 MPI processes<br>
type: seqaij<br>
rows=100, cols=100, bs=2<br>
total: nonzeros=976, allocated nonzeros=100000<br>
total number of mallocs used during MatSetValues calls =0<br>
using I-node routines: found 50 nodes, limit used is 5<br>
PC Object: (sub_1_) 1 MPI processes<br>
type: hypre<br>
HYPRE Pilut preconditioning<br>
HYPRE Pilut: maximum number of iterations 1<br>
HYPRE Pilut: drop tolerance 0.1<br>
HYPRE Pilut: default factor row size <br>
linear system matrix = precond matrix:<br>
Mat Object: 1 MPI processes<br>
type: seqaij<br>
rows=100, cols=100, bs=2<br>
total: nonzeros=976, allocated nonzeros=100000<br>
total number of mallocs used during MatSetValues calls =0<br>
using I-node routines: found 50 nodes, limit used is 5<br>
<br>
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 <br>
<br>
I've attached a new version of testmain2.c that runs your solver and also my version.<br>
<br>
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<br>
<br>
Let me know if you have any questions,<br>
<br>
Barry<br>
<br>
<br>
<br>
<br>
> On Jan 23, 2017, at 8:25 AM, Sébastien Loisel <<a href="mailto:sloisel@gmail.com" target="_blank">sloisel@gmail.com</a>> wrote:<br>
> <br>
> Hi Barry,<br>
> <br>
> 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. :)<br>
> <br>
> On Mon, Jan 23, 2017 at 2:26 AM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
> 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.<br>
> <br>
> Barry<br>
> <br>
> <br>
> <br>
> > On Jan 20, 2017, at 5:48 PM, Sébastien Loisel <<a href="mailto:sloisel@gmail.com" target="_blank">sloisel@gmail.com</a>> wrote:<br>
> ><br>
> > OK I'm attaching the prototype.<br>
> ><br>
> > 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.<br>
> ><br>
> > 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!<br>
> ><br>
> > On Fri, Jan 20, 2017 at 11:29 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
> ><br>
> > Sure, email the PCSHELL version, best case scenario I only need change out the PCSHELL and it takes me 5 minutes :-)<br>
> ><br>
> ><br>
> > > On Jan 20, 2017, at 5:07 PM, Sébastien Loisel <<a href="mailto:sloisel@gmail.com" target="_blank">sloisel@gmail.com</a>> wrote:<br>
> > ><br>
> > > Hi all,<br>
> > ><br>
> > > 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.<br>
> > ><br>
> > > Thanks,<br>
> > ><br>
> > > On Fri, Jan 20, 2017 at 9:50 PM, Robert Annewandter <<a href="mailto:robert.annewandter@opengosim.com" target="_blank">robert.annewandter@opengosim.com</a>> wrote:<br>
> > > Indeed that would be very helpful! We're very happy to support to test things out, provide feedback etc<br>
> > ><br>
> > > Many thanks!<br>
> > > Robert<br>
> > ><br>
> > ><br>
> > ><br>
> > > On 20/01/17 21:22, Hammond, Glenn E wrote:<br>
> > >> 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.<br>
> > >><br>
> > >> Thanks,<br>
> > >><br>
> > >> Glenn<br>
> > >><br>
> > >><br>
> > >>> -----Original Message-----<br>
> > >>> From: Barry Smith [<br>
> > >>> mailto:<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a><br>
> > >>> ]<br>
> > >>> Sent: Friday, January 20, 2017 12:58 PM<br>
> > >>> To: Hammond, Glenn E<br>
> > >>> <<a href="mailto:gehammo@sandia.gov" target="_blank">gehammo@sandia.gov</a>><br>
> > >>><br>
> > >>> Cc: Sébastien Loisel<br>
> > >>> <<a href="mailto:sloisel@gmail.com" target="_blank">sloisel@gmail.com</a>><br>
> > >>> ; Robert Annewandter<br>
> > >>><br>
> > >>> <<a href="mailto:robert.annewandter@opengosim.com" target="_blank">robert.annewandter@opengosim.com</a>>; Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>><br>
> > >>> ;<br>
> > >>> Paolo Orsini<br>
> > >>> <<a href="mailto:paolo.orsini@opengosim.com" target="_blank">paolo.orsini@opengosim.com</a>><br>
> > >>> ; Matthew Knepley<br>
> > >>><br>
> > >>> <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br>
> > >>><br>
> > >>> Subject: Re: [EXTERNAL] CPR preconditioning<br>
> > >>><br>
> > >>><br>
> > >>> Glenn,<br>
> > >>><br>
> > >>> Sorry about the delay in processing this, too much going on ...<br>
> > >>><br>
> > >>> I think the best thing is for us (the PETSc developers) to implement a CPR<br>
> > >>> preconditioner "directly" as its own PC and then have you guys try it out. I am<br>
> > >>> planning to do this.<br>
> > >>><br>
> > >>> Barry<br>
> > >>><br>
> > >>><br>
> > >>>> On Jan 20, 2017, at 2:50 PM, Hammond, Glenn E <<a href="mailto:gehammo@sandia.gov" target="_blank">gehammo@sandia.gov</a>><br>
> > >>> wrote:<br>
> > >>><br>
> > >>>> Barry, Jed or Matt,<br>
> > >>>><br>
> > >>>> Do you have any suggestions for how to address the limitations of<br>
> > >>>><br>
> > >>> PetscFieldSplit() discussed below. Will they need to manipulate the matrices<br>
> > >>> manually?<br>
> > >>><br>
> > >>>> Thanks,<br>
> > >>>><br>
> > >>>> Glenn<br>
> > >>>><br>
> > >>>> From: Sébastien Loisel [<br>
> > >>>> mailto:<a href="mailto:sloisel@gmail.com" target="_blank">sloisel@gmail.com</a><br>
> > >>>> ]<br>
> > >>>> Sent: Wednesday, January 11, 2017 3:33 AM<br>
> > >>>> To: Barry Smith<br>
> > >>>> <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>><br>
> > >>>> ; Robert Annewandter<br>
> > >>>><br>
> > >>>> <<a href="mailto:robert.annewandter@opengosim.com" target="_blank">robert.annewandter@opengosim.com</a>><br>
> > >>>> ; Hammond, Glenn E<br>
> > >>>><br>
> > >>>> <<a href="mailto:gehammo@sandia.gov" target="_blank">gehammo@sandia.gov</a>>; Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>><br>
> > >>>> ; Paolo Orsini<br>
> > >>>><br>
> > >>>> <<a href="mailto:paolo.orsini@opengosim.com" target="_blank">paolo.orsini@opengosim.com</a>><br>
> > >>>><br>
> > >>>> Subject: [EXTERNAL] CPR preconditioning<br>
> > >>>><br>
> > >>>> Dear Friends,<br>
> > >>>><br>
> > >>>> Paolo has asked me to write this email to clarify issues surrounding<br>
> > >>>> the CPR preconditioner that is widely used in multiphase flow. I know<br>
> > >>>> Barry from a long time ago but we only met once when I was a PhD<br>
> > >>>> student so I would be shocked if he remembered me. :)<br>
> > >>>><br>
> > >>>> I'm a math assistant professor and one of my areas of specialization is linear<br>
> > >>>><br>
> > >>> algebra and preconditioning.<br>
> > >>><br>
> > >>>> The main issue that is useful to clarify is the following. There was a proposal<br>
> > >>>><br>
> > >>> to use PetSC's PETSCFIELDSPLIT in conjunction with PCCOMPOSITE in order to<br>
> > >>> implement CPR preconditioning. Although this is morally the right idea, this<br>
> > >>> seems to be currently impossible because PETSCFIELDSPLIT lacks the<br>
> > >>> capabilities it would need to implement the T1 preconditioner. This is due to<br>
> > >>> limitations in the API exposed by PETSCFIELDSPLIT (and no doubt limitations<br>
> > >>> in the underlying implementation).<br>
> > >>><br>
> > >>>> In order to be as clear as possible, please allow me to describe<br>
> > >>>><br>
> > >>> unambiguously the first of the two parts of the CPR preconditioner using<br>
> > >>> some MATLAB snippets. Let I denote the N by N identity, and Z the N by N<br>
> > >>> zero matrix. Put WT = [I I I] and C = [Z;Z;I]. The pressure matrix is Ap =<br>
> > >>> WT*A*C, and the T1 preconditioner is C*inv(Ap)*WT, where inv(Ap) is to be<br>
> > >>> implemented with AMG.<br>
> > >>><br>
> > >>>> This T1 preconditioner is the one that would have to be implemented by<br>
> > >>>><br>
> > >>> PETSCFIELDSPLIT. The following limitations in PETSCFIELDSPLIT prevents one<br>
> > >>> to implement T1:<br>
> > >>><br>
> > >>>> • One must select the "pressure" by specifying an IS object to<br>
> > >>>><br>
> > >>> PCFieldSplitSetIS(). However, since WT = [I I I], the pressure is obtained by<br>
> > >>> summing over the three fields. As far as I can tell, an IS object does not allow<br>
> > >>> one to sum over several entries to obtain the pressure field.<br>
> > >>><br>
> > >>>> • The pressure matrix is Ap = WT*A*C; note that the matrix WT on<br>
> > >>>><br>
> > >>> the left is different from the matrix C on the right. However, PCFIELDSPLIT<br>
> > >>> has no notion of a "left-IS" vs "right-IS"; morally, only diagonal blocks of A can<br>
> > >>> be used by PCFIELDSPLIT.<br>
> > >>><br>
> > >>>> • PCFIELDSPLIT offers a range of hard-coded block structures for the<br>
> > >>>><br>
> > >>> final preconditioner, but the choice T1 = C*inv(Ap)*WT is not one of these<br>
> > >>> choices. Indeed, this first stage CPR preconditioner T1 is *singular*, but there<br>
> > >>> is no obvious way for PCFIELDSPLIT to produce a singular preconditioner.<br>
> > >>><br>
> > >>>> Note that the documentation for PETSCFIELDSPLIT says that "The<br>
> > >>>><br>
> > >>> Constrained Pressure Preconditioner (CPR) does not appear to be currently<br>
> > >>> implementable directly with PCFIELDSPLIT". Unless there are very significant<br>
> > >>> capabilities that are not documented, I don't see how CPR can be<br>
> > >>> implemented with PETSCFIELDSPLIT.<br>
> > >>><br>
> > >>>> Elsewhere, someone proposed putting the two preconditioners T1 and T2<br>
> > >>>> on each side of A, e.g. T1*A*T2. That is a very bad idea because T1 is<br>
> > >>>> singular and hence T1*A*T2 is also singular. The correct CPR<br>
> > >>>> preconditioner is nonsingular despite the fact that T1 is singular,<br>
> > >>>> and MCPR is given by the formula MCPR = T2*(I-A*T1)+T1, where T2 =<br>
> > >>>> BILU(0) of A. (There is a proof, due to Felix Kwok, that BILU(0) works<br>
> > >>>> even though ILU(0) craps out on vanishing diagonal entries.)<br>
> > >>>><br>
> > >>>> I'm also attaching a sample MATLAB code that runs the CPR preconditioner<br>
> > >>>><br>
> > >>> on some fabricated random matrix A. I emphasize that this is not a realistic<br>
> > >>> matrix, but it still illustrates that the algorithm works, and that MCPR is better<br>
> > >>> than T2 alone. Note that T1 cannot be used alone since it is singular. Further<br>
> > >>> gains are expected when the Robert's realistic code with correct physics will<br>
> > >>> come online.<br>
> > >>><br>
> > >>>> <image003.jpg><br>
> > >>>> I hope this clarifies some things.<br>
> > >>>><br>
> > >>>><br>
> > >>>> Sébastien Loisel<br>
> > >>>> Assistant Professor<br>
> > >>>> Department of Mathematics, Heriot-Watt University Riccarton, EH14 4AS,<br>
> > >>>> United Kingdom<br>
> > >>>> web:<br>
> > >>>> <a href="http://www.ma.hw.ac.uk/~loisel/" rel="noreferrer" target="_blank">http://www.ma.hw.ac.uk/~loisel/</a><br>
> > >>>><br>
> > >>>> email: S.Loisel at<br>
> > >>>> <a href="http://hw.ac.uk" rel="noreferrer" target="_blank">hw.ac.uk</a><br>
> > >>>><br>
> > >>>> phone:<br>
> > >>>> +44 131 451 3234<br>
> > >>>><br>
> > >>>> fax:<br>
> > >>>> +44 131 451 3249<br>
> > ><br>
> > ><br>
> > ><br>
> > ><br>
> > > --<br>
> > > Sébastien Loisel<br>
> > > Assistant Professor<br>
> > > Department of Mathematics, Heriot-Watt University<br>
> > > Riccarton, EH14 4AS, United Kingdom<br>
> > > web: <a href="http://www.ma.hw.ac.uk/~loisel/" rel="noreferrer" target="_blank">http://www.ma.hw.ac.uk/~loisel/</a><br>
> > > email: S.Loisel at <a href="http://hw.ac.uk" rel="noreferrer" target="_blank">hw.ac.uk</a><br>
> > > phone: +44 131 451 3234<br>
> > > fax: +44 131 451 3249<br>
> > ><br>
> ><br>
> ><br>
> ><br>
> ><br>
> > --<br>
> > Sébastien Loisel<br>
> > Assistant Professor<br>
> > Department of Mathematics, Heriot-Watt University<br>
> > Riccarton, EH14 4AS, United Kingdom<br>
> > web: <a href="http://www.ma.hw.ac.uk/~loisel/" rel="noreferrer" target="_blank">http://www.ma.hw.ac.uk/~loisel/</a><br>
> > email: S.Loisel at <a href="http://hw.ac.uk" rel="noreferrer" target="_blank">hw.ac.uk</a><br>
> > phone: +44 131 451 3234<br>
> > fax: +44 131 451 3249<br>
> ><br>
> > <petsc.zip><br>
> <br>
> <br>
> <br>
> <br>
> -- <br>
> Sébastien Loisel<br>
> Assistant Professor<br>
> Department of Mathematics, Heriot-Watt University<br>
> Riccarton, EH14 4AS, United Kingdom<br>
> web: <a href="http://www.ma.hw.ac.uk/~loisel/" rel="noreferrer" target="_blank">http://www.ma.hw.ac.uk/~loisel/</a><br>
> email: S.Loisel at <a href="http://hw.ac.uk" rel="noreferrer" target="_blank">hw.ac.uk</a><br>
> phone: +44 131 451 3234<br>
> fax: +44 131 451 3249<br>
> <br>
<br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>