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