xref: /petsc/src/dm/tutorials/ex15.c (revision 38f6737523616bfeabdd577adb30ea0990bf5fe7) 
1 
2 static char help[] = "Tests VecView() functionality with DMDA objects when using:"
3                      "(i) a PetscViewer binary with MPI-IO support; and (ii) when the binary header is skipped.\n\n";
4 
5 #include <petscdm.h>
6 #include <petscdmda.h>
7 
8 #define DMDA_I 5
9 #define DMDA_J 4
10 #define DMDA_K 6
11 
12 const PetscReal dmda_i_val[] = {1.10, 2.3006, 2.32444, 3.44006, 66.9009};
13 const PetscReal dmda_j_val[] = {0.0, 0.25, 0.5, 0.75};
14 const PetscReal dmda_k_val[] = {0.0, 1.1, 2.2, 3.3, 4.4, 5.5};
15 
16 PetscErrorCode MyVecDump(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x) {
17   MPI_Comm    comm;
18   PetscViewer viewer;
19   PetscBool   ismpiio, isskip;
20 
21   PetscFunctionBeginUser;
22   PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
23 
24   PetscCall(PetscViewerCreate(comm, &viewer));
25   PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
26   if (skippheader) PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE));
27   PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
28   if (usempiio) PetscCall(PetscViewerBinarySetUseMPIIO(viewer, PETSC_TRUE));
29   PetscCall(PetscViewerFileSetName(viewer, fname));
30 
31   PetscCall(VecView(x, viewer));
32 
33   PetscCall(PetscViewerBinaryGetUseMPIIO(viewer, &ismpiio));
34   if (ismpiio) PetscCall(PetscPrintf(comm, "*** PetscViewer[write] using MPI-IO ***\n"));
35   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &isskip));
36   if (isskip) PetscCall(PetscPrintf(comm, "*** PetscViewer[write] skipping header ***\n"));
37 
38   PetscCall(PetscViewerDestroy(&viewer));
39   PetscFunctionReturn(0);
40 }
41 
42 PetscErrorCode MyVecLoad(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x) {
43   MPI_Comm    comm;
44   PetscViewer viewer;
45   PetscBool   ismpiio, isskip;
46 
47   PetscFunctionBeginUser;
48   PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
49 
50   PetscCall(PetscViewerCreate(comm, &viewer));
51   PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
52   if (skippheader) PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE));
53   PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
54   if (usempiio) PetscCall(PetscViewerBinarySetUseMPIIO(viewer, PETSC_TRUE));
55   PetscCall(PetscViewerFileSetName(viewer, fname));
56 
57   PetscCall(VecLoad(x, viewer));
58 
59   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &isskip));
60   if (isskip) PetscCall(PetscPrintf(comm, "*** PetscViewer[load] skipping header ***\n"));
61   PetscCall(PetscViewerBinaryGetUseMPIIO(viewer, &ismpiio));
62   if (ismpiio) PetscCall(PetscPrintf(comm, "*** PetscViewer[load] using MPI-IO ***\n"));
63 
64   PetscCall(PetscViewerDestroy(&viewer));
65   PetscFunctionReturn(0);
66 }
67 
68 PetscErrorCode DMDAVecGenerateEntries(DM dm, Vec a) {
69   PetscScalar ****LA_v;
70   PetscInt        i, j, k, l, si, sj, sk, ni, nj, nk, M, N, dof;
71 
72   PetscFunctionBeginUser;
73   PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
74   PetscCall(DMDAGetCorners(dm, &si, &sj, &sk, &ni, &nj, &nk));
75   PetscCall(DMDAVecGetArrayDOF(dm, a, &LA_v));
76   for (k = sk; k < sk + nk; k++) {
77     for (j = sj; j < sj + nj; j++) {
78       for (i = si; i < si + ni; i++) {
79         PetscScalar test_value_s;
80 
81         test_value_s = dmda_i_val[i] * ((PetscScalar)i) + dmda_j_val[j] * ((PetscScalar)(i + j * M)) + dmda_k_val[k] * ((PetscScalar)(i + j * M + k * M * N));
82         for (l = 0; l < dof; l++) { LA_v[k][j][i][l] = (PetscScalar)dof * test_value_s + (PetscScalar)l; }
83       }
84     }
85   }
86   PetscCall(DMDAVecRestoreArrayDOF(dm, a, &LA_v));
87   PetscFunctionReturn(0);
88 }
89 
90 PetscErrorCode HeaderlessBinaryReadCheck(DM dm, const char name[]) {
91   int         fdes;
92   PetscScalar buffer[DMDA_I * DMDA_J * DMDA_K * 10];
93   PetscInt    len, d, i, j, k, M, N, dof;
94   PetscMPIInt rank;
95   PetscBool   dataverified = PETSC_TRUE;
96 
97   PetscFunctionBeginUser;
98   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
99   PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
100   len = DMDA_I * DMDA_J * DMDA_K * dof;
101   if (rank == 0) {
102     PetscCall(PetscBinaryOpen(name, FILE_MODE_READ, &fdes));
103     PetscCall(PetscBinaryRead(fdes, buffer, len, NULL, PETSC_SCALAR));
104     PetscCall(PetscBinaryClose(fdes));
105 
106     for (k = 0; k < DMDA_K; k++) {
107       for (j = 0; j < DMDA_J; j++) {
108         for (i = 0; i < DMDA_I; i++) {
109           for (d = 0; d < dof; d++) {
110             PetscScalar v, test_value_s, test_value;
111             PetscInt    index;
112 
113             test_value_s = dmda_i_val[i] * ((PetscScalar)i) + dmda_j_val[j] * ((PetscScalar)(i + j * M)) + dmda_k_val[k] * ((PetscScalar)(i + j * M + k * M * N));
114             test_value   = (PetscScalar)dof * test_value_s + (PetscScalar)d;
115 
116             index = dof * (i + j * M + k * M * N) + d;
117             v     = PetscAbsScalar(test_value - buffer[index]);
118 #if defined(PETSC_USE_COMPLEX)
119             if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
120               PetscCall(PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference > 1.0e-10 occurred (delta = (%+1.12e,%+1.12e) [loc %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "(%" PetscInt_FMT ")])\n", (double)PetscRealPart(test_value), (double)PetscImaginaryPart(test_value), i, j, k, d));
121               dataverified = PETSC_FALSE;
122             }
123 #else
124             if (PetscRealPart(v) > 1.0e-10) {
125               PetscCall(PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "(%" PetscInt_FMT ")])\n", (double)PetscRealPart(test_value), i, j, k, d));
126               dataverified = PETSC_FALSE;
127             }
128 #endif
129           }
130         }
131       }
132     }
133     if (dataverified) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Headerless read of data verified for: %s\n", name));
134   }
135   PetscFunctionReturn(0);
136 }
137 
138 PetscErrorCode VecCompare(Vec a, Vec b) {
139   PetscInt  locmin[2], locmax[2];
140   PetscReal min[2], max[2];
141   Vec       ref;
142 
143   PetscFunctionBeginUser;
144   PetscCall(VecMin(a, &locmin[0], &min[0]));
145   PetscCall(VecMax(a, &locmax[0], &max[0]));
146 
147   PetscCall(VecMin(b, &locmin[1], &min[1]));
148   PetscCall(VecMax(b, &locmax[1], &max[1]));
149 
150   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecCompare\n"));
151   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(a)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)min[0], locmin[0]));
152   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  max(a)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)max[0], locmax[0]));
153 
154   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(b)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)min[1], locmin[1]));
155   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  max(b)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)max[1], locmax[1]));
156 
157   PetscCall(VecDuplicate(a, &ref));
158   PetscCall(VecCopy(a, ref));
159   PetscCall(VecAXPY(ref, -1.0, b));
160   PetscCall(VecMin(ref, &locmin[0], &min[0]));
161   if (PetscAbsReal(min[0]) > 1.0e-10) {
162     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  ERROR: min(a-b) > 1.0e-10\n"));
163     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(a-b) = %+1.10e\n", (double)PetscAbsReal(min[0])));
164   } else {
165     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(a-b) < 1.0e-10\n"));
166   }
167   PetscCall(VecDestroy(&ref));
168   PetscFunctionReturn(0);
169 }
170 
171 PetscErrorCode TestDMDAVec(PetscBool usempiio) {
172   DM        dm;
173   Vec       x_ref, x_test;
174   PetscBool skipheader = PETSC_TRUE;
175 
176   PetscFunctionBeginUser;
177   if (!usempiio) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s\n", PETSC_FUNCTION_NAME));
178   else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s [using mpi-io]\n", PETSC_FUNCTION_NAME));
179   PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, DMDA_I, DMDA_J, DMDA_K, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 3, 2, NULL, NULL, NULL, &dm));
180   PetscCall(DMSetFromOptions(dm));
181   PetscCall(DMSetUp(dm));
182 
183   PetscCall(DMCreateGlobalVector(dm, &x_ref));
184   PetscCall(DMDAVecGenerateEntries(dm, x_ref));
185 
186   if (!usempiio) {
187     PetscCall(MyVecDump("dmda.pbvec", skipheader, PETSC_FALSE, x_ref));
188   } else {
189     PetscCall(MyVecDump("dmda-mpiio.pbvec", skipheader, PETSC_TRUE, x_ref));
190   }
191 
192   PetscCall(DMCreateGlobalVector(dm, &x_test));
193 
194   if (!usempiio) {
195     PetscCall(MyVecLoad("dmda.pbvec", skipheader, usempiio, x_test));
196   } else {
197     PetscCall(MyVecLoad("dmda-mpiio.pbvec", skipheader, usempiio, x_test));
198   }
199 
200   PetscCall(VecCompare(x_ref, x_test));
201 
202   if (!usempiio) {
203     PetscCall(HeaderlessBinaryReadCheck(dm, "dmda.pbvec"));
204   } else {
205     PetscCall(HeaderlessBinaryReadCheck(dm, "dmda-mpiio.pbvec"));
206   }
207   PetscCall(VecDestroy(&x_ref));
208   PetscCall(VecDestroy(&x_test));
209   PetscCall(DMDestroy(&dm));
210   PetscFunctionReturn(0);
211 }
212 
213 int main(int argc, char **args) {
214   PetscBool usempiio = PETSC_FALSE;
215 
216   PetscFunctionBeginUser;
217   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
218   PetscCall(PetscOptionsGetBool(NULL, NULL, "-usempiio", &usempiio, NULL));
219   if (!usempiio) {
220     PetscCall(TestDMDAVec(PETSC_FALSE));
221   } else {
222 #if defined(PETSC_HAVE_MPIIO)
223     PetscCall(TestDMDAVec(PETSC_TRUE));
224 #else
225     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n"));
226 #endif
227   }
228   PetscCall(PetscFinalize());
229   return 0;
230 }
231 
232 /*TEST
233 
234    test:
235 
236    test:
237       suffix: 2
238       nsize: 12
239 
240    test:
241       suffix: 3
242       nsize: 12
243       requires: defined(PETSC_HAVE_MPIIO)
244       args: -usempiio
245 
246 TEST*/
247