1*c4762a1bSJed Brown 2*c4762a1bSJed Brown static char help[] = "\n\n"; 3*c4762a1bSJed Brown 4*c4762a1bSJed Brown /* 5*c4762a1bSJed Brown Demonstrates using DM_BOUNDARY_GHOSTED how to handle a rotated boundary conditions where one edge 6*c4762a1bSJed Brown is connected to its immediate neighbor 7*c4762a1bSJed Brown 8*c4762a1bSJed Brown Consider the domain (with natural numbering) 9*c4762a1bSJed Brown 10*c4762a1bSJed Brown 6 7 8 11*c4762a1bSJed Brown 3 4 5 12*c4762a1bSJed Brown 0 1 2 13*c4762a1bSJed Brown 14*c4762a1bSJed Brown The ghost points along the bottom (directly below the three columns above) should be 0 3 and 6 15*c4762a1bSJed Brown while the ghost points along the left side should be 0 1 2 16*c4762a1bSJed Brown 17*c4762a1bSJed Brown Note that the ghosted local vectors extend in both the x and y directions so, for example if we have a 18*c4762a1bSJed Brown single MPI process the ghosted vector has (in the original natural numbering) 19*c4762a1bSJed Brown 20*c4762a1bSJed Brown x x x x x 21*c4762a1bSJed Brown 2 6 7 8 x 22*c4762a1bSJed Brown 1 3 4 5 x 23*c4762a1bSJed Brown 0 0 1 2 x 24*c4762a1bSJed Brown x 0 3 6 x 25*c4762a1bSJed Brown 26*c4762a1bSJed Brown where x indicates a location that is not updated by the communication and should be used. 27*c4762a1bSJed Brown 28*c4762a1bSJed Brown For this to make sense the number of grid points in the x and y directions must be the same 29*c4762a1bSJed Brown 30*c4762a1bSJed Brown This ghost point mapping was suggested by: Wenbo Zhao <zhaowenbo.npic@gmail.com> 31*c4762a1bSJed Brown */ 32*c4762a1bSJed Brown 33*c4762a1bSJed Brown #include <petscdm.h> 34*c4762a1bSJed Brown #include <petscdmda.h> 35*c4762a1bSJed Brown 36*c4762a1bSJed Brown int main(int argc,char **argv) 37*c4762a1bSJed Brown { 38*c4762a1bSJed Brown PetscInt M = 6; 39*c4762a1bSJed Brown PetscErrorCode ierr; 40*c4762a1bSJed Brown DM da; 41*c4762a1bSJed Brown Vec local,global,natural; 42*c4762a1bSJed Brown PetscInt i,start,end,*ifrom,x,y,xm,ym; 43*c4762a1bSJed Brown PetscScalar *xnatural; 44*c4762a1bSJed Brown IS from,to; 45*c4762a1bSJed Brown AO ao; 46*c4762a1bSJed Brown VecScatter scatter1,scatter2; 47*c4762a1bSJed Brown PetscViewer subviewer; 48*c4762a1bSJed Brown 49*c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 50*c4762a1bSJed Brown 51*c4762a1bSJed Brown /* Create distributed array and get vectors */ 52*c4762a1bSJed Brown ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_GHOSTED,DM_BOUNDARY_GHOSTED,DMDA_STENCIL_STAR,M,M,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr); 53*c4762a1bSJed Brown ierr = DMSetUp(da);CHKERRQ(ierr); 54*c4762a1bSJed Brown ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr); 55*c4762a1bSJed Brown ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr); 56*c4762a1bSJed Brown 57*c4762a1bSJed Brown /* construct global to local scatter for the left side of the domain to the ghost on the bottom */ 58*c4762a1bSJed Brown ierr = DMDAGetCorners(da,&x,&y,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 59*c4762a1bSJed Brown if (!y) { /* only processes on the bottom of the domain fill up the ghost locations */ 60*c4762a1bSJed Brown ierr = ISCreateStride(PETSC_COMM_SELF,xm,1,1,&to);CHKERRQ(ierr); 61*c4762a1bSJed Brown } else { 62*c4762a1bSJed Brown ierr = ISCreateStride(PETSC_COMM_SELF,0,0,0,&to);CHKERRQ(ierr); 63*c4762a1bSJed Brown } 64*c4762a1bSJed Brown ierr = PetscMalloc1(xm,&ifrom);CHKERRQ(ierr); 65*c4762a1bSJed Brown for (i=x;i<x+xm;i++) { 66*c4762a1bSJed Brown ifrom[i-x] = M*i; 67*c4762a1bSJed Brown } 68*c4762a1bSJed Brown ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 69*c4762a1bSJed Brown ierr = AOApplicationToPetsc(ao,xm,ifrom);CHKERRQ(ierr); 70*c4762a1bSJed Brown if (!y) { 71*c4762a1bSJed Brown ierr = ISCreateGeneral(PETSC_COMM_WORLD,xm,ifrom,PETSC_OWN_POINTER,&from);CHKERRQ(ierr); 72*c4762a1bSJed Brown } else { 73*c4762a1bSJed Brown ierr = PetscFree(ifrom);CHKERRQ(ierr); 74*c4762a1bSJed Brown ierr = ISCreateGeneral(PETSC_COMM_WORLD,0,NULL,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 75*c4762a1bSJed Brown } 76*c4762a1bSJed Brown ierr = VecScatterCreate(global,from,local,to,&scatter1);CHKERRQ(ierr); 77*c4762a1bSJed Brown ierr = ISDestroy(&to);CHKERRQ(ierr); 78*c4762a1bSJed Brown ierr = ISDestroy(&from);CHKERRQ(ierr); 79*c4762a1bSJed Brown 80*c4762a1bSJed Brown /* construct global to local scatter for the bottom side of the domain to the ghost on the right */ 81*c4762a1bSJed Brown if (!x) { /* only processes on the left side of the domain fill up the ghost locations */ 82*c4762a1bSJed Brown ierr = ISCreateStride(PETSC_COMM_SELF,ym,xm+2,xm+2,&to);CHKERRQ(ierr); 83*c4762a1bSJed Brown } else { 84*c4762a1bSJed Brown ierr = ISCreateStride(PETSC_COMM_SELF,0,0,0,&to);CHKERRQ(ierr); 85*c4762a1bSJed Brown } 86*c4762a1bSJed Brown ierr = PetscMalloc1(ym,&ifrom);CHKERRQ(ierr); 87*c4762a1bSJed Brown for (i=y;i<y+ym;i++) { 88*c4762a1bSJed Brown ifrom[i-y] = i; 89*c4762a1bSJed Brown } 90*c4762a1bSJed Brown ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 91*c4762a1bSJed Brown ierr = AOApplicationToPetsc(ao,ym,ifrom);CHKERRQ(ierr); 92*c4762a1bSJed Brown if (!x) { 93*c4762a1bSJed Brown ierr = ISCreateGeneral(PETSC_COMM_WORLD,ym,ifrom,PETSC_OWN_POINTER,&from);CHKERRQ(ierr); 94*c4762a1bSJed Brown } else { 95*c4762a1bSJed Brown ierr = PetscFree(ifrom);CHKERRQ(ierr); 96*c4762a1bSJed Brown ierr = ISCreateGeneral(PETSC_COMM_WORLD,0,NULL,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 97*c4762a1bSJed Brown } 98*c4762a1bSJed Brown ierr = VecScatterCreate(global,from,local,to,&scatter2);CHKERRQ(ierr); 99*c4762a1bSJed Brown ierr = ISDestroy(&to);CHKERRQ(ierr); 100*c4762a1bSJed Brown ierr = ISDestroy(&from);CHKERRQ(ierr); 101*c4762a1bSJed Brown 102*c4762a1bSJed Brown /* 103*c4762a1bSJed Brown fill the global vector with the natural global numbering for each local entry 104*c4762a1bSJed Brown this is only done for testing purposes since it is easy to see if the scatter worked correctly 105*c4762a1bSJed Brown */ 106*c4762a1bSJed Brown ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 107*c4762a1bSJed Brown ierr = VecGetOwnershipRange(natural,&start,&end);CHKERRQ(ierr); 108*c4762a1bSJed Brown ierr = VecGetArray(natural,&xnatural);CHKERRQ(ierr); 109*c4762a1bSJed Brown for (i=start; i<end; i++) { 110*c4762a1bSJed Brown xnatural[i-start] = i; 111*c4762a1bSJed Brown } 112*c4762a1bSJed Brown ierr = VecRestoreArray(natural,&xnatural);CHKERRQ(ierr); 113*c4762a1bSJed Brown ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,global);CHKERRQ(ierr); 114*c4762a1bSJed Brown ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,global);CHKERRQ(ierr); 115*c4762a1bSJed Brown ierr = VecDestroy(&natural);CHKERRQ(ierr); 116*c4762a1bSJed Brown 117*c4762a1bSJed Brown /* scatter from global to local */ 118*c4762a1bSJed Brown ierr = VecScatterBegin(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 119*c4762a1bSJed Brown ierr = VecScatterEnd(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 120*c4762a1bSJed Brown ierr = VecScatterBegin(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 121*c4762a1bSJed Brown ierr = VecScatterEnd(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 122*c4762a1bSJed Brown /* 123*c4762a1bSJed Brown normally here you would also call 124*c4762a1bSJed Brown ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr); 125*c4762a1bSJed Brown ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr); 126*c4762a1bSJed Brown to update all the interior ghost cells between neighboring processes. 127*c4762a1bSJed Brown We don't do it here since this is only a test of "special" ghost points. 128*c4762a1bSJed Brown */ 129*c4762a1bSJed Brown 130*c4762a1bSJed Brown /* view each local ghosted vector */ 131*c4762a1bSJed Brown ierr = PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);CHKERRQ(ierr); 132*c4762a1bSJed Brown ierr = VecView(local,subviewer);CHKERRQ(ierr); 133*c4762a1bSJed Brown ierr = PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);CHKERRQ(ierr); 134*c4762a1bSJed Brown 135*c4762a1bSJed Brown ierr = VecScatterDestroy(&scatter1);CHKERRQ(ierr); 136*c4762a1bSJed Brown ierr = VecScatterDestroy(&scatter2);CHKERRQ(ierr); 137*c4762a1bSJed Brown ierr = VecDestroy(&local);CHKERRQ(ierr); 138*c4762a1bSJed Brown ierr = VecDestroy(&global);CHKERRQ(ierr); 139*c4762a1bSJed Brown ierr = DMDestroy(&da);CHKERRQ(ierr); 140*c4762a1bSJed Brown ierr = PetscFinalize(); 141*c4762a1bSJed Brown return ierr; 142*c4762a1bSJed Brown } 143*c4762a1bSJed Brown 144*c4762a1bSJed Brown 145*c4762a1bSJed Brown 146*c4762a1bSJed Brown /*TEST 147*c4762a1bSJed Brown 148*c4762a1bSJed Brown test: 149*c4762a1bSJed Brown 150*c4762a1bSJed Brown test: 151*c4762a1bSJed Brown suffix: 2 152*c4762a1bSJed Brown nsize: 2 153*c4762a1bSJed Brown 154*c4762a1bSJed Brown test: 155*c4762a1bSJed Brown suffix: 4 156*c4762a1bSJed Brown nsize: 4 157*c4762a1bSJed Brown 158*c4762a1bSJed Brown test: 159*c4762a1bSJed Brown suffix: 9 160*c4762a1bSJed Brown nsize: 9 161*c4762a1bSJed Brown 162*c4762a1bSJed Brown TEST*/ 163