xref: /petsc/src/dm/tests/ex23.c (revision d5b43468fb8780a8feea140ccd6fa3e6a50411cc)
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 
7 int main(int argc, char **argv)
8 {
9   PetscMPIInt     size;
10   PetscInt        N = 6, m = PETSC_DECIDE, n = PETSC_DECIDE, p = PETSC_DECIDE, M = 8, dof = 1, stencil_width = 1, P = 5, pt = 0, st = 0;
11   PetscBool       flg2, flg3, native = PETSC_FALSE;
12   DMBoundaryType  bx = DM_BOUNDARY_NONE, by = DM_BOUNDARY_NONE, bz = DM_BOUNDARY_NONE;
13   DMDAStencilType stencil_type = DMDA_STENCIL_STAR;
14   DM              da;
15   Vec             global1, global2, global3, global4;
16   PetscScalar     mone = -1.0;
17   PetscReal       norm;
18   PetscViewer     viewer;
19   PetscRandom     rdm;
20 
21   PetscFunctionBeginUser;
22   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
23   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
24   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
25   PetscCall(PetscOptionsGetInt(NULL, NULL, "-P", &P, NULL));
26   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
27   PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &stencil_width, NULL));
28   PetscCall(PetscOptionsGetInt(NULL, NULL, "-periodic", &pt, NULL));
29   PetscCall(PetscOptionsGetBool(NULL, NULL, "-native", &native, NULL));
30   if (pt == 1) bx = DM_BOUNDARY_PERIODIC;
31   if (pt == 2) by = DM_BOUNDARY_PERIODIC;
32   if (pt == 3) {
33     bx = DM_BOUNDARY_PERIODIC;
34     by = DM_BOUNDARY_PERIODIC;
35   }
36   if (pt == 4) bz = DM_BOUNDARY_PERIODIC;
37 
38   PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_type", &st, NULL));
39   stencil_type = (DMDAStencilType)st;
40 
41   PetscCall(PetscOptionsHasName(NULL, NULL, "-one", &flg2));
42   PetscCall(PetscOptionsHasName(NULL, NULL, "-two", &flg2));
43   PetscCall(PetscOptionsHasName(NULL, NULL, "-three", &flg3));
44   if (flg2) {
45     PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bx, by, stencil_type, M, N, m, n, dof, stencil_width, 0, 0, &da));
46   } else if (flg3) {
47     PetscCall(DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, stencil_type, M, N, P, m, n, p, dof, stencil_width, 0, 0, 0, &da));
48   } else {
49     PetscCall(DMDACreate1d(PETSC_COMM_WORLD, bx, M, dof, stencil_width, NULL, &da));
50   }
51   PetscCall(DMSetFromOptions(da));
52   PetscCall(DMSetUp(da));
53 
54   PetscCall(DMCreateGlobalVector(da, &global1));
55   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
56   PetscCall(PetscRandomSetFromOptions(rdm));
57   PetscCall(DMCreateGlobalVector(da, &global2));
58   PetscCall(DMCreateGlobalVector(da, &global3));
59   PetscCall(DMCreateGlobalVector(da, &global4));
60 
61   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_WRITE, &viewer));
62   if (native) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
63   PetscCall(VecSetRandom(global1, rdm));
64   PetscCall(VecView(global1, viewer));
65   PetscCall(VecSetRandom(global3, rdm));
66   PetscCall(VecView(global3, viewer));
67   if (native) PetscCall(PetscViewerPopFormat(viewer));
68   PetscCall(PetscViewerDestroy(&viewer));
69 
70   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer));
71   if (native) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
72   PetscCall(VecLoad(global2, viewer));
73   PetscCall(VecLoad(global4, viewer));
74   if (native) PetscCall(PetscViewerPopFormat(viewer));
75   PetscCall(PetscViewerDestroy(&viewer));
76 
77   if (native) {
78     Vec       filenative;
79     PetscBool same;
80     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "temp", FILE_MODE_READ, &viewer));
81     PetscCall(DMDACreateNaturalVector(da, &filenative));
82     /* DMDA "natural" Vec does not commandeer VecLoad.  The following load will only work when run on the same process
83      * layout, where as the standard VecView/VecLoad (using DMDA and not PETSC_VIEWER_NATIVE) can be read on a different
84      * number of processors. */
85     PetscCall(VecLoad(filenative, viewer));
86     PetscCall(VecEqual(global2, filenative, &same));
87     if (!same) {
88       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex23: global vector does not match contents of file\n"));
89       PetscCall(VecView(global2, 0));
90       PetscCall(VecView(filenative, 0));
91     }
92     PetscCall(PetscViewerDestroy(&viewer));
93     PetscCall(VecDestroy(&filenative));
94   }
95 
96   PetscCall(VecAXPY(global2, mone, global1));
97   PetscCall(VecNorm(global2, NORM_MAX, &norm));
98   if (norm != 0.0) {
99     PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
100     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex23: Norm of difference %g should be zero\n", (double)norm));
101     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Number of processors %d\n", size));
102     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  M,N,P,dof %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N, P, dof));
103     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  stencil_width %" PetscInt_FMT " stencil_type %d periodic %d\n", stencil_width, (int)stencil_type, (int)bx));
104     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  dimension %d\n", 1 + (int)flg2 + (int)flg3));
105   }
106   PetscCall(VecAXPY(global4, mone, global3));
107   PetscCall(VecNorm(global4, NORM_MAX, &norm));
108   if (norm != 0.0) {
109     PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
110     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex23: Norm of difference %g should be zero\n", (double)norm));
111     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Number of processors %d\n", size));
112     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  M,N,P,dof %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N, P, dof));
113     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  stencil_width %" PetscInt_FMT " stencil_type %d periodic %d\n", stencil_width, (int)stencil_type, (int)bx));
114     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  dimension %d\n", 1 + (int)flg2 + (int)flg3));
115   }
116 
117   PetscCall(PetscRandomDestroy(&rdm));
118   PetscCall(DMDestroy(&da));
119   PetscCall(VecDestroy(&global1));
120   PetscCall(VecDestroy(&global2));
121   PetscCall(VecDestroy(&global3));
122   PetscCall(VecDestroy(&global4));
123   PetscCall(PetscFinalize());
124   return 0;
125 }
126 
127 /*TEST
128 
129    test:
130       nsize: {{1  3}}
131       args: -one -dof {{1 2 3}} -stencil_type {{0 1}}
132 
133    test:
134       suffix: 3
135       nsize: {{2 4}}
136       args: -two -dof {{1 3}} -stencil_type {{0 1}}
137 
138    test:
139       suffix: 4
140       nsize: {{1 4}}
141       args: -three -dof {{2 3}} -stencil_type {{0 1}}
142 
143    test:
144       suffix: 2
145       nsize: 2
146       args: -two -native
147 
148 TEST*/
149