Actual source code: ex20.c

  2: static char help[] = "Tests SDALocalToLocalxxx().\n\n";

 4:  #include petscda.h

  6: /*
  7:          For testing purposes this example also creates a 
  8:    DA context. Actually codes using SDA routines will probably 
  9:    not also work with DA contexts.
 10: */

 14: int main(int argc,char **argv)
 15: {
 16:   PetscMPIInt    size,rank;
 17:   PetscInt       M=8,dof=1,stencil_width=1,i,start,end,P=5,N = 6,m=PETSC_DECIDE,n=PETSC_DECIDE,p=PETSC_DECIDE,pt = 0,st = 0;
 19:   PetscTruth     flg2,flg3,flg;
 20:   DAPeriodicType periodic = DA_NONPERIODIC;
 21:   DAStencilType  stencil_type = DA_STENCIL_STAR;
 22:   DA             da;
 23:   SDA            sda;
 24:   Vec            local,global,local_copy;
 25:   PetscScalar    value,*in,*out;
 26:   PetscReal      norm,work;
 27:   PetscViewer    viewer;
 28:   char           filename[PETSC_MAX_PATH_LEN];
 29:   FILE           *file;


 32:   PetscInitialize(&argc,&argv,(char*)0,help);

 34:   PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
 35:   PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
 36:   PetscOptionsGetInt(PETSC_NULL,"-P",&P,PETSC_NULL);
 37:   PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);
 38:   PetscOptionsGetInt(PETSC_NULL,"-stencil_width",&stencil_width,PETSC_NULL);
 39:   PetscOptionsGetInt(PETSC_NULL,"-periodic",&pt,PETSC_NULL);
 40:   periodic = (DAPeriodicType) pt;
 41:   PetscOptionsGetInt(PETSC_NULL,"-stencil_type",&st,PETSC_NULL);
 42:   stencil_type = (DAStencilType) st;

 44:   PetscOptionsHasName(PETSC_NULL,"-1d",&flg2);
 45:   PetscOptionsHasName(PETSC_NULL,"-2d",&flg2);
 46:   PetscOptionsHasName(PETSC_NULL,"-3d",&flg3);
 47:   if (flg2) {
 48:     DACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,m,n,dof,stencil_width,0,0,&da);
 49:     SDACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,m,n,dof,stencil_width,0,0,&sda);
 50:   } else if (flg3) {
 51:     DACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,P,m,n,p,dof,stencil_width,0,0,0,&da);
 52:     SDACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,P,m,n,p,dof,stencil_width,0,0,0,&sda);
 53:   }
 54:   else {
 55:     DACreate1d(PETSC_COMM_WORLD,periodic,M,dof,stencil_width,PETSC_NULL,&da);
 56:     SDACreate1d(PETSC_COMM_WORLD,periodic,M,dof,stencil_width,PETSC_NULL,&sda);
 57:   }

 59:   DACreateGlobalVector(da,&global);
 60:   DACreateLocalVector(da,&local);
 61:   VecDuplicate(local,&local_copy);
 62: 
 63:   /* zero out vectors so that ghostpoints are zero */
 64:   value = 0;
 65:   VecSet(local,value);
 66:   VecSet(local_copy,value);

 68:   VecGetOwnershipRange(global,&start,&end);
 69:   for (i=start; i<end; i++) {
 70:     value = i + 1;
 71:     VecSetValues(global,1,&i,&value,INSERT_VALUES);
 72:   }
 73:   VecAssemblyBegin(global);
 74:   VecAssemblyEnd(global);

 76:   DAGlobalToLocalBegin(da,global,INSERT_VALUES,local);
 77:   DAGlobalToLocalEnd(da,global,INSERT_VALUES,local);

 79:   flg  = PETSC_FALSE;
 80:   PetscOptionsGetTruth(PETSC_NULL,"-same_array",&flg,PETSC_NULL);
 81:   if (flg) {
 82:     /* test the case where the input and output array is the same */
 83:     VecCopy(local,local_copy);
 84:     VecGetArray(local_copy,&in);
 85:     VecRestoreArray(local_copy,PETSC_NULL);
 86:     SDALocalToLocalBegin(sda,in,INSERT_VALUES,in);
 87:     SDALocalToLocalEnd(sda,in,INSERT_VALUES,in);
 88:   } else {
 89:     VecGetArray(local,&out);
 90:     VecRestoreArray(local,PETSC_NULL);
 91:     VecGetArray(local_copy,&in);
 92:     VecRestoreArray(local_copy,PETSC_NULL);
 93:     SDALocalToLocalBegin(sda,out,INSERT_VALUES,in);
 94:     SDALocalToLocalEnd(sda,out,INSERT_VALUES,in);
 95:   }

 97:   flg  = PETSC_FALSE;
 98:   PetscOptionsGetTruth(PETSC_NULL,"-save",&flg,PETSC_NULL);
 99:   if (flg) {
100:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
101:     sprintf(filename,"local.%d",rank);
102:     PetscViewerASCIIOpen(PETSC_COMM_SELF,filename,&viewer);
103:     PetscViewerASCIIGetPointer(viewer,&file);
104:     VecView(local,viewer);
105:     fprintf(file,"Vector with correct ghost points\n");
106:     VecView(local_copy,viewer);
107:     PetscViewerDestroy(viewer);
108:   }

110:   VecAXPY(local_copy,-1.0,local);
111:   VecNorm(local_copy,NORM_MAX,&work);
112:   MPI_Allreduce(&work,&norm,1,MPIU_REAL,MPI_MAX,PETSC_COMM_WORLD);
113:   if (norm != 0.0) {
114:     MPI_Comm_size(PETSC_COMM_WORLD,&size);
115:     PetscPrintf(PETSC_COMM_WORLD,"Norm of difference %G should be zero\n",norm);
116:     PetscPrintf(PETSC_COMM_WORLD,"  Number of processors %d\n",size);
117:     PetscPrintf(PETSC_COMM_WORLD,"  M,N,P,dof %D %D %D %D\n",M,N,P,dof);
118:     PetscPrintf(PETSC_COMM_WORLD,"  stencil_width %D stencil_type %d periodic %d\n",stencil_width,(int)stencil_type,(int)periodic);
119:     PetscPrintf(PETSC_COMM_WORLD,"  dimension %d\n",1 + (int) flg2 + (int) flg3);
120:   }
121: 
122:   DADestroy(da);
123:   SDADestroy(sda);
124:   VecDestroy(local_copy);
125:   VecDestroy(local);
126:   VecDestroy(global);
127:   PetscFinalize();
128:   return 0;
129: }