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
main(int argc,char ** argv)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, NULL, 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