1 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
2 #include <petscdraw.h>
3
DMView_DA_2d(DM da,PetscViewer viewer)4 static PetscErrorCode DMView_DA_2d(DM da, PetscViewer viewer)
5 {
6 PetscMPIInt rank;
7 PetscBool isascii, isdraw, isglvis, isbinary;
8 DM_DA *dd = (DM_DA *)da->data;
9 #if defined(PETSC_HAVE_MATLAB)
10 PetscBool ismatlab;
11 #endif
12
13 PetscFunctionBegin;
14 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank));
15
16 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
17 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
18 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
19 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
20 #if defined(PETSC_HAVE_MATLAB)
21 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab));
22 #endif
23 if (isascii) {
24 PetscViewerFormat format;
25
26 PetscCall(PetscViewerGetFormat(viewer, &format));
27 if (format == PETSC_VIEWER_LOAD_BALANCE) {
28 PetscInt i, nmax = 0, nmin = PETSC_INT_MAX, navg = 0, *nz, nzlocal;
29 DMDALocalInfo info;
30 PetscMPIInt size;
31 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
32 PetscCall(DMDAGetLocalInfo(da, &info));
33 nzlocal = info.xm * info.ym;
34 PetscCall(PetscMalloc1(size, &nz));
35 PetscCallMPI(MPI_Allgather(&nzlocal, 1, MPIU_INT, nz, 1, MPIU_INT, PetscObjectComm((PetscObject)da)));
36 for (i = 0; i < size; i++) {
37 nmax = PetscMax(nmax, nz[i]);
38 nmin = PetscMin(nmin, nz[i]);
39 navg += nz[i];
40 }
41 PetscCall(PetscFree(nz));
42 navg = navg / size;
43 PetscCall(PetscViewerASCIIPrintf(viewer, " Load Balance - Grid Points: Min %" PetscInt_FMT " avg %" PetscInt_FMT " max %" PetscInt_FMT "\n", nmin, navg, nmax));
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46 if (format != PETSC_VIEWER_ASCII_GLVIS) {
47 DMDALocalInfo info;
48 PetscCall(DMDAGetLocalInfo(da, &info));
49 PetscCall(PetscViewerASCIIPushSynchronized(viewer));
50 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Processor [%d] M %" PetscInt_FMT " N %" PetscInt_FMT " m %" PetscInt_FMT " n %" PetscInt_FMT " w %" PetscInt_FMT " s %" PetscInt_FMT "\n", rank, dd->M, dd->N, dd->m, dd->n, dd->w, dd->s));
51 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "X range of indices: %" PetscInt_FMT " %" PetscInt_FMT ", Y range of indices: %" PetscInt_FMT " %" PetscInt_FMT "\n", info.xs, info.xs + info.xm, info.ys, info.ys + info.ym));
52 PetscCall(PetscViewerFlush(viewer));
53 PetscCall(PetscViewerASCIIPopSynchronized(viewer));
54 } else if (format == PETSC_VIEWER_ASCII_GLVIS) PetscCall(DMView_DA_GLVis(da, viewer));
55 } else if (isdraw) {
56 PetscDraw draw;
57 double ymin = -1 * dd->s - 1, ymax = dd->N + dd->s;
58 double xmin = -1 * dd->s - 1, xmax = dd->M + dd->s;
59 double x, y;
60 PetscInt base;
61 const PetscInt *idx;
62 char node[10];
63 PetscBool isnull;
64
65 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
66 PetscCall(PetscDrawIsNull(draw, &isnull));
67 if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
68
69 PetscCall(PetscDrawCheckResizedWindow(draw));
70 PetscCall(PetscDrawClear(draw));
71 PetscCall(PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax));
72
73 PetscDrawCollectiveBegin(draw);
74 /* first processor draw all node lines */
75 if (rank == 0) {
76 ymin = 0.0;
77 ymax = dd->N - 1;
78 for (xmin = 0; xmin < dd->M; xmin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_BLACK));
79 xmin = 0.0;
80 xmax = dd->M - 1;
81 for (ymin = 0; ymin < dd->N; ymin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK));
82 }
83 PetscDrawCollectiveEnd(draw);
84 PetscCall(PetscDrawFlush(draw));
85 PetscCall(PetscDrawPause(draw));
86
87 PetscDrawCollectiveBegin(draw);
88 /* draw my box */
89 xmin = dd->xs / dd->w;
90 xmax = (dd->xe - 1) / dd->w;
91 ymin = dd->ys;
92 ymax = dd->ye - 1;
93 PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED));
94 PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED));
95 PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED));
96 PetscCall(PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED));
97 /* put in numbers */
98 base = (dd->base) / dd->w;
99 for (y = ymin; y <= ymax; y++) {
100 for (x = xmin; x <= xmax; x++) {
101 PetscCall(PetscSNPrintf(node, sizeof(node), "%" PetscInt_FMT, base++));
102 PetscCall(PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node));
103 }
104 }
105 PetscDrawCollectiveEnd(draw);
106 PetscCall(PetscDrawFlush(draw));
107 PetscCall(PetscDrawPause(draw));
108
109 PetscDrawCollectiveBegin(draw);
110 /* overlay ghost numbers, useful for error checking */
111 PetscCall(ISLocalToGlobalMappingGetBlockIndices(da->ltogmap, &idx));
112 base = 0;
113 xmin = dd->Xs;
114 xmax = dd->Xe;
115 ymin = dd->Ys;
116 ymax = dd->Ye;
117 for (y = ymin; y < ymax; y++) {
118 for (x = xmin; x < xmax; x++) {
119 if ((base % dd->w) == 0) {
120 PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)(idx[base / dd->w])));
121 PetscCall(PetscDrawString(draw, x / dd->w, y, PETSC_DRAW_BLUE, node));
122 }
123 base++;
124 }
125 }
126 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx));
127 PetscDrawCollectiveEnd(draw);
128 PetscCall(PetscDrawFlush(draw));
129 PetscCall(PetscDrawPause(draw));
130 PetscCall(PetscDrawSave(draw));
131 } else if (isglvis) {
132 PetscCall(DMView_DA_GLVis(da, viewer));
133 } else if (isbinary) {
134 PetscCall(DMView_DA_Binary(da, viewer));
135 #if defined(PETSC_HAVE_MATLAB)
136 } else if (ismatlab) {
137 PetscCall(DMView_DA_Matlab(da, viewer));
138 #endif
139 }
140 PetscFunctionReturn(PETSC_SUCCESS);
141 }
142
143 #if defined(new)
144 /*
145 DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix-free matrix where local
146 function lives on a DMDA
147
148 y ~= (F(u + ha) - F(u))/h,
149 where F = nonlinear function, as set by SNESSetFunction()
150 u = current iterate
151 h = difference interval
152 */
DMDAGetDiagonal_MFFD(DM da,Vec U,Vec a)153 PetscErrorCode DMDAGetDiagonal_MFFD(DM da, Vec U, Vec a)
154 {
155 PetscScalar h, *aa, *ww, v;
156 PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
157 PetscInt gI, nI;
158 MatStencil stencil;
159 DMDALocalInfo info;
160
161 PetscFunctionBegin;
162 PetscCall((*ctx->func)(0, U, a, ctx->funcctx));
163 PetscCall((*ctx->funcisetbase)(U, ctx->funcctx));
164
165 PetscCall(VecGetArray(U, &ww));
166 PetscCall(VecGetArray(a, &aa));
167
168 nI = 0;
169 h = ww[gI];
170 if (h == 0.0) h = 1.0;
171 if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
172 else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
173 h *= epsilon;
174
175 ww[gI] += h;
176 PetscCall((*ctx->funci)(i, w, &v, ctx->funcctx));
177 aa[nI] = (v - aa[nI]) / h;
178 ww[gI] -= h;
179 nI++;
180
181 PetscCall(VecRestoreArray(U, &ww));
182 PetscCall(VecRestoreArray(a, &aa));
183 PetscFunctionReturn(PETSC_SUCCESS);
184 }
185 #endif
186
DMSetUp_DA_2D(DM da)187 PetscErrorCode DMSetUp_DA_2D(DM da)
188 {
189 DM_DA *dd = (DM_DA *)da->data;
190 const PetscInt M = dd->M;
191 const PetscInt N = dd->N;
192 PetscMPIInt m, n;
193 const PetscInt dof = dd->w;
194 const PetscInt s = dd->s;
195 DMBoundaryType bx = dd->bx;
196 DMBoundaryType by = dd->by;
197 DMDAStencilType stencil_type = dd->stencil_type;
198 PetscInt *lx = dd->lx;
199 PetscInt *ly = dd->ly;
200 MPI_Comm comm;
201 PetscMPIInt rank, size, n0, n1, n2, n3, n5, n6, n7, n8;
202 PetscInt xs, xe, ys, ye, x, y, Xs, Xe, Ys, Ye, IXs, IXe, IYs, IYe;
203 PetscInt up, down, left, right, i, *idx, nn;
204 PetscInt xbase, *bases, *ldims, j, x_t, y_t, s_t, base, count;
205 PetscInt s_x, s_y; /* s proportionalized to w */
206 PetscMPIInt sn0 = 0, sn2 = 0, sn6 = 0, sn8 = 0;
207 Vec local, global;
208 VecScatter gtol;
209 IS to, from;
210
211 PetscFunctionBegin;
212 PetscCheck(stencil_type != DMDA_STENCIL_BOX || (bx != DM_BOUNDARY_MIRROR && by != DM_BOUNDARY_MIRROR), PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Mirror boundary and box stencil");
213 PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
214 #if !defined(PETSC_USE_64BIT_INDICES)
215 PetscCheck(((PetscInt64)M) * ((PetscInt64)N) * ((PetscInt64)dof) <= (PetscInt64)PETSC_MPI_INT_MAX, comm, PETSC_ERR_INT_OVERFLOW, "Mesh of %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " (dof) is too large for 32-bit indices", M, N, dof);
216 #endif
217 PetscCall(PetscMPIIntCast(dd->m, &m));
218 PetscCall(PetscMPIIntCast(dd->n, &n));
219
220 PetscCallMPI(MPI_Comm_size(comm, &size));
221 PetscCallMPI(MPI_Comm_rank(comm, &rank));
222
223 dd->p = 1;
224 if (m != PETSC_DECIDE) {
225 PetscCheck(m >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in X direction: %d", m);
226 PetscCheck(m <= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in X direction: %d %d", m, size);
227 }
228 if (n != PETSC_DECIDE) {
229 PetscCheck(n >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Y direction: %d", n);
230 PetscCheck(n <= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Y direction: %d %d", n, size);
231 }
232
233 if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
234 if (n != PETSC_DECIDE) {
235 m = size / n;
236 } else if (m != PETSC_DECIDE) {
237 n = size / m;
238 } else {
239 /* try for squarish distribution */
240 m = (PetscMPIInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N)));
241 if (!m) m = 1;
242 while (m > 0) {
243 n = size / m;
244 if (m * n == size) break;
245 m--;
246 }
247 if (M > N && m < n) {
248 PetscMPIInt _m = m;
249 m = n;
250 n = _m;
251 }
252 }
253 PetscCheck(m * n == size, comm, PETSC_ERR_PLIB, "Unable to create partition, check the size of the communicator and input m and n ");
254 } else PetscCheck(m * n == size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Given Bad partition");
255
256 PetscCheck(M >= m, comm, PETSC_ERR_ARG_OUTOFRANGE, "Partition in x direction is too fine! %" PetscInt_FMT " %d", M, m);
257 PetscCheck(N >= n, comm, PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %d", N, n);
258
259 /*
260 Determine locally owned region
261 xs is the first local node number, x is the number of local nodes
262 */
263 if (!lx) {
264 PetscCall(PetscMalloc1(m, &dd->lx));
265 lx = dd->lx;
266 for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > i);
267 }
268 x = lx[rank % m];
269 xs = 0;
270 for (i = 0; i < (rank % m); i++) xs += lx[i];
271 if (PetscDefined(USE_DEBUG)) {
272 left = xs;
273 for (i = (rank % m); i < m; i++) left += lx[i];
274 PetscCheck(left == M, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Sum of lx across processors not equal to M: %" PetscInt_FMT " %" PetscInt_FMT, left, M);
275 }
276
277 /*
278 Determine locally owned region
279 ys is the first local node number, y is the number of local nodes
280 */
281 if (!ly) {
282 PetscCall(PetscMalloc1(n, &dd->ly));
283 ly = dd->ly;
284 for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > i);
285 }
286 y = ly[rank / m];
287 ys = 0;
288 for (i = 0; i < (rank / m); i++) ys += ly[i];
289 if (PetscDefined(USE_DEBUG)) {
290 left = ys;
291 for (i = (rank / m); i < n; i++) left += ly[i];
292 PetscCheck(left == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Sum of ly across processors not equal to N: %" PetscInt_FMT " %" PetscInt_FMT, left, N);
293 }
294
295 /*
296 check if the scatter requires more than one process neighbor or wraps around
297 the domain more than once
298 */
299 PetscCheck((x >= s) || ((m <= 1) && (bx != DM_BOUNDARY_PERIODIC)), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local x-width of domain x %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, x, s);
300 PetscCheck((y >= s) || ((n <= 1) && (by != DM_BOUNDARY_PERIODIC)), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local y-width of domain y %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, y, s);
301 PetscCheck((x > s) || (bx != DM_BOUNDARY_MIRROR), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local x-width of domain x %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT " with mirror", x, s);
302 PetscCheck((y > s) || (by != DM_BOUNDARY_MIRROR), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local y-width of domain y %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT " with mirror", y, s);
303 xe = xs + x;
304 ye = ys + y;
305
306 /* determine ghost region (Xs) and region scattered into (IXs) */
307 if (xs - s > 0) {
308 Xs = xs - s;
309 IXs = xs - s;
310 } else {
311 if (bx) {
312 Xs = xs - s;
313 } else {
314 Xs = 0;
315 }
316 IXs = 0;
317 }
318 if (xe + s <= M) {
319 Xe = xe + s;
320 IXe = xe + s;
321 } else {
322 if (bx) {
323 Xs = xs - s;
324 Xe = xe + s;
325 } else {
326 Xe = M;
327 }
328 IXe = M;
329 }
330
331 if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
332 IXs = xs - s;
333 IXe = xe + s;
334 Xs = xs - s;
335 Xe = xe + s;
336 }
337
338 if (ys - s > 0) {
339 Ys = ys - s;
340 IYs = ys - s;
341 } else {
342 if (by) {
343 Ys = ys - s;
344 } else {
345 Ys = 0;
346 }
347 IYs = 0;
348 }
349 if (ye + s <= N) {
350 Ye = ye + s;
351 IYe = ye + s;
352 } else {
353 if (by) {
354 Ye = ye + s;
355 } else {
356 Ye = N;
357 }
358 IYe = N;
359 }
360
361 if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
362 IYs = ys - s;
363 IYe = ye + s;
364 Ys = ys - s;
365 Ye = ye + s;
366 }
367
368 /* stencil length in each direction */
369 s_x = s;
370 s_y = s;
371
372 /* determine starting point of each processor */
373 nn = x * y;
374 PetscCall(PetscMalloc2(size + 1, &bases, size, &ldims));
375 PetscCallMPI(MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm));
376 bases[0] = 0;
377 for (i = 1; i <= size; i++) bases[i] = ldims[i - 1];
378 for (i = 1; i <= size; i++) bases[i] += bases[i - 1];
379 base = bases[rank] * dof;
380
381 /* allocate the base parallel and sequential vectors */
382 dd->Nlocal = x * y * dof;
383 PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global));
384 dd->nlocal = (Xe - Xs) * (Ye - Ys) * dof;
385 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local));
386
387 /* generate global to local vector scatter and local to global mapping*/
388
389 /* global to local must include ghost points within the domain,
390 but not ghost points outside the domain that aren't periodic */
391 PetscCall(PetscMalloc1((IXe - IXs) * (IYe - IYs), &idx));
392 if (stencil_type == DMDA_STENCIL_BOX) {
393 left = IXs - Xs;
394 right = left + (IXe - IXs);
395 down = IYs - Ys;
396 up = down + (IYe - IYs);
397 count = 0;
398 for (i = down; i < up; i++) {
399 for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs);
400 }
401 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));
402
403 } else {
404 /* must drop into cross shape region */
405 /* ---------|
406 | top |
407 |--- ---| up
408 | middle |
409 | |
410 ---- ---- down
411 | bottom |
412 -----------
413 Xs xs xe Xe */
414 left = xs - Xs;
415 right = left + x;
416 down = ys - Ys;
417 up = down + y;
418 count = 0;
419 /* bottom */
420 for (i = (IYs - Ys); i < down; i++) {
421 for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs);
422 }
423 /* middle */
424 for (i = down; i < up; i++) {
425 for (j = (IXs - Xs); j < (IXe - Xs); j++) idx[count++] = j + i * (Xe - Xs);
426 }
427 /* top */
428 for (i = up; i < up + IYe - ye; i++) {
429 for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs);
430 }
431 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to));
432 }
433
434 /* determine who lies on each side of us stored in n6 n7 n8
435 n3 n5
436 n0 n1 n2
437 */
438
439 /* Assume the Non-Periodic Case */
440 n1 = rank - m;
441 if (rank % m) {
442 n0 = n1 - 1;
443 } else {
444 n0 = -1;
445 }
446 if ((rank + 1) % m) {
447 n2 = n1 + 1;
448 n5 = rank + 1;
449 n8 = rank + m + 1;
450 if (n8 >= m * n) n8 = -1;
451 } else {
452 n2 = -1;
453 n5 = -1;
454 n8 = -1;
455 }
456 if (rank % m) {
457 n3 = rank - 1;
458 n6 = n3 + m;
459 if (n6 >= m * n) n6 = -1;
460 } else {
461 n3 = -1;
462 n6 = -1;
463 }
464 n7 = rank + m;
465 if (n7 >= m * n) n7 = -1;
466
467 if (bx == DM_BOUNDARY_PERIODIC && by == DM_BOUNDARY_PERIODIC) {
468 /* Modify for Periodic Cases */
469 /* Handle all four corners */
470 if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m - 1;
471 if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
472 if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size - m;
473 if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size - 1;
474
475 /* Handle Top and Bottom Sides */
476 if (n1 < 0) n1 = rank + m * (n - 1);
477 if (n7 < 0) n7 = rank - m * (n - 1);
478 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
479 if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1;
480 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
481 if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1;
482
483 /* Handle Left and Right Sides */
484 if (n3 < 0) n3 = rank + (m - 1);
485 if (n5 < 0) n5 = rank - (m - 1);
486 if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1;
487 if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1;
488 if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1;
489 if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1;
490 } else if (by == DM_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */
491 if (n1 < 0) n1 = rank + m * (n - 1);
492 if (n7 < 0) n7 = rank - m * (n - 1);
493 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
494 if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1;
495 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
496 if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1;
497 } else if (bx == DM_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */
498 if (n3 < 0) n3 = rank + (m - 1);
499 if (n5 < 0) n5 = rank - (m - 1);
500 if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1;
501 if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1;
502 if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1;
503 if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1;
504 }
505
506 PetscCall(PetscMalloc1(9, &dd->neighbors));
507
508 dd->neighbors[0] = n0;
509 dd->neighbors[1] = n1;
510 dd->neighbors[2] = n2;
511 dd->neighbors[3] = n3;
512 dd->neighbors[4] = rank;
513 dd->neighbors[5] = n5;
514 dd->neighbors[6] = n6;
515 dd->neighbors[7] = n7;
516 dd->neighbors[8] = n8;
517
518 if (stencil_type == DMDA_STENCIL_STAR) {
519 /* save corner processor numbers */
520 sn0 = n0;
521 sn2 = n2;
522 sn6 = n6;
523 sn8 = n8;
524 n0 = n2 = n6 = n8 = -1;
525 }
526
527 PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys), &idx));
528
529 nn = 0;
530 xbase = bases[rank];
531 for (i = 1; i <= s_y; i++) {
532 if (n0 >= 0) { /* left below */
533 x_t = lx[n0 % m];
534 y_t = ly[n0 / m];
535 s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x;
536 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
537 }
538
539 if (n1 >= 0) { /* directly below */
540 x_t = x;
541 y_t = ly[n1 / m];
542 s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t;
543 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
544 } else if (by == DM_BOUNDARY_MIRROR) {
545 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j;
546 }
547
548 if (n2 >= 0) { /* right below */
549 x_t = lx[n2 % m];
550 y_t = ly[n2 / m];
551 s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t;
552 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
553 }
554 }
555
556 for (i = 0; i < y; i++) {
557 if (n3 >= 0) { /* directly left */
558 x_t = lx[n3 % m];
559 /* y_t = y; */
560 s_t = bases[n3] + (i + 1) * x_t - s_x;
561 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
562 } else if (bx == DM_BOUNDARY_MIRROR) {
563 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j;
564 }
565
566 for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */
567
568 if (n5 >= 0) { /* directly right */
569 x_t = lx[n5 % m];
570 /* y_t = y; */
571 s_t = bases[n5] + (i)*x_t;
572 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
573 } else if (bx == DM_BOUNDARY_MIRROR) {
574 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j;
575 }
576 }
577
578 for (i = 1; i <= s_y; i++) {
579 if (n6 >= 0) { /* left above */
580 x_t = lx[n6 % m];
581 /* y_t = ly[n6 / m]; */
582 s_t = bases[n6] + (i)*x_t - s_x;
583 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
584 }
585
586 if (n7 >= 0) { /* directly above */
587 x_t = x;
588 /* y_t = ly[n7 / m]; */
589 s_t = bases[n7] + (i - 1) * x_t;
590 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
591 } else if (by == DM_BOUNDARY_MIRROR) {
592 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j;
593 }
594
595 if (n8 >= 0) { /* right above */
596 x_t = lx[n8 % m];
597 /* y_t = ly[n8 / m]; */
598 s_t = bases[n8] + (i - 1) * x_t;
599 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
600 }
601 }
602
603 PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from));
604 PetscCall(VecScatterCreate(global, from, local, to, >ol));
605 PetscCall(ISDestroy(&to));
606 PetscCall(ISDestroy(&from));
607
608 if (stencil_type == DMDA_STENCIL_STAR) {
609 n0 = sn0;
610 n2 = sn2;
611 n6 = sn6;
612 n8 = sn8;
613 }
614
615 if ((stencil_type == DMDA_STENCIL_STAR) || (bx && bx != DM_BOUNDARY_PERIODIC) || (by && by != DM_BOUNDARY_PERIODIC)) {
616 /*
617 Recompute the local to global mappings, this time keeping the
618 information about the cross corner processor numbers and any ghosted
619 but not periodic indices.
620 */
621 nn = 0;
622 xbase = bases[rank];
623 for (i = 1; i <= s_y; i++) {
624 if (n0 >= 0) { /* left below */
625 x_t = lx[n0 % m];
626 y_t = ly[n0 / m];
627 s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x;
628 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
629 } else if (xs - Xs > 0 && ys - Ys > 0) {
630 for (j = 0; j < s_x; j++) idx[nn++] = -1;
631 }
632 if (n1 >= 0) { /* directly below */
633 x_t = x;
634 y_t = ly[n1 / m];
635 s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t;
636 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
637 } else if (ys - Ys > 0) {
638 if (by == DM_BOUNDARY_MIRROR) {
639 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j;
640 } else {
641 for (j = 0; j < x; j++) idx[nn++] = -1;
642 }
643 }
644 if (n2 >= 0) { /* right below */
645 x_t = lx[n2 % m];
646 y_t = ly[n2 / m];
647 s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t;
648 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
649 } else if (Xe - xe > 0 && ys - Ys > 0) {
650 for (j = 0; j < s_x; j++) idx[nn++] = -1;
651 }
652 }
653
654 for (i = 0; i < y; i++) {
655 if (n3 >= 0) { /* directly left */
656 x_t = lx[n3 % m];
657 /* y_t = y; */
658 s_t = bases[n3] + (i + 1) * x_t - s_x;
659 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
660 } else if (xs - Xs > 0) {
661 if (bx == DM_BOUNDARY_MIRROR) {
662 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j;
663 } else {
664 for (j = 0; j < s_x; j++) idx[nn++] = -1;
665 }
666 }
667
668 for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */
669
670 if (n5 >= 0) { /* directly right */
671 x_t = lx[n5 % m];
672 /* y_t = y; */
673 s_t = bases[n5] + (i)*x_t;
674 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
675 } else if (Xe - xe > 0) {
676 if (bx == DM_BOUNDARY_MIRROR) {
677 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j;
678 } else {
679 for (j = 0; j < s_x; j++) idx[nn++] = -1;
680 }
681 }
682 }
683
684 for (i = 1; i <= s_y; i++) {
685 if (n6 >= 0) { /* left above */
686 x_t = lx[n6 % m];
687 /* y_t = ly[n6 / m]; */
688 s_t = bases[n6] + (i)*x_t - s_x;
689 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
690 } else if (xs - Xs > 0 && Ye - ye > 0) {
691 for (j = 0; j < s_x; j++) idx[nn++] = -1;
692 }
693 if (n7 >= 0) { /* directly above */
694 x_t = x;
695 /* y_t = ly[n7 / m]; */
696 s_t = bases[n7] + (i - 1) * x_t;
697 for (j = 0; j < x_t; j++) idx[nn++] = s_t++;
698 } else if (Ye - ye > 0) {
699 if (by == DM_BOUNDARY_MIRROR) {
700 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j;
701 } else {
702 for (j = 0; j < x; j++) idx[nn++] = -1;
703 }
704 }
705 if (n8 >= 0) { /* right above */
706 x_t = lx[n8 % m];
707 /* y_t = ly[n8 / m]; */
708 s_t = bases[n8] + (i - 1) * x_t;
709 for (j = 0; j < s_x; j++) idx[nn++] = s_t++;
710 } else if (Xe - xe > 0 && Ye - ye > 0) {
711 for (j = 0; j < s_x; j++) idx[nn++] = -1;
712 }
713 }
714 }
715 /*
716 Set the local to global ordering in the global vector, this allows use
717 of VecSetValuesLocal().
718 */
719 PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap));
720
721 PetscCall(PetscFree2(bases, ldims));
722 dd->m = m;
723 dd->n = n;
724 /* note PETSc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
725 dd->xs = xs * dof;
726 dd->xe = xe * dof;
727 dd->ys = ys;
728 dd->ye = ye;
729 dd->zs = 0;
730 dd->ze = 1;
731 dd->Xs = Xs * dof;
732 dd->Xe = Xe * dof;
733 dd->Ys = Ys;
734 dd->Ye = Ye;
735 dd->Zs = 0;
736 dd->Ze = 1;
737
738 PetscCall(VecDestroy(&local));
739 PetscCall(VecDestroy(&global));
740
741 dd->gtol = gtol;
742 dd->base = base;
743 da->ops->view = DMView_DA_2d;
744 dd->ltol = NULL;
745 dd->ao = NULL;
746 PetscFunctionReturn(PETSC_SUCCESS);
747 }
748
749 /*@
750 DMDACreate2d - Creates an object that will manage the communication of two-dimensional
751 regular array data that is distributed across one or more MPI processes.
752
753 Collective
754
755 Input Parameters:
756 + comm - MPI communicator
757 . bx - type of ghost nodes the x array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
758 . by - type of ghost nodes the y array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`.
759 . stencil_type - stencil type. Use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`.
760 . M - global dimension in x direction of the array
761 . N - global dimension in y direction of the array
762 . m - corresponding number of processors in x dimension (or `PETSC_DECIDE` to have calculated)
763 . n - corresponding number of processors in y dimension (or `PETSC_DECIDE` to have calculated)
764 . dof - number of degrees of freedom per node
765 . s - stencil width
766 . lx - arrays containing the number of nodes in each cell along the x coordinates, or `NULL`.
767 - ly - arrays containing the number of nodes in each cell along the y coordinates, or `NULL`.
768
769 Output Parameter:
770 . da - the resulting distributed array object
771
772 Options Database Keys:
773 + -dm_view - Calls `DMView()` at the conclusion of `DMDACreate2d()`
774 . -da_grid_x <nx> - number of grid points in x direction
775 . -da_grid_y <ny> - number of grid points in y direction
776 . -da_processors_x <nx> - number of processors in x direction
777 . -da_processors_y <ny> - number of processors in y direction
778 . -da_bd_x <bx> - boundary type in x direction
779 . -da_bd_y <by> - boundary type in y direction
780 . -da_bd_all <bt> - boundary type in all directions
781 . -da_refine_x <rx> - refinement ratio in x direction
782 . -da_refine_y <ry> - refinement ratio in y direction
783 - -da_refine <n> - refine the `DMDA` n times before creating
784
785 Level: beginner
786
787 Notes:
788 If `lx` or `ly` are non-null, these must be of length as `m` and `n`, and the corresponding
789 `m` and `n` cannot be `PETSC_DECIDE`. The sum of the `lx` entries must be `M`, and the sum of
790 the `ly` entries must be `N`.
791
792 The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the
793 standard 5-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes
794 the standard 9-pt stencil.
795
796 The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects;
797 The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()`
798 and DMCreateLocalVector() and calls to `VecDuplicate()` if more are needed.
799
800 You must call `DMSetUp()` after this call before using this `DM`.
801
802 To use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call
803 but before `DMSetUp()`.
804
805 .seealso: [](sec_struct), `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate3d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`,
806 `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`,
807 `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`,
808 `DMStagCreate2d()`, `DMBoundaryType`
809 @*/
DMDACreate2d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMDAStencilType stencil_type,PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM * da)810 PetscErrorCode DMDACreate2d(MPI_Comm comm, DMBoundaryType bx, DMBoundaryType by, DMDAStencilType stencil_type, PetscInt M, PetscInt N, PetscInt m, PetscInt n, PetscInt dof, PetscInt s, const PetscInt lx[], const PetscInt ly[], DM *da)
811 {
812 PetscFunctionBegin;
813 PetscCall(DMDACreate(comm, da));
814 PetscCall(DMSetDimension(*da, 2));
815 PetscCall(DMDASetSizes(*da, M, N, 1));
816 PetscCall(DMDASetNumProcs(*da, m, n, PETSC_DECIDE));
817 PetscCall(DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE));
818 PetscCall(DMDASetDof(*da, dof));
819 PetscCall(DMDASetStencilType(*da, stencil_type));
820 PetscCall(DMDASetStencilWidth(*da, s));
821 PetscCall(DMDASetOwnershipRanges(*da, lx, ly, NULL));
822 PetscFunctionReturn(PETSC_SUCCESS);
823 }
824