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