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