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