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   isascii, 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, &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;

     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;

     65:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
     66:     PetscCall(PetscDrawIsNull(draw, &isnull));
     67:     if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

     69:     PetscCall(PetscDrawCheckResizedWindow(draw));
     70:     PetscCall(PetscDrawClear(draw));
     71:     PetscCall(PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax));

     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));

     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));

    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: }

    143: #if defined(new)
    144: /*
    145:   DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix-free matrix where local
    146:     function lives on a DMDA

    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: */
    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;

    161:   PetscFunctionBegin;
    162:   PetscCall((*ctx->func)(0, U, a, ctx->funcctx));
    163:   PetscCall((*ctx->funcisetbase)(U, ctx->funcctx));

    165:   PetscCall(VecGetArray(U, &ww));
    166:   PetscCall(VecGetArray(a, &aa));

    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;

    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++;

    181:   PetscCall(VecRestoreArray(U, &ww));
    182:   PetscCall(VecRestoreArray(a, &aa));
    183:   PetscFunctionReturn(PETSC_SUCCESS);
    184: }
    185: #endif

    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;

    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));

    220:   PetscCallMPI(MPI_Comm_size(comm, &size));
    221:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    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:   }

    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");

    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);

    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:   }

    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:   }

    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;

    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:   }

    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:   }

    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:   }

    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:   }

    368:   /* stencil length in each direction */
    369:   s_x = s;
    370:   s_y = s;

    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;

    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));

    387:   /* generate global to local vector scatter and local to global mapping*/

    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));

    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:   }

    434:   /* determine who lies on each side of us stored in    n6 n7 n8
    435:                                                         n3    n5
    436:                                                         n0 n1 n2
    437:   */

    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;

    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;

    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;

    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:   }

    506:   PetscCall(PetscMalloc1(9, &dd->neighbors));

    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;

    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:   }

    527:   PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys), &idx));

    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:     }

    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:     }

    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:   }

    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:     }

    566:     for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */

    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:   }

    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:     }

    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:     }

    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:   }

    603:   PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from));
    604:   PetscCall(VecScatterCreate(global, from, local, to, &gtol));
    605:   PetscCall(ISDestroy(&to));
    606:   PetscCall(ISDestroy(&from));

    608:   if (stencil_type == DMDA_STENCIL_STAR) {
    609:     n0 = sn0;
    610:     n2 = sn2;
    611:     n6 = sn6;
    612:     n8 = sn8;
    613:   }

    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:     }

    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:       }

    668:       for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */

    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:     }

    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));

    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;

    738:   PetscCall(VecDestroy(&local));
    739:   PetscCall(VecDestroy(&global));

    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: }

    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.

    753:   Collective

    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 the x direction of the array
    761: . N            - global dimension in the y direction of the array
    762: . m            - corresponding number of processors in the x dimension (or `PETSC_DECIDE` to have calculated)
    763: . n            - corresponding number of processors in the 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`.

    769:   Output Parameter:
    770: . da - the resulting distributed array object

    772:   Options Database Keys:
    773: + -dm_view            - Calls `DMView()` at the conclusion of `DMDACreate2d()`
    774: . -da_grid_x nx       - number of grid points in the x direction
    775: . -da_grid_y ny       - number of grid points in the y direction
    776: . -da_processors_x nx - number of processors in the x direction
    777: . -da_processors_y ny - number of processors in the y direction
    778: . -da_bd_x bx         - boundary type in the x direction
    779: . -da_bd_y by         - boundary type in the y direction
    780: . -da_bd_all bt       - boundary type in all directions
    781: . -da_refine_x rx     - refinement ratio in the x direction
    782: . -da_refine_y ry     - refinement ratio in the y direction
    783: - -da_refine n        - refine the `DMDA` `n` times before creating

    785:   Level: beginner

    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`.

    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.

    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.

    800:   You must call `DMSetUp()` after this call before using this `DM`.

    802:   To use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call
    803:   but before `DMSetUp()`.

    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: @*/
    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: }