[petsc-users] DMCreateGlobalVector in fortran.

Sun Feb 5 12:13:26 CST 2017

  Are your -1 index locations used for "boundary conditions", that is extra variables that are used to update the other variables but don't get updated themselves?

  If so you should look at using a DMBoundaryType of  DM_BOUNDARY_GHOSTED, this will result in the arrays returned from DMDAVecGetArrayF90() starting with -1. 

   But I urge you to start with a simple example in order to understand how DMDA and its vectors work. I think you may misunderstand it.

For example:

> call DMCreateGlobalVector(da,globalv,ierr)
> >         call DMCreateLocalVector(da,localv,ierr)
> >
> >         call DMGlobalToLocalBegin(da,globalv,INSERT_VALUES,localv,ierr)
> >         call DMGlobalToLocalEnd(da,globalv,INSERT_VALUES,localv,ierr)
> >

  You create a global and local array which have all zero entries (PETSc vectors by default have all zeros) and then you scatter the zero entries from the global into the local. This scatter serves no purpose, unless you are just checking that it does not crash.

   Barry

> On Feb 4, 2017, at 10:19 PM, Manuel Valera <mvalera at mail.sdsu.edu> wrote:
> 
> Ok, let me try to clarify a little bit,
> 
> I have a u0 3d array with the mentioned indices (starting from -1), im creating a dmda with the sizes in the constructor, no problem with that.
> 
> Next, im trying to distribute this u0 array into global and local vectors, apparently i can do that since no error comes out. 
> 
> Now, when i get the array back with VecGetArrayF90(), the indices start now from 0, is that right? so the whole array have been shifted in its indices.
> 
> This changes things for me because my idea was to use the previous serial code to operate on the arrays, but this is written with indices starting from -1. There are dozens of arrays with the same structure, but they also operate with other arrays with indices starting from 1, at the same time, usually in triple nested loops which code the stencil used within them, to update and correct the velocities after solving for pressure, per example.
> 
> So it would be much easier to shift the indices of these arrays by one, rather than changing all the code to match the new indices starting from zero,
> 
> Is there a way to do this?
> 
> Thanks,
> 
> 
> 
> 
> On Sat, Feb 4, 2017 at 5:46 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
>    For DMDAVecGetArrayF90 you need to declare the "array" arguments as Fortran pointers you don't declare them like
> 
> > u0(-1:IMax+2,-1:JMax+1,-1:KMax+1)
> 
> 
> 
> > On Feb 4, 2017, at 7:34 PM, Matthew Knepley <knepley at gmail.com> wrote:
> >
> > On Sat, Feb 4, 2017 at 7:30 PM, Manuel Valera <mvalera at mail.sdsu.edu> wrote:
> > Thanks, that helped a lot,
> >
> > Now im a little confused with my next error, in my serial code i have several arrays with dimensions declared as:
> >
> > u0(-1:IMax+2,-1:JMax+1,-1:KMax+1)
> >
> > and i create my DMDA for this array as:
> >
> > call CDA3(daub,IMax+4,JMax+3,KMax+3)
> > call DAs(daub,u0p,u0pl,u0)
> >
> > With the subroutines as:
> >
> >  SUBROUTINE CDA3(da3d,dim3dx,dim3dy,dim3dz)
> >
> >
> >
> >       DM, intent(inout)         :: da3d
> >       PetscInt, intent(in)      :: dim3dx,dim3dy,dim3dz
> >       DMBoundaryType            :: bx,by,bz
> >
> >
> >         bx = DM_BOUNDARY_NONE
> >         by = DM_BOUNDARY_NONE
> >         bz = DM_BOUNDARY_NONE
> >
> >       call DMDACreate3d(PETSC_COMM_WORLD,bx,by,bz,DMDA_STENCIL_BOX,dim3dx,dim3dy,dim3dz,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da3d,ierr)
> >
> >       call DMSetFromOptions(da3d,ierr)
> >       call DMSetUp(da3d,ierr)
> >
> >
> >  END SUBROUTINE CDA3
> >
> >  SUBROUTINE DAs(da,globalv,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 either 1d,2d or 3d, already created and setup
> >        ! globalv = global vector to be operated into
> >        ! localv = local chunk each processor works in.
> >        ! array = the array to be petscified. ONLY 3D ARRAYS rn.
> >        ! just in case.
> >
> >         Vec,intent(inout)             :: globalv,localv
> >
> >
> >         PetscScalar,dimension(:,:,:)           :: array
> >         DM,intent(inout)              :: da
> >
> >
> >
> >
> >
> >
> >         call DMCreateGlobalVector(da,globalv,ierr)
> >         call DMCreateLocalVector(da,localv,ierr)
> >
> >         call DMGlobalToLocalBegin(da,globalv,INSERT_VALUES,localv,ierr)
> >         call DMGlobalToLocalEnd(da,globalv,INSERT_VALUES,localv,ierr)
> >
> >
> >
> >  END SUBROUTINE DAs
> >
> > Which i later open into an array to operate into it, but when i try to access u0(:,j-2,:) it gives me the following error:
> >
> > Fortran runtime error: Index '-1' of dimension 2 of array 'u0' below lower bound of 0
> >
> > So i'm guessing the DMDA is creating the arrays with dimensions starting at 0 instead of -1 as it should, right? what can be done in this case?
> >
> > 1) The DMDA only takes sizes, not index bounds, in the constructor.
> >
> > 2) Its not the DMDA you must be talking about, but the array you get from DMDAVecGetArrayF90() or VecGetArrayF90()
> >
> > 3) The VecGetArrayF90() always start from 0 since it is indexed in local indices
> >
> > 4) The DMDVecGetArrayF90() arrays are indexed by global indices for the patch your process owns (including ghosts if its a local vector)
> >
> >   Thanks,
> >
> >      Matt
> >
> > Thanks,
> > Manuel
> >
> > On Sat, Feb 4, 2017 at 4:18 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> >    See:  http://www.mcs.anl.gov/petsc/documentation/changes/dev.html  search for DMDACreate  and update your code.
> >
> >   Barry
> >
> >
> > > On Feb 4, 2017, at 6:02 PM, Manuel Valera <mvalera at mail.sdsu.edu> wrote:
> > >
> > > Hello devs,
> > >
> > > Ive been creating the framework to distribute the arrays from my model using the global/local vectors and the DMDAs, for this i created a small subroutine that is responsible for making the necessary calls:
> > >
> > >  SUBROUTINE DAs(da,globalv,localv)
> > >
> > >        ! 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 either 1d,2d or 3d, already created and setup
> > >         ! globalv = global vector to be operated into
> > >         ! localv = local chunk each processor works in.
> > >
> > >
> > >
> > >         Vec,intent(inout)             :: globalv,localv
> > >
> > >
> > >
> > >         DM,intent(inout)              :: da
> > >
> > >         call DMCreateGlobalVector(da,globalv,ierr)
> > >         call DMCreateLocalVector(da,localv,ierr)
> > >
> > >         call DMGlobalToLocalBegin(da,globalv,INSERT_VALUES,localv,ierr)
> > >         call DMGlobalToLocalEnd(da,globalv,INSERT_VALUES,localv,ierr)
> > >
> > >
> > >
> > > END SUBROUTINE DAs
> > >
> > >
> > > I call this from another program where i have declared a vector as:
> > >
> > > #include "petsc/finclude/petscdm.h"
> > > #include "petsc/finclude/petscdmda.h"
> > >
> > > use petsc
> > > use petscdm
> > > use petscdmda
> > >
> > > IMPLICIT NONE
> > >
> > > Vec                                        :: xpdmda,yp,zp,xcp,ycp,zcp
> > > DM                                         :: da,dac,dau,dav,daw
> > >
> > > CONTAINS
> > > SUBROUTINE InitGridP
> > >
> > > IMPLICIT NONE
> > > !PetscInt     :: ierr
> > >
> > > ! CREATE DMDA:
> > >
> > > call CDA3(da,IMax-1,JMax-1,KMax-1)
> > >
> > > !xpdmda,yp,zp have the same dimensions than x,y,z
> > > !xpl,ypl,zpl are the local versions of xp,yp,zp
> > >
> > > !We create the global and local vectors and communicate their values:
> > > call DAs(da,xpdmda,xpl,x)
> > >
> > > And at this call it crashes with the following message:
> > >
> > > [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> > > [0]PETSC ERROR: Null argument, when expecting valid pointer
> > > [0]PETSC ERROR: Null Object: Parameter # 2
> > > [0]PETSC ERROR: See
> > > http://www.mcs.anl.gov/petsc/documentation/faq.html
> > >  for trouble shooting.
> > > [0]PETSC ERROR: Petsc Development GIT revision: v3.7.5-2949-gc3c607c  GIT Date: 2017-01-20 17:34:16 -0600
> > > [0]PETSC ERROR: ./ucmsLEP on a arch-linux2-c-debug named ocean by valera Sat Feb  4 15:53:53 2017
> > > [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich --download-hdf5 --download-netcdf --download-hypre --download-metis --download-parmetis --download-trillinos --with-debugging=1
> > > [0]PETSC ERROR: #3 VecSetLocalToGlobalMapping() line 79 in /usr/dataC/home/valera/petsc/src/vec/vec/interface/vector.c
> > > [0]PETSC ERROR: #4 DMCreateGlobalVector_DA() line 41 in /usr/dataC/home/valera/petsc/src/dm/impls/da/dadist.c
> > > [0]PETSC ERROR: #5 DMCreateGlobalVector() line 840 in /usr/dataC/home/valera/petsc/src/dm/interface/dm.c
> > >
> > >
> > >
> > > So my question is, do i have to setup the vectors in any additional way? i followed the fortran DMDA examples closely and i didnt notice anything extra,
> > >
> > > Thanks,
> > >
> > > 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
> 
> 