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