1 2 /* 3 Plots vectors obtained with DMDACreate2d() 4 */ 5 6 #include <private/daimpl.h> /*I "petscdmda.h" I*/ 7 #include <private/vecimpl.h> 8 9 /* 10 The data that is passed into the graphics callback 11 */ 12 typedef struct { 13 PetscInt m,n,step,k; 14 PetscReal min,max,scale; 15 PetscScalar *xy,*v; 16 PetscBool showgrid; 17 } ZoomCtx; 18 19 /* 20 This does the drawing for one particular field 21 in one particular set of coordinates. It is a callback 22 called from PetscDrawZoom() 23 */ 24 #undef __FUNCT__ 25 #define __FUNCT__ "VecView_MPI_Draw_DA2d_Zoom" 26 PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw,void *ctx) 27 { 28 ZoomCtx *zctx = (ZoomCtx*)ctx; 29 PetscErrorCode ierr; 30 PetscInt m,n,i,j,k,step,id,c1,c2,c3,c4; 31 PetscReal s,min,x1,x2,x3,x4,y_1,y2,y3,y4; 32 PetscScalar *v,*xy; 33 34 PetscFunctionBegin; 35 m = zctx->m; 36 n = zctx->n; 37 step = zctx->step; 38 k = zctx->k; 39 v = zctx->v; 40 xy = zctx->xy; 41 s = zctx->scale; 42 min = zctx->min; 43 44 /* PetscDraw the contour plot patch */ 45 for (j=0; j<n-1; j++) { 46 for (i=0; i<m-1; i++) { 47 #if !defined(PETSC_USE_COMPLEX) 48 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)); 49 id = i+j*m+1; x2 = xy[2*id];y2 = y_1; c2 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 50 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)); 51 id = i+j*m+m; x4 = x1; y4 = y3; c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 52 #else 53 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)); 54 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)); 55 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)); 56 id = i+j*m+m; x4 = x1; y4 = y3; c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min)); 57 #endif 58 ierr = PetscDrawTriangle(draw,x1,y_1,x2,y2,x3,y3,c1,c2,c3);CHKERRQ(ierr); 59 ierr = PetscDrawTriangle(draw,x1,y_1,x3,y3,x4,y4,c1,c3,c4);CHKERRQ(ierr); 60 if (zctx->showgrid) { 61 ierr = PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);CHKERRQ(ierr); 62 ierr = PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);CHKERRQ(ierr); 63 ierr = PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);CHKERRQ(ierr); 64 ierr = PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);CHKERRQ(ierr); 65 } 66 } 67 } 68 PetscFunctionReturn(0); 69 } 70 71 #undef __FUNCT__ 72 #define __FUNCT__ "VecView_MPI_Draw_DA2d" 73 PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer) 74 { 75 DM da,dac,dag; 76 PetscErrorCode ierr; 77 PetscMPIInt rank; 78 PetscInt N,s,M,w; 79 const PetscInt *lx,*ly; 80 PetscReal coors[4],ymin,ymax,xmin,xmax; 81 PetscDraw draw,popup; 82 PetscBool isnull,useports = PETSC_FALSE; 83 MPI_Comm comm; 84 Vec xlocal,xcoor,xcoorl; 85 DMDABoundaryType bx,by; 86 DMDAStencilType st; 87 ZoomCtx zctx; 88 PetscDrawViewPorts *ports = PETSC_NULL; 89 PetscViewerFormat format; 90 PetscInt *displayfields; 91 PetscInt ndisplayfields,i,nbounds; 92 PetscBool flg; 93 const PetscReal *bounds; 94 95 PetscFunctionBegin; 96 zctx.showgrid = PETSC_FALSE; 97 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 98 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 99 ierr = PetscViewerDrawGetBounds(viewer,&nbounds,&bounds);CHKERRQ(ierr); 100 101 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 102 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 103 104 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 105 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 106 107 ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&st);CHKERRQ(ierr); 108 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,PETSC_NULL);CHKERRQ(ierr); 109 110 /* 111 Obtain a sequential vector that is going to contain the local values plus ONE layer of 112 ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will 113 update the local values pluse ONE layer of ghost values. 114 */ 115 ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr); 116 if (!xlocal) { 117 if (bx != DMDA_BOUNDARY_NONE || by != DMDA_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 118 /* 119 if original da is not of stencil width one, or periodic or not a box stencil then 120 create a special DMDA to handle one level of ghost points for graphics 121 */ 122 ierr = DMDACreate2d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,w,1,lx,ly,&dac);CHKERRQ(ierr); 123 ierr = PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");CHKERRQ(ierr); 124 } else { 125 /* otherwise we can use the da we already have */ 126 dac = da; 127 } 128 /* create local vector for holding ghosted values used in graphics */ 129 ierr = DMCreateLocalVector(dac,&xlocal);CHKERRQ(ierr); 130 if (dac != da) { 131 /* don't keep any public reference of this DMDA, is is only available through xlocal */ 132 ierr = PetscObjectDereference((PetscObject)dac);CHKERRQ(ierr); 133 } else { 134 /* remove association between xlocal and da, because below we compose in the opposite 135 direction and if we left this connect we'd get a loop, so the objects could 136 never be destroyed */ 137 ierr = PetscObjectRemoveReference((PetscObject)xlocal,"DM");CHKERRQ(ierr); 138 } 139 ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr); 140 ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr); 141 } else { 142 if (bx != DMDA_BOUNDARY_NONE || by != DMDA_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) { 143 ierr = PetscObjectQuery((PetscObject)xlocal,"DM",(PetscObject*)&dac);CHKERRQ(ierr); 144 } else { 145 dac = da; 146 } 147 } 148 149 /* 150 Get local (ghosted) values of vector 151 */ 152 ierr = DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 153 ierr = DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 154 ierr = VecGetArray(xlocal,&zctx.v);CHKERRQ(ierr); 155 156 /* get coordinates of nodes */ 157 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 158 if (!xcoor) { 159 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);CHKERRQ(ierr); 160 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 161 } 162 163 /* 164 Determine the min and max coordinates in plot 165 */ 166 ierr = VecStrideMin(xcoor,0,PETSC_NULL,&xmin);CHKERRQ(ierr); 167 ierr = VecStrideMax(xcoor,0,PETSC_NULL,&xmax);CHKERRQ(ierr); 168 ierr = VecStrideMin(xcoor,1,PETSC_NULL,&ymin);CHKERRQ(ierr); 169 ierr = VecStrideMax(xcoor,1,PETSC_NULL,&ymax);CHKERRQ(ierr); 170 coors[0] = xmin - .05*(xmax- xmin); coors[2] = xmax + .05*(xmax - xmin); 171 coors[1] = ymin - .05*(ymax- ymin); coors[3] = ymax + .05*(ymax - ymin); 172 ierr = PetscInfo4(da,"Preparing DMDA 2d contour plot coordinates %G %G %G %G\n",coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 173 174 /* 175 get local ghosted version of coordinates 176 */ 177 ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr); 178 if (!xcoorl) { 179 /* create DMDA to get local version of graphics */ 180 ierr = DMDACreate2d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,2,1,lx,ly,&dag);CHKERRQ(ierr); 181 ierr = PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");CHKERRQ(ierr); 182 ierr = DMCreateLocalVector(dag,&xcoorl);CHKERRQ(ierr); 183 ierr = PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);CHKERRQ(ierr); 184 ierr = PetscObjectDereference((PetscObject)dag);CHKERRQ(ierr); 185 ierr = PetscObjectDereference((PetscObject)xcoorl);CHKERRQ(ierr); 186 } else { 187 ierr = PetscObjectQuery((PetscObject)xcoorl,"DM",(PetscObject*)&dag);CHKERRQ(ierr); 188 } 189 ierr = DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 190 ierr = DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 191 ierr = VecGetArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 192 193 /* 194 Get information about size of area each processor must do graphics for 195 */ 196 ierr = DMDAGetInfo(dac,0,&M,&N,0,0,0,0,&zctx.step,0,&bx,&by,0,0);CHKERRQ(ierr); 197 ierr = DMDAGetGhostCorners(dac,0,0,0,&zctx.m,&zctx.n,0);CHKERRQ(ierr); 198 199 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_contour_grid",&zctx.showgrid,PETSC_NULL);CHKERRQ(ierr); 200 ierr = PetscMalloc(zctx.step*sizeof(PetscInt),&displayfields);CHKERRQ(ierr); 201 for (i=0; i<zctx.step; i++) displayfields[i] = i; 202 ndisplayfields = zctx.step; 203 ierr = PetscOptionsGetIntArray(PETSC_NULL,"-draw_fields",displayfields,&ndisplayfields,&flg);CHKERRQ(ierr); 204 if (!flg) ndisplayfields = zctx.step; 205 206 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 207 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_ports",&useports,PETSC_NULL);CHKERRQ(ierr); 208 if (useports || format == PETSC_VIEWER_DRAW_PORTS){ 209 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 210 ierr = PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);CHKERRQ(ierr); 211 } 212 213 /* 214 Loop over each field; drawing each in a different window 215 */ 216 for (i=0; i<ndisplayfields; i++) { 217 zctx.k = displayfields[i]; 218 if (useports) { 219 ierr = PetscDrawViewPortsSet(ports,i);CHKERRQ(ierr); 220 } else { 221 ierr = PetscViewerDrawGetDraw(viewer,i,&draw);CHKERRQ(ierr); 222 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 223 } 224 225 /* 226 Determine the min and max color in plot 227 */ 228 ierr = VecStrideMin(xin,zctx.k,PETSC_NULL,&zctx.min);CHKERRQ(ierr); 229 ierr = VecStrideMax(xin,zctx.k,PETSC_NULL,&zctx.max);CHKERRQ(ierr); 230 if (zctx.k < nbounds) { 231 zctx.min = PetscMin(zctx.min,bounds[2*zctx.k]); 232 zctx.max = PetscMax(zctx.max,bounds[2*zctx.k+1]); 233 } 234 if (zctx.min == zctx.max) { 235 zctx.min -= 1.e-12; 236 zctx.max += 1.e-12; 237 } 238 239 if (!rank) { 240 const char *title; 241 242 ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr); 243 if (title) { 244 ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr); 245 } 246 } 247 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 248 ierr = PetscInfo2(da,"DMDA 2d contour plot min %G max %G\n",zctx.min,zctx.max);CHKERRQ(ierr); 249 250 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 251 if (popup) {ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);} 252 253 zctx.scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/(zctx.max - zctx.min); 254 255 ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr); 256 } 257 ierr = PetscFree(displayfields);CHKERRQ(ierr); 258 ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr); 259 260 ierr = VecRestoreArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 261 ierr = VecRestoreArray(xlocal,&zctx.v);CHKERRQ(ierr); 262 PetscFunctionReturn(0); 263 } 264 265 266 #if defined(PETSC_HAVE_HDF5) 267 #undef __FUNCT__ 268 #define __FUNCT__ "VecView_MPI_HDF5_DA" 269 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 270 { 271 PetscErrorCode ierr; 272 DM dm; 273 DM_DA *da; 274 hsize_t dim,dims[5]; 275 hsize_t count[5]; 276 hsize_t offset[5]; 277 PetscInt cnt = 0; 278 PetscScalar *x; 279 const char *vecname; 280 hid_t filespace; /* file dataspace identifier */ 281 hid_t plist_id; /* property list identifier */ 282 hid_t dset_id; /* dataset identifier */ 283 hid_t memspace; /* memory dataspace identifier */ 284 hid_t file_id; 285 herr_t status; 286 287 PetscFunctionBegin; 288 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 289 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&dm);CHKERRQ(ierr); 290 if (!dm) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 291 da = (DM_DA*)dm->data; 292 293 /* Create the dataspace for the dataset */ 294 dim = PetscHDF5IntCast(da->dim + ((da->w == 1) ? 0 : 1)); 295 if (da->dim == 3) dims[cnt++] = PetscHDF5IntCast(da->P); 296 if (da->dim > 1) dims[cnt++] = PetscHDF5IntCast(da->N); 297 dims[cnt++] = PetscHDF5IntCast(da->M); 298 if (da->w > 1) dims[cnt++] = PetscHDF5IntCast(da->w); 299 #if defined(PETSC_USE_COMPLEX) 300 dim++; 301 dims[cnt++] = 2; 302 #endif 303 filespace = H5Screate_simple(dim, dims, NULL); 304 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 305 306 /* Create the dataset with default properties and close filespace */ 307 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 308 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 309 dset_id = H5Dcreate2(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); 310 #else 311 dset_id = H5Dcreate(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT); 312 #endif 313 if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()"); 314 status = H5Sclose(filespace);CHKERRQ(status); 315 316 /* Each process defines a dataset and writes it to the hyperslab in the file */ 317 cnt = 0; 318 if (da->dim == 3) offset[cnt++] = PetscHDF5IntCast(da->zs); 319 if (da->dim > 1) offset[cnt++] = PetscHDF5IntCast(da->ys); 320 offset[cnt++] = PetscHDF5IntCast(da->xs/da->w); 321 if (da->w > 1) offset[cnt++] = 0; 322 #if defined(PETSC_USE_COMPLEX) 323 offset[cnt++] = 0; 324 #endif 325 cnt = 0; 326 if (da->dim == 3) count[cnt++] = PetscHDF5IntCast(da->ze - da->zs); 327 if (da->dim > 1) count[cnt++] = PetscHDF5IntCast(da->ye - da->ys); 328 count[cnt++] = PetscHDF5IntCast((da->xe - da->xs)/da->w); 329 if (da->w > 1) count[cnt++] = PetscHDF5IntCast(da->w); 330 #if defined(PETSC_USE_COMPLEX) 331 count[cnt++] = 2; 332 #endif 333 memspace = H5Screate_simple(dim, count, NULL); 334 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 335 336 337 filespace = H5Dget_space(dset_id); 338 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()"); 339 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 340 341 /* Create property list for collective dataset write */ 342 plist_id = H5Pcreate(H5P_DATASET_XFER); 343 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 344 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 345 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 346 #endif 347 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 348 349 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 350 status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 351 status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status); 352 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 353 354 /* Close/release resources */ 355 status = H5Pclose(plist_id);CHKERRQ(status); 356 status = H5Sclose(filespace);CHKERRQ(status); 357 status = H5Sclose(memspace);CHKERRQ(status); 358 status = H5Dclose(dset_id);CHKERRQ(status); 359 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 360 PetscFunctionReturn(0); 361 } 362 #endif 363 364 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 365 366 #if defined(PETSC_HAVE_MPIIO) 367 #undef __FUNCT__ 368 #define __FUNCT__ "DMDAArrayMPIIO" 369 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 370 { 371 PetscErrorCode ierr; 372 MPI_File mfdes; 373 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 374 MPI_Datatype view; 375 const PetscScalar *array; 376 MPI_Offset off; 377 MPI_Aint ub,ul; 378 PetscInt type,rows,vecrows,tr[2]; 379 DM_DA *dd = (DM_DA*)da->data; 380 381 PetscFunctionBegin; 382 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 383 if (!write) { 384 /* Read vector header. */ 385 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 386 type = tr[0]; 387 rows = tr[1]; 388 if (type != VEC_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not vector next in file"); 389 if (rows != vecrows) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 390 } else { 391 tr[0] = VEC_FILE_CLASSID; 392 tr[1] = vecrows; 393 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 394 } 395 396 dof = PetscMPIIntCast(dd->w); 397 gsizes[0] = dof; gsizes[1] = PetscMPIIntCast(dd->M); gsizes[2] = PetscMPIIntCast(dd->N); gsizes[3] = PetscMPIIntCast(dd->P); 398 lsizes[0] = dof;lsizes[1] = PetscMPIIntCast((dd->xe-dd->xs)/dof); lsizes[2] = PetscMPIIntCast(dd->ye-dd->ys); lsizes[3] = PetscMPIIntCast(dd->ze-dd->zs); 399 lstarts[0] = 0; lstarts[1] = PetscMPIIntCast(dd->xs/dof); lstarts[2] = PetscMPIIntCast(dd->ys); lstarts[3] = PetscMPIIntCast(dd->zs); 400 ierr = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 401 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 402 403 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 404 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 405 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char *)"native",MPI_INFO_NULL);CHKERRQ(ierr); 406 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 407 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 408 if (write) { 409 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 410 } else { 411 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 412 } 413 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 414 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 415 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 416 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 417 PetscFunctionReturn(0); 418 } 419 #endif 420 421 EXTERN_C_BEGIN 422 #undef __FUNCT__ 423 #define __FUNCT__ "VecView_MPI_DA" 424 PetscErrorCode VecView_MPI_DA(Vec xin,PetscViewer viewer) 425 { 426 DM da; 427 PetscErrorCode ierr; 428 PetscInt dim; 429 Vec natural; 430 PetscBool isdraw,isvtk; 431 #if defined(PETSC_HAVE_HDF5) 432 PetscBool ishdf5; 433 #endif 434 const char *prefix,*name; 435 436 PetscFunctionBegin; 437 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 438 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 439 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 440 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 441 #if defined(PETSC_HAVE_HDF5) 442 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 443 #endif 444 if (isdraw) { 445 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 446 if (dim == 1) { 447 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 448 } else if (dim == 2) { 449 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 450 } else { 451 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 452 } 453 } else if (isvtk) { /* Duplicate the Vec and hold a reference to the DM */ 454 Vec Y; 455 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 456 ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr); 457 ierr = VecCopy(xin,Y);CHKERRQ(ierr); 458 ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,(PetscObject)Y);CHKERRQ(ierr); 459 #if defined(PETSC_HAVE_HDF5) 460 } else if (ishdf5) { 461 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 462 #endif 463 } else { 464 #if defined(PETSC_HAVE_MPIIO) 465 PetscBool isbinary,isMPIIO; 466 467 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 468 if (isbinary) { 469 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 470 if (isMPIIO) { 471 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 472 PetscFunctionReturn(0); 473 } 474 } 475 #endif 476 477 /* call viewer on natural ordering */ 478 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 479 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 480 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 481 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 482 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 483 ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr); 484 ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr); 485 ierr = VecView(natural,viewer);CHKERRQ(ierr); 486 ierr = VecDestroy(&natural);CHKERRQ(ierr); 487 } 488 PetscFunctionReturn(0); 489 } 490 EXTERN_C_END 491 492 #if defined(PETSC_HAVE_HDF5) 493 #undef __FUNCT__ 494 #define __FUNCT__ "VecLoad_HDF5_DA" 495 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 496 { 497 DM da; 498 PetscErrorCode ierr; 499 hsize_t dim; 500 hsize_t count[5]; 501 hsize_t offset[5]; 502 PetscInt cnt = 0; 503 PetscScalar *x; 504 const char *vecname; 505 hid_t filespace; /* file dataspace identifier */ 506 hid_t plist_id; /* property list identifier */ 507 hid_t dset_id; /* dataset identifier */ 508 hid_t memspace; /* memory dataspace identifier */ 509 hid_t file_id; 510 herr_t status; 511 DM_DA *dd; 512 513 PetscFunctionBegin; 514 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 515 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 516 dd = (DM_DA*)da->data; 517 518 /* Create the dataspace for the dataset */ 519 dim = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1)); 520 #if defined(PETSC_USE_COMPLEX) 521 dim++; 522 #endif 523 524 /* Create the dataset with default properties and close filespace */ 525 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 526 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 527 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 528 #else 529 dset_id = H5Dopen(file_id, vecname); 530 #endif 531 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 532 filespace = H5Dget_space(dset_id); 533 534 /* Each process defines a dataset and reads it from the hyperslab in the file */ 535 cnt = 0; 536 if (dd->dim == 3) offset[cnt++] = PetscHDF5IntCast(dd->zs); 537 if (dd->dim > 1) offset[cnt++] = PetscHDF5IntCast(dd->ys); 538 offset[cnt++] = PetscHDF5IntCast(dd->xs/dd->w); 539 if (dd->w > 1) offset[cnt++] = 0; 540 #if defined(PETSC_USE_COMPLEX) 541 offset[cnt++] = 0; 542 #endif 543 cnt = 0; 544 if (dd->dim == 3) count[cnt++] = PetscHDF5IntCast(dd->ze - dd->zs); 545 if (dd->dim > 1) count[cnt++] = PetscHDF5IntCast(dd->ye - dd->ys); 546 count[cnt++] = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w); 547 if (dd->w > 1) count[cnt++] = PetscHDF5IntCast(dd->w); 548 #if defined(PETSC_USE_COMPLEX) 549 count[cnt++] = 2; 550 #endif 551 memspace = H5Screate_simple(dim, count, NULL); 552 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 553 554 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 555 556 /* Create property list for collective dataset write */ 557 plist_id = H5Pcreate(H5P_DATASET_XFER); 558 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 559 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 560 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 561 #endif 562 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 563 564 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 565 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 566 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 567 568 /* Close/release resources */ 569 status = H5Pclose(plist_id);CHKERRQ(status); 570 status = H5Sclose(filespace);CHKERRQ(status); 571 status = H5Sclose(memspace);CHKERRQ(status); 572 status = H5Dclose(dset_id);CHKERRQ(status); 573 PetscFunctionReturn(0); 574 } 575 #endif 576 577 #undef __FUNCT__ 578 #define __FUNCT__ "VecLoad_Binary_DA" 579 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 580 { 581 DM da; 582 PetscErrorCode ierr; 583 Vec natural; 584 const char *prefix; 585 PetscInt bs; 586 PetscBool flag; 587 DM_DA *dd; 588 #if defined(PETSC_HAVE_MPIIO) 589 PetscBool isMPIIO; 590 #endif 591 592 PetscFunctionBegin; 593 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 594 dd = (DM_DA*)da->data; 595 #if defined(PETSC_HAVE_MPIIO) 596 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 597 if (isMPIIO) { 598 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 599 PetscFunctionReturn(0); 600 } 601 #endif 602 603 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 604 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 605 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 606 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 607 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 608 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 609 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 610 ierr = VecDestroy(&natural);CHKERRQ(ierr); 611 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 612 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 613 if (flag && bs != dd->w) { 614 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 615 } 616 PetscFunctionReturn(0); 617 } 618 619 EXTERN_C_BEGIN 620 #undef __FUNCT__ 621 #define __FUNCT__ "VecLoad_Default_DA" 622 PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer) 623 { 624 PetscErrorCode ierr; 625 DM da; 626 PetscBool isbinary; 627 #if defined(PETSC_HAVE_HDF5) 628 PetscBool ishdf5; 629 #endif 630 631 PetscFunctionBegin; 632 ierr = PetscObjectQuery((PetscObject)xin,"DM",(PetscObject*)&da);CHKERRQ(ierr); 633 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 634 635 #if defined(PETSC_HAVE_HDF5) 636 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 637 #endif 638 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 639 640 if (isbinary) { 641 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 642 #if defined(PETSC_HAVE_HDF5) 643 } else if (ishdf5) { 644 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 645 #endif 646 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 647 PetscFunctionReturn(0); 648 } 649 EXTERN_C_END 650