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 - the lower left corner of the requested patch 29 . upper - the upper right corner of the requested patch 30 - commz - the new communicator for the patch, MPI_COMM_NULL indicates that the given rank will not own a patch 31 32 Output Parameters: 33 + dmz - the patch DM 34 . sfz - the PetscSF mapping the patch+halo to the zoomed version (optional) 35 - sfzr - the PetscSF mapping the patch to the restricted zoomed version 36 37 Level: intermediate 38 39 .seealso: `DMPatchSolve()`, `DMDACreatePatchIS()` 40 @*/ 41 PetscErrorCode DMPatchZoom(DM dm, MatStencil lower, MatStencil upper, MPI_Comm commz, DM *dmz, PetscSF *sfz, PetscSF *sfzr) 42 { 43 DMDAStencilType st; 44 MatStencil blower, bupper, loclower, locupper; 45 IS is; 46 const PetscInt *ranges, *indices; 47 PetscInt *localPoints = NULL; 48 PetscSFNode *remotePoints = NULL; 49 PetscInt dim, dof; 50 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; 51 PetscMPIInt size; 52 PetscBool patchis_offproc = PETSC_TRUE; 53 Vec X; 54 55 PetscFunctionBegin; 56 if (!sfz) halo = 0; 57 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size)); 58 /* Create patch DM */ 59 PetscCall(DMDAGetInfo(dm, &dim, &M, &N, &P, NULL,NULL,NULL, &dof, NULL,NULL,NULL,NULL, &st)); 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 PetscCall(DMDACreate(commz, dmz)); 71 PetscCall(DMSetDimension(*dmz, dim)); 72 PetscCall(DMDASetSizes(*dmz, rM, rN, rP)); 73 PetscCall(DMDASetNumProcs(*dmz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE)); 74 PetscCall(DMDASetBoundaryType(*dmz, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE)); 75 PetscCall(DMDASetDof(*dmz, dof)); 76 PetscCall(DMDASetStencilType(*dmz, st)); 77 PetscCall(DMDASetStencilWidth(*dmz, 0)); 78 PetscCall(DMDASetOwnershipRanges(*dmz, NULL, NULL, NULL)); 79 PetscCall(DMSetFromOptions(*dmz)); 80 PetscCall(DMSetUp(*dmz)); 81 PetscCall(DMDAGetCorners(*dmz, &sxb, &syb, &szb, &mxb, &myb, &mzb)); 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 PetscCall(PetscMalloc2(dof*rM*rN*PetscMax(rP,1),&localPoints,dof*rM*rN*PetscMax(rP,1),&remotePoints)); 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 PetscCall(DMCreateGlobalVector(dm,&X)); 95 PetscCall(VecGetOwnershipRanges(X,&ranges)); 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 PetscCall(DMDACreatePatchIS(dm, &loclower, &locupper, &is, patchis_offproc)); 102 PetscCall(ISGetIndices(is, &indices)); 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 PetscCall(PetscFindInt(indices[q], size+1, ranges, &r)); 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 PetscCall(ISRestoreIndices(is, &indices)); 127 PetscCall(ISDestroy(&is)); 128 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), sfzr)); 129 PetscCall(PetscObjectSetName((PetscObject) *sfzr, "Restricted Map")); 130 PetscCall(PetscSFSetGraph(*sfzr, dof*M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES)); 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 PetscCall(DMDACreatePatchIS(dm, &loclower, &locupper, &is, patchis_offproc)); 139 PetscCall(ISGetIndices(is, &indices)); 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 PetscCall(PetscFindInt(indices[q], size+1, ranges, &r)); 149 remotePoints[q].rank = r < 0 ? -(r+1) - 1 : r; 150 remotePoints[q].index = indices[q] - ranges[remotePoints[q].rank]; 151 } 152 } 153 } 154 PetscCall(ISRestoreIndices(is, &indices)); 155 PetscCall(ISDestroy(&is)); 156 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), sfz)); 157 PetscCall(PetscObjectSetName((PetscObject) *sfz, "Buffered Map")); 158 PetscCall(PetscSFSetGraph(*sfz, M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES)); 159 } 160 161 PetscCall(VecDestroy(&X)); 162 PetscCall(PetscFree2(localPoints, remotePoints)); 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 180 PetscFunctionBegin; 181 PetscCall(PetscObjectGetComm((PetscObject)dm,&comm)); 182 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 183 PetscCallMPI(MPI_Comm_size(comm, &size)); 184 PetscCall(DMPatchGetCoarse(dm, &dmc)); 185 PetscCall(DMPatchGetPatchSize(dm, &patchSize)); 186 PetscCall(DMPatchGetCommSize(dm, &commSize)); 187 PetscCall(DMPatchGetCommSize(dm, &commSize)); 188 PetscCall(DMGetGlobalVector(dmc, &XC)); 189 PetscCall(DMDAGetInfo(dmc, NULL, &M, &N, &P, &l, &m, &n, NULL,NULL,NULL,NULL,NULL,NULL)); 190 M = PetscMax(M, 1); l = PetscMax(l, 1); 191 N = PetscMax(N, 1); m = PetscMax(m, 1); 192 P = PetscMax(P, 1); n = PetscMax(n, 1); 193 194 gridRank.i = rank % l; 195 gridRank.j = rank/l % m; 196 gridRank.k = rank/(l*m) % n; 197 198 if (commSize.i*commSize.j*commSize.k == size || commSize.i*commSize.j*commSize.k == 0) { 199 commSize.i = l; commSize.j = m; commSize.k = n; 200 commz = comm; 201 } else if (commSize.i*commSize.j*commSize.k == 1) { 202 commz = PETSC_COMM_SELF; 203 } else { 204 const PetscMPIInt newComm = ((gridRank.k/commSize.k)*(m/commSize.j) + gridRank.j/commSize.j)*(l/commSize.i) + (gridRank.i/commSize.i); 205 const PetscMPIInt newRank = ((gridRank.k%commSize.k)*commSize.j + gridRank.j%commSize.j)*commSize.i + (gridRank.i%commSize.i); 206 207 PetscCallMPI(MPI_Comm_split(comm, newComm, newRank, &commz)); 208 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Rank %d color %d key %d commz %p\n", rank, newComm, newRank, (void*)(MPI_Aint)commz)); 209 } 210 /* 211 Assumptions: 212 - patchSize divides gridSize 213 - commSize divides gridSize 214 - commSize divides l,m,n 215 Ignore multiple patches per rank for now 216 217 Multiple ranks per patch: 218 - l,m,n divides patchSize 219 - commSize divides patchSize 220 */ 221 for (k = 0; k < P; k += PetscMax(patchSize.k, 1)) { 222 for (j = 0; j < N; j += PetscMax(patchSize.j, 1)) { 223 for (i = 0; i < M; i += PetscMax(patchSize.i, 1), ++p) { 224 MPI_Comm commp = MPI_COMM_NULL; 225 DM dmz = NULL; 226 #if 0 227 DM dmf = NULL; 228 Mat interpz = NULL; 229 #endif 230 Vec XZ = NULL; 231 PetscScalar *xcarray = NULL; 232 PetscScalar *xzarray = NULL; 233 234 if ((gridRank.k/commSize.k == p/(l/commSize.i * m/commSize.j) % n/commSize.k) && 235 (gridRank.j/commSize.j == p/(l/commSize.i) % m/commSize.j) && 236 (gridRank.i/commSize.i == p % l/commSize.i)) { 237 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Rank %d is accepting Patch %" PetscInt_FMT "\n", rank, p)); 238 commp = commz; 239 } 240 /* Zoom to coarse patch */ 241 lower.i = i; lower.j = j; lower.k = k; 242 upper.i = i + patchSize.i; upper.j = j + patchSize.j; upper.k = k + patchSize.k; 243 PetscCall(DMPatchZoom(dmc, lower, upper, commp, &dmz, &sfz, &sfzr)); 244 lower.c = 0; /* initialize member, otherwise compiler issues warnings */ 245 upper.c = 0; /* initialize member, otherwise compiler issues warnings */ 246 if (debug) PetscCall(PetscPrintf(comm, "Patch %" PetscInt_FMT ": (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")--(%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", p, lower.i, lower.j, lower.k, upper.i, upper.j, upper.k)); 247 if (dmz) PetscCall(DMView(dmz, PETSC_VIEWER_STDOUT_(commz))); 248 PetscCall(PetscSFView(sfz, PETSC_VIEWER_STDOUT_(comm))); 249 PetscCall(PetscSFView(sfzr, PETSC_VIEWER_STDOUT_(comm))); 250 /* Scatter Xcoarse -> Xzoom */ 251 if (dmz) PetscCall(DMGetGlobalVector(dmz, &XZ)); 252 if (XZ) PetscCall(VecGetArray(XZ, &xzarray)); 253 PetscCall(VecGetArray(XC, &xcarray)); 254 PetscCall(PetscSFBcastBegin(sfz, MPIU_SCALAR, xcarray, xzarray,MPI_REPLACE)); 255 PetscCall(PetscSFBcastEnd(sfz, MPIU_SCALAR, xcarray, xzarray,MPI_REPLACE)); 256 PetscCall(VecRestoreArray(XC, &xcarray)); 257 if (XZ) PetscCall(VecRestoreArray(XZ, &xzarray)); 258 #if 0 259 /* Interpolate Xzoom -> Xfine, note that this may be on subcomms */ 260 PetscCall(DMRefine(dmz, MPI_COMM_NULL, &dmf)); 261 PetscCall(DMCreateInterpolation(dmz, dmf, &interpz, NULL)); 262 PetscCall(DMInterpolate(dmz, interpz, dmf)); 263 /* Smooth Xfine using two-step smoother, normal smoother plus Kaczmarz---moves back and forth from dmzoom to dmfine */ 264 /* Compute residual Rfine */ 265 /* Restrict Rfine to Rzoom_restricted */ 266 #endif 267 /* Scatter Rzoom_restricted -> Rcoarse_restricted */ 268 if (XZ) PetscCall(VecGetArray(XZ, &xzarray)); 269 PetscCall(VecGetArray(XC, &xcarray)); 270 PetscCall(PetscSFReduceBegin(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM)); 271 PetscCall(PetscSFReduceEnd(sfzr, MPIU_SCALAR, xzarray, xcarray, MPIU_SUM)); 272 PetscCall(VecRestoreArray(XC, &xcarray)); 273 if (XZ) PetscCall(VecRestoreArray(XZ, &xzarray)); 274 if (dmz) PetscCall(DMRestoreGlobalVector(dmz, &XZ)); 275 /* Compute global residual Rcoarse */ 276 /* TauCoarse = Rcoarse - Rcoarse_restricted */ 277 278 PetscCall(PetscSFDestroy(&sfz)); 279 PetscCall(PetscSFDestroy(&sfzr)); 280 PetscCall(DMDestroy(&dmz)); 281 } 282 } 283 } 284 PetscCall(DMRestoreGlobalVector(dmc, &XC)); 285 PetscFunctionReturn(0); 286 } 287 288 PetscErrorCode DMPatchView_ASCII(DM dm, PetscViewer viewer) 289 { 290 DM_Patch *mesh = (DM_Patch*) dm->data; 291 PetscViewerFormat format; 292 const char *name; 293 294 PetscFunctionBegin; 295 PetscCall(PetscViewerGetFormat(viewer, &format)); 296 /* if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) */ 297 PetscCall(PetscObjectGetName((PetscObject) dm, &name)); 298 PetscCall(PetscViewerASCIIPrintf(viewer, "Patch DM %s\n", name)); 299 PetscCall(PetscViewerASCIIPushTab(viewer)); 300 PetscCall(PetscViewerASCIIPrintf(viewer, "Coarse DM\n")); 301 PetscCall(DMView(mesh->dmCoarse, viewer)); 302 PetscCall(PetscViewerASCIIPopTab(viewer)); 303 PetscFunctionReturn(0); 304 } 305 306 PetscErrorCode DMView_Patch(DM dm, PetscViewer viewer) 307 { 308 PetscBool iascii, isbinary; 309 310 PetscFunctionBegin; 311 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 312 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 313 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii)); 314 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary)); 315 if (iascii) { 316 PetscCall(DMPatchView_ASCII(dm, viewer)); 317 } 318 PetscFunctionReturn(0); 319 } 320 321 PetscErrorCode DMDestroy_Patch(DM dm) 322 { 323 DM_Patch *mesh = (DM_Patch*) dm->data; 324 325 PetscFunctionBegin; 326 if (--mesh->refct > 0) PetscFunctionReturn(0); 327 PetscCall(DMDestroy(&mesh->dmCoarse)); 328 /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */ 329 PetscCall(PetscFree(mesh)); 330 PetscFunctionReturn(0); 331 } 332 333 PetscErrorCode DMSetUp_Patch(DM dm) 334 { 335 DM_Patch *mesh = (DM_Patch*) dm->data; 336 337 PetscFunctionBegin; 338 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 339 PetscCall(DMSetUp(mesh->dmCoarse)); 340 PetscFunctionReturn(0); 341 } 342 343 PetscErrorCode DMCreateGlobalVector_Patch(DM dm, Vec *g) 344 { 345 DM_Patch *mesh = (DM_Patch*) dm->data; 346 347 PetscFunctionBegin; 348 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 349 PetscCall(DMCreateGlobalVector(mesh->dmCoarse, g)); 350 PetscFunctionReturn(0); 351 } 352 353 PetscErrorCode DMCreateLocalVector_Patch(DM dm, Vec *l) 354 { 355 DM_Patch *mesh = (DM_Patch*) dm->data; 356 357 PetscFunctionBegin; 358 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 359 PetscCall(DMCreateLocalVector(mesh->dmCoarse, l)); 360 PetscFunctionReturn(0); 361 } 362 363 PetscErrorCode DMCreateSubDM_Patch(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm) 364 { 365 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Tell me to code this"); 366 } 367 368 PetscErrorCode DMPatchGetCoarse(DM dm, DM *dmCoarse) 369 { 370 DM_Patch *mesh = (DM_Patch*) dm->data; 371 372 PetscFunctionBegin; 373 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 374 *dmCoarse = mesh->dmCoarse; 375 PetscFunctionReturn(0); 376 } 377 378 PetscErrorCode DMPatchGetPatchSize(DM dm, MatStencil *patchSize) 379 { 380 DM_Patch *mesh = (DM_Patch*) dm->data; 381 382 PetscFunctionBegin; 383 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 384 PetscValidPointer(patchSize, 2); 385 *patchSize = mesh->patchSize; 386 PetscFunctionReturn(0); 387 } 388 389 PetscErrorCode DMPatchSetPatchSize(DM dm, MatStencil patchSize) 390 { 391 DM_Patch *mesh = (DM_Patch*) dm->data; 392 393 PetscFunctionBegin; 394 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 395 mesh->patchSize = patchSize; 396 PetscFunctionReturn(0); 397 } 398 399 PetscErrorCode DMPatchGetCommSize(DM dm, MatStencil *commSize) 400 { 401 DM_Patch *mesh = (DM_Patch*) dm->data; 402 403 PetscFunctionBegin; 404 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 405 PetscValidPointer(commSize, 2); 406 *commSize = mesh->commSize; 407 PetscFunctionReturn(0); 408 } 409 410 PetscErrorCode DMPatchSetCommSize(DM dm, MatStencil commSize) 411 { 412 DM_Patch *mesh = (DM_Patch*) dm->data; 413 414 PetscFunctionBegin; 415 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 416 mesh->commSize = commSize; 417 PetscFunctionReturn(0); 418 } 419