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