1 2 static char help[] = "Solves the 1-dimensional wave equation.\n\n"; 3 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 #include <petscdraw.h> 7 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; 19 20 PetscFunctionBeginUser; 21 PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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 115 TEST*/ 116