Actual source code: ex10.c
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
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);
9: The file storage of the vector is independent of the number of processes used.
10: */
12: #include <petscdm.h>
13: #include <petscdmda.h>
14: #include <petscsys.h>
15: #include <petscviewerhdf5.h>
17: static char help[] = "Test to write HDF5 file from PETSc DMDA Vec.\n\n";
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;
33: dx = (xp - xm) / (Nx - 1);
34: dy = (yp - ym) / (Ny - 1);
36: /* Initialize the Petsc context */
38: PetscInitialize(&argc, &argv, (char *)0, help);
39: 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: DMSetFromOptions(da2D);
41: DMSetUp(da2D);
43: /* Set the coordinates */
44: DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
46: /* Declare gauss as a DMDA component */
47: DMCreateGlobalVector(da2D, &gauss);
48: PetscObjectSetName((PetscObject)gauss, "pressure");
50: /* Initialize vector gauss with a constant value (=1) */
51: VecSet(gauss, value);
53: /* Get the coordinates of the corners for each process */
54: DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0);
56: /* Build the gaussian profile (exp(-x^2-y^2)) */
57: 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: DMDAVecRestoreArray(da2D, gauss, &gauss_ptr);
63: /* Create the HDF5 viewer */
64: PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_WRITE, &H5viewer);
65: PetscViewerSetFromOptions(H5viewer);
67: /* Write the H5 file */
68: VecView(gauss, H5viewer);
70: /* Close the viewer */
71: PetscViewerDestroy(&H5viewer);
73: VecDuplicate(gauss, &input);
74: PetscObjectGetName((PetscObject)gauss, &vecname);
75: PetscObjectSetName((PetscObject)input, vecname);
77: /* Create the HDF5 viewer for reading */
78: PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_READ, &H5viewer);
79: PetscViewerSetFromOptions(H5viewer);
80: VecLoad(input, H5viewer);
81: PetscViewerDestroy(&H5viewer);
83: VecAXPY(input, -1.0, gauss);
84: VecNorm(input, NORM_2, &norm);
87: VecDestroy(&input);
88: VecDestroy(&gauss);
89: DMDestroy(&da2D);
90: PetscFinalize();
91: return 0;
92: }
94: /*TEST
96: build:
97: requires: hdf5 !defined(PETSC_USE_CXXCOMPLEX)
99: test:
100: nsize: 4
102: test:
103: nsize: 4
104: suffix: 2
105: args: -viewer_hdf5_base_dimension2
106: output_file: output/ex10_1.out
108: test:
109: nsize: 4
110: suffix: 3
111: args: -viewer_hdf5_sp_output
112: output_file: output/ex10_1.out
114: TEST*/