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