Actual source code: ex12.c


  2: /*
  3:    Simple example to show how PETSc programs can be run from MATLAB.
  4:   See ex12.m.
  5: */

  7: static char help[] = "Solves the one dimensional heat equation.\n\n";

  9: #include <petscdm.h>
 10: #include <petscdmda.h>

 12: int main(int argc, char **argv)
 13: {
 14:   PetscMPIInt  rank, size;
 15:   PetscInt     M = 14, time_steps = 20, w = 1, s = 1, localsize, j, i, mybase, myend, globalsize;
 16:   DM           da;
 17:   Vec          global, local;
 18:   PetscScalar *globalptr, *localptr;
 19:   PetscReal    h, k;

 22:   PetscInitialize(&argc, &argv, (char *)0, help);
 23:   PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
 24:   PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL);

 26:   /* Set up the array */
 27:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, M, w, s, NULL, &da);
 28:   DMSetFromOptions(da);
 29:   DMSetUp(da);
 30:   DMCreateGlobalVector(da, &global);
 31:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 32:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 34:   /* Make copy of local array for doing updates */
 35:   DMCreateLocalVector(da, &local);

 37:   /* determine starting point of each processor */
 38:   VecGetOwnershipRange(global, &mybase, &myend);

 40:   /* Initialize the Array */
 41:   VecGetLocalSize(global, &globalsize);
 42:   VecGetArray(global, &globalptr);

 44:   for (i = 0; i < globalsize; i++) {
 45:     j = i + mybase;

 47:     globalptr[i] = PetscSinReal((PETSC_PI * j * 6) / ((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI * j * 2) / ((PetscReal)M))) * 4 + 4;
 48:   }

 50:   VecRestoreArray(global, &localptr);

 52:   /* Assign Parameters */
 53:   h = 1.0 / M;
 54:   k = h * h / 2.2;
 55:   VecGetLocalSize(local, &localsize);

 57:   for (j = 0; j < time_steps; j++) {
 58:     /* Global to Local */
 59:     DMGlobalToLocalBegin(da, global, INSERT_VALUES, local);
 60:     DMGlobalToLocalEnd(da, global, INSERT_VALUES, local);

 62:     /*Extract local array */
 63:     VecGetArray(local, &localptr);
 64:     VecGetArray(global, &globalptr);

 66:     /* Update Locally - Make array of new values */
 67:     /* Note: I don't do anything for the first and last entry */
 68:     for (i = 1; i < localsize - 1; i++) globalptr[i - 1] = localptr[i] + (k / (h * h)) * (localptr[i + 1] - 2.0 * localptr[i] + localptr[i - 1]);

 70:     VecRestoreArray(global, &globalptr);
 71:     VecRestoreArray(local, &localptr);

 73:     /* View Wave */
 74:     /* Set Up Display to Show Heat Graph */
 75: #if defined(PETSC_USE_SOCKET_VIEWER)
 76:     VecView(global, PETSC_VIEWER_SOCKET_WORLD);
 77: #endif
 78:   }

 80:   VecDestroy(&local);
 81:   VecDestroy(&global);
 82:   DMDestroy(&da);
 83:   PetscFinalize();
 84:   return 0;
 85: }