1 static char help[] = "Tests VecView()/VecLoad() for DMDA vectors (this tests DMDAGlobalToNatural()).\n\n"; 2 3 #include <petscdm.h> 4 #include <petscdmda.h> 5 #include <petscviewerhdf5.h> 6 7 int main(int argc, char **argv) 8 { 9 PetscMPIInt rank, size; 10 PetscInt N = 6, M = 8, P = 5, dof = 1; 11 PetscInt stencil_width = 1, pt = 0, st = 0; 12 PetscBool flg2, flg3, isbinary, mpiio; 13 DMBoundaryType bx = DM_BOUNDARY_NONE, by = DM_BOUNDARY_NONE, bz = DM_BOUNDARY_NONE; 14 DMDAStencilType stencil_type = DMDA_STENCIL_STAR; 15 DM da, da2; 16 Vec global1, global2; 17 PetscScalar mone = -1.0; 18 PetscReal norm; 19 PetscViewer viewer; 20 PetscRandom rdm; 21 #if defined(PETSC_HAVE_HDF5) 22 PetscBool ishdf5; 23 #endif 24 25 PetscFunctionBeginUser; 26 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 27 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 28 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 29 30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 31 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL)); 32 PetscCall(PetscOptionsGetInt(NULL, NULL, "-P", &P, NULL)); 33 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL)); 34 PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &stencil_width, NULL)); 35 PetscCall(PetscOptionsGetInt(NULL, NULL, "-periodic", &pt, NULL)); 36 if (pt == 1) bx = DM_BOUNDARY_PERIODIC; 37 if (pt == 2) by = DM_BOUNDARY_PERIODIC; 38 if (pt == 4) { 39 bx = DM_BOUNDARY_PERIODIC; 40 by = DM_BOUNDARY_PERIODIC; 41 } 42 43 PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_type", &st, NULL)); 44 stencil_type = (DMDAStencilType)st; 45 46 PetscCall(PetscOptionsHasName(NULL, NULL, "-oned", &flg2)); 47 PetscCall(PetscOptionsHasName(NULL, NULL, "-twod", &flg2)); 48 PetscCall(PetscOptionsHasName(NULL, NULL, "-threed", &flg3)); 49 50 PetscCall(PetscOptionsHasName(NULL, NULL, "-binary", &isbinary)); 51 #if defined(PETSC_HAVE_HDF5) 52 PetscCall(PetscOptionsHasName(NULL, NULL, "-hdf5", &ishdf5)); 53 #endif 54 PetscCall(PetscOptionsHasName(NULL, NULL, "-mpiio", &mpiio)); 55 if (flg2) { 56 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bx, by, stencil_type, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, 0, 0, &da)); 57 } else if (flg3) { 58 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, stencil_type, M, N, P, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, 0, 0, 0, &da)); 59 } else { 60 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, bx, M, dof, stencil_width, NULL, &da)); 61 } 62 PetscCall(DMSetFromOptions(da)); 63 PetscCall(DMSetUp(da)); 64 65 PetscCall(DMCreateGlobalVector(da, &global1)); 66 PetscCall(PetscObjectSetName((PetscObject)global1, "Test_Vec")); 67 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm)); 68 PetscCall(PetscRandomSetFromOptions(rdm)); 69 PetscCall(VecSetRandom(global1, rdm)); 70 if (isbinary) { 71 if (mpiio) PetscCall(PetscOptionsSetValue(NULL, "-viewer_binary_mpiio", "")); 72 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_WRITE, &viewer)); 73 #if defined(PETSC_HAVE_HDF5) 74 } else if (ishdf5) { 75 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "temp", FILE_MODE_WRITE, &viewer)); 76 #endif 77 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Invalid Viewer : Run with -binary or -hdf5 option"); 78 PetscCall(VecView(global1, viewer)); 79 PetscCall(PetscViewerDestroy(&viewer)); 80 81 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Global vector written to temp file is \n")); 82 PetscCall(VecView(global1, PETSC_VIEWER_STDOUT_WORLD)); 83 84 if (flg2) { 85 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bx, by, stencil_type, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, 0, 0, &da2)); 86 } else if (flg3) { 87 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, stencil_type, M, N, P, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, 0, 0, 0, &da2)); 88 } else { 89 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, bx, M, dof, stencil_width, NULL, &da2)); 90 } 91 PetscCall(DMSetFromOptions(da2)); 92 PetscCall(DMSetUp(da2)); 93 94 if (isbinary) { 95 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer)); 96 #if defined(PETSC_HAVE_HDF5) 97 } else if (ishdf5) { 98 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer)); 99 #endif 100 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Invalid Viewer : Run with -binary or -hdf5 option"); 101 102 PetscCall(DMCreateGlobalVector(da2, &global2)); 103 PetscCall(PetscObjectSetName((PetscObject)global2, "Test_Vec")); 104 PetscCall(VecLoad(global2, viewer)); 105 PetscCall(PetscViewerDestroy(&viewer)); 106 107 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Global vector read from temp file is \n")); 108 PetscCall(VecView(global2, PETSC_VIEWER_STDOUT_WORLD)); 109 PetscCall(VecAXPY(global2, mone, global1)); 110 PetscCall(VecNorm(global2, NORM_MAX, &norm)); 111 if (norm != 0.0) { 112 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex23: Norm of difference %g should be zero\n", (double)norm)); 113 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Number of processors %d\n", size)); 114 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " M,N,P,dof %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N, P, dof)); 115 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " stencil_width %" PetscInt_FMT " stencil_type %d periodic %d\n", stencil_width, (int)stencil_type, (int)pt)); 116 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " dimension %d\n", 1 + (int)flg2 + (int)flg3)); 117 } 118 119 PetscCall(PetscRandomDestroy(&rdm)); 120 PetscCall(DMDestroy(&da)); 121 PetscCall(DMDestroy(&da2)); 122 PetscCall(VecDestroy(&global1)); 123 PetscCall(VecDestroy(&global2)); 124 PetscCall(PetscFinalize()); 125 return 0; 126 } 127