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