1 2 static char help[] = "Solves the 1-dimensional wave equation.\n\n"; 3 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 #include <petscdraw.h> 7 8 int main(int argc,char **argv) 9 { 10 PetscMPIInt rank,size; 11 PetscErrorCode ierr; 12 PetscInt M = 60,time_steps = 100, localsize,j,i,mybase,myend,width,xbase,*localnodes = NULL; 13 DM da; 14 PetscViewer viewer,viewer_private; 15 PetscDraw draw; 16 Vec local,global; 17 PetscScalar *localptr,*globalptr; 18 PetscReal a,h,k; 19 PetscBool flg = PETSC_FALSE; 20 21 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 22 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 23 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 24 25 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 26 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL)); 27 /* 28 Test putting two nodes on each processor, exact last processor gets the rest 29 */ 30 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-distribute",&flg,NULL)); 31 if (flg) { 32 CHKERRQ(PetscMalloc1(size,&localnodes)); 33 for (i=0; i<size-1; i++) localnodes[i] = 2; 34 localnodes[size-1] = M - 2*(size-1); 35 } 36 37 /* Set up the array */ 38 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,1,localnodes,&da)); 39 CHKERRQ(DMSetFromOptions(da)); 40 CHKERRQ(DMSetUp(da)); 41 CHKERRQ(PetscFree(localnodes)); 42 CHKERRQ(DMCreateGlobalVector(da,&global)); 43 CHKERRQ(DMCreateLocalVector(da,&local)); 44 45 /* Set up display to show combined wave graph */ 46 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"Entire Solution",20,480,800,200,&viewer)); 47 CHKERRQ(PetscViewerDrawGetDraw(viewer,0,&draw)); 48 CHKERRQ(PetscDrawSetDoubleBuffer(draw)); 49 50 /* determine starting point of each processor */ 51 CHKERRQ(VecGetOwnershipRange(global,&mybase,&myend)); 52 53 /* set up display to show my portion of the wave */ 54 xbase = (int)((mybase)*((800.0 - 4.0*size)/M) + 4.0*rank); 55 width = (int)((myend-mybase)*800./M); 56 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Local Portion of Solution",xbase,200,width,200,&viewer_private)); 57 CHKERRQ(PetscViewerDrawGetDraw(viewer_private,0,&draw)); 58 CHKERRQ(PetscDrawSetDoubleBuffer(draw)); 59 60 /* Initialize the array */ 61 CHKERRQ(VecGetLocalSize(local,&localsize)); 62 CHKERRQ(VecGetArray(global,&globalptr)); 63 64 for (i=1; i<localsize-1; i++) { 65 j = (i-1)+mybase; 66 globalptr[i-1] = PetscSinReal((PETSC_PI*j*6)/((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI*j*2)/((PetscReal)M))) * 2; 67 } 68 69 CHKERRQ(VecRestoreArray(global,&globalptr)); 70 71 /* Assign Parameters */ 72 a= 1.0; 73 h= 1.0/M; 74 k= h; 75 76 for (j=0; j<time_steps; j++) { 77 78 /* Global to Local */ 79 CHKERRQ(DMGlobalToLocalBegin(da,global,INSERT_VALUES,local)); 80 CHKERRQ(DMGlobalToLocalEnd(da,global,INSERT_VALUES,local)); 81 82 /*Extract local array */ 83 CHKERRQ(VecGetArray(local,&localptr)); 84 CHKERRQ(VecGetArray(global,&globalptr)); 85 86 /* Update Locally - Make array of new values */ 87 /* Note: I don't do anything for the first and last entry */ 88 for (i=1; i< localsize-1; i++) { 89 globalptr[i-1] = .5*(localptr[i+1]+localptr[i-1]) - (k / (2.0*a*h)) * (localptr[i+1] - localptr[i-1]); 90 } 91 CHKERRQ(VecRestoreArray(global,&globalptr)); 92 CHKERRQ(VecRestoreArray(local,&localptr)); 93 94 /* View my part of Wave */ 95 CHKERRQ(VecView(global,viewer_private)); 96 97 /* View global Wave */ 98 CHKERRQ(VecView(global,viewer)); 99 } 100 101 CHKERRQ(DMDestroy(&da)); 102 CHKERRQ(PetscViewerDestroy(&viewer)); 103 CHKERRQ(PetscViewerDestroy(&viewer_private)); 104 CHKERRQ(VecDestroy(&local)); 105 CHKERRQ(VecDestroy(&global)); 106 107 ierr = PetscFinalize(); 108 return ierr; 109 } 110 111 /*TEST 112 113 test: 114 nsize: 3 115 args: -time 50 -nox 116 requires: x 117 118 TEST*/ 119