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