Actual source code: ex15.c
petsc-3.9.4 2018-09-11
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 dim = 2;
14: PetscInt numFields = 1;
15: PetscInt numBC = 0;
16: PetscInt numComp[1] = {2};
17: PetscInt numDof[3] = {2, 0, 0};
18: PetscInt bcFields[1] = {0};
19: IS bcPoints[1] = {NULL};
20: PetscSection section;
21: PetscReal norm;
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: PetscOptionsInt("-dim","the dimension of the problem","",2,&dim,NULL);
31: PetscOptionsEnd();
33: DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, &dm);
34: DMGetDimension(dm, &dim);
35: numDof[0] = dim;
36: DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, NULL, NULL, §ion);
37: DMSetDefaultSection(dm, section);
38: PetscSectionDestroy(§ion);
39: DMSetUseNatural(dm, PETSC_TRUE);
40: {
41: PetscPartitioner part;
42: DM dmDist;
44: DMPlexGetPartitioner(dm,&part);
45: PetscPartitionerSetFromOptions(part);
46: DMPlexDistribute(dm, 0, NULL, &dmDist);
47: if (dmDist) {
48: DMDestroy(&dm);
49: dm = dmDist;
50: }
51: }
53: DMCreateGlobalVector(dm, &v);
54: PetscObjectSetName((PetscObject) v, "V");
55: DMGetCoordinates(dm, &coord);
56: VecCopy(coord, v);
58: if (verbose) {
59: PetscInt size, bs;
61: VecGetSize(v, &size);
62: VecGetBlockSize(v, &bs);
63: PetscPrintf(PETSC_COMM_WORLD, "==== original V in global ordering. size==%d\tblock size=%d\n", size, bs);
64: VecView(v, PETSC_VIEWER_STDOUT_WORLD);
65: }
67: DMCreateGlobalVector(dm, &nv);
68: PetscObjectSetName((PetscObject) nv, "NV");
69: DMPlexGlobalToNaturalBegin(dm, v, nv);
70: DMPlexGlobalToNaturalEnd(dm, v, nv);
72: if (verbose) {
73: PetscInt size, bs;
75: VecGetSize(nv, &size);
76: VecGetBlockSize(nv, &bs);
77: PetscPrintf(PETSC_COMM_WORLD, "==== V in natural ordering. size==%d\tblock size=%d\n", size, bs);
78: VecView(nv, PETSC_VIEWER_STDOUT_WORLD);
79: }
81: VecViewFromOptions(v, NULL, "-global_vec_view");
83: if (test_read) {
84: DMCreateGlobalVector(dm, &rv);
85: PetscObjectSetName((PetscObject) rv, "V");
86: /* Test native read */
87: PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
88: PetscViewerPushFormat(hdf5Viewer, PETSC_VIEWER_NATIVE);
89: VecLoad(rv, hdf5Viewer);
90: PetscViewerPopFormat(hdf5Viewer);
91: PetscViewerDestroy(&hdf5Viewer);
92: if (verbose) {
93: PetscInt size, bs;
95: VecGetSize(rv, &size);
96: VecGetBlockSize(rv, &bs);
97: PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
98: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
99: }
100: VecEqual(rv, v, &flg);
101: if (flg) {
102: PetscPrintf(PETSC_COMM_WORLD, "V and RV are equal\n");
103: } else {
104: PetscPrintf(PETSC_COMM_WORLD, "V and RV are not equal\n\n");
105: VecAXPY(rv, -1.0, v);
106: VecNorm(rv, NORM_INFINITY, &norm);
107: PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm);
108: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
109: }
110: /* Test raw read */
111: PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
112: VecLoad(rv, hdf5Viewer);
113: PetscViewerDestroy(&hdf5Viewer);
114: if (verbose) {
115: PetscInt size, bs;
117: VecGetSize(rv, &size);
118: VecGetBlockSize(rv, &bs);
119: PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
120: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
121: }
122: VecEqual(rv, nv, &flg);
123: if (flg) {
124: PetscPrintf(PETSC_COMM_WORLD, "NV and RV are equal\n");
125: } else {
126: PetscPrintf(PETSC_COMM_WORLD, "NV and RV are not equal\n\n");
127: VecAXPY(rv, -1.0, v);
128: VecNorm(rv, NORM_INFINITY, &norm);
129: PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double) norm);
130: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
131: }
132: VecDestroy(&rv);
133: }
134: VecDestroy(&nv);
135: VecDestroy(&v);
136: DMDestroy(&dm);
137: PetscFinalize();
138: return ierr;
139: }
141: /*TEST
142: build:
143: requires: triangle hdf5
144: test:
145: suffix: 0
146: requires: triangle hdf5
147: nsize: 2
148: args: -petscpartitioner_type simple -verbose -globaltonatural_sf_view
149: test:
150: suffix: 1
151: requires: triangle hdf5
152: nsize: 2
153: args: -petscpartitioner_type simple -verbose -global_vec_view hdf5:V.h5:native -test_read
155: TEST*/