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