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