Actual source code: ex10.c
petsc-3.4.5 2014-06-29
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 <petscdmda.h>
13: #include <petscsys.h>
14: #include <petscviewerhdf5.h>
16: static char help[] = "Test to write HDF5 file from PETSc DMDA Vec.\n\n";
18: int main(int argc,char **argv)
19: {
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;
29: PetscScalar **gauss_ptr;
31: dx=(xp-xm)/(Nx-1);
32: dy=(yp-ym)/(Ny-1);
34: /* Initialize the Petsc context */
35: PetscInitialize(&argc,&argv,(char*)0,help);
37: /* Build of the DMDA */
38: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,Nx,Ny,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da2D);
40: /* Set the coordinates */
41: DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
43: /* Declare gauss as a DMDA component */
44: DMCreateGlobalVector(da2D,&gauss);
45: PetscObjectSetName((PetscObject) gauss, "pressure");
47: /* Initialize vector gauss with a constant value (=1) */
48: VecSet(gauss,value);
50: /* Get the coordinates of the corners for each process */
51: DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0);
53: /* Build the gaussian profile (exp(-x^2-y^2)) */
54: DMDAVecGetArray(da2D,gauss,&gauss_ptr);
55: for (j=iys; j<iys+iym; j++) {
56: for (i=ixs; i<ixs+ixm; i++) {
57: gauss_ptr[j][i]=exp(-(xm+i*dx)*(xm+i*dx)-(ym+j*dy)*(ym+j*dy));
58: }
59: }
60: DMDAVecRestoreArray(da2D,gauss,&gauss_ptr);
62: /* Create the HDF5 viewer */
63: PetscViewerHDF5Open(PETSC_COMM_WORLD,"gauss.h5",FILE_MODE_WRITE,&H5viewer);
65: /* Write the H5 file */
66: VecView(gauss,H5viewer);
68: /* Cleaning stage */
69: PetscViewerDestroy(&H5viewer);
70: VecDestroy(&gauss);
71: DMDestroy(&da2D);
72: PetscFinalize();
73: return 0;
74: }