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