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