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, PetscCtx 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) PetscCall(VecGetDM(xlocal, &dac));
188: else dac = da;
189: }
191: /*
192: Get local (ghosted) values of vector
193: */
194: PetscCall(DMGlobalToLocalBegin(dac, xin, INSERT_VALUES, xlocal));
195: PetscCall(DMGlobalToLocalEnd(dac, xin, INSERT_VALUES, xlocal));
196: PetscCall(VecGetArrayRead(xlocal, &zctx.v));
198: /*
199: Get coordinates of nodes
200: */
201: PetscCall(DMGetCoordinates(da, &xcoor));
202: if (!xcoor) {
203: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
204: PetscCall(DMGetCoordinates(da, &xcoor));
205: }
207: /*
208: Determine the min and max coordinates in plot
209: */
210: PetscCall(VecStrideMin(xcoor, 0, NULL, &zctx.xmin));
211: PetscCall(VecStrideMax(xcoor, 0, NULL, &zctx.xmax));
212: PetscCall(VecStrideMin(xcoor, 1, NULL, &zctx.ymin));
213: PetscCall(VecStrideMax(xcoor, 1, NULL, &zctx.ymax));
214: PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_contour_axis", &zctx.showaxis, NULL));
215: if (zctx.showaxis) {
216: coors[0] = zctx.xmin - .05 * (zctx.xmax - zctx.xmin);
217: coors[1] = zctx.ymin - .05 * (zctx.ymax - zctx.ymin);
218: coors[2] = zctx.xmax + .05 * (zctx.xmax - zctx.xmin);
219: coors[3] = zctx.ymax + .05 * (zctx.ymax - zctx.ymin);
220: } else {
221: coors[0] = zctx.xmin;
222: coors[1] = zctx.ymin;
223: coors[2] = zctx.xmax;
224: coors[3] = zctx.ymax;
225: }
226: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-draw_coordinates", coors, &ncoors, NULL));
227: 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]));
229: /*
230: Get local ghosted version of coordinates
231: */
232: PetscCall(PetscObjectQuery((PetscObject)da, "GraphicsCoordinateGhosted", (PetscObject *)&xcoorl));
233: if (!xcoorl) {
234: /* create DMDA to get local version of graphics */
235: PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, M, N, zctx.m, zctx.n, 2, 1, lx, ly, &dag));
236: PetscCall(DMSetUp(dag));
237: PetscCall(PetscInfo(dag, "Creating auxiliary DMDA for managing graphics coordinates ghost points\n"));
238: PetscCall(DMCreateLocalVector(dag, &xcoorl));
239: PetscCall(PetscObjectCompose((PetscObject)da, "GraphicsCoordinateGhosted", (PetscObject)xcoorl));
240: PetscCall(PetscObjectDereference((PetscObject)dag));
241: PetscCall(PetscObjectDereference((PetscObject)xcoorl));
242: } else PetscCall(VecGetDM(xcoorl, &dag));
243: PetscCall(DMGlobalToLocalBegin(dag, xcoor, INSERT_VALUES, xcoorl));
244: PetscCall(DMGlobalToLocalEnd(dag, xcoor, INSERT_VALUES, xcoorl));
245: PetscCall(VecGetArrayRead(xcoorl, &zctx.xy));
246: PetscCall(DMDAGetCoordinateName(da, 0, &zctx.name0));
247: PetscCall(DMDAGetCoordinateName(da, 1, &zctx.name1));
249: /*
250: Get information about size of area each processor must do graphics for
251: */
252: PetscCall(DMDAGetInfo(dac, NULL, &M, &N, NULL, NULL, NULL, NULL, &zctx.dof, NULL, &bx, &by, NULL, NULL));
253: PetscCall(DMDAGetGhostCorners(dac, NULL, NULL, NULL, &zctx.m, &zctx.n, NULL));
254: PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_contour_grid", &zctx.showgrid, NULL));
256: PetscCall(DMDASelectFields(da, &ndisplayfields, &displayfields));
257: PetscCall(PetscViewerGetFormat(viewer, &format));
258: PetscCall(PetscOptionsGetBool(NULL, NULL, "-draw_ports", &useports, NULL));
259: if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
260: if (useports) {
261: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
262: PetscCall(PetscDrawCheckResizedWindow(draw));
263: PetscCall(PetscDrawClear(draw));
264: PetscCall(PetscDrawViewPortsCreate(draw, ndisplayfields, &ports));
265: }
267: /*
268: Loop over each field; drawing each in a different window
269: */
270: for (i = 0; i < ndisplayfields; i++) {
271: zctx.k = displayfields[i];
273: /* determine the min and max value in plot */
274: PetscCall(VecStrideMin(xin, zctx.k, NULL, &zctx.min));
275: PetscCall(VecStrideMax(xin, zctx.k, NULL, &zctx.max));
276: if (zctx.k < nbounds) {
277: zctx.min = bounds[2 * zctx.k];
278: zctx.max = bounds[2 * zctx.k + 1];
279: }
280: if (zctx.min == zctx.max) {
281: zctx.min -= 1.e-12;
282: zctx.max += 1.e-12;
283: }
284: PetscCall(PetscInfo(da, "DMDA 2d contour plot min %g max %g\n", (double)zctx.min, (double)zctx.max));
286: if (useports) {
287: PetscCall(PetscDrawViewPortsSet(ports, i));
288: } else {
289: const char *title;
290: PetscCall(PetscViewerDrawGetDraw(viewer, i, &draw));
291: PetscCall(DMDAGetFieldName(da, zctx.k, &title));
292: if (title) PetscCall(PetscDrawSetTitle(draw, title));
293: }
295: PetscCall(PetscDrawGetPopup(draw, &popup));
296: PetscCall(PetscDrawScalePopup(popup, zctx.min, zctx.max));
297: PetscCall(PetscDrawSetCoordinates(draw, coors[0], coors[1], coors[2], coors[3]));
298: PetscCall(PetscDrawZoom(draw, VecView_MPI_Draw_DA2d_Zoom, &zctx));
299: if (!useports) PetscCall(PetscDrawSave(draw));
300: }
301: if (useports) {
302: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
303: PetscCall(PetscDrawSave(draw));
304: }
306: PetscCall(PetscDrawViewPortsDestroy(ports));
307: PetscCall(PetscFree(displayfields));
308: PetscCall(VecRestoreArrayRead(xcoorl, &zctx.xy));
309: PetscCall(VecRestoreArrayRead(xlocal, &zctx.v));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: #if PetscDefined(HAVE_HDF5)
314: static PetscErrorCode VecGetHDF5ChunkSize(DM_DA *da, Vec xin, PetscInt dimension, PetscInt timestep, hsize_t *chunkDims)
315: {
316: PetscMPIInt comm_size;
317: hsize_t chunk_size, target_size, dim;
318: hsize_t vec_size = sizeof(PetscScalar) * da->M * da->N * da->P * da->w;
319: hsize_t avg_local_vec_size, KiB = 1024, MiB = KiB * KiB, GiB = MiB * KiB, min_size = MiB;
320: hsize_t max_chunks = 64 * KiB; /* HDF5 internal limitation */
321: hsize_t max_chunk_size = 4 * GiB; /* HDF5 internal limitation */
322: hsize_t zslices = da->p, yslices = da->n, xslices = da->m;
324: PetscFunctionBegin;
325: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &comm_size));
326: avg_local_vec_size = (hsize_t)PetscCeilInt64(vec_size, comm_size); /* we will attempt to use this as the chunk size */
328: 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))));
329: /* following line uses sizeof(PetscReal) instead of sizeof(PetscScalar) because the last dimension of chunkDims[] captures the 2* when complex numbers are being used */
330: 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);
332: /*
333: if size/rank > max_chunk_size, we need radical measures: even going down to
334: avg_local_vec_size is not enough, so we simply use chunk size of 4 GiB no matter
335: what, composed in the most efficient way possible.
336: N.B. this minimises the number of chunks, which may or may not be the optimal
337: solution. In a BG, for example, the optimal solution is probably to make # chunks = #
338: IO nodes involved, but this author has no access to a BG to figure out how to
339: reliably find the right number. And even then it may or may not be enough.
340: */
341: if (avg_local_vec_size > max_chunk_size) {
342: hsize_t zslices_full = da->P;
343: hsize_t yslices_full = da->N;
345: /* check if we can just split local z-axis: is that enough? */
346: zslices = PetscCeilInt64(vec_size, da->p * max_chunk_size) * zslices;
347: if (zslices > zslices_full) {
348: /* lattice is too large in xy-directions, splitting z only is not enough */
349: zslices = zslices_full;
350: yslices = PetscCeilInt64(vec_size, zslices * da->n * max_chunk_size) * yslices;
351: if (yslices > yslices_full) {
352: /* lattice is too large in x-direction, splitting along z, y is not enough */
353: yslices = yslices_full;
354: xslices = PetscCeilInt64(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 PetscDefined(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 (PetscDefined(USE_COMPLEX)) {
483: dims[dim] = 2;
484: maxDims[dim] = dims[dim];
485: chunkDims[dim] = dims[dim];
486: ++dim;
487: }
489: PetscCall(VecGetHDF5ChunkSize(da, xin, dimension, timestep, chunkDims));
491: PetscCallHDF5Return(filespace, H5Screate_simple, (int)dim, dims, maxDims);
493: #if PetscDefined(USE_REAL_SINGLE)
494: memscalartype = H5T_NATIVE_FLOAT;
495: filescalartype = H5T_NATIVE_FLOAT;
496: #elif PetscDefined(USE_REAL___FLOAT128)
497: #error "HDF5 output with 128 bit floats not supported."
498: #elif PetscDefined(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 (PetscDefined(USE_COMPLEX)) offset[dim++] = 0;
531: dim = 0;
532: if (timestep >= 0) {
533: count[dim] = 1;
534: ++dim;
535: }
536: if (dimension == 3) PetscCall(PetscHDF5IntCast(da->ze - da->zs, count + dim++));
537: if (dimension > 1) PetscCall(PetscHDF5IntCast(da->ye - da->ys, count + dim++));
538: PetscCall(PetscHDF5IntCast((da->xe - da->xs) / da->w, count + dim++));
539: if (da->w > 1 || dim2) PetscCall(PetscHDF5IntCast(da->w, count + dim++));
540: if (PetscDefined(USE_COMPLEX)) count[dim++] = 2;
541: PetscCallHDF5Return(memspace, H5Screate_simple, (int)dim, count, NULL);
542: PetscCallHDF5Return(filespace, H5Dget_space, dset_id);
543: PetscCallHDF5(H5Sselect_hyperslab, filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
545: PetscCall(VecGetArrayRead(xin, &x));
546: PetscCallHDF5(H5Dwrite, dset_id, memscalartype, memspace, filespace, hdf5->dxpl_id, x);
547: PetscCallHDF5(H5Fflush, file_id, H5F_SCOPE_GLOBAL);
548: PetscCall(VecRestoreArrayRead(xin, &x));
550: if (PetscDefined(USE_COMPLEX)) {
551: PetscBool tru = PETSC_TRUE;
552: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "complex", PETSC_BOOL, &tru));
553: }
554: if (timestepping) PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "timestepping", PETSC_BOOL, ×tepping));
556: /* Close/release resources */
557: if (group != file_id) PetscCallHDF5(H5Gclose, group);
558: PetscCallHDF5(H5Sclose, filespace);
559: PetscCallHDF5(H5Sclose, memspace);
560: PetscCallHDF5(H5Dclose, dset_id);
561: PetscCall(PetscInfo(xin, "Wrote Vec object with name %s\n", vecname));
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
564: #endif
566: extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec, PetscViewer);
568: #if PetscDefined(HAVE_MPIIO)
569: static PetscErrorCode DMDAArrayMPIIO(DM da, PetscViewer viewer, Vec xin, PetscBool write)
570: {
571: MPI_File mfdes;
572: PetscMPIInt gsizes[4], lsizes[4], lstarts[4], asiz, dof;
573: MPI_Datatype view;
574: const PetscScalar *array;
575: MPI_Offset off;
576: MPI_Aint ub, ul;
577: PetscInt type, rows, vecrows, tr[2];
578: DM_DA *dd = (DM_DA *)da->data;
579: PetscBool skipheader;
581: PetscFunctionBegin;
582: PetscCall(VecGetSize(xin, &vecrows));
583: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipheader));
584: if (!write) {
585: /* Read vector header. */
586: if (!skipheader) {
587: PetscCall(PetscViewerBinaryRead(viewer, tr, 2, NULL, PETSC_INT));
588: type = tr[0];
589: rows = tr[1];
590: PetscCheck(type == VEC_FILE_CLASSID, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Not vector next in file");
591: PetscCheck(rows == vecrows, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Vector in file not same size as DMDA vector");
592: }
593: } else {
594: tr[0] = VEC_FILE_CLASSID;
595: tr[1] = vecrows;
596: if (!skipheader) PetscCall(PetscViewerBinaryWrite(viewer, tr, 2, PETSC_INT));
597: }
599: PetscCall(PetscMPIIntCast(dd->w, &dof));
600: gsizes[0] = dof;
601: PetscCall(PetscMPIIntCast(dd->M, gsizes + 1));
602: PetscCall(PetscMPIIntCast(dd->N, gsizes + 2));
603: PetscCall(PetscMPIIntCast(dd->P, gsizes + 3));
604: lsizes[0] = dof;
605: PetscCall(PetscMPIIntCast((dd->xe - dd->xs) / dof, lsizes + 1));
606: PetscCall(PetscMPIIntCast(dd->ye - dd->ys, lsizes + 2));
607: PetscCall(PetscMPIIntCast(dd->ze - dd->zs, lsizes + 3));
608: lstarts[0] = 0;
609: PetscCall(PetscMPIIntCast(dd->xs / dof, lstarts + 1));
610: PetscCall(PetscMPIIntCast(dd->ys, lstarts + 2));
611: PetscCall(PetscMPIIntCast(dd->zs, lstarts + 3));
612: PetscCallMPI(MPI_Type_create_subarray((PetscMPIInt)(da->dim + 1), gsizes, lsizes, lstarts, MPI_ORDER_FORTRAN, MPIU_SCALAR, &view));
613: PetscCallMPI(MPI_Type_commit(&view));
615: PetscCall(PetscViewerBinaryGetMPIIODescriptor(viewer, &mfdes));
616: PetscCall(PetscViewerBinaryGetMPIIOOffset(viewer, &off));
617: PetscCallMPI(MPI_File_set_view(mfdes, off, MPIU_SCALAR, view, (char *)"native", MPI_INFO_NULL));
618: PetscCall(VecGetArrayRead(xin, &array));
619: asiz = lsizes[1] * (lsizes[2] > 0 ? lsizes[2] : 1) * (lsizes[3] > 0 ? lsizes[3] : 1) * dof;
620: if (write) PetscCall(MPIU_File_write_all(mfdes, (PetscScalar *)array, asiz, MPIU_SCALAR, MPI_STATUS_IGNORE));
621: else PetscCall(MPIU_File_read_all(mfdes, (PetscScalar *)array, asiz, MPIU_SCALAR, MPI_STATUS_IGNORE));
622: PetscCallMPI(MPI_Type_get_extent(view, &ul, &ub));
623: PetscCall(PetscViewerBinaryAddMPIIOOffset(viewer, ub));
624: PetscCall(VecRestoreArrayRead(xin, &array));
625: PetscCallMPI(MPI_Type_free(&view));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
628: #endif
630: PetscErrorCode VecView_MPI_DA(Vec xin, PetscViewer viewer)
631: {
632: DM da;
633: PetscInt dim;
634: Vec natural;
635: PetscBool isdraw, isvtk, isglvis;
636: #if PetscDefined(HAVE_HDF5)
637: PetscBool ishdf5;
638: #endif
639: const char *prefix, *name;
640: PetscViewerFormat format;
642: PetscFunctionBegin;
643: PetscCall(VecGetDM(xin, &da));
644: PetscCheck(da, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA");
645: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
646: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
647: #if PetscDefined(HAVE_HDF5)
648: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
649: #endif
650: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
651: if (isdraw) {
652: PetscCall(DMDAGetInfo(da, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
653: if (dim == 1) {
654: PetscCall(VecView_MPI_Draw_DA1d(xin, viewer));
655: } else if (dim == 2) {
656: PetscCall(VecView_MPI_Draw_DA2d(xin, viewer));
657: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Cannot graphically view vector associated with this dimensional DMDA %" PetscInt_FMT, dim);
658: } else if (isvtk) { /* Duplicate the Vec */
659: Vec Y;
660: PetscCall(VecDuplicate(xin, &Y));
661: if (((PetscObject)xin)->name) {
662: /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */
663: PetscCall(PetscObjectSetName((PetscObject)Y, ((PetscObject)xin)->name));
664: }
665: PetscCall(VecCopy(xin, Y));
666: {
667: PetscObject dmvtk;
668: PetscBool compatible, compatibleSet;
669: PetscCall(PetscViewerVTKGetDM(viewer, &dmvtk));
670: if (dmvtk) {
672: PetscCall(DMGetCompatibility(da, (DM)dmvtk, &compatible, &compatibleSet));
673: 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.");
674: }
675: PetscCall(PetscViewerVTKAddField(viewer, (PetscObject)da, DMDAVTKWriteAll, PETSC_DEFAULT, PETSC_VTK_POINT_FIELD, PETSC_FALSE, (PetscObject)Y));
676: }
677: #if PetscDefined(HAVE_HDF5)
678: } else if (ishdf5) {
679: PetscCall(VecView_MPI_HDF5_DA(xin, viewer));
680: #endif
681: } else if (isglvis) {
682: PetscCall(VecView_GLVis(xin, viewer));
683: } else {
684: #if PetscDefined(HAVE_MPIIO)
685: PetscBool isbinary, isMPIIO;
687: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
688: if (isbinary) {
689: PetscCall(PetscViewerBinaryGetUseMPIIO(viewer, &isMPIIO));
690: if (isMPIIO) {
691: PetscCall(DMDAArrayMPIIO(da, viewer, xin, PETSC_TRUE));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
694: }
695: #endif
697: /* call viewer on natural ordering */
698: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)xin, &prefix));
699: PetscCall(DMDACreateNaturalVector(da, &natural));
700: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)natural, prefix));
701: PetscCall(DMDAGlobalToNaturalBegin(da, xin, INSERT_VALUES, natural));
702: PetscCall(DMDAGlobalToNaturalEnd(da, xin, INSERT_VALUES, natural));
703: PetscCall(PetscObjectGetName((PetscObject)xin, &name));
704: PetscCall(PetscObjectSetName((PetscObject)natural, name));
706: PetscCall(PetscViewerGetFormat(viewer, &format));
707: if (format == PETSC_VIEWER_BINARY_MATLAB) {
708: /* temporarily remove viewer format so it won't trigger in the VecView() */
709: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_DEFAULT));
710: }
712: ((PetscObject)natural)->donotPetscObjectPrintClassNamePrefixType = PETSC_TRUE;
713: PetscCall(VecView(natural, viewer));
714: ((PetscObject)natural)->donotPetscObjectPrintClassNamePrefixType = PETSC_FALSE;
716: if (format == PETSC_VIEWER_BINARY_MATLAB) {
717: MPI_Comm comm;
718: FILE *info;
719: const char *fieldname;
720: char fieldbuf[256];
721: PetscInt dim, ni, nj, nk, pi, pj, pk, dof, n;
723: /* set the viewer format back into the viewer */
724: PetscCall(PetscViewerPopFormat(viewer));
725: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
726: PetscCall(PetscViewerBinaryGetInfoPointer(viewer, &info));
727: PetscCall(DMDAGetInfo(da, &dim, &ni, &nj, &nk, &pi, &pj, &pk, &dof, NULL, NULL, NULL, NULL, NULL));
728: PetscCall(PetscFPrintf(comm, info, "#--- begin code written by PetscViewerBinary for MATLAB format ---#\n"));
729: PetscCall(PetscFPrintf(comm, info, "#$$ tmp = PetscBinaryRead(fd); \n"));
730: if (dim == 1) PetscCall(PetscFPrintf(comm, info, "#$$ tmp = reshape(tmp,%" PetscInt_FMT ",%" PetscInt_FMT ");\n", dof, ni));
731: if (dim == 2) PetscCall(PetscFPrintf(comm, info, "#$$ tmp = reshape(tmp,%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ");\n", dof, ni, nj));
732: if (dim == 3) PetscCall(PetscFPrintf(comm, info, "#$$ tmp = reshape(tmp,%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ");\n", dof, ni, nj, nk));
734: for (n = 0; n < dof; n++) {
735: PetscCall(DMDAGetFieldName(da, n, &fieldname));
736: if (!fieldname || !fieldname[0]) {
737: PetscCall(PetscSNPrintf(fieldbuf, sizeof(fieldbuf), "field%" PetscInt_FMT, n));
738: fieldname = fieldbuf;
739: }
740: if (dim == 1) PetscCall(PetscFPrintf(comm, info, "#$$ Set.%s.%s = squeeze(tmp(%" PetscInt_FMT ",:))';\n", name, fieldname, n + 1));
741: if (dim == 2) PetscCall(PetscFPrintf(comm, info, "#$$ Set.%s.%s = squeeze(tmp(%" PetscInt_FMT ",:,:))';\n", name, fieldname, n + 1));
742: if (dim == 3) PetscCall(PetscFPrintf(comm, info, "#$$ Set.%s.%s = permute(squeeze(tmp(%" PetscInt_FMT ",:,:,:)),[2 1 3]);\n", name, fieldname, n + 1));
743: }
744: PetscCall(PetscFPrintf(comm, info, "#$$ clear tmp; \n"));
745: PetscCall(PetscFPrintf(comm, info, "#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n"));
746: }
748: PetscCall(VecDestroy(&natural));
749: }
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: #if PetscDefined(HAVE_HDF5)
754: static PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer)
755: {
756: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
757: DM da;
758: int dim, rdim;
759: hsize_t dims[6] = {0}, count[6] = {0}, offset[6] = {0};
760: PetscBool dim2 = PETSC_FALSE, timestepping = PETSC_FALSE;
761: PetscInt dimension, timestep = PETSC_INT_MIN, dofInd;
762: PetscScalar *x;
763: const char *vecname;
764: hid_t filespace; /* file dataspace identifier */
765: hid_t dset_id; /* dataset identifier */
766: hid_t memspace; /* memory dataspace identifier */
767: hid_t file_id, group;
768: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
769: DM_DA *dd;
771: PetscFunctionBegin;
772: #if PetscDefined(USE_REAL_SINGLE)
773: scalartype = H5T_NATIVE_FLOAT;
774: #elif PetscDefined(USE_REAL___FLOAT128)
775: #error "HDF5 output with 128 bit floats not supported."
776: #elif PetscDefined(USE_REAL___FP16)
777: #error "HDF5 output with 16 bit floats not supported."
778: #else
779: scalartype = H5T_NATIVE_DOUBLE;
780: #endif
782: PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group));
783: PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
784: PetscCall(PetscViewerHDF5CheckTimestepping_Internal(viewer, vecname));
785: PetscCall(PetscViewerHDF5IsTimestepping(viewer, ×tepping));
786: if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, ×tep));
787: PetscCall(VecGetDM(xin, &da));
788: dd = (DM_DA *)da->data;
789: PetscCall(DMGetDimension(da, &dimension));
791: /* Open dataset */
792: PetscCallHDF5Return(dset_id, H5Dopen2, group, vecname, H5P_DEFAULT);
794: /* Retrieve the dataspace for the dataset */
795: PetscCallHDF5Return(filespace, H5Dget_space, dset_id);
796: PetscCallHDF5Return(rdim, H5Sget_simple_extent_dims, filespace, dims, NULL);
798: /* Expected dimension for holding the dof's */
799: dofInd = PetscDefined(USE_COMPLEX) ? rdim - 2 : rdim - 1;
801: /* The expected number of dimensions, assuming basedimension2 = false */
802: dim = (int)dimension;
803: if (dd->w > 1) ++dim;
804: if (timestep >= 0) ++dim;
805: if (PetscDefined(USE_COMPLEX)) ++dim;
807: /* In this case the input dataset have one extra, unexpected dimension. */
808: if (rdim == dim + 1) {
809: /* In this case the block size unity */
810: if (dd->w == 1 && dims[dofInd] == 1) dim2 = PETSC_TRUE;
812: /* Special error message for the case where dof does not match the input file */
813: 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);
815: /* Other cases where rdim != dim cannot be handled currently */
816: } 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);
818: /* Set up the hyperslab size */
819: dim = 0;
820: if (timestep >= 0) {
821: offset[dim] = timestep;
822: count[dim] = 1;
823: ++dim;
824: }
825: if (dimension == 3) {
826: PetscCall(PetscHDF5IntCast(dd->zs, offset + dim));
827: PetscCall(PetscHDF5IntCast(dd->ze - dd->zs, count + dim));
828: ++dim;
829: }
830: if (dimension > 1) {
831: PetscCall(PetscHDF5IntCast(dd->ys, offset + dim));
832: PetscCall(PetscHDF5IntCast(dd->ye - dd->ys, count + dim));
833: ++dim;
834: }
835: PetscCall(PetscHDF5IntCast(dd->xs / dd->w, offset + dim));
836: PetscCall(PetscHDF5IntCast((dd->xe - dd->xs) / dd->w, count + dim));
837: ++dim;
838: if (dd->w > 1 || dim2) {
839: offset[dim] = 0;
840: PetscCall(PetscHDF5IntCast(dd->w, count + dim));
841: ++dim;
842: }
843: if (PetscDefined(USE_COMPLEX)) {
844: offset[dim] = 0;
845: count[dim] = 2;
846: ++dim;
847: }
849: /* Create the memory and filespace */
850: PetscCallHDF5Return(memspace, H5Screate_simple, dim, count, NULL);
851: PetscCallHDF5(H5Sselect_hyperslab, filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
853: PetscCall(VecGetArray(xin, &x));
854: PetscCallHDF5(H5Dread, dset_id, scalartype, memspace, filespace, hdf5->dxpl_id, x);
855: PetscCall(VecRestoreArray(xin, &x));
857: /* Close/release resources */
858: if (group != file_id) PetscCallHDF5(H5Gclose, group);
859: PetscCallHDF5(H5Sclose, filespace);
860: PetscCallHDF5(H5Sclose, memspace);
861: PetscCallHDF5(H5Dclose, dset_id);
862: PetscFunctionReturn(PETSC_SUCCESS);
863: }
864: #endif
866: static PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer)
867: {
868: DM da;
869: Vec natural;
870: const char *prefix;
871: PetscInt bs;
872: PetscBool flag;
873: DM_DA *dd;
874: #if PetscDefined(HAVE_MPIIO)
875: PetscBool isMPIIO;
876: #endif
878: PetscFunctionBegin;
879: PetscCall(VecGetDM(xin, &da));
880: dd = (DM_DA *)da->data;
881: #if PetscDefined(HAVE_MPIIO)
882: PetscCall(PetscViewerBinaryGetUseMPIIO(viewer, &isMPIIO));
883: if (isMPIIO) {
884: PetscCall(DMDAArrayMPIIO(da, viewer, xin, PETSC_FALSE));
885: PetscFunctionReturn(PETSC_SUCCESS);
886: }
887: #endif
889: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)xin, &prefix));
890: PetscCall(DMDACreateNaturalVector(da, &natural));
891: PetscCall(PetscObjectSetName((PetscObject)natural, ((PetscObject)xin)->name));
892: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)natural, prefix));
893: PetscCall(VecLoad(natural, viewer));
894: PetscCall(DMDANaturalToGlobalBegin(da, natural, INSERT_VALUES, xin));
895: PetscCall(DMDANaturalToGlobalEnd(da, natural, INSERT_VALUES, xin));
896: PetscCall(VecDestroy(&natural));
897: PetscCall(PetscInfo(xin, "Loading vector from natural ordering into DMDA\n"));
898: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)xin)->prefix, "-vecload_block_size", &bs, &flag));
899: 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));
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: PetscErrorCode VecLoad_Default_DA(Vec xin, PetscViewer viewer)
904: {
905: DM da;
906: PetscBool isbinary;
907: #if PetscDefined(HAVE_HDF5)
908: PetscBool ishdf5;
909: #endif
911: PetscFunctionBegin;
912: PetscCall(VecGetDM(xin, &da));
913: PetscCheck(da, PetscObjectComm((PetscObject)xin), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMDA");
915: #if PetscDefined(HAVE_HDF5)
916: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
917: #endif
918: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
920: if (isbinary) {
921: PetscCall(VecLoad_Binary_DA(xin, viewer));
922: #if PetscDefined(HAVE_HDF5)
923: } else if (ishdf5) {
924: PetscCall(VecLoad_HDF5_DA(xin, viewer));
925: #endif
926: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name);
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }