Actual source code: stag1d.c
1: /*
2: Functions specific to the 1-dimensional implementation of DMStag
3: */
4: #include <petsc/private/dmstagimpl.h>
6: /*@C
7: DMStagCreate1d - Create an object to manage data living on the elements and vertices of a parallelized regular 1D grid.
9: Collective
11: Input Parameters:
12: + comm - MPI communicator
13: . bndx - boundary type: `DM_BOUNDARY_NONE`, `DM_BOUNDARY_PERIODIC`, or `DM_BOUNDARY_GHOSTED`
14: . M - global number of elements
15: . dof0 - number of degrees of freedom per vertex/0-cell
16: . dof1 - number of degrees of freedom per element/1-cell
17: . stencilType - ghost/halo region type: `DMSTAG_STENCIL_BOX` or `DMSTAG_STENCIL_NONE`
18: . stencilWidth - width, in elements, of halo/ghost region
19: - lx - array of local sizes, of length equal to the comm size, summing to M
21: Output Parameter:
22: . dm - the new DMStag object
24: Options Database Keys:
25: + -dm_view - calls `DMViewFromOptions()` at the conclusion of `DMSetUp()`
26: . -stag_grid_x <nx> - number of elements in the x direction
27: . -stag_ghost_stencil_width - width of ghost region, in elements
28: - -stag_boundary_type_x <none,ghosted,periodic> - `DMBoundaryType` value
30: Level: beginner
32: Notes:
33: You must call `DMSetUp()` after this call before using the `DM`.
34: If you wish to use the options database (see the keys above) to change values in the `DMSTAG`, you must call
35: `DMSetFromOptions()` after this function but before `DMSetUp()`.
37: .seealso: [](ch_stag), `DMSTAG`, `DMStagCreate2d()`, `DMStagCreate3d()`, `DMDestroy()`, `DMView()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMLocalToGlobalBegin()`, `DMDACreate1d()`
38: @*/
39: PETSC_EXTERN PetscErrorCode DMStagCreate1d(MPI_Comm comm, DMBoundaryType bndx, PetscInt M, PetscInt dof0, PetscInt dof1, DMStagStencilType stencilType, PetscInt stencilWidth, const PetscInt lx[], DM *dm)
40: {
41: PetscMPIInt size;
43: PetscFunctionBegin;
44: PetscCallMPI(MPI_Comm_size(comm, &size));
45: PetscCall(DMCreate(comm, dm));
46: PetscCall(DMSetDimension(*dm, 1));
47: PetscCall(DMStagInitialize(bndx, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, M, 0, 0, size, 0, 0, dof0, dof1, 0, 0, stencilType, stencilWidth, lx, NULL, NULL, *dm));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: PETSC_INTERN PetscErrorCode DMStagRestrictSimple_1d(DM dmf, Vec xf_local, DM dmc, Vec xc_local)
52: {
53: PetscScalar **LA_xf, **LA_xc;
54: PetscInt i, start, n, nextra, N;
55: PetscInt d, dof[2];
56: PetscInt slot_left_coarse, slot_element_coarse, slot_left_fine, slot_element_fine;
58: PetscFunctionBegin;
59: PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], NULL, NULL));
60: PetscCall(DMStagGetCorners(dmc, &start, NULL, NULL, &n, NULL, NULL, &nextra, NULL, NULL));
61: PetscCall(DMStagGetGlobalSizes(dmc, &N, NULL, NULL));
62: if (PetscDefined(USE_DEBUG)) {
63: PetscInt dof_check[2], n_fine, start_fine;
65: PetscCall(DMStagGetDOF(dmf, &dof_check[0], &dof_check[1], NULL, NULL));
66: PetscCall(DMStagGetCorners(dmf, &start_fine, NULL, NULL, &n_fine, NULL, NULL, NULL, NULL, NULL));
67: for (d = 0; d < 2; ++d)
68: PetscCheck(dof_check[d] == dof[d], PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_INCOMP, "Cannot transfer between DMStag objects with different dof on each stratum. Stratum %" PetscInt_FMT " has %" PetscInt_FMT " dof (fine) but %" PetscInt_FMT " dof (coarse)", d, dof_check[d], dof[d]);
69: PetscCheck(n_fine == 2 * n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot transfer between DMStag objects unless there is a 2-1 coarsening. The fine DM has %" PetscInt_FMT " local elements and the coarse DM has %" PetscInt_FMT "", n_fine, n);
70: PetscCheck(start_fine == 2 * start, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot transfer between DMStag objects unless there is a 2-1 coarsening. The fine DM starts at element %" PetscInt_FMT " and the coarse DM starts at %" PetscInt_FMT "", start_fine, start);
71: {
72: PetscInt size_local, entries_local;
74: PetscCall(DMStagGetEntriesLocal(dmf, &entries_local));
75: PetscCall(VecGetLocalSize(xf_local, &size_local));
76: PetscCheck(entries_local == size_local, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Fine vector must be a local vector of size %" PetscInt_FMT ", but a vector of size %" PetscInt_FMT " was supplied", entries_local, size_local);
77: }
78: {
79: PetscInt size_local, entries_local;
81: PetscCall(DMStagGetEntriesLocal(dmc, &entries_local));
82: PetscCall(VecGetLocalSize(xc_local, &size_local));
83: PetscCheck(entries_local == size_local, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Coarse vector must be a local vector of size %" PetscInt_FMT ", but a vector of size %" PetscInt_FMT " was supplied", entries_local, size_local);
84: }
85: }
86: PetscCall(VecZeroEntries(xc_local));
87: PetscCall(DMStagVecGetArray(dmf, xf_local, &LA_xf));
88: PetscCall(DMStagVecGetArray(dmc, xc_local, &LA_xc));
89: PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_LEFT, 0, &slot_left_fine));
90: PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_ELEMENT, 0, &slot_element_fine));
91: PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_LEFT, 0, &slot_left_coarse));
92: PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_ELEMENT, 0, &slot_element_coarse));
93: for (i = start; i < start + n + nextra; ++i) {
94: for (d = 0; d < dof[0]; ++d) LA_xc[i][slot_left_coarse + d] = LA_xf[2 * i][slot_left_fine + d];
95: if (i < N) {
96: for (d = 0; d < dof[1]; ++d) LA_xc[i][slot_element_coarse + d] = 0.5 * (LA_xf[2 * i][slot_element_fine + d] + LA_xf[2 * i + 1][slot_element_fine + d]);
97: }
98: }
99: PetscCall(DMStagVecRestoreArray(dmf, xf_local, &LA_xf));
100: PetscCall(DMStagVecRestoreArray(dmc, xc_local, &LA_xc));
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_1d(DM dm, PetscReal xmin, PetscReal xmax)
105: {
106: DM_Stag *stagCoord;
107: DM dmCoord;
108: Vec coordLocal;
109: PetscReal h, min;
110: PetscScalar **arr;
111: PetscInt start_ghost, n_ghost, s;
112: PetscInt ileft, ielement;
114: PetscFunctionBegin;
115: PetscCall(DMGetCoordinateDM(dm, &dmCoord));
116: stagCoord = (DM_Stag *)dmCoord->data;
117: for (s = 0; s < 2; ++s) {
118: PetscCheck(stagCoord->dof[s] == 0 || stagCoord->dof[s] == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Coordinate DM in 1 dimensions must have 0 or 1 dof on each stratum, but stratum %" PetscInt_FMT " has %" PetscInt_FMT " dof", s,
119: stagCoord->dof[s]);
120: }
121: PetscCall(DMCreateLocalVector(dmCoord, &coordLocal));
123: PetscCall(DMStagVecGetArray(dmCoord, coordLocal, &arr));
124: if (stagCoord->dof[0]) PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_LEFT, 0, &ileft));
125: if (stagCoord->dof[1]) PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_ELEMENT, 0, &ielement));
126: PetscCall(DMStagGetGhostCorners(dmCoord, &start_ghost, NULL, NULL, &n_ghost, NULL, NULL));
128: min = xmin;
129: h = (xmax - xmin) / stagCoord->N[0];
131: for (PetscInt ind = start_ghost; ind < start_ghost + n_ghost; ++ind) {
132: if (stagCoord->dof[0]) {
133: const PetscReal off = 0.0;
134: arr[ind][ileft] = min + ((PetscReal)ind + off) * h;
135: }
136: if (stagCoord->dof[1]) {
137: const PetscReal off = 0.5;
138: arr[ind][ielement] = min + ((PetscReal)ind + off) * h;
139: }
140: }
141: PetscCall(DMStagVecRestoreArray(dmCoord, coordLocal, &arr));
142: PetscCall(DMSetCoordinatesLocal(dm, coordLocal));
143: PetscCall(VecDestroy(&coordLocal));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /* Helper functions used in DMSetUp_Stag() */
148: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM);
150: PETSC_INTERN PetscErrorCode DMSetUp_Stag_1d(DM dm)
151: {
152: DM_Stag *const stag = (DM_Stag *)dm->data;
153: PetscMPIInt size, rank;
154: MPI_Comm comm;
155: PetscInt j;
157: PetscFunctionBegin;
158: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
159: PetscCallMPI(MPI_Comm_size(comm, &size));
160: PetscCallMPI(MPI_Comm_rank(comm, &rank));
162: /* Check Global size */
163: PetscCheck(stag->N[0] >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Global grid size of %" PetscInt_FMT " < 1 specified", stag->N[0]);
165: /* Local sizes */
166: PetscCheck(stag->N[0] >= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "More ranks (%d) than elements (%" PetscInt_FMT ") specified", size, stag->N[0]);
167: if (!stag->l[0]) {
168: /* Divide equally, giving an extra elements to higher ranks */
169: PetscCall(PetscMalloc1(stag->nRanks[0], &stag->l[0]));
170: for (j = 0; j < stag->nRanks[0]; ++j) stag->l[0][j] = stag->N[0] / stag->nRanks[0] + (stag->N[0] % stag->nRanks[0] > j ? 1 : 0);
171: }
172: {
173: PetscInt Nchk = 0;
174: for (j = 0; j < size; ++j) Nchk += stag->l[0][j];
175: PetscCheck(Nchk == stag->N[0], comm, PETSC_ERR_ARG_OUTOFRANGE, "Sum of specified local sizes (%" PetscInt_FMT ") is not equal to global size (%" PetscInt_FMT ")", Nchk, stag->N[0]);
176: }
177: stag->n[0] = stag->l[0][rank];
179: /* Rank (trivial in 1d) */
180: stag->rank[0] = rank;
181: stag->firstRank[0] = (PetscBool)(rank == 0);
182: stag->lastRank[0] = (PetscBool)(rank == size - 1);
184: /* Local (unghosted) numbers of entries */
185: stag->entriesPerElement = stag->dof[0] + stag->dof[1];
186: switch (stag->boundaryType[0]) {
187: case DM_BOUNDARY_NONE:
188: case DM_BOUNDARY_GHOSTED:
189: stag->entries = stag->n[0] * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
190: break;
191: case DM_BOUNDARY_PERIODIC:
192: stag->entries = stag->n[0] * stag->entriesPerElement;
193: break;
194: default:
195: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
196: }
198: /* Starting element */
199: stag->start[0] = 0;
200: for (j = 0; j < stag->rank[0]; ++j) stag->start[0] += stag->l[0][j];
202: /* Local/ghosted size and starting element */
203: switch (stag->boundaryType[0]) {
204: case DM_BOUNDARY_NONE:
205: switch (stag->stencilType) {
206: case DMSTAG_STENCIL_NONE: /* Only dummy cells on the right */
207: stag->startGhost[0] = stag->start[0];
208: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
209: break;
210: case DMSTAG_STENCIL_STAR:
211: case DMSTAG_STENCIL_BOX:
212: stag->startGhost[0] = stag->firstRank[0] ? stag->start[0] : stag->start[0] - stag->stencilWidth;
213: stag->nGhost[0] = stag->n[0];
214: stag->nGhost[0] += stag->firstRank[0] ? 0 : stag->stencilWidth;
215: stag->nGhost[0] += stag->lastRank[0] ? 1 : stag->stencilWidth;
216: break;
217: default:
218: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
219: }
220: break;
221: case DM_BOUNDARY_GHOSTED:
222: switch (stag->stencilType) {
223: case DMSTAG_STENCIL_NONE:
224: stag->startGhost[0] = stag->start[0];
225: stag->nGhost[0] = stag->n[0] + (stag->lastRank[0] ? 1 : 0);
226: break;
227: case DMSTAG_STENCIL_STAR:
228: case DMSTAG_STENCIL_BOX:
229: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
230: stag->nGhost[0] = stag->n[0] + 2 * stag->stencilWidth + (stag->lastRank[0] && stag->stencilWidth == 0 ? 1 : 0);
231: break;
232: default:
233: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
234: }
235: break;
236: case DM_BOUNDARY_PERIODIC:
237: switch (stag->stencilType) {
238: case DMSTAG_STENCIL_NONE:
239: stag->startGhost[0] = stag->start[0];
240: stag->nGhost[0] = stag->n[0];
241: break;
242: case DMSTAG_STENCIL_STAR:
243: case DMSTAG_STENCIL_BOX:
244: stag->startGhost[0] = stag->start[0] - stag->stencilWidth; /* This value may be negative */
245: stag->nGhost[0] = stag->n[0] + 2 * stag->stencilWidth;
246: break;
247: default:
248: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
249: }
250: break;
251: default:
252: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
253: }
255: /* Total size of ghosted/local representation */
256: stag->entriesGhost = stag->nGhost[0] * stag->entriesPerElement;
258: /* Define neighbors */
259: PetscCall(PetscMalloc1(3, &stag->neighbors));
260: if (stag->firstRank[0]) {
261: switch (stag->boundaryType[0]) {
262: case DM_BOUNDARY_GHOSTED:
263: case DM_BOUNDARY_NONE:
264: stag->neighbors[0] = -1;
265: break;
266: case DM_BOUNDARY_PERIODIC:
267: stag->neighbors[0] = stag->nRanks[0] - 1;
268: break;
269: default:
270: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
271: }
272: } else {
273: stag->neighbors[0] = stag->rank[0] - 1;
274: }
275: stag->neighbors[1] = stag->rank[0];
276: if (stag->lastRank[0]) {
277: switch (stag->boundaryType[0]) {
278: case DM_BOUNDARY_GHOSTED:
279: case DM_BOUNDARY_NONE:
280: stag->neighbors[2] = -1;
281: break;
282: case DM_BOUNDARY_PERIODIC:
283: stag->neighbors[2] = 0;
284: break;
285: default:
286: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
287: }
288: } else {
289: stag->neighbors[2] = stag->rank[0] + 1;
290: }
292: PetscCheck(stag->n[0] >= stag->stencilWidth, PETSC_COMM_SELF, PETSC_ERR_SUP, "DMStag 1d setup does not support local sizes (%" PetscInt_FMT ") smaller than the elementwise stencil width (%" PetscInt_FMT ")", stag->n[0], stag->stencilWidth);
294: /* Create global->local VecScatter and ISLocalToGlobalMapping */
295: {
296: PetscInt *idxLocal, *idxGlobal, *idxGlobalAll;
297: PetscInt i, iLocal, d, entriesToTransferTotal, ghostOffsetStart, ghostOffsetEnd, nNonDummyGhost;
298: IS isLocal, isGlobal;
300: /* The offset on the right (may not be equal to the stencil width, as we
301: always have at least one ghost element, to account for the boundary
302: point, and may with ghosted boundaries), and the number of non-dummy ghost elements */
303: ghostOffsetStart = stag->start[0] - stag->startGhost[0];
304: ghostOffsetEnd = stag->startGhost[0] + stag->nGhost[0] - (stag->start[0] + stag->n[0]);
305: nNonDummyGhost = stag->nGhost[0] - (stag->lastRank[0] ? ghostOffsetEnd : 0) - (stag->firstRank[0] ? ghostOffsetStart : 0);
307: /* Compute the number of non-dummy entries in the local representation
308: This is equal to the number of non-dummy elements in the local (ghosted) representation,
309: plus some extra entries on the right boundary on the last rank*/
310: switch (stag->boundaryType[0]) {
311: case DM_BOUNDARY_GHOSTED:
312: case DM_BOUNDARY_NONE:
313: entriesToTransferTotal = nNonDummyGhost * stag->entriesPerElement + (stag->lastRank[0] ? stag->dof[0] : 0);
314: break;
315: case DM_BOUNDARY_PERIODIC:
316: entriesToTransferTotal = stag->entriesGhost; /* No dummy points */
317: break;
318: default:
319: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
320: }
322: PetscCall(PetscMalloc1(entriesToTransferTotal, &idxLocal));
323: PetscCall(PetscMalloc1(entriesToTransferTotal, &idxGlobal));
324: PetscCall(PetscMalloc1(stag->entriesGhost, &idxGlobalAll));
325: if (stag->boundaryType[0] == DM_BOUNDARY_NONE) {
326: PetscInt count = 0, countAll = 0;
327: /* Left ghost points and native points */
328: for (i = stag->startGhost[0], iLocal = 0; iLocal < nNonDummyGhost; ++i, ++iLocal) {
329: for (d = 0; d < stag->entriesPerElement; ++d, ++count, ++countAll) {
330: idxLocal[count] = iLocal * stag->entriesPerElement + d;
331: idxGlobal[count] = i * stag->entriesPerElement + d;
332: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
333: }
334: }
335: /* Ghost points on the right
336: Special case for last (partial dummy) element on the last rank */
337: if (stag->lastRank[0]) {
338: i = stag->N[0];
339: iLocal = (stag->nGhost[0] - ghostOffsetEnd);
340: /* Only vertex (0-cell) dofs in global representation */
341: for (d = 0; d < stag->dof[0]; ++d, ++count, ++countAll) {
342: idxGlobal[count] = i * stag->entriesPerElement + d;
343: idxLocal[count] = iLocal * stag->entriesPerElement + d;
344: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
345: }
346: for (d = stag->dof[0]; d < stag->entriesPerElement; ++d, ++countAll) { /* Additional dummy entries */
347: idxGlobalAll[countAll] = -1;
348: }
349: }
350: } else if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC) {
351: PetscInt count = 0, iLocal = 0; /* No dummy points, so idxGlobal and idxGlobalAll are identical */
352: const PetscInt iMin = stag->firstRank[0] ? stag->start[0] : stag->startGhost[0];
353: const PetscInt iMax = stag->lastRank[0] ? stag->startGhost[0] + stag->nGhost[0] - stag->stencilWidth : stag->startGhost[0] + stag->nGhost[0];
354: /* Ghost points on the left */
355: if (stag->firstRank[0]) {
356: for (i = stag->N[0] - stag->stencilWidth; iLocal < stag->stencilWidth; ++i, ++iLocal) {
357: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
358: idxGlobal[count] = i * stag->entriesPerElement + d;
359: idxLocal[count] = iLocal * stag->entriesPerElement + d;
360: idxGlobalAll[count] = idxGlobal[count];
361: }
362: }
363: }
364: /* Native points */
365: for (i = iMin; i < iMax; ++i, ++iLocal) {
366: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
367: idxGlobal[count] = i * stag->entriesPerElement + d;
368: idxLocal[count] = iLocal * stag->entriesPerElement + d;
369: idxGlobalAll[count] = idxGlobal[count];
370: }
371: }
372: /* Ghost points on the right */
373: if (stag->lastRank[0]) {
374: for (i = 0; iLocal < stag->nGhost[0]; ++i, ++iLocal) {
375: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
376: idxGlobal[count] = i * stag->entriesPerElement + d;
377: idxLocal[count] = iLocal * stag->entriesPerElement + d;
378: idxGlobalAll[count] = idxGlobal[count];
379: }
380: }
381: }
382: } else if (stag->boundaryType[0] == DM_BOUNDARY_GHOSTED) {
383: PetscInt count = 0, countAll = 0;
384: /* Dummy elements on the left, on the first rank */
385: if (stag->firstRank[0]) {
386: for (iLocal = 0; iLocal < ghostOffsetStart; ++iLocal) {
387: /* Complete elements full of dummy entries */
388: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
389: }
390: i = 0; /* nonDummy entries start with global entry 0 */
391: } else {
392: /* nonDummy entries start as usual */
393: i = stag->startGhost[0];
394: iLocal = 0;
395: }
397: /* non-Dummy entries */
398: {
399: PetscInt iLocalNonDummyMax = stag->firstRank[0] ? nNonDummyGhost + ghostOffsetStart : nNonDummyGhost;
400: for (; iLocal < iLocalNonDummyMax; ++i, ++iLocal) {
401: for (d = 0; d < stag->entriesPerElement; ++d, ++count, ++countAll) {
402: idxLocal[count] = iLocal * stag->entriesPerElement + d;
403: idxGlobal[count] = i * stag->entriesPerElement + d;
404: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
405: }
406: }
407: }
409: /* (partial) dummy elements on the right, on the last rank */
410: if (stag->lastRank[0]) {
411: /* First one is partial dummy */
412: i = stag->N[0];
413: iLocal = (stag->nGhost[0] - ghostOffsetEnd);
414: for (d = 0; d < stag->dof[0]; ++d, ++count, ++countAll) { /* Only vertex (0-cell) dofs in global representation */
415: idxLocal[count] = iLocal * stag->entriesPerElement + d;
416: idxGlobal[count] = i * stag->entriesPerElement + d;
417: idxGlobalAll[countAll] = i * stag->entriesPerElement + d;
418: }
419: for (d = stag->dof[0]; d < stag->entriesPerElement; ++d, ++countAll) { /* Additional dummy entries */
420: idxGlobalAll[countAll] = -1;
421: }
422: for (iLocal = stag->nGhost[0] - ghostOffsetEnd + 1; iLocal < stag->nGhost[0]; ++iLocal) {
423: /* Additional dummy elements */
424: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
425: }
426: }
427: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported x boundary type %s", DMBoundaryTypes[stag->boundaryType[0]]);
429: /* Create Local IS (transferring pointer ownership) */
430: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxLocal, PETSC_OWN_POINTER, &isLocal));
432: /* Create Global IS (transferring pointer ownership) */
433: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxGlobal, PETSC_OWN_POINTER, &isGlobal));
435: /* Create stag->gtol, which doesn't include dummy entries */
436: {
437: Vec local, global;
438: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global));
439: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, stag->entriesPerElement, stag->entriesGhost, NULL, &local));
440: PetscCall(VecScatterCreate(global, isGlobal, local, isLocal, &stag->gtol));
441: PetscCall(VecDestroy(&global));
442: PetscCall(VecDestroy(&local));
443: }
445: /* In special cases, create a dedicated injective local-to-global map */
446: if (stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) PetscCall(DMStagPopulateLocalToGlobalInjective(dm));
448: /* Destroy ISs */
449: PetscCall(ISDestroy(&isLocal));
450: PetscCall(ISDestroy(&isGlobal));
452: /* Create local-to-global map (transferring pointer ownership) */
453: PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stag->entriesGhost, idxGlobalAll, PETSC_OWN_POINTER, &dm->ltogmap));
454: }
456: /* Precompute location offsets */
457: PetscCall(DMStagComputeLocationOffsets_1d(dm));
459: /* View from Options */
460: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: static PetscErrorCode DMStagComputeLocationOffsets_1d(DM dm)
466: {
467: DM_Stag *const stag = (DM_Stag *)dm->data;
468: const PetscInt epe = stag->entriesPerElement;
470: PetscFunctionBegin;
471: PetscCall(PetscMalloc1(DMSTAG_NUMBER_LOCATIONS, &stag->locationOffsets));
472: stag->locationOffsets[DMSTAG_LEFT] = 0;
473: stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[0];
474: stag->locationOffsets[DMSTAG_RIGHT] = stag->locationOffsets[DMSTAG_LEFT] + epe;
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_1d(DM dm)
479: {
480: DM_Stag *const stag = (DM_Stag *)dm->data;
481: PetscInt *idxLocal, *idxGlobal;
482: PetscInt i, iLocal, d, count;
483: IS isLocal, isGlobal;
485: PetscFunctionBegin;
486: PetscCall(PetscMalloc1(stag->entries, &idxLocal));
487: PetscCall(PetscMalloc1(stag->entries, &idxGlobal));
488: count = 0;
489: iLocal = stag->start[0] - stag->startGhost[0];
490: for (i = stag->start[0]; i < stag->start[0] + stag->n[0]; ++i, ++iLocal) {
491: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
492: idxGlobal[count] = i * stag->entriesPerElement + d;
493: idxLocal[count] = iLocal * stag->entriesPerElement + d;
494: }
495: }
496: if (stag->lastRank[0] && stag->boundaryType[0] != DM_BOUNDARY_PERIODIC) {
497: i = stag->start[0] + stag->n[0];
498: iLocal = stag->start[0] - stag->startGhost[0] + stag->n[0];
499: for (d = 0; d < stag->dof[0]; ++d, ++count) {
500: idxGlobal[count] = i * stag->entriesPerElement + d;
501: idxLocal[count] = iLocal * stag->entriesPerElement + d;
502: }
503: }
504: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxLocal, PETSC_OWN_POINTER, &isLocal));
505: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxGlobal, PETSC_OWN_POINTER, &isGlobal));
506: {
507: Vec local, global;
508: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global));
509: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, stag->entriesPerElement, stag->entriesGhost, NULL, &local));
510: PetscCall(VecScatterCreate(local, isLocal, global, isGlobal, &stag->ltog_injective));
511: PetscCall(VecDestroy(&global));
512: PetscCall(VecDestroy(&local));
513: }
514: PetscCall(ISDestroy(&isLocal));
515: PetscCall(ISDestroy(&isGlobal));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: PETSC_INTERN PetscErrorCode DMCreateMatrix_Stag_1D_AIJ_Assemble(DM dm, Mat A)
520: {
521: DMStagStencilType stencil_type;
522: PetscInt dof[2], start, n, n_extra, stencil_width, N, epe;
523: DMBoundaryType boundary_type_x;
525: PetscFunctionBegin;
526: PetscCall(DMStagGetDOF(dm, &dof[0], &dof[1], NULL, NULL));
527: PetscCall(DMStagGetStencilType(dm, &stencil_type));
528: PetscCall(DMStagGetStencilWidth(dm, &stencil_width));
529: PetscCall(DMStagGetCorners(dm, &start, NULL, NULL, &n, NULL, NULL, &n_extra, NULL, NULL));
530: PetscCall(DMStagGetGlobalSizes(dm, &N, NULL, NULL));
531: PetscCall(DMStagGetEntriesPerElement(dm, &epe));
532: PetscCall(DMStagGetBoundaryTypes(dm, &boundary_type_x, NULL, NULL));
533: if (stencil_type == DMSTAG_STENCIL_NONE) {
534: /* Couple all DOF at each location to each other */
535: DMStagStencil *row_vertex, *row_element;
537: PetscCall(PetscMalloc1(dof[0], &row_vertex));
538: for (PetscInt c = 0; c < dof[0]; ++c) {
539: row_vertex[c].loc = DMSTAG_LEFT;
540: row_vertex[c].c = c;
541: }
543: PetscCall(PetscMalloc1(dof[1], &row_element));
544: for (PetscInt c = 0; c < dof[1]; ++c) {
545: row_element[c].loc = DMSTAG_ELEMENT;
546: row_element[c].c = c;
547: }
549: for (PetscInt e = start; e < start + n + n_extra; ++e) {
550: {
551: for (PetscInt c = 0; c < dof[0]; ++c) row_vertex[c].i = e;
552: PetscCall(DMStagMatSetValuesStencil(dm, A, dof[0], row_vertex, dof[0], row_vertex, NULL, INSERT_VALUES));
553: }
554: if (e < N) {
555: for (PetscInt c = 0; c < dof[1]; ++c) row_element[c].i = e;
556: PetscCall(DMStagMatSetValuesStencil(dm, A, dof[1], row_element, dof[1], row_element, NULL, INSERT_VALUES));
557: }
558: }
559: PetscCall(PetscFree(row_vertex));
560: PetscCall(PetscFree(row_element));
561: } else if (stencil_type == DMSTAG_STENCIL_STAR || stencil_type == DMSTAG_STENCIL_BOX) {
562: DMStagStencil *col, *row;
564: PetscCall(PetscMalloc1(epe, &row));
565: {
566: PetscInt nrows = 0;
567: for (PetscInt c = 0; c < dof[0]; ++c) {
568: row[nrows].c = c;
569: row[nrows].loc = DMSTAG_LEFT;
570: ++nrows;
571: }
572: for (PetscInt c = 0; c < dof[1]; ++c) {
573: row[nrows].c = c;
574: row[nrows].loc = DMSTAG_ELEMENT;
575: ++nrows;
576: }
577: }
578: PetscCall(PetscMalloc1(epe, &col));
579: {
580: PetscInt ncols = 0;
581: for (PetscInt c = 0; c < dof[0]; ++c) {
582: col[ncols].c = c;
583: col[ncols].loc = DMSTAG_LEFT;
584: ++ncols;
585: }
586: for (PetscInt c = 0; c < dof[1]; ++c) {
587: col[ncols].c = c;
588: col[ncols].loc = DMSTAG_ELEMENT;
589: ++ncols;
590: }
591: }
592: for (PetscInt e = start; e < start + n + n_extra; ++e) {
593: for (PetscInt i = 0; i < epe; ++i) row[i].i = e;
594: for (PetscInt offset = -stencil_width; offset <= stencil_width; ++offset) {
595: const PetscInt e_offset = e + offset;
597: /* Only set values corresponding to elements which can have non-dummy entries,
598: meaning those that map to unknowns in the global representation. In the periodic
599: case, this is the entire stencil, but in all other cases, only includes a single
600: "extra" element which is partially outside the physical domain (those points in the
601: global representation */
602: if (boundary_type_x == DM_BOUNDARY_PERIODIC || (e_offset < N + 1 && e_offset >= 0)) {
603: for (PetscInt i = 0; i < epe; ++i) col[i].i = e_offset;
604: PetscCall(DMStagMatSetValuesStencil(dm, A, epe, row, epe, col, NULL, INSERT_VALUES));
605: }
606: }
607: }
608: PetscCall(PetscFree(row));
609: PetscCall(PetscFree(col));
610: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported stencil type %s", DMStagStencilTypes[stencil_type]);
611: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
612: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }