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