Actual source code: gr2.c

      1: #define PETSCDM_DLL

      3: /*
      4:    Plots vectors obtained with DACreate2d()
      5: */

     7:  #include src/dm/da/daimpl.h
     8:  #include vecimpl.h

     10: #if defined(PETSC_HAVE_PNETCDF)
     12: #include "pnetcdf.h"
     14: #endif

     17: /*
     18:         The data that is passed into the graphics callback
     19: */
     20: typedef struct {
     21:   PetscInt     m,n,step,k;
     22:   PetscReal    min,max,scale;
     23:   PetscScalar  *xy,*v;
     24:   PetscTruth   showgrid;
     25: } ZoomCtx;

     27: /*
     28:        This does the drawing for one particular field
     29:     in one particular set of coordinates. It is a callback
     30:     called from PetscDrawZoom()
     31: */
     34: PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw,void *ctx)
     35: {
     36:   ZoomCtx        *zctx = (ZoomCtx*)ctx;
     38:   PetscInt       m,n,i,j,k,step,id,c1,c2,c3,c4;
     39:   PetscReal      s,min,x1,x2,x3,x4,y_1,y2,y3,y4;
     40:   PetscScalar   *v,*xy;

     43:   m    = zctx->m;
     44:   n    = zctx->n;
     45:   step = zctx->step;
     46:   k    = zctx->k;
     47:   v    = zctx->v;
     48:   xy   = zctx->xy;
     49:   s    = zctx->scale;
     50:   min  = zctx->min;
     51:
     52:   /* PetscDraw the contour plot patch */
     53:   for (j=0; j<n-1; j++) {
     54:     for (i=0; i<m-1; i++) {
     55: #if !defined(PETSC_USE_COMPLEX)
     56:       id = i+j*m;    x1 = xy[2*id];y_1 = xy[2*id+1];c1 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min));
     57:       id = i+j*m+1;  x2 = xy[2*id];y2  = y_1;       c2 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min));
     58:       id = i+j*m+1+m;x3 = x2;      y3  = xy[2*id+1];c3 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min));
     59:       id = i+j*m+m;  x4 = x1;      y4  = y3;        c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min));
     60: #else
     61:       id = i+j*m;    x1 = PetscRealPart(xy[2*id]);y_1 = PetscRealPart(xy[2*id+1]);c1 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min));
     62:       id = i+j*m+1;  x2 = PetscRealPart(xy[2*id]);y2  = y_1;       c2 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min));
     63:       id = i+j*m+1+m;x3 = x2;      y3  = PetscRealPart(xy[2*id+1]);c3 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min));
     64:       id = i+j*m+m;  x4 = x1;      y4  = y3;        c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min));
     65: #endif
     66:       PetscDrawTriangle(draw,x1,y_1,x2,y2,x3,y3,c1,c2,c3);
     67:       PetscDrawTriangle(draw,x1,y_1,x3,y3,x4,y4,c1,c3,c4);
     68:       if (zctx->showgrid) {
     69:         PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);
     70:         PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);
     71:         PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);
     72:         PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);
     73:       }
     74:     }
     75:   }
     76:   return(0);
     77: }

     81: PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer)
     82: {
     83:   DA                 da,dac,dag;
     84:   PetscErrorCode     ierr;
     85:   PetscMPIInt        rank;
     86:   PetscInt           igstart,N,s,M,istart,isize,jgstart,*lx,*ly,w;
     87:   PetscReal          coors[4],ymin,ymax,xmin,xmax;
     88:   PetscDraw          draw,popup;
     89:   PetscTruth         isnull,useports;
     90:   MPI_Comm           comm;
     91:   Vec                xlocal,xcoor,xcoorl;
     92:   DAPeriodicType     periodic;
     93:   DAStencilType      st;
     94:   ZoomCtx            zctx;
     95:   PetscDrawViewPorts *ports;
     96:   PetscViewerFormat  format;

     99:   PetscViewerDrawGetDraw(viewer,0,&draw);
    100:   PetscDrawIsNull(draw,&isnull); if (isnull) return(0);

    102:   PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);
    103:   if (!da) SETERRQ(PETSC_ERR_ARG_WRONG,"Vector not generated from a DA");

    105:   PetscObjectGetComm((PetscObject)xin,&comm);
    106:   MPI_Comm_rank(comm,&rank);

    108:   DAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&periodic,&st);
    109:   DAGetOwnershipRange(da,&lx,&ly,PETSC_NULL);

    111:   /*
    112:         Obtain a sequential vector that is going to contain the local values plus ONE layer of
    113:      ghosted values to draw the graphics from. We also need its corresponding DA (dac) that will
    114:      update the local values pluse ONE layer of ghost values.
    115:   */
    116:   PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);
    117:   if (!xlocal) {
    118:     if (periodic != DA_NONPERIODIC || s != 1 || st != DA_STENCIL_BOX) {
    119:       /*
    120:          if original da is not of stencil width one, or periodic or not a box stencil then
    121:          create a special DA to handle one level of ghost points for graphics
    122:       */
    123:       DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,M,N,zctx.m,zctx.n,w,1,lx,ly,&dac);
    124:       PetscLogInfo((da,"VecView_MPI_Draw_DA2d:Creating auxilary DA for managing graphics ghost points\n"));
    125:     } else {
    126:       /* otherwise we can use the da we already have */
    127:       dac = da;
    128:     }
    129:     /* create local vector for holding ghosted values used in graphics */
    130:     DACreateLocalVector(dac,&xlocal);
    131:     if (dac != da) {
    132:       /* don't keep any public reference of this DA, is is only available through xlocal */
    133:       DADestroy(dac);
    134:     } else {
    135:       /* remove association between xlocal and da, because below we compose in the opposite
    136:          direction and if we left this connect we'd get a loop, so the objects could
    137:          never be destroyed */
    138:       PetscObjectCompose((PetscObject)xlocal,"DA",0);
    139:     }
    140:     PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);
    141:     PetscObjectDereference((PetscObject)xlocal);
    142:   } else {
    143:     if (periodic == DA_NONPERIODIC && s == 1 && st == DA_STENCIL_BOX) {
    144:       dac = da;
    145:     } else {
    146:       PetscObjectQuery((PetscObject)xlocal,"DA",(PetscObject*)&dac);
    147:     }
    148:   }

    150:   /*
    151:       Get local (ghosted) values of vector
    152:   */
    153:   DAGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);
    154:   DAGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);
    155:   VecGetArray(xlocal,&zctx.v);

    157:   /* get coordinates of nodes */
    158:   DAGetCoordinates(da,&xcoor);
    159:   if (!xcoor) {
    160:     DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);
    161:     DAGetCoordinates(da,&xcoor);
    162:   }

    164:   /*
    165:       Determine the min and max  coordinates in plot
    166:   */
    167:   VecStrideMin(xcoor,0,PETSC_NULL,&xmin);
    168:   VecStrideMax(xcoor,0,PETSC_NULL,&xmax);
    169:   VecStrideMin(xcoor,1,PETSC_NULL,&ymin);
    170:   VecStrideMax(xcoor,1,PETSC_NULL,&ymax);
    171:   coors[0] = xmin - .05*(xmax- xmin); coors[2] = xmax + .05*(xmax - xmin);
    172:   coors[1] = ymin - .05*(ymax- ymin); coors[3] = ymax + .05*(ymax - ymin);
    173:   PetscLogInfo((da,"VecView_MPI_Draw_DA2d:Preparing DA 2d contour plot coordinates %g %g %g %g\n",coors[0],coors[1],coors[2],coors[3]));

    175:   /*
    176:        get local ghosted version of coordinates
    177:   */
    178:   PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);
    179:   if (!xcoorl) {
    180:     /* create DA to get local version of graphics */
    181:     DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,M,N,zctx.m,zctx.n,2,1,lx,ly,&dag);
    182:     PetscLogInfo((dag,"VecView_MPI_Draw_DA2d:Creating auxilary DA for managing graphics coordinates ghost points\n"));
    183:     DACreateLocalVector(dag,&xcoorl);
    184:     PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);
    185:     DADestroy(dag);/* dereference dag */
    186:     PetscObjectDereference((PetscObject)xcoorl);
    187:   } else {
    188:     PetscObjectQuery((PetscObject)xcoorl,"DA",(PetscObject*)&dag);
    189:   }
    190:   DAGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);
    191:   DAGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);
    192:   VecGetArray(xcoorl,&zctx.xy);
    193:
    194:   /*
    195:         Get information about size of area each processor must do graphics for
    196:   */
    197:   DAGetInfo(dac,0,&M,&N,0,0,0,0,&zctx.step,0,&periodic,0);
    198:   DAGetGhostCorners(dac,&igstart,&jgstart,0,&zctx.m,&zctx.n,0);
    199:   DAGetCorners(dac,&istart,0,0,&isize,0,0);

    201:   PetscOptionsHasName(PETSC_NULL,"-draw_contour_grid",&zctx.showgrid);

    203:   PetscViewerGetFormat(viewer,&format);
    204:   PetscOptionsHasName(PETSC_NULL,"-draw_ports",&useports);
    205:   if (useports || format == PETSC_VIEWER_DRAW_PORTS){
    206:     PetscDrawSynchronizedClear(draw);
    207:     PetscDrawViewPortsCreate(draw,zctx.step,&ports);
    208:   }
    209:   /*
    210:      Loop over each field; drawing each in a different window
    211:   */
    212:   for (zctx.k=0; zctx.k<zctx.step; zctx.k++) {
    213:     if (useports) {
    214:       PetscDrawViewPortsSet(ports,zctx.k);
    215:     } else {
    216:       PetscViewerDrawGetDraw(viewer,zctx.k,&draw);
    217:       PetscDrawSynchronizedClear(draw);
    218:     }

    220:     /*
    221:         Determine the min and max color in plot
    222:     */
    223:     VecStrideMin(xin,zctx.k,PETSC_NULL,&zctx.min);
    224:     VecStrideMax(xin,zctx.k,PETSC_NULL,&zctx.max);
    225:     if (zctx.min == zctx.max) {
    226:       zctx.min -= 1.e-12;
    227:       zctx.max += 1.e-12;
    228:     }

    230:     if (!rank) {
    231:       char *title;

    233:       DAGetFieldName(da,zctx.k,&title);
    234:       if (title) {
    235:         PetscDrawSetTitle(draw,title);
    236:       }
    237:     }
    238:     PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);
    239:     PetscLogInfo((da,"VecView_MPI_Draw_DA2d:DA 2d contour plot min %g max %g\n",zctx.min,zctx.max));

    241:     PetscDrawGetPopup(draw,&popup);
    242:     if (popup) {PetscDrawScalePopup(popup,zctx.min,zctx.max);}

    244:     zctx.scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/(zctx.max - zctx.min);

    246:     PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);
    247:   }
    248:   if (useports){
    249:     PetscDrawViewPortsDestroy(ports);
    250:   }

    252:   VecRestoreArray(xcoorl,&zctx.xy);
    253:   VecRestoreArray(xlocal,&zctx.v);
    254:   return(0);
    255: }

    257: EXTERN PetscErrorCode VecView_MPI_HDF4_Ex(Vec X, PetscViewer viewer, PetscInt d, PetscInt *dims);

    259: #if defined(PETSC_HAVE_HDF4)
    262: PetscErrorCode VecView_MPI_HDF4_DA2d(Vec xin,PetscViewer viewer)
    263: {
    265:   PetscInt       dims[2];
    266:   DA             da;
    267:   Vec            natural;

    271:   PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);
    272:   if (!da) SETERRQ(PETSC_ERR_ARG_WRONG,"Vector not generated from a DA");

    274:   dims[0] = da->M;
    275:   dims[1] = da->N;

    277:   DACreateNaturalVector(da,&natural);
    278:   DAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);
    279:   DAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);
    280:   VecView_MPI_HDF4_Ex(natural, viewer, 2, dims);
    281:   VecDestroy(natural);
    282:   return(0);
    283: }
    284: #endif

    286: #if defined(PETSC_HAVE_PNETCDF)
    289: PetscErrorCode VecView_MPI_Netcdf_DA(Vec xin,PetscViewer viewer)
    290: {
    292:   PetscInt       ncid,xstart,xdim_num=1;
    293:   PetscInt       i,j,len,dim,m,n,p,dof,swidth,M,N,P;
    294:   PetscInt       xin_dim,xin_id,xin_n,xin_N,xyz_dim,xyz_id,xyz_n,xyz_N;
    295:   PetscInt       *lx,*ly,*lz;
    296:   PetscScalar    *xarray;
    297:   DA             da,dac;
    298:   Vec            natural,xyz;
    299:   DAStencilType  stencil;
    300:   DAPeriodicType periodic;
    301:   MPI_Comm       comm;

    304:   PetscObjectGetComm((PetscObject)xin,&comm);
    305:   PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);
    306:   if (!da) SETERRQ(PETSC_ERR_ARG_WRONG,,"Vector not generated from a DA");
    307:   DAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&dof,&swidth,&periodic,&stencil);

    309:   /* create the appropriate DA to map the coordinates to natural ordering */
    310:   DAGetOwnershipRange(da,&lx,&ly,&lz);
    311:   if (dim == 1) {
    312:     DACreate1d(comm,DA_NONPERIODIC,m,dim,0,lx,&dac);
    313:   } else if (dim == 2) {
    314:     DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,m,n,M,N,dim,0,lx,ly,&dac);
    315:   } else if (dim == 3) {
    316:     DACreate3d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,m,n,p,M,N,P,dim,0,lx,ly,lz,&dac);
    317:   } else {
    318:     SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Dimension is not 1 2 or 3: %D\n",dim);
    319:   }
    320:   DACreateNaturalVector(dac,&xyz);
    321:   PetscObjectSetOptionsPrefix((PetscObject)xyz,"coor_");
    322:   DAGlobalToNaturalBegin(dac,da->coordinates,INSERT_VALUES,xyz);
    323:   DAGlobalToNaturalEnd(dac,da->coordinates,INSERT_VALUES,xyz);
    324:   /* Create the DA vector in natural ordering */
    325:   DACreateNaturalVector(da,&natural);
    326:   DAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);
    327:   DAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);
    328:   /* Write the netCDF dataset */
    329:   PetscViewerNetcdfGetID(viewer,&ncid);
    330:   if (ncid < 0) SETERRQ(PETSC_ERR_ORDER,"First call PetscViewerNetcdfOpen to create NetCDF dataset");
    331:   /* define dimensions */
    332:   VecGetSize(xin,&xin_N);
    333:   VecGetLocalSize(xin,&xin_n);
    334:   ncmpi_def_dim(ncid,"PETSc_DA_Vector_Global_Size",xin_N,&xin_dim);
    335:   VecGetSize(xyz,&xyz_N);
    336:   VecGetLocalSize(xyz,&xyz_n);
    337:   ncmpi_def_dim(ncid,"PETSc_DA_Coordinate_Vector_Global_Size",xyz_N,&xyz_dim);
    338:   /* define variables */
    339:   ncmpi_def_var(ncid,"PETSc_DA_Vector",NC_DOUBLE,xdim_num,&xin_dim,&xin_id);
    340:   ncmpi_def_var(ncid,"PETSc_DA_Coordinate_Vector",NC_DOUBLE,xdim_num,&xyz_dim,&xyz_id);
    341:   /* leave define mode */
    342:   ncmpi_enddef(ncid);
    343:   /* store the vector */
    344:   VecGetArray(xin,&xarray);
    345:   VecGetOwnershipRange(xin,&xstart,PETSC_NULL);
    346:   ncmpi_put_vara_double_all(ncid,xin_id,(const size_t*)&xstart,(const size_t*)&xin_n,xarray);
    347:   VecRestoreArray(xin,&xarray);
    348:   /* store the coordinate vector */
    349:   VecGetArray(xyz,&xarray);
    350:   VecGetOwnershipRange(xyz,&xstart,PETSC_NULL);
    351:   ncmpi_put_vara_double_all(ncid,xyz_id,(const size_t*)&xstart,(const size_t*)&xyz_n,xarray);
    352:   VecRestoreArray(xyz,&xarray);
    353:   /* destroy the vectors and da */
    354:   VecDestroy(natural);
    355:   VecDestroy(xyz);
    356:   DADestroy(dac);
    357:   return(0);
    358: }
    359: #endif

    361: EXTERN PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer);

    366: PetscErrorCode PETSCDM_DLLEXPORT VecView_MPI_DA(Vec xin,PetscViewer viewer)
    367: {
    368:   DA             da;
    370:   PetscInt       dim;
    371:   Vec            natural;
    372:   PetscTruth     isdraw;
    373: #if defined(PETSC_HAVE_HDF4)
    374:   PetscTruth     ishdf4;
    375: #endif
    376: #if defined(PETSC_HAVE_PNETCDF)
    377:   PetscTruth     isnetcdf;
    378: #endif
    379:   const char     *prefix;

    382:   PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);
    383:   if (!da) SETERRQ(PETSC_ERR_ARG_WRONG,"Vector not generated from a DA");
    384:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
    385: #if defined(PETSC_HAVE_HDF4)
    386:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_HDF4,&ishdf4);
    387: #endif
    388: #if defined(PETSC_HAVE_PNETCDF)
    389:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_NETCDF,&isnetcdf);
    390: #endif
    391:   if (isdraw) {
    392:     DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);
    393:     if (dim == 1) {
    394:       VecView_MPI_Draw_DA1d(xin,viewer);
    395:     } else if (dim == 2) {
    396:       VecView_MPI_Draw_DA2d(xin,viewer);
    397:     } else {
    398:       SETERRQ1(PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DA %D",dim);
    399:     }
    400: #if defined(PETSC_HAVE_HDF4)
    401:   } else if (ishdf4) {
    402:     DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);
    403:     switch (dim) {
    404:     case 2:
    405:       VecView_MPI_HDF4_DA2d(xin,viewer);
    406:       break;
    407:     default:
    408:       SETERRQ1(PETSC_ERR_SUP,"Cannot view HDF4 vector associated with this dimensional DA %D",dim);
    409:     }
    410: #endif
    411: #if defined(PETSC_HAVE_PNETCDF)
    412:   } else if (isnetcdf) {
    413:     VecView_MPI_Netcdf_DA(xin,viewer);
    414: #endif
    415:   } else {
    416:     /* call viewer on natural ordering */
    417:     PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);
    418:     DACreateNaturalVector(da,&natural);
    419:     PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);
    420:     DAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);
    421:     DAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);
    422:     PetscObjectName((PetscObject)xin);
    423:     PetscObjectSetName((PetscObject)natural,xin->name);
    424:     VecView(natural,viewer);
    425:     VecDestroy(natural);
    426:   }
    427:   return(0);
    428: }

    434: PetscErrorCode PETSCDM_DLLEXPORT VecLoadIntoVector_Binary_DA(PetscViewer viewer,Vec xin)
    435: {
    436:   DA             da;
    438:   Vec            natural;
    439:   const char     *prefix;

    442:   PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);
    443:   if (!da) SETERRQ(PETSC_ERR_ARG_WRONG,"Vector not generated from a DA");
    444:   PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);
    445:   DACreateNaturalVector(da,&natural);
    446:   PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);
    447:   VecLoadIntoVector(viewer,natural);
    448:   DANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);
    449:   DANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);
    450:   VecDestroy(natural);
    451:   PetscLogInfo((xin,"VecLoadIntoVector_Binary_DA:Loading vector from natural ordering into DA\n"));
    452:   return(0);
    453: }