Actual source code: ex3.c
2: static char help[] = "Solves the 1-dimensional wave equation.\n\n";
4: #include <petscdm.h>
5: #include <petscdmda.h>
6: #include <petscdraw.h>
8: int main(int argc, char **argv)
9: {
10: PetscMPIInt rank, size;
11: PetscInt M = 60, time_steps = 100, localsize, j, i, mybase, myend, width, xbase, *localnodes = NULL;
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;
21: PetscInitialize(&argc, &argv, (char *)0, help);
22: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
23: MPI_Comm_size(PETSC_COMM_WORLD, &size);
25: PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
26: PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL);
27: /*
28: Test putting two nodes on each processor, exact last processor gets the rest
29: */
30: PetscOptionsGetBool(NULL, NULL, "-distribute", &flg, NULL);
31: if (flg) {
32: PetscMalloc1(size, &localnodes);
33: for (i = 0; i < size - 1; i++) localnodes[i] = 2;
34: localnodes[size - 1] = M - 2 * (size - 1);
35: }
37: /* Set up the array */
38: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, M, 1, 1, localnodes, &da);
39: DMSetFromOptions(da);
40: DMSetUp(da);
41: PetscFree(localnodes);
42: DMCreateGlobalVector(da, &global);
43: DMCreateLocalVector(da, &local);
45: /* Set up display to show combined wave graph */
46: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "Entire Solution", 20, 480, 800, 200, &viewer);
47: PetscViewerDrawGetDraw(viewer, 0, &draw);
48: PetscDrawSetDoubleBuffer(draw);
50: /* determine starting point of each processor */
51: VecGetOwnershipRange(global, &mybase, &myend);
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: PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "Local Portion of Solution", xbase, 200, width, 200, &viewer_private);
57: PetscViewerDrawGetDraw(viewer_private, 0, &draw);
58: PetscDrawSetDoubleBuffer(draw);
60: /* Initialize the array */
61: VecGetLocalSize(local, &localsize);
62: VecGetArray(global, &globalptr);
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: }
69: VecRestoreArray(global, &globalptr);
71: /* Assign Parameters */
72: a = 1.0;
73: h = 1.0 / M;
74: k = h;
76: for (j = 0; j < time_steps; j++) {
77: /* Global to Local */
78: DMGlobalToLocalBegin(da, global, INSERT_VALUES, local);
79: DMGlobalToLocalEnd(da, global, INSERT_VALUES, local);
81: /*Extract local array */
82: VecGetArray(local, &localptr);
83: VecGetArray(global, &globalptr);
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: VecRestoreArray(global, &globalptr);
89: VecRestoreArray(local, &localptr);
91: /* View my part of Wave */
92: VecView(global, viewer_private);
94: /* View global Wave */
95: VecView(global, viewer);
96: }
98: DMDestroy(&da);
99: PetscViewerDestroy(&viewer);
100: PetscViewerDestroy(&viewer_private);
101: VecDestroy(&local);
102: VecDestroy(&global);
104: PetscFinalize();
105: return 0;
106: }
108: /*TEST
110: test:
111: nsize: 3
112: args: -time 50 -nox
113: requires: x
115: TEST*/