1 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/ 2 #include <petscds.h> 3 4 // Greatest common divisor of two nonnegative integers 5 PetscInt PetscGCD(PetscInt a, PetscInt b) 6 { 7 while (b != 0) { 8 PetscInt tmp = b; 9 b = a % b; 10 a = tmp; 11 } 12 return a; 13 } 14 15 PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm, Vec *vec) 16 { 17 PetscSection gSection; 18 PetscInt localSize, bs, blockSize = -1, pStart, pEnd, p; 19 PetscInt in[2], out[2]; 20 21 PetscFunctionBegin; 22 PetscCall(DMGetGlobalSection(dm, &gSection)); 23 PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd)); 24 for (p = pStart; p < pEnd; ++p) { 25 PetscInt dof, cdof; 26 27 PetscCall(PetscSectionGetDof(gSection, p, &dof)); 28 PetscCall(PetscSectionGetConstraintDof(gSection, p, &cdof)); 29 30 if (dof - cdof > 0) { 31 if (blockSize < 0) { 32 /* set blockSize */ 33 blockSize = dof - cdof; 34 } else { 35 blockSize = PetscGCD(dof - cdof, blockSize); 36 } 37 } 38 } 39 40 in[0] = blockSize < 0 ? PETSC_MIN_INT : -blockSize; 41 in[1] = blockSize; 42 PetscCall(MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm))); 43 /* -out[0] = min(blockSize), out[1] = max(blockSize) */ 44 if (-out[0] == out[1]) { 45 bs = out[1]; 46 } else bs = 1; 47 48 if (bs < 0) { /* Everyone was empty */ 49 blockSize = 1; 50 bs = 1; 51 } 52 53 PetscCall(PetscSectionGetConstrainedStorageSize(gSection, &localSize)); 54 PetscCheck(localSize % blockSize == 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %" PetscInt_FMT " and local storage size %" PetscInt_FMT, blockSize, localSize); 55 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), vec)); 56 PetscCall(VecSetSizes(*vec, localSize, PETSC_DETERMINE)); 57 PetscCall(VecSetBlockSize(*vec, bs)); 58 PetscCall(VecSetType(*vec, dm->vectype)); 59 PetscCall(VecSetDM(*vec, dm)); 60 /* PetscCall(VecSetLocalToGlobalMapping(*vec, dm->ltogmap)); */ 61 PetscFunctionReturn(PETSC_SUCCESS); 62 } 63 64 PetscErrorCode DMCreateLocalVector_Section_Private(DM dm, Vec *vec) 65 { 66 PetscSection section; 67 PetscInt localSize, blockSize = -1, pStart, pEnd, p; 68 69 PetscFunctionBegin; 70 PetscCall(DMGetLocalSection(dm, §ion)); 71 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 72 for (p = pStart; p < pEnd; ++p) { 73 PetscInt dof; 74 75 PetscCall(PetscSectionGetDof(section, p, &dof)); 76 if ((blockSize < 0) && (dof > 0)) blockSize = dof; 77 if (dof > 0) blockSize = PetscGCD(dof, blockSize); 78 } 79 PetscCall(PetscSectionGetStorageSize(section, &localSize)); 80 PetscCall(VecCreate(PETSC_COMM_SELF, vec)); 81 PetscCall(VecSetSizes(*vec, localSize, localSize)); 82 PetscCall(VecSetBlockSize(*vec, blockSize)); 83 PetscCall(VecSetType(*vec, dm->vectype)); 84 PetscCall(VecSetDM(*vec, dm)); 85 PetscFunctionReturn(PETSC_SUCCESS); 86 } 87 88 /*@C 89 DMCreateSectionSubDM - Returns an `IS` and subDM+subSection encapsulating a subproblem defined by the fields in a `PetscSection` in the `DM`. 90 91 Not Collective 92 93 Input Parameters: 94 + dm - The `DM` object 95 . numFields - The number of fields in this subproblem 96 - fields - The field numbers of the selected fields 97 98 Output Parameters: 99 + is - The global indices for the subproblem 100 - subdm - The `DM` for the subproblem, which must already have be cloned from `dm` 101 102 Level: intermediate 103 104 Note: 105 This handles all information in the `DM` class and the `PetscSection`. This is used as the basis for creating subDMs in specialized classes, 106 such as `DMPLEX` and `DMFOREST` 107 108 .seealso `DMCreateSubDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()` 109 @*/ 110 PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm) 111 { 112 PetscSection section, sectionGlobal; 113 PetscInt *subIndices; 114 PetscInt subSize = 0, subOff = 0, Nf, f, pStart, pEnd, p; 115 116 PetscFunctionBegin; 117 if (!numFields) PetscFunctionReturn(PETSC_SUCCESS); 118 PetscCall(DMGetLocalSection(dm, §ion)); 119 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 120 PetscCheck(section, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 121 PetscCheck(sectionGlobal, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 122 PetscCall(PetscSectionGetNumFields(section, &Nf)); 123 PetscCheck(numFields <= Nf, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Number of requested fields %" PetscInt_FMT " greater than number of DM fields %" PetscInt_FMT, numFields, Nf); 124 if (is) { 125 PetscInt bs, bsLocal[2], bsMinMax[2]; 126 127 for (f = 0, bs = 0; f < numFields; ++f) { 128 PetscInt Nc; 129 130 PetscCall(PetscSectionGetFieldComponents(section, fields[f], &Nc)); 131 bs += Nc; 132 } 133 PetscCall(PetscSectionGetChart(sectionGlobal, &pStart, &pEnd)); 134 for (p = pStart; p < pEnd; ++p) { 135 PetscInt gdof, pSubSize = 0; 136 137 PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof)); 138 if (gdof > 0) { 139 for (f = 0; f < numFields; ++f) { 140 PetscInt fdof, fcdof; 141 142 PetscCall(PetscSectionGetFieldDof(section, p, fields[f], &fdof)); 143 PetscCall(PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof)); 144 pSubSize += fdof - fcdof; 145 } 146 subSize += pSubSize; 147 if (pSubSize && bs != pSubSize) { 148 /* Layout does not admit a pointwise block size */ 149 bs = 1; 150 } 151 } 152 } 153 /* Must have same blocksize on all procs (some might have no points) */ 154 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; 155 bsLocal[1] = bs; 156 PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)dm), bsLocal, bsMinMax)); 157 if (bsMinMax[0] != bsMinMax[1]) { 158 bs = 1; 159 } else { 160 bs = bsMinMax[0]; 161 } 162 PetscCall(PetscMalloc1(subSize, &subIndices)); 163 for (p = pStart; p < pEnd; ++p) { 164 PetscInt gdof, goff; 165 166 PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof)); 167 if (gdof > 0) { 168 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 169 for (f = 0; f < numFields; ++f) { 170 PetscInt fdof, fcdof, fc, f2, poff = 0; 171 172 /* Can get rid of this loop by storing field information in the global section */ 173 for (f2 = 0; f2 < fields[f]; ++f2) { 174 PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof)); 175 PetscCall(PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof)); 176 poff += fdof - fcdof; 177 } 178 PetscCall(PetscSectionGetFieldDof(section, p, fields[f], &fdof)); 179 PetscCall(PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof)); 180 for (fc = 0; fc < fdof - fcdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc; 181 } 182 } 183 } 184 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is)); 185 if (bs > 1) { 186 /* We need to check that the block size does not come from non-contiguous fields */ 187 PetscInt i, j, set = 1, rset = 1; 188 for (i = 0; i < subSize; i += bs) { 189 for (j = 0; j < bs; ++j) { 190 if (subIndices[i + j] != subIndices[i] + j) { 191 set = 0; 192 break; 193 } 194 } 195 } 196 PetscCallMPI(MPI_Allreduce(&set, &rset, 1, MPIU_INT, MPI_PROD, PetscObjectComm((PetscObject)dm))); 197 if (rset) PetscCall(ISSetBlockSize(*is, bs)); 198 } 199 } 200 if (subdm) { 201 PetscSection subsection; 202 PetscBool haveNull = PETSC_FALSE; 203 PetscInt f, nf = 0, of = 0; 204 205 PetscCall(PetscSectionCreateSubsection(section, numFields, fields, &subsection)); 206 PetscCall(DMSetLocalSection(*subdm, subsection)); 207 PetscCall(PetscSectionDestroy(&subsection)); 208 for (f = 0; f < numFields; ++f) { 209 (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]]; 210 if ((*subdm)->nullspaceConstructors[f]) { 211 haveNull = PETSC_TRUE; 212 nf = f; 213 of = fields[f]; 214 } 215 } 216 if (dm->probs) { 217 PetscCall(DMSetNumFields(*subdm, numFields)); 218 for (f = 0; f < numFields; ++f) { 219 PetscObject disc; 220 221 PetscCall(DMGetField(dm, fields[f], NULL, &disc)); 222 PetscCall(DMSetField(*subdm, f, NULL, disc)); 223 } 224 PetscCall(DMCreateDS(*subdm)); 225 if (numFields == 1 && is) { 226 PetscObject disc, space, pmat; 227 228 PetscCall(DMGetField(*subdm, 0, NULL, &disc)); 229 PetscCall(PetscObjectQuery(disc, "nullspace", &space)); 230 if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", space)); 231 PetscCall(PetscObjectQuery(disc, "nearnullspace", &space)); 232 if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nearnullspace", space)); 233 PetscCall(PetscObjectQuery(disc, "pmat", &pmat)); 234 if (pmat) PetscCall(PetscObjectCompose((PetscObject)*is, "pmat", pmat)); 235 } 236 /* Check if DSes record their DM fields */ 237 if (dm->probs[0].fields) { 238 PetscInt d, e; 239 240 for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) { 241 const PetscInt Nf = dm->probs[d].ds->Nf; 242 const PetscInt *fld; 243 PetscInt f, g; 244 245 PetscCall(ISGetIndices(dm->probs[d].fields, &fld)); 246 for (f = 0; f < Nf; ++f) { 247 for (g = 0; g < numFields; ++g) 248 if (fld[f] == fields[g]) break; 249 if (g < numFields) break; 250 } 251 PetscCall(ISRestoreIndices(dm->probs[d].fields, &fld)); 252 if (f == Nf) continue; 253 PetscCall(PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds)); 254 PetscCall(PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds)); 255 /* Translate DM fields to DS fields */ 256 { 257 IS infields, dsfields; 258 const PetscInt *fld, *ofld; 259 PetscInt *fidx; 260 PetscInt onf, nf, f, g; 261 262 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields)); 263 PetscCall(ISIntersect(infields, dm->probs[d].fields, &dsfields)); 264 PetscCall(ISDestroy(&infields)); 265 PetscCall(ISGetLocalSize(dsfields, &nf)); 266 PetscCheck(nf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields"); 267 PetscCall(ISGetIndices(dsfields, &fld)); 268 PetscCall(ISGetLocalSize(dm->probs[d].fields, &onf)); 269 PetscCall(ISGetIndices(dm->probs[d].fields, &ofld)); 270 PetscCall(PetscMalloc1(nf, &fidx)); 271 for (f = 0, g = 0; f < onf && g < nf; ++f) { 272 if (ofld[f] == fld[g]) fidx[g++] = f; 273 } 274 PetscCall(ISRestoreIndices(dm->probs[d].fields, &ofld)); 275 PetscCall(ISRestoreIndices(dsfields, &fld)); 276 PetscCall(ISDestroy(&dsfields)); 277 PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds)); 278 PetscCall(PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds)); 279 PetscCall(PetscFree(fidx)); 280 } 281 ++e; 282 } 283 } else { 284 PetscCall(PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds)); 285 PetscCall(PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds)); 286 PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds)); 287 PetscCall(PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds)); 288 } 289 } 290 if (haveNull && is) { 291 MatNullSpace nullSpace; 292 293 PetscCall((*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace)); 294 PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace)); 295 PetscCall(MatNullSpaceDestroy(&nullSpace)); 296 } 297 if (dm->coarseMesh) PetscCall(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh)); 298 } 299 PetscFunctionReturn(PETSC_SUCCESS); 300 } 301 302 /*@C 303 DMCreateSectionSuperDM - Returns an arrays of `IS` and DM+Section encapsulating a superproblem defined by the DM+Sections passed in. 304 305 Not Collective 306 307 Input Parameters: 308 + dms - The `DM` objects 309 - len - The number of `DM`s 310 311 Output Parameters: 312 + is - The global indices for the subproblem, or `NULL` 313 - superdm - The `DM` for the superproblem, which must already have be cloned 314 315 Level: intermediate 316 317 Note: 318 This handles all information in the `DM` class and the `PetscSection`. This is used as the basis for creating subDMs in specialized classes, 319 such as `DMPLEX` and `DMFOREST` 320 321 .seealso `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()` 322 @*/ 323 PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm) 324 { 325 MPI_Comm comm; 326 PetscSection supersection, *sections, *sectionGlobals; 327 PetscInt *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i; 328 PetscBool haveNull = PETSC_FALSE; 329 330 PetscFunctionBegin; 331 PetscCall(PetscObjectGetComm((PetscObject)dms[0], &comm)); 332 /* Pull out local and global sections */ 333 PetscCall(PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals)); 334 for (i = 0; i < len; ++i) { 335 PetscCall(DMGetLocalSection(dms[i], §ions[i])); 336 PetscCall(DMGetGlobalSection(dms[i], §ionGlobals[i])); 337 PetscCheck(sections[i], comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 338 PetscCheck(sectionGlobals[i], comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 339 PetscCall(PetscSectionGetNumFields(sections[i], &Nfs[i])); 340 Nf += Nfs[i]; 341 } 342 /* Create the supersection */ 343 PetscCall(PetscSectionCreateSupersection(sections, len, &supersection)); 344 PetscCall(DMSetLocalSection(*superdm, supersection)); 345 /* Create ISes */ 346 if (is) { 347 PetscSection supersectionGlobal; 348 PetscInt bs = -1, startf = 0; 349 350 PetscCall(PetscMalloc1(len, is)); 351 PetscCall(DMGetGlobalSection(*superdm, &supersectionGlobal)); 352 for (i = 0; i < len; startf += Nfs[i], ++i) { 353 PetscInt *subIndices; 354 PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy; 355 356 PetscCall(PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd)); 357 PetscCall(PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize)); 358 PetscCall(PetscMalloc1(subSize, &subIndices)); 359 for (p = pStart, subOff = 0; p < pEnd; ++p) { 360 PetscInt gdof, gcdof, gtdof, d; 361 362 PetscCall(PetscSectionGetDof(sectionGlobals[i], p, &gdof)); 363 PetscCall(PetscSectionGetConstraintDof(sections[i], p, &gcdof)); 364 gtdof = gdof - gcdof; 365 if (gdof > 0 && gtdof) { 366 if (bs < 0) { 367 bs = gtdof; 368 } else if (bs != gtdof) { 369 bs = 1; 370 } 371 PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy)); 372 PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf + Nfs[i] - 1, &dummy, &end)); 373 PetscCheck(end - start == gtdof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid number of global dofs %" PetscInt_FMT " != %" PetscInt_FMT " for dm %" PetscInt_FMT " on point %" PetscInt_FMT, end - start, gtdof, i, p); 374 for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d; 375 } 376 } 377 PetscCall(ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i])); 378 /* Must have same blocksize on all procs (some might have no points) */ 379 { 380 PetscInt bs = -1, bsLocal[2], bsMinMax[2]; 381 382 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; 383 bsLocal[1] = bs; 384 PetscCall(PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax)); 385 if (bsMinMax[0] != bsMinMax[1]) { 386 bs = 1; 387 } else { 388 bs = bsMinMax[0]; 389 } 390 PetscCall(ISSetBlockSize((*is)[i], bs)); 391 } 392 } 393 } 394 /* Preserve discretizations */ 395 if (len && dms[0]->probs) { 396 PetscCall(DMSetNumFields(*superdm, Nf)); 397 for (i = 0, supf = 0; i < len; ++i) { 398 for (f = 0; f < Nfs[i]; ++f, ++supf) { 399 PetscObject disc; 400 401 PetscCall(DMGetField(dms[i], f, NULL, &disc)); 402 PetscCall(DMSetField(*superdm, supf, NULL, disc)); 403 } 404 } 405 PetscCall(DMCreateDS(*superdm)); 406 } 407 /* Preserve nullspaces */ 408 for (i = 0, supf = 0; i < len; ++i) { 409 for (f = 0; f < Nfs[i]; ++f, ++supf) { 410 (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f]; 411 if ((*superdm)->nullspaceConstructors[supf]) { 412 haveNull = PETSC_TRUE; 413 nullf = supf; 414 oldf = f; 415 } 416 } 417 } 418 /* Attach nullspace to IS */ 419 if (haveNull && is) { 420 MatNullSpace nullSpace; 421 422 PetscCall((*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace)); 423 PetscCall(PetscObjectCompose((PetscObject)(*is)[nullf], "nullspace", (PetscObject)nullSpace)); 424 PetscCall(MatNullSpaceDestroy(&nullSpace)); 425 } 426 PetscCall(PetscSectionDestroy(&supersection)); 427 PetscCall(PetscFree3(Nfs, sections, sectionGlobals)); 428 PetscFunctionReturn(PETSC_SUCCESS); 429 } 430