Actual source code: ex6.c

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

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

      8:     Consider the domain (with natural numbering)

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

     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

     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)

     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

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

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

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

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

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

     49:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

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

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

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

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

    117:   /* scatter from global to local */
    118:   VecScatterBegin(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);
    119:   VecScatterEnd(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);
    120:   VecScatterBegin(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);
    121:   VecScatterEnd(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);
    122:   /*
    123:      normally here you would also call
    124:   DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);
    125:   DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);
    126:     to update all the interior ghost cells between neighboring processes.
    127:     We don't do it here since this is only a test of "special" ghost points.
    128:   */

    130:   /* view each local ghosted vector */
    131:   PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);
    132:   VecView(local,subviewer);
    133:   PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);

    135:   VecScatterDestroy(&scatter1);
    136:   VecScatterDestroy(&scatter2);
    137:   VecDestroy(&local);
    138:   VecDestroy(&global);
    139:   DMDestroy(&da);
    140:   PetscFinalize();
    141:   return ierr;
    142: }

    146: /*TEST

    148:    test:

    150:    test:
    151:       suffix: 2
    152:       nsize: 2

    154:    test:
    155:       suffix: 4
    156:       nsize: 4

    158:    test:
    159:       suffix: 9
    160:       nsize: 9

    162: TEST*/