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