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 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, (char *)0, 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