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);
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: DMPlexDistributeSetDefault(dm, PETSC_FALSE);
35: DMSetFromOptions(dm);
36: DMViewFromOptions(dm, NULL, "-dm_view");
37: DMGetDimension(dm, &dim);
38: numDof[0] = dim;
39: DMSetNumFields(dm, numFields);
40: DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcFields, bcPoints, NULL, NULL, §ion);
41: DMSetLocalSection(dm, section);
42: PetscSectionDestroy(§ion);
43: DMSetUseNatural(dm, PETSC_TRUE);
44: {
45: PetscPartitioner part;
46: DM dmDist;
48: DMPlexGetPartitioner(dm, &part);
49: PetscPartitionerSetFromOptions(part);
50: DMPlexDistribute(dm, 0, NULL, &dmDist);
51: if (dmDist) {
52: DMDestroy(&dm);
53: dm = dmDist;
54: }
55: }
57: DMCreateGlobalVector(dm, &v);
58: PetscObjectSetName((PetscObject)v, "V");
59: DMGetCoordinates(dm, &coord);
60: VecCopy(coord, v);
62: if (verbose) {
63: PetscInt size, bs;
65: VecGetSize(v, &size);
66: VecGetBlockSize(v, &bs);
67: PetscPrintf(PETSC_COMM_WORLD, "==== original V in global ordering. size==%d\tblock size=%d\n", size, bs);
68: VecView(v, PETSC_VIEWER_STDOUT_WORLD);
69: }
71: DMPlexCreateNaturalVector(dm, &nv);
72: PetscObjectSetName((PetscObject)nv, "NV");
73: DMPlexGlobalToNaturalBegin(dm, v, nv);
74: DMPlexGlobalToNaturalEnd(dm, v, nv);
76: if (verbose) {
77: PetscInt size, bs;
79: VecGetSize(nv, &size);
80: VecGetBlockSize(nv, &bs);
81: PetscPrintf(PETSC_COMM_WORLD, "==== V in natural ordering. size==%d\tblock size=%d\n", size, bs);
82: VecView(nv, PETSC_VIEWER_STDOUT_WORLD);
83: }
85: VecViewFromOptions(v, NULL, "-global_vec_view");
87: if (0 && test_read) {
88: DMCreateGlobalVector(dm, &rv);
89: PetscObjectSetName((PetscObject)rv, "V");
90: /* Test native read */
91: PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
92: PetscViewerPushFormat(hdf5Viewer, PETSC_VIEWER_NATIVE);
93: VecLoad(rv, hdf5Viewer);
94: PetscViewerPopFormat(hdf5Viewer);
95: PetscViewerDestroy(&hdf5Viewer);
96: if (verbose) {
97: PetscInt size, bs;
99: VecGetSize(rv, &size);
100: VecGetBlockSize(rv, &bs);
101: PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
102: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
103: }
104: VecEqual(rv, v, &flg);
105: if (flg) {
106: PetscPrintf(PETSC_COMM_WORLD, "V and RV are equal\n");
107: } else {
108: PetscPrintf(PETSC_COMM_WORLD, "V and RV are not equal\n\n");
109: VecAXPY(rv, -1.0, v);
110: VecNorm(rv, NORM_INFINITY, &norm);
111: PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double)norm);
112: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
113: }
114: /* Test raw read */
115: PetscViewerHDF5Open(comm, "V.h5", FILE_MODE_READ, &hdf5Viewer);
116: VecLoad(rv, hdf5Viewer);
117: PetscViewerDestroy(&hdf5Viewer);
118: if (verbose) {
119: PetscInt size, bs;
121: VecGetSize(rv, &size);
122: VecGetBlockSize(rv, &bs);
123: PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%d\tblock size=%d\n", size, bs);
124: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
125: }
126: VecEqual(rv, nv, &flg);
127: if (flg) {
128: PetscPrintf(PETSC_COMM_WORLD, "NV and RV are equal\n");
129: } else {
130: PetscPrintf(PETSC_COMM_WORLD, "NV and RV are not equal\n\n");
131: VecAXPY(rv, -1.0, v);
132: VecNorm(rv, NORM_INFINITY, &norm);
133: PetscPrintf(PETSC_COMM_WORLD, "diff norm is = %g\n", (double)norm);
134: VecView(rv, PETSC_VIEWER_STDOUT_WORLD);
135: }
136: VecDestroy(&rv);
137: }
138: VecDestroy(&nv);
139: VecDestroy(&v);
140: DMDestroy(&dm);
141: PetscFinalize();
142: return 0;
143: }
145: /*TEST
146: build:
147: requires: triangle hdf5
148: test:
149: suffix: 0
150: requires: triangle hdf5
151: nsize: 2
152: args: -petscpartitioner_type simple -verbose -globaltonatural_sf_view
153: test:
154: suffix: 1
155: TODO: broken
156: requires: triangle hdf5
157: nsize: 2
158: args: -petscpartitioner_type simple -verbose -global_vec_view hdf5:V.h5:native -test_read
160: TEST*/