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