Re: [petsc-dev] Note about CPR preconditioner in PCFIELDSPLIT man page?

[Mills, Richard Tran via petsc-dev]

On 3/3/19 8:21 PM, Jed Brown wrote:

"Mills, Richard Tran"  writes:



Alright, I see the source of my confusion: I've now seen three different that 
are all called a "Constrained Pressure Residual" preconditioner. What I have 
generally seen is exactly what Matt describes: R is just a restriction to a 
subset of DoFs (pressure, and sometimes saturation), and P is just R' (using 
Matlab notation here). A CPR preconditioner in which P differs from R' is not 
something I've seen until now.



Please read the code again.  P = R' for this case, but R is restriction to the 
sum of field components, not a subset of the variables.

Oops -- yes, I see. Thanks, Jed. Elsewhere that I've seen a CPR preconditioner 
used, the restriction is only to a subset of the variables. In that case, the 
preconditioner can be constructed with only PCFIELDSPLIT and PCCOMPOSITE, but 
in the case described in the email messages to Barry, I do see that PCGALERKIN 
is needed.

--Richard





I'll see if I can think of a succinct change to the documentation for 
PCFIELDSPLIT to describe the two cases and submit a pull request.

--Richard

On 3/3/19 8:38 AM, Jed Brown wrote:

Matthew Knepley 

 writes:



On Sun, Mar 3, 2019 at 12:58 AM Jed Brown via petsc-dev <
petsc-dev@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.



Look at the source.  These are not subsets:

 /*
 Apply the restriction operator for the Galkerin problem
 */
 PetscErrorCode ApplyR(Mat A, Vec x,Vec y)
 {
   PetscErrorCode ierr;
   PetscInt   b;
   PetscFunctionBegin;
   ierr = VecGetBlockSize(x,&b);CHKERRQ(ierr);
   ierr = VecStrideGather(x,0,y,INSERT_VALUES);CHKERRQ(ierr);
   for (PetscInt k=1;k

Re: [petsc-dev] Note about CPR preconditioner in PCFIELDSPLIT man page?

[Jed Brown via petsc-dev]

"Mills, Richard Tran"  writes:

> Alright, I see the source of my confusion: I've now seen three different that 
> are all called a "Constrained Pressure Residual" preconditioner. What I have 
> generally seen is exactly what Matt describes: R is just a restriction to a 
> subset of DoFs (pressure, and sometimes saturation), and P is just R' (using 
> Matlab notation here). A CPR preconditioner in which P differs from R' is not 
> something I've seen until now.

Please read the code again.  P = R' for this case, but R is restriction to the 
sum of field components, not a subset of the variables.

> I'll see if I can think of a succinct change to the documentation for 
> PCFIELDSPLIT to describe the two cases and submit a pull request.
>
> --Richard
>
> On 3/3/19 8:38 AM, Jed Brown wrote:
>
> Matthew Knepley  writes:
>
>
>
> On Sun, Mar 3, 2019 at 12:58 AM Jed Brown via petsc-dev <
> petsc-dev@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.
>
>
>
> Look at the source.  These are not subsets:
>
>  /*
>  Apply the restriction operator for the Galkerin problem
>  */
>  PetscErrorCode ApplyR(Mat A, Vec x,Vec y)
>  {
>PetscErrorCode ierr;
>PetscInt   b;
>PetscFunctionBegin;
>ierr = VecGetBlockSize(x,&b);CHKERRQ(ierr);
>ierr = VecStrideGather(x,0,y,INSERT_VALUES);CHKERRQ(ierr);
>for (PetscInt k=1;k VecStrideGather(x,k,y,ADD_VALUES);CHKERRQ(ierr);}
>PetscFunctionReturn(0);
>  }
>
>  /*
>  Apply the interpolation operator for the Galerkin problem
>  */
>  PetscErrorCode ApplyP(Mat A, Vec x,Vec y)
>  {
>PetscErrorCode ierr;
>PetscInt   offset = 1;
>PetscFunctionBegin;
>ierr = VecStrideScatter(x,offset,y,INSERT_VALUES);CHKERRQ(ierr);
>PetscFunctionReturn(0);
>  }

Re: [petsc-dev] Note about CPR preconditioner in PCFIELDSPLIT man page?

[Mills, Richard Tran via petsc-dev]

Alright, I see the source of my confusion: I've now seen three different that 
are all called a "Constrained Pressure Residual" preconditioner. What I have 
generally seen is exactly what Matt describes: R is just a restriction to a 
subset of DoFs (pressure, and sometimes saturation), and P is just R' (using 
Matlab notation here). A CPR preconditioner in which P differs from R' is not 
something I've seen until now.

