[petsc-users] PCDiagonalScale

[Klaij, Christiaan]

Barry,

Thanks for your reply. Good to know it's indeed a missing check.
I'm not following your suggestion to modify PETSc; it means I
would have to modify it on a number of machines, some of which I
can't even access.

For the time being, I'm considering scaling the system in my code
before passing it to PETSc. Besides, this would allow me to use
FGMRES in case GCR gives trouble.

Chris


Date: Wed, 2 Mar 2011 12:57:38 -0600
From: Barry Smith <bsmith at mcs.anl.gov>
Subject: Re: [petsc-users] PCDiagonalScale
To: PETSc users list <petsc-users at mcs.anl.gov>
Message-ID: <F3851033-54D8-492E-B521-D8932AC5BF91 at mcs.anl.gov>
Content-Type: text/plain; charset=us-ascii


On Mar 1, 2011, at 2:28 AM, Klaij, Christiaan wrote:

> Barry,
>
> I do not want to use it for the reason given in the manual.
>
> Instead I'm trying to use SIMPLE as a preconditioner for a colocated, cell-centered FV discretization of the incompressible Navier-Stokes equations. I'm using a matfree approach for the matrix and the (right) preconditioner. I would like to scale the system as D A M D^{-1} z = D b but I get the feeling that PCDiagonalScaleSet with GCR doesn't do anything.
>
> How about Jed's comment:
>
> "GCR also tolerates a variable preconditioner and it does not have such a check. I don't know if that means it can use diagonal scaling or just that someone forgot to check, but you could try it."

   KSPSolve_GCR() is not coded for using that extra diagonal scaling and does not check to see if you have requested it.  So, as you saw, if you turn it on it has no effect.

   I have added an error check in petsc-dev so it will tell you that it does not use it.  Thank you for letting us know about this bug.

   You may be able to add support for it.  The code is in src/ksp/ksp/impls/gcr/grc.c  Note that the current code has the lines

  ierr = PCApply( pc, r, s ); CHKERRQ(ierr); /* s = B^{-1} r */
    ierr = MatMult( A, s, v ); CHKERRQ(ierr);  /* v = A s */

  to fix the code to do the diagonal scaling you need to do a few things

1) replace the two lines above with a single call to KSP_PCApplyBAorAB() this will incorporate the scaling in the application of    A*B so instead it provides D A B D^{-1}

2) this is more difficult. Understand the code well enough to fix it so that it properly takes into account the diagonal scaling if the initial guess is nonzero and also the fact that you are solving D A B D^{-1} y = D b instead of A B z = b  and hence  A x = b where x = B z so x = B *D^{-1} y

  Barry


dr. ir. Christiaan Klaij
CFD Researcher
Research & Development
E mailto:C.Klaij at marin.nl
T +31 317 49 33 44

MARIN
2, Haagsteeg, P.O. Box 28, 6700 AA Wageningen, The Netherlands
T +31 317 49 39 11, F +31 317 49 32 45, I www.marin.nl