1 static char help[] = "Tests DMCreateDomainDecomposition.\n\n"; 2 3 /* 4 Use the options 5 -da_grid_x <nx> - number of grid points in x direction, if M < 0 6 -da_grid_y <ny> - number of grid points in y direction, if N < 0 7 -da_processors_x <MX> number of processors in x directio 8 -da_processors_y <MY> number of processors in x direction 9 */ 10 11 #include <petscdm.h> 12 #include <petscdmda.h> 13 14 PetscErrorCode FillLocalSubdomain(DM da, Vec gvec) 15 { 16 DMDALocalInfo info; 17 PetscMPIInt rank; 18 PetscInt i, j, k, l; 19 20 PetscFunctionBeginUser; 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 22 PetscCall(DMDAGetLocalInfo(da, &info)); 23 24 if (info.dim == 3) { 25 PetscScalar ***g; 26 PetscCall(DMDAVecGetArray(da, gvec, &g)); 27 /* loop over ghosts */ 28 for (k = info.zs; k < info.zs + info.zm; k++) { 29 for (j = info.ys; j < info.ys + info.ym; j++) { 30 for (i = info.xs; i < info.xs + info.xm; i++) { 31 g[k][j][info.dof * i + 0] = i; 32 g[k][j][info.dof * i + 1] = j; 33 g[k][j][info.dof * i + 2] = k; 34 } 35 } 36 } 37 PetscCall(DMDAVecRestoreArray(da, gvec, &g)); 38 } 39 if (info.dim == 2) { 40 PetscScalar **g; 41 PetscCall(DMDAVecGetArray(da, gvec, &g)); 42 /* loop over ghosts */ 43 for (j = info.ys; j < info.ys + info.ym; j++) { 44 for (i = info.xs; i < info.xs + info.xm; i++) { 45 for (l = 0; l < info.dof; l++) { 46 g[j][info.dof * i + 0] = i; 47 g[j][info.dof * i + 1] = j; 48 g[j][info.dof * i + 2] = rank; 49 } 50 } 51 } 52 PetscCall(DMDAVecRestoreArray(da, gvec, &g)); 53 } 54 PetscFunctionReturn(PETSC_SUCCESS); 55 } 56 57 int main(int argc, char **argv) 58 { 59 DM da, *subda; 60 PetscInt i, dim = 3, dof = 3; 61 PetscInt M = 25, N = 25, P = 25; 62 PetscMPIInt size, rank; 63 Vec v; 64 Vec slvec, sgvec; 65 IS *ois, *iis; 66 VecScatter oscata; 67 VecScatter *iscat, *oscat, *gscat; 68 DMDALocalInfo info; 69 PetscBool patchis_offproc = PETSC_TRUE; 70 71 PetscFunctionBeginUser; 72 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 73 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL)); 74 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL)); 75 76 /* Create distributed array and get vectors */ 77 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 78 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 79 if (dim == 2) { 80 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, &da)); 81 } else if (dim == 3) { 82 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, dof, 1, NULL, NULL, NULL, &da)); 83 } 84 PetscCall(DMSetFromOptions(da)); 85 PetscCall(DMSetUp(da)); 86 PetscCall(DMDAGetLocalInfo(da, &info)); 87 88 PetscCall(DMCreateDomainDecomposition(da, NULL, NULL, &iis, &ois, &subda)); 89 PetscCall(DMCreateDomainDecompositionScatters(da, 1, subda, &iscat, &oscat, &gscat)); 90 91 /* TODO: broken for dim=2 */ 92 if (dim == 3) { 93 DMDALocalInfo subinfo; 94 MatStencil lower, upper; 95 IS patchis; 96 Vec smallvec; 97 Vec largevec; 98 VecScatter patchscat; 99 100 PetscCall(DMDAGetLocalInfo(subda[0], &subinfo)); 101 102 lower.i = info.xs; 103 lower.j = info.ys; 104 lower.k = info.zs; 105 upper.i = info.xs + info.xm; 106 upper.j = info.ys + info.ym; 107 upper.k = info.zs + info.zm; 108 109 /* test the patch IS as a thing to scatter to/from */ 110 PetscCall(DMDACreatePatchIS(da, &lower, &upper, &patchis, patchis_offproc)); 111 PetscCall(DMGetGlobalVector(da, &largevec)); 112 113 PetscCall(VecCreate(PETSC_COMM_SELF, &smallvec)); 114 PetscCall(VecSetSizes(smallvec, info.dof * (upper.i - lower.i) * (upper.j - lower.j) * (upper.k - lower.k), PETSC_DECIDE)); 115 PetscCall(VecSetFromOptions(smallvec)); 116 PetscCall(VecScatterCreate(smallvec, NULL, largevec, patchis, &patchscat)); 117 118 PetscCall(FillLocalSubdomain(subda[0], smallvec)); 119 PetscCall(VecSet(largevec, 0)); 120 121 PetscCall(VecScatterBegin(patchscat, smallvec, largevec, ADD_VALUES, SCATTER_FORWARD)); 122 PetscCall(VecScatterEnd(patchscat, smallvec, largevec, ADD_VALUES, SCATTER_FORWARD)); 123 PetscCall(ISView(patchis, PETSC_VIEWER_STDOUT_WORLD)); 124 PetscCall(VecScatterView(patchscat, PETSC_VIEWER_STDOUT_WORLD)); 125 126 for (i = 0; i < size; i++) { 127 if (i == rank) PetscCall(VecView(smallvec, PETSC_VIEWER_STDOUT_SELF)); 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