1 static char help[] = "Solves the 1-dimensional wave equation.\n\n";
2
3 #include <petscdm.h>
4 #include <petscdmda.h>
5 #include <petscdraw.h>
6
main(int argc,char ** argv)7 int main(int argc, char **argv)
8 {
9 PetscMPIInt rank, size;
10 PetscInt M = 60, time_steps = 100, localsize, j, i, mybase, myend, *localnodes = NULL;
11 int width, xbase;
12 DM da;
13 PetscViewer viewer, viewer_private;
14 PetscDraw draw;
15 Vec local, global;
16 PetscScalar *localptr, *globalptr;
17 PetscReal a, h, k;
18 PetscBool flg = PETSC_FALSE;
19
20 PetscFunctionBeginUser;
21 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
22 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
24
25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL));
27 /*
28 Test putting two nodes on each processor, exact last processor gets the rest
29 */
30 PetscCall(PetscOptionsGetBool(NULL, NULL, "-distribute", &flg, NULL));
31 if (flg) {
32 PetscCall(PetscMalloc1(size, &localnodes));
33 for (i = 0; i < size - 1; i++) localnodes[i] = 2;
34 localnodes[size - 1] = M - 2 * (size - 1);
35 }
36
37 /* Set up the array */
38 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, M, 1, 1, localnodes, &da));
39 PetscCall(DMSetFromOptions(da));
40 PetscCall(DMSetUp(da));
41 PetscCall(PetscFree(localnodes));
42 PetscCall(DMCreateGlobalVector(da, &global));
43 PetscCall(DMCreateLocalVector(da, &local));
44
45 /* Set up display to show combined wave graph */
46 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "Entire Solution", 20, 480, 800, 200, &viewer));
47 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
48 PetscCall(PetscDrawSetDoubleBuffer(draw));
49
50 /* determine starting point of each processor */
51 PetscCall(VecGetOwnershipRange(global, &mybase, &myend));
52
53 /* set up display to show my portion of the wave */
54 xbase = (int)((mybase) * ((800.0 - 4.0 * size) / M) + 4.0 * rank);
55 width = (int)((myend - mybase) * 800. / M);
56 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "Local Portion of Solution", xbase, 200, width, 200, &viewer_private));
57 PetscCall(PetscViewerDrawGetDraw(viewer_private, 0, &draw));
58 PetscCall(PetscDrawSetDoubleBuffer(draw));
59
60 /* Initialize the array */
61 PetscCall(VecGetLocalSize(local, &localsize));
62 PetscCall(VecGetArray(global, &globalptr));
63
64 for (i = 1; i < localsize - 1; i++) {
65 j = (i - 1) + mybase;
66 globalptr[i - 1] = PetscSinReal((PETSC_PI * j * 6) / ((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI * j * 2) / ((PetscReal)M))) * 2;
67 }
68
69 PetscCall(VecRestoreArray(global, &globalptr));
70
71 /* Assign Parameters */
72 a = 1.0;
73 h = 1.0 / M;
74 k = h;
75
76 for (j = 0; j < time_steps; j++) {
77 /* Global to Local */
78 PetscCall(DMGlobalToLocalBegin(da, global, INSERT_VALUES, local));
79 PetscCall(DMGlobalToLocalEnd(da, global, INSERT_VALUES, local));
80
81 /*Extract local array */
82 PetscCall(VecGetArray(local, &localptr));
83 PetscCall(VecGetArray(global, &globalptr));
84
85 /* Update Locally - Make array of new values */
86 /* Note: I don't do anything for the first and last entry */
87 for (i = 1; i < localsize - 1; i++) globalptr[i - 1] = .5 * (localptr[i + 1] + localptr[i - 1]) - (k / (2.0 * a * h)) * (localptr[i + 1] - localptr[i - 1]);
88 PetscCall(VecRestoreArray(global, &globalptr));
89 PetscCall(VecRestoreArray(local, &localptr));
90
91 /* View my part of Wave */
92 PetscCall(VecView(global, viewer_private));
93
94 /* View global Wave */
95 PetscCall(VecView(global, viewer));
96 }
97
98 PetscCall(DMDestroy(&da));
99 PetscCall(PetscViewerDestroy(&viewer));
100 PetscCall(PetscViewerDestroy(&viewer_private));
101 PetscCall(VecDestroy(&local));
102 PetscCall(VecDestroy(&global));
103
104 PetscCall(PetscFinalize());
105 return 0;
106 }
107
108 /*TEST
109
110 test:
111 nsize: 3
112 args: -time 50 -nox
113 requires: x
114 output_file: output/empty.out
115
116 TEST*/
117