Actual source code: dmi.c

  1: #include <petsc/private/dmimpl.h>
  2: #include <petscds.h>

  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: }

 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];

 21:   PetscFunctionBegin;
 22:   PetscCall(DMGetGlobalSection(dm, &gSection));
 23:   PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd));
 24:   for (p = pStart; p < pEnd; ++p) {
 25:     PetscInt dof, cdof;

 27:     PetscCall(PetscSectionGetDof(gSection, p, &dof));
 28:     PetscCall(PetscSectionGetConstraintDof(gSection, p, &cdof));

 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:   }

 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;

 48:   if (bs < 0) { /* Everyone was empty */
 49:     blockSize = 1;
 50:     bs        = 1;
 51:   }

 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: }

 64: PetscErrorCode DMCreateLocalVector_Section_Private(DM dm, Vec *vec)
 65: {
 66:   PetscSection section;
 67:   PetscInt     localSize, blockSize = -1, pStart, pEnd, p;

 69:   PetscFunctionBegin;
 70:   PetscCall(DMGetLocalSection(dm, &section));
 71:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
 72:   for (p = pStart; p < pEnd; ++p) {
 73:     PetscInt dof;

 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: }

 88: /*@C
 89:   DMCreateSectionSubDM - Returns an `IS` and subDM+subSection encapsulating a subproblem defined by the fields in a `PetscSection` in the `DM`.

 91:   Not Collective

 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

 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`

102:   Level: intermediate

104:   Note:
105:   This handles all information in the `DM` class and the `PetscSection`. This is used as the basis for creating subDMs in specialized classes,
106:   such as `DMPLEX` and `DMFOREST`

108: .seealso `DMCreateSubDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
109: @*/
110: PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
111: {
112:   PetscSection section, sectionGlobal;
113:   PetscInt    *subIndices;
114:   PetscInt     subSize = 0, subOff = 0, Nf, f, pStart, pEnd, p;

116:   PetscFunctionBegin;
117:   if (!numFields) PetscFunctionReturn(PETSC_SUCCESS);
118:   PetscCall(DMGetLocalSection(dm, &section));
119:   PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
120:   PetscCheck(section, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
121:   PetscCheck(sectionGlobal, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
122:   PetscCall(PetscSectionGetNumFields(section, &Nf));
123:   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);
124:   if (is) {
125:     PetscInt bs, bsLocal[2], bsMinMax[2];

127:     for (f = 0, bs = 0; f < numFields; ++f) {
128:       PetscInt Nc;

130:       PetscCall(PetscSectionGetFieldComponents(section, fields[f], &Nc));
131:       bs += Nc;
132:     }
133:     PetscCall(PetscSectionGetChart(sectionGlobal, &pStart, &pEnd));
134:     for (p = pStart; p < pEnd; ++p) {
135:       PetscInt gdof, pSubSize = 0;

137:       PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
138:       if (gdof > 0) {
139:         for (f = 0; f < numFields; ++f) {
140:           PetscInt fdof, fcdof;

142:           PetscCall(PetscSectionGetFieldDof(section, p, fields[f], &fdof));
143:           PetscCall(PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof));
144:           pSubSize += fdof - fcdof;
145:         }
146:         subSize += pSubSize;
147:         if (pSubSize && bs != pSubSize) {
148:           /* Layout does not admit a pointwise block size */
149:           bs = 1;
150:         }
151:       }
152:     }
153:     /* Must have same blocksize on all procs (some might have no points) */
154:     bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs;
155:     bsLocal[1] = bs;
156:     PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)dm), bsLocal, bsMinMax));
157:     if (bsMinMax[0] != bsMinMax[1]) {
158:       bs = 1;
159:     } else {
160:       bs = bsMinMax[0];
161:     }
162:     PetscCall(PetscMalloc1(subSize, &subIndices));
163:     for (p = pStart; p < pEnd; ++p) {
164:       PetscInt gdof, goff;

166:       PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
167:       if (gdof > 0) {
168:         PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
169:         for (f = 0; f < numFields; ++f) {
170:           PetscInt fdof, fcdof, fc, f2, poff = 0;

172:           /* Can get rid of this loop by storing field information in the global section */
173:           for (f2 = 0; f2 < fields[f]; ++f2) {
174:             PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof));
175:             PetscCall(PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof));
176:             poff += fdof - fcdof;
177:           }
178:           PetscCall(PetscSectionGetFieldDof(section, p, fields[f], &fdof));
179:           PetscCall(PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof));
180:           for (fc = 0; fc < fdof - fcdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
181:         }
182:       }
183:     }
184:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is));
185:     if (bs > 1) {
186:       /* We need to check that the block size does not come from non-contiguous fields */
187:       PetscInt i, j, set = 1, rset = 1;
188:       for (i = 0; i < subSize; i += bs) {
189:         for (j = 0; j < bs; ++j) {
190:           if (subIndices[i + j] != subIndices[i] + j) {
191:             set = 0;
192:             break;
193:           }
194:         }
195:       }
196:       PetscCall(MPIU_Allreduce(&set, &rset, 1, MPIU_INT, MPI_PROD, PetscObjectComm((PetscObject)dm)));
197:       if (rset) PetscCall(ISSetBlockSize(*is, bs));
198:     }
199:   }
200:   if (subdm) {
201:     PetscSection subsection;
202:     PetscBool    haveNull = PETSC_FALSE;
203:     PetscInt     f, nf = 0, of = 0;

205:     PetscCall(PetscSectionCreateSubsection(section, numFields, fields, &subsection));
206:     PetscCall(DMSetLocalSection(*subdm, subsection));
207:     PetscCall(PetscSectionDestroy(&subsection));
208:     for (f = 0; f < numFields; ++f) {
209:       (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
210:       if ((*subdm)->nullspaceConstructors[f]) {
211:         haveNull = PETSC_TRUE;
212:         nf       = f;
213:         of       = fields[f];
214:       }
215:     }
216:     if (dm->probs) {
217:       PetscCall(DMSetNumFields(*subdm, numFields));
218:       for (f = 0; f < numFields; ++f) {
219:         PetscObject disc;

221:         PetscCall(DMGetField(dm, fields[f], NULL, &disc));
222:         PetscCall(DMSetField(*subdm, f, NULL, disc));
223:       }
224:       PetscCall(DMCreateDS(*subdm));
225:       if (numFields == 1 && is) {
226:         PetscObject disc, space, pmat;

228:         PetscCall(DMGetField(*subdm, 0, NULL, &disc));
229:         PetscCall(PetscObjectQuery(disc, "nullspace", &space));
230:         if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", space));
231:         PetscCall(PetscObjectQuery(disc, "nearnullspace", &space));
232:         if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nearnullspace", space));
233:         PetscCall(PetscObjectQuery(disc, "pmat", &pmat));
234:         if (pmat) PetscCall(PetscObjectCompose((PetscObject)*is, "pmat", pmat));
235:       }
236:       /* Check if DSes record their DM fields */
237:       if (dm->probs[0].fields) {
238:         PetscInt d, e;

240:         for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) {
241:           const PetscInt  Nf = dm->probs[d].ds->Nf;
242:           const PetscInt *fld;
243:           PetscInt        f, g;

245:           PetscCall(ISGetIndices(dm->probs[d].fields, &fld));
246:           for (f = 0; f < Nf; ++f) {
247:             for (g = 0; g < numFields; ++g)
248:               if (fld[f] == fields[g]) break;
249:             if (g < numFields) break;
250:           }
251:           PetscCall(ISRestoreIndices(dm->probs[d].fields, &fld));
252:           if (f == Nf) continue;
253:           PetscCall(PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds));
254:           PetscCall(PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds));
255:           /* Translate DM fields to DS fields */
256:           {
257:             IS              infields, dsfields;
258:             const PetscInt *fld, *ofld;
259:             PetscInt       *fidx;
260:             PetscInt        onf, nf, f, g;

262:             PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields));
263:             PetscCall(ISIntersect(infields, dm->probs[d].fields, &dsfields));
264:             PetscCall(ISDestroy(&infields));
265:             PetscCall(ISGetLocalSize(dsfields, &nf));
266:             PetscCheck(nf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields");
267:             PetscCall(ISGetIndices(dsfields, &fld));
268:             PetscCall(ISGetLocalSize(dm->probs[d].fields, &onf));
269:             PetscCall(ISGetIndices(dm->probs[d].fields, &ofld));
270:             PetscCall(PetscMalloc1(nf, &fidx));
271:             for (f = 0, g = 0; f < onf && g < nf; ++f) {
272:               if (ofld[f] == fld[g]) fidx[g++] = f;
273:             }
274:             PetscCall(ISRestoreIndices(dm->probs[d].fields, &ofld));
275:             PetscCall(ISRestoreIndices(dsfields, &fld));
276:             PetscCall(ISDestroy(&dsfields));
277:             PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds));
278:             PetscCall(PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds));
279:             PetscCall(PetscFree(fidx));
280:           }
281:           ++e;
282:         }
283:       } else {
284:         PetscCall(PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds));
285:         PetscCall(PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds));
286:         PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds));
287:         PetscCall(PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds));
288:       }
289:     }
290:     if (haveNull && is) {
291:       MatNullSpace nullSpace;

293:       PetscCall((*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace));
294:       PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace));
295:       PetscCall(MatNullSpaceDestroy(&nullSpace));
296:     }
297:     if (dm->coarseMesh) PetscCall(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh));
298:   }
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

