Actual source code: gr1.c

      1: /*
      2:    Plots vectors obtained with DMDACreate1d()
      3: */

      5: #include <petsc/private/dmdaimpl.h>

      7: /*@
      8:   DMDASetUniformCoordinates - Sets a `DMDA` coordinates to be a uniform grid

     10:   Collective

     12:   Input Parameters:
     13: + da   - the `DMDA` object
     14: . xmin - min extreme in the x direction
     15: . xmax - max extreme in the x direction
     16: . ymin - min extreme in the y direction (value ignored for 1 dimensional problems)
     17: . ymax - max extreme in the y direction (value ignored for 1 dimensional problems)
     18: . zmin - min extreme in the z direction (value ignored for 1 or 2 dimensional problems)
     19: - zmax - max extreme in the z direction (value ignored for 1 or 2 dimensional problems)

     21:   Level: beginner

     23: .seealso: [](sec_struct), `DM`, `DMDA`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMStagSetUniformCoordinates()`
     24: @*/
     25: PetscErrorCode DMDASetUniformCoordinates(DM da, PetscReal xmin, PetscReal xmax, PetscReal ymin, PetscReal ymax, PetscReal zmin, PetscReal zmax)
     26: {
     27:   MPI_Comm       comm;
     28:   DM             cda;
     29:   DM_DA         *dd = (DM_DA *)da->data;
     30:   DMBoundaryType bx, by, bz;
     31:   Vec            xcoor;
     32:   PetscScalar   *coors;
     33:   PetscReal      hx = 0., hy = 0., hz_ = 0.;
     34:   PetscInt       i, j, k, M, N, P, istart, isize, jstart, jsize, kstart, ksize, dim, cnt;

     36:   PetscFunctionBegin;
     38:   PetscCheck(dd->gtol, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Cannot set coordinates until after DMDA has been setup");
     39:   PetscCall(DMDAGetInfo(da, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, &bx, &by, &bz, NULL));
     40:   PetscCheck(xmax >= xmin, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_INCOMP, "xmax must be larger than xmin %g %g", (double)xmin, (double)xmax);
     41:   PetscCheck(!(dim > 1) || !(ymax < ymin), PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_INCOMP, "ymax must be larger than ymin %g %g", (double)ymin, (double)ymax);
     42:   PetscCheck(!(dim > 2) || !(zmax < zmin), PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_INCOMP, "zmax must be larger than zmin %g %g", (double)zmin, (double)zmax);
     43:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
     44:   PetscCall(DMDAGetCorners(da, &istart, &jstart, &kstart, &isize, &jsize, &ksize));
     45:   PetscCall(DMGetCoordinateDM(da, &cda));
     46:   PetscCall(DMCreateGlobalVector(cda, &xcoor));
     47:   if (dim == 1) {
     48:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax - xmin) / M;
     49:     else hx = (xmax - xmin) / (M - 1);
     50:     PetscCall(VecGetArray(xcoor, &coors));
     51:     for (i = 0; i < isize; i++) coors[i] = xmin + hx * (i + istart);
     52:     PetscCall(VecRestoreArray(xcoor, &coors));
     53:   } else if (dim == 2) {
     54:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax - xmin) / (M);
     55:     else hx = (xmax - xmin) / (M - 1);
     56:     if (by == DM_BOUNDARY_PERIODIC) hy = (ymax - ymin) / (N);
     57:     else hy = (ymax - ymin) / (N - 1);
     58:     PetscCall(VecGetArray(xcoor, &coors));
     59:     cnt = 0;
     60:     for (j = 0; j < jsize; j++) {
     61:       for (i = 0; i < isize; i++) {
     62:         coors[cnt++] = xmin + hx * (i + istart);
     63:         coors[cnt++] = ymin + hy * (j + jstart);
     64:       }
     65:     }
     66:     PetscCall(VecRestoreArray(xcoor, &coors));
     67:   } else if (dim == 3) {
     68:     if (bx == DM_BOUNDARY_PERIODIC) hx = (xmax - xmin) / (M);
     69:     else hx = (xmax - xmin) / (M - 1);
     70:     if (by == DM_BOUNDARY_PERIODIC) hy = (ymax - ymin) / (N);
     71:     else hy = (ymax - ymin) / (N - 1);
     72:     if (bz == DM_BOUNDARY_PERIODIC) hz_ = (zmax - zmin) / (P);
     73:     else hz_ = (zmax - zmin) / (P - 1);
     74:     PetscCall(VecGetArray(xcoor, &coors));
     75:     cnt = 0;
     76:     for (k = 0; k < ksize; k++) {
     77:       for (j = 0; j < jsize; j++) {
     78:         for (i = 0; i < isize; i++) {
     79:           coors[cnt++] = xmin + hx * (i + istart);
     80:           coors[cnt++] = ymin + hy * (j + jstart);
     81:           coors[cnt++] = zmin + hz_ * (k + kstart);
     82:         }
     83:       }
     84:     }
     85:     PetscCall(VecRestoreArray(xcoor, &coors));
     86:   } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Cannot create uniform coordinates for this dimension %" PetscInt_FMT, dim);
     87:   PetscCall(DMSetCoordinates(da, xcoor));
     88:   PetscCall(VecDestroy(&xcoor));
     89:   // Handle periodicity
     90:   if (bx == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_PERIODIC || bz == DM_BOUNDARY_PERIODIC) {
     91:     PetscReal maxCell[3] = {hx, hy, hz_};
     92:     PetscReal Lstart[3]  = {xmin, ymin, zmin};
     93:     PetscReal L[3]       = {xmax - xmin, ymax - ymin, zmax - zmin};

     95:     PetscCall(DMSetPeriodicity(da, maxCell, Lstart, L));
     96:   }
     97:   PetscFunctionReturn(PETSC_SUCCESS);
     98: }

    100: /*
    101:     Allows a user to select a subset of the fields to be drawn by VecView() when the vector comes from a DMDA
    102: */
    103: PetscErrorCode DMDASelectFields(DM da, PetscInt *outfields, PetscInt **fields)
    104: {
    105:   PetscInt  step, ndisplayfields, *displayfields, k, j;
    106:   PetscBool flg;

    108:   PetscFunctionBegin;
    109:   PetscCall(DMDAGetInfo(da, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &step, NULL, NULL, NULL, NULL, NULL));
    110:   PetscCall(PetscMalloc1(step, &displayfields));
    111:   for (k = 0; k < step; k++) displayfields[k] = k;
    112:   ndisplayfields = step;
    113:   PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-draw_fields", displayfields, &ndisplayfields, &flg));
    114:   if (!ndisplayfields) ndisplayfields = step;
    115:   if (!flg) {
    116:     char      **fields;
    117:     const char *fieldname;
    118:     PetscInt    nfields = step;
    119:     PetscCall(PetscMalloc1(step, &fields));
    120:     PetscCall(PetscOptionsGetStringArray(NULL, NULL, "-draw_fields_by_name", fields, &nfields, &flg));
    121:     if (flg) {
    122:       ndisplayfields = 0;
    123:       for (k = 0; k < nfields; k++) {
    124:         for (j = 0; j < step; j++) {
    125:           PetscCall(DMDAGetFieldName(da, j, &fieldname));
    126:           PetscCall(PetscStrcmp(fieldname, fields[k], &flg));
    127:           if (flg) goto found;
    128:         }
    129:         SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Unknown fieldname %s", fields[k]);
    130:       found:
    131:         displayfields[ndisplayfields++] = j;
    132:       }
    133:     }
    134:     for (k = 0; k < nfields; k++) PetscCall(PetscFree(fields[k]));
    135:     PetscCall(PetscFree(fields));
    136:   }
    137:   *fields    = displayfields;
    138:   *outfields = ndisplayfields;
    139:   PetscFunctionReturn(PETSC_SUCCESS);
    140: }

    142: #include <petscdraw.h>

    144: PetscErrorCode VecView_MPI_Draw_DA1d(Vec xin, PetscViewer v)
    145: {
    146:   DM                  da;
    147:   PetscMPIInt         rank, size, tag;
    148:   PetscInt            i, n, N, dof, istart, isize, j, nbounds;
    149:   MPI_Status          status;
    150:   PetscReal           min, max, xmin = 0.0, xmax = 0.0, tmp = 0.0, xgtmp = 0.0;
    151:   const PetscScalar  *array, *xg;
    152:   PetscDraw           draw;
    153:   PetscBool           isnull, useports = PETSC_FALSE, showmarkers = PETSC_FALSE;
    154:   MPI_Comm            comm;
    155:   PetscDrawAxis       axis;
    156:   Vec                 xcoor;
    157:   DMBoundaryType      bx;
    158:   const char         *tlabel = NULL, *xlabel = NULL;
    159:   const PetscReal    *bounds;
    160:   PetscInt           *displayfields;
    161:   PetscInt            k, ndisplayfields;
    162:   PetscBool           hold;
    163:   PetscDrawViewPorts *ports = NULL;
    164:   PetscViewerFormat   format;

    166:   PetscFunctionBegin;
    167:   PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
    168:   PetscCall(PetscDrawIsNull(draw, &isnull));
    169:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
    170:   PetscCall(PetscViewerDrawGetBounds(v, &nbounds, &bounds));

    172:   PetscCall(VecGetDM(xin, &da));
    173:   PetscCheck(da, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA");
    174:   PetscCall(PetscObjectGetComm((PetscObject)xin, &comm));
    175:   PetscCallMPI(MPI_Comm_size(comm, &size));
    176:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    178:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_vec_use_markers", &showmarkers, NULL));

    180:   PetscCall(DMDAGetInfo(da, NULL, &N, NULL, NULL, NULL, NULL, NULL, &dof, NULL, &bx, NULL, NULL, NULL));
    181:   PetscCall(DMDAGetCorners(da, &istart, NULL, NULL, &isize, NULL, NULL));
    182:   PetscCall(VecGetArrayRead(xin, &array));
    183:   PetscCall(VecGetLocalSize(xin, &n));
    184:   n = n / dof;

    186:   /* Get coordinates of nodes */
    187:   PetscCall(DMGetCoordinates(da, &xcoor));
    188:   if (!xcoor) {
    189:     PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0));
    190:     PetscCall(DMGetCoordinates(da, &xcoor));
    191:   }
    192:   PetscCall(VecGetArrayRead(xcoor, &xg));
    193:   PetscCall(DMDAGetCoordinateName(da, 0, &xlabel));

    195:   /* Determine the min and max coordinate in plot */
    196:   if (rank == 0) xmin = PetscRealPart(xg[0]);
    197:   if (rank == size - 1) xmax = PetscRealPart(xg[n - 1]);
    198:   PetscCallMPI(MPI_Bcast(&xmin, 1, MPIU_REAL, 0, comm));
    199:   PetscCallMPI(MPI_Bcast(&xmax, 1, MPIU_REAL, size - 1, comm));

    201:   PetscCall(DMDASelectFields(da, &ndisplayfields, &displayfields));
    202:   PetscCall(PetscViewerGetFormat(v, &format));
    203:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_ports", &useports, NULL));
    204:   if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
    205:   if (useports) {
    206:     PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
    207:     PetscCall(PetscViewerDrawGetDrawAxis(v, 0, &axis));
    208:     PetscCall(PetscDrawCheckResizedWindow(draw));
    209:     PetscCall(PetscDrawClear(draw));
    210:     PetscCall(PetscDrawViewPortsCreate(draw, ndisplayfields, &ports));
    211:   }

    213:   /* Loop over each field; drawing each in a different window */
    214:   for (k = 0; k < ndisplayfields; k++) {
    215:     j = displayfields[k];

    217:     /* determine the min and max value in plot */
    218:     PetscCall(VecStrideMin(xin, j, NULL, &min));
    219:     PetscCall(VecStrideMax(xin, j, NULL, &max));
    220:     if (j < nbounds) {
    221:       min = PetscMin(min, bounds[2 * j]);
    222:       max = PetscMax(max, bounds[2 * j + 1]);
    223:     }
    224:     if (min == max) {
    225:       min -= 1.e-5;
    226:       max += 1.e-5;
    227:     }

    229:     if (useports) {
    230:       PetscCall(PetscDrawViewPortsSet(ports, k));
    231:       PetscCall(DMDAGetFieldName(da, j, &tlabel));
    232:     } else {
    233:       const char *title;
    234:       PetscCall(PetscViewerDrawGetHold(v, &hold));
    235:       PetscCall(PetscViewerDrawGetDraw(v, k, &draw));
    236:       PetscCall(PetscViewerDrawGetDrawAxis(v, k, &axis));
    237:       PetscCall(DMDAGetFieldName(da, j, &title));
    238:       if (title) PetscCall(PetscDrawSetTitle(draw, title));
    239:       PetscCall(PetscDrawCheckResizedWindow(draw));
    240:       if (!hold) PetscCall(PetscDrawClear(draw));
    241:     }
    242:     PetscCall(PetscDrawAxisSetLabels(axis, tlabel, xlabel, NULL));
    243:     PetscCall(PetscDrawAxisSetLimits(axis, xmin, xmax, min, max));
    244:     PetscCall(PetscDrawAxisDraw(axis));

    246:     /* draw local part of vector */
    247:     PetscCall(PetscObjectGetNewTag((PetscObject)xin, &tag));
    248:     if (rank < size - 1) { /*send value to right */
    249:       PetscCallMPI(MPI_Send((void *)&xg[n - 1], 1, MPIU_REAL, rank + 1, tag, comm));
    250:       PetscCallMPI(MPI_Send((void *)&array[j + (n - 1) * dof], 1, MPIU_REAL, rank + 1, tag, comm));
    251:     }
    252:     if (rank) { /* receive value from left */
    253:       PetscCallMPI(MPI_Recv(&xgtmp, 1, MPIU_REAL, rank - 1, tag, comm, &status));
    254:       PetscCallMPI(MPI_Recv(&tmp, 1, MPIU_REAL, rank - 1, tag, comm, &status));
    255:     }
    256:     PetscDrawCollectiveBegin(draw);
    257:     if (rank) {
    258:       PetscCall(PetscDrawLine(draw, xgtmp, tmp, PetscRealPart(xg[0]), PetscRealPart(array[j]), PETSC_DRAW_RED));
    259:       if (showmarkers) PetscCall(PetscDrawPoint(draw, xgtmp, tmp, PETSC_DRAW_BLACK));
    260:     }
    261:     for (i = 1; i < n; i++) {
    262:       PetscCall(PetscDrawLine(draw, PetscRealPart(xg[i - 1]), PetscRealPart(array[j + dof * (i - 1)]), PetscRealPart(xg[i]), PetscRealPart(array[j + dof * i]), PETSC_DRAW_RED));
    263:       if (showmarkers) PetscCall(PetscDrawMarker(draw, PetscRealPart(xg[i - 1]), PetscRealPart(array[j + dof * (i - 1)]), PETSC_DRAW_BLACK));
    264:     }
    265:     if (rank == size - 1) {
    266:       if (showmarkers) PetscCall(PetscDrawMarker(draw, PetscRealPart(xg[n - 1]), PetscRealPart(array[j + dof * (n - 1)]), PETSC_DRAW_BLACK));
    267:     }
    268:     PetscDrawCollectiveEnd(draw);
    269:     PetscCall(PetscDrawFlush(draw));
    270:     PetscCall(PetscDrawPause(draw));
    271:     if (!useports) PetscCall(PetscDrawSave(draw));
    272:   }
    273:   if (useports) {
    274:     PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
    275:     PetscCall(PetscDrawSave(draw));
    276:   }

    278:   PetscCall(PetscDrawViewPortsDestroy(ports));
    279:   PetscCall(PetscFree(displayfields));
    280:   PetscCall(VecRestoreArrayRead(xcoor, &xg));
    281:   PetscCall(VecRestoreArrayRead(xin, &array));
    282:   PetscFunctionReturn(PETSC_SUCCESS);
    283: }