1 static char help[] = "Tests DMLocalToLocalxxx() for DMDA.\n\n";
2
3 #include <petscdmda.h>
4
main(int argc,char ** argv)5 int main(int argc, char **argv)
6 {
7 PetscMPIInt rank;
8 PetscInt M = 8, dof = 1, stencil_width = 1, i, start, end, P = 5, N = 6, m = PETSC_DECIDE, n = PETSC_DECIDE, p = PETSC_DECIDE, pt = 0, st = 0;
9 PetscBool flg = PETSC_FALSE, flg2, flg3;
10 DMBoundaryType periodic;
11 DMDAStencilType stencil_type;
12 DM da;
13 Vec local, global, local_copy;
14 PetscScalar value;
15 PetscReal norm, work;
16 PetscViewer viewer;
17 char filename[64];
18 FILE *file;
19
20 PetscFunctionBeginUser;
21 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
22 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &stencil_width, NULL));
26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-periodic", &pt, NULL));
27
28 periodic = (DMBoundaryType)pt;
29
30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_type", &st, NULL));
31
32 stencil_type = (DMDAStencilType)st;
33
34 PetscCall(PetscOptionsHasName(NULL, NULL, "-grid2d", &flg2));
35 PetscCall(PetscOptionsHasName(NULL, NULL, "-grid3d", &flg3));
36 if (flg2) {
37 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, periodic, periodic, stencil_type, M, N, m, n, dof, stencil_width, NULL, NULL, &da));
38 } else if (flg3) {
39 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, periodic, periodic, periodic, stencil_type, M, N, P, m, n, p, dof, stencil_width, NULL, NULL, NULL, &da));
40 } else {
41 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, periodic, M, dof, stencil_width, NULL, &da));
42 }
43 PetscCall(DMSetFromOptions(da));
44 PetscCall(DMSetUp(da));
45
46 PetscCall(DMCreateGlobalVector(da, &global));
47 PetscCall(DMCreateLocalVector(da, &local));
48 PetscCall(VecDuplicate(local, &local_copy));
49
50 /* zero out vectors so that ghostpoints are zero */
51 value = 0;
52 PetscCall(VecSet(local, value));
53 PetscCall(VecSet(local_copy, value));
54
55 PetscCall(VecGetOwnershipRange(global, &start, &end));
56 for (i = start; i < end; i++) {
57 value = i + 1;
58 PetscCall(VecSetValues(global, 1, &i, &value, INSERT_VALUES));
59 }
60 PetscCall(VecAssemblyBegin(global));
61 PetscCall(VecAssemblyEnd(global));
62
63 PetscCall(DMGlobalToLocalBegin(da, global, INSERT_VALUES, local));
64 PetscCall(DMGlobalToLocalEnd(da, global, INSERT_VALUES, local));
65
66 PetscCall(DMLocalToLocalBegin(da, local, INSERT_VALUES, local_copy));
67 PetscCall(DMLocalToLocalEnd(da, local, INSERT_VALUES, local_copy));
68
69 PetscCall(PetscOptionsGetBool(NULL, NULL, "-save", &flg, NULL));
70 if (flg) {
71 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
72 PetscCall(PetscSNPrintf(filename, PETSC_STATIC_ARRAY_LENGTH(filename), "local.%d", rank));
73 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_SELF, filename, &viewer));
74 PetscCall(PetscViewerASCIIGetPointer(viewer, &file));
75 PetscCall(VecView(local, viewer));
76 fprintf(file, "Vector with correct ghost points\n");
77 PetscCall(VecView(local_copy, viewer));
78 PetscCall(PetscViewerDestroy(&viewer));
79 }
80
81 PetscCall(VecAXPY(local_copy, -1.0, local));
82 PetscCall(VecNorm(local_copy, NORM_MAX, &work));
83 PetscCallMPI(MPIU_Allreduce(&work, &norm, 1, MPIU_REAL, MPIU_MAX, PETSC_COMM_WORLD));
84 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of difference %g should be zero\n", (double)norm));
85
86 PetscCall(VecDestroy(&local_copy));
87 PetscCall(VecDestroy(&local));
88 PetscCall(VecDestroy(&global));
89 PetscCall(DMDestroy(&da));
90 PetscCall(PetscFinalize());
91 return 0;
92 }
93
94 /*TEST
95
96 test:
97 nsize: 8
98 args: -dof 3 -stencil_width 2 -M 50 -N 50 -periodic
99
100 test:
101 suffix: 2
102 nsize: 8
103 args: -dof 3 -stencil_width 2 -M 50 -N 50 -periodic -grid2d
104 output_file: output/ex7_1.out
105
106 test:
107 suffix: 3
108 nsize: 8
109 args: -dof 3 -stencil_width 2 -M 50 -N 50 -periodic -grid3d
110 output_file: output/ex7_1.out
111
112 TEST*/
113