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: }