1 2 static char help[] = "Tests DMCreateDomainDecomposition.\n\n"; 3 4 /* 5 Use the options 6 -da_grid_x <nx> - number of grid points in x direction, if M < 0 7 -da_grid_y <ny> - number of grid points in y direction, if N < 0 8 -da_processors_x <MX> number of processors in x directio 9 -da_processors_y <MY> number of processors in x direction 10 */ 11 12 #include <petscdm.h> 13 #include <petscdmda.h> 14 15 PetscErrorCode FillLocalSubdomain(DM da, Vec gvec) 16 { 17 DMDALocalInfo info; 18 PetscMPIInt rank; 19 PetscInt i,j,k,l; 20 21 PetscFunctionBeginUser; 22 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 23 PetscCall(DMDAGetLocalInfo(da,&info)); 24 25 if (info.dim == 3) { 26 PetscScalar ***g; 27 PetscCall(DMDAVecGetArray(da,gvec,&g)); 28 /* loop over ghosts */ 29 for (k=info.zs; k<info.zs+info.zm; k++) { 30 for (j=info.ys; j<info.ys+info.ym; j++) { 31 for (i=info.xs; i<info.xs+info.xm; i++) { 32 g[k][j][info.dof*i+0] = i; 33 g[k][j][info.dof*i+1] = j; 34 g[k][j][info.dof*i+2] = k; 35 } 36 } 37 } 38 PetscCall(DMDAVecRestoreArray(da,gvec,&g)); 39 } 40 if (info.dim == 2) { 41 PetscScalar **g; 42 PetscCall(DMDAVecGetArray(da,gvec,&g)); 43 /* loop over ghosts */ 44 for (j=info.ys; j<info.ys+info.ym; j++) { 45 for (i=info.xs; i<info.xs+info.xm; i++) { 46 for (l = 0;l<info.dof;l++) { 47 g[j][info.dof*i+0] = i; 48 g[j][info.dof*i+1] = j; 49 g[j][info.dof*i+2] = rank; 50 } 51 } 52 } 53 PetscCall(DMDAVecRestoreArray(da,gvec,&g)); 54 } 55 PetscFunctionReturn(0); 56 } 57 58 int main(int argc,char **argv) 59 { 60 DM da,*subda; 61 PetscInt i,dim = 3; 62 PetscInt M = 25, N = 25, P = 25; 63 PetscMPIInt size,rank; 64 Vec v; 65 Vec slvec,sgvec; 66 IS *ois,*iis; 67 VecScatter oscata; 68 VecScatter *iscat,*oscat,*gscat; 69 DMDALocalInfo info; 70 PetscBool patchis_offproc = PETSC_TRUE; 71 72 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 73 PetscCall(PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL)); 74 75 /* Create distributed array and get vectors */ 76 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 77 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 78 if (dim == 2) { 79 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da)); 80 } else if (dim == 3) { 81 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,NULL,&da)); 82 } 83 PetscCall(DMSetFromOptions(da)); 84 PetscCall(DMSetUp(da)); 85 PetscCall(DMDAGetLocalInfo(da,&info)); 86 87 PetscCall(DMCreateDomainDecomposition(da,NULL,NULL,&iis,&ois,&subda)); 88 PetscCall(DMCreateDomainDecompositionScatters(da,1,subda,&iscat,&oscat,&gscat)); 89 90 { 91 DMDALocalInfo subinfo; 92 MatStencil lower,upper; 93 IS patchis; 94 Vec smallvec; 95 Vec largevec; 96 VecScatter patchscat; 97 98 PetscCall(DMDAGetLocalInfo(subda[0],&subinfo)); 99 100 lower.i = info.xs; 101 lower.j = info.ys; 102 lower.k = info.zs; 103 upper.i = info.xs+info.xm; 104 upper.j = info.ys+info.ym; 105 upper.k = info.zs+info.zm; 106 107 /* test the patch IS as a thing to scatter to/from */ 108 PetscCall(DMDACreatePatchIS(da,&lower,&upper,&patchis,patchis_offproc)); 109 PetscCall(DMGetGlobalVector(da,&largevec)); 110 111 PetscCall(VecCreate(PETSC_COMM_SELF,&smallvec)); 112 PetscCall(VecSetSizes(smallvec,info.dof*(upper.i - lower.i)*(upper.j - lower.j)*(upper.k - lower.k),PETSC_DECIDE)); 113 PetscCall(VecSetFromOptions(smallvec)); 114 PetscCall(VecScatterCreate(smallvec,NULL,largevec,patchis,&patchscat)); 115 116 PetscCall(FillLocalSubdomain(subda[0],smallvec)); 117 PetscCall(VecSet(largevec,0)); 118 119 PetscCall(VecScatterBegin(patchscat,smallvec,largevec,ADD_VALUES,SCATTER_FORWARD)); 120 PetscCall(VecScatterEnd(patchscat,smallvec,largevec,ADD_VALUES,SCATTER_FORWARD)); 121 PetscCall(ISView(patchis,PETSC_VIEWER_STDOUT_WORLD)); 122 PetscCall(VecScatterView(patchscat,PETSC_VIEWER_STDOUT_WORLD)); 123 124 for (i = 0; i < size; i++) { 125 if (i == rank) { 126 PetscCall(VecView(smallvec,PETSC_VIEWER_STDOUT_SELF)); 127 } 128 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 129 } 130 131 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 132 PetscCall(VecView(largevec,PETSC_VIEWER_STDOUT_WORLD)); 133 134 PetscCall(VecDestroy(&smallvec)); 135 PetscCall(DMRestoreGlobalVector(da,&largevec)); 136 PetscCall(ISDestroy(&patchis)); 137 PetscCall(VecScatterDestroy(&patchscat)); 138 } 139 140 /* view the various parts */ 141 { 142 for (i = 0; i < size; i++) { 143 if (i == rank) { 144 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Processor %d: \n",i)); 145 PetscCall(DMView(subda[0],PETSC_VIEWER_STDOUT_SELF)); 146 } 147 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 148 } 149 150 PetscCall(DMGetLocalVector(subda[0],&slvec)); 151 PetscCall(DMGetGlobalVector(subda[0],&sgvec)); 152 PetscCall(DMGetGlobalVector(da,&v)); 153 154 /* test filling outer between the big DM and the small ones with the IS scatter*/ 155 PetscCall(VecScatterCreate(v,ois[0],sgvec,NULL,&oscata)); 156 157 PetscCall(FillLocalSubdomain(subda[0],sgvec)); 158 159 PetscCall(VecScatterBegin(oscata,sgvec,v,ADD_VALUES,SCATTER_REVERSE)); 160 PetscCall(VecScatterEnd(oscata,sgvec,v,ADD_VALUES,SCATTER_REVERSE)); 161 162 /* test the local-to-local scatter */ 163 164 /* fill up the local subdomain and then add them together */ 165 PetscCall(FillLocalSubdomain(da,v)); 166 167 PetscCall(VecScatterBegin(gscat[0],v,slvec,ADD_VALUES,SCATTER_FORWARD)); 168 PetscCall(VecScatterEnd(gscat[0],v,slvec,ADD_VALUES,SCATTER_FORWARD)); 169 170 PetscCall(VecView(v,PETSC_VIEWER_STDOUT_WORLD)); 171 172 /* test ghost scattering backwards */ 173 174 PetscCall(VecSet(v,0)); 175 176 PetscCall(VecScatterBegin(gscat[0],slvec,v,ADD_VALUES,SCATTER_REVERSE)); 177 PetscCall(VecScatterEnd(gscat[0],slvec,v,ADD_VALUES,SCATTER_REVERSE)); 178 179 PetscCall(VecView(v,PETSC_VIEWER_STDOUT_WORLD)); 180 181 /* test overlap scattering backwards */ 182 183 PetscCall(DMLocalToGlobalBegin(subda[0],slvec,ADD_VALUES,sgvec)); 184 PetscCall(DMLocalToGlobalEnd(subda[0],slvec,ADD_VALUES,sgvec)); 185 186 PetscCall(VecSet(v,0)); 187 188 PetscCall(VecScatterBegin(oscat[0],sgvec,v,ADD_VALUES,SCATTER_REVERSE)); 189 PetscCall(VecScatterEnd(oscat[0],sgvec,v,ADD_VALUES,SCATTER_REVERSE)); 190 191 PetscCall(VecView(v,PETSC_VIEWER_STDOUT_WORLD)); 192 193 /* test interior scattering backwards */ 194 195 PetscCall(VecSet(v,0)); 196 197 PetscCall(VecScatterBegin(iscat[0],sgvec,v,ADD_VALUES,SCATTER_REVERSE)); 198 PetscCall(VecScatterEnd(iscat[0],sgvec,v,ADD_VALUES,SCATTER_REVERSE)); 199 200 PetscCall(VecView(v,PETSC_VIEWER_STDOUT_WORLD)); 201 202 /* test matrix allocation */ 203 for (i = 0; i < size; i++) { 204 if (i == rank) { 205 Mat m; 206 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Processor %d: \n",i)); 207 PetscCall(DMSetMatType(subda[0],MATAIJ)); 208 PetscCall(DMCreateMatrix(subda[0],&m)); 209 PetscCall(MatView(m,PETSC_VIEWER_STDOUT_SELF)); 210 PetscCall(MatDestroy(&m)); 211 } 212 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 213 } 214 PetscCall(DMRestoreLocalVector(subda[0],&slvec)); 215 PetscCall(DMRestoreGlobalVector(subda[0],&sgvec)); 216 PetscCall(DMRestoreGlobalVector(da,&v)); 217 } 218 219 PetscCall(DMDestroy(&subda[0])); 220 PetscCall(ISDestroy(&ois[0])); 221 PetscCall(ISDestroy(&iis[0])); 222 223 PetscCall(VecScatterDestroy(&iscat[0])); 224 PetscCall(VecScatterDestroy(&oscat[0])); 225 PetscCall(VecScatterDestroy(&gscat[0])); 226 PetscCall(VecScatterDestroy(&oscata)); 227 228 PetscCall(PetscFree(iscat)); 229 PetscCall(PetscFree(oscat)); 230 PetscCall(PetscFree(gscat)); 231 PetscCall(PetscFree(oscata)); 232 233 PetscCall(PetscFree(subda)); 234 PetscCall(PetscFree(ois)); 235 PetscCall(PetscFree(iis)); 236 237 PetscCall(DMDestroy(&da)); 238 PetscCall(PetscFinalize()); 239 return 0; 240 } 241