1 static char help[] = "An example of writing a global Vec from a DMPlex with HDF5 format.\n";
2
3 #include <petscdmplex.h>
4 #include <petscviewerhdf5.h>
5
main(int argc,char ** argv)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;
22
23 PetscFunctionBeginUser;
24 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
25 comm = PETSC_COMM_WORLD;
26
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();
31
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;
47
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 }
56
57 PetscCall(DMCreateGlobalVector(dm, &v));
58 PetscCall(PetscObjectSetName((PetscObject)v, "V"));
59 PetscCall(DMGetCoordinates(dm, &coord));
60 PetscCall(VecCopy(coord, v));
61
62 if (verbose) {
63 PetscInt size, bs;
64
65 PetscCall(VecGetSize(v, &size));
66 PetscCall(VecGetBlockSize(v, &bs));
67 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== original V in global ordering. size==%" PetscInt_FMT "\tblock size=%" PetscInt_FMT "\n", size, bs));
68 PetscCall(VecView(v, PETSC_VIEWER_STDOUT_WORLD));
69 }
70
71 PetscCall(DMPlexCreateNaturalVector(dm, &nv));
72 PetscCall(PetscObjectSetName((PetscObject)nv, "NV"));
73 PetscCall(DMPlexGlobalToNaturalBegin(dm, v, nv));
74 PetscCall(DMPlexGlobalToNaturalEnd(dm, v, nv));
75
76 if (verbose) {
77 PetscInt size, bs;
78
79 PetscCall(VecGetSize(nv, &size));
80 PetscCall(VecGetBlockSize(nv, &bs));
81 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== V in natural ordering. size==%" PetscInt_FMT "\tblock size=%" PetscInt_FMT "\n", size, bs));
82 PetscCall(VecView(nv, PETSC_VIEWER_STDOUT_WORLD));
83 }
84
85 PetscCall(VecViewFromOptions(v, NULL, "-global_vec_view"));
86
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;
98
99 PetscCall(VecGetSize(rv, &size));
100 PetscCall(VecGetBlockSize(rv, &bs));
101 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%" PetscInt_FMT "\tblock size=%" PetscInt_FMT "\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;
120
121 PetscCall(VecGetSize(rv, &size));
122 PetscCall(VecGetBlockSize(rv, &bs));
123 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "==== Vector from file. size==%" PetscInt_FMT "\tblock size=%" PetscInt_FMT "\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 }
144
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
158
159 TEST*/
160