302: /*@C
303:   DMCreateSectionSuperDM - Returns an arrays of `IS` and DM+Section encapsulating a superproblem defined by the DM+Sections passed in.

305:   Not Collective

307:   Input Parameters:
308: + dms - The `DM` objects
309: - len - The number of `DM`s

311:   Output Parameters:
312: + is - The global indices for the subproblem, or `NULL`
313: - superdm - The `DM` for the superproblem, which must already have be cloned

315:   Level: intermediate

317:   Note:
318:   This handles all information in the `DM` class and the `PetscSection`. This is used as the basis for creating subDMs in specialized classes,
319:   such as `DMPLEX` and `DMFOREST`

321: .seealso `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
322: @*/
323: PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
324: {
325:   MPI_Comm     comm;
326:   PetscSection supersection, *sections, *sectionGlobals;
327:   PetscInt    *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i;
328:   PetscBool    haveNull = PETSC_FALSE;

330:   PetscFunctionBegin;
331:   PetscCall(PetscObjectGetComm((PetscObject)dms[0], &comm));
332:   /* Pull out local and global sections */
333:   PetscCall(PetscMalloc3(len, &Nfs, len, &sections, len, &sectionGlobals));
334:   for (i = 0; i < len; ++i) {
335:     PetscCall(DMGetLocalSection(dms[i], &sections[i]));
336:     PetscCall(DMGetGlobalSection(dms[i], &sectionGlobals[i]));
337:     PetscCheck(sections[i], comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
338:     PetscCheck(sectionGlobals[i], comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
339:     PetscCall(PetscSectionGetNumFields(sections[i], &Nfs[i]));
340:     Nf += Nfs[i];
341:   }
342:   /* Create the supersection */
343:   PetscCall(PetscSectionCreateSupersection(sections, len, &supersection));
344:   PetscCall(DMSetLocalSection(*superdm, supersection));
345:   /* Create ISes */
346:   if (is) {
347:     PetscSection supersectionGlobal;
348:     PetscInt     bs = -1, startf = 0;

350:     PetscCall(PetscMalloc1(len, is));
351:     PetscCall(DMGetGlobalSection(*superdm, &supersectionGlobal));
352:     for (i = 0; i < len; startf += Nfs[i], ++i) {
353:       PetscInt *subIndices;
354:       PetscInt  subSize, subOff, pStart, pEnd, p, start, end, dummy;

356:       PetscCall(PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd));
357:       PetscCall(PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize));
358:       PetscCall(PetscMalloc1(subSize, &subIndices));
359:       for (p = pStart, subOff = 0; p < pEnd; ++p) {
360:         PetscInt gdof, gcdof, gtdof, d;

362:         PetscCall(PetscSectionGetDof(sectionGlobals[i], p, &gdof));
363:         PetscCall(PetscSectionGetConstraintDof(sections[i], p, &gcdof));
364:         gtdof = gdof - gcdof;
365:         if (gdof > 0 && gtdof) {
366:           if (bs < 0) {
367:             bs = gtdof;
368:           } else if (bs != gtdof) {
369:             bs = 1;
370:           }
371:           PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy));
372:           PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf + Nfs[i] - 1, &dummy, &end));
373:           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);
374:           for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d;
375:         }
376:       }
377:       PetscCall(ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]));
378:       /* Must have same blocksize on all procs (some might have no points) */
379:       {
380:         PetscInt bs = -1, bsLocal[2], bsMinMax[2];

382:         bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs;
383:         bsLocal[1] = bs;
384:         PetscCall(PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax));
385:         if (bsMinMax[0] != bsMinMax[1]) {
386:           bs = 1;
387:         } else {
388:           bs = bsMinMax[0];
389:         }
390:         PetscCall(ISSetBlockSize((*is)[i], bs));
391:       }
392:     }
393:   }
394:   /* Preserve discretizations */
395:   if (len && dms[0]->probs) {
396:     PetscCall(DMSetNumFields(*superdm, Nf));
397:     for (i = 0, supf = 0; i < len; ++i) {
398:       for (f = 0; f < Nfs[i]; ++f, ++supf) {
399:         PetscObject disc;

401:         PetscCall(DMGetField(dms[i], f, NULL, &disc));
402:         PetscCall(DMSetField(*superdm, supf, NULL, disc));
403:       }
404:     }
405:     PetscCall(DMCreateDS(*superdm));
406:   }
407:   /* Preserve nullspaces */
408:   for (i = 0, supf = 0; i < len; ++i) {
409:     for (f = 0; f < Nfs[i]; ++f, ++supf) {
410:       (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f];
411:       if ((*superdm)->nullspaceConstructors[supf]) {
412:         haveNull = PETSC_TRUE;
413:         nullf    = supf;
414:         oldf     = f;
415:       }
416:     }
417:   }
418:   /* Attach nullspace to IS */
419:   if (haveNull && is) {
420:     MatNullSpace nullSpace;

422:     PetscCall((*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace));
423:     PetscCall(PetscObjectCompose((PetscObject)(*is)[nullf], "nullspace", (PetscObject)nullSpace));
424:     PetscCall(MatNullSpaceDestroy(&nullSpace));
425:   }
426:   PetscCall(PetscSectionDestroy(&supersection));
427:   PetscCall(PetscFree3(Nfs, sections, sectionGlobals));
428:   PetscFunctionReturn(PETSC_SUCCESS);
429: }