1 /*
2 Simple example to show how PETSc programs can be run from MATLAB.
3 See ex12.m.
4 */
5
6 static char help[] = "Solves the one dimensional heat equation.\n\n";
7
8 #include <petscdm.h>
9 #include <petscdmda.h>
10
main(int argc,char ** argv)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;
19
20 PetscFunctionBeginUser;
21 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
22 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL));
24
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));
32
33 /* Make copy of local array for doing updates */
34 PetscCall(DMCreateLocalVector(da, &local));
35
36 /* determine starting point of each processor */
37 PetscCall(VecGetOwnershipRange(global, &mybase, &myend));
38
39 /* Initialize the Array */
40 PetscCall(VecGetLocalSize(global, &globalsize));
41 PetscCall(VecGetArray(global, &globalptr));
42
43 for (i = 0; i < globalsize; i++) {
44 j = i + mybase;
45
46 globalptr[i] = PetscSinReal((PETSC_PI * j * 6) / ((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI * j * 2) / ((PetscReal)M))) * 4 + 4;
47 }
48
49 PetscCall(VecRestoreArray(global, &localptr));
50
51 /* Assign Parameters */
52 h = 1.0 / M;
53 k = h * h / 2.2;
54 PetscCall(VecGetLocalSize(local, &localsize));
55
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));
60
61 /*Extract local array */
62 PetscCall(VecGetArray(local, &localptr));
63 PetscCall(VecGetArray(global, &globalptr));
64
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]);
68
69 PetscCall(VecRestoreArray(global, &globalptr));
70 PetscCall(VecRestoreArray(local, &localptr));
71
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 }
78
79 PetscCall(VecDestroy(&local));
80 PetscCall(VecDestroy(&global));
81 PetscCall(DMDestroy(&da));
82 PetscCall(PetscFinalize());
83 return 0;
84 }
85