> On 21 Jul 2016, at 09:41, domenico lahaye <domenico_lah...@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
signature.asc
Description: Message signed with OpenPGP using GPGMail