Actual source code: gr2.c
1: /*
2: Plots vectors obtained with DMDACreate2d()
3: */
5: #include <petsc/private/dmdaimpl.h>
6: #include <petsc/private/glvisvecimpl.h>
7: #include <petsc/private/viewerhdf5impl.h>
8: #include <petscdraw.h>
10: /*
11: The data that is passed into the graphics callback
12: */
13: typedef struct {
14: PetscMPIInt rank;
15: PetscInt m, n, dof, k;
16: PetscReal xmin, xmax, ymin, ymax, min, max;
17: const PetscScalar *xy, *v;
18: PetscBool showaxis, showgrid;
19: const char *name0, *name1;
20: } ZoomCtx;
22: /*
23: This does the drawing for one particular field
24: in one particular set of coordinates. It is a callback
25: called from PetscDrawZoom()
26: */
27: static PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw, void *ctx)
28: {
29: ZoomCtx *zctx = (ZoomCtx *)ctx;
30: PetscInt m, n, i, j, k, dof, id;
31: int c1, c2, c3, c4;
32: PetscReal min, max, x1, x2, x3, x4, y_1, y2, y3, y4;
33: const PetscScalar *xy, *v;
35: PetscFunctionBegin;
36: m = zctx->m;
37: n = zctx->n;
38: dof = zctx->dof;
39: k = zctx->k;
40: xy = zctx->xy;
41: v = zctx->v;
42: min = zctx->min;
43: max = zctx->max;
45: /* PetscDraw the contour plot patch */
46: PetscDrawCollectiveBegin(draw);
47: for (j = 0; j < n - 1; j++) {
48: for (i = 0; i < m - 1; i++) {
49: id = i + j * m;
50: x1 = PetscRealPart(xy[2 * id]);
51: y_1 = PetscRealPart(xy[2 * id + 1]);
52: c1 = PetscDrawRealToColor(PetscRealPart(v[k + dof * id]), min, max);
54: id = i + j * m + 1;
55: x2 = PetscRealPart(xy[2 * id]);
56: y2 = PetscRealPart(xy[2 * id + 1]);
57: c2 = PetscDrawRealToColor(PetscRealPart(v[k + dof * id]), min, max);
59: id = i + j * m + 1 + m;
60: x3 = PetscRealPart(xy[2 * id]);
61: y3 = PetscRealPart(xy[2 * id + 1]);
62: c3 = PetscDrawRealToColor(PetscRealPart(v[k + dof * id]), min, max);
64: id = i + j * m + m;
65: x4 = PetscRealPart(xy[2 * id]);
66: y4 = PetscRealPart(xy[2 * id + 1]);
67: c4 = PetscDrawRealToColor(PetscRealPart(v[k + dof * id]), min, max);
69: PetscCall(PetscDrawTriangle(draw, x1, y_1, x2, y2, x3, y3, c1, c2, c3));
70: PetscCall(PetscDrawTriangle(draw, x1, y_1, x3, y3, x4, y4, c1, c3, c4));
71: if (zctx->showgrid) {
72: PetscCall(PetscDrawLine(draw, x1, y_1, x2, y2, PETSC_DRAW_BLACK));
73: PetscCall(PetscDrawLine(draw, x2, y2, x3, y3, PETSC_DRAW_BLACK));
74: PetscCall(PetscDrawLine(draw, x3, y3, x4, y4, PETSC_DRAW_BLACK));
75: PetscCall(PetscDrawLine(draw, x4, y4, x1, y_1, PETSC_DRAW_BLACK));
76: }
77: }
78: }
79: if (zctx->showaxis && !zctx->rank) {
80: if (zctx->name0 || zctx->name1) {
81: PetscReal xl, yl, xr, yr, x, y;
82: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
83: x = xl + .30 * (xr - xl);
84: xl = xl + .01 * (xr - xl);
85: y = yr - .30 * (yr - yl);
86: yl = yl + .01 * (yr - yl);
87: if (zctx->name0) PetscCall(PetscDrawString(draw, x, yl, PETSC_DRAW_BLACK, zctx->name0));
88: if (zctx->name1) PetscCall(PetscDrawStringVertical(draw, xl, y, PETSC_DRAW_BLACK, zctx->name1));
89: }
90: /*
91: Ideally we would use the PetscDrawAxis object to manage displaying the coordinate limits
92: but that may require some refactoring.
93: */
94: {
95: double xmin = (double)zctx->xmin, ymin = (double)zctx->ymin;
96: double xmax = (double)zctx->xmax, ymax = (double)zctx->ymax;
97: char value[16];
98: size_t len;
99: PetscReal w;
100: PetscCall(PetscSNPrintf(value, 16, "%0.2e", xmin));
101: PetscCall(PetscDrawString(draw, xmin, ymin - .05 * (ymax - ymin), PETSC_DRAW_BLACK, value));
102: PetscCall(PetscSNPrintf(value, 16, "%0.2e", xmax));
103: PetscCall(PetscStrlen(value, &len));
104: PetscCall(PetscDrawStringGetSize(draw, &w, NULL));
105: PetscCall(PetscDrawString(draw, xmax - len * w, ymin - .05 * (ymax - ymin), PETSC_DRAW_BLACK, value));
106: PetscCall(PetscSNPrintf(value, 16, "%0.2e", ymin));
107: PetscCall(PetscDrawString(draw, xmin - .05 * (xmax - xmin), ymin, PETSC_DRAW_BLACK, value));
108: PetscCall(PetscSNPrintf(value, 16, "%0.2e", ymax));
109: PetscCall(PetscDrawString(draw, xmin - .05 * (xmax - xmin), ymax, PETSC_DRAW_BLACK, value));
110: }
111: }
112: PetscDrawCollectiveEnd(draw);
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin, PetscViewer viewer)
117: {
118: DM da, dac, dag;
119: PetscInt N, s, M, w, ncoors = 4;
120: const PetscInt *lx, *ly;
121: PetscReal coors[4];
122: PetscDraw draw, popup;
123: PetscBool isnull, useports = PETSC_FALSE;
124: MPI_Comm comm;
125: Vec xlocal, xcoor, xcoorl;
126: DMBoundaryType bx, by;
127: DMDAStencilType st;
128: ZoomCtx zctx;
129: PetscDrawViewPorts *ports = NULL;
130: PetscViewerFormat format;
131: PetscInt *displayfields;
132: PetscInt ndisplayfields, i, nbounds;
133: const PetscReal *bounds;
135: PetscFunctionBegin;
136: zctx.showgrid = PETSC_FALSE;
137: zctx.showaxis = PETSC_TRUE;
139: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
140: PetscCall(PetscDrawIsNull(draw, &isnull));
141: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
143: PetscCall(PetscViewerDrawGetBounds(viewer, &nbounds, &bounds));
145: PetscCall(VecGetDM(xin, &da));
146: PetscCheck(da, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA");
148: PetscCall(PetscObjectGetComm((PetscObject)xin, &comm));
149: PetscCallMPI(MPI_Comm_rank(comm, &zctx.rank));
151: PetscCall(DMDAGetInfo(da, NULL, &M, &N, NULL, &zctx.m, &zctx.n, NULL, &w, &s, &bx, &by, NULL, &st));
152: PetscCall(DMDAGetOwnershipRanges(da, &lx, &ly, NULL));
154: /*
155: Obtain a sequential vector that is going to contain the local values plus ONE layer of
156: ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will
157: update the local values plus ONE layer of ghost values.
158: */
159: PetscCall(PetscObjectQuery((PetscObject)da, "GraphicsGhosted", (PetscObject *)&xlocal));
160: if (!xlocal) {
161: if (bx != DM_BOUNDARY_NONE || by != DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) {
162: /*
163: if original da is not of stencil width one, or periodic or not a box stencil then
164: create a special DMDA to handle one level of ghost points for graphics
165: */
166: PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, M, N, zctx.m, zctx.n, w, 1, lx, ly, &dac));
167: PetscCall(DMSetUp(dac));
168: PetscCall(PetscInfo(da, "Creating auxiliary DMDA for managing graphics ghost points\n"));
169: } else {
170: /* otherwise we can use the da we already have */
171: dac = da;
172: }
173: /* create local vector for holding ghosted values used in graphics */
174: PetscCall(DMCreateLocalVector(dac, &xlocal));
175: if (dac != da) {
176: /* don't keep any public reference of this DMDA, it is only available through xlocal */
177: PetscCall(PetscObjectDereference((PetscObject)dac));
178: } else {
179: /* remove association between xlocal and da, because below we compose in the opposite
180: direction and if we left this connect we'd get a loop, so the objects could
181: never be destroyed */
182: PetscCall(PetscObjectRemoveReference((PetscObject)xlocal, "__PETSc_dm"));
183: }
184: PetscCall(PetscObjectCompose((PetscObject)da, "GraphicsGhosted", (PetscObject)xlocal));
185: PetscCall(PetscObjectDereference((PetscObject)xlocal));
186: } else {
187: if (bx != DM_BOUNDARY_NONE || by != DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) {
188: PetscCall(VecGetDM(xlocal, &dac));
189: } else {
190: dac = da;
191: }
192: }
194: /*
195: Get local (ghosted) values of vector
196: */
197: PetscCall(DMGlobalToLocalBegin(dac, xin, INSERT_VALUES, xlocal));
198: PetscCall(DMGlobalToLocalEnd(dac, xin, INSERT_VALUES, xlocal));
199: PetscCall(VecGetArrayRead(xlocal, &zctx.v));
201: /*
202: Get coordinates of nodes
203: */
204: PetscCall(DMGetCoordinates(da, &xcoor));
205: if (!xcoor) {
206: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
207: PetscCall(DMGetCoordinates(da, &xcoor));
208: }
210: /*
211: Determine the min and max coordinates in plot
212: */
213: PetscCall(VecStrideMin(xcoor, 0, NULL, &zctx.xmin));
214: PetscCall(VecStrideMax(xcoor, 0, NULL, &zctx.xmax));
215: PetscCall(VecStrideMin(xcoor, 1, NULL, &zctx.ymin));
216: PetscCall(VecStrideMax(xcoor, 1, NULL, &zctx.ymax));
217: PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_contour_axis", &zctx.showaxis, NULL));
218: if (zctx.showaxis) {
219: coors[0] = zctx.xmin - .05 * (zctx.xmax - zctx.xmin);
220: coors[1] = zctx.ymin - .05 * (zctx.ymax - zctx.ymin);
221: coors[2] = zctx.xmax + .05 * (zctx.xmax - zctx.xmin);
222: coors[3] = zctx.ymax + .05 * (zctx.ymax - zctx.ymin);
223: } else {
224: coors[0] = zctx.xmin;
225: coors[1] = zctx.ymin;
226: coors[2] = zctx.xmax;
227: coors[3] = zctx.ymax;
228: }
229: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-draw_coordinates", coors, &ncoors, NULL));
230: PetscCall(PetscInfo(da, "Preparing DMDA 2d contour plot coordinates %g %g %g %g\n", (double)coors[0], (double)coors[1], (double)coors[2], (double)coors[3]));
232: /*
233: Get local ghosted version of coordinates
234: */
235: PetscCall(PetscObjectQuery((PetscObject)da, "GraphicsCoordinateGhosted", (PetscObject *)&xcoorl));
236: if (!xcoorl) {
237: /* create DMDA to get local version of graphics */
238: PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, M, N, zctx.m, zctx.n, 2, 1, lx, ly, &dag));
239: PetscCall(DMSetUp(dag));
240: PetscCall(PetscInfo(dag, "Creating auxiliary DMDA for managing graphics coordinates ghost points\n"));
241: PetscCall(DMCreateLocalVector(dag, &xcoorl));
242: PetscCall(PetscObjectCompose((PetscObject)da, "GraphicsCoordinateGhosted", (PetscObject)xcoorl));
243: PetscCall(PetscObjectDereference((PetscObject)dag));
244: PetscCall(PetscObjectDereference((PetscObject)xcoorl));
245: } else PetscCall(VecGetDM(xcoorl, &dag));
246: PetscCall(DMGlobalToLocalBegin(dag, xcoor, INSERT_VALUES, xcoorl));
247: PetscCall(DMGlobalToLocalEnd(dag, xcoor, INSERT_VALUES, xcoorl));
248: PetscCall(VecGetArrayRead(xcoorl, &zctx.xy));
249: PetscCall(DMDAGetCoordinateName(da, 0, &zctx.name0));
250: PetscCall(DMDAGetCoordinateName(da, 1, &zctx.name1));
252: /*
253: Get information about size of area each processor must do graphics for
254: */
255: PetscCall(DMDAGetInfo(dac, NULL, &M, &N, NULL, NULL, NULL, NULL, &zctx.dof, NULL, &bx, &by, NULL, NULL));
256: PetscCall(DMDAGetGhostCorners(dac, NULL, NULL, NULL, &zctx.m, &zctx.n, NULL));
257: PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_contour_grid", &zctx.showgrid, NULL));
259: PetscCall(DMDASelectFields(da, &ndisplayfields, &displayfields));
260: PetscCall(PetscViewerGetFormat(viewer, &format));
261: PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_ports", &useports, NULL));
262: if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
263: if (useports) {
264: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
265: PetscCall(PetscDrawCheckResizedWindow(draw));
266: PetscCall(PetscDrawClear(draw));
267: PetscCall(PetscDrawViewPortsCreate(draw, ndisplayfields, &ports));
268: }
270: /*
271: Loop over each field; drawing each in a different window
272: */
273: for (i = 0; i < ndisplayfields; i++) {
274: zctx.k = displayfields[i];
276: /* determine the min and max value in plot */
277: PetscCall(VecStrideMin(xin, zctx.k, NULL, &zctx.min));
278: PetscCall(VecStrideMax(xin, zctx.k, NULL, &zctx.max));
279: if (zctx.k < nbounds) {
280: zctx.min = bounds[2 * zctx.k];
281: zctx.max = bounds[2 * zctx.k + 1];
282: }
283: if (zctx.min == zctx.max) {
284: zctx.min -= 1.e-12;
285: zctx.max += 1.e-12;
286: }
287: PetscCall(PetscInfo(da, "DMDA 2d contour plot min %g max %g\n", (double)zctx.min, (double)zctx.max));
289: if (useports) {
290: PetscCall(PetscDrawViewPortsSet(ports, i));
291: } else {
292: const char *title;
293: PetscCall(PetscViewerDrawGetDraw(viewer, i, &draw));
294: PetscCall(DMDAGetFieldName(da, zctx.k, &title));
295: if (title) PetscCall(PetscDrawSetTitle(draw, title));
296: }
298: PetscCall(PetscDrawGetPopup(draw, &popup));
299: PetscCall(PetscDrawScalePopup(popup, zctx.min, zctx.max));
300: PetscCall(PetscDrawSetCoordinates(draw, coors[0], coors[1], coors[2], coors[3]));
301: PetscCall(PetscDrawZoom(draw, VecView_MPI_Draw_DA2d_Zoom, &zctx));
302: if (!useports) PetscCall(PetscDrawSave(draw));
303: }
304: if (useports) {
305: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
306: PetscCall(PetscDrawSave(draw));
307: }
309: PetscCall(PetscDrawViewPortsDestroy(ports));
310: PetscCall(PetscFree(displayfields));
311: PetscCall(VecRestoreArrayRead(xcoorl, &zctx.xy));
312: PetscCall(VecRestoreArrayRead(xlocal, &zctx.v));
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: #if defined(PETSC_HAVE_HDF5)
317: static PetscErrorCode VecGetHDF5ChunkSize(DM_DA *da, Vec xin, PetscInt dimension, PetscInt timestep, hsize_t *chunkDims)
318: {
319: PetscMPIInt comm_size;
320: hsize_t chunk_size, target_size, dim;
321: hsize_t vec_size = sizeof(PetscScalar) * da->M * da->N * da->P * da->w;
322: hsize_t avg_local_vec_size, KiB = 1024, MiB = KiB * KiB, GiB = MiB * KiB, min_size = MiB;
323: hsize_t max_chunks = 64 * KiB; /* HDF5 internal limitation */
324: hsize_t max_chunk_size = 4 * GiB; /* HDF5 internal limitation */
325: PetscInt zslices = da->p, yslices = da->n, xslices = da->m;
327: PetscFunctionBegin;
328: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &comm_size));
329: avg_local_vec_size = (hsize_t)PetscCeilInt(vec_size, comm_size); /* we will attempt to use this as the chunk size */
331: target_size = (hsize_t)PetscMin((PetscInt64)vec_size, PetscMin((PetscInt64)max_chunk_size, PetscMax((PetscInt64)avg_local_vec_size, PetscMax(PetscCeilInt64(vec_size, max_chunks), (PetscInt64)min_size))));
332: /* following line uses sizeof(PetscReal) instead of sizeof(PetscScalar) because the last dimension of chunkDims[] captures the 2* when complex numbers are being used */
333: chunk_size = (hsize_t)PetscMax(1, chunkDims[0]) * PetscMax(1, chunkDims[1]) * PetscMax(1, chunkDims[2]) * PetscMax(1, chunkDims[3]) * PetscMax(1, chunkDims[4]) * PetscMax(1, chunkDims[5]) * sizeof(PetscReal);
335: /*
336: if size/rank > max_chunk_size, we need radical measures: even going down to
337: avg_local_vec_size is not enough, so we simply use chunk size of 4 GiB no matter
338: what, composed in the most efficient way possible.
339: N.B. this minimises the number of chunks, which may or may not be the optimal
340: solution. In a BG, for example, the optimal solution is probably to make # chunks = #
341: IO nodes involved, but this author has no access to a BG to figure out how to
342: reliably find the right number. And even then it may or may not be enough.
343: */
344: if (avg_local_vec_size > max_chunk_size) {
345: /* check if we can just split local z-axis: is that enough? */
346: zslices = PetscCeilInt(vec_size, da->p * max_chunk_size) * zslices;
347: if (zslices > da->P) {
348: /* lattice is too large in xy-directions, splitting z only is not enough */
349: zslices = da->P;
350: yslices = PetscCeilInt(vec_size, zslices * da->n * max_chunk_size) * yslices;
351: if (yslices > da->N) {
352: /* lattice is too large in x-direction, splitting along z, y is not enough */
353: yslices = da->N;
354: xslices = PetscCeilInt(vec_size, zslices * yslices * da->m * max_chunk_size) * xslices;
355: }
356: }
357: dim = 0;
358: if (timestep >= 0) ++dim;
359: /* prefer to split z-axis, even down to planar slices */
360: if (dimension == 3) {
361: chunkDims[dim++] = (hsize_t)da->P / zslices;
362: chunkDims[dim++] = (hsize_t)da->N / yslices;
363: chunkDims[dim++] = (hsize_t)da->M / xslices;
364: } else {
365: /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */
366: chunkDims[dim++] = (hsize_t)da->N / yslices;
367: chunkDims[dim++] = (hsize_t)da->M / xslices;
368: }
369: chunk_size = (hsize_t)PetscMax(1, chunkDims[0]) * PetscMax(1, chunkDims[1]) * PetscMax(1, chunkDims[2]) * PetscMax(1, chunkDims[3]) * PetscMax(1, chunkDims[4]) * PetscMax(1, chunkDims[5]) * sizeof(double);
370: } else {
371: if (target_size < chunk_size) {
372: /* only change the defaults if target_size < chunk_size */
373: dim = 0;
374: if (timestep >= 0) ++dim;
375: /* prefer to split z-axis, even down to planar slices */
376: if (dimension == 3) {
377: /* try splitting the z-axis to core-size bits, i.e. divide chunk size by # comm_size in z-direction */
378: if (target_size >= chunk_size / da->p) {
379: /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */
380: chunkDims[dim] = (hsize_t)PetscCeilInt(da->P, da->p);
381: } else {
382: /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be
383: radical and let everyone write all they've got */
384: chunkDims[dim++] = (hsize_t)PetscCeilInt(da->P, da->p);
385: chunkDims[dim++] = (hsize_t)PetscCeilInt(da->N, da->n);
386: chunkDims[dim++] = (hsize_t)PetscCeilInt(da->M, da->m);
387: }
388: } else {
389: /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */
390: if (target_size >= chunk_size / da->n) {
391: /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */
392: chunkDims[dim] = (hsize_t)PetscCeilInt(da->N, da->n);
393: } else {
394: /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be
395: radical and let everyone write all they've got */
396: chunkDims[dim++] = (hsize_t)PetscCeilInt(da->N, da->n);
397: chunkDims[dim++] = (hsize_t)PetscCeilInt(da->M, da->m);
398: }
399: }
400: chunk_size = (hsize_t)PetscMax(1, chunkDims[0]) * PetscMax(1, chunkDims[1]) * PetscMax(1, chunkDims[2]) * PetscMax(1, chunkDims[3]) * PetscMax(1, chunkDims[4]) * PetscMax(1, chunkDims[5]) * sizeof(double);
401: } else {
402: /* precomputed chunks are fine, we don't need to do anything */
403: }
404: }
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
407: #endif
409: #if defined(PETSC_HAVE_HDF5)
410: static PetscErrorCode VecView_MPI_HDF5_DA(Vec xin, PetscViewer viewer)
411: {
412: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
413: DM dm;
414: DM_DA *da;
415: hid_t filespace; /* file dataspace identifier */
416: hid_t chunkspace; /* chunk dataset property identifier */
417: hid_t dset_id; /* dataset identifier */
418: hid_t memspace; /* memory dataspace identifier */
419: hid_t file_id;
420: hid_t group;
421: hid_t memscalartype; /* scalar type for mem (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
422: hid_t filescalartype; /* scalar type for file (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
423: hsize_t dim;
424: hsize_t maxDims[6] = {0}, dims[6] = {0}, chunkDims[6] = {0}, count[6] = {0}, offset[6] = {0}; /* we depend on these being sane later on */
425: PetscBool timestepping = PETSC_FALSE, dim2, spoutput;
426: PetscInt timestep = PETSC_INT_MIN, dimension;
427: const PetscScalar *x;
428: const char *vecname;
430: PetscFunctionBegin;
431: PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group));
432: PetscCall(PetscViewerHDF5IsTimestepping(viewer, ×tepping));
433: if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, ×tep));
434: PetscCall(PetscViewerHDF5GetBaseDimension2(viewer, &dim2));
435: PetscCall(PetscViewerHDF5GetSPOutput(viewer, &spoutput));
437: PetscCall(VecGetDM(xin, &dm));
438: PetscCheck(dm, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA");
439: da = (DM_DA *)dm->data;
440: PetscCall(DMGetDimension(dm, &dimension));
442: /* Create the dataspace for the dataset.
443: *
444: * dims - holds the current dimensions of the dataset
445: *
446: * maxDims - holds the maximum dimensions of the dataset (unlimited
447: * for the number of time steps with the current dimensions for the
448: * other dimensions; so only additional time steps can be added).
449: *
450: * chunkDims - holds the size of a single time step (required to
451: * permit extending dataset).
452: */
453: dim = 0;
454: if (timestep >= 0) {
455: dims[dim] = timestep + 1;
456: maxDims[dim] = H5S_UNLIMITED;
457: chunkDims[dim] = 1;
458: ++dim;
459: }
460: if (dimension == 3) {
461: PetscCall(PetscHDF5IntCast(da->P, dims + dim));
462: maxDims[dim] = dims[dim];
463: chunkDims[dim] = dims[dim];
464: ++dim;
465: }
466: if (dimension > 1) {
467: PetscCall(PetscHDF5IntCast(da->N, dims + dim));
468: maxDims[dim] = dims[dim];
469: chunkDims[dim] = dims[dim];
470: ++dim;
471: }
472: PetscCall(PetscHDF5IntCast(da->M, dims + dim));
473: maxDims[dim] = dims[dim];
474: chunkDims[dim] = dims[dim];
475: ++dim;
476: if (da->w > 1 || dim2) {
477: PetscCall(PetscHDF5IntCast(da->w, dims + dim));
478: maxDims[dim] = dims[dim];
479: chunkDims[dim] = dims[dim];
480: ++dim;
481: }
482: #if defined(PETSC_USE_COMPLEX)
483: dims[dim] = 2;
484: maxDims[dim] = dims[dim];
485: chunkDims[dim] = dims[dim];
486: ++dim;
487: #endif
489: PetscCall(VecGetHDF5ChunkSize(da, xin, dimension, timestep, chunkDims));
491: PetscCallHDF5Return(filespace, H5Screate_simple, ((int)dim, dims, maxDims));
493: #if defined(PETSC_USE_REAL_SINGLE)
494: memscalartype = H5T_NATIVE_FLOAT;
495: filescalartype = H5T_NATIVE_FLOAT;
496: #elif defined(PETSC_USE_REAL___FLOAT128)
497: #error "HDF5 output with 128 bit floats not supported."
498: #elif defined(PETSC_USE_REAL___FP16)
499: #error "HDF5 output with 16 bit floats not supported."
500: #else
501: memscalartype = H5T_NATIVE_DOUBLE;
502: if (spoutput == PETSC_TRUE) filescalartype = H5T_NATIVE_FLOAT;
503: else filescalartype = H5T_NATIVE_DOUBLE;
504: #endif
506: /* Create the dataset with default properties and close filespace */
507: PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
508: if (!H5Lexists(group, vecname, H5P_DEFAULT)) {
509: /* Create chunk */
510: PetscCallHDF5Return(chunkspace, H5Pcreate, (H5P_DATASET_CREATE));
511: PetscCallHDF5(H5Pset_chunk, (chunkspace, (int)dim, chunkDims));
513: PetscCallHDF5Return(dset_id, H5Dcreate2, (group, vecname, filescalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
514: } else {
515: PetscCallHDF5Return(dset_id, H5Dopen2, (group, vecname, H5P_DEFAULT));
516: PetscCallHDF5(H5Dset_extent, (dset_id, dims));
517: }
518: PetscCallHDF5(H5Sclose, (filespace));
520: /* Each process defines a dataset and writes it to the hyperslab in the file */
521: dim = 0;
522: if (timestep >= 0) {
523: offset[dim] = timestep;
524: ++dim;
525: }
526: if (dimension == 3) PetscCall(PetscHDF5IntCast(da->zs, offset + dim++));
527: if (dimension > 1) PetscCall(PetscHDF5IntCast(da->ys, offset + dim++));
528: PetscCall(PetscHDF5IntCast(da->xs / da->w, offset + dim++));
529: if (da->w > 1 || dim2) offset[dim++] = 0;
530: #if defined(PETSC_USE_COMPLEX)
531: offset[dim++] = 0;
532: #endif
533: dim = 0;
534: if (timestep >= 0) {
535: count[dim] = 1;
536: ++dim;
537: }
538: if (dimension == 3) PetscCall(PetscHDF5IntCast(da->ze - da->zs, count + dim++));
539: if (dimension > 1) PetscCall(PetscHDF5IntCast(da->ye - da->ys, count + dim++));
540: PetscCall(PetscHDF5IntCast((da->xe - da->xs) / da->w, count + dim++));
541: if (da->w > 1 || dim2) PetscCall(PetscHDF5IntCast(da->w, count + dim++));
542: #if defined(PETSC_USE_COMPLEX)
543: count[dim++] = 2;
544: #endif
545: PetscCallHDF5Return(memspace, H5Screate_simple, ((int)dim, count, NULL));
546: PetscCallHDF5Return(filespace, H5Dget_space, (dset_id));
547: PetscCallHDF5(H5Sselect_hyperslab, (filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
549: PetscCall(VecGetArrayRead(xin, &x));
550: PetscCallHDF5(H5Dwrite, (dset_id, memscalartype, memspace, filespace, hdf5->dxpl_id, x));
551: PetscCallHDF5(H5Fflush, (file_id, H5F_SCOPE_GLOBAL));
552: PetscCall(VecRestoreArrayRead(xin, &x));
554: #if defined(PETSC_USE_COMPLEX)
555: {
556: PetscBool tru = PETSC_TRUE;
557: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "complex", PETSC_BOOL, &tru));
558: }
559: #endif
560: if (timestepping) PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "timestepping", PETSC_BOOL, ×tepping));
562: /* Close/release resources */
563: if (group != file_id) PetscCallHDF5(H5Gclose, (group));
564: PetscCallHDF5(H5Sclose, (filespace));
565: PetscCallHDF5(H5Sclose, (memspace));
566: PetscCallHDF5(H5Dclose, (dset_id));
567: PetscCall(PetscInfo(xin, "Wrote Vec object with name %s\n", vecname));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
570: #endif
572: extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec, PetscViewer);
574: #if defined(PETSC_HAVE_MPIIO)
575: static PetscErrorCode DMDAArrayMPIIO(DM da, PetscViewer viewer, Vec xin, PetscBool write)
576: {
577: MPI_File mfdes;
578: PetscMPIInt gsizes[4], lsizes[4], lstarts[4], asiz, dof;
579: MPI_Datatype view;
580: const PetscScalar *array;
581: MPI_Offset off;
582: MPI_Aint ub, ul;
583: PetscInt type, rows, vecrows, tr[2];
584: DM_DA *dd = (DM_DA *)da->data;
585: PetscBool skipheader;
587: PetscFunctionBegin;
588: PetscCall(VecGetSize(xin, &vecrows));
589: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipheader));
590: if (!write) {
591: /* Read vector header. */
592: if (!skipheader) {
593: PetscCall(PetscViewerBinaryRead(viewer, tr, 2, NULL, PETSC_INT));
594: type = tr[0];
595: rows = tr[1];
596: PetscCheck(type == VEC_FILE_CLASSID, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Not vector next in file");
597: PetscCheck(rows == vecrows, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Vector in file not same size as DMDA vector");
598: }
599: } else {
600: tr[0] = VEC_FILE_CLASSID;
601: tr[1] = vecrows;
602: if (!skipheader) PetscCall(PetscViewerBinaryWrite(viewer, tr, 2, PETSC_INT));
603: }
605: PetscCall(PetscMPIIntCast(dd->w, &dof));
606: gsizes[0] = dof;
607: PetscCall(PetscMPIIntCast(dd->M, gsizes + 1));
608: PetscCall(PetscMPIIntCast(dd->N, gsizes + 2));
609: PetscCall(PetscMPIIntCast(dd->P, gsizes + 3));
610: lsizes[0] = dof;
611: PetscCall(PetscMPIIntCast((dd->xe - dd->xs) / dof, lsizes + 1));
612: PetscCall(PetscMPIIntCast(dd->ye - dd->ys, lsizes + 2));
613: PetscCall(PetscMPIIntCast(dd->ze - dd->zs, lsizes + 3));
614: lstarts[0] = 0;
615: PetscCall(PetscMPIIntCast(dd->xs / dof, lstarts + 1));
616: PetscCall(PetscMPIIntCast(dd->ys, lstarts + 2));
617: PetscCall(PetscMPIIntCast(dd->zs, lstarts + 3));
618: PetscCallMPI(MPI_Type_create_subarray((PetscMPIInt)(da->dim + 1), gsizes, lsizes, lstarts, MPI_ORDER_FORTRAN, MPIU_SCALAR, &view));
619: PetscCallMPI(MPI_Type_commit(&view));
621: PetscCall(PetscViewerBinaryGetMPIIODescriptor(viewer, &mfdes));
622: PetscCall(PetscViewerBinaryGetMPIIOOffset(viewer, &off));
623: PetscCallMPI(MPI_File_set_view(mfdes, off, MPIU_SCALAR, view, (char *)"native", MPI_INFO_NULL));
624: PetscCall(VecGetArrayRead(xin, &array));
625: asiz = lsizes[1] * (lsizes[2] > 0 ? lsizes[2] : 1) * (lsizes[3] > 0 ? lsizes[3] : 1) * dof;
626: if (write) PetscCall(MPIU_File_write_all(mfdes, (PetscScalar *)array, asiz, MPIU_SCALAR, MPI_STATUS_IGNORE));
627: else PetscCall(MPIU_File_read_all(mfdes, (PetscScalar *)array, asiz, MPIU_SCALAR, MPI_STATUS_IGNORE));
628: PetscCallMPI(MPI_Type_get_extent(view, &ul, &ub));
629: PetscCall(PetscViewerBinaryAddMPIIOOffset(viewer, ub));
630: PetscCall(VecRestoreArrayRead(xin, &array));
631: PetscCallMPI(MPI_Type_free(&view));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
634: #endif
636: PetscErrorCode VecView_MPI_DA(Vec xin, PetscViewer viewer)
637: {
638: DM da;
639: PetscInt dim;
640: Vec natural;
641: PetscBool isdraw, isvtk, isglvis;
642: #if defined(PETSC_HAVE_HDF5)
643: PetscBool ishdf5;
644: #endif
645: const char *prefix, *name;
646: PetscViewerFormat format;
648: PetscFunctionBegin;
649: PetscCall(VecGetDM(xin, &da));
650: PetscCheck(da, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA");
651: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
652: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
653: #if defined(PETSC_HAVE_HDF5)
654: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
655: #endif
656: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
657: if (isdraw) {
658: PetscCall(DMDAGetInfo(da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
659: if (dim == 1) {
660: PetscCall(VecView_MPI_Draw_DA1d(xin, viewer));
661: } else if (dim == 2) {
662: PetscCall(VecView_MPI_Draw_DA2d(xin, viewer));
663: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Cannot graphically view vector associated with this dimensional DMDA %" PetscInt_FMT, dim);
664: } else if (isvtk) { /* Duplicate the Vec */
665: Vec Y;
666: PetscCall(VecDuplicate(xin, &Y));
667: if (((PetscObject)xin)->name) {
668: /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */
669: PetscCall(PetscObjectSetName((PetscObject)Y, ((PetscObject)xin)->name));
670: }
671: PetscCall(VecCopy(xin, Y));
672: {
673: PetscObject dmvtk;
674: PetscBool compatible, compatibleSet;
675: PetscCall(PetscViewerVTKGetDM(viewer, &dmvtk));
676: if (dmvtk) {
678: PetscCall(DMGetCompatibility(da, (DM)dmvtk, &compatible, &compatibleSet));
679: PetscCheck(compatibleSet && compatible, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_INCOMP, "Cannot confirm compatibility of DMs associated with Vecs viewed in the same VTK file. Check that grids are the same.");
680: }
681: PetscCall(PetscViewerVTKAddField(viewer, (PetscObject)da, DMDAVTKWriteAll, PETSC_DEFAULT, PETSC_VTK_POINT_FIELD, PETSC_FALSE, (PetscObject)Y));
682: }
683: #if defined(PETSC_HAVE_HDF5)
684: } else if (ishdf5) {
685: PetscCall(VecView_MPI_HDF5_DA(xin, viewer));
686: #endif
687: } else if (isglvis) {
688: PetscCall(VecView_GLVis(xin, viewer));
689: } else {
690: #if defined(PETSC_HAVE_MPIIO)
691: PetscBool isbinary, isMPIIO;
693: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
694: if (isbinary) {
695: PetscCall(PetscViewerBinaryGetUseMPIIO(viewer, &isMPIIO));
696: if (isMPIIO) {
697: PetscCall(DMDAArrayMPIIO(da, viewer, xin, PETSC_TRUE));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
700: }
701: #endif
703: /* call viewer on natural ordering */
704: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)xin, &prefix));
705: PetscCall(DMDACreateNaturalVector(da, &natural));
706: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)natural, prefix));
707: PetscCall(DMDAGlobalToNaturalBegin(da, xin, INSERT_VALUES, natural));
708: PetscCall(DMDAGlobalToNaturalEnd(da, xin, INSERT_VALUES, natural));
709: PetscCall(PetscObjectGetName((PetscObject)xin, &name));
710: PetscCall(PetscObjectSetName((PetscObject)natural, name));
712: PetscCall(PetscViewerGetFormat(viewer, &format));
713: if (format == PETSC_VIEWER_BINARY_MATLAB) {
714: /* temporarily remove viewer format so it won't trigger in the VecView() */
715: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_DEFAULT));
716: }
718: ((PetscObject)natural)->donotPetscObjectPrintClassNamePrefixType = PETSC_TRUE;
719: PetscCall(VecView(natural, viewer));
720: ((PetscObject)natural)->donotPetscObjectPrintClassNamePrefixType = PETSC_FALSE;
722: if (format == PETSC_VIEWER_BINARY_MATLAB) {
723: MPI_Comm comm;
724: FILE *info;
725: const char *fieldname;
726: char fieldbuf[256];
727: PetscInt dim, ni, nj, nk, pi, pj, pk, dof, n;
729: /* set the viewer format back into the viewer */
730: PetscCall(PetscViewerPopFormat(viewer));
731: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
732: PetscCall(PetscViewerBinaryGetInfoPointer(viewer, &info));
733: PetscCall(DMDAGetInfo(da, &dim, &ni, &nj, &nk, &pi, &pj, &pk, &dof, NULL, NULL, NULL, NULL, NULL));
734: PetscCall(PetscFPrintf(comm, info, "#--- begin code written by PetscViewerBinary for MATLAB format ---#\n"));
735: PetscCall(PetscFPrintf(comm, info, "#$$ tmp = PetscBinaryRead(fd); \n"));
736: if (dim == 1) PetscCall(PetscFPrintf(comm, info, "#$$ tmp = reshape(tmp,%" PetscInt_FMT ",%" PetscInt_FMT ");\n", dof, ni));
737: if (dim == 2) PetscCall(PetscFPrintf(comm, info, "#$$ tmp = reshape(tmp,%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ");\n", dof, ni, nj));
738: if (dim == 3) PetscCall(PetscFPrintf(comm, info, "#$$ tmp = reshape(tmp,%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ");\n", dof, ni, nj, nk));
740: for (n = 0; n < dof; n++) {
741: PetscCall(DMDAGetFieldName(da, n, &fieldname));
742: if (!fieldname || !fieldname[0]) {
743: PetscCall(PetscSNPrintf(fieldbuf, sizeof fieldbuf, "field%" PetscInt_FMT, n));
744: fieldname = fieldbuf;
745: }
746: if (dim == 1) PetscCall(PetscFPrintf(comm, info, "#$$ Set.%s.%s = squeeze(tmp(%" PetscInt_FMT ",:))';\n", name, fieldname, n + 1));
747: if (dim == 2) PetscCall(PetscFPrintf(comm, info, "#$$ Set.%s.%s = squeeze(tmp(%" PetscInt_FMT ",:,:))';\n", name, fieldname, n + 1));
748: if (dim == 3) PetscCall(PetscFPrintf(comm, info, "#$$ Set.%s.%s = permute(squeeze(tmp(%" PetscInt_FMT ",:,:,:)),[2 1 3]);\n", name, fieldname, n + 1));
749: }
750: PetscCall(PetscFPrintf(comm, info, "#$$ clear tmp; \n"));
751: PetscCall(PetscFPrintf(comm, info, "#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n"));
752: }
754: PetscCall(VecDestroy(&natural));
755: }
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: #if defined(PETSC_HAVE_HDF5)
760: static PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer)
761: {
762: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
763: DM da;
764: int dim, rdim;
765: hsize_t dims[6] = {0}, count[6] = {0}, offset[6] = {0};
766: PetscBool dim2 = PETSC_FALSE, timestepping = PETSC_FALSE;
767: PetscInt dimension, timestep = PETSC_INT_MIN, dofInd;
768: PetscScalar *x;
769: const char *vecname;
770: hid_t filespace; /* file dataspace identifier */
771: hid_t dset_id; /* dataset identifier */
772: hid_t memspace; /* memory dataspace identifier */
773: hid_t file_id, group;
774: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
775: DM_DA *dd;
777: PetscFunctionBegin;
778: #if defined(PETSC_USE_REAL_SINGLE)
779: scalartype = H5T_NATIVE_FLOAT;
780: #elif defined(PETSC_USE_REAL___FLOAT128)
781: #error "HDF5 output with 128 bit floats not supported."
782: #elif defined(PETSC_USE_REAL___FP16)
783: #error "HDF5 output with 16 bit floats not supported."
784: #else
785: scalartype = H5T_NATIVE_DOUBLE;
786: #endif
788: PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group));
789: PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
790: PetscCall(PetscViewerHDF5CheckTimestepping_Internal(viewer, vecname));
791: PetscCall(PetscViewerHDF5IsTimestepping(viewer, ×tepping));
792: if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, ×tep));
793: PetscCall(VecGetDM(xin, &da));
794: dd = (DM_DA *)da->data;
795: PetscCall(DMGetDimension(da, &dimension));
797: /* Open dataset */
798: PetscCallHDF5Return(dset_id, H5Dopen2, (group, vecname, H5P_DEFAULT));
800: /* Retrieve the dataspace for the dataset */
801: PetscCallHDF5Return(filespace, H5Dget_space, (dset_id));
802: PetscCallHDF5Return(rdim, H5Sget_simple_extent_dims, (filespace, dims, NULL));
804: /* Expected dimension for holding the dof's */
805: #if defined(PETSC_USE_COMPLEX)
806: dofInd = rdim - 2;
807: #else
808: dofInd = rdim - 1;
809: #endif
811: /* The expected number of dimensions, assuming basedimension2 = false */
812: dim = (int)dimension;
813: if (dd->w > 1) ++dim;
814: if (timestep >= 0) ++dim;
815: #if defined(PETSC_USE_COMPLEX)
816: ++dim;
817: #endif
819: /* In this case the input dataset have one extra, unexpected dimension. */
820: if (rdim == dim + 1) {
821: /* In this case the block size unity */
822: if (dd->w == 1 && dims[dofInd] == 1) dim2 = PETSC_TRUE;
824: /* Special error message for the case where dof does not match the input file */
825: else PetscCheck(dd->w == (PetscInt)dims[dofInd], PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Number of dofs in file is %" PetscInt_FMT ", not %" PetscInt_FMT " as expected", (PetscInt)dims[dofInd], dd->w);
827: /* Other cases where rdim != dim cannot be handled currently */
828: } else PetscCheck(rdim == dim, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file is %d, not %d as expected with dof = %" PetscInt_FMT, rdim, dim, dd->w);
830: /* Set up the hyperslab size */
831: dim = 0;
832: if (timestep >= 0) {
833: offset[dim] = timestep;
834: count[dim] = 1;
835: ++dim;
836: }
837: if (dimension == 3) {
838: PetscCall(PetscHDF5IntCast(dd->zs, offset + dim));
839: PetscCall(PetscHDF5IntCast(dd->ze - dd->zs, count + dim));
840: ++dim;
841: }
842: if (dimension > 1) {
843: PetscCall(PetscHDF5IntCast(dd->ys, offset + dim));
844: PetscCall(PetscHDF5IntCast(dd->ye - dd->ys, count + dim));
845: ++dim;
846: }
847: PetscCall(PetscHDF5IntCast(dd->xs / dd->w, offset + dim));
848: PetscCall(PetscHDF5IntCast((dd->xe - dd->xs) / dd->w, count + dim));
849: ++dim;
850: if (dd->w > 1 || dim2) {
851: offset[dim] = 0;
852: PetscCall(PetscHDF5IntCast(dd->w, count + dim));
853: ++dim;
854: }
855: #if defined(PETSC_USE_COMPLEX)
856: offset[dim] = 0;
857: count[dim] = 2;
858: ++dim;
859: #endif
861: /* Create the memory and filespace */
862: PetscCallHDF5Return(memspace, H5Screate_simple, (dim, count, NULL));
863: PetscCallHDF5(H5Sselect_hyperslab, (filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
865: PetscCall(VecGetArray(xin, &x));
866: PetscCallHDF5(H5Dread, (dset_id, scalartype, memspace, filespace, hdf5->dxpl_id, x));
867: PetscCall(VecRestoreArray(xin, &x));
869: /* Close/release resources */
870: if (group != file_id) PetscCallHDF5(H5Gclose, (group));
871: PetscCallHDF5(H5Sclose, (filespace));
872: PetscCallHDF5(H5Sclose, (memspace));
873: PetscCallHDF5(H5Dclose, (dset_id));
874: PetscFunctionReturn(PETSC_SUCCESS);
875: }
876: #endif
878: static PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer)
879: {
880: DM da;
881: Vec natural;
882: const char *prefix;
883: PetscInt bs;
884: PetscBool flag;
885: DM_DA *dd;
886: #if defined(PETSC_HAVE_MPIIO)
887: PetscBool isMPIIO;
888: #endif
890: PetscFunctionBegin;
891: PetscCall(VecGetDM(xin, &da));
892: dd = (DM_DA *)da->data;
893: #if defined(PETSC_HAVE_MPIIO)
894: PetscCall(PetscViewerBinaryGetUseMPIIO(viewer, &isMPIIO));
895: if (isMPIIO) {
896: PetscCall(DMDAArrayMPIIO(da, viewer, xin, PETSC_FALSE));
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
899: #endif
901: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)xin, &prefix));
902: PetscCall(DMDACreateNaturalVector(da, &natural));
903: PetscCall(PetscObjectSetName((PetscObject)natural, ((PetscObject)xin)->name));
904: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)natural, prefix));
905: PetscCall(VecLoad(natural, viewer));
906: PetscCall(DMDANaturalToGlobalBegin(da, natural, INSERT_VALUES, xin));
907: PetscCall(DMDANaturalToGlobalEnd(da, natural, INSERT_VALUES, xin));
908: PetscCall(VecDestroy(&natural));
909: PetscCall(PetscInfo(xin, "Loading vector from natural ordering into DMDA\n"));
910: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)xin)->prefix, "-vecload_block_size", &bs, &flag));
911: if (flag && bs != dd->w) PetscCall(PetscInfo(xin, "Block size in file %" PetscInt_FMT " not equal to DMDA's dof %" PetscInt_FMT "\n", bs, dd->w));
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer)
916: {
917: DM da;
918: PetscBool isbinary;
919: #if defined(PETSC_HAVE_HDF5)
920: PetscBool ishdf5;
921: #endif
923: PetscFunctionBegin;
924: PetscCall(VecGetDM(xin, &da));
925: PetscCheck(da, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA");
927: #if defined(PETSC_HAVE_HDF5)
928: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
929: #endif
930: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
932: if (isbinary) {
933: PetscCall(VecLoad_Binary_DA(xin, viewer));
934: #if defined(PETSC_HAVE_HDF5)
935: } else if (ishdf5) {
936: PetscCall(VecLoad_HDF5_DA(xin, viewer));
937: #endif
938: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name);
939: PetscFunctionReturn(PETSC_SUCCESS);
940: }