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