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