1c4762a1bSJed Brown /*
2c4762a1bSJed Brown Simple example to show how PETSc programs can be run from MATLAB.
3c4762a1bSJed Brown See ex12.m.
4c4762a1bSJed Brown */
5c4762a1bSJed Brown
6c4762a1bSJed Brown static char help[] = "Solves the one dimensional heat equation.\n\n";
7c4762a1bSJed Brown
8c4762a1bSJed Brown #include <petscdm.h>
9c4762a1bSJed Brown #include <petscdmda.h>
10c4762a1bSJed Brown
main(int argc,char ** argv)11d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
12d71ae5a4SJacob Faibussowitsch {
13c4762a1bSJed Brown PetscMPIInt rank, size;
14c4762a1bSJed Brown PetscInt M = 14, time_steps = 20, w = 1, s = 1, localsize, j, i, mybase, myend, globalsize;
15c4762a1bSJed Brown DM da;
16c4762a1bSJed Brown Vec global, local;
17c4762a1bSJed Brown PetscScalar *globalptr, *localptr;
18c4762a1bSJed Brown PetscReal h, k;
19c4762a1bSJed Brown
20327415f7SBarry Smith PetscFunctionBeginUser;
21*c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &argv, NULL, help));
229566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
239566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL));
24c4762a1bSJed Brown
25c4762a1bSJed Brown /* Set up the array */
269566063dSJacob Faibussowitsch PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, M, w, s, NULL, &da));
279566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(da));
289566063dSJacob Faibussowitsch PetscCall(DMSetUp(da));
299566063dSJacob Faibussowitsch PetscCall(DMCreateGlobalVector(da, &global));
309566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
319566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
32c4762a1bSJed Brown
33c4762a1bSJed Brown /* Make copy of local array for doing updates */
349566063dSJacob Faibussowitsch PetscCall(DMCreateLocalVector(da, &local));
35c4762a1bSJed Brown
36c4762a1bSJed Brown /* determine starting point of each processor */
379566063dSJacob Faibussowitsch PetscCall(VecGetOwnershipRange(global, &mybase, &myend));
38c4762a1bSJed Brown
39c4762a1bSJed Brown /* Initialize the Array */
409566063dSJacob Faibussowitsch PetscCall(VecGetLocalSize(global, &globalsize));
419566063dSJacob Faibussowitsch PetscCall(VecGetArray(global, &globalptr));
42c4762a1bSJed Brown
43c4762a1bSJed Brown for (i = 0; i < globalsize; i++) {
44c4762a1bSJed Brown j = i + mybase;
45c4762a1bSJed Brown
46c4762a1bSJed Brown globalptr[i] = PetscSinReal((PETSC_PI * j * 6) / ((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI * j * 2) / ((PetscReal)M))) * 4 + 4;
47c4762a1bSJed Brown }
48c4762a1bSJed Brown
499566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(global, &localptr));
50c4762a1bSJed Brown
51c4762a1bSJed Brown /* Assign Parameters */
52c4762a1bSJed Brown h = 1.0 / M;
53c4762a1bSJed Brown k = h * h / 2.2;
549566063dSJacob Faibussowitsch PetscCall(VecGetLocalSize(local, &localsize));
55c4762a1bSJed Brown
56c4762a1bSJed Brown for (j = 0; j < time_steps; j++) {
57c4762a1bSJed Brown /* Global to Local */
589566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalBegin(da, global, INSERT_VALUES, local));
599566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalEnd(da, global, INSERT_VALUES, local));
60c4762a1bSJed Brown
61c4762a1bSJed Brown /*Extract local array */
629566063dSJacob Faibussowitsch PetscCall(VecGetArray(local, &localptr));
639566063dSJacob Faibussowitsch PetscCall(VecGetArray(global, &globalptr));
64c4762a1bSJed Brown
65c4762a1bSJed Brown /* Update Locally - Make array of new values */
66c4762a1bSJed Brown /* Note: I don't do anything for the first and last entry */
67ad540459SPierre Jolivet 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]);
68c4762a1bSJed Brown
699566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(global, &globalptr));
709566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(local, &localptr));
71c4762a1bSJed Brown
72c4762a1bSJed Brown /* View Wave */
73c4762a1bSJed Brown /* Set Up Display to Show Heat Graph */
74c4762a1bSJed Brown #if defined(PETSC_USE_SOCKET_VIEWER)
759566063dSJacob Faibussowitsch PetscCall(VecView(global, PETSC_VIEWER_SOCKET_WORLD));
76c4762a1bSJed Brown #endif
77c4762a1bSJed Brown }
78c4762a1bSJed Brown
799566063dSJacob Faibussowitsch PetscCall(VecDestroy(&local));
809566063dSJacob Faibussowitsch PetscCall(VecDestroy(&global));
819566063dSJacob Faibussowitsch PetscCall(DMDestroy(&da));
829566063dSJacob Faibussowitsch PetscCall(PetscFinalize());
83b122ec5aSJacob Faibussowitsch return 0;
84c4762a1bSJed Brown }
85