Actual source code: stag2d.c
1: /* Functions specific to the 2-dimensional implementation of DMStag */
2: #include <petsc/private/dmstagimpl.h>
4: /*@C
5: DMStagCreate2d - Create an object to manage data living on the elements, faces, and vertices of a parallelized regular 2D grid.
7: Collective
9: Input Parameters:
10: + comm - MPI communicator
11: . bndx - x boundary type, `DM_BOUNDARY_NONE`, `DM_BOUNDARY_PERIODIC`, or
12: `DM_BOUNDARY_GHOSTED`
13: . bndy - y boundary type, `DM_BOUNDARY_NONE`, `DM_BOUNDARY_PERIODIC`, or `DM_BOUNDARY_GHOSTED`
14: . M - global number of elements in x direction
15: . N - global number of elements in y direction
16: . m - number of ranks in the x direction (may be `PETSC_DECIDE`)
17: . n - number of ranks in the y direction (may be `PETSC_DECIDE`)
18: . dof0 - number of degrees of freedom per vertex/0-cell
19: . dof1 - number of degrees of freedom per face/1-cell
20: . dof2 - number of degrees of freedom per element/2-cell
21: . stencilType - ghost/halo region type: `DMSTAG_STENCIL_NONE`, `DMSTAG_STENCIL_BOX`, or `DMSTAG_STENCIL_STAR`
22: . stencilWidth - width, in elements, of halo/ghost region
23: . lx - array of local x element counts, of length equal to `m`, summing to `M`
24: - ly - array of local y element counts, of length equal to `n`, summing to `N`
26: Output Parameter:
27: . dm - the new `DMSTAG` object
29: Options Database Keys:
30: + -dm_view - calls `DMViewFromOptions()` at the conclusion of `DMSetUp()`
31: . -stag_grid_x <nx> - number of elements in the x direction
32: . -stag_grid_y <ny> - number of elements in the y direction
33: . -stag_ranks_x <rx> - number of ranks in the x direction
34: . -stag_ranks_y <ry> - number of ranks in the y direction
35: . -stag_ghost_stencil_width - width of ghost region, in elements
36: . -stag_boundary_type_x <none,ghosted,periodic> - `DMBoundaryType` value
37: - -stag_boundary_type_y <none,ghosted,periodic> - `DMBoundaryType` value
39: Level: beginner
41: Notes:
42: You must call `DMSetUp()` after this call, before using the `DM`.
43: If you wish to use the options database (see the keys above) to change values in the `DMSTAG`, you must call
44: `DMSetFromOptions()` after this function but before `DMSetUp()`.
46: .seealso: [](ch_stag), `DMSTAG`, `DMStagCreate1d()`, `DMStagCreate3d()`, `DMDestroy()`, `DMView()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMLocalToGlobalBegin()`, `DMDACreate2d()`
47: @*/
48: PETSC_EXTERN PetscErrorCode DMStagCreate2d(MPI_Comm comm, DMBoundaryType bndx, DMBoundaryType bndy, PetscInt M, PetscInt N, PetscInt m, PetscInt n, PetscInt dof0, PetscInt dof1, PetscInt dof2, DMStagStencilType stencilType, PetscInt stencilWidth, const PetscInt lx[], const PetscInt ly[], DM *dm)
49: {
50: PetscFunctionBegin;
51: PetscCall(DMCreate(comm, dm));
52: PetscCall(DMSetDimension(*dm, 2));
53: PetscCall(DMStagInitialize(bndx, bndy, DM_BOUNDARY_NONE, M, N, 0, m, n, 0, dof0, dof1, dof2, 0, stencilType, stencilWidth, lx, ly, NULL, *dm));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: PETSC_INTERN PetscErrorCode DMStagRestrictSimple_2d(DM dmf, Vec xf_local, DM dmc, Vec xc_local)
58: {
59: PetscInt Mf, Nf, Mc, Nc, factorx, factory, dof[3];
60: PetscInt xc, yc, mc, nc, nExtraxc, nExtrayc, i, j, d;
61: PetscInt idownleftf, ileftf, idownf, ielemf, idownleftc, ileftc, idownc, ielemc;
62: const PetscScalar ***arrf;
63: PetscScalar ***arrc;
65: PetscFunctionBegin;
66: PetscCall(DMStagGetGlobalSizes(dmf, &Mf, &Nf, NULL));
67: PetscCall(DMStagGetGlobalSizes(dmc, &Mc, &Nc, NULL));
68: factorx = Mf / Mc;
69: factory = Nf / Nc;
70: PetscCall(DMStagGetDOF(dmc, &dof[0], &dof[1], &dof[2], NULL));
72: PetscCall(DMStagGetCorners(dmc, &xc, &yc, NULL, &mc, &nc, NULL, &nExtraxc, &nExtrayc, NULL));
73: PetscCall(VecZeroEntries(xc_local));
74: PetscCall(DMStagVecGetArray(dmf, xf_local, &arrf));
75: PetscCall(DMStagVecGetArray(dmc, xc_local, &arrc));
76: PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_DOWN_LEFT, 0, &idownleftf));
77: PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_LEFT, 0, &ileftf));
78: PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_DOWN, 0, &idownf));
79: PetscCall(DMStagGetLocationSlot(dmf, DMSTAG_ELEMENT, 0, &ielemf));
80: PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_DOWN_LEFT, 0, &idownleftc));
81: PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_LEFT, 0, &ileftc));
82: PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_DOWN, 0, &idownc));
83: PetscCall(DMStagGetLocationSlot(dmc, DMSTAG_ELEMENT, 0, &ielemc));
85: for (d = 0; d < dof[0]; ++d)
86: for (j = yc; j < yc + nc + nExtrayc; ++j)
87: for (i = xc; i < xc + mc + nExtraxc; ++i) {
88: const PetscInt ii = factorx * i, jj = factory * j;
90: arrc[j][i][idownleftc + d] = arrf[jj][ii][idownleftf + d];
91: }
93: for (d = 0; d < dof[1]; ++d)
94: for (j = yc; j < yc + nc; ++j)
95: for (i = xc; i < xc + mc + nExtraxc; ++i) {
96: const PetscInt ii = factorx * i, jj = factory * j + factory / 2;
98: if (factory % 2 == 0) arrc[j][i][ileftc + d] = 0.5 * (arrf[jj - 1][ii][ileftf + d] + arrf[jj][ii][ileftf + d]);
99: else arrc[j][i][ileftc + d] = arrf[jj][ii][ileftf + d];
100: }
102: for (d = 0; d < dof[1]; ++d)
103: for (j = yc; j < yc + nc + nExtrayc; ++j)
104: for (i = xc; i < xc + mc; ++i) {
105: const PetscInt ii = factorx * i + factorx / 2, jj = factory * j;
107: if (factorx % 2 == 0) arrc[j][i][idownc + d] = 0.5 * (arrf[jj][ii - 1][idownf + d] + arrf[jj][ii][idownf + d]);
108: else arrc[j][i][idownc + d] = arrf[jj][ii][idownf + d];
109: }
111: for (d = 0; d < dof[2]; ++d)
112: for (j = yc; j < yc + nc; ++j)
113: for (i = xc; i < xc + mc; ++i) {
114: const PetscInt ii = factorx * i + factorx / 2, jj = factory * j + factory / 2;
116: if (factorx % 2 == 0 && factory % 2 == 0) arrc[j][i][ielemc + d] = 0.25 * (arrf[jj - 1][ii - 1][ielemf + d] + arrf[jj][ii - 1][ielemf + d] + arrf[jj - 1][ii][ielemf + d] + arrf[jj][ii][ielemf + d]);
117: else if (factorx % 2 == 0) arrc[j][i][ielemc + d] = 0.5 * (arrf[jj - 1][ii - 1][ielemf + d] + arrf[jj][ii - 1][ielemf + d]);
118: else if (factory % 2 == 0) arrc[j][i][ielemc + d] = 0.5 * (arrf[jj - 1][ii - 1][ielemf + d] + arrf[jj - 1][ii][ielemf + d]);
119: else arrc[j][i][ielemc + d] = arrf[jj][ii][ielemf + d];
120: }
122: PetscCall(DMStagVecRestoreArray(dmf, xf_local, &arrf));
123: PetscCall(DMStagVecRestoreArray(dmc, xc_local, &arrc));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_2d(DM dm, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax)
128: {
129: DM_Stag *stagCoord;
130: DM dmCoord;
131: Vec coordLocal;
132: PetscReal h[2], min[2];
133: PetscScalar ***arr;
134: PetscInt ind[2], start_ghost[2], n_ghost[2], s, c;
135: PetscInt idownleft, idown, ileft, ielement;
137: PetscFunctionBegin;
138: PetscCall(DMGetCoordinateDM(dm, &dmCoord));
139: stagCoord = (DM_Stag *)dmCoord->data;
140: for (s = 0; s < 3; ++s) {
141: PetscCheck(stagCoord->dof[s] == 0 || stagCoord->dof[s] == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Coordinate DM in 2 dimensions must have 0 or 2 dof on each stratum, but stratum %" PetscInt_FMT " has %" PetscInt_FMT " dof", s,
142: stagCoord->dof[s]);
143: }
144: PetscCall(DMCreateLocalVector(dmCoord, &coordLocal));
146: PetscCall(DMStagVecGetArray(dmCoord, coordLocal, &arr));
147: if (stagCoord->dof[0]) PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_DOWN_LEFT, 0, &idownleft));
148: if (stagCoord->dof[1]) {
149: PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_DOWN, 0, &idown));
150: PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_LEFT, 0, &ileft));
151: }
152: if (stagCoord->dof[2]) PetscCall(DMStagGetLocationSlot(dmCoord, DMSTAG_ELEMENT, 0, &ielement));
153: PetscCall(DMStagGetGhostCorners(dmCoord, &start_ghost[0], &start_ghost[1], NULL, &n_ghost[0], &n_ghost[1], NULL));
155: min[0] = xmin;
156: min[1] = ymin;
157: h[0] = (xmax - xmin) / stagCoord->N[0];
158: h[1] = (ymax - ymin) / stagCoord->N[1];
160: for (ind[1] = start_ghost[1]; ind[1] < start_ghost[1] + n_ghost[1]; ++ind[1]) {
161: for (ind[0] = start_ghost[0]; ind[0] < start_ghost[0] + n_ghost[0]; ++ind[0]) {
162: if (stagCoord->dof[0]) {
163: const PetscReal offs[2] = {0.0, 0.0};
164: for (c = 0; c < 2; ++c) arr[ind[1]][ind[0]][idownleft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
165: }
166: if (stagCoord->dof[1]) {
167: const PetscReal offs[2] = {0.5, 0.0};
168: for (c = 0; c < 2; ++c) arr[ind[1]][ind[0]][idown + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
169: }
170: if (stagCoord->dof[1]) {
171: const PetscReal offs[2] = {0.0, 0.5};
172: for (c = 0; c < 2; ++c) arr[ind[1]][ind[0]][ileft + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
173: }
174: if (stagCoord->dof[2]) {
175: const PetscReal offs[2] = {0.5, 0.5};
176: for (c = 0; c < 2; ++c) arr[ind[1]][ind[0]][ielement + c] = min[c] + ((PetscReal)ind[c] + offs[c]) * h[c];
177: }
178: }
179: }
180: PetscCall(DMStagVecRestoreArray(dmCoord, coordLocal, &arr));
181: PetscCall(DMSetCoordinatesLocal(dm, coordLocal));
182: PetscCall(VecDestroy(&coordLocal));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /* Helper functions used in DMSetUp_Stag() */
187: static PetscErrorCode DMStagSetUpBuildRankGrid_2d(DM);
188: static PetscErrorCode DMStagSetUpBuildNeighbors_2d(DM);
189: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_2d(DM, PetscInt **);
190: static PetscErrorCode DMStagComputeLocationOffsets_2d(DM);
192: PETSC_INTERN PetscErrorCode DMSetUp_Stag_2d(DM dm)
193: {
194: DM_Stag *const stag = (DM_Stag *)dm->data;
195: PetscMPIInt size, rank;
196: PetscInt i, j, d, entriesPerElementRowGhost, entriesPerCorner, entriesPerFace, entriesPerElementRow;
197: MPI_Comm comm;
198: PetscInt *globalOffsets;
199: PetscBool star, dummyStart[2], dummyEnd[2];
200: const PetscInt dim = 2;
202: PetscFunctionBegin;
203: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
204: PetscCallMPI(MPI_Comm_size(comm, &size));
205: PetscCallMPI(MPI_Comm_rank(comm, &rank));
207: /* Rank grid sizes (populates stag->nRanks) */
208: PetscCall(DMStagSetUpBuildRankGrid_2d(dm));
210: /* Determine location of rank in grid (these get extra boundary points on the last element)
211: Order is x-fast, as usual */
212: stag->rank[0] = rank % stag->nRanks[0];
213: stag->rank[1] = rank / stag->nRanks[0];
214: for (i = 0; i < dim; ++i) {
215: stag->firstRank[i] = PetscNot(stag->rank[i]);
216: stag->lastRank[i] = (PetscBool)(stag->rank[i] == stag->nRanks[i] - 1);
217: }
219: /* Determine Locally owned region
221: Divide equally, giving lower ranks in each dimension and extra element if needbe.
223: Note that this uses O(P) storage. If this ever becomes an issue, this could
224: be refactored to not keep this data around. */
225: for (i = 0; i < dim; ++i) {
226: if (!stag->l[i]) {
227: const PetscInt Ni = stag->N[i], nRanksi = stag->nRanks[i];
228: PetscCall(PetscMalloc1(stag->nRanks[i], &stag->l[i]));
229: for (j = 0; j < stag->nRanks[i]; ++j) stag->l[i][j] = Ni / nRanksi + ((Ni % nRanksi) > j);
230: }
231: }
233: /* Retrieve local size in stag->n */
234: for (i = 0; i < dim; ++i) stag->n[i] = stag->l[i][stag->rank[i]];
235: if (PetscDefined(USE_DEBUG)) {
236: for (i = 0; i < dim; ++i) {
237: PetscInt Ncheck, j;
238: Ncheck = 0;
239: for (j = 0; j < stag->nRanks[i]; ++j) Ncheck += stag->l[i][j];
240: PetscCheck(Ncheck == stag->N[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local sizes in dimension %" PetscInt_FMT " don't add up. %" PetscInt_FMT " != %" PetscInt_FMT, i, Ncheck, stag->N[i]);
241: }
242: }
244: /* Compute starting elements */
245: for (i = 0; i < dim; ++i) {
246: stag->start[i] = 0;
247: for (j = 0; j < stag->rank[i]; ++j) stag->start[i] += stag->l[i][j];
248: }
250: /* Determine ranks of neighbors, using DMDA's convention
252: n6 n7 n8
253: n3 n5
254: n0 n1 n2 */
255: PetscCall(DMStagSetUpBuildNeighbors_2d(dm));
257: /* Determine whether the ghost region includes dummies or not. This is currently
258: equivalent to having a non-periodic boundary. If not, then
259: ghostOffset{Start,End}[d] elements correspond to elements on the neighbor.
260: If true, then
261: - at the start, there are ghostOffsetStart[d] ghost elements
262: - at the end, there is a layer of extra "physical" points inside a layer of
263: ghostOffsetEnd[d] ghost elements
264: Note that this computation should be updated if any boundary types besides
265: NONE, GHOSTED, and PERIODIC are supported. */
266: for (d = 0; d < 2; ++d) dummyStart[d] = (PetscBool)(stag->firstRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
267: for (d = 0; d < 2; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
269: /* Define useful sizes */
270: stag->entriesPerElement = stag->dof[0] + 2 * stag->dof[1] + stag->dof[2];
271: entriesPerFace = stag->dof[0] + stag->dof[1];
272: entriesPerCorner = stag->dof[0];
273: entriesPerElementRow = stag->n[0] * stag->entriesPerElement + (dummyEnd[0] ? entriesPerFace : 0);
274: stag->entries = stag->n[1] * entriesPerElementRow + (dummyEnd[1] ? stag->n[0] * entriesPerFace : 0) + (dummyEnd[0] && dummyEnd[1] ? entriesPerCorner : 0);
276: /* Compute offsets for each rank into global vectors
277: This again requires O(P) storage, which could be replaced with some global
278: communication. */
279: PetscCall(DMStagSetUpBuildGlobalOffsets_2d(dm, &globalOffsets));
281: for (d = 0; d < dim; ++d)
282: PetscCheck(stag->boundaryType[d] == DM_BOUNDARY_NONE || stag->boundaryType[d] == DM_BOUNDARY_PERIODIC || stag->boundaryType[d] == DM_BOUNDARY_GHOSTED, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported boundary type");
284: /* Define ghosted/local sizes */
285: if (stag->stencilType != DMSTAG_STENCIL_NONE && (stag->n[0] < stag->stencilWidth || stag->n[1] < stag->stencilWidth)) {
286: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "DMStag 2d setup does not support local sizes (%" PetscInt_FMT " x %" PetscInt_FMT ") smaller than the elementwise stencil width (%" PetscInt_FMT ")", stag->n[0], stag->n[1], stag->stencilWidth);
287: }
288: for (d = 0; d < dim; ++d) {
289: switch (stag->boundaryType[d]) {
290: case DM_BOUNDARY_NONE:
291: /* Note: for a elements-only DMStag, the extra elements on the faces aren't necessary but we include them anyway */
292: switch (stag->stencilType) {
293: case DMSTAG_STENCIL_NONE: /* only the extra one on the right/top faces */
294: stag->nGhost[d] = stag->n[d];
295: stag->startGhost[d] = stag->start[d];
296: if (stag->lastRank[d]) stag->nGhost[d] += 1;
297: break;
298: case DMSTAG_STENCIL_STAR: /* allocate the corners but don't use them */
299: case DMSTAG_STENCIL_BOX:
300: stag->nGhost[d] = stag->n[d];
301: stag->startGhost[d] = stag->start[d];
302: if (!stag->firstRank[d]) {
303: stag->nGhost[d] += stag->stencilWidth; /* add interior ghost elements */
304: stag->startGhost[d] -= stag->stencilWidth;
305: }
306: if (!stag->lastRank[d]) {
307: stag->nGhost[d] += stag->stencilWidth; /* add interior ghost elements */
308: } else {
309: stag->nGhost[d] += 1; /* one element on the boundary to complete blocking */
310: }
311: break;
312: default:
313: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
314: }
315: break;
316: case DM_BOUNDARY_GHOSTED:
317: switch (stag->stencilType) {
318: case DMSTAG_STENCIL_NONE:
319: stag->startGhost[d] = stag->start[d];
320: stag->nGhost[d] = stag->n[d] + (stag->lastRank[d] ? 1 : 0);
321: break;
322: case DMSTAG_STENCIL_STAR:
323: case DMSTAG_STENCIL_BOX:
324: stag->startGhost[d] = stag->start[d] - stag->stencilWidth; /* This value may be negative */
325: stag->nGhost[d] = stag->n[d] + 2 * stag->stencilWidth + (stag->lastRank[d] && stag->stencilWidth == 0 ? 1 : 0);
326: break;
327: default:
328: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
329: }
330: break;
331: case DM_BOUNDARY_PERIODIC:
332: switch (stag->stencilType) {
333: case DMSTAG_STENCIL_NONE: /* only the extra one on the right/top faces */
334: stag->nGhost[d] = stag->n[d];
335: stag->startGhost[d] = stag->start[d];
336: break;
337: case DMSTAG_STENCIL_STAR:
338: case DMSTAG_STENCIL_BOX:
339: stag->nGhost[d] = stag->n[d] + 2 * stag->stencilWidth;
340: stag->startGhost[d] = stag->start[d] - stag->stencilWidth;
341: break;
342: default:
343: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unrecognized ghost stencil type %d", stag->stencilType);
344: }
345: break;
346: default:
347: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported boundary type in dimension %" PetscInt_FMT, d);
348: }
349: }
350: stag->entriesGhost = stag->nGhost[0] * stag->nGhost[1] * stag->entriesPerElement;
351: entriesPerElementRowGhost = stag->nGhost[0] * stag->entriesPerElement;
353: /* Create global-->local VecScatter and local->global ISLocalToGlobalMapping
355: We iterate over all local points twice. First, we iterate over each neighbor, populating
356: 1. idxLocal[] : the subset of points, in local numbering ("S" from 0 on all points including ghosts), which correspond to global points. That is, the set of all non-dummy points in the ghosted representation
357: 2. idxGlobal[]: the corresponding global points, in global numbering (Nested "S"s - ranks then non-ghost points in each rank)
359: Next, we iterate over all points in the local ordering, populating
360: 3. idxGlobalAll[] : entry i is the global point corresponding to local point i, or -1 if local point i is a dummy.
362: Note further here that the local/ghosted vectors:
363: - Are always an integral number of elements-worth of points, in all directions.
364: - Contain three flavors of points:
365: 1. Points which "live here" in the global representation
366: 2. Ghost points which correspond to points on other ranks in the global representation
367: 3. Ghost points, which we call "dummy points," which do not correspond to any point in the global representation
369: Dummy ghost points arise in at least three ways:
370: 1. As padding for the right, top, and front physical boundaries, to complete partial elements
371: 2. As unused space in the "corners" on interior ranks when using a star stencil
372: 3. As additional work space on all physical boundaries, when DM_BOUNDARY_GHOSTED is used
374: Note that, because of the boundary dummies,
375: with a stencil width of zero, on 1 rank, local and global vectors
376: are still different!
378: We assume that the size on each rank is greater than or equal to the
379: stencil width.
380: */
382: /* Check stencil type */
383: PetscCheck(stag->stencilType == DMSTAG_STENCIL_NONE || stag->stencilType == DMSTAG_STENCIL_BOX || stag->stencilType == DMSTAG_STENCIL_STAR, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported stencil type %s", DMStagStencilTypes[stag->stencilType]);
384: star = (PetscBool)(stag->stencilType == DMSTAG_STENCIL_STAR || stag->stencilType == DMSTAG_STENCIL_NONE);
386: {
387: PetscInt *idxLocal, *idxGlobal, *idxGlobalAll;
388: PetscInt count, countAll, entriesToTransferTotal, i, j, d, ghostOffsetStart[2], ghostOffsetEnd[2];
389: IS isLocal, isGlobal;
390: PetscInt jghost, ighost;
391: PetscInt nNeighbors[9][2];
392: PetscBool nextToDummyEnd[2];
394: /* Compute numbers of elements on each neighbor */
395: for (i = 0; i < 9; ++i) {
396: const PetscInt neighborRank = stag->neighbors[i];
397: if (neighborRank >= 0) { /* note we copy the values for our own rank (neighbor 4) */
398: nNeighbors[i][0] = stag->l[0][neighborRank % stag->nRanks[0]];
399: nNeighbors[i][1] = stag->l[1][neighborRank / stag->nRanks[0]];
400: } else {
401: nNeighbors[i][0] = 0;
402: nNeighbors[i][1] = 0;
403: }
404: }
406: /* These offsets should always be non-negative, and describe how many
407: ghost elements exist at each boundary. These are not always equal to the stencil width,
408: because we may have different numbers of ghost elements at the boundaries. In particular,
409: we always have at least one ghost (dummy) element at the right/top/front. */
410: for (d = 0; d < 2; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
411: for (d = 0; d < 2; ++d) ghostOffsetEnd[d] = stag->startGhost[d] + stag->nGhost[d] - (stag->start[d] + stag->n[d]);
413: /* Compute whether the next rank has an extra point (only used in x direction) */
414: for (d = 0; d < 2; ++d) nextToDummyEnd[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC && stag->rank[d] == stag->nRanks[d] - 2);
416: /* Compute the number of local entries which correspond to any global entry */
417: {
418: PetscInt nNonDummyGhost[2];
419: for (d = 0; d < 2; ++d) nNonDummyGhost[d] = stag->nGhost[d] - (dummyStart[d] ? ghostOffsetStart[d] : 0) - (dummyEnd[d] ? ghostOffsetEnd[d] : 0);
420: if (star) {
421: entriesToTransferTotal = (nNonDummyGhost[0] * stag->n[1] + stag->n[0] * nNonDummyGhost[1] - stag->n[0] * stag->n[1]) * stag->entriesPerElement + (dummyEnd[0] ? nNonDummyGhost[1] * entriesPerFace : 0) + (dummyEnd[1] ? nNonDummyGhost[0] * entriesPerFace : 0) + (dummyEnd[0] && dummyEnd[1] ? entriesPerCorner : 0);
422: } else {
423: entriesToTransferTotal = nNonDummyGhost[0] * nNonDummyGhost[1] * stag->entriesPerElement + (dummyEnd[0] ? nNonDummyGhost[1] * entriesPerFace : 0) + (dummyEnd[1] ? nNonDummyGhost[0] * entriesPerFace : 0) + (dummyEnd[0] && dummyEnd[1] ? entriesPerCorner : 0);
424: }
425: }
427: /* Allocate arrays to populate */
428: PetscCall(PetscMalloc1(entriesToTransferTotal, &idxLocal));
429: PetscCall(PetscMalloc1(entriesToTransferTotal, &idxGlobal));
431: /* Counts into idxLocal/idxGlobal */
432: count = 0;
434: /* Here and below, we work with (i,j) describing element numbers within a neighboring rank's global ordering,
435: to be offset by that rank's global offset,
436: and (ighost,jghost) referring to element numbers within this ranks local (ghosted) ordering */
438: /* Neighbor 0 (down left) */
439: if (!star && !dummyStart[0] && !dummyStart[1]) {
440: const PetscInt neighbor = 0;
441: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
442: const PetscInt *const nNeighbor = nNeighbors[neighbor];
443: const PetscInt entriesPerElementRowNeighbor = stag->entriesPerElement * nNeighbor[0];
444: for (jghost = 0; jghost < ghostOffsetStart[1]; ++jghost) {
445: const PetscInt j = nNeighbor[1] - ghostOffsetStart[1] + jghost;
446: for (ighost = 0; ighost < ghostOffsetStart[0]; ++ighost) {
447: const PetscInt i = nNeighbor[0] - ghostOffsetStart[0] + ighost;
448: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
449: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
450: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
451: }
452: }
453: }
454: }
456: /* Neighbor 1 (down) */
457: if (!dummyStart[1]) {
458: /* We may be a ghosted boundary in x, in which case the neighbor is also */
459: const PetscInt neighbor = 1;
460: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
461: const PetscInt *const nNeighbor = nNeighbors[neighbor];
462: const PetscInt entriesPerElementRowNeighbor = entriesPerElementRow; /* same as here */
463: for (jghost = 0; jghost < ghostOffsetStart[1]; ++jghost) {
464: const PetscInt j = nNeighbor[1] - ghostOffsetStart[1] + jghost;
465: for (ighost = ghostOffsetStart[0]; ighost < stag->nGhost[0] - ghostOffsetEnd[0]; ++ighost) {
466: const PetscInt i = ighost - ghostOffsetStart[0];
467: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
468: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
469: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
470: }
471: }
472: if (dummyEnd[0]) {
473: const PetscInt ighost = stag->nGhost[0] - ghostOffsetEnd[0];
474: const PetscInt i = stag->n[0];
475: for (d = 0; d < stag->dof[0]; ++d, ++count) { /* Vertex */
476: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
477: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
478: }
479: for (d = 0; d < stag->dof[1]; ++d, ++count) { /* Face */
480: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + stag->dof[0] + d;
481: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + stag->dof[0] + stag->dof[1] + d;
482: }
483: }
484: }
485: }
487: /* Neighbor 2 (down right) */
488: if (!star && !dummyEnd[0] && !dummyStart[1]) {
489: const PetscInt neighbor = 2;
490: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
491: const PetscInt *const nNeighbor = nNeighbors[neighbor];
492: const PetscInt entriesPerElementRowNeighbor = nNeighbor[0] * stag->entriesPerElement + (nextToDummyEnd[0] ? entriesPerFace : 0);
493: for (jghost = 0; jghost < ghostOffsetStart[1]; ++jghost) {
494: const PetscInt j = nNeighbor[1] - ghostOffsetStart[1] + jghost;
495: for (i = 0; i < ghostOffsetEnd[0]; ++i) {
496: const PetscInt ighost = stag->nGhost[0] - ghostOffsetEnd[0] + i;
497: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
498: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
499: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
500: }
501: }
502: }
503: }
505: /* Neighbor 3 (left) */
506: if (!dummyStart[0]) {
507: /* Our neighbor is never a ghosted boundary in x, but we may be
508: Here, we may be a ghosted boundary in y and thus so will our neighbor be */
509: const PetscInt neighbor = 3;
510: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
511: const PetscInt *const nNeighbor = nNeighbors[neighbor];
512: const PetscInt entriesPerElementRowNeighbor = nNeighbor[0] * stag->entriesPerElement;
513: for (jghost = ghostOffsetStart[1]; jghost < stag->nGhost[1] - ghostOffsetEnd[1]; ++jghost) {
514: const PetscInt j = jghost - ghostOffsetStart[1];
515: for (ighost = 0; ighost < ghostOffsetStart[0]; ++ighost) {
516: const PetscInt i = nNeighbor[0] - ghostOffsetStart[0] + ighost;
517: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
518: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
519: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
520: }
521: }
522: }
523: if (dummyEnd[1]) {
524: const PetscInt jghost = stag->nGhost[1] - ghostOffsetEnd[1];
525: const PetscInt j = stag->n[1];
526: for (ighost = 0; ighost < ghostOffsetStart[0]; ++ighost) {
527: const PetscInt i = nNeighbor[0] - ghostOffsetStart[0] + ighost;
528: for (d = 0; d < entriesPerFace; ++d, ++count) { /* only vertices and horizontal face (which are the first dof) */
529: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * entriesPerFace + d; /* i moves by face here */
530: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
531: }
532: }
533: }
534: }
536: /* Interior/Resident-here-in-global elements ("Neighbor 4" - same rank)
537: *including* entries from boundary dummy elements */
538: {
539: const PetscInt neighbor = 4;
540: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
541: for (j = 0; j < stag->n[1]; ++j) {
542: const PetscInt jghost = j + ghostOffsetStart[1];
543: for (i = 0; i < stag->n[0]; ++i) {
544: const PetscInt ighost = i + ghostOffsetStart[0];
545: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
546: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * stag->entriesPerElement + d;
547: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
548: }
549: }
550: if (dummyEnd[0]) {
551: const PetscInt ighost = i + ghostOffsetStart[0];
552: i = stag->n[0];
553: for (d = 0; d < stag->dof[0]; ++d, ++count) { /* vertex first */
554: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * stag->entriesPerElement + d;
555: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
556: }
557: for (d = 0; d < stag->dof[1]; ++d, ++count) { /* then left edge (skipping bottom face) */
558: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * stag->entriesPerElement + stag->dof[0] + d;
559: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + stag->dof[0] + stag->dof[1] + d;
560: }
561: }
562: }
563: if (dummyEnd[1]) {
564: const PetscInt jghost = j + ghostOffsetStart[1];
565: j = stag->n[1];
566: for (i = 0; i < stag->n[0]; ++i) {
567: const PetscInt ighost = i + ghostOffsetStart[0];
568: for (d = 0; d < entriesPerFace; ++d, ++count) { /* vertex and bottom face (which are the first entries) */
569: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * entriesPerFace + d; /* note i increment by entriesPerFace */
570: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
571: }
572: }
573: if (dummyEnd[0]) {
574: const PetscInt ighost = i + ghostOffsetStart[0];
575: i = stag->n[0];
576: for (d = 0; d < entriesPerCorner; ++d, ++count) { /* vertex only */
577: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * entriesPerFace + d; /* note i increment by entriesPerFace */
578: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
579: }
580: }
581: }
582: }
584: /* Neighbor 5 (right) */
585: if (!dummyEnd[0]) {
586: /* We can never be right boundary, but we may be a top boundary, along with the right neighbor */
587: const PetscInt neighbor = 5;
588: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
589: const PetscInt *const nNeighbor = nNeighbors[neighbor];
590: const PetscInt entriesPerElementRowNeighbor = nNeighbor[0] * stag->entriesPerElement + (nextToDummyEnd[0] ? entriesPerFace : 0);
591: for (jghost = ghostOffsetStart[1]; jghost < stag->nGhost[1] - ghostOffsetEnd[1]; ++jghost) {
592: const PetscInt j = jghost - ghostOffsetStart[1];
593: for (i = 0; i < ghostOffsetEnd[0]; ++i) {
594: const PetscInt ighost = stag->nGhost[0] - ghostOffsetEnd[0] + i;
595: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
596: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
597: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
598: }
599: }
600: }
601: if (dummyEnd[1]) {
602: const PetscInt jghost = stag->nGhost[1] - ghostOffsetEnd[1];
603: const PetscInt j = nNeighbor[1];
604: for (i = 0; i < ghostOffsetEnd[0]; ++i) {
605: const PetscInt ighost = stag->nGhost[0] - ghostOffsetEnd[0] + i;
606: for (d = 0; d < entriesPerFace; ++d, ++count) { /* only vertices and horizontal face (which are the first dof) */
607: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * entriesPerFace + d; /* Note i increment by entriesPerFace */
608: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
609: }
610: }
611: }
612: }
614: /* Neighbor 6 (up left) */
615: if (!star && !dummyStart[0] && !dummyEnd[1]) {
616: /* We can never be a top boundary, but our neighbor may be
617: We may be a right boundary, but our neighbor cannot be */
618: const PetscInt neighbor = 6;
619: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
620: const PetscInt *const nNeighbor = nNeighbors[neighbor];
621: const PetscInt entriesPerElementRowNeighbor = nNeighbor[0] * stag->entriesPerElement;
622: for (j = 0; j < ghostOffsetEnd[1]; ++j) {
623: const PetscInt jghost = stag->nGhost[1] - ghostOffsetEnd[1] + j;
624: for (ighost = 0; ighost < ghostOffsetStart[0]; ++ighost) {
625: const PetscInt i = nNeighbor[0] - ghostOffsetStart[0] + ighost;
626: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
627: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
628: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
629: }
630: }
631: }
632: }
634: /* Neighbor 7 (up) */
635: if (!dummyEnd[1]) {
636: /* We cannot be the last rank in y, though our neighbor may be
637: We may be the last rank in x, in which case our neighbor is also */
638: const PetscInt neighbor = 7;
639: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
640: const PetscInt *const nNeighbor = nNeighbors[neighbor];
641: const PetscInt entriesPerElementRowNeighbor = entriesPerElementRow; /* same as here */
642: for (j = 0; j < ghostOffsetEnd[1]; ++j) {
643: const PetscInt jghost = stag->nGhost[1] - ghostOffsetEnd[1] + j;
644: for (ighost = ghostOffsetStart[0]; ighost < stag->nGhost[0] - ghostOffsetEnd[0]; ++ighost) {
645: const PetscInt i = ighost - ghostOffsetStart[0];
646: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
647: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
648: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
649: }
650: }
651: if (dummyEnd[0]) {
652: const PetscInt ighost = stag->nGhost[0] - ghostOffsetEnd[0];
653: const PetscInt i = nNeighbor[0];
654: for (d = 0; d < stag->dof[0]; ++d, ++count) { /* Vertex */
655: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
656: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
657: }
658: for (d = 0; d < stag->dof[1]; ++d, ++count) { /* Face */
659: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + stag->dof[0] + d;
660: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + stag->dof[0] + stag->dof[1] + d;
661: }
662: }
663: }
664: }
666: /* Neighbor 8 (up right) */
667: if (!star && !dummyEnd[0] && !dummyEnd[1]) {
668: /* We can never be a ghosted boundary
669: Our neighbor may be a top boundary, a right boundary, or both */
670: const PetscInt neighbor = 8;
671: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
672: const PetscInt *const nNeighbor = nNeighbors[neighbor];
673: const PetscInt entriesPerElementRowNeighbor = nNeighbor[0] * stag->entriesPerElement + (nextToDummyEnd[0] ? entriesPerFace : 0);
674: for (j = 0; j < ghostOffsetEnd[1]; ++j) {
675: const PetscInt jghost = stag->nGhost[1] - ghostOffsetEnd[1] + j;
676: for (i = 0; i < ghostOffsetEnd[0]; ++i) {
677: const PetscInt ighost = stag->nGhost[0] - ghostOffsetEnd[0] + i;
678: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
679: idxGlobal[count] = globalOffset + j * entriesPerElementRowNeighbor + i * stag->entriesPerElement + d;
680: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
681: }
682: }
683: }
684: }
686: PetscCheck(count == entriesToTransferTotal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of entries computed in gtol (%" PetscInt_FMT ") is not as expected (%" PetscInt_FMT ")", count, entriesToTransferTotal);
688: /* Create Local and Global ISs (transferring pointer ownership) */
689: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxLocal, PETSC_OWN_POINTER, &isLocal));
690: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), entriesToTransferTotal, idxGlobal, PETSC_OWN_POINTER, &isGlobal));
692: /* Create stag->gtol. The order is computed as PETSc ordering, and doesn't include dummy entries */
693: {
694: Vec local, global;
695: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global));
696: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, stag->entriesPerElement, stag->entriesGhost, NULL, &local));
697: PetscCall(VecScatterCreate(global, isGlobal, local, isLocal, &stag->gtol));
698: PetscCall(VecDestroy(&global));
699: PetscCall(VecDestroy(&local));
700: }
702: /* Destroy ISs */
703: PetscCall(ISDestroy(&isLocal));
704: PetscCall(ISDestroy(&isGlobal));
706: /* Next, we iterate over the local entries again, in local order, recording the global entry to which each maps,
707: or -1 if there is none */
708: PetscCall(PetscMalloc1(stag->entriesGhost, &idxGlobalAll));
710: countAll = 0;
712: /* Loop over rows 1/3 : down */
713: if (!dummyStart[1]) {
714: for (jghost = 0; jghost < ghostOffsetStart[1]; ++jghost) {
715: /* Loop over columns 1/3 : down left */
716: if (!star && !dummyStart[0]) {
717: const PetscInt neighbor = 0;
718: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
719: const PetscInt *const nNeighbor = nNeighbors[neighbor];
720: const PetscInt j = nNeighbor[1] - ghostOffsetStart[1] + jghost; /* Note: this is actually the same value for the whole row of ranks below, so recomputing it for the next two ranks is redundant, and one could even get rid of jghost entirely if desired */
721: const PetscInt eprNeighbor = nNeighbor[0] * stag->entriesPerElement;
722: for (i = nNeighbor[0] - ghostOffsetStart[0]; i < nNeighbor[0]; ++i) {
723: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + d;
724: }
725: } else {
726: /* Down Left dummies */
727: for (ighost = 0; ighost < ghostOffsetStart[0]; ++ighost) {
728: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
729: }
730: }
732: /* Loop over columns 2/3 : down middle */
733: {
734: const PetscInt neighbor = 1;
735: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
736: const PetscInt *const nNeighbor = nNeighbors[neighbor];
737: const PetscInt j = nNeighbor[1] - ghostOffsetStart[1] + jghost;
738: const PetscInt eprNeighbor = entriesPerElementRow; /* same as here */
739: for (i = 0; i < nNeighbor[0]; ++i) {
740: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + d;
741: }
742: }
744: /* Loop over columns 3/3 : down right */
745: if (!star && !dummyEnd[0]) {
746: const PetscInt neighbor = 2;
747: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
748: const PetscInt *const nNeighbor = nNeighbors[neighbor];
749: const PetscInt j = nNeighbor[1] - ghostOffsetStart[1] + jghost;
750: const PetscInt eprNeighbor = nNeighbor[0] * stag->entriesPerElement + (nextToDummyEnd[0] ? entriesPerFace : 0);
751: for (i = 0; i < ghostOffsetEnd[0]; ++i) {
752: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + d;
753: }
754: } else if (dummyEnd[0]) {
755: /* Down right partial dummy elements, living on the *down* rank */
756: const PetscInt neighbor = 1;
757: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
758: const PetscInt *const nNeighbor = nNeighbors[neighbor];
759: const PetscInt j = nNeighbor[1] - ghostOffsetStart[1] + jghost;
760: const PetscInt eprNeighbor = entriesPerElementRow; /* same as here */
761: PetscInt dGlobal;
762: i = nNeighbor[0];
763: for (d = 0, dGlobal = 0; d < stag->dof[0]; ++d, ++dGlobal, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + dGlobal;
764: for (; d < stag->dof[0] + stag->dof[1]; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy down face point */ }
765: for (; d < stag->dof[0] + 2 * stag->dof[1]; ++d, ++dGlobal, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + dGlobal;
766: for (; d < stag->entriesPerElement; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy element point */ }
767: ++i;
768: for (; i < nNeighbor[0] + ghostOffsetEnd[0]; ++i) {
769: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
770: }
771: } else {
772: /* Down Right dummies */
773: for (ighost = 0; ighost < ghostOffsetEnd[0]; ++ighost) {
774: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
775: }
776: }
777: }
778: } else {
779: /* Down dummies row */
780: for (jghost = 0; jghost < ghostOffsetStart[1]; ++jghost) {
781: for (ighost = 0; ighost < stag->nGhost[0]; ++ighost) {
782: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
783: }
784: }
785: }
787: /* Loop over rows 2/3 : center */
788: for (j = 0; j < stag->n[1]; ++j) {
789: /* Loop over columns 1/3 : left */
790: if (!dummyStart[0]) {
791: const PetscInt neighbor = 3;
792: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
793: const PetscInt *const nNeighbor = nNeighbors[neighbor];
794: const PetscInt eprNeighbor = nNeighbor[0] * stag->entriesPerElement;
795: for (i = nNeighbor[0] - ghostOffsetStart[0]; i < nNeighbor[0]; ++i) {
796: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + d;
797: }
798: } else {
799: /* (Middle) Left dummies */
800: for (ighost = 0; ighost < ghostOffsetStart[0]; ++ighost) {
801: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
802: }
803: }
805: /* Loop over columns 2/3 : here (the "neighbor" is ourselves, here) */
806: {
807: const PetscInt neighbor = 4;
808: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
809: const PetscInt eprNeighbor = entriesPerElementRow; /* same as here (obviously) */
810: for (i = 0; i < stag->n[0]; ++i) {
811: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + d;
812: }
813: }
815: /* Loop over columns 3/3 : right */
816: if (!dummyEnd[0]) {
817: const PetscInt neighbor = 5;
818: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
819: const PetscInt *const nNeighbor = nNeighbors[neighbor];
820: const PetscInt eprNeighbor = nNeighbor[0] * stag->entriesPerElement + (nextToDummyEnd[0] ? entriesPerFace : 0);
821: for (i = 0; i < ghostOffsetEnd[0]; ++i) {
822: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + d;
823: }
824: } else {
825: /* -1's for right layer of partial dummies, living on *this* rank */
826: const PetscInt neighbor = 4;
827: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
828: const PetscInt *const nNeighbor = nNeighbors[neighbor];
829: const PetscInt eprNeighbor = entriesPerElementRow; /* same as here (obviously) */
830: PetscInt dGlobal;
831: i = nNeighbor[0];
832: for (d = 0, dGlobal = 0; d < stag->dof[0]; ++d, ++dGlobal, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + dGlobal;
833: for (; d < stag->dof[0] + stag->dof[1]; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy down face point */ }
834: for (; d < stag->dof[0] + 2 * stag->dof[1]; ++d, ++dGlobal, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + dGlobal;
835: for (; d < stag->entriesPerElement; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy element point */ }
836: ++i;
837: for (; i < nNeighbor[0] + ghostOffsetEnd[0]; ++i) {
838: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
839: }
840: }
841: }
843: /* Loop over rows 3/3 : up */
844: if (!dummyEnd[1]) {
845: for (j = 0; j < ghostOffsetEnd[1]; ++j) {
846: /* Loop over columns 1/3 : up left */
847: if (!star && !dummyStart[0]) {
848: const PetscInt neighbor = 6;
849: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
850: const PetscInt *const nNeighbor = nNeighbors[neighbor];
851: const PetscInt eprNeighbor = nNeighbor[0] * stag->entriesPerElement;
852: for (i = nNeighbor[0] - ghostOffsetStart[0]; i < nNeighbor[0]; ++i) {
853: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + d;
854: }
855: } else {
856: /* Up Left dummies */
857: for (ighost = 0; ighost < ghostOffsetStart[0]; ++ighost) {
858: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
859: }
860: }
862: /* Loop over columns 2/3 : up */
863: {
864: const PetscInt neighbor = 7;
865: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
866: const PetscInt *const nNeighbor = nNeighbors[neighbor];
867: const PetscInt eprNeighbor = entriesPerElementRow; /* Same as here */
868: for (i = 0; i < nNeighbor[0]; ++i) {
869: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + d;
870: }
871: }
873: /* Loop over columns 3/3 : up right */
874: if (!star && !dummyEnd[0]) {
875: const PetscInt neighbor = 8;
876: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
877: const PetscInt *const nNeighbor = nNeighbors[neighbor];
878: const PetscInt eprNeighbor = nNeighbor[0] * stag->entriesPerElement + (nextToDummyEnd[0] ? entriesPerFace : 0);
879: for (i = 0; i < ghostOffsetEnd[0]; ++i) {
880: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + d;
881: }
882: } else if (dummyEnd[0]) {
883: /* -1's for right layer of partial dummies, living on rank above */
884: const PetscInt neighbor = 7;
885: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
886: const PetscInt *const nNeighbor = nNeighbors[neighbor];
887: const PetscInt eprNeighbor = entriesPerElementRow; /* Same as here */
888: PetscInt dGlobal;
889: i = nNeighbor[0];
890: for (d = 0, dGlobal = 0; d < stag->dof[0]; ++d, ++dGlobal, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + dGlobal;
891: for (; d < stag->dof[0] + stag->dof[1]; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy down face point */ }
892: for (; d < stag->dof[0] + 2 * stag->dof[1]; ++d, ++dGlobal, ++countAll) idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * stag->entriesPerElement + dGlobal;
893: for (; d < stag->entriesPerElement; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy element point */ }
894: ++i;
895: for (; i < nNeighbor[0] + ghostOffsetEnd[0]; ++i) {
896: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
897: }
898: } else {
899: /* Up Right dummies */
900: for (ighost = 0; ighost < ghostOffsetEnd[0]; ++ighost) {
901: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
902: }
903: }
904: }
905: } else {
906: j = stag->n[1];
907: /* Top layer of partial dummies */
909: /* up left partial dummies layer : Loop over columns 1/3 : living on *left* neighbor */
910: if (!dummyStart[0]) {
911: const PetscInt neighbor = 3;
912: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
913: const PetscInt *const nNeighbor = nNeighbors[neighbor];
914: const PetscInt eprNeighbor = nNeighbor[0] * stag->entriesPerElement;
915: for (i = nNeighbor[0] - ghostOffsetStart[0]; i < nNeighbor[0]; ++i) {
916: for (d = 0; d < stag->dof[0] + stag->dof[1]; ++d, ++countAll) { idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * entriesPerFace + d; /* Note entriesPerFace here */ }
917: for (; d < stag->entriesPerElement; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy left face and element points */ }
918: }
919: } else {
920: for (ighost = 0; ighost < ghostOffsetStart[0]; ++ighost) {
921: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
922: }
923: }
925: /* up partial dummies layer : Loop over columns 2/3 : living on *this* rank */
926: {
927: const PetscInt neighbor = 4;
928: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
929: const PetscInt eprNeighbor = entriesPerElementRow; /* same as here (obviously) */
930: for (i = 0; i < stag->n[0]; ++i) {
931: for (d = 0; d < stag->dof[0] + stag->dof[1]; ++d, ++countAll) { idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * entriesPerFace + d; /* Note entriesPerFace here */ }
932: for (; d < stag->entriesPerElement; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy left face and element points */ }
933: }
934: }
936: if (!dummyEnd[0]) {
937: /* up right partial dummies layer : Loop over columns 3/3 : living on *right* neighbor */
938: const PetscInt neighbor = 5;
939: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
940: const PetscInt *const nNeighbor = nNeighbors[neighbor];
941: const PetscInt eprNeighbor = nNeighbor[0] * stag->entriesPerElement + (nextToDummyEnd[0] ? entriesPerFace : 0);
942: for (i = 0; i < ghostOffsetEnd[0]; ++i) {
943: for (d = 0; d < stag->dof[0] + stag->dof[1]; ++d, ++countAll) { idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * entriesPerFace + d; /* Note entriesPerFace here */ }
944: for (; d < stag->entriesPerElement; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy left face and element points */ }
945: }
946: } else {
947: /* Top partial dummies layer : Loop over columns 3/3 : right, living *here* */
948: const PetscInt neighbor = 4;
949: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
950: const PetscInt eprNeighbor = entriesPerElementRow; /* same as here (obviously) */
951: i = stag->n[0];
952: for (d = 0; d < stag->dof[0]; ++d, ++countAll) { /* Note just the vertex here */
953: idxGlobalAll[countAll] = globalOffset + j * eprNeighbor + i * entriesPerFace + d; /* Note entriesPerFace here */
954: }
955: for (; d < stag->entriesPerElement; ++d, ++countAll) { idxGlobalAll[countAll] = -1; /* dummy bottom face, left face and element points */ }
956: ++i;
957: for (; i < stag->n[0] + ghostOffsetEnd[0]; ++i) {
958: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
959: }
960: }
961: ++j;
962: /* Additional top dummy layers */
963: for (; j < stag->n[1] + ghostOffsetEnd[1]; ++j) {
964: for (ighost = 0; ighost < stag->nGhost[0]; ++ighost) {
965: for (d = 0; d < stag->entriesPerElement; ++d, ++countAll) idxGlobalAll[countAll] = -1;
966: }
967: }
968: }
970: /* Create local-to-global map (in local ordering, includes maps to -1 for dummy points) */
971: PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stag->entriesGhost, idxGlobalAll, PETSC_OWN_POINTER, &dm->ltogmap));
972: }
974: /* In special cases, create a dedicated injective local-to-global map */
975: if ((stag->boundaryType[0] == DM_BOUNDARY_PERIODIC && stag->nRanks[0] == 1) || (stag->boundaryType[1] == DM_BOUNDARY_PERIODIC && stag->nRanks[1] == 1)) PetscCall(DMStagPopulateLocalToGlobalInjective(dm));
977: /* Free global offsets */
978: PetscCall(PetscFree(globalOffsets));
980: /* Precompute location offsets */
981: PetscCall(DMStagComputeLocationOffsets_2d(dm));
983: /* View from Options */
984: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: /* adapted from da2.c */
989: static PetscErrorCode DMStagSetUpBuildRankGrid_2d(DM dm)
990: {
991: DM_Stag *const stag = (DM_Stag *)dm->data;
992: PetscInt m, n;
993: PetscMPIInt rank, size;
994: const PetscInt M = stag->N[0];
995: const PetscInt N = stag->N[1];
997: PetscFunctionBegin;
998: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
999: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1000: m = stag->nRanks[0];
1001: n = stag->nRanks[1];
1002: if (m != PETSC_DECIDE) {
1003: PetscCheck(m >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of ranks in X direction: %" PetscInt_FMT, m);
1004: PetscCheck(m <= size, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Too many ranks in X direction: %" PetscInt_FMT " %d", m, size);
1005: }
1006: if (n != PETSC_DECIDE) {
1007: PetscCheck(n >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of ranks in Y direction: %" PetscInt_FMT, n);
1008: PetscCheck(n <= size, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Too many ranks in Y direction: %" PetscInt_FMT " %d", n, size);
1009: }
1010: if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
1011: if (n != PETSC_DECIDE) {
1012: m = size / n;
1013: } else if (m != PETSC_DECIDE) {
1014: n = size / m;
1015: } else {
1016: /* try for squarish distribution */
1017: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N)));
1018: if (!m) m = 1;
1019: while (m > 0) {
1020: n = size / m;
1021: if (m * n == size) break;
1022: m--;
1023: }
1024: if (M > N && m < n) {
1025: PetscInt _m = m;
1026: m = n;
1027: n = _m;
1028: }
1029: }
1030: PetscCheck(m * n == size, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Unable to create partition, check the size of the communicator and input m and n ");
1031: } else PetscCheck(m * n == size, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Given Bad partition. Product of sizes (%" PetscInt_FMT ") does not equal communicator size (%d)", m * n, size);
1032: PetscCheck(M >= m, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Partition in x direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, M, m);
1033: PetscCheck(N >= n, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, N, n);
1034: stag->nRanks[0] = m;
1035: stag->nRanks[1] = n;
1036: PetscFunctionReturn(PETSC_SUCCESS);
1037: }
1039: static PetscErrorCode DMStagSetUpBuildNeighbors_2d(DM dm)
1040: {
1041: DM_Stag *const stag = (DM_Stag *)dm->data;
1042: PetscInt d, i;
1043: PetscBool per[2], first[2], last[2];
1044: PetscInt neighborRank[9][2], r[2], n[2];
1045: const PetscInt dim = 2;
1047: PetscFunctionBegin;
1048: for (d = 0; d < dim; ++d)
1049: PetscCheck(stag->boundaryType[d] == DM_BOUNDARY_NONE || stag->boundaryType[d] == DM_BOUNDARY_PERIODIC || stag->boundaryType[d] == DM_BOUNDARY_GHOSTED, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Neighbor determination not implemented for %s",
1050: DMBoundaryTypes[stag->boundaryType[d]]);
1052: /* Assemble some convenience variables */
1053: for (d = 0; d < dim; ++d) {
1054: per[d] = (PetscBool)(stag->boundaryType[d] == DM_BOUNDARY_PERIODIC);
1055: first[d] = stag->firstRank[d];
1056: last[d] = stag->lastRank[d];
1057: r[d] = stag->rank[d];
1058: n[d] = stag->nRanks[d];
1059: }
1061: /* First, compute the position in the rank grid for all neighbors */
1062: neighborRank[0][0] = first[0] ? (per[0] ? n[0] - 1 : -1) : r[0] - 1; /* left down */
1063: neighborRank[0][1] = first[1] ? (per[1] ? n[1] - 1 : -1) : r[1] - 1;
1065: neighborRank[1][0] = r[0]; /* down */
1066: neighborRank[1][1] = first[1] ? (per[1] ? n[1] - 1 : -1) : r[1] - 1;
1068: neighborRank[2][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right down */
1069: neighborRank[2][1] = first[1] ? (per[1] ? n[1] - 1 : -1) : r[1] - 1;
1071: neighborRank[3][0] = first[0] ? (per[0] ? n[0] - 1 : -1) : r[0] - 1; /* left */
1072: neighborRank[3][1] = r[1];
1074: neighborRank[4][0] = r[0]; /* */
1075: neighborRank[4][1] = r[1];
1077: neighborRank[5][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right */
1078: neighborRank[5][1] = r[1];
1080: neighborRank[6][0] = first[0] ? (per[0] ? n[0] - 1 : -1) : r[0] - 1; /* left up */
1081: neighborRank[6][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
1083: neighborRank[7][0] = r[0]; /* up */
1084: neighborRank[7][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
1086: neighborRank[8][0] = last[0] ? (per[0] ? 0 : -1) : r[0] + 1; /* right up */
1087: neighborRank[8][1] = last[1] ? (per[1] ? 0 : -1) : r[1] + 1;
1089: /* Then, compute the rank of each in the linear ordering */
1090: PetscCall(PetscMalloc1(9, &stag->neighbors));
1091: for (i = 0; i < 9; ++i) {
1092: if (neighborRank[i][0] >= 0 && neighborRank[i][1] >= 0) {
1093: stag->neighbors[i] = neighborRank[i][0] + n[0] * neighborRank[i][1];
1094: } else {
1095: stag->neighbors[i] = -1;
1096: }
1097: }
1098: PetscFunctionReturn(PETSC_SUCCESS);
1099: }
1101: static PetscErrorCode DMStagSetUpBuildGlobalOffsets_2d(DM dm, PetscInt **pGlobalOffsets)
1102: {
1103: const DM_Stag *const stag = (DM_Stag *)dm->data;
1104: PetscInt *globalOffsets;
1105: PetscInt i, j, d, entriesPerFace, count;
1106: PetscMPIInt size;
1107: PetscBool extra[2];
1109: PetscFunctionBegin;
1110: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
1111: for (d = 0; d < 2; ++d) extra[d] = (PetscBool)(stag->boundaryType[d] != DM_BOUNDARY_PERIODIC); /* Extra points in global rep */
1112: entriesPerFace = stag->dof[0] + stag->dof[1];
1113: PetscCall(PetscMalloc1(size, pGlobalOffsets));
1114: globalOffsets = *pGlobalOffsets;
1115: globalOffsets[0] = 0;
1116: count = 1; /* note the count is offset by 1 here. We add the size of the previous rank */
1117: for (j = 0; j < stag->nRanks[1] - 1; ++j) {
1118: const PetscInt nnj = stag->l[1][j];
1119: for (i = 0; i < stag->nRanks[0] - 1; ++i) {
1120: const PetscInt nni = stag->l[0][i];
1121: globalOffsets[count] = globalOffsets[count - 1] + nnj * nni * stag->entriesPerElement; /* No right/top/front boundaries */
1122: ++count;
1123: }
1124: {
1125: /* i = stag->nRanks[0]-1; */
1126: const PetscInt nni = stag->l[0][i];
1127: globalOffsets[count] = globalOffsets[count - 1] + nnj * nni * stag->entriesPerElement + (extra[0] ? nnj * entriesPerFace : 0); /* Extra faces on the right */
1128: ++count;
1129: }
1130: }
1131: {
1132: /* j = stag->nRanks[1]-1; */
1133: const PetscInt nnj = stag->l[1][j];
1134: for (i = 0; i < stag->nRanks[0] - 1; ++i) {
1135: const PetscInt nni = stag->l[0][i];
1136: globalOffsets[count] = globalOffsets[count - 1] + nni * nnj * stag->entriesPerElement + (extra[1] ? nni * entriesPerFace : 0); /* Extra faces on the top */
1137: ++count;
1138: }
1139: /* Don't need to compute entries in last element */
1140: }
1141: PetscFunctionReturn(PETSC_SUCCESS);
1142: }
1144: static PetscErrorCode DMStagComputeLocationOffsets_2d(DM dm)
1145: {
1146: DM_Stag *const stag = (DM_Stag *)dm->data;
1147: const PetscInt epe = stag->entriesPerElement;
1148: const PetscInt epr = stag->nGhost[0] * epe;
1150: PetscFunctionBegin;
1151: PetscCall(PetscMalloc1(DMSTAG_NUMBER_LOCATIONS, &stag->locationOffsets));
1152: stag->locationOffsets[DMSTAG_DOWN_LEFT] = 0;
1153: stag->locationOffsets[DMSTAG_DOWN] = stag->locationOffsets[DMSTAG_DOWN_LEFT] + stag->dof[0];
1154: stag->locationOffsets[DMSTAG_DOWN_RIGHT] = stag->locationOffsets[DMSTAG_DOWN_LEFT] + epe;
1155: stag->locationOffsets[DMSTAG_LEFT] = stag->locationOffsets[DMSTAG_DOWN] + stag->dof[1];
1156: stag->locationOffsets[DMSTAG_ELEMENT] = stag->locationOffsets[DMSTAG_LEFT] + stag->dof[1];
1157: stag->locationOffsets[DMSTAG_RIGHT] = stag->locationOffsets[DMSTAG_LEFT] + epe;
1158: stag->locationOffsets[DMSTAG_UP_LEFT] = stag->locationOffsets[DMSTAG_DOWN_LEFT] + epr;
1159: stag->locationOffsets[DMSTAG_UP] = stag->locationOffsets[DMSTAG_DOWN] + epr;
1160: stag->locationOffsets[DMSTAG_UP_RIGHT] = stag->locationOffsets[DMSTAG_UP_LEFT] + epe;
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToGlobalInjective_2d(DM dm)
1165: {
1166: DM_Stag *const stag = (DM_Stag *)dm->data;
1167: PetscInt *idxLocal, *idxGlobal, *globalOffsetsRecomputed;
1168: const PetscInt *globalOffsets;
1169: PetscInt i, j, d, count, entriesPerCorner, entriesPerFace, entriesPerElementRowGhost, entriesPerElementRow, ghostOffsetStart[2];
1170: IS isLocal, isGlobal;
1171: PetscBool dummyEnd[2];
1173: PetscFunctionBegin;
1174: PetscCall(DMStagSetUpBuildGlobalOffsets_2d(dm, &globalOffsetsRecomputed)); /* note that we don't actually use all of these. An available optimization is to pass them, when available */
1175: globalOffsets = globalOffsetsRecomputed;
1176: PetscCall(PetscMalloc1(stag->entries, &idxLocal));
1177: PetscCall(PetscMalloc1(stag->entries, &idxGlobal));
1178: for (d = 0; d < 2; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
1179: entriesPerCorner = stag->dof[0];
1180: entriesPerFace = stag->dof[0] + stag->dof[1];
1181: entriesPerElementRow = stag->n[0] * stag->entriesPerElement + (dummyEnd[0] ? entriesPerFace : 0);
1182: entriesPerElementRowGhost = stag->nGhost[0] * stag->entriesPerElement;
1183: count = 0;
1184: for (d = 0; d < 2; ++d) ghostOffsetStart[d] = stag->start[d] - stag->startGhost[d];
1185: {
1186: const PetscInt neighbor = 4;
1187: const PetscInt globalOffset = globalOffsets[stag->neighbors[neighbor]];
1188: for (j = 0; j < stag->n[1]; ++j) {
1189: const PetscInt jghost = j + ghostOffsetStart[1];
1190: for (i = 0; i < stag->n[0]; ++i) {
1191: const PetscInt ighost = i + ghostOffsetStart[0];
1192: for (d = 0; d < stag->entriesPerElement; ++d, ++count) {
1193: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * stag->entriesPerElement + d;
1194: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
1195: }
1196: }
1197: if (dummyEnd[0]) {
1198: const PetscInt ighost = i + ghostOffsetStart[0];
1199: i = stag->n[0];
1200: for (d = 0; d < stag->dof[0]; ++d, ++count) { /* vertex first */
1201: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * stag->entriesPerElement + d;
1202: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
1203: }
1204: for (d = 0; d < stag->dof[1]; ++d, ++count) { /* then left edge (skipping bottom face) */
1205: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * stag->entriesPerElement + stag->dof[0] + d;
1206: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + stag->dof[0] + stag->dof[1] + d;
1207: }
1208: }
1209: }
1210: if (dummyEnd[1]) {
1211: const PetscInt jghost = j + ghostOffsetStart[1];
1212: j = stag->n[1];
1213: for (i = 0; i < stag->n[0]; ++i) {
1214: const PetscInt ighost = i + ghostOffsetStart[0];
1215: for (d = 0; d < entriesPerFace; ++d, ++count) { /* vertex and bottom face (which are the first entries) */
1216: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * entriesPerFace + d; /* note i increment by entriesPerFace */
1217: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
1218: }
1219: }
1220: if (dummyEnd[0]) {
1221: const PetscInt ighost = i + ghostOffsetStart[0];
1222: i = stag->n[0];
1223: for (d = 0; d < entriesPerCorner; ++d, ++count) { /* vertex only */
1224: idxGlobal[count] = globalOffset + j * entriesPerElementRow + i * entriesPerFace + d; /* note i increment by entriesPerFace */
1225: idxLocal[count] = jghost * entriesPerElementRowGhost + ighost * stag->entriesPerElement + d;
1226: }
1227: }
1228: }
1229: }
1230: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxLocal, PETSC_OWN_POINTER, &isLocal));
1231: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), stag->entries, idxGlobal, PETSC_OWN_POINTER, &isGlobal));
1232: {
1233: Vec local, global;
1234: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm), 1, stag->entries, PETSC_DECIDE, NULL, &global));
1235: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, stag->entriesPerElement, stag->entriesGhost, NULL, &local));
1236: PetscCall(VecScatterCreate(local, isLocal, global, isGlobal, &stag->ltog_injective));
1237: PetscCall(VecDestroy(&global));
1238: PetscCall(VecDestroy(&local));
1239: }
1240: PetscCall(ISDestroy(&isLocal));
1241: PetscCall(ISDestroy(&isGlobal));
1242: if (globalOffsetsRecomputed) PetscCall(PetscFree(globalOffsetsRecomputed));
1243: PetscFunctionReturn(PETSC_SUCCESS);
1244: }
1246: PETSC_INTERN PetscErrorCode DMStagPopulateLocalToLocal2d_Internal(DM dm)
1247: {
1248: DM_Stag *const stag = (DM_Stag *)dm->data;
1249: PetscInt *idxRemap;
1250: PetscBool dummyEnd[2];
1251: PetscInt i, j, d, count, leftGhostElements, downGhostElements, entriesPerRowGhost, iOffset, jOffset;
1252: PetscInt dOffset[4] = {0};
1254: PetscFunctionBegin;
1255: PetscCall(VecScatterCopy(stag->gtol, &stag->ltol));
1256: PetscCall(PetscMalloc1(stag->entries, &idxRemap));
1258: for (d = 0; d < 2; ++d) dummyEnd[d] = (PetscBool)(stag->lastRank[d] && stag->boundaryType[d] != DM_BOUNDARY_PERIODIC);
1259: leftGhostElements = stag->start[0] - stag->startGhost[0];
1260: downGhostElements = stag->start[1] - stag->startGhost[1];
1261: entriesPerRowGhost = stag->nGhost[0] * stag->entriesPerElement;
1262: dOffset[1] = dOffset[0] + stag->dof[0];
1263: dOffset[2] = dOffset[1] + stag->dof[1];
1264: dOffset[3] = dOffset[2] + stag->dof[1];
1266: count = 0;
1267: for (j = 0; j < stag->n[1]; ++j) {
1268: jOffset = entriesPerRowGhost * (downGhostElements + j);
1269: for (i = 0; i < stag->n[0]; ++i) {
1270: iOffset = stag->entriesPerElement * (leftGhostElements + i);
1271: // all
1272: for (d = 0; d < stag->entriesPerElement; ++d) idxRemap[count++] = jOffset + iOffset + d;
1273: }
1274: if (dummyEnd[0]) {
1275: iOffset = stag->entriesPerElement * (leftGhostElements + stag->n[0]);
1276: // down left, left
1277: for (d = 0; d < stag->dof[0]; ++d) idxRemap[count++] = jOffset + iOffset + dOffset[0] + d;
1278: for (d = 0; d < stag->dof[1]; ++d) idxRemap[count++] = jOffset + iOffset + dOffset[2] + d;
1279: }
1280: }
1281: if (dummyEnd[1]) {
1282: jOffset = entriesPerRowGhost * (downGhostElements + stag->n[1]);
1283: for (i = 0; i < stag->n[0]; ++i) {
1284: iOffset = stag->entriesPerElement * (leftGhostElements + i);
1285: // down left, down
1286: for (d = 0; d < stag->dof[0]; ++d) idxRemap[count++] = jOffset + iOffset + dOffset[0] + d;
1287: for (d = 0; d < stag->dof[1]; ++d) idxRemap[count++] = jOffset + iOffset + dOffset[1] + d;
1288: }
1289: if (dummyEnd[0]) {
1290: iOffset = stag->entriesPerElement * (leftGhostElements + stag->n[0]);
1291: // down left
1292: for (d = 0; d < stag->dof[0]; ++d) idxRemap[count++] = jOffset + iOffset + dOffset[0] + d;
1293: }
1294: }
1296: PetscCheck(count == stag->entries, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of entries computed in ltol (%" PetscInt_FMT ") is not as expected (%" PetscInt_FMT ")", count, stag->entries);
1298: PetscCall(VecScatterRemap(stag->ltol, idxRemap, NULL));
1299: PetscCall(PetscFree(idxRemap));
1300: PetscFunctionReturn(PETSC_SUCCESS);
1301: }
1303: PETSC_INTERN PetscErrorCode DMCreateMatrix_Stag_2D_AIJ_Assemble(DM dm, Mat A)
1304: {
1305: PetscInt entries, dof[DMSTAG_MAX_STRATA], epe, stencil_width, N[2], start[2], n[2], n_extra[2];
1306: DMStagStencilType stencil_type;
1307: DMBoundaryType boundary_type[2];
1309: PetscFunctionBegin;
1310: PetscCall(DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], NULL));
1311: PetscCall(DMStagGetStencilType(dm, &stencil_type));
1312: PetscCall(DMStagGetStencilWidth(dm, &stencil_width));
1313: PetscCall(DMStagGetEntries(dm, &entries));
1314: PetscCall(DMStagGetEntriesPerElement(dm, &epe));
1315: PetscCall(DMStagGetCorners(dm, &start[0], &start[1], NULL, &n[0], &n[1], NULL, &n_extra[0], &n_extra[1], NULL));
1316: PetscCall(DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL));
1317: PetscCall(DMStagGetBoundaryTypes(dm, &boundary_type[0], &boundary_type[1], NULL));
1319: if (stencil_type == DMSTAG_STENCIL_NONE) {
1320: /* Couple all DOF at each location to each other */
1321: DMStagStencil *row_vertex, *row_face_down, *row_face_left, *row_element;
1323: PetscCall(PetscMalloc1(dof[0], &row_vertex));
1324: for (PetscInt c = 0; c < dof[0]; ++c) {
1325: row_vertex[c].loc = DMSTAG_DOWN_LEFT;
1326: row_vertex[c].c = c;
1327: }
1329: PetscCall(PetscMalloc1(dof[1], &row_face_down));
1330: for (PetscInt c = 0; c < dof[1]; ++c) {
1331: row_face_down[c].loc = DMSTAG_DOWN;
1332: row_face_down[c].c = c;
1333: }
1335: PetscCall(PetscMalloc1(dof[1], &row_face_left));
1336: for (PetscInt c = 0; c < dof[1]; ++c) {
1337: row_face_left[c].loc = DMSTAG_LEFT;
1338: row_face_left[c].c = c;
1339: }
1341: PetscCall(PetscMalloc1(dof[2], &row_element));
1342: for (PetscInt c = 0; c < dof[2]; ++c) {
1343: row_element[c].loc = DMSTAG_ELEMENT;
1344: row_element[c].c = c;
1345: }
1347: for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
1348: for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
1349: {
1350: for (PetscInt c = 0; c < dof[0]; ++c) {
1351: row_vertex[c].i = ex;
1352: row_vertex[c].j = ey;
1353: }
1354: PetscCall(DMStagMatSetValuesStencil(dm, A, dof[0], row_vertex, dof[0], row_vertex, NULL, INSERT_VALUES));
1355: }
1356: if (ex < N[0]) {
1357: for (PetscInt c = 0; c < dof[1]; ++c) {
1358: row_face_down[c].i = ex;
1359: row_face_down[c].j = ey;
1360: }
1361: PetscCall(DMStagMatSetValuesStencil(dm, A, dof[1], row_face_down, dof[1], row_face_down, NULL, INSERT_VALUES));
1362: }
1363: if (ey < N[1]) {
1364: for (PetscInt c = 0; c < dof[1]; ++c) {
1365: row_face_left[c].i = ex;
1366: row_face_left[c].j = ey;
1367: }
1368: PetscCall(DMStagMatSetValuesStencil(dm, A, dof[1], row_face_left, dof[1], row_face_left, NULL, INSERT_VALUES));
1369: }
1370: if (ex < N[0] && ey < N[1]) {
1371: for (PetscInt c = 0; c < dof[2]; ++c) {
1372: row_element[c].i = ex;
1373: row_element[c].j = ey;
1374: }
1375: PetscCall(DMStagMatSetValuesStencil(dm, A, dof[2], row_element, dof[2], row_element, NULL, INSERT_VALUES));
1376: }
1377: }
1378: }
1379: PetscCall(PetscFree(row_vertex));
1380: PetscCall(PetscFree(row_face_left));
1381: PetscCall(PetscFree(row_face_down));
1382: PetscCall(PetscFree(row_element));
1383: } else if (stencil_type == DMSTAG_STENCIL_STAR || stencil_type == DMSTAG_STENCIL_BOX) {
1384: DMStagStencil *col, *row;
1386: PetscCall(PetscMalloc1(epe, &row));
1387: {
1388: PetscInt nrows = 0;
1390: for (PetscInt c = 0; c < dof[0]; ++c) {
1391: row[nrows].c = c;
1392: row[nrows].loc = DMSTAG_DOWN_LEFT;
1393: ++nrows;
1394: }
1395: for (PetscInt c = 0; c < dof[1]; ++c) {
1396: row[nrows].c = c;
1397: row[nrows].loc = DMSTAG_LEFT;
1398: ++nrows;
1399: }
1400: for (PetscInt c = 0; c < dof[1]; ++c) {
1401: row[nrows].c = c;
1402: row[nrows].loc = DMSTAG_DOWN;
1403: ++nrows;
1404: }
1405: for (PetscInt c = 0; c < dof[2]; ++c) {
1406: row[nrows].c = c;
1407: row[nrows].loc = DMSTAG_ELEMENT;
1408: ++nrows;
1409: }
1410: }
1412: PetscCall(PetscMalloc1(epe, &col));
1413: {
1414: PetscInt ncols = 0;
1416: for (PetscInt c = 0; c < dof[0]; ++c) {
1417: col[ncols].c = c;
1418: col[ncols].loc = DMSTAG_DOWN_LEFT;
1419: ++ncols;
1420: }
1421: for (PetscInt c = 0; c < dof[1]; ++c) {
1422: col[ncols].c = c;
1423: col[ncols].loc = DMSTAG_LEFT;
1424: ++ncols;
1425: }
1426: for (PetscInt c = 0; c < dof[1]; ++c) {
1427: col[ncols].c = c;
1428: col[ncols].loc = DMSTAG_DOWN;
1429: ++ncols;
1430: }
1431: for (PetscInt c = 0; c < dof[2]; ++c) {
1432: col[ncols].c = c;
1433: col[ncols].loc = DMSTAG_ELEMENT;
1434: ++ncols;
1435: }
1436: }
1438: for (PetscInt ey = start[1]; ey < start[1] + n[1] + n_extra[1]; ++ey) {
1439: for (PetscInt ex = start[0]; ex < start[0] + n[0] + n_extra[0]; ++ex) {
1440: for (PetscInt i = 0; i < epe; ++i) {
1441: row[i].i = ex;
1442: row[i].j = ey;
1443: }
1444: for (PetscInt offset_y = -stencil_width; offset_y <= stencil_width; ++offset_y) {
1445: const PetscInt ey_offset = ey + offset_y;
1446: for (PetscInt offset_x = -stencil_width; offset_x <= stencil_width; ++offset_x) {
1447: const PetscInt ex_offset = ex + offset_x;
1448: /* Only set values corresponding to elements which can have non-dummy entries,
1449: meaning those that map to unknowns in the global representation. In the periodic
1450: case, this is the entire stencil, but in all other cases, only includes a single
1451: "extra" element which is partially outside the physical domain (those points in the
1452: global representation */
1453: if ((stencil_type == DMSTAG_STENCIL_BOX || offset_x == 0 || offset_y == 0) && (boundary_type[0] == DM_BOUNDARY_PERIODIC || (ex_offset < N[0] + 1 && ex_offset >= 0)) && (boundary_type[1] == DM_BOUNDARY_PERIODIC || (ey_offset < N[1] + 1 && ey_offset >= 0))) {
1454: for (PetscInt i = 0; i < epe; ++i) {
1455: col[i].i = ex_offset;
1456: col[i].j = ey_offset;
1457: }
1458: PetscCall(DMStagMatSetValuesStencil(dm, A, epe, row, epe, col, NULL, INSERT_VALUES));
1459: }
1460: }
1461: }
1462: }
1463: }
1464: PetscCall(PetscFree(row));
1465: PetscCall(PetscFree(col));
1466: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported stencil type %s", DMStagStencilTypes[stencil_type]);
1467: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1468: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1469: PetscFunctionReturn(PETSC_SUCCESS);
1470: }