1 static char help[] = "Tests various 3-dimensional DMDA routines.\n\n"; 2 3 #include <petscdm.h> 4 #include <petscdmda.h> 5 #include <petscao.h> 6 7 int main(int argc,char **argv) 8 { 9 PetscMPIInt rank; 10 PetscInt M = 3,N = 5,P=3,s=1,w=2,nloc,l,i,j,k,kk,m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE; 11 PetscInt Xs,Xm,Ys,Ym,Zs,Zm,iloc,*iglobal; 12 const PetscInt *ltog; 13 PetscInt *lx = NULL,*ly = NULL,*lz = NULL; 14 PetscBool test_order = PETSC_FALSE; 15 DM da; 16 PetscViewer viewer; 17 Vec local,global; 18 PetscScalar value; 19 DMBoundaryType bx = DM_BOUNDARY_NONE,by = DM_BOUNDARY_NONE,bz = DM_BOUNDARY_NONE; 20 DMDAStencilType stencil_type = DMDA_STENCIL_BOX; 21 AO ao; 22 PetscBool flg = PETSC_FALSE; 23 24 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 25 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",300,0,400,300,&viewer)); 26 27 /* Read options */ 28 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-NX",&M,NULL)); 29 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-NY",&N,NULL)); 30 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-NZ",&P,NULL)); 31 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 32 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 33 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL)); 34 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-s",&s,NULL)); 35 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-w",&w,NULL)); 36 flg = PETSC_FALSE; 37 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-star",&flg,NULL)); 38 if (flg) stencil_type = DMDA_STENCIL_STAR; 39 flg = PETSC_FALSE; 40 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-box",&flg,NULL)); 41 if (flg) stencil_type = DMDA_STENCIL_BOX; 42 43 flg = PETSC_FALSE; 44 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-xperiodic",&flg,NULL)); 45 if (flg) bx = DM_BOUNDARY_PERIODIC; 46 flg = PETSC_FALSE; 47 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-xghosted",&flg,NULL)); 48 if (flg) bx = DM_BOUNDARY_GHOSTED; 49 flg = PETSC_FALSE; 50 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-xnonghosted",&flg,NULL)); 51 52 flg = PETSC_FALSE; 53 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-yperiodic",&flg,NULL)); 54 if (flg) by = DM_BOUNDARY_PERIODIC; 55 flg = PETSC_FALSE; 56 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-yghosted",&flg,NULL)); 57 if (flg) by = DM_BOUNDARY_GHOSTED; 58 flg = PETSC_FALSE; 59 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-ynonghosted",&flg,NULL)); 60 61 flg = PETSC_FALSE; 62 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-zperiodic",&flg,NULL)); 63 if (flg) bz = DM_BOUNDARY_PERIODIC; 64 flg = PETSC_FALSE; 65 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-zghosted",&flg,NULL)); 66 if (flg) bz = DM_BOUNDARY_GHOSTED; 67 flg = PETSC_FALSE; 68 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-znonghosted",&flg,NULL)); 69 70 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testorder",&test_order,NULL)); 71 72 flg = PETSC_FALSE; 73 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-distribute",&flg,NULL)); 74 if (flg) { 75 PetscCheckFalse(m == PETSC_DECIDE,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must set -m option with -distribute option"); 76 CHKERRQ(PetscMalloc1(m,&lx)); 77 for (i=0; i<m-1; i++) lx[i] = 4; 78 lx[m-1] = M - 4*(m-1); 79 PetscCheckFalse(n == PETSC_DECIDE,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must set -n option with -distribute option"); 80 CHKERRQ(PetscMalloc1(n,&ly)); 81 for (i=0; i<n-1; i++) ly[i] = 2; 82 ly[n-1] = N - 2*(n-1); 83 PetscCheckFalse(p == PETSC_DECIDE,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must set -p option with -distribute option"); 84 CHKERRQ(PetscMalloc1(p,&lz)); 85 for (i=0; i<p-1; i++) lz[i] = 2; 86 lz[p-1] = P - 2*(p-1); 87 } 88 89 /* Create distributed array and get vectors */ 90 CHKERRQ(DMDACreate3d(PETSC_COMM_WORLD,bx,by,bz,stencil_type,M,N,P,m,n,p,w,s,lx,ly,lz,&da)); 91 CHKERRQ(DMSetFromOptions(da)); 92 CHKERRQ(DMSetUp(da)); 93 CHKERRQ(PetscFree(lx)); 94 CHKERRQ(PetscFree(ly)); 95 CHKERRQ(PetscFree(lz)); 96 CHKERRQ(DMView(da,viewer)); 97 CHKERRQ(DMCreateGlobalVector(da,&global)); 98 CHKERRQ(DMCreateLocalVector(da,&local)); 99 100 /* Set global vector; send ghost points to local vectors */ 101 value = 1; 102 CHKERRQ(VecSet(global,value)); 103 CHKERRQ(DMGlobalToLocalBegin(da,global,INSERT_VALUES,local)); 104 CHKERRQ(DMGlobalToLocalEnd(da,global,INSERT_VALUES,local)); 105 106 /* Scale local vectors according to processor rank; pass to global vector */ 107 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 108 value = rank; 109 CHKERRQ(VecScale(local,value)); 110 CHKERRQ(DMLocalToGlobalBegin(da,local,INSERT_VALUES,global)); 111 CHKERRQ(DMLocalToGlobalEnd(da,local,INSERT_VALUES,global)); 112 113 if (!test_order) { /* turn off printing when testing ordering mappings */ 114 if (M*N*P<40) { 115 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nGlobal Vector:\n")); 116 CHKERRQ(VecView(global,PETSC_VIEWER_STDOUT_WORLD)); 117 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n")); 118 } 119 } 120 121 /* Send ghost points to local vectors */ 122 CHKERRQ(DMGlobalToLocalBegin(da,global,INSERT_VALUES,local)); 123 CHKERRQ(DMGlobalToLocalEnd(da,global,INSERT_VALUES,local)); 124 125 flg = PETSC_FALSE; 126 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-local_print",&flg,NULL)); 127 if (flg) { 128 PetscViewer sviewer; 129 CHKERRQ(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 130 CHKERRQ(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"\nLocal Vector: processor %d\n",rank)); 131 CHKERRQ(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&sviewer)); 132 CHKERRQ(VecView(local,sviewer)); 133 CHKERRQ(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&sviewer)); 134 CHKERRQ(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 135 CHKERRQ(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 136 } 137 138 /* Tests mappings between application/PETSc orderings */ 139 if (test_order) { 140 ISLocalToGlobalMapping ltogm; 141 142 CHKERRQ(DMGetLocalToGlobalMapping(da,<ogm)); 143 CHKERRQ(ISLocalToGlobalMappingGetSize(ltogm,&nloc)); 144 CHKERRQ(ISLocalToGlobalMappingGetIndices(ltogm,<og)); 145 146 CHKERRQ(DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,&Xm,&Ym,&Zm)); 147 CHKERRQ(DMDAGetAO(da,&ao)); 148 /* CHKERRQ(AOView(ao,PETSC_VIEWER_STDOUT_WORLD)); */ 149 CHKERRQ(PetscMalloc1(nloc,&iglobal)); 150 151 /* Set iglobal to be global indices for each processor's local and ghost nodes, 152 using the DMDA ordering of grid points */ 153 kk = 0; 154 for (k=Zs; k<Zs+Zm; k++) { 155 for (j=Ys; j<Ys+Ym; j++) { 156 for (i=Xs; i<Xs+Xm; i++) { 157 iloc = w*((k-Zs)*Xm*Ym + (j-Ys)*Xm + i-Xs); 158 for (l=0; l<w; l++) { 159 iglobal[kk++] = ltog[iloc+l]; 160 } 161 } 162 } 163 } 164 165 /* Map this to the application ordering (which for DMDAs is just the natural ordering 166 that would be used for 1 processor, numbering most rapidly by x, then y, then z) */ 167 CHKERRQ(AOPetscToApplication(ao,nloc,iglobal)); 168 169 /* Then map the application ordering back to the PETSc DMDA ordering */ 170 CHKERRQ(AOApplicationToPetsc(ao,nloc,iglobal)); 171 172 /* Verify the mappings */ 173 kk=0; 174 for (k=Zs; k<Zs+Zm; k++) { 175 for (j=Ys; j<Ys+Ym; j++) { 176 for (i=Xs; i<Xs+Xm; i++) { 177 iloc = w*((k-Zs)*Xm*Ym + (j-Ys)*Xm + i-Xs); 178 for (l=0; l<w; l++) { 179 if (iglobal[kk] != ltog[iloc+l]) { 180 CHKERRQ(PetscPrintf(MPI_COMM_WORLD,"[%D] Problem with mapping: z=%D, j=%D, i=%D, l=%D, petsc1=%D, petsc2=%D\n",rank,k,j,i,l,ltog[iloc+l],iglobal[kk])); 181 } 182 kk++; 183 } 184 } 185 } 186 } 187 CHKERRQ(PetscFree(iglobal)); 188 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(ltogm,<og)); 189 } 190 191 /* Free memory */ 192 CHKERRQ(PetscViewerDestroy(&viewer)); 193 CHKERRQ(VecDestroy(&local)); 194 CHKERRQ(VecDestroy(&global)); 195 CHKERRQ(DMDestroy(&da)); 196 CHKERRQ(PetscFinalize()); 197 return 0; 198 } 199 200 /*TEST 201 202 test: 203 args: -testorder -nox 204 205 TEST*/ 206