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 dim, i, d;
264: PetscInt *l[DMSTAG_MAX_DIM];
266: PetscFunctionBegin;
267: PetscCall(DMStagDuplicateWithoutSetup(dm, comm, dmc));
268: PetscCall(DMSetOptionsPrefix(*dmc, ((PetscObject)dm)->prefix));
269: PetscCall(DMGetDimension(dm, &dim));
270: 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");
271: PetscCall(DMStagSetGlobalSizes(*dmc, stag->N[0] / stag->refineFactor[0], stag->N[1] / stag->refineFactor[1], stag->N[2] / stag->refineFactor[2]));
272: for (d = 0; d < dim; ++d) {
273: PetscCall(PetscMalloc1(stag->nRanks[d], &l[d]));
274: for (i = 0; i < stag->nRanks[d]; ++i) {
275: 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");
276: l[d][i] = stag->l[d][i] / stag->refineFactor[d]; /* Just divide everything */
277: }
278: }
279: PetscCall(DMStagSetOwnershipRanges(*dmc, l[0], l[1], l[2]));
280: for (d = 0; d < dim; ++d) PetscCall(PetscFree(l[d]));
281: PetscCall(DMSetUp(*dmc));
283: if (dm->coordinates[0].dm) { /* Note that with product coordinates, dm->coordinates = NULL, so we check the DM */
284: DM coordinate_dm, coordinate_dmc;
285: PetscBool isstag, isprod;
287: PetscCall(DMGetCoordinateDM(dm, &coordinate_dm));
288: PetscCall(PetscObjectTypeCompare((PetscObject)coordinate_dm, DMSTAG, &isstag));
289: PetscCall(PetscObjectTypeCompare((PetscObject)coordinate_dm, DMPRODUCT, &isprod));
290: if (isstag) {
291: PetscCall(DMStagSetUniformCoordinatesExplicit(*dmc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); /* Coordinates will be overwritten */
292: PetscCall(DMGetCoordinateDM(*dmc, &coordinate_dmc));
293: PetscCall(DMStagRestrictSimple(coordinate_dm, dm->coordinates[0].x, coordinate_dmc, (*dmc)->coordinates[0].x));
294: } else if (isprod) {
295: PetscCall(DMStagSetUniformCoordinatesProduct(*dmc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); /* Coordinates will be overwritten */
296: PetscCall(DMGetCoordinateDM(*dmc, &coordinate_dmc));
297: for (d = 0; d < dim; ++d) {
298: DM subdm_coarse, subdm_coord_coarse, subdm_fine, subdm_coord_fine;
300: PetscCall(DMProductGetDM(coordinate_dm, d, &subdm_fine));
301: PetscCall(DMGetCoordinateDM(subdm_fine, &subdm_coord_fine));
302: PetscCall(DMProductGetDM(coordinate_dmc, d, &subdm_coarse));
303: PetscCall(DMGetCoordinateDM(subdm_coarse, &subdm_coord_coarse));
304: PetscCall(DMStagRestrictSimple(subdm_coord_fine, subdm_fine->coordinates[0].xl, subdm_coord_coarse, subdm_coarse->coordinates[0].xl));
305: }
306: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown coordinate DM type");
307: }
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: static PetscErrorCode DMRefine_Stag(DM dm, MPI_Comm comm, DM *dmf)
312: {
313: const DM_Stag *const stag = (DM_Stag *)dm->data;
314: PetscInt dim, i, d;
315: PetscInt *l[DMSTAG_MAX_DIM];
317: PetscFunctionBegin;
318: PetscCall(DMStagDuplicateWithoutSetup(dm, comm, dmf));
319: PetscCall(DMSetOptionsPrefix(*dmf, ((PetscObject)dm)->prefix));
320: PetscCall(DMStagSetGlobalSizes(*dmf, stag->N[0] * stag->refineFactor[0], stag->N[1] * stag->refineFactor[1], stag->N[2] * stag->refineFactor[2]));
321: PetscCall(DMGetDimension(dm, &dim));
322: for (d = 0; d < dim; ++d) {
323: PetscCall(PetscMalloc1(stag->nRanks[d], &l[d]));
324: for (i = 0; i < stag->nRanks[d]; ++i) l[d][i] = stag->l[d][i] * stag->refineFactor[d]; /* Just multiply everything */
325: }
326: PetscCall(DMStagSetOwnershipRanges(*dmf, l[0], l[1], l[2]));
327: for (d = 0; d < dim; ++d) PetscCall(PetscFree(l[d]));
328: PetscCall(DMSetUp(*dmf));
329: /* Note: For now, we do not refine coordinates */
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode DMDestroy_Stag(DM dm)
334: {
335: DM_Stag *stag;
336: PetscInt i;
338: PetscFunctionBegin;
339: stag = (DM_Stag *)dm->data;
340: for (i = 0; i < DMSTAG_MAX_DIM; ++i) PetscCall(PetscFree(stag->l[i]));
341: PetscCall(VecScatterDestroy(&stag->gtol));
342: PetscCall(VecScatterDestroy(&stag->ltog_injective));
343: PetscCall(VecScatterDestroy(&stag->ltol));
344: PetscCall(PetscFree(stag->neighbors));
345: PetscCall(PetscFree(stag->locationOffsets));
346: PetscCall(PetscFree(stag->coordinateDMType));
347: PetscCall(PetscFree(stag));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: static PetscErrorCode DMCreateGlobalVector_Stag(DM dm, Vec *vec)
352: {
353: DM_Stag *const stag = (DM_Stag *)dm->data;
355: PetscFunctionBegin;
356: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
357: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), vec));
358: PetscCall(VecSetSizes(*vec, stag->entries, PETSC_DETERMINE));
359: PetscCall(VecSetType(*vec, dm->vectype));
360: PetscCall(VecSetDM(*vec, dm));
361: /* Could set some ops, as DMDA does */
362: PetscCall(VecSetLocalToGlobalMapping(*vec, dm->ltogmap));
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: static PetscErrorCode DMCreateLocalVector_Stag(DM dm, Vec *vec)
367: {
368: DM_Stag *const stag = (DM_Stag *)dm->data;
370: PetscFunctionBegin;
371: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
372: PetscCall(VecCreate(PETSC_COMM_SELF, vec));
373: PetscCall(VecSetSizes(*vec, stag->entriesGhost, PETSC_DETERMINE));
374: PetscCall(VecSetType(*vec, dm->vectype));
375: PetscCall(VecSetBlockSize(*vec, stag->entriesPerElement));
376: PetscCall(VecSetDM(*vec, dm));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /* Helper function to check for the limited situations for which interpolation
381: and restriction functions are implemented */
382: static PetscErrorCode CheckTransferOperatorRequirements_Private(DM dmc, DM dmf)
383: {
384: PetscInt dim, stencilWidthc, stencilWidthf, nf[DMSTAG_MAX_DIM], nc[DMSTAG_MAX_DIM], doff[DMSTAG_MAX_STRATA], dofc[DMSTAG_MAX_STRATA];
386: PetscFunctionBegin;
387: PetscCall(DMGetDimension(dmc, &dim));
388: PetscCall(DMStagGetStencilWidth(dmc, &stencilWidthc));
389: PetscCheck(stencilWidthc >= 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "DMCreateRestriction not implemented for coarse grid stencil width < 1");
390: PetscCall(DMStagGetStencilWidth(dmf, &stencilWidthf));
391: PetscCheck(stencilWidthf >= 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "DMCreateRestriction not implemented for fine grid stencil width < 1");
392: PetscCall(DMStagGetLocalSizes(dmf, &nf[0], &nf[1], &nf[2]));
393: PetscCall(DMStagGetLocalSizes(dmc, &nc[0], &nc[1], &nc[2]));
394: 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");
395: PetscCall(DMStagGetDOF(dmc, &dofc[0], &dofc[1], &dofc[2], &dofc[3]));
396: PetscCall(DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]));
397: for (PetscInt d = 0; d < dim + 1; ++d)
398: 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]);
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /* Since the interpolation uses MATMAIJ for dof > 0 we convert requests for non-MATAIJ baseded matrices to MATAIJ.
403: This is a bit of a hack; the reason for it is partially because -dm_mat_type defines the
404: matrix type for both the operator matrices and the interpolation matrices so that users
405: can select matrix types of base MATAIJ for accelerators
407: Note: The ConvertToAIJ() code below *has been copied from dainterp.c*! ConvertToAIJ() should perhaps be placed somewhere
408: in mat/utils to avoid code duplication, but then the DMStag and DMDA code would need to include the private Mat headers.
409: Since it is only used in two places, I have simply duplicated the code to avoid the need to exposure the private
410: Mat routines in parts of DM. If we find a need for ConvertToAIJ() elsewhere, then we should consolidate it to one
411: place in mat/utils.
412: */
413: static PetscErrorCode ConvertToAIJ(MatType intype, MatType *outtype)
414: {
415: PetscInt i;
416: char const *types[3] = {MATAIJ, MATSEQAIJ, MATMPIAIJ};
417: PetscBool flg;
419: PetscFunctionBegin;
420: *outtype = MATAIJ;
421: for (i = 0; i < 3; i++) {
422: PetscCall(PetscStrbeginswith(intype, types[i], &flg));
423: if (flg) {
424: *outtype = intype;
425: break;
426: }
427: }
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: static PetscErrorCode DMCreateInterpolation_Stag(DM dmc, DM dmf, Mat *A, Vec *vec)
432: {
433: PetscInt dim, entriesf, entriesc;
434: ISLocalToGlobalMapping ltogmf, ltogmc;
435: MatType mattype;
437: PetscFunctionBegin;
438: PetscCall(CheckTransferOperatorRequirements_Private(dmc, dmf));
440: PetscCall(DMStagGetEntries(dmf, &entriesf));
441: PetscCall(DMStagGetEntries(dmc, &entriesc));
442: PetscCall(DMGetLocalToGlobalMapping(dmf, <ogmf));
443: PetscCall(DMGetLocalToGlobalMapping(dmc, <ogmc));
445: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), A));
446: PetscCall(MatSetSizes(*A, entriesf, entriesc, PETSC_DECIDE, PETSC_DECIDE));
447: PetscCall(ConvertToAIJ(dmc->mattype, &mattype));
448: PetscCall(MatSetType(*A, mattype));
449: PetscCall(MatSetLocalToGlobalMapping(*A, ltogmf, ltogmc));
451: PetscCall(DMGetDimension(dmc, &dim));
452: if (dim == 1) PetscCall(DMStagPopulateInterpolation1d_Internal(dmc, dmf, *A));
453: else if (dim == 2) PetscCall(DMStagPopulateInterpolation2d_Internal(dmc, dmf, *A));
454: else if (dim == 3) PetscCall(DMStagPopulateInterpolation3d_Internal(dmc, dmf, *A));
455: else SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
456: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
457: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
459: if (vec) *vec = NULL;
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: static PetscErrorCode DMCreateRestriction_Stag(DM dmc, DM dmf, Mat *A)
464: {
465: PetscInt dim, entriesf, entriesc, doff[DMSTAG_MAX_STRATA];
466: ISLocalToGlobalMapping ltogmf, ltogmc;
467: MatType mattype;
469: PetscFunctionBegin;
470: PetscCall(CheckTransferOperatorRequirements_Private(dmc, dmf));
472: PetscCall(DMStagGetEntries(dmf, &entriesf));
473: PetscCall(DMStagGetEntries(dmc, &entriesc));
474: PetscCall(DMGetLocalToGlobalMapping(dmf, <ogmf));
475: PetscCall(DMGetLocalToGlobalMapping(dmc, <ogmc));
477: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), A));
478: PetscCall(MatSetSizes(*A, entriesc, entriesf, PETSC_DECIDE, PETSC_DECIDE)); /* Note transpose wrt interpolation */
479: PetscCall(ConvertToAIJ(dmc->mattype, &mattype));
480: PetscCall(MatSetType(*A, mattype));
481: PetscCall(MatSetLocalToGlobalMapping(*A, ltogmc, ltogmf)); /* Note transpose wrt interpolation */
483: PetscCall(DMGetDimension(dmc, &dim));
484: PetscCall(DMStagGetDOF(dmf, &doff[0], &doff[1], &doff[2], &doff[3]));
485: if (dim == 1) PetscCall(DMStagPopulateRestriction1d_Internal(dmc, dmf, *A));
486: else if (dim == 2) PetscCall(DMStagPopulateRestriction2d_Internal(dmc, dmf, *A));
487: else if (dim == 3) PetscCall(DMStagPopulateRestriction3d_Internal(dmc, dmf, *A));
488: else SETERRQ(PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
490: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
491: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: static PetscErrorCode DMCreateMatrix_Stag(DM dm, Mat *mat)
496: {
497: MatType mat_type;
498: PetscBool is_shell, is_aij;
499: PetscInt dim, entries;
500: ISLocalToGlobalMapping ltogmap;
502: PetscFunctionBegin;
503: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This function must be called after DMSetUp()");
504: PetscCall(DMGetDimension(dm, &dim));
505: PetscCall(DMGetMatType(dm, &mat_type));
506: PetscCall(DMStagGetEntries(dm, &entries));
507: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), mat));
508: PetscCall(MatSetSizes(*mat, entries, entries, PETSC_DETERMINE, PETSC_DETERMINE));
509: PetscCall(MatSetType(*mat, mat_type));
510: PetscCall(MatSetUp(*mat));
511: PetscCall(DMGetLocalToGlobalMapping(dm, <ogmap));
512: PetscCall(MatSetLocalToGlobalMapping(*mat, ltogmap, ltogmap));
513: PetscCall(MatSetDM(*mat, dm));
515: /* Compare to similar and perhaps superior logic in DMCreateMatrix_DA, which creates
516: the matrix first and then performs this logic by checking for preallocation functions */
517: PetscCall(PetscObjectBaseTypeCompare((PetscObject)*mat, MATAIJ, &is_aij));
518: if (!is_aij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)*mat, MATSEQAIJ, &is_aij));
519: if (!is_aij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)*mat, MATMPIAIJ, &is_aij));
520: PetscCall(PetscStrcmp(mat_type, MATSHELL, &is_shell));
521: if (is_aij) {
522: Mat preallocator;
523: PetscInt m, n;
524: const PetscBool fill_with_zeros = PETSC_FALSE;
526: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), &preallocator));
527: PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
528: PetscCall(MatGetLocalSize(*mat, &m, &n));
529: PetscCall(MatSetSizes(preallocator, m, n, PETSC_DECIDE, PETSC_DECIDE));
530: PetscCall(MatSetLocalToGlobalMapping(preallocator, ltogmap, ltogmap));
531: PetscCall(MatSetUp(preallocator));
532: switch (dim) {
533: case 1:
534: PetscCall(DMCreateMatrix_Stag_1D_AIJ_Assemble(dm, preallocator));
535: break;
536: case 2:
537: PetscCall(DMCreateMatrix_Stag_2D_AIJ_Assemble(dm, preallocator));
538: break;
539: case 3:
540: PetscCall(DMCreateMatrix_Stag_3D_AIJ_Assemble(dm, preallocator));
541: break;
542: default:
543: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
544: }
545: PetscCall(MatPreallocatorPreallocate(preallocator, fill_with_zeros, *mat));
546: PetscCall(MatDestroy(&preallocator));
548: if (!dm->prealloc_only) {
549: /* Bind to CPU before assembly, to prevent unnecessary copies of zero entries from CPU to GPU */
550: PetscCall(MatBindToCPU(*mat, PETSC_TRUE));
551: switch (dim) {
552: case 1:
553: PetscCall(DMCreateMatrix_Stag_1D_AIJ_Assemble(dm, *mat));
554: break;
555: case 2:
556: PetscCall(DMCreateMatrix_Stag_2D_AIJ_Assemble(dm, *mat));
557: break;
558: case 3:
559: PetscCall(DMCreateMatrix_Stag_3D_AIJ_Assemble(dm, *mat));
560: break;
561: default:
562: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
563: }
564: PetscCall(MatBindToCPU(*mat, PETSC_FALSE));
565: }
566: } else if (is_shell) {
567: /* nothing more to do */
568: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for Mattype %s", mat_type);
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: static PetscErrorCode DMGetCompatibility_Stag(DM dm, DM dm2, PetscBool *compatible, PetscBool *set)
573: {
574: const DM_Stag *const stag = (DM_Stag *)dm->data;
575: const DM_Stag *const stag2 = (DM_Stag *)dm2->data;
576: PetscInt dim, dim2, i;
577: MPI_Comm comm;
578: PetscMPIInt sameComm;
579: DMType type2;
580: PetscBool sameType;
582: PetscFunctionBegin;
583: PetscCall(DMGetType(dm2, &type2));
584: PetscCall(PetscStrcmp(DMSTAG, type2, &sameType));
585: if (!sameType) {
586: PetscCall(PetscInfo((PetscObject)dm, "DMStag compatibility check not implemented with DM of type %s\n", type2));
587: *set = PETSC_FALSE;
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
592: PetscCallMPI(MPI_Comm_compare(comm, PetscObjectComm((PetscObject)dm2), &sameComm));
593: if (sameComm != MPI_IDENT) {
594: 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)));
595: *set = PETSC_FALSE;
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
598: PetscCall(DMGetDimension(dm, &dim));
599: PetscCall(DMGetDimension(dm2, &dim2));
600: if (dim != dim2) {
601: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different dimensions\n"));
602: *set = PETSC_TRUE;
603: *compatible = PETSC_FALSE;
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
606: for (i = 0; i < dim; ++i) {
607: if (stag->N[i] != stag2->N[i]) {
608: 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]));
609: *set = PETSC_TRUE;
610: *compatible = PETSC_FALSE;
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
613: if (stag->n[i] != stag2->n[i]) {
614: 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]));
615: *set = PETSC_TRUE;
616: *compatible = PETSC_FALSE;
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
619: if (stag->boundaryType[i] != stag2->boundaryType[i]) {
620: 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]]));
621: *set = PETSC_TRUE;
622: *compatible = PETSC_FALSE;
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
625: }
626: /* Note: we include stencil type and width in the notion of compatibility, as this affects
627: the "atlas" (local subdomains). This might be irritating in legitimate cases
628: of wanting to transfer between two other-wise compatible DMs with different
629: stencil characteristics. */
630: if (stag->stencilType != stag2->stencilType) {
631: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different ghost stencil types: %s != %s\n", DMStagStencilTypes[stag->stencilType], DMStagStencilTypes[stag2->stencilType]));
632: *set = PETSC_TRUE;
633: *compatible = PETSC_FALSE;
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
636: if (stag->stencilWidth != stag2->stencilWidth) {
637: PetscCall(PetscInfo((PetscObject)dm, "DMStag objects have different ghost stencil widths: %" PetscInt_FMT " != %" PetscInt_FMT "\n", stag->stencilWidth, stag->stencilWidth));
638: *set = PETSC_TRUE;
639: *compatible = PETSC_FALSE;
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
642: *set = PETSC_TRUE;
643: *compatible = PETSC_TRUE;
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }
647: static PetscErrorCode DMHasCreateInjection_Stag(DM dm, PetscBool *flg)
648: {
649: PetscFunctionBegin;
651: PetscAssertPointer(flg, 2);
652: *flg = PETSC_FALSE;
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: /*
657: Note there are several orderings in play here.
658: 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).
659: Also in all cases, only subdomains which are the last in their dimension have partial elements.
661: 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.
662: 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.
663: 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.
664: */
666: static PetscErrorCode DMLocalToGlobalBegin_Stag(DM dm, Vec l, InsertMode mode, Vec g)
667: {
668: DM_Stag *const stag = (DM_Stag *)dm->data;
670: PetscFunctionBegin;
671: if (mode == ADD_VALUES) {
672: PetscCall(VecScatterBegin(stag->gtol, l, g, mode, SCATTER_REVERSE));
673: } else if (mode == INSERT_VALUES) {
674: if (stag->ltog_injective) {
675: PetscCall(VecScatterBegin(stag->ltog_injective, l, g, mode, SCATTER_FORWARD));
676: } else {
677: PetscCall(VecScatterBegin(stag->gtol, l, g, mode, SCATTER_REVERSE_LOCAL));
678: }
679: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported InsertMode");
680: PetscFunctionReturn(PETSC_SUCCESS);
681: }
683: static PetscErrorCode DMLocalToGlobalEnd_Stag(DM dm, Vec l, InsertMode mode, Vec g)
684: {
685: DM_Stag *const stag = (DM_Stag *)dm->data;
687: PetscFunctionBegin;
688: if (mode == ADD_VALUES) {
689: PetscCall(VecScatterEnd(stag->gtol, l, g, mode, SCATTER_REVERSE));
690: } else if (mode == INSERT_VALUES) {
691: if (stag->ltog_injective) {
692: PetscCall(VecScatterEnd(stag->ltog_injective, l, g, mode, SCATTER_FORWARD));
693: } else {
694: PetscCall(VecScatterEnd(stag->gtol, l, g, mode, SCATTER_REVERSE_LOCAL));
695: }
696: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported InsertMode");
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: static PetscErrorCode DMGlobalToLocalBegin_Stag(DM dm, Vec g, InsertMode mode, Vec l)
701: {
702: DM_Stag *const stag = (DM_Stag *)dm->data;
704: PetscFunctionBegin;
705: PetscCall(VecScatterBegin(stag->gtol, g, l, mode, SCATTER_FORWARD));
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: static PetscErrorCode DMGlobalToLocalEnd_Stag(DM dm, Vec g, InsertMode mode, Vec l)
710: {
711: DM_Stag *const stag = (DM_Stag *)dm->data;
713: PetscFunctionBegin;
714: PetscCall(VecScatterEnd(stag->gtol, g, l, mode, SCATTER_FORWARD));
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
718: static PetscErrorCode DMLocalToLocalBegin_Stag(DM dm, Vec g, InsertMode mode, Vec l)
719: {
720: DM_Stag *const stag = (DM_Stag *)dm->data;
722: PetscFunctionBegin;
723: if (!stag->ltol) {
724: PetscInt dim;
725: PetscCall(DMGetDimension(dm, &dim));
726: switch (dim) {
727: case 1:
728: PetscCall(DMStagPopulateLocalToLocal1d_Internal(dm));
729: break;
730: case 2:
731: PetscCall(DMStagPopulateLocalToLocal2d_Internal(dm));
732: break;
733: case 3:
734: PetscCall(DMStagPopulateLocalToLocal3d_Internal(dm));
735: break;
736: default:
737: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
738: }
739: }
740: PetscCall(VecScatterBegin(stag->ltol, g, l, mode, SCATTER_FORWARD));
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: static PetscErrorCode DMLocalToLocalEnd_Stag(DM dm, Vec g, InsertMode mode, Vec l)
745: {
746: DM_Stag *const stag = (DM_Stag *)dm->data;
748: PetscFunctionBegin;
749: PetscCall(VecScatterEnd(stag->ltol, g, l, mode, SCATTER_FORWARD));
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: /*
754: If a stratum is active (non-zero dof), make it active in the coordinate DM.
755: */
756: static PetscErrorCode DMCreateCoordinateDM_Stag(DM dm, DM *dmc)
757: {
758: DM_Stag *const stag = (DM_Stag *)dm->data;
759: PetscInt dim;
760: PetscBool isstag, isproduct;
761: const char *prefix;
763: PetscFunctionBegin;
764: PetscCheck(stag->coordinateDMType, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Before creating a coordinate DM, a type must be specified with DMStagSetCoordinateDMType()");
766: PetscCall(DMGetDimension(dm, &dim));
767: PetscCall(PetscStrcmp(stag->coordinateDMType, DMSTAG, &isstag));
768: PetscCall(PetscStrcmp(stag->coordinateDMType, DMPRODUCT, &isproduct));
769: if (isstag) {
770: 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));
771: } else if (isproduct) {
772: PetscCall(DMCreate(PETSC_COMM_WORLD, dmc));
773: PetscCall(DMSetType(*dmc, DMPRODUCT));
774: PetscCall(DMSetDimension(*dmc, dim));
775: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate DM type %s", stag->coordinateDMType);
776: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
777: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dmc, prefix));
778: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*dmc, "cdm_"));
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
782: static PetscErrorCode DMGetNeighbors_Stag(DM dm, PetscInt *nRanks, const PetscMPIInt *ranks[])
783: {
784: DM_Stag *const stag = (DM_Stag *)dm->data;
785: PetscInt dim;
787: PetscFunctionBegin;
788: PetscCall(DMGetDimension(dm, &dim));
789: switch (dim) {
790: case 1:
791: *nRanks = 3;
792: break;
793: case 2:
794: *nRanks = 9;
795: break;
796: case 3:
797: *nRanks = 27;
798: break;
799: default:
800: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Get neighbors not implemented for dim = %" PetscInt_FMT, dim);
801: }
802: *ranks = stag->neighbors;
803: PetscFunctionReturn(PETSC_SUCCESS);
804: }
806: static PetscErrorCode DMView_Stag(DM dm, PetscViewer viewer)
807: {
808: DM_Stag *const stag = (DM_Stag *)dm->data;
809: PetscBool isascii, viewAllRanks;
810: PetscMPIInt rank, size;
811: PetscInt dim, maxRanksToView, i;
813: PetscFunctionBegin;
814: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
815: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
816: PetscCall(DMGetDimension(dm, &dim));
817: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
818: if (isascii) {
819: PetscCall(PetscViewerASCIIPrintf(viewer, "Dimension: %" PetscInt_FMT "\n", dim));
820: switch (dim) {
821: case 1:
822: PetscCall(PetscViewerASCIIPrintf(viewer, "Global size: %" PetscInt_FMT "\n", stag->N[0]));
823: break;
824: case 2:
825: PetscCall(PetscViewerASCIIPrintf(viewer, "Global sizes: %" PetscInt_FMT " x %" PetscInt_FMT "\n", stag->N[0], stag->N[1]));
826: PetscCall(PetscViewerASCIIPrintf(viewer, "Parallel decomposition: %d x %d ranks\n", stag->nRanks[0], stag->nRanks[1]));
827: break;
828: case 3:
829: PetscCall(PetscViewerASCIIPrintf(viewer, "Global sizes: %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "\n", stag->N[0], stag->N[1], stag->N[2]));
830: PetscCall(PetscViewerASCIIPrintf(viewer, "Parallel decomposition: %d x %d x %d ranks\n", stag->nRanks[0], stag->nRanks[1], stag->nRanks[2]));
831: break;
832: default:
833: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim==%" PetscInt_FMT, dim);
834: }
835: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary ghosting:"));
836: for (i = 0; i < dim; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, " %s", DMBoundaryTypes[stag->boundaryType[i]]));
837: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
838: PetscCall(PetscViewerASCIIPrintf(viewer, "Elementwise ghost stencil: %s", DMStagStencilTypes[stag->stencilType]));
839: if (stag->stencilType != DMSTAG_STENCIL_NONE) {
840: PetscCall(PetscViewerASCIIPrintf(viewer, ", width %" PetscInt_FMT "\n", stag->stencilWidth));
841: } else {
842: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
843: }
844: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per vertex (0D)\n", stag->dof[0]));
845: if (dim == 3) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per edge (1D)\n", stag->dof[1]));
846: if (dim > 1) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per face (%" PetscInt_FMT "D)\n", stag->dof[dim - 1], dim - 1));
847: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " DOF per element (%" PetscInt_FMT "D)\n", stag->dof[dim], dim));
848: if (dm->coordinates[0].dm) PetscCall(PetscViewerASCIIPrintf(viewer, "Has coordinate DM\n"));
849: maxRanksToView = 16;
850: viewAllRanks = (PetscBool)(size <= maxRanksToView);
851: if (viewAllRanks) {
852: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
853: switch (dim) {
854: case 1:
855: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local elements : %" PetscInt_FMT " (%" PetscInt_FMT " with ghosts)\n", rank, stag->n[0], stag->nGhost[0]));
856: break;
857: case 2:
858: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Rank coordinates (%d,%d)\n", rank, stag->rank[0], stag->rank[1]));
859: 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]));
860: break;
861: case 3:
862: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Rank coordinates (%d,%d,%d)\n", rank, stag->rank[0], stag->rank[1], stag->rank[2]));
863: 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],
864: stag->n[2], stag->nGhost[0], stag->nGhost[1], stag->nGhost[2]));
865: break;
866: default:
867: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "not implemented for dim==%" PetscInt_FMT, dim);
868: }
869: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local native entries: %" PetscInt_FMT "\n", rank, stag->entries));
870: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local entries total : %" PetscInt_FMT "\n", rank, stag->entriesGhost));
871: PetscCall(PetscViewerFlush(viewer));
872: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
873: } else {
874: PetscCall(PetscViewerASCIIPrintf(viewer, "(Per-rank information omitted since >%" PetscInt_FMT " ranks used)\n", maxRanksToView));
875: }
876: }
877: PetscFunctionReturn(PETSC_SUCCESS);
878: }
880: static PetscErrorCode DMSetFromOptions_Stag(DM dm, PetscOptionItems *PetscOptionsObject)
881: {
882: DM_Stag *const stag = (DM_Stag *)dm->data;
883: PetscInt dim, nRefine = 0, refineFactorTotal[DMSTAG_MAX_DIM], i, d;
885: PetscFunctionBegin;
886: PetscCall(DMGetDimension(dm, &dim));
887: PetscOptionsHeadBegin(PetscOptionsObject, "DMStag Options");
888: PetscCall(PetscOptionsInt("-stag_grid_x", "Number of grid points in x direction", "DMStagSetGlobalSizes", stag->N[0], &stag->N[0], NULL));
889: if (dim > 1) PetscCall(PetscOptionsInt("-stag_grid_y", "Number of grid points in y direction", "DMStagSetGlobalSizes", stag->N[1], &stag->N[1], NULL));
890: if (dim > 2) PetscCall(PetscOptionsInt("-stag_grid_z", "Number of grid points in z direction", "DMStagSetGlobalSizes", stag->N[2], &stag->N[2], NULL));
891: PetscCall(PetscOptionsMPIInt("-stag_ranks_x", "Number of ranks in x direction", "DMStagSetNumRanks", stag->nRanks[0], &stag->nRanks[0], NULL));
892: if (dim > 1) PetscCall(PetscOptionsMPIInt("-stag_ranks_y", "Number of ranks in y direction", "DMStagSetNumRanks", stag->nRanks[1], &stag->nRanks[1], NULL));
893: if (dim > 2) PetscCall(PetscOptionsMPIInt("-stag_ranks_z", "Number of ranks in z direction", "DMStagSetNumRanks", stag->nRanks[2], &stag->nRanks[2], NULL));
894: PetscCall(PetscOptionsInt("-stag_stencil_width", "Elementwise stencil width", "DMStagSetStencilWidth", stag->stencilWidth, &stag->stencilWidth, NULL));
895: PetscCall(PetscOptionsEnum("-stag_stencil_type", "Elementwise stencil stype", "DMStagSetStencilType", DMStagStencilTypes, (PetscEnum)stag->stencilType, (PetscEnum *)&stag->stencilType, NULL));
896: PetscCall(PetscOptionsEnum("-stag_boundary_type_x", "Treatment of (physical) boundaries in x direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[0], (PetscEnum *)&stag->boundaryType[0], NULL));
897: PetscCall(PetscOptionsEnum("-stag_boundary_type_y", "Treatment of (physical) boundaries in y direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[1], (PetscEnum *)&stag->boundaryType[1], NULL));
898: PetscCall(PetscOptionsEnum("-stag_boundary_type_z", "Treatment of (physical) boundaries in z direction", "DMStagSetBoundaryTypes", DMBoundaryTypes, (PetscEnum)stag->boundaryType[2], (PetscEnum *)&stag->boundaryType[2], NULL));
899: PetscCall(PetscOptionsInt("-stag_dof_0", "Number of dof per 0-cell (vertex)", "DMStagSetDOF", stag->dof[0], &stag->dof[0], NULL));
900: 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));
901: 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));
902: PetscCall(PetscOptionsInt("-stag_dof_3", "Number of dof per 3-cell (element in 3D)", "DMStagSetDOF", stag->dof[3], &stag->dof[3], NULL));
903: PetscCall(PetscOptionsBoundedInt("-stag_refine_x", "Refinement factor in x-direction", "DMStagSetRefinementFactor", stag->refineFactor[0], &stag->refineFactor[0], NULL, 1));
904: if (dim > 1) PetscCall(PetscOptionsBoundedInt("-stag_refine_y", "Refinement factor in y-direction", "DMStagSetRefinementFactor", stag->refineFactor[1], &stag->refineFactor[1], NULL, 1));
905: if (dim > 2) PetscCall(PetscOptionsBoundedInt("-stag_refine_z", "Refinement factor in z-direction", "DMStagSetRefinementFactor", stag->refineFactor[2], &stag->refineFactor[2], NULL, 1));
906: PetscCall(PetscOptionsBoundedInt("-stag_refine", "Refine grid one or more times", "None", nRefine, &nRefine, NULL, 0));
907: PetscOptionsHeadEnd();
909: for (d = 0; d < dim; ++d) refineFactorTotal[d] = 1;
910: while (nRefine--)
911: for (d = 0; d < dim; ++d) refineFactorTotal[d] *= stag->refineFactor[d];
912: for (d = 0; d < dim; ++d) {
913: stag->N[d] *= refineFactorTotal[d];
914: if (stag->l[d])
915: for (i = 0; i < stag->nRanks[d]; ++i) stag->l[d][i] *= refineFactorTotal[d];
916: }
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: /*MC
921: DMSTAG - `"stag"` - A `DM` object representing a "staggered grid" or a structured cell complex.
923: Level: beginner
925: Notes:
926: This implementation parallels the `DMDA` implementation in many ways, but allows degrees of freedom
927: to be associated with all "strata" in a logically-rectangular grid.
929: Each stratum can be characterized by the dimension of the entities ("points", to borrow the `DMPLEX`
930: terminology), from 0- to 3-dimensional.
932: In some cases this numbering is used directly, for example with `DMStagGetDOF()`.
933: To allow easier reading and to some extent more similar code between different-dimensional implementations
934: of the same problem, we associate canonical names for each type of point, for each dimension of DMStag.
936: * 1-dimensional `DMSTAG` objects have vertices (0D) and elements (1D).
937: * 2-dimensional `DMSTAG` objects have vertices (0D), faces (1D), and elements (2D).
938: * 3-dimensional `DMSTAG` objects have vertices (0D), edges (1D), faces (2D), and elements (3D).
940: This naming is reflected when viewing a `DMSTAG` object with `DMView()`, and in forming
941: convenient options prefixes when creating a decomposition with `DMCreateFieldDecomposition()`.
943: .seealso: [](ch_stag), `DM`, `DMPRODUCT`, `DMDA`, `DMPLEX`, `DMStagCreate1d()`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMType`, `DMCreate()`,
944: `DMSetType()`, `DMStagVecSplitToDMDA()`
945: M*/
947: PETSC_EXTERN PetscErrorCode DMCreate_Stag(DM dm)
948: {
949: DM_Stag *stag;
950: PetscInt i, dim;
952: PetscFunctionBegin;
953: PetscAssertPointer(dm, 1);
954: PetscCall(PetscNew(&stag));
955: dm->data = stag;
957: stag->gtol = NULL;
958: stag->ltog_injective = NULL;
959: stag->ltol = NULL;
960: for (i = 0; i < DMSTAG_MAX_STRATA; ++i) stag->dof[i] = 0;
961: for (i = 0; i < DMSTAG_MAX_DIM; ++i) stag->l[i] = NULL;
962: stag->stencilType = DMSTAG_STENCIL_NONE;
963: stag->stencilWidth = 0;
964: for (i = 0; i < DMSTAG_MAX_DIM; ++i) stag->nRanks[i] = -1;
965: stag->coordinateDMType = NULL;
966: for (i = 0; i < DMSTAG_MAX_DIM; ++i) stag->refineFactor[i] = 2;
968: PetscCall(DMGetDimension(dm, &dim));
969: 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");
971: PetscCall(PetscMemzero(dm->ops, sizeof(*dm->ops)));
972: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Stag;
973: dm->ops->createglobalvector = DMCreateGlobalVector_Stag;
974: dm->ops->createlocalvector = DMCreateLocalVector_Stag;
975: dm->ops->creatematrix = DMCreateMatrix_Stag;
976: dm->ops->hascreateinjection = DMHasCreateInjection_Stag;
977: dm->ops->refine = DMRefine_Stag;
978: dm->ops->coarsen = DMCoarsen_Stag;
979: dm->ops->createinterpolation = DMCreateInterpolation_Stag;
980: dm->ops->createrestriction = DMCreateRestriction_Stag;
981: dm->ops->destroy = DMDestroy_Stag;
982: dm->ops->getneighbors = DMGetNeighbors_Stag;
983: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Stag;
984: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Stag;
985: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Stag;
986: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Stag;
987: dm->ops->localtolocalbegin = DMLocalToLocalBegin_Stag;
988: dm->ops->localtolocalend = DMLocalToLocalEnd_Stag;
989: dm->ops->setfromoptions = DMSetFromOptions_Stag;
990: switch (dim) {
991: case 1:
992: dm->ops->setup = DMSetUp_Stag_1d;
993: break;
994: case 2:
995: dm->ops->setup = DMSetUp_Stag_2d;
996: break;
997: case 3:
998: dm->ops->setup = DMSetUp_Stag_3d;
999: break;
1000: default:
1001: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
1002: }
1003: dm->ops->clone = DMClone_Stag;
1004: dm->ops->view = DMView_Stag;
1005: dm->ops->getcompatibility = DMGetCompatibility_Stag;
1006: dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Stag;
1007: PetscFunctionReturn(PETSC_SUCCESS);
1008: }