1*c4762a1bSJed Brown 2*c4762a1bSJed Brown /* 3*c4762a1bSJed Brown Simple example to show how PETSc programs can be run from MATLAB. 4*c4762a1bSJed Brown See ex12.m. 5*c4762a1bSJed Brown */ 6*c4762a1bSJed Brown 7*c4762a1bSJed Brown static char help[] = "Solves the one dimensional heat equation.\n\n"; 8*c4762a1bSJed Brown 9*c4762a1bSJed Brown #include <petscdm.h> 10*c4762a1bSJed Brown #include <petscdmda.h> 11*c4762a1bSJed Brown 12*c4762a1bSJed Brown int main(int argc,char **argv) 13*c4762a1bSJed Brown { 14*c4762a1bSJed Brown PetscMPIInt rank,size; 15*c4762a1bSJed Brown PetscInt M = 14,time_steps = 20,w=1,s=1,localsize,j,i,mybase,myend,globalsize; 16*c4762a1bSJed Brown PetscErrorCode ierr; 17*c4762a1bSJed Brown DM da; 18*c4762a1bSJed Brown Vec global,local; 19*c4762a1bSJed Brown PetscScalar *globalptr,*localptr; 20*c4762a1bSJed Brown PetscReal h,k; 21*c4762a1bSJed Brown 22*c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 23*c4762a1bSJed Brown ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 24*c4762a1bSJed Brown ierr = PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);CHKERRQ(ierr); 25*c4762a1bSJed Brown 26*c4762a1bSJed Brown /* Set up the array */ 27*c4762a1bSJed Brown ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,M,w,s,NULL,&da);CHKERRQ(ierr); 28*c4762a1bSJed Brown ierr = DMSetFromOptions(da);CHKERRQ(ierr); 29*c4762a1bSJed Brown ierr = DMSetUp(da);CHKERRQ(ierr); 30*c4762a1bSJed Brown ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr); 31*c4762a1bSJed Brown ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 32*c4762a1bSJed Brown ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 33*c4762a1bSJed Brown 34*c4762a1bSJed Brown /* Make copy of local array for doing updates */ 35*c4762a1bSJed Brown ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr); 36*c4762a1bSJed Brown 37*c4762a1bSJed Brown 38*c4762a1bSJed Brown /* determine starting point of each processor */ 39*c4762a1bSJed Brown ierr = VecGetOwnershipRange(global,&mybase,&myend);CHKERRQ(ierr); 40*c4762a1bSJed Brown 41*c4762a1bSJed Brown /* Initialize the Array */ 42*c4762a1bSJed Brown ierr = VecGetLocalSize (global,&globalsize);CHKERRQ(ierr); 43*c4762a1bSJed Brown ierr = VecGetArray (global,&globalptr);CHKERRQ(ierr); 44*c4762a1bSJed Brown 45*c4762a1bSJed Brown 46*c4762a1bSJed Brown for (i=0; i<globalsize; i++) { 47*c4762a1bSJed Brown j = i + mybase; 48*c4762a1bSJed Brown 49*c4762a1bSJed Brown globalptr[i] = PetscSinReal((PETSC_PI*j*6)/((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI*j*2)/((PetscReal)M))) * 4+4; 50*c4762a1bSJed Brown } 51*c4762a1bSJed Brown 52*c4762a1bSJed Brown ierr = VecRestoreArray(global,&localptr);CHKERRQ(ierr); 53*c4762a1bSJed Brown 54*c4762a1bSJed Brown /* Assign Parameters */ 55*c4762a1bSJed Brown h= 1.0/M; 56*c4762a1bSJed Brown k= h*h/2.2; 57*c4762a1bSJed Brown ierr = VecGetLocalSize(local,&localsize);CHKERRQ(ierr); 58*c4762a1bSJed Brown 59*c4762a1bSJed Brown for (j=0; j<time_steps; j++) { 60*c4762a1bSJed Brown 61*c4762a1bSJed Brown /* Global to Local */ 62*c4762a1bSJed Brown ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr); 63*c4762a1bSJed Brown ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr); 64*c4762a1bSJed Brown 65*c4762a1bSJed Brown /*Extract local array */ 66*c4762a1bSJed Brown ierr = VecGetArray(local,&localptr);CHKERRQ(ierr); 67*c4762a1bSJed Brown ierr = VecGetArray (global,&globalptr);CHKERRQ(ierr); 68*c4762a1bSJed Brown 69*c4762a1bSJed Brown /* Update Locally - Make array of new values */ 70*c4762a1bSJed Brown /* Note: I don't do anything for the first and last entry */ 71*c4762a1bSJed Brown for (i=1; i< localsize-1; i++) { 72*c4762a1bSJed Brown globalptr[i-1] = localptr[i] + (k/(h*h)) * (localptr[i+1]-2.0*localptr[i]+localptr[i-1]); 73*c4762a1bSJed Brown } 74*c4762a1bSJed Brown 75*c4762a1bSJed Brown ierr = VecRestoreArray (global,&globalptr);CHKERRQ(ierr); 76*c4762a1bSJed Brown ierr = VecRestoreArray(local,&localptr);CHKERRQ(ierr); 77*c4762a1bSJed Brown 78*c4762a1bSJed Brown /* View Wave */ 79*c4762a1bSJed Brown /* Set Up Display to Show Heat Graph */ 80*c4762a1bSJed Brown #if defined(PETSC_USE_SOCKET_VIEWER) 81*c4762a1bSJed Brown ierr = VecView(global,PETSC_VIEWER_SOCKET_WORLD);CHKERRQ(ierr); 82*c4762a1bSJed Brown #endif 83*c4762a1bSJed Brown } 84*c4762a1bSJed Brown 85*c4762a1bSJed Brown ierr = VecDestroy(&local);CHKERRQ(ierr); 86*c4762a1bSJed Brown ierr = VecDestroy(&global);CHKERRQ(ierr); 87*c4762a1bSJed Brown ierr = DMDestroy(&da);CHKERRQ(ierr); 88*c4762a1bSJed Brown ierr = PetscFinalize(); 89*c4762a1bSJed Brown return ierr; 90*c4762a1bSJed Brown } 91*c4762a1bSJed Brown 92