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