[petsc-users] Extract subdomain values from DMDA

[Barry Smith]

  Tina,

   Is the matrix also obtained by extracting it from a larger matrix on the larger grid or do you generate the matrix for the smaller problem directly on the smaller grid?

   In order to "pull out" the part you want, you need to create two index sets that indicate the values you are taking from the larger array and where you are putting them in the smaller array. This is easier to do in the natural ordering than the PETSc ordering on the DMDA so I would do the following create a DMDA for both the large grid and the smaller grid (I would write a test code in 2d first because it is much easier to debug in 2d), then you can create Vecs to hold the vectors in  "natural ordering" for each grid with DMDACreateNaturalVector().
You can use DMDAGlobalToNaturalBegin()/DMDAGlobalToNaturalEnd() to take the larger vector into natural ordering, then use a VecScatterBegin/End() to select the parts you want for the smaller grid, then use DMDANaturalToGlobalBegin()/DMDANaturalToGlobalEnd() to move the results to the DMDA ordering on the smaller grid, solve the system and then do the reverse process to get the values back onto the larger grid. The only slightly tricky part is generating the two IS you need to create the VecScatter that you use to move from the natural ordering on the larger grid to the natural ordering on the smaller grid. But since you can create the natural ordering numbering from the i,j (in 2d or i,j,k in 3d) location of the grid points on the larger grid to the smaller grid this is not too difficult; but I cannot tell you how to do this since it depends on what locations are you taking from the larger grid to the smaller grid). Feel free to ask more questions as you work on it.

   Barry


> On Apr 22, 2017, at 8:37 PM, Tina Patel <tinap89 at yahoo.com> wrote:
> 
> Hello everyone,
> 
> I want to manipulate a global vector that has values only on the subgrid of a larger grid obtained by DMDACreate3d(). For example, I'd like to extract a 6x6x6 chunk from a 7x7x7 grid. 
> 
> At some point, I need to put that 6x6x6 chunk from the logical xyz cartesian coordinates into a natural vector "b" to solve Ax=b at every iteration, then map it back to xyz coordinate system to do further calculations with ghost values.
> 
> The closest thing I can compare it to is MPI_Type_create_subarray. Do I just use this or is there a better way?
> 
> 
> Thanks for your time,
> Tina