Actual source code: ex15.c
1: static char help[] = "An example of writing a global Vec from a DMPlex with HDF5 format.\n";
3: #include <petscdmplex.h>
4: #include <petscviewerhdf5.h>
6: int main(int argc, char **argv)
7: {
8: MPI_Comm comm;
9: DM dm;
10: Vec v, nv, rv, coord;
11: PetscBool test_read = PETSC_FALSE, verbose = PETSC_FALSE, flg;
12: PetscViewer hdf5Viewer;
13: PetscInt numFields = 1;
14: PetscInt numBC = 0;
15: PetscInt numComp[1] = {2};
16: PetscInt numDof[3] = {2, 0, 0};
17: PetscInt bcFields[1] = {0};
18: IS bcPoints[1] = {NULL};
19: PetscSection section;
20: PetscReal norm;
21: PetscInt dim;
24: PetscInitialize(&argc, &argv, (char *) 0, help);if (ierr) return ierr;
25: comm = PETSC_COMM_WORLD;
27: PetscOptionsBegin(PETSC_COMM_WORLD,"","Test Options","none");
28: PetscOptionsBool("-test_read","Test reading from the HDF5 file","",PETSC_FALSE,&test_read,NULL);
29: PetscOptionsBool("-verbose","print the Vecs","",PETSC_FALSE,&verbose,NULL);
30: PetscOptionsEnd();
32: DMCreate(comm, &dm);
33: DMSetType(dm, DMPLEX);
34: DMSetFromOptions(dm);
35: DMViewFromOptions(dm, NULL, "-dm_view");
36: DMGetDimension(dm, &dim);
37: numDof[0] = dim;
38: DMSetNumFields(dm, numFields);
39: DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcFields, bcPoints, NULL, NULL, §ion);
40: DMSetLocalSection(dm, section);
41: PetscSectionDestroy(§ion);
42: DMSetUseNatural(dm, PETSC_TRUE);
43: {
44: PetscPartitioner part;
45: DM dmDist;
47: DMPlexGetPartitioner(dm,&part);
48: PetscPartitionerSetFromOptions(part);
49: DMPlexDistribute(dm, 0, NULL, &dmDist);
50: if (dmDist) {
51: DMDestroy(&dm);
52: dm = dmDist;
53: }
54: }
56: DMCreateGlobalVector(dm, &v);
57: PetscObjectSetName((PetscObject) v, "V");
58: DMGetCoordinates(dm, &coord);
59: VecCopy(coord, v);
61: if (verbose) {
62: PetscInt size, bs;
64: VecGetSize(v, &size);
65: VecGetBlockSize(v, &bs);
66: PetscPrintf(PETSC_COMM_WORLD, "==== original V in global ordering. size==%d\tblock size=%d\n", size, bs);
67: VecView(v, PETSC_VIEWER_STDOUT_WORLD);
68: }
70: DMCreateGlobalVector(dm, &nv);
71: PetscObjectSetName((PetscObject) nv, "NV");
72: DMPlexGlobalToNaturalBegin(dm, v, nv);
73: DMPlexGlobalToNaturalEnd(dm, v, nv);
75: if (verbose) {
76: PetscInt size, bs;
78: VecGetSize(nv, &size);
79: VecGetBlockSize(nv, &bs);
80: PetscPrintf(PETSC_COMM_WORLD, "==== V in natural ordering. size==%d\tblock size=%d\n", size, bs);
81: VecView(nv, PETSC_VIEWER_STDOUT_WORLD);
82: }
84: VecViewFromOptions(v, NULL, "-global_vec_view");
86: if (test_read) {
87: DMCreateGlobalVector(dm, &rv);
88: PetscObjectSetName((PetscObject) rv, "V");
89: /* Test native read */
90: PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
91: PetscViewerPushFormat(hdf5Viewer, PETSC_VIEWER_NATIVE);
92: VecLoad(rv, hdf5Viewer);
93: PetscViewerPopFormat(hdf5Viewer);
94: PetscViewerDestroy(&hdf5Viewer);
95: if (verbose) {
96: PetscInt size, bs;
98: VecGetSize(rv, &size);
99: VecGetBlockSize(rv, &bs);
100: PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
101: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
102: }
103: VecEqual(rv, v, &flg);
104: if (flg) {
105: PetscPrintf(PETSC_COMM_WORLD, "V and RV are equal\n");
106: } else {
107: PetscPrintf(PETSC_COMM_WORLD, "V and RV are not equal\n\n");
108: VecAXPY(rv, -1.0, v);
109: VecNorm(rv, NORM_INFINITY, &norm);
110: PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm);
111: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
112: }
113: /* Test raw read */
114: PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
115: VecLoad(rv, hdf5Viewer);
116: PetscViewerDestroy(&hdf5Viewer);
117: if (verbose) {
118: PetscInt size, bs;
120: VecGetSize(rv, &size);
121: VecGetBlockSize(rv, &bs);
122: PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
123: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
124: }
125: VecEqual(rv, nv, &flg);
126: if (flg) {
127: PetscPrintf(PETSC_COMM_WORLD, "NV and RV are equal\n");
128: } else {
129: PetscPrintf(PETSC_COMM_WORLD, "NV and RV are not equal\n\n");
130: VecAXPY(rv, -1.0, v);
131: VecNorm(rv, NORM_INFINITY, &norm);
132: PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm);
133: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
134: }
135: VecDestroy(&rv);
136: }
137: VecDestroy(&nv);
138: VecDestroy(&v);
139: DMDestroy(&dm);
140: PetscFinalize();
141: return ierr;
142: }
144: /*TEST
145: build:
146: requires: triangle hdf5
147: test:
148: suffix: 0
149: requires: triangle hdf5
150: nsize: 2
151: args: -petscpartitioner_type simple -verbose -globaltonatural_sf_view
152: test:
153: suffix: 1
154: requires: triangle hdf5
155: nsize: 2
156: args: -petscpartitioner_type simple -verbose -global_vec_view hdf5:V.h5:native -test_read
158: TEST*/