1 /*
2 Demonstrates using the HDF5 viewer with a DMDA Vec
3 - create a global vector containing a gauss profile (exp(-x^2-y^2))
4 - write the global vector in a hdf5 file
5
6 The resulting file gauss.h5 can be viewed with Visit (an open-source visualization package)
7 Or with some versions of MATLAB with data=hdfread('gauss.h5','pressure'); mesh(data);
8
9 The file storage of the vector is independent of the number of processes used.
10 */
11
12 #include <petscdm.h>
13 #include <petscdmda.h>
14 #include <petscsys.h>
15 #include <petscviewerhdf5.h>
16
17 static char help[] = "Test to write HDF5 file from PETSc DMDA Vec.\n\n";
18
main(int argc,char ** argv)19 int main(int argc, char **argv)
20 {
21 DM da2D;
22 PetscInt i, j, ixs, ixm, iys, iym;
23 PetscViewer H5viewer;
24 PetscScalar xm = -1.0, xp = 1.0;
25 PetscScalar ym = -1.0, yp = 1.0;
26 PetscScalar value = 1.0, dx, dy;
27 PetscInt Nx = 40, Ny = 40;
28 Vec gauss, input;
29 PetscScalar **gauss_ptr;
30 PetscReal norm;
31 const char *vecname;
32
33 dx = (xp - xm) / (Nx - 1);
34 dy = (yp - ym) / (Ny - 1);
35
36 /* Initialize the PETSc context */
37 PetscFunctionBeginUser;
38 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
39 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, Nx, Ny, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da2D));
40 PetscCall(DMSetFromOptions(da2D));
41 PetscCall(DMSetUp(da2D));
42
43 /* Set the coordinates */
44 PetscCall(DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
45
46 /* Declare gauss as a DMDA component */
47 PetscCall(DMCreateGlobalVector(da2D, &gauss));
48 PetscCall(PetscObjectSetName((PetscObject)gauss, "pressure"));
49
50 /* Initialize vector gauss with a constant value (=1) */
51 PetscCall(VecSet(gauss, value));
52
53 /* Get the coordinates of the corners for each process */
54 PetscCall(DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0));
55
56 /* Build the gaussian profile (exp(-x^2-y^2)) */
57 PetscCall(DMDAVecGetArray(da2D, gauss, &gauss_ptr));
58 for (j = iys; j < iys + iym; j++) {
59 for (i = ixs; i < ixs + ixm; i++) gauss_ptr[j][i] = PetscExpScalar(-(xm + i * dx) * (xm + i * dx) - (ym + j * dy) * (ym + j * dy));
60 }
61 PetscCall(DMDAVecRestoreArray(da2D, gauss, &gauss_ptr));
62
63 /* Create the HDF5 viewer */
64 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_WRITE, &H5viewer));
65 PetscCall(PetscViewerSetFromOptions(H5viewer));
66
67 /* Write the H5 file */
68 PetscCall(VecView(gauss, H5viewer));
69
70 /* Close the viewer */
71 PetscCall(PetscViewerDestroy(&H5viewer));
72
73 PetscCall(VecDuplicate(gauss, &input));
74 PetscCall(PetscObjectGetName((PetscObject)gauss, &vecname));
75 PetscCall(PetscObjectSetName((PetscObject)input, vecname));
76
77 /* Create the HDF5 viewer for reading */
78 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_READ, &H5viewer));
79 PetscCall(PetscViewerSetFromOptions(H5viewer));
80 PetscCall(VecLoad(input, H5viewer));
81 PetscCall(PetscViewerDestroy(&H5viewer));
82
83 PetscCall(VecAXPY(input, -1.0, gauss));
84 PetscCall(VecNorm(input, NORM_2, &norm));
85 PetscCheck(norm <= 1.e-6, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Vec read in does not match vector written out");
86
87 PetscCall(VecDestroy(&input));
88 PetscCall(VecDestroy(&gauss));
89 PetscCall(DMDestroy(&da2D));
90 PetscCall(PetscFinalize());
91 return 0;
92 }
93
94 /*TEST
95
96 build:
97 requires: hdf5 !defined(PETSC_USE_CXXCOMPLEX)
98
99 test:
100 nsize: 4
101 output_file: output/empty.out
102
103 test:
104 nsize: 4
105 suffix: 2
106 args: -viewer_hdf5_base_dimension2
107 output_file: output/empty.out
108
109 test:
110 nsize: 4
111 suffix: 3
112 args: -viewer_hdf5_sp_output
113 output_file: output/empty.out
114
115 TEST*/
116