[petsc-users] How to extract a slice at a given coordinate and view it ?

[Barry Smith]

  I wrote a routine DMDAGetRay() that pulls a 1 dimensional slice out of a 2d DMDA and puts it on process 0. It uses AOApplicationToPetsc() so is not truly scalable but perhaps you could take a look at that. Since you say  "(I have lots of points in phi, so this can happen easily)" it may be ok for you to just stick the 2d slice on process 0 and then save it?

  Barry

Without using AOApplicationToPetsc() or something similar yes it is in general a nightmare.


> On Aug 28, 2015, at 10:15 PM, Timothée Nicolas <timothee.nicolas at gmail.com> wrote:
> 
> Hi,
> 
> I have been thinking for several hours about this problem and can't find an efficient solution, however I imagine this must be possible somehow with Petsc.
> 
> My problem is the following :
> 
> I work in 3D (R,Z,Phi) which makes my data quite heavy and I don't want to save all the data of all my fields, even just once in a while. Instead, I would like to save in a binary file a slice at a given angle, say phi=0.
> 
> As I did not find if it's natively possible in Petsc, I considered creating a second 2D DMDA, on which I can create 2D vectors and view them with the binary viewer. So far so good. However, upon creating the 2D DMDA, naturally the distribution of processors does not correspond to the distribution of the 3D DMDA. So I was considering creating global arrays, filling them with the data of the 3D array in phi=0, then doing an MPI_allgather to give the information to all the processors, to be able to read the array and fill the 2D Petsc Vector with it. So the code would be something along the lines of :
> 
>   PetscScalar, pointer :: gX2D(:,:,:)
>   PetscScalar, pointer :: gX(:,:,:,:)
>   ! LocalArray is locally filled                                                                                                                                                                       
>   ! It is transmitted to GlobalArray via MPI_Allgather                                                                                                                                                 
>   real(8)          :: LocalArray(user%dof,user%mr,user%mz)
>   real(8)          :: GlobalArray(user%dof,user%mr,user%mz)
> 
>   call DMDAVecGetArrayF90(da_phi0,X2D,gX2D,ierr)
>   call DMDAVecGetArrayF90(da,X,gX,ierr)
> 
>   do k = user%phis,user%phie
>      do j = user%zs,user%ze
>         do i = user%rs,user%re
>            do l=1,user%dof
>               if (k.eq.phi_print) then
>                  ! Numbering obtained with DMDAGetArrayF90 differs from usual                                                                                                                          
>                  LocalArray(l,i,j) = gX(l-1,i-1,j-1,k-1)
>               end if
>            end do
>         end do
>      end do
>   end do
> 
>   nvals = user%dof*user%rm*user%zm
> 
>   call MPI_AllGather(LocalArray(1,user%rs,user%zs),        &
>        &             nvals,MPI_REAL,                        &
>        &             GlobalArray,       &
>        &             nvals,MPI_REAL,MPI_COMM_WORLD,ierr)
> 
>   do j = zs2D,ze2D
>      do i = rs2D,re2D
>         do l=1,user%dof
>            gX2D(l-1,i-1,j-1) = GlobalArray(l,i,j)
>         end do
>      end do
>   end do
> 
>  call DMDAVecRestoreArrayF90(da_phi0,X2D,gX2D,ierr)
>   call DMDAVecRestoreArrayF90(da,X,gX,ierr)
> 
> The problem is that MPI_allgather is not at all that simple. Exchanging array information is much more complicated that I had anticipated ! See this long post on stackoverflow : 
> 
> http://stackoverflow.com/questions/17508647/sending-2d-arrays-in-fortran-with-mpi-gather
> 
> I could probably get it to work eventually, but it's pretty complicated, and I was wondering if there was not a simpler alternative I could not see. Besides, I am concerned about what could happen if the number of processors is so large that the 2D Vector gets less than 2 points per processor (I have lots of points in phi, so this can happen easily). Then Petsc would complain.
> 
> Does anybody have ideas ?
> 
> Best
> 
> Timothée