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, §ion));
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, §ion));
119: PetscCall(DMGetGlobalSection(dm, §ionGlobal));
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, §ions, len, §ionGlobals));
334: for (i = 0; i < len; ++i) {
335: PetscCall(DMGetLocalSection(dms[i], §ions[i]));
336: PetscCall(DMGetGlobalSection(dms[i], §ionGlobals[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: }