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 PetscFunctionBeginUser; 21 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 22 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 24 25 PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 26 PetscCall(PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL)); 27 /* 28 Test putting two nodes on each processor, exact last processor gets the rest 29 */ 30 PetscCall(PetscOptionsGetBool(NULL,NULL,"-distribute",&flg,NULL)); 31 if (flg) { 32 PetscCall(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 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,1,localnodes,&da)); 39 PetscCall(DMSetFromOptions(da)); 40 PetscCall(DMSetUp(da)); 41 PetscCall(PetscFree(localnodes)); 42 PetscCall(DMCreateGlobalVector(da,&global)); 43 PetscCall(DMCreateLocalVector(da,&local)); 44 45 /* Set up display to show combined wave graph */ 46 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"Entire Solution",20,480,800,200,&viewer)); 47 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 48 PetscCall(PetscDrawSetDoubleBuffer(draw)); 49 50 /* determine starting point of each processor */ 51 PetscCall(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 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Local Portion of Solution",xbase,200,width,200,&viewer_private)); 57 PetscCall(PetscViewerDrawGetDraw(viewer_private,0,&draw)); 58 PetscCall(PetscDrawSetDoubleBuffer(draw)); 59 60 /* Initialize the array */ 61 PetscCall(VecGetLocalSize(local,&localsize)); 62 PetscCall(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 PetscCall(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 PetscCall(DMGlobalToLocalBegin(da,global,INSERT_VALUES,local)); 80 PetscCall(DMGlobalToLocalEnd(da,global,INSERT_VALUES,local)); 81 82 /*Extract local array */ 83 PetscCall(VecGetArray(local,&localptr)); 84 PetscCall(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 PetscCall(VecRestoreArray(global,&globalptr)); 92 PetscCall(VecRestoreArray(local,&localptr)); 93 94 /* View my part of Wave */ 95 PetscCall(VecView(global,viewer_private)); 96 97 /* View global Wave */ 98 PetscCall(VecView(global,viewer)); 99 } 100 101 PetscCall(DMDestroy(&da)); 102 PetscCall(PetscViewerDestroy(&viewer)); 103 PetscCall(PetscViewerDestroy(&viewer_private)); 104 PetscCall(VecDestroy(&local)); 105 PetscCall(VecDestroy(&global)); 106 107 PetscCall(PetscFinalize()); 108 return 0; 109 } 110 111 /*TEST 112 113 test: 114 nsize: 3 115 args: -time 50 -nox 116 requires: x 117 118 TEST*/ 119