Actual source code: ex6.c

      1: static char help[] = "\n\n";

      3: /*
      4:      Demonstrates using DM_BOUNDARY_GHOSTED how to handle a rotated boundary conditions where one edge
      5:     is connected to its immediate neighbor

      7:     Consider the domain (with natural numbering)

      9:      6   7   8
     10:      3   4   5
     11:      0   1   2

     13:     The ghost points along the bottom (directly below the three columns above) should be 0 3 and 6
     14:     while the ghost points along the left side should be 0 1 2

     16:     Note that the ghosted local vectors extend in both the x and y directions so, for example if we have a
     17:     single MPI process the ghosted vector has (in the original natural numbering)

     19:      x  x  x  x  x
     20:      2  6  7  8  x
     21:      1  3  4  5  x
     22:      0  0  1  2  x
     23:      x  0  3  6  x

     25:     where x indicates a location that is not updated by the communication and should be used.

     27:     For this to make sense the number of grid points in the x and y directions must be the same

     29:     This ghost point mapping was suggested by: Wenbo Zhao <zhaowenbo.npic@gmail.com>
     30: */

     32: #include <petscdm.h>
     33: #include <petscdmda.h>

     35: int main(int argc, char **argv)
     36: {
     37:   PetscInt     M = 6;
     38:   DM           da;
     39:   Vec          local, global, natural;
     40:   PetscInt     i, start, end, *ifrom, x, y, xm, ym;
     41:   PetscScalar *xnatural;
     42:   IS           from, to;
     43:   AO           ao;
     44:   VecScatter   scatter1, scatter2;
     45:   PetscViewer  subviewer;

     47:   PetscFunctionBeginUser;
     48:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

     50:   /* Create distributed array and get vectors */
     51:   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));
     52:   PetscCall(DMSetUp(da));
     53:   PetscCall(DMCreateGlobalVector(da, &global));
     54:   PetscCall(DMCreateLocalVector(da, &local));

     56:   /* construct global to local scatter for the left side of the domain to the ghost on the bottom */
     57:   PetscCall(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:     PetscCall(ISCreateStride(PETSC_COMM_SELF, xm, 1, 1, &to));
     60:   } else {
     61:     PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 0, &to));
     62:   }
     63:   PetscCall(PetscMalloc1(xm, &ifrom));
     64:   for (i = x; i < x + xm; i++) ifrom[i - x] = M * i;
     65:   PetscCall(DMDAGetAO(da, &ao));
     66:   PetscCall(AOApplicationToPetsc(ao, xm, ifrom));
     67:   if (!y) {
     68:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, xm, ifrom, PETSC_OWN_POINTER, &from));
     69:   } else {
     70:     PetscCall(PetscFree(ifrom));
     71:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_COPY_VALUES, &from));
     72:   }
     73:   PetscCall(VecScatterCreate(global, from, local, to, &scatter1));
     74:   PetscCall(ISDestroy(&to));
     75:   PetscCall(ISDestroy(&from));

     77:   /* construct global to local scatter for the bottom side of the domain to the ghost on the right */
     78:   if (!x) { /* only processes on the left side of the domain fill up the ghost locations */
     79:     PetscCall(ISCreateStride(PETSC_COMM_SELF, ym, xm + 2, xm + 2, &to));
     80:   } else {
     81:     PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 0, &to));
     82:   }
     83:   PetscCall(PetscMalloc1(ym, &ifrom));
     84:   for (i = y; i < y + ym; i++) ifrom[i - y] = i;
     85:   PetscCall(DMDAGetAO(da, &ao));
     86:   PetscCall(AOApplicationToPetsc(ao, ym, ifrom));
     87:   if (!x) {
     88:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ym, ifrom, PETSC_OWN_POINTER, &from));
     89:   } else {
     90:     PetscCall(PetscFree(ifrom));
     91:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_COPY_VALUES, &from));
     92:   }
     93:   PetscCall(VecScatterCreate(global, from, local, to, &scatter2));
     94:   PetscCall(ISDestroy(&to));
     95:   PetscCall(ISDestroy(&from));

     97:   /*
     98:      fill the global vector with the natural global numbering for each local entry
     99:      this is only done for testing purposes since it is easy to see if the scatter worked correctly
    100:   */
    101:   PetscCall(DMDACreateNaturalVector(da, &natural));
    102:   PetscCall(VecGetOwnershipRange(natural, &start, &end));
    103:   PetscCall(VecGetArray(natural, &xnatural));
    104:   for (i = start; i < end; i++) xnatural[i - start] = i;
    105:   PetscCall(VecRestoreArray(natural, &xnatural));
    106:   PetscCall(DMDANaturalToGlobalBegin(da, natural, INSERT_VALUES, global));
    107:   PetscCall(DMDANaturalToGlobalEnd(da, natural, INSERT_VALUES, global));
    108:   PetscCall(VecDestroy(&natural));

    110:   /* scatter from global to local */
    111:   PetscCall(VecScatterBegin(scatter1, global, local, INSERT_VALUES, SCATTER_FORWARD));
    112:   PetscCall(VecScatterEnd(scatter1, global, local, INSERT_VALUES, SCATTER_FORWARD));
    113:   PetscCall(VecScatterBegin(scatter2, global, local, INSERT_VALUES, SCATTER_FORWARD));
    114:   PetscCall(VecScatterEnd(scatter2, global, local, INSERT_VALUES, SCATTER_FORWARD));
    115:   /*
    116:      normally here you would also call
    117:   PetscCall(DMGlobalToLocalBegin(da,global,INSERT_VALUES,local));
    118:   PetscCall(DMGlobalToLocalEnd(da,global,INSERT_VALUES,local));
    119:     to update all the interior ghost cells between neighboring processes.
    120:     We don't do it here since this is only a test of "special" ghost points.
    121:   */

    123:   /* view each local ghosted vector */
    124:   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &subviewer));
    125:   PetscCall(VecView(local, subviewer));
    126:   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &subviewer));

    128:   PetscCall(VecScatterDestroy(&scatter1));
    129:   PetscCall(VecScatterDestroy(&scatter2));
    130:   PetscCall(VecDestroy(&local));
    131:   PetscCall(VecDestroy(&global));
    132:   PetscCall(DMDestroy(&da));
    133:   PetscCall(PetscFinalize());
    134:   return 0;
    135: }

    137: /*TEST

    139:    test:

    141:    test:
    142:       suffix: 2
    143:       nsize: 2

    145:    test:
    146:       suffix: 4
    147:       nsize: 4

    149:    test:
    150:       suffix: 9
    151:       nsize: 9

    153: TEST*/