[petsc-users] PCMG with matrix-free operators accessing DMDA

[Barry Smith]

> On Jul 28, 2016, at 12:48 PM, Valeria Mele <valeria.mele at unina.it> wrote:
> 
> Hi everyone,
> this time I am using PETSc to do something that is more complicated than my usual and I want to do it at the highest possible abstraction level.
> 
> To put it in a nutshell, my intent is to build a parallel multigrid to solve a linear system via DM, KSP and PCMG (I would like to use DMMG but probably I should have the same problems or more).
   
   DMMG doesn't exist anymore. It was refactored away many years ago, its functionality is handled by PCMG and DM.


> 
> I created the distributed object, da, with DMDACreate3d, even if it is distributed (as yet) only in the x-dimension and has 3 dof.
> Then I create the KSP (type KSPRICHARDSON) and set the nonzero initial guess and PCMG as preconditioner. Here I start to tune the MG. 
> 
> The point is that I need to define all the operators as matrix-free, since they will do several operations on x to obtain y, and I am not familiar with the way to access all the elements or informations in the two levels involved and/or among the processors with a so-high level interface.
>  
> So please (please please please please) tell me if I correctly understand the mechanism or I am on the wrong way and clear my doubts.
> 
> That is, let's say that my operation for the shell are:
> 	• A_mult(Mat mat,Vec x, Vec y) //coefficients matrix
> in this case the level is only one but should I write it taking into account only the local data (I think so) and accessing them via the informations in da? 

   Yes. You can use VecGetDM(x,&da) to get the DMDA object

> 
> For example, if I use DMDAVecGetArray, DMDAVecGetCorners (or DMDAGetGhostCorners) and DMDAVecRestoreArray, will they retrieve informations from the right level each time (I am pretty sure that in some official examples it is done in this way)?
> 
> Or should I handle just Vecs as local structures with their usual indices (through VecGetArray and VecRestoreArray)?

   No, no, no because then you would need to mange all the structured grid information yourself, since the DMDA manages it for you you should use it.

> 	• P_mult(Mat mat,Vec x, Vec y) //interpolation matrix that is NOT conceptually the traspose of Restriction
> in this case x and y will be from two different levels (respectively L and L+1), so, if I retrieve informations from the da... how can I access the two at different levels?

   Use VecGetDM(x,  and VecGetDM(y to get access to both DMDA.

> 
> I am sorry if it seems that they are trivial questions, and I will be grateful to anyone will help me. 

  Additional information. Since the PCMG will be requesting the matrices and the interpolation/restriction operations (rather than you setting them into each level of multigrid) you will need to use DMShellSetCreateMatrix() and DMShellSetCreateInterpolation() and DMShellSetCreateRestriction() to provide the routines that will create the Shell matrices you need to represent the operators on the levels and the restriction and interpolation (Even though you are using a DMDA you can still call these routines).

   Barry

> 
> Thanks a lot,
> Valeria
> 
> 
> 
> 
> 
> 
> ---------------------------------------------------------------------------------------------
> PhD Valeria Mele
> 
> University of Naples Federico II
> Department of Mathematics and Applications "R. Caccioppoli"
> Complesso Universitario M.S. Angelo, Via Cinthia
> 80126 Naples 
> ---------------------------------------------------------------------------------------------