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