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 PetscFunctionBeginUser; 49 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 50 51 /* Create distributed array and get vectors */ 52 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DMDA_STENCIL_STAR, M, M, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da)); 53 PetscCall(DMSetUp(da)); 54 PetscCall(DMCreateGlobalVector(da, &global)); 55 PetscCall(DMCreateLocalVector(da, &local)); 56 57 /* construct global to local scatter for the left side of the domain to the ghost on the bottom */ 58 PetscCall(DMDAGetCorners(da, &x, &y, NULL, &xm, &ym, NULL)); 59 if (!y) { /* only processes on the bottom of the domain fill up the ghost locations */ 60 PetscCall(ISCreateStride(PETSC_COMM_SELF, xm, 1, 1, &to)); 61 } else { 62 PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 0, &to)); 63 } 64 PetscCall(PetscMalloc1(xm, &ifrom)); 65 for (i = x; i < x + xm; i++) ifrom[i - x] = M * i; 66 PetscCall(DMDAGetAO(da, &ao)); 67 PetscCall(AOApplicationToPetsc(ao, xm, ifrom)); 68 if (!y) { 69 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, xm, ifrom, PETSC_OWN_POINTER, &from)); 70 } else { 71 PetscCall(PetscFree(ifrom)); 72 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_COPY_VALUES, &from)); 73 } 74 PetscCall(VecScatterCreate(global, from, local, to, &scatter1)); 75 PetscCall(ISDestroy(&to)); 76 PetscCall(ISDestroy(&from)); 77 78 /* construct global to local scatter for the bottom side of the domain to the ghost on the right */ 79 if (!x) { /* only processes on the left side of the domain fill up the ghost locations */ 80 PetscCall(ISCreateStride(PETSC_COMM_SELF, ym, xm + 2, xm + 2, &to)); 81 } else { 82 PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 0, &to)); 83 } 84 PetscCall(PetscMalloc1(ym, &ifrom)); 85 for (i = y; i < y + ym; i++) ifrom[i - y] = i; 86 PetscCall(DMDAGetAO(da, &ao)); 87 PetscCall(AOApplicationToPetsc(ao, ym, ifrom)); 88 if (!x) { 89 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ym, ifrom, PETSC_OWN_POINTER, &from)); 90 } else { 91 PetscCall(PetscFree(ifrom)); 92 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_COPY_VALUES, &from)); 93 } 94 PetscCall(VecScatterCreate(global, from, local, to, &scatter2)); 95 PetscCall(ISDestroy(&to)); 96 PetscCall(ISDestroy(&from)); 97 98 /* 99 fill the global vector with the natural global numbering for each local entry 100 this is only done for testing purposes since it is easy to see if the scatter worked correctly 101 */ 102 PetscCall(DMDACreateNaturalVector(da, &natural)); 103 PetscCall(VecGetOwnershipRange(natural, &start, &end)); 104 PetscCall(VecGetArray(natural, &xnatural)); 105 for (i = start; i < end; i++) xnatural[i - start] = i; 106 PetscCall(VecRestoreArray(natural, &xnatural)); 107 PetscCall(DMDANaturalToGlobalBegin(da, natural, INSERT_VALUES, global)); 108 PetscCall(DMDANaturalToGlobalEnd(da, natural, INSERT_VALUES, global)); 109 PetscCall(VecDestroy(&natural)); 110 111 /* scatter from global to local */ 112 PetscCall(VecScatterBegin(scatter1, global, local, INSERT_VALUES, SCATTER_FORWARD)); 113 PetscCall(VecScatterEnd(scatter1, global, local, INSERT_VALUES, SCATTER_FORWARD)); 114 PetscCall(VecScatterBegin(scatter2, global, local, INSERT_VALUES, SCATTER_FORWARD)); 115 PetscCall(VecScatterEnd(scatter2, global, local, INSERT_VALUES, SCATTER_FORWARD)); 116 /* 117 normally here you would also call 118 PetscCall(DMGlobalToLocalBegin(da,global,INSERT_VALUES,local)); 119 PetscCall(DMGlobalToLocalEnd(da,global,INSERT_VALUES,local)); 120 to update all the interior ghost cells between neighboring processes. 121 We don't do it here since this is only a test of "special" ghost points. 122 */ 123 124 /* view each local ghosted vector */ 125 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &subviewer)); 126 PetscCall(VecView(local, subviewer)); 127 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &subviewer)); 128 129 PetscCall(VecScatterDestroy(&scatter1)); 130 PetscCall(VecScatterDestroy(&scatter2)); 131 PetscCall(VecDestroy(&local)); 132 PetscCall(VecDestroy(&global)); 133 PetscCall(DMDestroy(&da)); 134 PetscCall(PetscFinalize()); 135 return 0; 136 } 137 138 /*TEST 139 140 test: 141 142 test: 143 suffix: 2 144 nsize: 2 145 146 test: 147 suffix: 4 148 nsize: 4 149 150 test: 151 suffix: 9 152 nsize: 9 153 154 TEST*/ 155