Actual source code: ex7.c
1: static char help[] = "Tests DMLocalToLocalxxx() for DMDA.\n\n";
3: #include <petscdmda.h>
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;
20: PetscFunctionBeginUser;
21: PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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));
28: periodic = (DMBoundaryType)pt;
30: PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_type", &st, NULL));
32: stencil_type = (DMDAStencilType)st;
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));
46: PetscCall(DMCreateGlobalVector(da, &global));
47: PetscCall(DMCreateLocalVector(da, &local));
48: PetscCall(VecDuplicate(local, &local_copy));
50: /* zero out vectors so that ghostpoints are zero */
51: value = 0;
52: PetscCall(VecSet(local, value));
53: PetscCall(VecSet(local_copy, value));
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));
63: PetscCall(DMGlobalToLocalBegin(da, global, INSERT_VALUES, local));
64: PetscCall(DMGlobalToLocalEnd(da, global, INSERT_VALUES, local));
66: PetscCall(DMLocalToLocalBegin(da, local, INSERT_VALUES, local_copy));
67: PetscCall(DMLocalToLocalEnd(da, local, INSERT_VALUES, local_copy));
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: }
81: PetscCall(VecAXPY(local_copy, -1.0, local));
82: PetscCall(VecNorm(local_copy, NORM_MAX, &work));
83: PetscCall(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));
86: PetscCall(VecDestroy(&local_copy));
87: PetscCall(VecDestroy(&local));
88: PetscCall(VecDestroy(&global));
89: PetscCall(DMDestroy(&da));
90: PetscCall(PetscFinalize());
91: return 0;
92: }
94: /*TEST
96: test:
97: nsize: 8
98: args: -dof 3 -stencil_width 2 -M 50 -N 50 -periodic
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
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
112: TEST*/