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