Actual source code: stag.c
1: /*
2: Implementation of DMStag, defining dimension-independent functions in the
3: DM API. stag1d.c, stag2d.c, and stag3d.c may include dimension-specific
4: implementations of DM API functions, and other files here contain additional
5: DMStag-specific API functions, as well as internal functions.
6: */
7: #include <petsc/private/dmstagimpl.h>
8: #include <petscsf.h>
10: static PetscErrorCode DMCreateFieldDecomposition_Stag(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
11: {
12: PetscInt f0, f1, f2, f3, dof0, dof1, dof2, dof3, n_entries, k, d, cnt, n_fields, dim;
13: DMStagStencil *stencil0, *stencil1, *stencil2, *stencil3;
15: PetscFunctionBegin;
16: PetscCall(DMGetDimension(dm, &dim));
17: PetscCall(DMStagGetDOF(dm, &dof0, &dof1, &dof2, &dof3));
18: PetscCall(DMStagGetEntriesPerElement(dm, &n_entries));
20: f0 = 1;
21: f1 = f2 = f3 = 0;
22: if (dim == 1) {
23: f1 = 1;
24: } else if (dim == 2) {
25: f1 = 2;
26: f2 = 1;
27: } else if (dim == 3) {
28: f1 = 3;
29: f2 = 3;
30: f3 = 1;
31: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
33: PetscCall(PetscCalloc1(f0 * dof0, &stencil0));
34: PetscCall(PetscCalloc1(f1 * dof1, &stencil1));
35: if (dim >= 2) PetscCall(PetscCalloc1(f2 * dof2, &stencil2));
36: if (dim >= 3) PetscCall(PetscCalloc1(f3 * dof3, &stencil3));
37: for (k = 0; k < f0; ++k) {
38: for (d = 0; d < dof0; ++d) {
39: stencil0[dof0 * k + d].i = 0;
40: stencil0[dof0 * k + d].j = 0;
41: stencil0[dof0 * k + d].j = 0;
42: }
43: }
44: for (k = 0; k < f1; ++k) {
45: for (d = 0; d < dof1; ++d) {
46: stencil1[dof1 * k + d].i = 0;
47: stencil1[dof1 * k + d].j = 0;
48: stencil1[dof1 * k + d].j = 0;
49: }
50: }
51: if (dim >= 2) {
52: for (k = 0; k < f2; ++k) {
53: for (d = 0; d < dof2; ++d) {
54: stencil2[dof2 * k + d].i = 0;
55: stencil2[dof2 * k + d].j = 0;
56: stencil2[dof2 * k + d].j = 0;
57: }
58: }
59: }
60: if (dim >= 3) {
61: for (k = 0; k < f3; ++k) {
62: for (d = 0; d < dof3; ++d) {
63: stencil3[dof3 * k + d].i = 0;
64: stencil3[dof3 * k + d].j = 0;
65: stencil3[dof3 * k + d].j = 0;
66: }
67: }
68: }
70: n_fields = 0;
71: if (dof0 != 0) ++n_fields;
72: if (dof1 != 0) ++n_fields;
73: if (dim >= 2 && dof2 != 0) ++n_fields;
74: if (dim >= 3 && dof3 != 0) ++n_fields;
75: if (len) *len = n_fields;
77: if (islist) {
78: PetscCall(PetscMalloc1(n_fields, islist));
80: if (dim == 1) {
81: /* face, element */
82: for (d = 0; d < dof0; ++d) {
83: stencil0[d].loc = DMSTAG_LEFT;
84: stencil0[d].c = d;
85: }
86: for (d = 0; d < dof1; ++d) {
87: stencil1[d].loc = DMSTAG_ELEMENT;
88: stencil1[d].c = d;
89: }
90: } else if (dim == 2) {
91: /* vertex, edge(down,left), element */
92: for (d = 0; d < dof0; ++d) {
93: stencil0[d].loc = DMSTAG_DOWN_LEFT;
94: stencil0[d].c = d;
95: }
96: /* edge */
97: cnt = 0;
98: for (d = 0; d < dof1; ++d) {
99: stencil1[cnt].loc = DMSTAG_DOWN;
100: stencil1[cnt].c = d;
101: ++cnt;
102: }
103: for (d = 0; d < dof1; ++d) {
104: stencil1[cnt].loc = DMSTAG_LEFT;
105: stencil1[cnt].c = d;
106: ++cnt;
107: }
108: /* element */
109: for (d = 0; d < dof2; ++d) {
110: stencil2[d].loc = DMSTAG_ELEMENT;
111: stencil2[d].c = d;
112: }
113: } else if (dim == 3) {
114: /* vertex, edge(down,left), face(down,left,back), element */
115: for (d = 0; d < dof0; ++d) {
116: stencil0[d].loc = DMSTAG_BACK_DOWN_LEFT;
117: stencil0[d].c = d;
118: }
119: /* edges */
120: cnt = 0;
121: for (d = 0; d < dof1; ++d) {
122: stencil1[cnt].loc = DMSTAG_BACK_DOWN;
123: stencil1[cnt].c = d;
124: ++cnt;
125: }
126: for (d = 0; d < dof1; ++d) {
127: stencil1[cnt].loc = DMSTAG_BACK_LEFT;
128: stencil1[cnt].c = d;
129: ++cnt;
130: }
131: for (d = 0; d < dof1; ++d) {
132: stencil1[cnt].loc = DMSTAG_DOWN_LEFT;
133: stencil1[cnt].c = d;
134: ++cnt;
135: }
136: /* faces */
137: cnt = 0;
138: for (d = 0; d < dof2; ++d) {
139: stencil2[cnt].loc = DMSTAG_BACK;
140: stencil2[cnt].c = d;
141: ++cnt;
142: }
143: for (d = 0; d < dof2; ++d) {
144: stencil2[cnt].loc = DMSTAG_DOWN;
145: stencil2[cnt].c = d;
146: ++cnt;
147: }
148: for (d = 0; d < dof2; ++d) {
149: stencil2[cnt].loc = DMSTAG_LEFT;
150: stencil2[cnt].c = d;
151: ++cnt;
152: }
153: /* elements */
154: for (d = 0; d < dof3; ++d) {
155: stencil3[d].loc = DMSTAG_ELEMENT;
156: stencil3[d].c = d;
157: }
158: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
160: cnt = 0;
161: if (dof0 != 0) {
162: PetscCall(DMStagCreateISFromStencils(dm, f0 * dof0, stencil0, &(*islist)[cnt]));
163: ++cnt;
164: }
165: if (dof1 != 0) {
166: PetscCall(DMStagCreateISFromStencils(dm, f1 * dof1, stencil1, &(*islist)[cnt]));
167: ++cnt;
168: }
169: if (dim >= 2 && dof2 != 0) {
170: PetscCall(DMStagCreateISFromStencils(dm, f2 * dof2, stencil2, &(*islist)[cnt]));
171: ++cnt;
172: }
173: if (dim >= 3 && dof3 != 0) {
174: PetscCall(DMStagCreateISFromStencils(dm, f3 * dof3, stencil3, &(*islist)[cnt]));
175: ++cnt;
176: }
177: }
179: if (namelist) {
180: PetscCall(PetscMalloc1(n_fields, namelist));
181: cnt = 0;
182: if (dim == 1) {
183: if (dof0 != 0) {
184: PetscCall(PetscStrallocpy("vertex", &(*namelist)[cnt]));
185: ++cnt;
186: }
187: if (dof1 != 0) {
188: PetscCall(PetscStrallocpy("element", &(*namelist)[cnt]));
189: ++cnt;
190: }
191: } else if (dim == 2) {
192: if (dof0 != 0) {
193: PetscCall(PetscStrallocpy("vertex", &(*namelist)[cnt]));
194: ++cnt;
195: }
196: if (dof1 != 0) {
197: PetscCall(PetscStrallocpy("face", &(*namelist)[cnt]));
198: ++cnt;
199: }
200: if (dof2 != 0) {
201: PetscCall(PetscStrallocpy("element", &(*namelist)[cnt]));
202: ++cnt;
203: }
204: } else if (dim == 3) {
205: if (dof0 != 0) {
206: PetscCall(PetscStrallocpy("vertex", &(*namelist)[cnt]));
207: ++cnt;
208: }
209: if (dof1 != 0) {
210: PetscCall(PetscStrallocpy("edge", &(*namelist)[cnt]));
211: ++cnt;
212: }
213: if (dof2 != 0) {
214: PetscCall(PetscStrallocpy("face", &(*namelist)[cnt]));
215: ++cnt;
216: }
217: if (dof3 != 0) {
218: PetscCall(PetscStrallocpy("element", &(*namelist)[cnt]));
219: ++cnt;
220: }
221: }
222: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
223: if (dmlist) {
224: PetscCall(PetscMalloc1(n_fields, dmlist));
225: cnt = 0;
226: if (dof0 != 0) {
227: PetscCall(DMStagCreateCompatibleDMStag(dm, dof0, 0, 0, 0, &(*dmlist)[cnt]));
228: ++cnt;
229: }
230: if (dof1 != 0) {
231: PetscCall(DMStagCreateCompatibleDMStag(dm, 0, dof1, 0, 0, &(*dmlist)[cnt]));
232: ++cnt;
233: }
234: if (dim >= 2 && dof2 != 0) {
235: PetscCall(DMStagCreateCompatibleDMStag(dm, 0, 0, dof2, 0, &(*dmlist)[cnt]));
236: ++cnt;
237: }
238: if (dim >= 3 && dof3 != 0) {
239: PetscCall(DMStagCreateCompatibleDMStag(dm, 0, 0, 0, dof3, &(*dmlist)[cnt]));
240: ++cnt;
241: }
242: }
243: PetscCall(PetscFree(stencil0));
244: PetscCall(PetscFree(stencil1));
245: if (dim >= 2) PetscCall(PetscFree(stencil2));
246: if (dim >= 3) PetscCall(PetscFree(stencil3));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: static PetscErrorCode DMClone_Stag(DM dm, DM *newdm)
251: {
252: PetscFunctionBegin;
253: /* Destroy the DM created by generic logic in DMClone() */
254: if (*newdm) PetscCall(DMDestroy(newdm));
255: PetscCall(DMStagDuplicateWithoutSetup(dm, PetscObjectComm((PetscObject)dm), newdm));
256: PetscCall(DMSetUp(*newdm));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: static PetscErrorCode DMCoarsen_Stag(DM dm, MPI_Comm comm, DM *dmc)
261: {
262: const DM_Stag *const stag = (DM_Stag *)dm->data;
263: PetscInt d, dim;
265: PetscFunctionBegin;
266: PetscCall(DMStagDuplicateWithoutSetup(dm, comm, dmc));
267: PetscCall(DMSetOptionsPrefix(*dmc, ((PetscObject)dm)->prefix));
268: PetscCall(DMGetDimension(dm, &dim));
269: for (d = 0; d < dim; ++d) PetscCheck(stag->N[d] % 2 == 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coarsening not supported except for even numbers of elements in each dimension ");
270: PetscCall(DMStagSetGlobalSizes(*dmc, stag->N[0] / 2, stag->N[1] / 2, stag->N[2] / 2));
271: {
272: PetscInt *l[DMSTAG_MAX_DIM];
273: for (d = 0; d < dim; ++d) {
274: PetscInt i;
275: PetscCall(PetscMalloc1(stag->nRanks[d], &l[d]));
276: for (i = 0; i < stag->nRanks[d]; ++i) {
277: PetscCheck(stag->l[d][i] % 2 == 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coarsening not supported except for an even number of elements in each direction on each rank");
278: l[d][i] = stag->l[d][i] / 2; /* Just halve everything */
279: }
280: }
281: PetscCall(DMStagSetOwnershipRanges(*dmc, l[0], l[1], l[2]));
282: for (d = 0; d < dim; ++d) PetscCall(PetscFree(l[d]));
283: }
284: PetscCall(DMSetUp(*dmc));
286: if (dm->coordinates[0].dm) { /* Note that with product coordinates, dm->coordinates = NULL, so we check the DM */
287: DM coordinate_dm, coordinate_dmc;
288: PetscBool isstag, isprod;
290: PetscCall(DMGetCoordinateDM(dm, &coordinate_dm));
291: PetscCall(PetscObjectTypeCompare((PetscObject)coordinate_dm, DMSTAG, &isstag));
292: PetscCall(PetscObjectTypeCompare((PetscObject)coordinate_dm, DMPRODUCT, &isprod));
293: if (isstag) {
294: PetscCall(DMStagSetUniformCoordinatesExplicit(*dmc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); /* Coordinates will be overwritten */
295: PetscCall(DMGetCoordinateDM(*dmc, &coordinate_dmc));
296: PetscCall(DMStagRestrictSimple(coordinate_dm, dm->coordinates[0].x, coordinate_dmc, (*dmc)->coordinates[0].x));
297: } else if (isprod) {
298: PetscCall(DMStagSetUniformCoordinatesProduct(*dmc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); /* Coordinates will be overwritten */
299: PetscCall(DMGetCoordinateDM(*dmc, &coordinate_dmc));
300: for (d = 0; d < dim; ++d) {
301: DM subdm_coarse, subdm_coord_coarse, subdm_fine, subdm_coord_fine;
303: PetscCall(DMProductGetDM(coordinate_dm, d, &subdm_fine));
304: PetscCall(DMGetCoordinateDM(subdm_fine, &subdm_coord_fine));
305: PetscCall(DMProductGetDM(coordinate_dmc, d, &subdm_coarse));
306: PetscCall(DMGetCoordinateDM(subdm_coarse, &subdm_coord_coarse));
307: PetscCall(DMStagRestrictSimple(subdm_coord_fine, subdm_fine->coordinates[0].xl, subdm_coord_coarse, subdm_coarse->coordinates[0].xl));
308: }
309: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown coordinate DM type");
310: }
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: static PetscErrorCode DMRefine_Stag(DM dm, MPI_Comm comm, DM *dmc)
315: {
316: const DM_Stag *const stag = (DM_Stag *)dm->data;
318: PetscFunctionBegin;
319: PetscCall(DMStagDuplicateWithoutSetup(dm, comm, dmc));
320: PetscCall(DMSetOptionsPrefix(*dmc, ((PetscObject)dm)->prefix));
321: PetscCall(DMStagSetGlobalSizes(*dmc, stag->N[0] * 2, stag->N[1] * 2, stag->N[2] * 2));
322: {
323: PetscInt dim, d;
324: PetscInt *l[DMSTAG_MAX_DIM];
325: PetscCall(DMGetDimension(dm, &dim));
326: for (d = 0; d < dim; ++d) {
327: PetscInt i;
328: PetscCall(PetscMalloc1(stag->nRanks[d], &l[d]));
329: for (i = 0; i < stag->nRanks[d]; ++i) { l[d][i] = stag->l[d][i] * 2; /* Just double everything */ }
330: }
331: PetscCall(DMStagSetOwnershipRanges(*dmc, l[0], l[1], l[2]));
332: for (d = 0; d < dim; ++d) PetscCall(PetscFree(l[d]));
333: }
334: PetscCall(DMSetUp(*dmc));
335: /* Note: For now, we do not refine coordinates */
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: static PetscErrorCode DMDestroy_Stag(DM dm)
340: {
341: DM_Stag *stag;
342: PetscInt i;
344: PetscFunctionBegin;
345: stag = (DM_Stag *)dm->data;
346: for (i = 0; i < DMSTAG_MAX_DIM; ++i) PetscCall(PetscFree(stag->l[i]));
347: PetscCall(VecScatterDestroy(&stag->gtol));
348: PetscCall(VecScatterDestroy(&stag->ltog_injective));
349: PetscCall(PetscFree(stag->neighbors));
350: PetscCall(PetscFree(stag->locationOffsets));
351: PetscCall(PetscFree(stag->coordinateDMType));
352: PetscCall(PetscFree(stag));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode DMCreateGlobalVector_Stag(DM dm, Vec *vec)
357: {
358: DM_Stag *const stag = (DM_Stag *)dm->data;
360: PetscFunctionBegin;
361: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
362: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), vec));
363: PetscCall(VecSetSizes(*vec, stag->entries, PETSC_DETERMINE));
364: PetscCall(VecSetType(*vec, dm->vectype));
365: PetscCall(VecSetDM(*vec, dm));
366: /* Could set some ops, as DMDA does */
367: PetscCall(VecSetLocalToGlobalMapping(*vec, dm->ltogmap));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode DMCreateLocalVector_Stag(DM dm, Vec *vec)
372: {
373: DM_Stag *const stag = (DM_Stag *)dm->data;
375: PetscFunctionBegin;
376: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
377: PetscCall(VecCreate(PETSC_COMM_SELF, vec));
378: PetscCall(VecSetSizes(*vec, stag->entriesGhost, PETSC_DETERMINE));
379: PetscCall(VecSetType(*vec, dm->vectype));
380: PetscCall(VecSetBlockSize(*vec, stag->entriesPerElement));
381: PetscCall(VecSetDM(*vec, dm));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /* Helper function to check for the limited situations for which interpolation
386: and restriction functions are implemented */
387: static PetscErrorCode CheckTransferOperatorRequirements_Private(DM dmc, DM dmf)
388: {
389: PetscInt dim, stencilWidthc, stencilWidthf, nf[DMSTAG_MAX_DIM], nc[DMSTAG_MAX_DIM], doff[DMSTAG_MAX_STRATA], dofc[DMSTAG_MAX_STRATA];
391: PetscFunctionBegin;
392: PetscCall(DMGetDimension(dmc, &dim));
393: PetscCall(DMStagGetStencilWidth(dmc, &stencilWidthc));
394: PetscCheck(stencilWidthc >= 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "DMCreateRestriction not implemented for coarse grid stencil width < 1");
395: PetscCall(DMStagGetStencilWidth(dmf, &stencilWidthf));
396: PetscCheck(stencilWidthf >= 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "DMCreateRestriction not implemented for fine grid stencil width < 1");
397: PetscCall(DMStagGetLocalSizes(dmf, &nf[0], &nf[1], &nf[2]));
398: PetscCall(DMStagGetLocalSizes(dmc, &nc[0], &nc[1], &nc[2]));
399: for (PetscInt d = 0; d < dim; ++d)
400: PetscCheck(nf[d] == 2 * nc[d], PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No support for fine to coarse ratio other than 2 (it is %" PetscInt_FMT " to %" PetscInt_FMT " in dimension %" PetscInt_FMT ")", nf[d], nc[d], d);
401: PetscCall(DMStagGetDOF(dmc, &dofc[0], &dofc[1], &dofc[2], &dofc[3]));
402: PetscCall(DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]));
403: for (PetscInt d = 0; d < dim + 1; ++d)
404: PetscCheck(dofc[d] == doff[d], PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No support for different numbers of dof per stratum between coarse and fine DMStag objects: dof%" PetscInt_FMT " is %" PetscInt_FMT " (fine) but %" PetscInt_FMT "(coarse))", d, doff[d], dofc[d]);
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /* Since the interpolation uses MATMAIJ for dof > 0 we convert requests for non-MATAIJ baseded matrices to MATAIJ.
409: This is a bit of a hack; the reason for it is partially because -dm_mat_type defines the
410: matrix type for both the operator matrices and the interpolation matrices so that users
411: can select matrix types of base MATAIJ for accelerators
413: Note: The ConvertToAIJ() code below *has been copied from dainterp.c*! ConvertToAIJ() should perhaps be placed somewhere
414: in mat/utils to avoid code duplication, but then the DMStag and DMDA code would need to include the private Mat headers.
415: Since it is only used in two places, I have simply duplicated the code to avoid the need to exposure the private
416: Mat routines in parts of DM. If we find a need for ConvertToAIJ() elsewhere, then we should consolidate it to one
417: place in mat/utils.
418: */
419: static PetscErrorCode ConvertToAIJ(MatType intype, MatType *outtype)
420: {
421: PetscInt i;
422: char const *types[3] = {MATAIJ, MATSEQAIJ, MATMPIAIJ};
423: PetscBool flg;
425: PetscFunctionBegin;
426: *outtype = MATAIJ;
427: for (i = 0; i < 3; i++) {
428: PetscCall(PetscStrbeginswith(intype, types[i], &flg));
429: if (flg) {
430: *outtype = intype;
431: break;
432: }
433: }
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: static PetscErrorCode DMCreateInterpolation_Stag(DM dmc, DM dmf, Mat *A, Vec *vec)
438: {
439: PetscInt dim, entriesf, entriesc, doff[DMSTAG_MAX_STRATA];
440: ISLocalToGlobalMapping ltogmf, ltogmc;
441: MatType mattype;
443: PetscFunctionBegin;
444: PetscCall(CheckTransferOperatorRequirements_Private(dmc, dmf));
446: PetscCall(DMStagGetEntries(dmf, &entriesf));
447: PetscCall(DMStagGetEntries(dmc, &entriesc));
448: PetscCall(DMGetLocalToGlobalMapping(dmf, <ogmf));
449: PetscCall(DMGetLocalToGlobalMapping(dmc, <ogmc));
451: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), A));
452: PetscCall(MatSetSizes(*A, entriesf, entriesc, PETSC_DECIDE, PETSC_DECIDE));
453: PetscCall(ConvertToAIJ(dmc->mattype, &mattype));
454: PetscCall(MatSetType(*A, mattype));
455: PetscCall(MatSetLocalToGlobalMapping(*A, ltogmf, ltogmc));
457: PetscCall(DMGetDimension(dmc, &dim));
458: PetscCall(DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]));
459: if (dim == 1) {
460: PetscCall(DMStagPopulateInterpolation1d_a_b_Private(dmc, dmf, *A));
461: } else if (dim == 2) {
462: if (doff[0] == 0) {
463: PetscCall(DMStagPopulateInterpolation2d_0_a_b_Private(dmc, dmf, *A));
464: } else
465: SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No default interpolation available between 2d DMStag objects with %" PetscInt_FMT " dof/vertex, %" PetscInt_FMT " dof/face and %" PetscInt_FMT " dof/element", doff[0], doff[1], doff[2]);
466: } else if (dim == 3) {
467: if (doff[0] == 0 && doff[1] == 0) {
468: PetscCall(DMStagPopulateInterpolation3d_0_0_a_b_Private(dmc, dmf, *A));
469: } else
470: SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No default interpolation available between 3d DMStag objects with %" PetscInt_FMT " dof/vertex, %" PetscInt_FMT " dof/edge, %" PetscInt_FMT " dof/face and %" PetscInt_FMT " dof/element", doff[0], doff[1], doff[2], doff[3]);
471: } else SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
472: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
473: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
475: if (vec) *vec = NULL;
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: static PetscErrorCode DMCreateRestriction_Stag(DM dmc, DM dmf, Mat *A)
480: {
481: PetscInt dim, entriesf, entriesc, doff[DMSTAG_MAX_STRATA];
482: ISLocalToGlobalMapping ltogmf, ltogmc;
483: MatType mattype;
485: PetscFunctionBegin;
486: PetscCall(CheckTransferOperatorRequirements_Private(dmc, dmf));
488: PetscCall(DMStagGetEntries(dmf, &entriesf));
489: PetscCall(DMStagGetEntries(dmc, &entriesc));
490: PetscCall(DMGetLocalToGlobalMapping(dmf, <ogmf));
491: PetscCall(DMGetLocalToGlobalMapping(dmc, <ogmc));
493: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), A));
494: PetscCall(MatSetSizes(*A, entriesc, entriesf, PETSC_DECIDE, PETSC_DECIDE)); /* Note transpose wrt interpolation */
495: PetscCall(ConvertToAIJ(dmc->mattype, &mattype));
496: PetscCall(MatSetType(*A, mattype));
497: PetscCall(MatSetLocalToGlobalMapping(*A, ltogmc, ltogmf)); /* Note transpose wrt interpolation */
499: PetscCall(DMGetDimension(dmc, &dim));
500: PetscCall(DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]));
501: if (dim == 1) {
502: PetscCall(DMStagPopulateRestriction1d_a_b_Private(dmc, dmf, *A));
503: } else if (dim == 2) {
504: if (doff[0] == 0) {
505: PetscCall(DMStagPopulateRestriction2d_0_a_b_Private(dmc, dmf, *A));
506: } else
507: SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No default restriction available between 2d DMStag objects with %" PetscInt_FMT " dof/vertex, %" PetscInt_FMT " dof/face and %" PetscInt_FMT " dof/element", doff[0], doff[1], doff[2]);
508: } else if (dim == 3) {
509: if (doff[0] == 0 && doff[0] == 0) {
510: PetscCall(DMStagPopulateRestriction3d_0_0_a_b_Private(dmc, dmf, *A));
511: } else
512: SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "No default restriction available between 3d DMStag objects with %" PetscInt_FMT " dof/vertex, %" PetscInt_FMT " dof/edge, %" PetscInt_FMT " dof/face and %" PetscInt_FMT " dof/element", doff[0], doff[1], doff[2], doff[3]);
513: } else SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
515: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
516: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: static PetscErrorCode DMCreateMatrix_Stag(DM dm, Mat *mat)
521: {
522: MatType mat_type;
523: PetscBool is_shell, is_aij;
524: PetscInt dim, entries;
525: ISLocalToGlobalMapping ltogmap;
527: PetscFunctionBegin;
528: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
529: PetscCall(DMGetDimension(dm, &dim));
530: PetscCall(DMGetMatType(dm, &mat_type));
531: PetscCall(DMStagGetEntries(dm, &entries));
532: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), mat));
533: PetscCall(MatSetSizes(*mat, entries, entries, PETSC_DETERMINE, PETSC_DETERMINE));
534: PetscCall(MatSetType(*mat, mat_type));
535: PetscCall(MatSetUp(*mat));
536: PetscCall(DMGetLocalToGlobalMapping(dm, <ogmap));
537: PetscCall(MatSetLocalToGlobalMapping(*mat, ltogmap, ltogmap));
538: PetscCall(MatSetDM(*mat, dm));
540: /* Compare to similar and perhaps superior logic in DMCreateMatrix_DA, which creates
541: the matrix first and then performs this logic by checking for preallocation functions */
542: PetscCall(PetscObjectBaseTypeCompare((PetscObject)*mat, MATAIJ, &is_aij));
543: if (!is_aij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)*mat, MATSEQAIJ, &is_aij));
544: if (!is_aij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)*mat, MATMPIAIJ, &is_aij));
545: PetscCall(PetscStrcmp(mat_type, MATSHELL, &is_shell));
546: if (is_aij) {
547: Mat preallocator;
548: PetscInt m, n;
549: const PetscBool fill_with_zeros = PETSC_FALSE;
551: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), &preallocator));
552: PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
553: PetscCall(MatGetLocalSize(*mat, &m, &n));
554: PetscCall(MatSetSizes(preallocator, m, n, PETSC_DECIDE, PETSC_DECIDE));
555: PetscCall(MatSetLocalToGlobalMapping(preallocator, ltogmap, ltogmap));
556: PetscCall(MatSetUp(preallocator));
557: switch (dim) {
558: case 1:
559: PetscCall(DMCreateMatrix_Stag_1D_AIJ_Assemble(dm, preallocator));
560: break;
561: case 2:
562: PetscCall(DMCreateMatrix_Stag_2D_AIJ_Assemble(dm, preallocator));
563: break;
564: case 3:
565: PetscCall(DMCreateMatrix_Stag_3D_AIJ_Assemble(dm, preallocator));
566: break;
567: default:
568: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
569: }
570: PetscCall(MatPreallocatorPreallocate(preallocator, fill_with_zeros, *mat));
571: PetscCall(MatDestroy(&preallocator));
573: if (!dm->prealloc_only) {
574: /* Bind to CPU before assembly, to prevent unnecessary copies of zero entries from CPU to GPU */
575: PetscCall(MatBindToCPU(*mat, PETSC_TRUE));
576: switch (dim) {
577: case 1:
578: PetscCall(DMCreateMatrix_Stag_1D_AIJ_Assemble(dm, *mat));
579: break;
580: case 2:
581: PetscCall(DMCreateMatrix_Stag_2D_AIJ_Assemble(dm, *mat));
582: break;
583: case 3:
584: PetscCall(DMCreateMatrix_Stag_3D_AIJ_Assemble(dm, *mat));
585: break;
586: default:
587: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
588: }
589: PetscCall(MatBindToCPU(*mat, PETSC_FALSE));
590: }
591: } else if (is_shell) {
592: /* nothing more to do */
593: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for Mattype %s", mat_type);
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: static PetscErrorCode DMGetCompatibility_Stag(DM dm, DM dm2, PetscBool *compatible, PetscBool *set)
598: {
599: const DM_Stag *const stag = (DM_Stag *)dm->data;
600: const DM_Stag *const stag2 = (DM_Stag *)dm2->data;
601: PetscInt dim, dim2, i;
602: MPI_Comm comm;
603: PetscMPIInt sameComm;
604: DMType type2;
605: PetscBool sameType;
607: PetscFunctionBegin;
608: PetscCall(DMGetType(dm2, &type2));
609: PetscCall(PetscStrcmp(DMSTAG, type2, &sameType));
610: if (!sameType) {
611: PetscCall(PetscInfo((PetscObject)dm, "DMStag compatibility check not implemented with DM of type %s\n", type2));
612: *set = PETSC_FALSE;
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
617: PetscCallMPI(MPI_Comm_compare(comm, PetscObjectComm((PetscObject)dm2), &sameComm));
618: if (sameComm != MPI_IDENT) {
619: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different communicators: %" PETSC_INTPTR_T_FMT " != %" PETSC_INTPTR_T_FMT "\n", (PETSC_INTPTR_T)comm, (PETSC_INTPTR_T)PetscObjectComm((PetscObject)dm2)));
620: *set = PETSC_FALSE;
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
623: PetscCall(DMGetDimension(dm, &dim));
624: PetscCall(DMGetDimension(dm2, &dim2));
625: if (dim != dim2) {
626: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different dimensions\n"));
627: *set = PETSC_TRUE;
628: *compatible = PETSC_FALSE;
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
631: for (i = 0; i < dim; ++i) {
632: if (stag->N[i] != stag2->N[i]) {
633: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different global numbers of elements in dimension %" PetscInt_FMT ": %" PetscInt_FMT " != %" PetscInt_FMT "\n", i, stag->n[i], stag2->n[i]));
634: *set = PETSC_TRUE;
635: *compatible = PETSC_FALSE;
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
638: if (stag->n[i] != stag2->n[i]) {
639: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different local numbers of elements in dimension %" PetscInt_FMT ": %" PetscInt_FMT " != %" PetscInt_FMT "\n", i, stag->n[i], stag2->n[i]));
640: *set = PETSC_TRUE;
641: *compatible = PETSC_FALSE;
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
644: if (stag->boundaryType[i] != stag2->boundaryType[i]) {
645: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different boundary types in dimension %" PetscInt_FMT ": %s != %s\n", i, DMBoundaryTypes[stag->boundaryType[i]], DMBoundaryTypes[stag2->boundaryType[i]]));
646: *set = PETSC_TRUE;
647: *compatible = PETSC_FALSE;
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
650: }
651: /* Note: we include stencil type and width in the notion of compatibility, as this affects
652: the "atlas" (local subdomains). This might be irritating in legitimate cases
653: of wanting to transfer between two other-wise compatible DMs with different
654: stencil characteristics. */
655: if (stag->stencilType != stag2->stencilType) {
656: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different ghost stencil types: %s != %s\n", DMStagStencilTypes[stag->stencilType], DMStagStencilTypes[stag2->stencilType]));
657: *set = PETSC_TRUE;
658: *compatible = PETSC_FALSE;
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
661: if (stag->stencilWidth != stag2->stencilWidth) {
662: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different ghost stencil widths: %" PetscInt_FMT " != %" PetscInt_FMT "\n", stag->stencilWidth, stag->stencilWidth));
663: *set = PETSC_TRUE;
664: *compatible = PETSC_FALSE;
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
667: *set = PETSC_TRUE;
668: *compatible = PETSC_TRUE;
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: static PetscErrorCode DMHasCreateInjection_Stag(DM dm, PetscBool *flg)
673: {
674: PetscFunctionBegin;
676: PetscAssertPointer(flg, 2);
677: *flg = PETSC_FALSE;
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: /*
682: Note there are several orderings in play here.
683: In all cases, non-element dof are associated with the element that they are below/left/behind, and the order in 2D proceeds vertex/bottom edge/left edge/element (with all dof on each together).
684: Also in all cases, only subdomains which are the last in their dimension have partial elements.
686: 1) "Natural" Ordering (not used). Number adding each full or partial (on the right or top) element, starting at the bottom left (i=0,j=0) and proceeding across the entire domain, row by row to get a global numbering.
687: 2) Global ("PETSc") ordering. The same as natural, but restricted to each domain. So, traverse all elements (again starting at the bottom left and going row-by-row) on rank 0, then continue numbering with rank 1, and so on.
688: 3) Local ordering. Including ghost elements (both interior and on the right/top/front to complete partial elements), use the same convention to create a local numbering.
689: */
691: static PetscErrorCode DMLocalToGlobalBegin_Stag(DM dm, Vec l, InsertMode mode, Vec g)
692: {
693: DM_Stag *const stag = (DM_Stag *)dm->data;
695: PetscFunctionBegin;
696: if (mode == ADD_VALUES) {
697: PetscCall(VecScatterBegin(stag->gtol, l, g, mode, SCATTER_REVERSE));
698: } else if (mode == INSERT_VALUES) {
699: if (stag->ltog_injective) {
700: PetscCall(VecScatterBegin(stag->ltog_injective, l, g, mode, SCATTER_FORWARD));
701: } else {
702: PetscCall(VecScatterBegin(stag->gtol, l, g, mode, SCATTER_REVERSE_LOCAL));
703: }
704: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported InsertMode");
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: static PetscErrorCode DMLocalToGlobalEnd_Stag(DM dm, Vec l, InsertMode mode, Vec g)
709: {
710: DM_Stag *const stag = (DM_Stag *)dm->data;
712: PetscFunctionBegin;
713: if (mode == ADD_VALUES) {
714: PetscCall(VecScatterEnd(stag->gtol, l, g, mode, SCATTER_REVERSE));
715: } else if (mode == INSERT_VALUES) {
716: if (stag->ltog_injective) {
717: PetscCall(VecScatterEnd(stag->ltog_injective, l, g, mode, SCATTER_FORWARD));
718: } else {
719: PetscCall(VecScatterEnd(stag->gtol, l, g, mode, SCATTER_REVERSE_LOCAL));
720: }
721: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported InsertMode");
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: static PetscErrorCode DMGlobalToLocalBegin_Stag(DM dm, Vec g, InsertMode mode, Vec l)
726: {
727: DM_Stag *const stag = (DM_Stag *)dm->data;
729: PetscFunctionBegin;
730: PetscCall(VecScatterBegin(stag->gtol, g, l, mode, SCATTER_FORWARD));
731: PetscFunctionReturn(PETSC_SUCCESS);
732: }
734: static PetscErrorCode DMGlobalToLocalEnd_Stag(DM dm, Vec g, InsertMode mode, Vec l)
735: {
736: DM_Stag *const stag = (DM_Stag *)dm->data;
738: PetscFunctionBegin;
739: PetscCall(VecScatterEnd(stag->gtol, g, l, mode, SCATTER_FORWARD));
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*
744: If a stratum is active (non-zero dof), make it active in the coordinate DM.
745: */
746: static PetscErrorCode DMCreateCoordinateDM_Stag(DM dm, DM *dmc)
747: {
748: DM_Stag *const stag = (DM_Stag *)dm->data;
749: PetscInt dim;
750: PetscBool isstag, isproduct;
751: const char *prefix;
753: PetscFunctionBegin;
755: PetscCheck(stag->coordinateDMType, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Before creating a coordinate DM, a type must be specified with DMStagSetCoordinateDMType()");
757: PetscCall(DMGetDimension(dm, &dim));
758: PetscCall(PetscStrcmp(stag->coordinateDMType, DMSTAG, &isstag));
759: PetscCall(PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &isproduct));
760: if (isstag) {
761: PetscCall(DMStagCreateCompatibleDMStag(dm, stag->dof[0] > 0 ? dim : 0, stag->dof[1] > 0 ? dim : 0, stag->dof[2] > 0 ? dim : 0, stag->dof[3] > 0 ? dim : 0, dmc));
762: } else if (isproduct) {
763: PetscCall(DMCreate(PETSC_COMM_WORLD, dmc));
764: PetscCall(DMSetType(*dmc, DMPRODUCT));
765: PetscCall(DMSetDimension(*dmc, dim));
766: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate DM type %s", stag->coordinateDMType);
767: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
768: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dmc, prefix));
769: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*dmc, "cdm_"));
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: static PetscErrorCode DMGetNeighbors_Stag(DM dm, PetscInt *nRanks, const PetscMPIInt *ranks[])
774: {
775: DM_Stag *const stag = (DM_Stag *)dm->data;
776: PetscInt dim;
778: PetscFunctionBegin;
779: PetscCall(DMGetDimension(dm, &dim));
780: switch (dim) {
781: case 1:
782: *nRanks = 3;
783: break;
784: case 2:
785: *nRanks = 9;
786: break;
787: case 3:
788: *nRanks = 27;
789: break;
790: default:
791: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Get neighbors not implemented for dim = %" PetscInt_FMT, dim);
792: }
793: *ranks = stag->neighbors;
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: static PetscErrorCode DMView_Stag(DM dm, PetscViewer viewer)
798: {
799: DM_Stag *const stag = (DM_Stag *)dm->data;
800: PetscBool isascii, viewAllRanks;
801: PetscMPIInt rank, size;
802: PetscInt dim, maxRanksToView, i;
804: PetscFunctionBegin;
805: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
806: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
807: PetscCall(DMGetDimension(dm, &dim));
808: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
809: if (isascii) {
810: PetscCall(PetscViewerASCIIPrintf(viewer, "Dimension: %" PetscInt_FMT "\n", dim));
811: switch (dim) {
812: case 1:
813: PetscCall(PetscViewerASCIIPrintf(viewer, "Global size: %" PetscInt_FMT "\n", stag->N[0]));
814: break;
815: case 2:
816: PetscCall(PetscViewerASCIIPrintf(viewer, "Global sizes: %" PetscInt_FMT " x %" PetscInt_FMT "\n", stag->N[0], stag->N[1]));
817: PetscCall(PetscViewerASCIIPrintf(viewer, "Parallel decomposition: %" PetscInt_FMT " x %" PetscInt_FMT " ranks\n", stag->nRanks[0], stag->nRanks[1]));
818: break;
819: case 3:
820: PetscCall(PetscViewerASCIIPrintf(viewer, "Global sizes: %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "\n", stag->N[0], stag->N[1], stag->N[2]));
821: PetscCall(PetscViewerASCIIPrintf(viewer, "Parallel decomposition: %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT " ranks\n", stag->nRanks[0], stag->nRanks[1], stag->nRanks[2]));
822: break;
823: default:
824: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim==%" PetscInt_FMT, dim);
825: }
826: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary ghosting:"));
827: for (i = 0; i < dim; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, " %s", DMBoundaryTypes[stag->boundaryType[i]]));
828: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
829: PetscCall(PetscViewerASCIIPrintf(viewer, "Elementwise ghost stencil: %s", DMStagStencilTypes[stag->stencilType]));
830: if (stag->stencilType != DMSTAG_STENCIL_NONE) {
831: PetscCall(PetscViewerASCIIPrintf(viewer, ", width %" PetscInt_FMT "\n", stag->stencilWidth));
832: } else {
833: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
834: }
835: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per vertex (0D)\n", stag->dof[0]));
836: if (dim == 3) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per edge (1D)\n", stag->dof[1]));
837: if (dim > 1) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per face (%" PetscInt_FMT "D)\n", stag->dof[dim - 1], dim - 1));
838: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per element (%" PetscInt_FMT "D)\n", stag->dof[dim], dim));
839: if (dm->coordinates[0].dm) PetscCall(PetscViewerASCIIPrintf(viewer, "Has coordinate DM\n"));
840: maxRanksToView = 16;
841: viewAllRanks = (PetscBool)(size <= maxRanksToView);
842: if (viewAllRanks) {
843: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
844: switch (dim) {
845: case 1:
846: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local elements : %" PetscInt_FMT " (%" PetscInt_FMT " with ghosts)\n", rank, stag->n[0], stag->nGhost[0]));
847: break;
848: case 2:
849: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Rank coordinates (%d,%d)\n", rank, stag->rank[0], stag->rank[1]));
850: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local elements : %" PetscInt_FMT " x %" PetscInt_FMT " (%" PetscInt_FMT " x %" PetscInt_FMT " with ghosts)\n", rank, stag->n[0], stag->n[1], stag->nGhost[0], stag->nGhost[1]));
851: break;
852: case 3:
853: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Rank coordinates (%d,%d,%d)\n", rank, stag->rank[0], stag->rank[1], stag->rank[2]));
854: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local elements : %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT " (%" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT " with ghosts)\n", rank, stag->n[0], stag->n[1],
855: stag->n[2], stag->nGhost[0], stag->nGhost[1], stag->nGhost[2]));
856: break;
857: default:
858: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim==%" PetscInt_FMT, dim);
859: }
860: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local native entries: %" PetscInt_FMT "\n", rank, stag->entries));
861: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local entries total : %" PetscInt_FMT "\n", rank, stag->entriesGhost));
862: PetscCall(PetscViewerFlush(viewer));
863: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
864: } else {
865: PetscCall(PetscViewerASCIIPrintf(viewer, "(Per-rank information omitted since >%" PetscInt_FMT " ranks used)\n", maxRanksToView));
866: }
867: }
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
871: static PetscErrorCode DMSetFromOptions_Stag(DM dm, PetscOptionItems *PetscOptionsObject)
872: {
873: DM_Stag *const stag = (DM_Stag *)dm->data;
874: PetscInt dim;
876: PetscFunctionBegin;
877: PetscCall(DMGetDimension(dm, &dim));
878: PetscOptionsHeadBegin(PetscOptionsObject, "DMStag Options");
879: PetscCall(PetscOptionsInt("-stag_grid_x", "Number of grid points in x direction", "DMStagSetGlobalSizes", stag->N[0], &stag->N[0], NULL));
880: if (dim > 1) PetscCall(PetscOptionsInt("-stag_grid_y", "Number of grid points in y direction", "DMStagSetGlobalSizes", stag->N[1], &stag->N[1], NULL));
881: if (dim > 2) PetscCall(PetscOptionsInt("-stag_grid_z", "Number of grid points in z direction", "DMStagSetGlobalSizes", stag->N[2], &stag->N[2], NULL));
882: PetscCall(PetscOptionsInt("-stag_ranks_x", "Number of ranks in x direction", "DMStagSetNumRanks", stag->nRanks[0], &stag->nRanks[0], NULL));
883: if (dim > 1) PetscCall(PetscOptionsInt("-stag_ranks_y", "Number of ranks in y direction", "DMStagSetNumRanks", stag->nRanks[1], &stag->nRanks[1], NULL));
884: if (dim > 2) PetscCall(PetscOptionsInt("-stag_ranks_z", "Number of ranks in z direction", "DMStagSetNumRanks", stag->nRanks[2], &stag->nRanks[2], NULL));
885: PetscCall(PetscOptionsInt("-stag_stencil_width", "Elementwise stencil width", "DMStagSetStencilWidth", stag->stencilWidth, &stag->stencilWidth, NULL));
886: PetscCall(PetscOptionsEnum("-stag_stencil_type", "Elementwise stencil stype", "DMStagSetStencilType", DMStagStencilTypes, (PetscEnum)stag->stencilType, (PetscEnum *)&stag->stencilType, NULL));
887: PetscCall(PetscOptionsEnum("-stag_boundary_type_x", "Treatment of (physical) boundaries in x direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[0], (PetscEnum *)&stag->boundaryType[0], NULL));
888: PetscCall(PetscOptionsEnum("-stag_boundary_type_y", "Treatment of (physical) boundaries in y direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[1], (PetscEnum *)&stag->boundaryType[1], NULL));
889: PetscCall(PetscOptionsEnum("-stag_boundary_type_z", "Treatment of (physical) boundaries in z direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[2], (PetscEnum *)&stag->boundaryType[2], NULL));
890: PetscCall(PetscOptionsInt("-stag_dof_0", "Number of dof per 0-cell (vertex)", "DMStagSetDOF", stag->dof[0], &stag->dof[0], NULL));
891: PetscCall(PetscOptionsInt("-stag_dof_1", "Number of dof per 1-cell (element in 1D, face in 2D, edge in 3D)", "DMStagSetDOF", stag->dof[1], &stag->dof[1], NULL));
892: PetscCall(PetscOptionsInt("-stag_dof_2", "Number of dof per 2-cell (element in 2D, face in 3D)", "DMStagSetDOF", stag->dof[2], &stag->dof[2], NULL));
893: PetscCall(PetscOptionsInt("-stag_dof_3", "Number of dof per 3-cell (element in 3D)", "DMStagSetDOF", stag->dof[3], &stag->dof[3], NULL));
894: PetscOptionsHeadEnd();
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: /*MC
899: DMSTAG - `"stag"` - A `DM` object representing a "staggered grid" or a structured cell complex.
901: Level: beginner
903: Notes:
904: This implementation parallels the `DMDA` implementation in many ways, but allows degrees of freedom
905: to be associated with all "strata" in a logically-rectangular grid.
907: Each stratum can be characterized by the dimension of the entities ("points", to borrow the `DMPLEX`
908: terminology), from 0- to 3-dimensional.
910: In some cases this numbering is used directly, for example with `DMStagGetDOF()`.
911: To allow easier reading and to some extent more similar code between different-dimensional implementations
912: of the same problem, we associate canonical names for each type of point, for each dimension of DMStag.
914: * 1-dimensional `DMSTAG` objects have vertices (0D) and elements (1D).
915: * 2-dimensional `DMSTAG` objects have vertices (0D), faces (1D), and elements (2D).
916: * 3-dimensional `DMSTAG` objects have vertices (0D), edges (1D), faces (2D), and elements (3D).
918: This naming is reflected when viewing a `DMSTAG` object with `DMView()`, and in forming
919: convenient options prefixes when creating a decomposition with `DMCreateFieldDecomposition()`.
921: .seealso: [](ch_stag), `DM`, `DMPRODUCT`, `DMDA`, `DMPLEX`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMType`, `DMCreate()`,
922: `DMSetType()`, `DMStagVecSplitToDMDA()`
923: M*/
925: PETSC_EXTERN PetscErrorCode DMCreate_Stag(DM dm)
926: {
927: DM_Stag *stag;
928: PetscInt i, dim;
930: PetscFunctionBegin;
931: PetscAssertPointer(dm, 1);
932: PetscCall(PetscNew(&stag));
933: dm->data = stag;
935: stag->gtol = NULL;
936: stag->ltog_injective = NULL;
937: for (i = 0; i < DMSTAG_MAX_STRATA; ++i) stag->dof[i] = 0;
938: for (i = 0; i < DMSTAG_MAX_DIM; ++i) stag->l[i] = NULL;
939: stag->stencilType = DMSTAG_STENCIL_NONE;
940: stag->stencilWidth = 0;
941: for (i = 0; i < DMSTAG_MAX_DIM; ++i) stag->nRanks[i] = -1;
942: stag->coordinateDMType = NULL;
944: PetscCall(DMGetDimension(dm, &dim));
945: PetscCheck(dim == 1 || dim == 2 || dim == 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMSetDimension() must be called to set a dimension with value 1, 2, or 3");
947: PetscCall(PetscMemzero(dm->ops, sizeof(*(dm->ops))));
948: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Stag;
949: dm->ops->createglobalvector = DMCreateGlobalVector_Stag;
950: dm->ops->createlocalvector = DMCreateLocalVector_Stag;
951: dm->ops->creatematrix = DMCreateMatrix_Stag;
952: dm->ops->hascreateinjection = DMHasCreateInjection_Stag;
953: dm->ops->refine = DMRefine_Stag;
954: dm->ops->coarsen = DMCoarsen_Stag;
955: dm->ops->createinterpolation = DMCreateInterpolation_Stag;
956: dm->ops->createrestriction = DMCreateRestriction_Stag;
957: dm->ops->destroy = DMDestroy_Stag;
958: dm->ops->getneighbors = DMGetNeighbors_Stag;
959: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Stag;
960: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Stag;
961: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Stag;
962: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Stag;
963: dm->ops->setfromoptions = DMSetFromOptions_Stag;
964: switch (dim) {
965: case 1:
966: dm->ops->setup = DMSetUp_Stag_1d;
967: break;
968: case 2:
969: dm->ops->setup = DMSetUp_Stag_2d;
970: break;
971: case 3:
972: dm->ops->setup = DMSetUp_Stag_3d;
973: break;
974: default:
975: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
976: }
977: dm->ops->clone = DMClone_Stag;
978: dm->ops->view = DMView_Stag;
979: dm->ops->getcompatibility = DMGetCompatibility_Stag;
980: dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Stag;
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }