Actual source code: ex15.c

  1: static char help[] = "Tests VecView() functionality with DMDA objects when using:"
  2:                      "(i) a PetscViewer binary with MPI-IO support; and (ii) when the binary header is skipped.\n\n";

  4: #include <petscdm.h>
  5: #include <petscdmda.h>

  7: #define DMDA_I 5
  8: #define DMDA_J 4
  9: #define DMDA_K 6

 11: const PetscReal dmda_i_val[] = {1.10, 2.3006, 2.32444, 3.44006, 66.9009};
 12: const PetscReal dmda_j_val[] = {0.0, 0.25, 0.5, 0.75};
 13: const PetscReal dmda_k_val[] = {0.0, 1.1, 2.2, 3.3, 4.4, 5.5};

 15: PetscErrorCode MyVecDump(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x)
 16: {
 17:   MPI_Comm    comm;
 18:   PetscViewer viewer;
 19:   PetscBool   ismpiio, isskip;

 21:   PetscFunctionBeginUser;
 22:   PetscCall(PetscObjectGetComm((PetscObject)x, &comm));

 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));

 31:   PetscCall(VecView(x, viewer));

 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"));

 38:   PetscCall(PetscViewerDestroy(&viewer));
 39:   PetscFunctionReturn(PETSC_SUCCESS);
 40: }

 42: PetscErrorCode MyVecLoad(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x)
 43: {
 44:   MPI_Comm    comm;
 45:   PetscViewer viewer;
 46:   PetscBool   ismpiio, isskip;

 48:   PetscFunctionBeginUser;
 49:   PetscCall(PetscObjectGetComm((PetscObject)x, &comm));

 51:   PetscCall(PetscViewerCreate(comm, &viewer));
 52:   PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
 53:   if (skippheader) PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE));
 54:   PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
 55:   if (usempiio) PetscCall(PetscViewerBinarySetUseMPIIO(viewer, PETSC_TRUE));
 56:   PetscCall(PetscViewerFileSetName(viewer, fname));

 58:   PetscCall(VecLoad(x, viewer));

 60:   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &isskip));
 61:   if (isskip) PetscCall(PetscPrintf(comm, "*** PetscViewer[load] skipping header ***\n"));
 62:   PetscCall(PetscViewerBinaryGetUseMPIIO(viewer, &ismpiio));
 63:   if (ismpiio) PetscCall(PetscPrintf(comm, "*** PetscViewer[load] using MPI-IO ***\n"));

 65:   PetscCall(PetscViewerDestroy(&viewer));
 66:   PetscFunctionReturn(PETSC_SUCCESS);
 67: }

 69: PetscErrorCode DMDAVecGenerateEntries(DM dm, Vec a)
 70: {
 71:   PetscScalar ****LA_v;
 72:   PetscInt        i, j, k, l, si, sj, sk, ni, nj, nk, M, N, dof;

 74:   PetscFunctionBeginUser;
 75:   PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
 76:   PetscCall(DMDAGetCorners(dm, &si, &sj, &sk, &ni, &nj, &nk));
 77:   PetscCall(DMDAVecGetArrayDOF(dm, a, &LA_v));
 78:   for (k = sk; k < sk + nk; k++) {
 79:     for (j = sj; j < sj + nj; j++) {
 80:       for (i = si; i < si + ni; i++) {
 81:         PetscScalar test_value_s;

 83:         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));
 84:         for (l = 0; l < dof; l++) LA_v[k][j][i][l] = (PetscScalar)dof * test_value_s + (PetscScalar)l;
 85:       }
 86:     }
 87:   }
 88:   PetscCall(DMDAVecRestoreArrayDOF(dm, a, &LA_v));
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: PetscErrorCode HeaderlessBinaryReadCheck(DM dm, const char name[])
 93: {
 94:   int         fdes;
 95:   PetscScalar buffer[DMDA_I * DMDA_J * DMDA_K * 10];
 96:   PetscInt    len, d, i, j, k, M, N, dof;
 97:   PetscMPIInt rank;
 98:   PetscBool   dataverified = PETSC_TRUE;

100:   PetscFunctionBeginUser;
101:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
102:   PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
103:   len = DMDA_I * DMDA_J * DMDA_K * dof;
104:   if (rank == 0) {
105:     PetscCall(PetscBinaryOpen(name, FILE_MODE_READ, &fdes));
106:     PetscCall(PetscBinaryRead(fdes, buffer, len, NULL, PETSC_SCALAR));
107:     PetscCall(PetscBinaryClose(fdes));

109:     for (k = 0; k < DMDA_K; k++) {
110:       for (j = 0; j < DMDA_J; j++) {
111:         for (i = 0; i < DMDA_I; i++) {
112:           for (d = 0; d < dof; d++) {
113:             PetscScalar v, test_value_s, test_value;
114:             PetscInt    index;

116:             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));
117:             test_value   = (PetscScalar)dof * test_value_s + (PetscScalar)d;

119:             index = dof * (i + j * M + k * M * N) + d;
120:             v     = PetscAbsScalar(test_value - buffer[index]);
121: #if defined(PETSC_USE_COMPLEX)
122:             if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
123:               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));
124:               dataverified = PETSC_FALSE;
125:             }
126: #else
127:             if (PetscRealPart(v) > 1.0e-10) {
128:               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));
129:               dataverified = PETSC_FALSE;
130:             }
131: #endif
132:           }
133:         }
134:       }
135:     }
136:     if (dataverified) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Headerless read of data verified for: %s\n", name));
137:   }
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: PetscErrorCode VecCompare(Vec a, Vec b)
142: {
143:   PetscInt  locmin[2], locmax[2];
144:   PetscReal min[2], max[2];
145:   Vec       ref;

147:   PetscFunctionBeginUser;
148:   PetscCall(VecMin(a, &locmin[0], &min[0]));
149:   PetscCall(VecMax(a, &locmax[0], &max[0]));

151:   PetscCall(VecMin(b, &locmin[1], &min[1]));
152:   PetscCall(VecMax(b, &locmax[1], &max[1]));

154:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecCompare\n"));
155:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(a)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)min[0], locmin[0]));
156:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  max(a)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)max[0], locmax[0]));

158:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(b)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)min[1], locmin[1]));
159:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  max(b)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)max[1], locmax[1]));

161:   PetscCall(VecDuplicate(a, &ref));
162:   PetscCall(VecCopy(a, ref));
163:   PetscCall(VecAXPY(ref, -1.0, b));
164:   PetscCall(VecMin(ref, &locmin[0], &min[0]));
165:   if (PetscAbsReal(min[0]) > 1.0e-10) {
166:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  ERROR: min(a-b) > 1.0e-10\n"));
167:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(a-b) = %+1.10e\n", (double)PetscAbsReal(min[0])));
168:   } else {
169:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(a-b) < 1.0e-10\n"));
170:   }
171:   PetscCall(VecDestroy(&ref));
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: PetscErrorCode TestDMDAVec(PetscBool usempiio)
176: {
177:   DM        dm;
178:   Vec       x_ref, x_test;
179:   PetscBool skipheader = PETSC_TRUE;

181:   PetscFunctionBeginUser;
182:   if (!usempiio) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s\n", PETSC_FUNCTION_NAME));
183:   else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s [using mpi-io]\n", PETSC_FUNCTION_NAME));
184:   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));
185:   PetscCall(DMSetFromOptions(dm));
186:   PetscCall(DMSetUp(dm));

188:   PetscCall(DMCreateGlobalVector(dm, &x_ref));
189:   PetscCall(DMDAVecGenerateEntries(dm, x_ref));

191:   if (!usempiio) {
192:     PetscCall(MyVecDump("dmda.pbvec", skipheader, PETSC_FALSE, x_ref));
193:   } else {
194:     PetscCall(MyVecDump("dmda-mpiio.pbvec", skipheader, PETSC_TRUE, x_ref));
195:   }

197:   PetscCall(DMCreateGlobalVector(dm, &x_test));

199:   if (!usempiio) {
200:     PetscCall(MyVecLoad("dmda.pbvec", skipheader, usempiio, x_test));
201:   } else {
202:     PetscCall(MyVecLoad("dmda-mpiio.pbvec", skipheader, usempiio, x_test));
203:   }

205:   PetscCall(VecCompare(x_ref, x_test));

207:   if (!usempiio) {
208:     PetscCall(HeaderlessBinaryReadCheck(dm, "dmda.pbvec"));
209:   } else {
210:     PetscCall(HeaderlessBinaryReadCheck(dm, "dmda-mpiio.pbvec"));
211:   }
212:   PetscCall(VecDestroy(&x_ref));
213:   PetscCall(VecDestroy(&x_test));
214:   PetscCall(DMDestroy(&dm));
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: }

218: int main(int argc, char **args)
219: {
220:   PetscBool usempiio = PETSC_FALSE;

222:   PetscFunctionBeginUser;
223:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
224:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-usempiio", &usempiio, NULL));
225:   if (!usempiio) {
226:     PetscCall(TestDMDAVec(PETSC_FALSE));
227:   } else {
228: #if defined(PETSC_HAVE_MPIIO)
229:     PetscCall(TestDMDAVec(PETSC_TRUE));
230: #else
231:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n"));
232: #endif
233:   }
234:   PetscCall(PetscFinalize());
235:   return 0;
236: }

238: /*TEST

240:    test:

242:    test:
243:       suffix: 2
244:       nsize: 12

246:    test:
247:       suffix: 3
248:       nsize: 12
249:       requires: defined(PETSC_HAVE_MPIIO)
250:       args: -usempiio

252: TEST*/