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*/