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, PetscInt rank, MPI_Comm commz, DM *dmz, PetscSF *sfz, PetscSF *sfzr) 44 { 45 const PetscInt *lx, *ly, *lz; 46 PetscInt dim, dof, sw; 47 PetscInt M, N, P, rM, rN, rP, halo; 48 DMDAStencilType st; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 if (commz == MPI_COMM_NULL) { 53 /* Split communicator */ 54 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not implemented"); 55 } 56 /* Create patch DM */ 57 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0, &dof, &sw,0,0,0, &st);CHKERRQ(ierr); 58 halo = sw; 59 60 /* Get piece for rank r, expanded by halo */ 61 ierr = DMDAGetOwnershipRanges(dm, &lx, &ly, &lz);CHKERRQ(ierr); 62 rM = PetscMin(M, lx[rank+1] - halo) - PetscMax(lx[rank] - halo, 0); 63 rN = PetscMin(N, ly[rank+1] - halo) - PetscMax(ly[rank] - halo, 0); 64 rP = PetscMin(P, lz[rank+1] - halo) - PetscMax(lz[rank] - halo, 0); 65 66 ierr = DMDACreate(commz, dmz);CHKERRQ(ierr); 67 ierr = DMDASetDim(*dmz, dim);CHKERRQ(ierr); 68 ierr = DMDASetSizes(*dmz, rM, rN, rP);CHKERRQ(ierr); 69 ierr = DMDASetNumProcs(*dmz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);CHKERRQ(ierr); 70 ierr = DMDASetBoundaryType(*dmz, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE);CHKERRQ(ierr); 71 ierr = DMDASetDof(*dmz, dof);CHKERRQ(ierr); 72 ierr = DMDASetStencilType(*dmz, st);CHKERRQ(ierr); 73 ierr = DMDASetStencilWidth(*dmz, sw);CHKERRQ(ierr); 74 ierr = DMDASetOwnershipRanges(*dmz, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 75 ierr = DMSetFromOptions(*dmz);CHKERRQ(ierr); 76 ierr = DMSetUp(*dmz);CHKERRQ(ierr); 77 /* Create SF for ghosted map */ 78 /* Create SF for restricted map */ 79 PetscFunctionReturn(0); 80 } 81 82 #undef __FUNCT__ 83 #define __FUNCT__ "DMPatchSolve" 84 PetscErrorCode DMPatchSolve(DM dm) 85 { 86 MPI_Comm comm = ((PetscObject) dm)->comm; 87 PetscSF sfz, sfr; 88 PetscInt r; 89 PetscMPIInt size; 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 94 for(r = 0; r < size; ++r) { 95 DM dmz; 96 97 ierr = DMPatchZoom(dm, r, comm, &dmz, &sfz, &sfr);CHKERRQ(ierr); 98 #if 0 99 DM dmf; 100 Mat interpz; 101 /* Scatter Xcoarse -> Xzoom */ 102 ierr = PetscSFBcastBegin(sfz, MPIU_SCALAR, xcarray, xzarray);CHKERRQ(ierr); 103 ierr = PetscSFBcastEnd(sfz, MPIU_SCALAR, xcarray, xzarray);CHKERRQ(ierr); 104 /* Interpolate Xzoom -> Xfine, note that this may be on subcomms */ 105 ierr = DMRefine(dmz, MPI_COMM_NULL, &dmf);CHKERRQ(ierr); 106 ierr = DMCreateInterpolation(dmz, dmf, &interpz, PETSC_NULL);CHKERRQ(ierr); 107 ierr = DMInterpolate(dmz, interpz, dmf);CHKERRQ(ierr); 108 /* Smooth Xfine using two-step smoother, normal smoother plus Kaczmarz---moves back and forth from dmzoom to dmfine */ 109 /* Compute residual Rfine */ 110 /* Restrict Rfine to Rzoom_restricted */ 111 /* Scatter Rzoom_restricted -> Rcoarse_restricted */ 112 ierr = PetscSFBcastBegin(sfzr, MPIU_SCALAR, xzarray, xcarray, MPI_SUM);CHKERRQ(ierr); 113 ierr = PetscSFBcastEnd(sfzr, MPIU_SCALAR, xzarray, xcarray, MPI_SUM);CHKERRQ(ierr); 114 /* Compute global residual Rcoarse */ 115 /* TauCoarse = Rcoarse - Rcoarse_restricted */ 116 #endif 117 } 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNCT__ 122 #define __FUNCT__ "DMPatchView_Ascii" 123 PetscErrorCode DMPatchView_Ascii(DM dm, PetscViewer viewer) 124 { 125 DM_Patch *mesh = (DM_Patch *) dm->data; 126 PetscViewerFormat format; 127 PetscInt p; 128 PetscErrorCode ierr; 129 130 PetscFunctionBegin; 131 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 132 /* if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) */ 133 ierr = PetscViewerASCIIPrintf(viewer, "%D patches\n", mesh->numPatches);CHKERRQ(ierr); 134 for(p = 0; p < mesh->numPatches; ++p) { 135 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 136 ierr = PetscViewerASCIIPrintf(viewer, "Patch %D\n", p);CHKERRQ(ierr); 137 ierr = DMView(mesh->patches[p], viewer);CHKERRQ(ierr); 138 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 139 } 140 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 141 ierr = PetscViewerASCIIPrintf(viewer, "Coarse DM\n");CHKERRQ(ierr); 142 ierr = DMView(mesh->dmCoarse, viewer);CHKERRQ(ierr); 143 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 147 #undef __FUNCT__ 148 #define __FUNCT__ "DMView_Patch" 149 PetscErrorCode DMView_Patch(DM dm, PetscViewer viewer) 150 { 151 PetscBool iascii, isbinary; 152 PetscErrorCode ierr; 153 154 PetscFunctionBegin; 155 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 156 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 157 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 158 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);CHKERRQ(ierr); 159 if (iascii) { 160 ierr = DMPatchView_Ascii(dm, viewer);CHKERRQ(ierr); 161 #if 0 162 } else if (isbinary) { 163 ierr = DMPatchView_Binary(dm, viewer);CHKERRQ(ierr); 164 #endif 165 } else SETERRQ1(((PetscObject)viewer)->comm,PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name); 166 PetscFunctionReturn(0); 167 } 168 169 #undef __FUNCT__ 170 #define __FUNCT__ "DMDestroy_Patch" 171 PetscErrorCode DMDestroy_Patch(DM dm) 172 { 173 DM_Patch *mesh = (DM_Patch *) dm->data; 174 PetscInt p; 175 PetscErrorCode ierr; 176 177 PetscFunctionBegin; 178 if (--mesh->refct > 0) {PetscFunctionReturn(0);} 179 for(p = 0; p < mesh->numPatches; ++p) { 180 ierr = DMDestroy(&mesh->patches[p]);CHKERRQ(ierr); 181 } 182 ierr = PetscFree(mesh->patches);CHKERRQ(ierr); 183 ierr = DMDestroy(&mesh->dmCoarse);CHKERRQ(ierr); 184 /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */ 185 ierr = PetscFree(mesh);CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 #undef __FUNCT__ 190 #define __FUNCT__ "DMSetUp_Patch" 191 PetscErrorCode DMSetUp_Patch(DM dm) 192 { 193 DM_Patch *mesh = (DM_Patch *) dm->data; 194 PetscInt p; 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 199 for(p = 0; p < mesh->numPatches; ++p) { 200 ierr = DMSetUp(mesh->patches[p]);CHKERRQ(ierr); 201 } 202 PetscFunctionReturn(0); 203 } 204 205 #undef __FUNCT__ 206 #define __FUNCT__ "DMCreateGlobalVector_Patch" 207 PetscErrorCode DMCreateGlobalVector_Patch(DM dm, Vec *g) 208 { 209 DM_Patch *mesh = (DM_Patch *) dm->data; 210 PetscErrorCode ierr; 211 212 PetscFunctionBegin; 213 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 214 ierr = DMCreateGlobalVector(mesh->patches[mesh->activePatch], g);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "DMCreateLocalVector_Patch" 220 PetscErrorCode DMCreateLocalVector_Patch(DM dm, Vec *l) 221 { 222 DM_Patch *mesh = (DM_Patch *) dm->data; 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 227 ierr = DMCreateLocalVector(mesh->patches[mesh->activePatch], l);CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "DMPatchGetActivePatch" 233 PetscErrorCode DMPatchGetActivePatch(DM dm, PetscInt *patch) 234 { 235 DM_Patch *mesh; 236 237 PetscFunctionBegin; 238 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 239 PetscValidPointer(patch, 2); 240 mesh = (DM_Patch *) dm->data; 241 *patch = mesh->activePatch; 242 PetscFunctionReturn(0); 243 } 244 245 #undef __FUNCT__ 246 #define __FUNCT__ "DMPatchSetActivePatch" 247 PetscErrorCode DMPatchSetActivePatch(DM dm, PetscInt patch) 248 { 249 DM_Patch *mesh; 250 251 PetscFunctionBegin; 252 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 253 mesh = (DM_Patch *) dm->data; 254 mesh->activePatch = patch; 255 PetscFunctionReturn(0); 256 } 257 258 #undef __FUNCT__ 259 #define __FUNCT__ "DMCreateSubDM_Patch" 260 PetscErrorCode DMCreateSubDM_Patch(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm) 261 { 262 PetscFunctionBegin; 263 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 264 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Tell me to code this"); 265 PetscFunctionReturn(0); 266 } 267