[petsc-users] Regarding ksp ex42 - Citations

[Lawrence Mitchell]

> On 21 Jul 2016, at 09:41, domenico lahaye <domenico_lahaye at yahoo.com> wrote:
> 
> Thanks.
> 
> KSPSetOperators() allows to precondition A^h with M^h.
> This is lovely and great as it allows to implement the shifted Laplace
> preconditioner for the Helmholtz equation.
> 
> Recently I managed to implement shifted Laplace using the DMDA
> infrastructure in 2D. This implementation avoids having to construct
> the hierarchy in Matlab as we did previously.
> 
> In next stage we would like to precondition A^H with M^H on a sequence
> of coarser grids. This is what Calandra does on two levels and what we do
> on multiple levels.
> 
> We currently have an implement in which we construct the hierarchy on A^h
> and M^h in Matlab, we read the hierarchy in PETSc, traverse the hierarchy and
> do SetOperators and do a lot more of dark magic and witch craft by combining
> preconditioners in a additive and multiplicative fashion.
> 
> It would be lovely to obtain a more readable piece of code.
> 
> I am not sure what kind of additional callbacks I need. My first guess here
> would be a multilevel extension of SetOperators allowing to define M^H
> a preconditioner for A^H on a sequence of coarser levels. But I currently
> fail to oversee the whole matter.
> 
> An alternative is to build a fragile code on top of DMDA first and get back
> to you with more informed guesses on what kind of call backs I precisely need.
> I think I prefer to go with this option.
> 
> Does this sound reasonable?


It sounds like what you need is that the coarse DM should have a way of building the operators via a callback.  I think this is already available.  Rather than doing KSPSetOperators.  You do KSPSetComputeOperators, providing the function to be called to build the operator.  Now, you need a way for the coarse grids to allocate the matrices that will be used for your operators.  If you have a DMDA, this is set up for you because the KSP calls DMCreateMatrix and the DMDA knows how to create a matrix.

One wrinkle here is that the interface doesn't currently support making separate matrices for A and M.  The code currently does (in KSPSetUp):

if (using_dm) {
    DMCreateMatrix(ksp->dm, &A);
    KSPSetOperators(ksp, A, A);
    ...
}

For your needs you'd need this to be:

if (using_dm) {
   DMCreateMatrices(ksp->dm, &A, &P);
   KSPSetOperators(ksp, A, P)
   ...
}

I think.

Adding this call should not be too hard, there have been discussions before about it.  See, for example, the thread here: http://lists.mcs.anl.gov/pipermail/petsc-dev/2015-March/017130.html

(which started here http://lists.mcs.anl.gov/pipermail/petsc-dev/2015-February/017008.html)

I note I never got round to making the suggested changes there.

Cheers,

Lawrence
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 455 bytes
Desc: Message signed with OpenPGP using GPGMail
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160721/69de7e61/attachment.pgp>