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 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 23 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 24 25 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 26 ierr = PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);CHKERRQ(ierr); 27 /* 28 Test putting two nodes on each processor, exact last processor gets the rest 29 */ 30 ierr = PetscOptionsGetBool(NULL,NULL,"-distribute",&flg,NULL);CHKERRQ(ierr); 31 if (flg) { 32 ierr = PetscMalloc1(size,&localnodes);CHKERRQ(ierr); 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 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,1,localnodes,&da);CHKERRQ(ierr); 39 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 40 ierr = DMSetUp(da);CHKERRQ(ierr); 41 ierr = PetscFree(localnodes);CHKERRQ(ierr); 42 ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr); 43 ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr); 44 45 /* Set up display to show combined wave graph */ 46 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"Entire Solution",20,480,800,200,&viewer);CHKERRQ(ierr); 47 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 48 ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr); 49 50 /* determine starting point of each processor */ 51 ierr = VecGetOwnershipRange(global,&mybase,&myend);CHKERRQ(ierr); 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 ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Local Portion of Solution",xbase,200,width,200,&viewer_private);CHKERRQ(ierr); 57 ierr = PetscViewerDrawGetDraw(viewer_private,0,&draw);CHKERRQ(ierr); 58 ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr); 59 60 61 62 /* Initialize the array */ 63 ierr = VecGetLocalSize(local,&localsize);CHKERRQ(ierr); 64 ierr = VecGetArray(global,&globalptr);CHKERRQ(ierr); 65 66 for (i=1; i<localsize-1; i++) { 67 j = (i-1)+mybase; 68 globalptr[i-1] = PetscSinReal((PETSC_PI*j*6)/((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI*j*2)/((PetscReal)M))) * 2; 69 } 70 71 ierr = VecRestoreArray(global,&globalptr);CHKERRQ(ierr); 72 73 /* Assign Parameters */ 74 a= 1.0; 75 h= 1.0/M; 76 k= h; 77 78 for (j=0; j<time_steps; j++) { 79 80 /* Global to Local */ 81 ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr); 82 ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr); 83 84 /*Extract local array */ 85 ierr = VecGetArray(local,&localptr);CHKERRQ(ierr); 86 ierr = VecGetArray(global,&globalptr);CHKERRQ(ierr); 87 88 /* Update Locally - Make array of new values */ 89 /* Note: I don't do anything for the first and last entry */ 90 for (i=1; i< localsize-1; i++) { 91 globalptr[i-1] = .5*(localptr[i+1]+localptr[i-1]) - (k / (2.0*a*h)) * (localptr[i+1] - localptr[i-1]); 92 } 93 ierr = VecRestoreArray(global,&globalptr);CHKERRQ(ierr); 94 ierr = VecRestoreArray(local,&localptr);CHKERRQ(ierr); 95 96 /* View my part of Wave */ 97 ierr = VecView(global,viewer_private);CHKERRQ(ierr); 98 99 /* View global Wave */ 100 ierr = VecView(global,viewer);CHKERRQ(ierr); 101 } 102 103 ierr = DMDestroy(&da);CHKERRQ(ierr); 104 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 105 ierr = PetscViewerDestroy(&viewer_private);CHKERRQ(ierr); 106 ierr = VecDestroy(&local);CHKERRQ(ierr); 107 ierr = VecDestroy(&global);CHKERRQ(ierr); 108 109 ierr = PetscFinalize(); 110 return ierr; 111 } 112 113 /*TEST 114 115 test: 116 nsize: 3 117 args: -time 50 -nox 118 requires: x 119 120 TEST*/ 121