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` containing a `PetscSection` that encapsulates a subproblem defined by a subset of 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 to incorporate into `subdm` 96 - fields - The field numbers of the selected fields 97 98 Output Parameters: 99 + is - The global indices for the subproblem or `NULL` 100 - subdm - The `DM` for the subproblem, which must already have be cloned from `dm` or `NULL` 101 102 Level: intermediate 103 104 Note: 105 If `is` and `subdm` are both `NULL` this does nothing 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 PetscCall(MPIU_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 `IS` and a `DM` containing a `PetscSection` that encapsulates a superproblem defined by the array of `DM` and their `PetscSection` 303 304 Not Collective 305 306 Input Parameters: 307 + dms - The `DM` objects, the must all have the same topology; for example obtained with `DMClone()` 308 - len - The number of `DM` in `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 and contain the same topology as the `dms` 313 314 Level: intermediate 315 316 .seealso: `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()` 317 @*/ 318 PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm) 319 { 320 MPI_Comm comm; 321 PetscSection supersection, *sections, *sectionGlobals; 322 PetscInt *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i; 323 PetscBool haveNull = PETSC_FALSE; 324 325 PetscFunctionBegin; 326 PetscCall(PetscObjectGetComm((PetscObject)dms[0], &comm)); 327 /* Pull out local and global sections */ 328 PetscCall(PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals)); 329 for (i = 0; i < len; ++i) { 330 PetscCall(DMGetLocalSection(dms[i], §ions[i])); 331 PetscCall(DMGetGlobalSection(dms[i], §ionGlobals[i])); 332 PetscCheck(sections[i], comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 333 PetscCheck(sectionGlobals[i], comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 334 PetscCall(PetscSectionGetNumFields(sections[i], &Nfs[i])); 335 Nf += Nfs[i]; 336 } 337 /* Create the supersection */ 338 PetscCall(PetscSectionCreateSupersection(sections, len, &supersection)); 339 PetscCall(DMSetLocalSection(*superdm, supersection)); 340 /* Create ISes */ 341 if (is) { 342 PetscSection supersectionGlobal; 343 PetscInt bs = -1, startf = 0; 344 345 PetscCall(PetscMalloc1(len, is)); 346 PetscCall(DMGetGlobalSection(*superdm, &supersectionGlobal)); 347 for (i = 0; i < len; startf += Nfs[i], ++i) { 348 PetscInt *subIndices; 349 PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy; 350 351 PetscCall(PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd)); 352 PetscCall(PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize)); 353 PetscCall(PetscMalloc1(subSize, &subIndices)); 354 for (p = pStart, subOff = 0; p < pEnd; ++p) { 355 PetscInt gdof, gcdof, gtdof, d; 356 357 PetscCall(PetscSectionGetDof(sectionGlobals[i], p, &gdof)); 358 PetscCall(PetscSectionGetConstraintDof(sections[i], p, &gcdof)); 359 gtdof = gdof - gcdof; 360 if (gdof > 0 && gtdof) { 361 if (bs < 0) { 362 bs = gtdof; 363 } else if (bs != gtdof) { 364 bs = 1; 365 } 366 PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy)); 367 PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf + Nfs[i] - 1, &dummy, &end)); 368 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); 369 for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d; 370 } 371 } 372 PetscCall(ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i])); 373 /* Must have same blocksize on all procs (some might have no points) */ 374 { 375 PetscInt bs = -1, bsLocal[2], bsMinMax[2]; 376 377 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; 378 bsLocal[1] = bs; 379 PetscCall(PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax)); 380 if (bsMinMax[0] != bsMinMax[1]) { 381 bs = 1; 382 } else { 383 bs = bsMinMax[0]; 384 } 385 PetscCall(ISSetBlockSize((*is)[i], bs)); 386 } 387 } 388 } 389 /* Preserve discretizations */ 390 if (len && dms[0]->probs) { 391 PetscCall(DMSetNumFields(*superdm, Nf)); 392 for (i = 0, supf = 0; i < len; ++i) { 393 for (f = 0; f < Nfs[i]; ++f, ++supf) { 394 PetscObject disc; 395 396 PetscCall(DMGetField(dms[i], f, NULL, &disc)); 397 PetscCall(DMSetField(*superdm, supf, NULL, disc)); 398 } 399 } 400 PetscCall(DMCreateDS(*superdm)); 401 } 402 /* Preserve nullspaces */ 403 for (i = 0, supf = 0; i < len; ++i) { 404 for (f = 0; f < Nfs[i]; ++f, ++supf) { 405 (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f]; 406 if ((*superdm)->nullspaceConstructors[supf]) { 407 haveNull = PETSC_TRUE; 408 nullf = supf; 409 oldf = f; 410 } 411 } 412 } 413 /* Attach nullspace to IS */ 414 if (haveNull && is) { 415 MatNullSpace nullSpace; 416 417 PetscCall((*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace)); 418 PetscCall(PetscObjectCompose((PetscObject)(*is)[nullf], "nullspace", (PetscObject)nullSpace)); 419 PetscCall(MatNullSpaceDestroy(&nullSpace)); 420 } 421 PetscCall(PetscSectionDestroy(&supersection)); 422 PetscCall(PetscFree3(Nfs, sections, sectionGlobals)); 423 PetscFunctionReturn(PETSC_SUCCESS); 424 } 425