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