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; 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 if (set) PetscCall(ISSetBlockSize(*is, bs)); 190 } 191 } 192 if (subdm) { 193 PetscSection subsection; 194 PetscBool haveNull = PETSC_FALSE; 195 PetscInt f, nf = 0, of = 0; 196 197 PetscCall(PetscSectionCreateSubsection(section, numFields, fields, &subsection)); 198 PetscCall(DMSetLocalSection(*subdm, subsection)); 199 PetscCall(PetscSectionDestroy(&subsection)); 200 for (f = 0; f < numFields; ++f) { 201 (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]]; 202 if ((*subdm)->nullspaceConstructors[f]) { 203 haveNull = PETSC_TRUE; 204 nf = f; 205 of = fields[f]; 206 } 207 } 208 if (dm->probs) { 209 PetscCall(DMSetNumFields(*subdm, numFields)); 210 for (f = 0; f < numFields; ++f) { 211 PetscObject disc; 212 213 PetscCall(DMGetField(dm, fields[f], NULL, &disc)); 214 PetscCall(DMSetField(*subdm, f, NULL, disc)); 215 } 216 PetscCall(DMCreateDS(*subdm)); 217 if (numFields == 1 && is) { 218 PetscObject disc, space, pmat; 219 220 PetscCall(DMGetField(*subdm, 0, NULL, &disc)); 221 PetscCall(PetscObjectQuery(disc, "nullspace", &space)); 222 if (space) PetscCall(PetscObjectCompose((PetscObject) *is, "nullspace", space)); 223 PetscCall(PetscObjectQuery(disc, "nearnullspace", &space)); 224 if (space) PetscCall(PetscObjectCompose((PetscObject) *is, "nearnullspace", space)); 225 PetscCall(PetscObjectQuery(disc, "pmat", &pmat)); 226 if (pmat) PetscCall(PetscObjectCompose((PetscObject) *is, "pmat", pmat)); 227 } 228 /* Check if DSes record their DM fields */ 229 if (dm->probs[0].fields) { 230 PetscInt d, e; 231 232 for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) { 233 const PetscInt Nf = dm->probs[d].ds->Nf; 234 const PetscInt *fld; 235 PetscInt f, g; 236 237 PetscCall(ISGetIndices(dm->probs[d].fields, &fld)); 238 for (f = 0; f < Nf; ++f) { 239 for (g = 0; g < numFields; ++g) if (fld[f] == fields[g]) break; 240 if (g < numFields) break; 241 } 242 PetscCall(ISRestoreIndices(dm->probs[d].fields, &fld)); 243 if (f == Nf) continue; 244 PetscCall(PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds)); 245 PetscCall(PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds)); 246 /* Translate DM fields to DS fields */ 247 { 248 IS infields, dsfields; 249 const PetscInt *fld, *ofld; 250 PetscInt *fidx; 251 PetscInt onf, nf, f, g; 252 253 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields)); 254 PetscCall(ISIntersect(infields, dm->probs[d].fields, &dsfields)); 255 PetscCall(ISDestroy(&infields)); 256 PetscCall(ISGetLocalSize(dsfields, &nf)); 257 PetscCheck(nf,PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields"); 258 PetscCall(ISGetIndices(dsfields, &fld)); 259 PetscCall(ISGetLocalSize(dm->probs[d].fields, &onf)); 260 PetscCall(ISGetIndices(dm->probs[d].fields, &ofld)); 261 PetscCall(PetscMalloc1(nf, &fidx)); 262 for (f = 0, g = 0; f < onf && g < nf; ++f) { 263 if (ofld[f] == fld[g]) fidx[g++] = f; 264 } 265 PetscCall(ISRestoreIndices(dm->probs[d].fields, &ofld)); 266 PetscCall(ISRestoreIndices(dsfields, &fld)); 267 PetscCall(ISDestroy(&dsfields)); 268 PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds)); 269 PetscCall(PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds)); 270 PetscCall(PetscFree(fidx)); 271 } 272 ++e; 273 } 274 } else { 275 PetscCall(PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds)); 276 PetscCall(PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds)); 277 PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds)); 278 PetscCall(PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds)); 279 } 280 } 281 if (haveNull && is) { 282 MatNullSpace nullSpace; 283 284 PetscCall((*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace)); 285 PetscCall(PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace)); 286 PetscCall(MatNullSpaceDestroy(&nullSpace)); 287 } 288 if (dm->coarseMesh) { 289 PetscCall(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh)); 290 } 291 } 292 PetscFunctionReturn(0); 293 } 294 295 /*@C 296 DMCreateSectionSuperDM - Returns an arrays of ISes and DM+Section encapsulating a superproblem defined by the DM+Sections passed in. 297 298 Not collective 299 300 Input Parameters: 301 + dms - The DM objects 302 - len - The number of DMs 303 304 Output Parameters: 305 + is - The global indices for the subproblem, or NULL 306 - superdm - The DM for the superproblem, which must already have be cloned 307 308 Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes, 309 such as Plex and Forest. 310 311 Level: intermediate 312 313 .seealso `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()` 314 @*/ 315 PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm) 316 { 317 MPI_Comm comm; 318 PetscSection supersection, *sections, *sectionGlobals; 319 PetscInt *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i; 320 PetscBool haveNull = PETSC_FALSE; 321 322 PetscFunctionBegin; 323 PetscCall(PetscObjectGetComm((PetscObject)dms[0], &comm)); 324 /* Pull out local and global sections */ 325 PetscCall(PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals)); 326 for (i = 0 ; i < len; ++i) { 327 PetscCall(DMGetLocalSection(dms[i], §ions[i])); 328 PetscCall(DMGetGlobalSection(dms[i], §ionGlobals[i])); 329 PetscCheck(sections[i],comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 330 PetscCheck(sectionGlobals[i],comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 331 PetscCall(PetscSectionGetNumFields(sections[i], &Nfs[i])); 332 Nf += Nfs[i]; 333 } 334 /* Create the supersection */ 335 PetscCall(PetscSectionCreateSupersection(sections, len, &supersection)); 336 PetscCall(DMSetLocalSection(*superdm, supersection)); 337 /* Create ISes */ 338 if (is) { 339 PetscSection supersectionGlobal; 340 PetscInt bs = -1, startf = 0; 341 342 PetscCall(PetscMalloc1(len, is)); 343 PetscCall(DMGetGlobalSection(*superdm, &supersectionGlobal)); 344 for (i = 0 ; i < len; startf += Nfs[i], ++i) { 345 PetscInt *subIndices; 346 PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy; 347 348 PetscCall(PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd)); 349 PetscCall(PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize)); 350 PetscCall(PetscMalloc1(subSize, &subIndices)); 351 for (p = pStart, subOff = 0; p < pEnd; ++p) { 352 PetscInt gdof, gcdof, gtdof, d; 353 354 PetscCall(PetscSectionGetDof(sectionGlobals[i], p, &gdof)); 355 PetscCall(PetscSectionGetConstraintDof(sections[i], p, &gcdof)); 356 gtdof = gdof - gcdof; 357 if (gdof > 0 && gtdof) { 358 if (bs < 0) {bs = gtdof;} 359 else if (bs != gtdof) {bs = 1;} 360 PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy)); 361 PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf+Nfs[i]-1, &dummy, &end)); 362 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); 363 for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d; 364 } 365 } 366 PetscCall(ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i])); 367 /* Must have same blocksize on all procs (some might have no points) */ 368 { 369 PetscInt bs = -1, bsLocal[2], bsMinMax[2]; 370 371 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs; 372 PetscCall(PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax)); 373 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;} 374 else {bs = bsMinMax[0];} 375 PetscCall(ISSetBlockSize((*is)[i], bs)); 376 } 377 } 378 } 379 /* Preserve discretizations */ 380 if (len && dms[0]->probs) { 381 PetscCall(DMSetNumFields(*superdm, Nf)); 382 for (i = 0, supf = 0; i < len; ++i) { 383 for (f = 0; f < Nfs[i]; ++f, ++supf) { 384 PetscObject disc; 385 386 PetscCall(DMGetField(dms[i], f, NULL, &disc)); 387 PetscCall(DMSetField(*superdm, supf, NULL, disc)); 388 } 389 } 390 PetscCall(DMCreateDS(*superdm)); 391 } 392 /* Preserve nullspaces */ 393 for (i = 0, supf = 0; i < len; ++i) { 394 for (f = 0; f < Nfs[i]; ++f, ++supf) { 395 (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f]; 396 if ((*superdm)->nullspaceConstructors[supf]) { 397 haveNull = PETSC_TRUE; 398 nullf = supf; 399 oldf = f; 400 } 401 } 402 } 403 /* Attach nullspace to IS */ 404 if (haveNull && is) { 405 MatNullSpace nullSpace; 406 407 PetscCall((*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace)); 408 PetscCall(PetscObjectCompose((PetscObject) (*is)[nullf], "nullspace", (PetscObject) nullSpace)); 409 PetscCall(MatNullSpaceDestroy(&nullSpace)); 410 } 411 PetscCall(PetscSectionDestroy(&supersection)); 412 PetscCall(PetscFree3(Nfs, sections, sectionGlobals)); 413 PetscFunctionReturn(0); 414 } 415