1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/dmlabelimpl.h> 3 #include <petsc/private/sectionimpl.h> 4 #include <petsc/private/sfimpl.h> 5 6 static PetscErrorCode DMTransferMaterialParameters(DM dm, PetscSF sf, DM odm) 7 { 8 Vec A; 9 10 PetscFunctionBegin; 11 /* TODO handle regions? */ 12 PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &A)); 13 if (A) { 14 DM dmAux, ocdm, odmAux; 15 Vec oA, oAt; 16 PetscSection sec, osec; 17 18 PetscCall(VecGetDM(A, &dmAux)); 19 PetscCall(DMClone(odm, &odmAux)); 20 PetscCall(DMGetCoordinateDM(odm, &ocdm)); 21 PetscCall(DMSetCoordinateDM(odmAux, ocdm)); 22 PetscCall(DMCopyDisc(dmAux, odmAux)); 23 24 PetscCall(DMGetLocalSection(dmAux, &sec)); 25 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)sec), &osec)); 26 PetscCall(VecCreate(PetscObjectComm((PetscObject)A), &oAt)); 27 PetscCall(DMPlexDistributeField(dmAux, sf, sec, A, osec, oAt)); 28 PetscCall(DMSetLocalSection(odmAux, osec)); 29 PetscCall(PetscSectionDestroy(&osec)); 30 PetscCall(DMCreateLocalVector(odmAux, &oA)); 31 PetscCall(DMDestroy(&odmAux)); 32 PetscCall(VecCopy(oAt, oA)); 33 PetscCall(VecDestroy(&oAt)); 34 PetscCall(DMSetAuxiliaryVec(odm, NULL, 0, 0, oA)); 35 PetscCall(VecDestroy(&oA)); 36 } 37 PetscFunctionReturn(PETSC_SUCCESS); 38 } 39 40 PetscErrorCode DMCreateDomainDecomposition_Plex(DM dm, PetscInt *nsub, char ***names, IS **innerises, IS **outerises, DM **dms) 41 { 42 DM odm; 43 PetscSF migrationSF = NULL, sectionSF; 44 PetscSection sec, tsec, ogsec, olsec; 45 PetscInt n, mh, ddovl = 0, pStart, pEnd, ni, no, nl; 46 PetscDS ds; 47 DMLabel label; 48 const char *oname = "__internal_plex_dd_ovl_"; 49 IS gi_is, li_is, go_is, gl_is, ll_is; 50 IS gis, lis; 51 PetscInt rst, ren, c, *gidxs, *lidxs, *tidxs; 52 Vec gvec; 53 54 PetscFunctionBegin; 55 n = 1; 56 if (nsub) *nsub = n; 57 if (names) PetscCall(PetscCalloc1(n, names)); 58 if (innerises) PetscCall(PetscCalloc1(n, innerises)); 59 if (outerises) PetscCall(PetscCalloc1(n, outerises)); 60 if (dms) PetscCall(PetscCalloc1(n, dms)); 61 62 PetscObjectOptionsBegin((PetscObject)dm); 63 PetscCall(PetscOptionsBoundedInt("-dm_plex_dd_overlap", "The size of the overlap halo for the subdomains", "DMCreateDomainDecomposition", ddovl, &ddovl, NULL, 0)); 64 PetscOptionsEnd(); 65 66 PetscCall(DMViewFromOptions(dm, NULL, "-dm_plex_dd_dm_view")); 67 PetscCall(DMPlexDistributeOverlap_Internal(dm, ddovl + 1, PETSC_COMM_SELF, oname, &migrationSF, &odm)); 68 if (!odm) PetscCall(DMClone(dm, &odm)); 69 if (migrationSF) PetscCall(PetscSFViewFromOptions(migrationSF, (PetscObject)dm, "-dm_plex_dd_sf_view")); 70 71 PetscCall(DMPlexGetMaxProjectionHeight(dm, &mh)); 72 PetscCall(DMPlexSetMaxProjectionHeight(odm, mh)); 73 74 /* share discretization */ 75 /* TODO Labels for regions may need to updated, 76 now it uses the original ones, not the ones for the odm. 77 Not sure what to do here */ 78 /* PetscCall(DMCopyDisc(dm, odm)); */ 79 PetscCall(DMGetDS(odm, &ds)); 80 if (!ds) { 81 PetscCall(DMGetLocalSection(dm, &sec)); 82 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)odm), &tsec)); 83 if (migrationSF) { 84 PetscCall(PetscSFDistributeSection(migrationSF, sec, NULL, tsec)); 85 } else { 86 PetscCall(PetscSectionCopy(sec, tsec)); 87 } 88 PetscCall(DMSetLocalSection(dm, tsec)); 89 PetscCall(PetscSectionDestroy(&tsec)); 90 } 91 /* TODO: what if these values changes? add to some DM hook? */ 92 PetscCall(DMTransferMaterialParameters(dm, migrationSF, odm)); 93 94 /* TODO: it would be nice to automatically translate filenames with wildcards: 95 name%r.vtu -> name${rank}.vtu 96 name%R.vtu -> name${PetscGlobalRank}.vtu 97 So that -dm_view vtk:name%R.vtu will automatically translate to the commented code below 98 */ 99 #if 0 100 { 101 char name[256]; 102 PetscViewer viewer; 103 104 PetscCall(PetscSNPrintf(name, sizeof(name), "ovl_mesh_%d.vtu", PetscGlobalRank)); 105 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)odm), name, FILE_MODE_WRITE, &viewer)); 106 PetscCall(DMView(odm, viewer)); 107 PetscCall(PetscViewerDestroy(&viewer)); 108 } 109 #endif 110 PetscCall(DMViewFromOptions(odm, (PetscObject)dm, "-dm_plex_dd_overlap_dm_view")); 111 112 /* propagate original global ordering to overlapping DM */ 113 PetscCall(DMGetSectionSF(dm, §ionSF)); 114 PetscCall(DMGetLocalSection(dm, &sec)); 115 PetscCall(PetscSectionGetStorageSize(sec, &nl)); 116 PetscCall(DMGetGlobalVector(dm, &gvec)); 117 PetscCall(VecGetOwnershipRange(gvec, &rst, &ren)); 118 PetscCall(DMRestoreGlobalVector(dm, &gvec)); 119 PetscCall(ISCreateStride(PETSC_COMM_SELF, ren - rst, rst, 1, &gi_is)); /* non-overlapping dofs */ 120 PetscCall(PetscMalloc1(nl, &lidxs)); 121 for (PetscInt i = 0; i < nl; i++) lidxs[i] = -1; 122 PetscCall(ISGetIndices(gi_is, (const PetscInt **)&gidxs)); 123 PetscCall(PetscSFBcastBegin(sectionSF, MPIU_INT, gidxs, lidxs, MPI_REPLACE)); 124 PetscCall(PetscSFBcastEnd(sectionSF, MPIU_INT, gidxs, lidxs, MPI_REPLACE)); 125 PetscCall(ISRestoreIndices(gi_is, (const PetscInt **)&gidxs)); 126 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), nl, lidxs, PETSC_OWN_POINTER, &lis)); 127 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)odm), &tsec)); 128 PetscCall(DMPlexDistributeFieldIS(dm, migrationSF, sec, lis, tsec, &gis)); 129 PetscCall(PetscSectionDestroy(&tsec)); 130 PetscCall(ISDestroy(&lis)); 131 PetscCall(PetscSFDestroy(&migrationSF)); 132 133 /* make dofs on the overlap boundary (not the global boundary) constrained */ 134 PetscCall(DMGetLabel(odm, oname, &label)); 135 PetscCall(DMPlexLabelComplete(odm, label)); 136 PetscCall(DMGetLocalSection(odm, &tsec)); 137 PetscCall(PetscSectionGetChart(tsec, &pStart, &pEnd)); 138 PetscCall(DMLabelCreateIndex(label, pStart, pEnd)); 139 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)tsec), &sec)); 140 PetscCall(PetscSectionCopy_Internal(tsec, sec, label->bt)); 141 PetscCall(DMSetLocalSection(odm, sec)); 142 PetscCall(PetscSectionDestroy(&sec)); 143 PetscCall(DMRemoveLabel(odm, oname, NULL)); 144 145 /* Create index sets for dofs in the overlap dm */ 146 PetscCall(DMGetSectionSF(odm, §ionSF)); 147 PetscCall(DMGetLocalSection(odm, &olsec)); 148 PetscCall(DMGetGlobalSection(odm, &ogsec)); 149 PetscCall(PetscSectionViewFromOptions(ogsec, (PetscObject)dm, "-dm_plex_dd_overlap_gsection_view")); 150 PetscCall(PetscSectionViewFromOptions(olsec, (PetscObject)dm, "-dm_plex_dd_overlap_lsection_view")); 151 ni = ren - rst; 152 PetscCall(PetscSectionGetConstrainedStorageSize(ogsec, &no)); /* dofs in the overlap */ 153 PetscCall(PetscSectionGetStorageSize(olsec, &nl)); /* local dofs in the overlap */ 154 PetscCall(PetscMalloc1(no, &gidxs)); 155 PetscCall(ISGetIndices(gis, (const PetscInt **)&lidxs)); 156 PetscCall(PetscSFReduceBegin(sectionSF, MPIU_INT, lidxs, gidxs, MPI_REPLACE)); 157 PetscCall(PetscSFReduceEnd(sectionSF, MPIU_INT, lidxs, gidxs, MPI_REPLACE)); 158 PetscCall(ISRestoreIndices(gis, (const PetscInt **)&lidxs)); 159 160 /* non-overlapping dofs */ 161 PetscCall(PetscMalloc1(no, &lidxs)); 162 c = 0; 163 for (PetscInt i = 0; i < no; i++) { 164 if (gidxs[i] >= rst && gidxs[i] < ren) lidxs[c++] = i; 165 } 166 PetscCheck(c == ni, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%" PetscInt_FMT " != %" PetscInt_FMT "\n", c, ni); 167 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)odm), ni, lidxs, PETSC_OWN_POINTER, &li_is)); 168 169 /* global dofs in the overlap */ 170 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), no, gidxs, PETSC_OWN_POINTER, &go_is)); 171 PetscCall(ISViewFromOptions(go_is, (PetscObject)dm, "-dm_plex_dd_overlap_gois_view")); 172 /* PetscCall(ISCreateStride(PetscObjectComm((PetscObject)odm), no, 0, 1, &lo_is)); */ 173 174 /* local dofs of the overlapping subdomain (we actually need only dofs on the boundary of the subdomain) */ 175 PetscCall(PetscMalloc1(nl, &lidxs)); 176 PetscCall(PetscMalloc1(nl, &gidxs)); 177 PetscCall(ISGetIndices(gis, (const PetscInt **)&tidxs)); 178 c = 0; 179 for (PetscInt i = 0; i < nl; i++) { 180 if (tidxs[i] < 0) continue; 181 lidxs[c] = i; 182 gidxs[c] = tidxs[i]; 183 c++; 184 } 185 PetscCall(ISRestoreIndices(gis, (const PetscInt **)&tidxs)); 186 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), c, gidxs, PETSC_OWN_POINTER, &gl_is)); 187 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)odm), c, lidxs, PETSC_OWN_POINTER, &ll_is)); 188 PetscCall(ISViewFromOptions(gl_is, (PetscObject)dm, "-dm_plex_dd_overlap_glis_view")); 189 190 PetscCall(PetscObjectCompose((PetscObject)odm, "__Plex_DD_IS_gi", (PetscObject)gi_is)); 191 PetscCall(PetscObjectCompose((PetscObject)odm, "__Plex_DD_IS_li", (PetscObject)li_is)); 192 PetscCall(PetscObjectCompose((PetscObject)odm, "__Plex_DD_IS_go", (PetscObject)go_is)); 193 PetscCall(PetscObjectCompose((PetscObject)odm, "__Plex_DD_IS_gl", (PetscObject)gl_is)); 194 PetscCall(PetscObjectCompose((PetscObject)odm, "__Plex_DD_IS_ll", (PetscObject)ll_is)); 195 196 if (innerises) (*innerises)[0] = gi_is; 197 else PetscCall(ISDestroy(&gi_is)); 198 if (outerises) (*outerises)[0] = go_is; 199 else PetscCall(ISDestroy(&go_is)); 200 if (dms) (*dms)[0] = odm; 201 else PetscCall(DMDestroy(&odm)); 202 PetscCall(ISDestroy(&li_is)); 203 PetscCall(ISDestroy(&gl_is)); 204 PetscCall(ISDestroy(&ll_is)); 205 PetscCall(ISDestroy(&gis)); 206 207 PetscFunctionReturn(PETSC_SUCCESS); 208 } 209 210 PetscErrorCode DMCreateDomainDecompositionScatters_Plex(DM dm, PetscInt n, DM *subdms, VecScatter **iscat, VecScatter **oscat, VecScatter **lscat) 211 { 212 Vec gvec, svec, lvec; 213 IS gi_is, li_is, go_is, gl_is, ll_is; 214 215 PetscFunctionBegin; 216 if (iscat) PetscCall(PetscMalloc1(n, iscat)); 217 if (oscat) PetscCall(PetscMalloc1(n, oscat)); 218 if (lscat) PetscCall(PetscMalloc1(n, lscat)); 219 220 PetscCall(DMGetGlobalVector(dm, &gvec)); 221 for (PetscInt i = 0; i < n; i++) { 222 PetscCall(DMGetGlobalVector(subdms[i], &svec)); 223 PetscCall(DMGetLocalVector(subdms[i], &lvec)); 224 PetscCall(PetscObjectQuery((PetscObject)subdms[i], "__Plex_DD_IS_gi", (PetscObject *)&gi_is)); 225 PetscCall(PetscObjectQuery((PetscObject)subdms[i], "__Plex_DD_IS_li", (PetscObject *)&li_is)); 226 PetscCall(PetscObjectQuery((PetscObject)subdms[i], "__Plex_DD_IS_go", (PetscObject *)&go_is)); 227 PetscCall(PetscObjectQuery((PetscObject)subdms[i], "__Plex_DD_IS_gl", (PetscObject *)&gl_is)); 228 PetscCall(PetscObjectQuery((PetscObject)subdms[i], "__Plex_DD_IS_ll", (PetscObject *)&ll_is)); 229 PetscCheck(gi_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "SubDM not obtained from DMCreateDomainDecomposition"); 230 PetscCheck(li_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "SubDM not obtained from DMCreateDomainDecomposition"); 231 PetscCheck(go_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "SubDM not obtained from DMCreateDomainDecomposition"); 232 PetscCheck(gl_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "SubDM not obtained from DMCreateDomainDecomposition"); 233 PetscCheck(ll_is, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "SubDM not obtained from DMCreateDomainDecomposition"); 234 if (iscat) PetscCall(VecScatterCreate(gvec, gi_is, svec, li_is, &(*iscat)[i])); 235 if (oscat) PetscCall(VecScatterCreate(gvec, go_is, svec, NULL, &(*oscat)[i])); 236 if (lscat) PetscCall(VecScatterCreate(gvec, gl_is, lvec, ll_is, &(*lscat)[i])); 237 PetscCall(DMRestoreGlobalVector(subdms[i], &svec)); 238 PetscCall(DMRestoreLocalVector(subdms[i], &lvec)); 239 } 240 PetscCall(DMRestoreGlobalVector(dm, &gvec)); 241 PetscFunctionReturn(PETSC_SUCCESS); 242 } 243 244 /* 245 DMCreateNeumannOverlap_Plex - Generates an IS, an unassembled (Neumann) Mat, a setup function, and the corresponding context to be used by PCHPDDM. 246 247 Input Parameter: 248 . dm - preconditioner context 249 250 Output Parameters: 251 + ovl - index set of overlapping subdomains 252 . J - unassembled (Neumann) local matrix 253 . setup - function for generating the matrix 254 - setup_ctx - context for setup 255 256 Options Database Keys: 257 + -dm_plex_view_neumann_original - view the DM without overlap 258 - -dm_plex_view_neumann_overlap - view the DM with overlap as needed by PCHPDDM 259 260 Level: advanced 261 262 .seealso: `DMCreate()`, `DM`, `MATIS`, `PCHPDDM`, `PCHPDDMSetAuxiliaryMat()` 263 */ 264 PetscErrorCode DMCreateNeumannOverlap_Plex(DM dm, IS *ovl, Mat *J, PetscErrorCode (**setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **setup_ctx) 265 { 266 DM odm; 267 Mat pJ; 268 PetscSF sf = NULL; 269 PetscSection sec, osec; 270 ISLocalToGlobalMapping l2g; 271 const PetscInt *idxs; 272 PetscInt n, mh; 273 274 PetscFunctionBegin; 275 *setup = NULL; 276 *setup_ctx = NULL; 277 *ovl = NULL; 278 *J = NULL; 279 280 /* Overlapped mesh 281 overlap is a little more generous, since it is not computed starting from the owned (Dirichlet) points, but from the locally owned cells */ 282 PetscCall(DMPlexDistributeOverlap(dm, 1, &sf, &odm)); 283 if (!odm) { 284 PetscCall(PetscSFDestroy(&sf)); 285 PetscFunctionReturn(PETSC_SUCCESS); 286 } 287 288 /* share discretization */ 289 PetscCall(DMGetLocalSection(dm, &sec)); 290 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)sec), &osec)); 291 PetscCall(PetscSFDistributeSection(sf, sec, NULL, osec)); 292 /* what to do here? using both is fine? */ 293 PetscCall(DMSetLocalSection(odm, osec)); 294 PetscCall(DMCopyDisc(dm, odm)); 295 PetscCall(DMPlexGetMaxProjectionHeight(dm, &mh)); 296 PetscCall(DMPlexSetMaxProjectionHeight(odm, mh)); 297 PetscCall(PetscSectionDestroy(&osec)); 298 299 /* material parameters */ 300 /* TODO: what if these values changes? add to some DM hook? */ 301 PetscCall(DMTransferMaterialParameters(dm, sf, odm)); 302 PetscCall(PetscSFDestroy(&sf)); 303 304 PetscCall(DMViewFromOptions(dm, NULL, "-dm_plex_view_neumann_original")); 305 PetscCall(PetscObjectSetName((PetscObject)odm, "OVL")); 306 PetscCall(DMViewFromOptions(odm, NULL, "-dm_plex_view_neumann_overlap")); 307 308 /* MATIS for the overlap region 309 the HPDDM interface wants local matrices, 310 we attach the global MATIS to the overlap IS, 311 since we need it to do assembly */ 312 PetscCall(DMSetMatType(odm, MATIS)); 313 PetscCall(DMCreateMatrix(odm, &pJ)); 314 PetscCall(MatISGetLocalMat(pJ, J)); 315 PetscCall(PetscObjectReference((PetscObject)*J)); 316 317 /* overlap IS */ 318 PetscCall(MatISGetLocalToGlobalMapping(pJ, &l2g, NULL)); 319 PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n)); 320 PetscCall(ISLocalToGlobalMappingGetIndices(l2g, &idxs)); 321 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)odm), n, idxs, PETSC_COPY_VALUES, ovl)); 322 PetscCall(ISLocalToGlobalMappingRestoreIndices(l2g, &idxs)); 323 PetscCall(PetscObjectCompose((PetscObject)*ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject)pJ)); 324 PetscCall(DMDestroy(&odm)); 325 PetscCall(MatDestroy(&pJ)); 326 327 /* special purpose setup function (composed in DMPlexSetSNESLocalFEM) */ 328 PetscCall(PetscObjectQueryFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", setup)); 329 if (*setup) PetscCall(PetscObjectCompose((PetscObject)*ovl, "_DM_Original_HPDDM", (PetscObject)dm)); 330 PetscFunctionReturn(PETSC_SUCCESS); 331 } 332