xref: /petsc/src/dm/tutorials/ex9.c (revision 6a98f8dc3f2c9149905a87dc2e9d0fedaf64e09a)
1 static char help[] = "Demonstrates HDF5 vector input/ouput\n\n";
2 
3 /*T
4    Concepts: viewers
5    Concepts: HDF5
6    Processors: n
7 T*/
8 #include <petscsys.h>
9 #include <petscdm.h>
10 #include <petscdmda.h>
11 #include <petscviewerhdf5.h>
12 
13 int main(int argc,char **argv)
14 {
15   PetscErrorCode ierr;
16   PetscViewer    viewer;
17   DM             da;
18   Vec            global,local,global2;
19   PetscMPIInt    rank;
20   PetscBool      flg;
21   PetscInt       ndof;
22 
23   /*
24     Every PETSc routine should begin with the PetscInitialize() routine.
25     argc, argv - These command line arguments are taken to extract the options
26                  supplied to PETSc and options supplied to MPI.
27     help       - When PETSc executable is invoked with the option -help,
28                  it prints the various options that can be applied at
29                  runtime.  The user can use the "help" variable place
30                  additional help messages in this printout.
31   */
32   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
33   /* Get number of DOF's from command line */
34   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"DMDA VecView/VecLoad example","");CHKERRQ(ierr);
35   {
36     ndof = 1;
37     PetscOptionsBoundedInt("-ndof","Number of DOF's in DMDA","",ndof,&ndof,NULL,1);
38   }
39   ierr = PetscOptionsEnd();CHKERRQ(ierr);
40 
41   /* Create a DMDA and an associated vector */
42   ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,100,90,PETSC_DECIDE,PETSC_DECIDE,ndof,1,NULL,NULL,&da);CHKERRQ(ierr);
43   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
44   ierr = DMSetUp(da);CHKERRQ(ierr);
45   ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);
46   ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);
47   ierr = VecSet(global,-1.0);CHKERRQ(ierr);
48   ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
49   ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
50   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
51   ierr = VecScale(local,rank+1);CHKERRQ(ierr);
52   ierr = DMLocalToGlobalBegin(da,local,ADD_VALUES,global);CHKERRQ(ierr);
53   ierr = DMLocalToGlobalEnd(da,local,ADD_VALUES,global);CHKERRQ(ierr);
54 
55   /* Create the HDF5 viewer for writing */
56   ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"hdf5output.h5",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
57   ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr);
58 
59   /* Write the Vec without one extra dimension for BS */
60   ierr = PetscViewerHDF5SetBaseDimension2(viewer, PETSC_FALSE);
61   ierr = PetscObjectSetName((PetscObject) global, "noBsDim");CHKERRQ(ierr);
62   ierr = VecView(global,viewer);CHKERRQ(ierr);
63 
64   /* Write the Vec with one extra, 1-sized, dimension for BS */
65   ierr = PetscViewerHDF5SetBaseDimension2(viewer, PETSC_TRUE);
66   ierr = PetscObjectSetName((PetscObject) global, "bsDim");CHKERRQ(ierr);
67   ierr = VecView(global,viewer);CHKERRQ(ierr);
68 
69   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
70   ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr);
71   ierr = VecDuplicate(global,&global2);CHKERRQ(ierr);
72 
73   /* Create the HDF5 viewer for reading */
74   ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"hdf5output.h5",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
75   ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr);
76 
77   /* Load the Vec without the BS dim and compare */
78   ierr = PetscObjectSetName((PetscObject) global2, "noBsDim");CHKERRQ(ierr);
79   ierr = VecLoad(global2,viewer);CHKERRQ(ierr);
80 
81   ierr = VecEqual(global,global2,&flg);CHKERRQ(ierr);
82   if (!flg) {
83     ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: Vectors are not equal\n");CHKERRQ(ierr);
84   }
85 
86   /* Load the Vec with one extra, 1-sized, BS dim and compare */
87   ierr = PetscObjectSetName((PetscObject) global2, "bsDim");CHKERRQ(ierr);
88   ierr = VecLoad(global2,viewer);CHKERRQ(ierr);
89   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
90 
91   ierr = VecEqual(global,global2,&flg);CHKERRQ(ierr);
92   if (!flg) {
93     ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: Vectors are not equal\n");CHKERRQ(ierr);
94   }
95 
96   /* clean up and exit */
97   ierr = VecDestroy(&local);CHKERRQ(ierr);
98   ierr = VecDestroy(&global);CHKERRQ(ierr);
99   ierr = VecDestroy(&global2);CHKERRQ(ierr);
100   ierr = DMDestroy(&da);CHKERRQ(ierr);
101 
102   ierr = PetscFinalize();
103   return ierr;
104 }
105 
106 
107 /*TEST
108 
109       build:
110          requires: hdf5
111 
112       test:
113          nsize: 4
114 
115       test:
116          nsize: 4
117          suffix: 2
118          args: -ndof 2
119          output_file: output/ex9_1.out
120 
121 TEST*/
122