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, &gtol));
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: }