Actual source code: ex12.c

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

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

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

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

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

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

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

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

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

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

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

 49:   PetscCall(VecRestoreArray(global, &localptr));

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

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

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

 65:     /* Update Locally - Make array of new values */
 66:     /* Note: I don't do anything for the first and last entry */
 67:     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]);

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

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

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