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