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 /* Initialize the array */ 61 ierr = VecGetLocalSize(local,&localsize);CHKERRQ(ierr); 62 ierr = VecGetArray(global,&globalptr);CHKERRQ(ierr); 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 ierr = VecRestoreArray(global,&globalptr);CHKERRQ(ierr); 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 ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr); 80 ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr); 81 82 /*Extract local array */ 83 ierr = VecGetArray(local,&localptr);CHKERRQ(ierr); 84 ierr = VecGetArray(global,&globalptr);CHKERRQ(ierr); 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 ierr = VecRestoreArray(global,&globalptr);CHKERRQ(ierr); 92 ierr = VecRestoreArray(local,&localptr);CHKERRQ(ierr); 93 94 /* View my part of Wave */ 95 ierr = VecView(global,viewer_private);CHKERRQ(ierr); 96 97 /* View global Wave */ 98 ierr = VecView(global,viewer);CHKERRQ(ierr); 99 } 100 101 ierr = DMDestroy(&da);CHKERRQ(ierr); 102 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 103 ierr = PetscViewerDestroy(&viewer_private);CHKERRQ(ierr); 104 ierr = VecDestroy(&local);CHKERRQ(ierr); 105 ierr = VecDestroy(&global);CHKERRQ(ierr); 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