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; 47 PetscInt *localPoints; 48 PetscSFNode *remotePoints; 49 PetscInt dim, dof; 50 PetscInt M, N, P, rM, rN, rP, halo = 1, sx, sy, sz, sxb, syb, szb, sxr, syr, szr, mxb, myb, mzb, mxr, myr, mzr, i, j, k, q; 51 PetscErrorCode ierr; 52 53 PetscFunctionBegin; 54 if (commz == MPI_COMM_NULL) { 55 /* Split communicator */ 56 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not implemented"); 57 } 58 /* Create patch DM */ 59 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0, &dof, 0,0,0,0, &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 ierr = DMDACreate(commz, dmz);CHKERRQ(ierr); 70 ierr = DMDASetDim(*dmz, dim);CHKERRQ(ierr); 71 ierr = DMDASetSizes(*dmz, rM, rN, rP);CHKERRQ(ierr); 72 ierr = DMDASetNumProcs(*dmz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);CHKERRQ(ierr); 73 ierr = DMDASetBoundaryType(*dmz, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE);CHKERRQ(ierr); 74 ierr = DMDASetDof(*dmz, dof);CHKERRQ(ierr); 75 ierr = DMDASetStencilType(*dmz, st);CHKERRQ(ierr); 76 ierr = DMDASetStencilWidth(*dmz, 0);CHKERRQ(ierr); 77 ierr = DMDASetOwnershipRanges(*dmz, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 78 ierr = DMSetFromOptions(*dmz);CHKERRQ(ierr); 79 ierr = DMSetUp(*dmz);CHKERRQ(ierr); 80 ierr = DMDAGetCorners(*dmz, &sx, &sy, &sz, &mxb, &myb, &mzb);CHKERRQ(ierr); 81 82 ierr = PetscMalloc2(rM*rN*rP,PetscInt,&localPoints,rM*rN*rP,PetscSFNode,&remotePoints);CHKERRQ(ierr); 83 84 /* Create SF for restricted map */ 85 q = 0; 86 szr = lower.k + sz; syr = lower.j + sy; sxr = lower.i + sx; 87 /* This is not quite right. Only works for serial patches */ 88 mzr = mzb - (bupper.k-upper.k) - (lower.k-blower.k); 89 myr = myb - (bupper.j-upper.j) - (lower.j-blower.j); 90 mxr = mxb - (bupper.i-upper.i) - (lower.i-blower.i); 91 for(k = szr; k < szr+mzr; ++k) { 92 for(j = syr; j < syr+myr; ++j) { 93 for(i = sxr; i < sxr+mxr; ++i, ++q) { 94 const PetscInt lp = ((k-szr)*rN + (j-syr))*rM + i-sxr; 95 const PetscInt gp = (k*N + j)*M + i; 96 97 localPoints[q] = lp; 98 /* Replace with a function that can figure this out */ 99 remotePoints[q].rank = 0; 100 remotePoints[q].index = gp; 101 } 102 } 103 } 104 ierr = PetscSFCreate(((PetscObject) dm)->comm, sfzr);CHKERRQ(ierr); 105 ierr = PetscSFSetGraph(*sfzr, M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);CHKERRQ(ierr); 106 107 /* Create SF for buffered map */ 108 q = 0; 109 szb = blower.k + sz; syb = blower.j + sy; sxb = blower.i + sx; 110 for(k = szb; k < szb+mzb; ++k) { 111 for(j = syb; j < syb+myb; ++j) { 112 for(i = sxb; i < sxb+mxb; ++i, ++q) { 113 const PetscInt lp = ((k-szb)*rN + (j-syb))*rM + i-sxb; 114 const PetscInt gp = (k*N + j)*M + i; 115 116 localPoints[q] = lp; 117 /* Replace with a function that can figure this out */ 118 remotePoints[q].rank = 0; 119 remotePoints[q].index = gp; 120 } 121 } 122 } 123 ierr = PetscSFCreate(((PetscObject) dm)->comm, sfz);CHKERRQ(ierr); 124 ierr = PetscSFSetGraph(*sfz, M*N*P, q, localPoints, PETSC_COPY_VALUES, remotePoints, PETSC_COPY_VALUES);CHKERRQ(ierr); 125 126 ierr = PetscFree2(localPoints, remotePoints);CHKERRQ(ierr); 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "DMPatchSolve" 132 PetscErrorCode DMPatchSolve(DM dm) 133 { 134 MPI_Comm comm = ((PetscObject) dm)->comm; 135 DM_Patch *mesh = (DM_Patch *) dm->data; 136 DM cdm = mesh->dmCoarse; 137 DM dmz; 138 PetscSF sfz, sfzr; 139 MatStencil patchSize, lower, upper; 140 PetscInt M, N, P, i, j, k, p = 0; 141 PetscErrorCode ierr; 142 143 PetscFunctionBegin; 144 ierr = DMPatchGetPatchSize(dm, &patchSize);CHKERRQ(ierr); 145 ierr = DMDAGetInfo(cdm, 0, &M, &N, &P, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 146 M = PetscMax(M, 1); 147 N = PetscMax(N, 1); 148 P = PetscMax(P, 1); 149 for(k = 0; k < P; k += PetscMax(patchSize.k, 1)) { 150 for(j = 0; j < N; j += PetscMax(patchSize.j, 1)) { 151 for(i = 0; i < M; i += PetscMax(patchSize.i, 1), ++p) { 152 DM dmf; 153 Mat interpz; 154 PetscScalar *xcarray, *xzarray; 155 156 /* Zoom to coarse patch */ 157 lower.i = i; lower.j = j; lower.k = k; 158 upper.i = i + patchSize.i; upper.j = j + patchSize.j; upper.k = k + patchSize.k; 159 ierr = DMPatchZoom(cdm, lower, upper, comm, &dmz, &sfz, &sfzr);CHKERRQ(ierr); 160 /* Debug */ 161 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); 162 ierr = DMView(dmz, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 163 ierr = PetscSFView(sfz, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 164 ierr = PetscSFView(sfzr, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 165 /* TODO: Need coarse and zoomed state vectors here: */ 166 xcarray = PETSC_NULL; 167 xzarray = PETSC_NULL; 168 /* Scatter Xcoarse -> Xzoom */ 169 ierr = PetscSFBcastBegin(sfz, MPIU_SCALAR, xcarray, xzarray);CHKERRQ(ierr); 170 ierr = PetscSFBcastEnd(sfz, MPIU_SCALAR, xcarray, xzarray);CHKERRQ(ierr); 171 /* Interpolate Xzoom -> Xfine, note that this may be on subcomms */ 172 ierr = DMRefine(dmz, MPI_COMM_NULL, &dmf);CHKERRQ(ierr); 173 ierr = DMCreateInterpolation(dmz, dmf, &interpz, PETSC_NULL);CHKERRQ(ierr); 174 ierr = DMInterpolate(dmz, interpz, dmf);CHKERRQ(ierr); 175 /* Smooth Xfine using two-step smoother, normal smoother plus Kaczmarz---moves back and forth from dmzoom to dmfine */ 176 /* Compute residual Rfine */ 177 /* Restrict Rfine to Rzoom_restricted */ 178 /* Scatter Rzoom_restricted -> Rcoarse_restricted */ 179 ierr = PetscSFReduceBegin(sfzr, MPIU_SCALAR, xzarray, xcarray, MPI_SUM);CHKERRQ(ierr); 180 ierr = PetscSFReduceEnd(sfzr, MPIU_SCALAR, xzarray, xcarray, MPI_SUM);CHKERRQ(ierr); 181 /* Compute global residual Rcoarse */ 182 /* TauCoarse = Rcoarse - Rcoarse_restricted */ 183 184 ierr = PetscSFDestroy(&sfz);CHKERRQ(ierr); 185 ierr = PetscSFDestroy(&sfzr);CHKERRQ(ierr); 186 ierr = DMDestroy(&dmz);CHKERRQ(ierr); 187 } 188 } 189 } 190 PetscFunctionReturn(0); 191 } 192 193 #undef __FUNCT__ 194 #define __FUNCT__ "DMPatchView_Ascii" 195 PetscErrorCode DMPatchView_Ascii(DM dm, PetscViewer viewer) 196 { 197 DM_Patch *mesh = (DM_Patch *) dm->data; 198 PetscViewerFormat format; 199 const char *name; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 204 /* if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) */ 205 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 206 ierr = PetscViewerASCIIPrintf(viewer, "Patch DM %s\n", name);CHKERRQ(ierr); 207 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 208 ierr = PetscViewerASCIIPrintf(viewer, "Coarse DM\n");CHKERRQ(ierr); 209 ierr = DMView(mesh->dmCoarse, viewer);CHKERRQ(ierr); 210 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 211 PetscFunctionReturn(0); 212 } 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "DMView_Patch" 216 PetscErrorCode DMView_Patch(DM dm, PetscViewer viewer) 217 { 218 PetscBool iascii, isbinary; 219 PetscErrorCode ierr; 220 221 PetscFunctionBegin; 222 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 223 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 224 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 225 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);CHKERRQ(ierr); 226 if (iascii) { 227 ierr = DMPatchView_Ascii(dm, viewer);CHKERRQ(ierr); 228 #if 0 229 } else if (isbinary) { 230 ierr = DMPatchView_Binary(dm, viewer);CHKERRQ(ierr); 231 #endif 232 } else SETERRQ1(((PetscObject)viewer)->comm,PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name); 233 PetscFunctionReturn(0); 234 } 235 236 #undef __FUNCT__ 237 #define __FUNCT__ "DMDestroy_Patch" 238 PetscErrorCode DMDestroy_Patch(DM dm) 239 { 240 DM_Patch *mesh = (DM_Patch *) dm->data; 241 PetscErrorCode ierr; 242 243 PetscFunctionBegin; 244 if (--mesh->refct > 0) {PetscFunctionReturn(0);} 245 ierr = DMDestroy(&mesh->dmCoarse);CHKERRQ(ierr); 246 /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */ 247 ierr = PetscFree(mesh);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 251 #undef __FUNCT__ 252 #define __FUNCT__ "DMSetUp_Patch" 253 PetscErrorCode DMSetUp_Patch(DM dm) 254 { 255 DM_Patch *mesh = (DM_Patch *) dm->data; 256 PetscErrorCode ierr; 257 258 PetscFunctionBegin; 259 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 260 ierr = DMSetUp(mesh->dmCoarse);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "DMCreateGlobalVector_Patch" 266 PetscErrorCode DMCreateGlobalVector_Patch(DM dm, Vec *g) 267 { 268 DM_Patch *mesh = (DM_Patch *) dm->data; 269 PetscErrorCode ierr; 270 271 PetscFunctionBegin; 272 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 273 ierr = DMCreateGlobalVector(mesh->dmCoarse, g);CHKERRQ(ierr); 274 PetscFunctionReturn(0); 275 } 276 277 #undef __FUNCT__ 278 #define __FUNCT__ "DMCreateLocalVector_Patch" 279 PetscErrorCode DMCreateLocalVector_Patch(DM dm, Vec *l) 280 { 281 DM_Patch *mesh = (DM_Patch *) dm->data; 282 PetscErrorCode ierr; 283 284 PetscFunctionBegin; 285 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 286 ierr = DMCreateLocalVector(mesh->dmCoarse, l);CHKERRQ(ierr); 287 PetscFunctionReturn(0); 288 } 289 290 #undef __FUNCT__ 291 #define __FUNCT__ "DMCreateSubDM_Patch" 292 PetscErrorCode DMCreateSubDM_Patch(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm) 293 { 294 PetscFunctionBegin; 295 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 296 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Tell me to code this"); 297 PetscFunctionReturn(0); 298 } 299 300 #undef __FUNCT__ 301 #define __FUNCT__ "DMPatchGetCoarse" 302 PetscErrorCode DMPatchGetCoarse(DM dm, DM *dmCoarse) 303 { 304 DM_Patch *mesh = (DM_Patch *) dm->data; 305 306 PetscFunctionBegin; 307 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 308 *dmCoarse = mesh->dmCoarse; 309 PetscFunctionReturn(0); 310 } 311 312 #undef __FUNCT__ 313 #define __FUNCT__ "DMPatchGetPatchSize" 314 PetscErrorCode DMPatchGetPatchSize(DM dm, MatStencil *patchSize) 315 { 316 DM_Patch *mesh = (DM_Patch *) dm->data; 317 318 PetscFunctionBegin; 319 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 320 PetscValidPointer(patchSize, 2); 321 *patchSize = mesh->patchSize; 322 PetscFunctionReturn(0); 323 } 324 325 #undef __FUNCT__ 326 #define __FUNCT__ "DMPatchSetPatchSize" 327 PetscErrorCode DMPatchSetPatchSize(DM dm, MatStencil patchSize) 328 { 329 DM_Patch *mesh = (DM_Patch *) dm->data; 330 331 PetscFunctionBegin; 332 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 333 mesh->patchSize = patchSize; 334 PetscFunctionReturn(0); 335 } 336