[petsc-users] A way to distribute 3D arrays.

[Matthew Knepley]

On Sat, Feb 18, 2017 at 2:44 PM, Manuel Valera <mvalera at mail.sdsu.edu>
wrote:

> thanks guys that helped a lot!
>
> I think i got it know, i copy the code i created in case you want to
> suggest something or maybe use it as example...
>
> I have a bonus question: Im going to operate next in several arrays at the
> same time, i created them using their slightly different layouts, one DA
> each, a peer told me different DA are not guaranteed to share the same
> corners, is this correct? if so, is there a way to enforce that? Thanks
>

I do not understand. Do these DAs have the same size, meaning number of
vertices in X, Y, and Z? If so, why are you
using separate arrays. You can just use 'dof' to put all your data there,
and strip out subvectors using VecStride
functions, such as


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecStrideScatter.html

If the DAs are different sizes, how would they ever have the same corners.

  Thanks,

     Matt


> SUBROUTINE DAs(da,veczero,localv,array)
>
>        ! use PetscObjects, only :: ierr
>
>        ! Umbrella program to update and communicate the arrays in a
>        ! distributed fashion using the DMDA objects from PETSc.
>        ! Manuel Valera 1/20/17
>
>        ! Arguments:
>        ! da = DMDA array (3D) already created and setup
>        ! veczero =
>        ! globalv =
>        ! localv = local chunk each processor works in.
>        ! array = the array to be petscified. ONLY 3D ARRAYS as now.
>        ! arrayp = the petsc version of the array, may be not needed.
>        ! cta = the scatter context to translate globalv to localv in DA
> indices
>
>         Vec,intent(inout)             :: veczero,localv
>         Vec                           :: natural,globalv
>         PetscScalar,dimension(:,:,:)  :: array
>         PetscScalar,pointer,dimension(:,:,:) :: arrayp
>         DM,intent(inout)              :: da
>         VecScatter                    :: cta
>         PetscInt                      ::xind,yind,zind,xwidth,
> ywidth,zwidth
>
>
>
>         !Debug:
>         !print*,"INFO BEFORE SCATTER:"
>         call  DAinf(da,xind,yind,zind,xwidth,ywidth,zwidth)
>         print*, SIZE(array,1),'x',SIZE(array,2),'x',SIZE(array,3)
>         !
>
>         call DMCreateLocalVector(da,localv,ierr)
>
>         !The following if-block is an attempt to scatter the arrays as
> vectors
>         !loaded from master node.
>
>       if (sizex/= 1)then
>         !PETSC Devs suggested:
>         !Buffer 3d array natural has diff ordering than DMDAs:
>         call DMDACreateNaturalVector(da,natural,ierr)
>         call DMDAVecGetArrayF90(da,natural,arrayp,ierr)
>         !**fill up veczero***:
>         !TODO
>         if(rank==0)then
>         arrayp = array(xind+1:xwidth,yind+1:ywidth,zind+1:zwidth)
>         endif
>         call DMDAVecRestoreArrayF90(da,natural,arrayp,ierr)
>         !call DMRestoreNaturalVector(da,natural,ierr)      !???? not
> needed?
>
>         !Distributing array:
>         call VecScatterCreateToZero(natural,cta,veczero,ierr)
>         call VecScatterBegin(cta,veczero,natural,INSERT_VALUES,SCATTER_
> REVERSE,ierr)
>         call VecScatterEnd(cta,veczero,natural,INSERT_VALUES,SCATTER_
> REVERSE,ierr)
>
>
>         call DMCreateGlobalVector(da,globalv,ierr)
>         call DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,globalv,
> ierr)
>         call DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,globalv,ierr)
>
>         call DMGlobalToLocalBegin(da,globalv,INSERT_VALUES,localv,ierr)
>         call DMGlobalToLocalEnd(da,globalv,INSERT_VALUES,localv,ierr)
>
>       else
>
>         print*, "This should be called in serial only."
>
>         call DMDAVecGetArrayF90(da,localv,arrayp,ierr)
>
>         arrayp = array
>
>         call DMDAVecRestoreArrayF90(da,localv,arrayp,ierr)
>         call DMRestoreLocalVector(da,localv,ierr)
>       endif
>
>       call  VecDestroy(globalv,ierr)
>       call  VecScatterDestroy(cta,ierr)
>       call  VecDestroy(natural,ierr)
>
>
>
>  END SUBROUTINE DAs
>
>
>
> On Wed, Feb 15, 2017 at 7:31 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>>
>>   You can do this in the style of
>>
>>  VecLoad_Binary_DA()
>>
>> first you take your sequential vector and make it parallel in the
>> "natural" ordering for 3d arrays
>>
>>   DMDACreateNaturalVector(da,&natural);
>>   VecScatterCreateToZero(natural,&scatter,&veczero);  /* veczero is of
>> full size on process 0 and has zero entries on all other processes*/
>>   /* fill up veczero */
>>    VecScatterBegin(scatter,veczero,natural,INSERT_VALUES,SCATT
>> ER_REVERSE);
>> VecScatterEnd(scatter,veczero,natural,INSERT_VALUES,SCATTER_REVERSE);
>>
>>  and then move it into the PETSc DMDA parallel ordering vector with
>>
>>   ierr = DMCreateGlobalVector(da,&xin);CHKERRQ(ierr);
>>   ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKER
>> RQ(ierr);
>>   ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ
>> (ierr);
>>
>>
>> > On Feb 15, 2017, at 7:16 PM, Manuel Valera <mvalera at mail.sdsu.edu>
>> wrote:
>> >
>> > Hello,
>> >
>> > My question this time is just if there is a way to distribute a 3D
>> array whos located at Zero rank over the processors, if possible using the
>> DMDAs, i'm trying not to do a lot of initialization I/O in parallel.
>> >
>> > Thanks for your time,
>> >
>> > Manuel
>>
>>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170218/324ef69f/attachment.html>