1 2 /* 3 Simple example to show how PETSc programs can be run from MATLAB. 4 See ex12.m. 5 */ 6 7 static char help[] = "Solves the one dimensional heat equation.\n\n"; 8 9 #include <petscdm.h> 10 #include <petscdmda.h> 11 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; 20 21 PetscFunctionBeginUser; 22 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL)); 25 26 /* Set up the array */ 27 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, M, w, s, NULL, &da)); 28 PetscCall(DMSetFromOptions(da)); 29 PetscCall(DMSetUp(da)); 30 PetscCall(DMCreateGlobalVector(da, &global)); 31 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 32 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 33 34 /* Make copy of local array for doing updates */ 35 PetscCall(DMCreateLocalVector(da, &local)); 36 37 /* determine starting point of each processor */ 38 PetscCall(VecGetOwnershipRange(global, &mybase, &myend)); 39 40 /* Initialize the Array */ 41 PetscCall(VecGetLocalSize(global, &globalsize)); 42 PetscCall(VecGetArray(global, &globalptr)); 43 44 for (i = 0; i < globalsize; i++) { 45 j = i + mybase; 46 47 globalptr[i] = PetscSinReal((PETSC_PI * j * 6) / ((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI * j * 2) / ((PetscReal)M))) * 4 + 4; 48 } 49 50 PetscCall(VecRestoreArray(global, &localptr)); 51 52 /* Assign Parameters */ 53 h = 1.0 / M; 54 k = h * h / 2.2; 55 PetscCall(VecGetLocalSize(local, &localsize)); 56 57 for (j = 0; j < time_steps; j++) { 58 /* Global to Local */ 59 PetscCall(DMGlobalToLocalBegin(da, global, INSERT_VALUES, local)); 60 PetscCall(DMGlobalToLocalEnd(da, global, INSERT_VALUES, local)); 61 62 /*Extract local array */ 63 PetscCall(VecGetArray(local, &localptr)); 64 PetscCall(VecGetArray(global, &globalptr)); 65 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]); 69 70 PetscCall(VecRestoreArray(global, &globalptr)); 71 PetscCall(VecRestoreArray(local, &localptr)); 72 73 /* View Wave */ 74 /* Set Up Display to Show Heat Graph */ 75 #if defined(PETSC_USE_SOCKET_VIEWER) 76 PetscCall(VecView(global, PETSC_VIEWER_SOCKET_WORLD)); 77 #endif 78 } 79 80 PetscCall(VecDestroy(&local)); 81 PetscCall(VecDestroy(&global)); 82 PetscCall(DMDestroy(&da)); 83 PetscCall(PetscFinalize()); 84 return 0; 85 } 86