I'll see if I can think of a succinct change to the documentation for 
PCFIELDSPLIT to describe the two cases and submit a pull request.

--Richard

On 3/3/19 8:38 AM, Jed Brown wrote:

Matthew Knepley  writes:



On Sun, Mar 3, 2019 at 12:58 AM Jed Brown via petsc-dev <
petsc-dev@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.



Look at the source.  These are not subsets:

 /*
 Apply the restriction operator for the Galkerin problem
 */
 PetscErrorCode ApplyR(Mat A, Vec x,Vec y)
 {
   PetscErrorCode ierr;
   PetscInt   b;
   PetscFunctionBegin;
   ierr = VecGetBlockSize(x,&b);CHKERRQ(ierr);
   ierr = VecStrideGather(x,0,y,INSERT_VALUES);CHKERRQ(ierr);
   for (PetscInt k=1;k

Re: [petsc-dev] Note about CPR preconditioner in PCFIELDSPLIT man page?

[Jed Brown via petsc-dev]

Matthew Knepley  writes:

> On Sun, Mar 3, 2019 at 12:58 AM Jed Brown via petsc-dev <
> petsc-dev@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.

Look at the source.  These are not subsets:

 /*
 Apply the restriction operator for the Galkerin problem
 */
 PetscErrorCode ApplyR(Mat A, Vec x,Vec y)
 {
   PetscErrorCode ierr;
   PetscInt   b;
   PetscFunctionBegin;
   ierr = VecGetBlockSize(x,&b);CHKERRQ(ierr);
   ierr = VecStrideGather(x,0,y,INSERT_VALUES);CHKERRQ(ierr);
   for (PetscInt k=1;k

Re: [petsc-dev] Note about CPR preconditioner in PCFIELDSPLIT man page?

[Jed Brown via petsc-dev]

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.

--- Begin Message ---

   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=1, initial guess is zero
  tolerances:  relative=1e-05, absolute=1e-50, divergence=1.
  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=1, initial guess is zero
tolerances:  relative=1e-05, absolute=1e-50, divergence=1.
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 coarse1
  HYPRE BoomerAMG: Relax down  symmetric-SOR/Jacobi
  HYPRE BoomerAMG: Relax upsymmetric-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 typelocal
  HYPRE BoomerAMG: Coarsen typeFalgout
  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=10
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=10
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_

[petsc-dev] Note about CPR preconditioner in PCFIELDSPLIT man page?

[Mills, Richard Tran via petsc-dev]

Folks,

There is a note in the PCFIELDSPLIT manual page about implementing a so-called 
"Constrained Pressure Residual" (CPR) preconditioner:

"The Constrained Pressure Preconditioner (CPR) can be implemented using 
PCCOMPOSITE
 with 
PCGALERKIN.
 CPR first solves an R A P subsystem, updates the residual on all variables 
(PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)),
 and then applies a simple ILU like preconditioner on all the variables."

I'm not sure why there is a reference to PCGALERKIN in there. Although CPR *is* 
using a Galerkin-type projection, I believe that this can all be set up on the 
command line by using PCFIELDSPLIT and PCCOMPOSITE. It seems that CPR 
preconditioners are defined in a few different ways, but I believe I can set up 
something like this in PFLOTRAN (where field 0 is pressure, 1 saturation, and 2 
is temperature) by doing something like this:

-flow_ksp_type fgmres -flow_pc_type composite \
-flow_pc_composite_type multiplicative -flow_pc_composite_pcs 
fieldsplit,bjacobi \
-flow_sub_0_ksp_type fgmres -flow_sub_0_pc_fieldsplit_type additive \
-flow_sub_0_pc_fieldsplit_0_fields 0 -flow_sub_0_pc_fieldsplit_1_fields 1,2 \
-flow_sub_0_fieldsplit_0_ksp_type richardson -flow_sub_0_fieldsplit_1_ksp_type 
preonly \
-flow_sub_0_fieldsplit_0_pc_type hypre -flow_sub_0_fieldsplit_0_pc_hypre_type 
boomeramg \
-flow_sub_0_fieldsplit_1_pc_type none -flow_sub_0_fieldsplit_1_sub_pc_type none 
\
-flow_sub_0_fieldsplit_1_ksp_max_it 10 -flow_sub_1_sub_pc_type ilu

Am I missing something? If the above setup is correct, it seems like the 
mention of PCGALERKIN is a bit confusing, since this is the PCFIELDSPLIT man 
page.

--Richard