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