1 static char help[] = "Tests various 2-dimensional DMDA routines.\n\n";
2
3 #include <petscdm.h>
4 #include <petscdmda.h>
5
main(int argc,char ** argv)6 int main(int argc, char **argv)
7 {
8 PetscMPIInt rank;
9 PetscInt M = 10, N = 8, m = PETSC_DECIDE;
10 PetscInt s = 2, w = 2, n = PETSC_DECIDE, nloc, l, i, j, kk;
11 PetscInt Xs, Xm, Ys, Ym, iloc, *iglobal;
12 const PetscInt *ltog;
13 PetscInt *lx = NULL, *ly = NULL;
14 PetscBool testorder = PETSC_FALSE, flg;
15 DMBoundaryType bx = DM_BOUNDARY_NONE, by = DM_BOUNDARY_NONE;
16 DM da;
17 PetscViewer viewer;
18 Vec local, global;
19 PetscScalar value;
20 DMDAStencilType st = DMDA_STENCIL_BOX;
21 AO ao;
22
23 PetscFunctionBeginUser;
24 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
25 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "", 300, 0, 400, 400, &viewer));
26
27 /* Readoptions */
28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-NX", &M, NULL));
29 PetscCall(PetscOptionsGetInt(NULL, NULL, "-NY", &N, NULL));
30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
31 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
32 PetscCall(PetscOptionsGetInt(NULL, NULL, "-s", &s, NULL));
33 PetscCall(PetscOptionsGetInt(NULL, NULL, "-w", &w, NULL));
34
35 flg = PETSC_FALSE;
36 PetscCall(PetscOptionsGetBool(NULL, NULL, "-xperiodic", &flg, NULL));
37 if (flg) bx = DM_BOUNDARY_PERIODIC;
38 flg = PETSC_FALSE;
39 PetscCall(PetscOptionsGetBool(NULL, NULL, "-yperiodic", &flg, NULL));
40 if (flg) by = DM_BOUNDARY_PERIODIC;
41 flg = PETSC_FALSE;
42 PetscCall(PetscOptionsGetBool(NULL, NULL, "-xghosted", &flg, NULL));
43 if (flg) bx = DM_BOUNDARY_GHOSTED;
44 flg = PETSC_FALSE;
45 PetscCall(PetscOptionsGetBool(NULL, NULL, "-yghosted", &flg, NULL));
46 if (flg) by = DM_BOUNDARY_GHOSTED;
47 flg = PETSC_FALSE;
48 PetscCall(PetscOptionsGetBool(NULL, NULL, "-star", &flg, NULL));
49 if (flg) st = DMDA_STENCIL_STAR;
50 flg = PETSC_FALSE;
51 PetscCall(PetscOptionsGetBool(NULL, NULL, "-box", &flg, NULL));
52 if (flg) st = DMDA_STENCIL_BOX;
53 flg = PETSC_FALSE;
54 PetscCall(PetscOptionsGetBool(NULL, NULL, "-testorder", &testorder, NULL));
55 /*
56 Test putting two nodes in x and y on each processor, exact last processor
57 in x and y gets the rest.
58 */
59 flg = PETSC_FALSE;
60 PetscCall(PetscOptionsGetBool(NULL, NULL, "-distribute", &flg, NULL));
61 if (flg) {
62 PetscCheck(m != PETSC_DECIDE, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must set -m option with -distribute option");
63 PetscCall(PetscMalloc1(m, &lx));
64 for (i = 0; i < m - 1; i++) lx[i] = 4;
65 lx[m - 1] = M - 4 * (m - 1);
66 PetscCheck(n != PETSC_DECIDE, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must set -n option with -distribute option");
67 PetscCall(PetscMalloc1(n, &ly));
68 for (i = 0; i < n - 1; i++) ly[i] = 2;
69 ly[n - 1] = N - 2 * (n - 1);
70 }
71
72 /* Create distributed array and get vectors */
73 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bx, by, st, M, N, m, n, w, s, lx, ly, &da));
74 PetscCall(DMSetFromOptions(da));
75 PetscCall(DMSetUp(da));
76 PetscCall(PetscFree(lx));
77 PetscCall(PetscFree(ly));
78
79 PetscCall(DMView(da, viewer));
80 PetscCall(DMCreateGlobalVector(da, &global));
81 PetscCall(DMCreateLocalVector(da, &local));
82
83 /* Set global vector; send ghost points to local vectors */
84 value = 1;
85 PetscCall(VecSet(global, value));
86 PetscCall(DMGlobalToLocalBegin(da, global, INSERT_VALUES, local));
87 PetscCall(DMGlobalToLocalEnd(da, global, INSERT_VALUES, local));
88
89 /* Scale local vectors according to processor rank; pass to global vector */
90 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
91 value = rank;
92 PetscCall(VecScale(local, value));
93 PetscCall(DMLocalToGlobalBegin(da, local, INSERT_VALUES, global));
94 PetscCall(DMLocalToGlobalEnd(da, local, INSERT_VALUES, global));
95
96 if (!testorder) { /* turn off printing when testing ordering mappings */
97 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nGlobal Vectors:\n"));
98 PetscCall(VecView(global, PETSC_VIEWER_STDOUT_WORLD));
99 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n"));
100 }
101
102 /* Send ghost points to local vectors */
103 PetscCall(DMGlobalToLocalBegin(da, global, INSERT_VALUES, local));
104 PetscCall(DMGlobalToLocalEnd(da, global, INSERT_VALUES, local));
105
106 flg = PETSC_FALSE;
107 PetscCall(PetscOptionsGetBool(NULL, NULL, "-local_print", &flg, NULL));
108 if (flg) {
109 PetscViewer sviewer;
110
111 PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
112 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\nLocal Vector: processor %d\n", rank));
113 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
114 PetscCall(VecView(local, sviewer));
115 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
116 PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
117 }
118
119 /* Tests mappings between application/PETSc orderings */
120 if (testorder) {
121 ISLocalToGlobalMapping ltogm;
122
123 PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
124 PetscCall(ISLocalToGlobalMappingGetSize(ltogm, &nloc));
125 PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, <og));
126 PetscCall(DMDAGetGhostCorners(da, &Xs, &Ys, NULL, &Xm, &Ym, NULL));
127 PetscCall(DMDAGetAO(da, &ao));
128 PetscCall(PetscMalloc1(nloc, &iglobal));
129
130 /* Set iglobal to be global indices for each processor's local and ghost nodes,
131 using the DMDA ordering of grid points */
132 kk = 0;
133 for (j = Ys; j < Ys + Ym; j++) {
134 for (i = Xs; i < Xs + Xm; i++) {
135 iloc = w * ((j - Ys) * Xm + i - Xs);
136 for (l = 0; l < w; l++) iglobal[kk++] = ltog[iloc + l];
137 }
138 }
139
140 /* Map this to the application ordering (which for DMDAs is just the natural ordering
141 that would be used for 1 processor, numbering most rapidly by x, then y) */
142 PetscCall(AOPetscToApplication(ao, nloc, iglobal));
143
144 /* Then map the application ordering back to the PETSc DMDA ordering */
145 PetscCall(AOApplicationToPetsc(ao, nloc, iglobal));
146
147 /* Verify the mappings */
148 kk = 0;
149 for (j = Ys; j < Ys + Ym; j++) {
150 for (i = Xs; i < Xs + Xm; i++) {
151 iloc = w * ((j - Ys) * Xm + i - Xs);
152 for (l = 0; l < w; l++) {
153 if (iglobal[kk] != ltog[iloc + l]) {
154 PetscCall(PetscFPrintf(PETSC_COMM_SELF, stdout, "[%d] Problem with mapping: j=%" PetscInt_FMT ", i=%" PetscInt_FMT ", l=%" PetscInt_FMT ", petsc1=%" PetscInt_FMT ", petsc2=%" PetscInt_FMT "\n", rank, j, i, l, ltog[iloc + l], iglobal[kk]));
155 }
156 kk++;
157 }
158 }
159 }
160 PetscCall(PetscFree(iglobal));
161 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, <og));
162 }
163
164 /* Free memory */
165 PetscCall(PetscViewerDestroy(&viewer));
166 PetscCall(VecDestroy(&local));
167 PetscCall(VecDestroy(&global));
168 PetscCall(DMDestroy(&da));
169
170 PetscCall(PetscFinalize());
171 return 0;
172 }
173
174 /*TEST
175
176 test:
177 nsize: 4
178 args: -nox
179 filter: grep -v -i Object
180 requires: x
181
182 test:
183 suffix: 2
184 args: -testorder -nox
185 requires: x
186 output_file: output/empty.out
187
188 TEST*/
189