1 #include <petsc/private/dmpatchimpl.h> /*I "petscdmpatch.h" I*/ 2 #include <petscdmda.h> 3 #include <petscsf.h> 4 5 /* 6 Solver loop to update \tau: 7 8 DMZoom(dmc, &dmz) 9 DMRefine(dmz, &dmf), 10 Scatter Xcoarse -> Xzoom, 11 Interpolate Xzoom -> Xfine (note that this may be on subcomms), 12 Smooth Xfine using two-step smoother 13 normal smoother plus Kaczmarz---moves back and forth from dmzoom to dmfine 14 Compute residual Rfine 15 Restrict Rfine to Rzoom_restricted 16 Scatter Rzoom_restricted -> Rcoarse_restricted 17 Compute global residual Rcoarse 18 TauCoarse = Rcoarse - Rcoarse_restricted 19 */ 20 21 /*@C 22 DMPatchZoom - Create patches of a DMDA on subsets of processes, indicated by commz 23 24 Collective on dm 25 26 Input Parameters: 27 + dm - the DM 28 . lower,upper - the upper right corner and the lower left corner of the requested patch 29 - commz - the new communicator for the patch, MPI_COMM_NULL indicates that the given rank will not own a patch 30 31 Output Parameters: 32 + dmz - the patch DM 33 . sfz - the PetscSF mapping the patch+halo to the zoomed version (optional) 34 - sfzr - the PetscSF mapping the patch to the restricted zoomed version 35 36 Level: intermediate 37 38 .seealso: DMPatchSolve(), DMDACreatePatchIS() 39 @*/ 40 PetscErrorCode DMPatchZoom(DM dm, MatStencil lower, MatStencil upper, MPI_Comm commz, DM *dmz, PetscSF *sfz, PetscSF *sfzr) 41 { 42 DMDAStencilType st; 43 MatStencil blower, bupper, loclower, locupper; 44 IS is; 45 const PetscInt *ranges, *indices; 46 PetscInt *localPoints = NULL; 47 PetscSFNode *remotePoints = NULL; 48 PetscInt dim, dof; 49 PetscInt M, N, P, rM, rN, rP, halo = 1, sxb, syb, szb, sxr, syr, szr, exr, eyr, ezr, mxb, myb, mzb, i, j, k, l, q; 50 PetscMPIInt size; 51 PetscBool patchis_offproc = PETSC_TRUE; 52 PetscErrorCode ierr; 53 Vec X; 54 55 PetscFunctionBegin; 56 if (!sfz) halo = 0; 57 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);CHKERRMPI(ierr); 58 /* Create patch DM */ 59 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, NULL,NULL,NULL, &dof, NULL,NULL,NULL,NULL, &st);CHKERRQ(ierr); 60 61 /* Get piece for rank r, expanded by halo */ 62 bupper.i = PetscMin(M, upper.i + halo); blower.i = PetscMax(lower.i - halo, 0); 63 bupper.j = PetscMin(N, upper.j + halo); blower.j = PetscMax(lower.j - halo, 0); 64 bupper.k = PetscMin(P, upper.k + halo); blower.k = PetscMax(lower.k - halo, 0); 65 rM = bupper.i - blower.i; 66 rN = bupper.j - blower.j; 67 rP = bupper.k - blower.k; 68 69 if (commz != MPI_COMM_NULL) { 70 ierr = DMDACreate(commz, dmz);CHKERRQ(ierr); 71 ierr = DMSetDimension(*dmz, dim);CHKERRQ(ierr); 72 ierr = DMDASetSizes(*dmz, rM, rN, rP);CHKERRQ(ierr); 73 ierr = DMDASetNumProcs(*dmz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);CHKERRQ(ierr); 74 ierr = DMDASetBoundaryType(*dmz, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE);CHKERRQ(ierr); 75 ierr = DMDASetDof(*dmz, dof);CHKERRQ(ierr); 76 ierr = DMDASetStencilType(*dmz, st);CHKERRQ(ierr); 77 ierr = DMDASetStencilWidth(*dmz, 0);CHKERRQ(ierr); 78 ierr = DMDASetOwnershipRanges(*dmz, NULL, NULL, NULL);CHKERRQ(ierr); 79 ierr = DMSetFromOptions(*dmz);CHKERRQ(ierr); 80 ierr = DMSetUp(*dmz);CHKERRQ(ierr); 81 ierr = DMDAGetCorners(*dmz, &sxb, &syb, &szb, &mxb, &myb, &mzb);CHKERRQ(ierr); 82 sxr = PetscMax(sxb, lower.i - blower.i); 83 syr = PetscMax(syb, lower.j - blower.j); 84 szr = PetscMax(szb, lower.k - blower.k); 85 exr = PetscMin(sxb+mxb, upper.i - blower.i); 86 eyr = PetscMin(syb+myb, upper.j - blower.j); 87 ezr = PetscMin(szb+mzb, upper.k - blower.k); 88 ierr = PetscMalloc2(dof*rM*rN*PetscMax(rP,1),&localPoints,dof*rM*rN*PetscMax(rP,1),&remotePoints);CHKERRQ(ierr); 89 } else { 90 sxr = syr = szr = exr = eyr = ezr = sxb = syb = szb = mxb = myb = mzb = 0; 91 } 92 93 /* Create SF for restricted map */ 94 ierr = DMCreateGlobalVector(dm,&X);CHKERRQ(ierr); 95 ierr = VecGetOwnershipRanges(X,&ranges);CHKERRQ(ierr); 96 97 loclower.i = blower.i + sxr; locupper.i = blower.i + exr; 98 loclower.j = blower.j + syr; locupper.j = blower.j + eyr; 99 loclower.k = blower.k + szr; locupper.k = blower.k + ezr; 100 101 ierr = DMDACreatePatchIS(dm, &loclower, &locupper, &is, patchis_offproc);CHKERRQ(ierr); 102 ierr = ISGetIndices(is, &indices);CHKERRQ(ierr); 103 104 if (dim < 3) {mzb = 1; ezr = 1;} 105 q = 0; 106 for (k = szb; k < szb+mzb; ++k) { 107 if ((k < szr) || (k >= ezr)) continue; 108 for (j = syb; j < syb+myb; ++j) { 109 if ((j < syr) || (j >= eyr)) continue; 110 for (i = sxb; i < sxb+mxb; ++i) { 111 for (l=0; l<dof; l++) { 112 const PetscInt lp = l + dof*(((k-szb)*rN + (j-syb))*rM + i-sxb); 113 PetscInt r; 114 115 if ((i < sxr) || (i >= exr)) continue; 116 localPoints[q] = lp; 117 ierr = PetscFindInt(indices[q], size+1, ranges, &r);CHKERRQ(ierr); 118 119 remotePoints[q].rank = r < 0 ? -(r+1) - 1 : r; 120 remotePoints[q].index = indices[q] - ranges[remotePoints[q].rank]; 121 ++q; 122 } 123 } 124 } 125 } 126 ierr = ISRestoreIndices(is, &indices);CHKERRQ(ierr); 127 ierr = ISDestroy(&is);CHKERRQ(ierr); 128 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfzr);CHKERRQ(ierr); 129 ierr = PetscObjectSetName((PetscObject) *sfzr, "Restricted Map");CHKERRQ(ierr); 130 ierr = PetscSFSetGraph(*sfzr, dof*M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);CHKERRQ(ierr); 131 132 if (sfz) { 133 /* Create SF for buffered map */ 134 loclower.i = blower.i + sxb; locupper.i = blower.i + sxb+mxb; 135 loclower.j = blower.j + syb; locupper.j = blower.j + syb+myb; 136 loclower.k = blower.k + szb; locupper.k = blower.k + szb+mzb; 137 138 ierr = DMDACreatePatchIS(dm, &loclower, &locupper, &is, patchis_offproc);CHKERRQ(ierr); 139 ierr = ISGetIndices(is, &indices);CHKERRQ(ierr); 140 141 q = 0; 142 for (k = szb; k < szb+mzb; ++k) { 143 for (j = syb; j < syb+myb; ++j) { 144 for (i = sxb; i < sxb+mxb; ++i, ++q) { 145 PetscInt r; 146 147 localPoints[q] = q; 148 ierr = PetscFindInt(indices[q], size+1, ranges, &r);CHKERRQ(ierr); 149 remotePoints[q].rank = r < 0 ? -(r+1) - 1 : r; 150 remotePoints[q].index = indices[q] - ranges[remotePoints[q].rank]; 151 } 152 } 153 } 154 ierr = ISRestoreIndices(is, &indices);CHKERRQ(ierr); 155 ierr = ISDestroy(&is);CHKERRQ(ierr); 156 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfz);CHKERRQ(ierr); 157 ierr = PetscObjectSetName((PetscObject) *sfz, "Buffered Map");CHKERRQ(ierr); 158 ierr = PetscSFSetGraph(*sfz, M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);CHKERRQ(ierr); 159 } 160 161 ierr = VecDestroy(&X);CHKERRQ(ierr); 162 ierr = PetscFree2(localPoints, remotePoints);CHKERRQ(ierr); 163 PetscFunctionReturn(0); 164 } 165 166 typedef enum {PATCH_COMM_TYPE_WORLD = 0, PATCH_COMM_TYPE_SELF = 1} PatchCommType; 167 168 PetscErrorCode DMPatchSolve(DM dm) 169 { 170 MPI_Comm comm; 171 MPI_Comm commz; 172 DM dmc; 173 PetscSF sfz, sfzr; 174 Vec XC; 175 MatStencil patchSize, commSize, gridRank, lower, upper; 176 PetscInt M, N, P, i, j, k, l, m, n, p = 0; 177 PetscMPIInt rank, size; 178 PetscInt debug = 0; 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 183 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 184 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 185 ierr = DMPatchGetCoarse(dm, &dmc);CHKERRQ(ierr); 186 ierr = DMPatchGetPatchSize(dm, &patchSize);CHKERRQ(ierr); 187 ierr = DMPatchGetCommSize(dm, &commSize);CHKERRQ(ierr); 188 ierr = DMPatchGetCommSize(dm, &commSize);CHKERRQ(ierr); 189 ierr = DMGetGlobalVector(dmc, &XC);CHKERRQ(ierr); 190 ierr = DMDAGetInfo(dmc, NULL, &M, &N, &P, &l, &m, &n, NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 191 M = PetscMax(M, 1); l = PetscMax(l, 1); 192 N = PetscMax(N, 1); m = PetscMax(m, 1); 193 P = PetscMax(P, 1); n = PetscMax(n, 1); 194 195 gridRank.i = rank % l; 196 gridRank.j = rank/l % m; 197 gridRank.k = rank/(l*m) % n; 198 199 if (commSize.i*commSize.j*commSize.k == size || commSize.i*commSize.j*commSize.k == 0) { 200 commSize.i = l; commSize.j = m; commSize.k = n; 201 commz = comm; 202 } else if (commSize.i*commSize.j*commSize.k == 1) { 203 commz = PETSC_COMM_SELF; 204 } else { 205 const PetscMPIInt newComm = ((gridRank.k/commSize.k)*(m/commSize.j) + gridRank.j/commSize.j)*(l/commSize.i) + (gridRank.i/commSize.i); 206 const PetscMPIInt newRank = ((gridRank.k%commSize.k)*commSize.j + gridRank.j%commSize.j)*commSize.i + (gridRank.i%commSize.i); 207 208 ierr = MPI_Comm_split(comm, newComm, newRank, &commz);CHKERRMPI(ierr); 209 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "Rank %d color %d key %d commz %d\n", rank, newComm, newRank, *((PetscMPIInt*) &commz));CHKERRQ(ierr);} 210 } 211 /* 212 Assumptions: 213 - patchSize divides gridSize 214 - commSize divides gridSize 215 - commSize divides l,m,n 216 Ignore multiple patches per rank for now 217 218 Multiple ranks per patch: 219 - l,m,n divides patchSize 220 - commSize divides patchSize 221 */ 222 for (k = 0; k < P; k += PetscMax(patchSize.k, 1)) { 223 for (j = 0; j < N; j += PetscMax(patchSize.j, 1)) { 224 for (i = 0; i < M; i += PetscMax(patchSize.i, 1), ++p) { 225 MPI_Comm commp = MPI_COMM_NULL; 226 DM dmz = NULL; 227 #if 0 228 DM dmf = NULL; 229 Mat interpz = NULL; 230 #endif 231 Vec XZ = NULL; 232 PetscScalar *xcarray = NULL; 233 PetscScalar *xzarray = NULL; 234 235 if ((gridRank.k/commSize.k == p/(l/commSize.i * m/commSize.j) % n/commSize.k) && 236 (gridRank.j/commSize.j == p/(l/commSize.i) % m/commSize.j) && 237 (gridRank.i/commSize.i == p % l/commSize.i)) { 238 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "Rank %d is accepting Patch %d\n", rank, p);CHKERRQ(ierr);} 239 commp = commz; 240 } 241 /* Zoom to coarse patch */ 242 lower.i = i; lower.j = j; lower.k = k; 243 upper.i = i + patchSize.i; upper.j = j + patchSize.j; upper.k = k + patchSize.k; 244 ierr = DMPatchZoom(dmc, lower, upper, commp, &dmz, &sfz, &sfzr);CHKERRQ(ierr); 245 lower.c = 0; /* initialize member, otherwise compiler issues warnings */ 246 upper.c = 0; /* initialize member, otherwise compiler issues warnings */ 247 if (debug) {ierr = PetscPrintf(comm, "Patch %d: (%d, %d, %d)--(%d, %d, %d)\n", p, lower.i, lower.j, lower.k, upper.i, upper.j, upper.k);CHKERRQ(ierr);} 248 if (dmz) {ierr = DMView(dmz, PETSC_VIEWER_STDOUT_(commz));CHKERRQ(ierr);} 249 ierr = PetscSFView(sfz, PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr); 250 ierr = PetscSFView(sfzr, PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr); 251 /* Scatter Xcoarse -> Xzoom */ 252 if (dmz) {ierr = DMGetGlobalVector(dmz, &XZ);CHKERRQ(ierr);} 253 if (XZ) {ierr = VecGetArray(XZ, &xzarray);CHKERRQ(ierr);} 254 ierr = VecGetArray(XC, &xcarray);CHKERRQ(ierr); 255 ierr = PetscSFBcastBegin(sfz, MPIU_SCALAR, xcarray, xzarray,MPI_REPLACE);CHKERRQ(ierr); 256 ierr = PetscSFBcastEnd(sfz, MPIU_SCALAR, xcarray, xzarray,MPI_REPLACE);CHKERRQ(ierr); 257 ierr = VecRestoreArray(XC, &xcarray);CHKERRQ(ierr); 258 if (XZ) {ierr = VecRestoreArray(XZ, &xzarray);CHKERRQ(ierr);} 259 #if 0 260 /* Interpolate Xzoom -> Xfine, note that this may be on subcomms */ 261 ierr = DMRefine(dmz, MPI_COMM_NULL, &dmf);CHKERRQ(ierr); 262 ierr = DMCreateInterpolation(dmz, dmf, &interpz, NULL);CHKERRQ(ierr); 263 ierr = DMInterpolate(dmz, interpz, dmf);CHKERRQ(ierr); 264 /* Smooth Xfine using two-step smoother, normal smoother plus Kaczmarz---moves back and forth from dmzoom to dmfine */ 265 /* Compute residual Rfine */ 266 /* Restrict Rfine to Rzoom_restricted */ 267 #endif 268 /* Scatter Rzoom_restricted -> Rcoarse_restricted */ 269 if (XZ) {ierr = VecGetArray(XZ, &xzarray);CHKERRQ(ierr);} 270 ierr = VecGetArray(XC, &xcarray);CHKERRQ(ierr); 271 ierr = PetscSFReduceBegin(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM);CHKERRQ(ierr); 272 ierr = PetscSFReduceEnd(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM);CHKERRQ(ierr); 273 ierr = VecRestoreArray(XC, &xcarray);CHKERRQ(ierr); 274 if (XZ) {ierr = VecRestoreArray(XZ, &xzarray);CHKERRQ(ierr);} 275 if (dmz) {ierr = DMRestoreGlobalVector(dmz, &XZ);CHKERRQ(ierr);} 276 /* Compute global residual Rcoarse */ 277 /* TauCoarse = Rcoarse - Rcoarse_restricted */ 278 279 ierr = PetscSFDestroy(&sfz);CHKERRQ(ierr); 280 ierr = PetscSFDestroy(&sfzr);CHKERRQ(ierr); 281 ierr = DMDestroy(&dmz);CHKERRQ(ierr); 282 } 283 } 284 } 285 ierr = DMRestoreGlobalVector(dmc, &XC);CHKERRQ(ierr); 286 PetscFunctionReturn(0); 287 } 288 289 PetscErrorCode DMPatchView_ASCII(DM dm, PetscViewer viewer) 290 { 291 DM_Patch *mesh = (DM_Patch*) dm->data; 292 PetscViewerFormat format; 293 const char *name; 294 PetscErrorCode ierr; 295 296 PetscFunctionBegin; 297 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 298 /* if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) */ 299 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 300 ierr = PetscViewerASCIIPrintf(viewer, "Patch DM %s\n", name);CHKERRQ(ierr); 301 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 302 ierr = PetscViewerASCIIPrintf(viewer, "Coarse DM\n");CHKERRQ(ierr); 303 ierr = DMView(mesh->dmCoarse, viewer);CHKERRQ(ierr); 304 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 305 PetscFunctionReturn(0); 306 } 307 308 PetscErrorCode DMView_Patch(DM dm, PetscViewer viewer) 309 { 310 PetscBool iascii, isbinary; 311 PetscErrorCode ierr; 312 313 PetscFunctionBegin; 314 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 315 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 316 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 317 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);CHKERRQ(ierr); 318 if (iascii) { 319 ierr = DMPatchView_ASCII(dm, viewer);CHKERRQ(ierr); 320 } 321 PetscFunctionReturn(0); 322 } 323 324 PetscErrorCode DMDestroy_Patch(DM dm) 325 { 326 DM_Patch *mesh = (DM_Patch*) dm->data; 327 PetscErrorCode ierr; 328 329 PetscFunctionBegin; 330 if (--mesh->refct > 0) PetscFunctionReturn(0); 331 ierr = DMDestroy(&mesh->dmCoarse);CHKERRQ(ierr); 332 /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */ 333 ierr = PetscFree(mesh);CHKERRQ(ierr); 334 PetscFunctionReturn(0); 335 } 336 337 PetscErrorCode DMSetUp_Patch(DM dm) 338 { 339 DM_Patch *mesh = (DM_Patch*) dm->data; 340 PetscErrorCode ierr; 341 342 PetscFunctionBegin; 343 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 344 ierr = DMSetUp(mesh->dmCoarse);CHKERRQ(ierr); 345 PetscFunctionReturn(0); 346 } 347 348 PetscErrorCode DMCreateGlobalVector_Patch(DM dm, Vec *g) 349 { 350 DM_Patch *mesh = (DM_Patch*) dm->data; 351 PetscErrorCode ierr; 352 353 PetscFunctionBegin; 354 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 355 ierr = DMCreateGlobalVector(mesh->dmCoarse, g);CHKERRQ(ierr); 356 PetscFunctionReturn(0); 357 } 358 359 PetscErrorCode DMCreateLocalVector_Patch(DM dm, Vec *l) 360 { 361 DM_Patch *mesh = (DM_Patch*) dm->data; 362 PetscErrorCode ierr; 363 364 PetscFunctionBegin; 365 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 366 ierr = DMCreateLocalVector(mesh->dmCoarse, l);CHKERRQ(ierr); 367 PetscFunctionReturn(0); 368 } 369 370 PetscErrorCode DMCreateSubDM_Patch(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm) 371 { 372 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tell me to code this"); 373 } 374 375 PetscErrorCode DMPatchGetCoarse(DM dm, DM *dmCoarse) 376 { 377 DM_Patch *mesh = (DM_Patch*) dm->data; 378 379 PetscFunctionBegin; 380 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 381 *dmCoarse = mesh->dmCoarse; 382 PetscFunctionReturn(0); 383 } 384 385 PetscErrorCode DMPatchGetPatchSize(DM dm, MatStencil *patchSize) 386 { 387 DM_Patch *mesh = (DM_Patch*) dm->data; 388 389 PetscFunctionBegin; 390 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 391 PetscValidPointer(patchSize, 2); 392 *patchSize = mesh->patchSize; 393 PetscFunctionReturn(0); 394 } 395 396 PetscErrorCode DMPatchSetPatchSize(DM dm, MatStencil patchSize) 397 { 398 DM_Patch *mesh = (DM_Patch*) dm->data; 399 400 PetscFunctionBegin; 401 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 402 mesh->patchSize = patchSize; 403 PetscFunctionReturn(0); 404 } 405 406 PetscErrorCode DMPatchGetCommSize(DM dm, MatStencil *commSize) 407 { 408 DM_Patch *mesh = (DM_Patch*) dm->data; 409 410 PetscFunctionBegin; 411 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 412 PetscValidPointer(commSize, 2); 413 *commSize = mesh->commSize; 414 PetscFunctionReturn(0); 415 } 416 417 PetscErrorCode DMPatchSetCommSize(DM dm, MatStencil commSize) 418 { 419 DM_Patch *mesh = (DM_Patch*) dm->data; 420 421 PetscFunctionBegin; 422 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 423 mesh->commSize = commSize; 424 PetscFunctionReturn(0); 425 } 426