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