Actual source code: ex6.c
petsc-3.4.2 2013-07-02

      2: static char help[] = "Demonstrates using 3 DMDA's to manage a slightly non-trivial grid";

      4: #include <petscdmda.h>
      5: #include <petscdraw.h>

      7: struct _p_FA {
      8:   MPI_Comm   comm[3];
      9:   PetscInt   xl[3],yl[3],ml[3],nl[3];    /* corners and sizes of local vector in DMDA */
     10:   PetscInt   xg[3],yg[3],mg[3],ng[3];    /* corners and sizes of global vector in DMDA */
     11:   PetscInt   offl[3],offg[3];            /* offset in local and global vector of region 1, 2 and 3 portions */
     12:   Vec        g,l;
     13:   VecScatter vscat;
     14:   PetscInt   p1,p2,r1,r2,r1g,r2g,sw;
     15: };
     16: typedef struct _p_FA *FA;

     18: typedef struct {
     19:   PetscScalar X;
     20:   PetscScalar Y;
     21: } Field;

     23: PetscErrorCode FAGetLocalCorners(FA fa,PetscInt j,PetscInt *x,PetscInt *y,PetscInt *m,PetscInt *n)
     24: {
     26:   if (fa->comm[j]) {
     27:     *x = fa->xl[j];
     28:     *y = fa->yl[j];
     29:     *m = fa->ml[j];
     30:     *n = fa->nl[j];
     31:   } else {
     32:     *x = *y = *m = *n = 0;
     33:   }
     34:   return(0);
     35: }

     37: PetscErrorCode FAGetGlobalCorners(FA fa,PetscInt j,PetscInt *x,PetscInt *y,PetscInt *m,PetscInt *n)
     38: {
     40:   if (fa->comm[j]) {
     41:     *x = fa->xg[j];
     42:     *y = fa->yg[j];
     43:     *m = fa->mg[j];
     44:     *n = fa->ng[j];
     45:   } else {
     46:     *x = *y = *m = *n = 0;
     47:   }
     48:   return(0);
     49: }

     51: PetscErrorCode FAGetLocalArray(FA fa,Vec v,PetscInt j,Field ***f)
     52: {
     54:   PetscScalar    *va;
     55:   PetscInt       i;
     56:   Field          **a;

     59:   if (fa->comm[j]) {
     60:     VecGetArray(v,&va);
     61:     PetscMalloc(fa->nl[j]*sizeof(Field*),&a);
     62:     for (i=0; i<fa->nl[j]; i++) (a)[i] = (Field*) (va + 2*fa->offl[j] + i*2*fa->ml[j] - 2*fa->xl[j]);
     63:     *f   = a - fa->yl[j];
     64:     VecRestoreArray(v,&va);
     65:   } else {
     66:     *f = 0;
     67:   }
     68:   return(0);
     69: }

     71: PetscErrorCode FARestoreLocalArray(FA fa,Vec v,PetscInt j,Field ***f)
     72: {
     74:   void           *dummy;

     77:   if (fa->comm[j]) {
     78:     dummy = *f + fa->yl[j];
     79:     PetscFree(dummy);
     80:   }
     81:   return(0);
     82: }

     84: PetscErrorCode FAGetGlobalArray(FA fa,Vec v,PetscInt j,Field ***f)
     85: {
     87:   PetscScalar    *va;
     88:   PetscInt       i;
     89:   Field          **a;

     92:   if (fa->comm[j]) {
     93:     VecGetArray(v,&va);
     94:     PetscMalloc(fa->ng[j]*sizeof(Field*),&a);
     95:     for (i=0; i<fa->ng[j]; i++) (a)[i] = (Field*) (va + 2*fa->offg[j] + i*2*fa->mg[j] - 2*fa->xg[j]);
     96:     *f   = a - fa->yg[j];
     97:     VecRestoreArray(v,&va);
     98:   } else {
     99:     *f = 0;
    100:   }
    101:   return(0);
    102: }

    104: PetscErrorCode FARestoreGlobalArray(FA fa,Vec v,PetscInt j,Field ***f)
    105: {
    107:   void           *dummy;

    110:   if (fa->comm[j]) {
    111:     dummy = *f + fa->yg[j];
    112:     PetscFree(dummy);
    113:   }
    114:   return(0);
    115: }

    117: PetscErrorCode FAGetGlobalVector(FA fa,Vec *v)
    118: {

    122:   VecDuplicate(fa->g,v);
    123:   return(0);
    124: }

    126: PetscErrorCode FAGetLocalVector(FA fa,Vec *v)
    127: {

    131:   VecDuplicate(fa->l,v);
    132:   return(0);
    133: }

    135: PetscErrorCode FAGlobalToLocal(FA fa,Vec g,Vec l)
    136: {

    140:   VecScatterBegin(fa->vscat,g,l,INSERT_VALUES,SCATTER_FORWARD);
    141:   VecScatterEnd(fa->vscat,g,l,INSERT_VALUES,SCATTER_FORWARD);
    142:   return(0);
    143: }

    145: PetscErrorCode FADestroy(FA *fa)
    146: {

    150:   VecDestroy(&(*fa)->g);
    151:   VecDestroy(&(*fa)->l);
    152:   VecScatterDestroy(&(*fa)->vscat);
    153:   PetscFree(*fa);
    154:   return(0);
    155: }

    157: PetscErrorCode FACreate(FA *infa)
    158: {
    159:   FA             fa;
    160:   PetscMPIInt    rank;
    161:   PetscInt       tonglobal,globalrstart,x,nx,y,ny,*tonatural,i,j,*to,*from,offt[3];
    162:   PetscInt       *fromnatural,fromnglobal,nscat,nlocal,cntl1,cntl2,cntl3,*indices;

    165:   /* Each DMDA manages the local vector for the portion of region 1, 2, and 3 for that processor
    166:      Each DMDA can belong on any subset (overlapping between DMDA's or not) of processors
    167:      For processes that a particular DMDA does not exist on, the corresponding comm should be set to zero
    168:   */
    169:   DM da1 = 0,da2 = 0,da3 = 0;
    170:   /*
    171:       v1, v2, v3 represent the local vector for a single DMDA
    172:   */
    173:   Vec vl1 = 0,vl2 = 0,vl3 = 0, vg1 = 0, vg2 = 0,vg3 = 0;

    175:   /*
    176:      globalvec and friends represent the global vectors that are used for the PETSc solvers
    177:      localvec represents the concatenation of the (up to) 3 local vectors; vl1, vl2, vl3

    179:      tovec and friends represent intermediate vectors that are ONLY used for setting up the
    180:      final communication patterns. Once this setup routine is complete they are destroyed.
    181:      The tovec  is like the globalvec EXCEPT it has redundant locations for the ghost points
    182:      between regions 2+3 and 1.
    183:   */
    184:   AO          toao,globalao;
    185:   IS          tois,globalis,is;
    186:   Vec         tovec,globalvec,localvec;
    187:   VecScatter  vscat;
    188:   PetscScalar *globalarray,*localarray,*toarray;

    190:   PetscNew(struct _p_FA,&fa);
    191:   /*
    192:       fa->sw is the stencil width

    194:       fa->p1 is the width of region 1, fa->p2 the width of region 2 (must be the same)
    195:       fa->r1 height of region 1
    196:       fa->r2 height of region 2

    198:       fa->r2 is also the height of region 3-4
    199:       (fa->p1 - fa->p2)/2 is the width of both region 3 and region 4
    200:   */
    201:   fa->p1  = 24;
    202:   fa->p2  = 15;
    203:   fa->r1  = 6;
    204:   fa->r2  = 6;
    205:   fa->sw  = 1;
    206:   fa->r1g = fa->r1 + fa->sw;
    207:   fa->r2g = fa->r2 + fa->sw;

    209:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

    211:   fa->comm[0] = PETSC_COMM_WORLD;
    212:   fa->comm[1] = PETSC_COMM_WORLD;
    213:   fa->comm[2] = PETSC_COMM_WORLD;
    214:   /* Test case with different communicators */
    215:   /* Normally one would use MPI_Comm routines to build MPI communicators on which you wish to partition the DMDAs*/
    216:   /*
    217:   if (!rank) {
    218:     fa->comm[0] = PETSC_COMM_SELF;
    219:     fa->comm[1] = 0;
    220:     fa->comm[2] = 0;
    221:   } else if (rank == 1) {
    222:     fa->comm[0] = 0;
    223:     fa->comm[1] = PETSC_COMM_SELF;
    224:     fa->comm[2] = 0;
    225:   } else {
    226:     fa->comm[0] = 0;
    227:     fa->comm[1] = 0;
    228:     fa->comm[2] = PETSC_COMM_SELF;
    229:   } */

    231:   if (fa->p2 > fa->p1 - 3) SETERRQ(PETSC_COMM_SELF,1,"Width of region fa->p2 must be at least 3 less then width of region 1");
    232:   if (!((fa->p2 - fa->p1) % 2)) SETERRQ(PETSC_COMM_SELF,1,"width of region 3 must NOT be divisible by 2");

    234:   if (fa->comm[1]) {
    235:     DMDACreate2d(fa->comm[1],DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,fa->p2,fa->r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa->sw,NULL,NULL,&da2);
    236:     DMGetLocalVector(da2,&vl2);
    237:     DMGetGlobalVector(da2,&vg2);
    238:   }
    239:   if (fa->comm[2]) {
    240:     DMDACreate2d(fa->comm[2],DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,fa->p1-fa->p2,fa->r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa->sw,NULL,NULL,&da3);
    241:     DMGetLocalVector(da3,&vl3);
    242:     DMGetGlobalVector(da3,&vg3);
    243:   }
    244:   if (fa->comm[0]) {
    245:     DMDACreate2d(fa->comm[0],DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,fa->p1,fa->r1g,PETSC_DECIDE,PETSC_DECIDE,1,fa->sw,NULL,NULL,&da1);
    246:     DMGetLocalVector(da1,&vl1);
    247:     DMGetGlobalVector(da1,&vg1);
    248:   }

    250:   /* count the number of unknowns owned on each processor and determine the starting point of each processors ownership
    251:      for global vector with redundancy */
    252:   tonglobal = 0;
    253:   if (fa->comm[1]) {
    254:     DMDAGetCorners(da2,&x,&y,0,&nx,&ny,0);
    255:     tonglobal += nx*ny;
    256:   }
    257:   if (fa->comm[2]) {
    258:     DMDAGetCorners(da3,&x,&y,0,&nx,&ny,0);
    259:     tonglobal += nx*ny;
    260:   }
    261:   if (fa->comm[0]) {
    262:     DMDAGetCorners(da1,&x,&y,0,&nx,&ny,0);
    263:     tonglobal += nx*ny;
    264:   }
    265:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Number of unknowns owned %d\n",rank,tonglobal);
    266:   PetscSynchronizedFlush(PETSC_COMM_WORLD);

    268:   /* Get tonatural number for each node */
    269:   PetscMalloc((tonglobal+1)*sizeof(PetscInt),&tonatural);
    270:   tonglobal = 0;
    271:   if (fa->comm[1]) {
    272:     DMDAGetCorners(da2,&x,&y,0,&nx,&ny,0);
    273:     for (j=0; j<ny; j++) {
    274:       for (i=0; i<nx; i++) {
    275:         tonatural[tonglobal++] = (fa->p1 - fa->p2)/2 + x + i + fa->p1*(y + j);
    276:       }
    277:     }
    278:   }
    279:   if (fa->comm[2]) {
    280:     DMDAGetCorners(da3,&x,&y,0,&nx,&ny,0);
    281:     for (j=0; j<ny; j++) {
    282:       for (i=0; i<nx; i++) {
    283:         if (x + i < (fa->p1 - fa->p2)/2) tonatural[tonglobal++] = x + i + fa->p1*(y + j);
    284:         else tonatural[tonglobal++] = fa->p2 + x + i + fa->p1*(y + j);
    285:       }
    286:     }
    287:   }
    288:   if (fa->comm[0]) {
    289:     DMDAGetCorners(da1,&x,&y,0,&nx,&ny,0);
    290:     for (j=0; j<ny; j++) {
    291:       for (i=0; i<nx; i++) {
    292:         tonatural[tonglobal++] = fa->p1*fa->r2g + x + i + fa->p1*(y + j);
    293:       }
    294:     }
    295:   }
    296:   /*  PetscIntView(tonglobal,tonatural,PETSC_VIEWER_STDOUT_WORLD); */
    297:   AOCreateBasic(PETSC_COMM_WORLD,tonglobal,tonatural,0,&toao);
    298:   PetscFree(tonatural);

    300:   /* count the number of unknowns owned on each processor and determine the starting point of each processors ownership
    301:      for global vector without redundancy */
    302:   fromnglobal = 0;
    303:   fa->offg[1] = 0;
    304:   offt[1]     = 0;
    305:   if (fa->comm[1]) {
    306:     DMDAGetCorners(da2,&x,&y,0,&nx,&ny,0);
    307:     offt[2] = nx*ny;
    308:     if (y+ny == fa->r2g) ny--;    /* includes the ghost points on the upper side */
    309:     fromnglobal += nx*ny;
    310:     fa->offg[2]  = fromnglobal;
    311:   } else {
    312:     offt[2]     = 0;
    313:     fa->offg[2] = 0;
    314:   }
    315:   if (fa->comm[2]) {
    316:     DMDAGetCorners(da3,&x,&y,0,&nx,&ny,0);
    317:     offt[0] = offt[2] + nx*ny;
    318:     if (y+ny == fa->r2g) ny--;    /* includes the ghost points on the upper side */
    319:     fromnglobal += nx*ny;
    320:     fa->offg[0]  = fromnglobal;
    321:   } else {
    322:     offt[0]     = offt[2];
    323:     fa->offg[0] = fromnglobal;
    324:   }
    325:   if (fa->comm[0]) {
    326:     DMDAGetCorners(da1,&x,&y,0,&nx,&ny,0);
    327:     if (y == 0) ny--;    /* includes the ghost points on the lower side */
    328:     fromnglobal += nx*ny;
    329:   }
    330:   MPI_Scan(&fromnglobal,&globalrstart,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
    331:   globalrstart -= fromnglobal;
    332:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Number of unknowns owned %d\n",rank,fromnglobal);
    333:   PetscSynchronizedFlush(PETSC_COMM_WORLD);

    335:   /* Get fromnatural number for each node */
    336:   PetscMalloc((fromnglobal+1)*sizeof(PetscInt),&fromnatural);
    337:   fromnglobal = 0;
    338:   if (fa->comm[1]) {
    339:     DMDAGetCorners(da2,&x,&y,0,&nx,&ny,0);
    340:     if (y+ny == fa->r2g) ny--;    /* includes the ghost points on the upper side */
    341:     fa->xg[1] = x; fa->yg[1] = y; fa->mg[1] = nx; fa->ng[1] = ny;
    342:     DMDAGetGhostCorners(da2,&fa->xl[1],&fa->yl[1],0,&fa->ml[1],&fa->nl[1],0);
    343:     for (j=0; j<ny; j++) {
    344:       for (i=0; i<nx; i++) {
    345:         fromnatural[fromnglobal++] = (fa->p1 - fa->p2)/2 + x + i + fa->p1*(y + j);
    346:       }
    347:     }
    348:   }
    349:   if (fa->comm[2]) {
    350:     DMDAGetCorners(da3,&x,&y,0,&nx,&ny,0);
    351:     if (y+ny == fa->r2g) ny--;    /* includes the ghost points on the upper side */
    352:     fa->xg[2] = x; fa->yg[2] = y; fa->mg[2] = nx; fa->ng[2] = ny;
    353:     DMDAGetGhostCorners(da3,&fa->xl[2],&fa->yl[2],0,&fa->ml[2],&fa->nl[2],0);
    354:     for (j=0; j<ny; j++) {
    355:       for (i=0; i<nx; i++) {
    356:         if (x + i < (fa->p1 - fa->p2)/2) fromnatural[fromnglobal++] = x + i + fa->p1*(y + j);
    357:         else fromnatural[fromnglobal++] = fa->p2 + x + i + fa->p1*(y + j);
    358:       }
    359:     }
    360:   }
    361:   if (fa->comm[0]) {
    362:     DMDAGetCorners(da1,&x,&y,0,&nx,&ny,0);
    363:     if (y == 0) ny--;    /* includes the ghost points on the lower side */
    364:     else y--;
    365:     fa->xg[0] = x; fa->yg[0] = y; fa->mg[0] = nx; fa->ng[0] = ny;
    366:     DMDAGetGhostCorners(da1,&fa->xl[0],&fa->yl[0],0,&fa->ml[0],&fa->nl[0],0);
    367:     for (j=0; j<ny; j++) {
    368:       for (i=0; i<nx; i++) {
    369:         fromnatural[fromnglobal++] = fa->p1*fa->r2 + x + i + fa->p1*(y + j);
    370:       }
    371:     }
    372:   }
    373:   /*PetscIntView(fromnglobal,fromnatural,PETSC_VIEWER_STDOUT_WORLD);*/
    374:   AOCreateBasic(PETSC_COMM_WORLD,fromnglobal,fromnatural,0,&globalao);
    375:   PetscFree(fromnatural);

    377:   /* ---------------------------------------------------*/
    378:   /* Create the scatter that updates 1 from 2 and 3 and 3 and 2 from 1 */
    379:   /* currently handles stencil width of 1 ONLY */
    380:   PetscMalloc(tonglobal*sizeof(PetscInt),&to);
    381:   PetscMalloc(tonglobal*sizeof(PetscInt),&from);
    382:   nscat = 0;
    383:   if (fa->comm[1]) {
    384:     DMDAGetCorners(da2,&x,&y,0,&nx,&ny,0);
    385:     for (j=0; j<ny; j++) {
    386:       for (i=0; i<nx; i++) {
    387:         to[nscat] = from[nscat] = (fa->p1 - fa->p2)/2 + x + i + fa->p1*(y + j);nscat++;
    388:       }
    389:     }
    390:   }
    391:   if (fa->comm[2]) {
    392:     DMDAGetCorners(da3,&x,&y,0,&nx,&ny,0);
    393:     for (j=0; j<ny; j++) {
    394:       for (i=0; i<nx; i++) {
    395:         if (x + i < (fa->p1 - fa->p2)/2) {
    396:           to[nscat] = from[nscat] = x + i + fa->p1*(y + j);nscat++;
    397:         } else {
    398:           to[nscat] = from[nscat] = fa->p2 + x + i + fa->p1*(y + j);nscat++;
    399:         }
    400:       }
    401:     }
    402:   }
    403:   if (fa->comm[0]) {
    404:     DMDAGetCorners(da1,&x,&y,0,&nx,&ny,0);
    405:     for (j=0; j<ny; j++) {
    406:       for (i=0; i<nx; i++) {
    407:         to[nscat]     = fa->p1*fa->r2g + x + i + fa->p1*(y + j);
    408:         from[nscat++] = fa->p1*(fa->r2 - 1) + x + i + fa->p1*(y + j);
    409:       }
    410:     }
    411:   }
    412:   AOApplicationToPetsc(toao,nscat,to);
    413:   AOApplicationToPetsc(globalao,nscat,from);
    414:   ISCreateGeneral(PETSC_COMM_WORLD,nscat,to,PETSC_COPY_VALUES,&tois);
    415:   ISCreateGeneral(PETSC_COMM_WORLD,nscat,from,PETSC_COPY_VALUES,&globalis);
    416:   PetscFree(to);
    417:   PetscFree(from);
    418:   VecCreateMPI(PETSC_COMM_WORLD,tonglobal,PETSC_DETERMINE,&tovec);
    419:   VecCreateMPI(PETSC_COMM_WORLD,fromnglobal,PETSC_DETERMINE,&globalvec);
    420:   VecScatterCreate(globalvec,globalis,tovec,tois,&vscat);
    421:   ISDestroy(&tois);
    422:   ISDestroy(&globalis);
    423:   AODestroy(&globalao);
    424:   AODestroy(&toao);

    426:   /* fill up global vector without redundant values with PETSc global numbering */
    427:   VecGetArray(globalvec,&globalarray);
    428:   for (i=0; i<fromnglobal; i++) {
    429:     globalarray[i] = globalrstart + i;
    430:   }
    431:   VecRestoreArray(globalvec,&globalarray);

    433:   /* scatter PETSc global indices to redundant valueed array */
    434:   VecScatterBegin(vscat,globalvec,tovec,INSERT_VALUES,SCATTER_FORWARD);
    435:   VecScatterEnd(vscat,globalvec,tovec,INSERT_VALUES,SCATTER_FORWARD);

    437:   /* Create local vector that is the concatenation of the local vectors */
    438:   nlocal = 0;
    439:   cntl1  = cntl2 = cntl3 = 0;
    440:   if (fa->comm[1]) {
    441:     VecGetSize(vl2,&cntl2);
    442:     nlocal += cntl2;
    443:   }
    444:   if (fa->comm[2]) {
    445:     VecGetSize(vl3,&cntl3);
    446:     nlocal += cntl3;
    447:   }
    448:   if (fa->comm[0]) {
    449:     VecGetSize(vl1,&cntl1);
    450:     nlocal += cntl1;
    451:   }
    452:   fa->offl[0] = cntl2 + cntl3;
    453:   fa->offl[1] = 0;
    454:   fa->offl[2] = cntl2;
    455:   VecCreateSeq(PETSC_COMM_SELF,nlocal,&localvec);

    457:   /* cheat so that  vl1, vl2, vl3 shared array memory with localvec */
    458:   VecGetArray(localvec,&localarray);
    459:   VecGetArray(tovec,&toarray);
    460:   if (fa->comm[1]) {
    461:     VecPlaceArray(vl2,localarray+fa->offl[1]);
    462:     VecPlaceArray(vg2,toarray+offt[1]);
    463:     DMGlobalToLocalBegin(da2,vg2,INSERT_VALUES,vl2);
    464:     DMGlobalToLocalEnd(da2,vg2,INSERT_VALUES,vl2);
    465:     DMRestoreGlobalVector(da2,&vg2);
    466:   }
    467:   if (fa->comm[2]) {
    468:     VecPlaceArray(vl3,localarray+fa->offl[2]);
    469:     VecPlaceArray(vg3,toarray+offt[2]);
    470:     DMGlobalToLocalBegin(da3,vg3,INSERT_VALUES,vl3);
    471:     DMGlobalToLocalEnd(da3,vg3,INSERT_VALUES,vl3);
    472:     DMRestoreGlobalVector(da3,&vg3);
    473:   }
    474:   if (fa->comm[0]) {
    475:     VecPlaceArray(vl1,localarray+fa->offl[0]);
    476:     VecPlaceArray(vg1,toarray+offt[0]);
    477:     DMGlobalToLocalBegin(da1,vg1,INSERT_VALUES,vl1);
    478:     DMGlobalToLocalEnd(da1,vg1,INSERT_VALUES,vl1);
    479:     DMRestoreGlobalVector(da1,&vg1);
    480:   }
    481:   VecRestoreArray(localvec,&localarray);
    482:   VecRestoreArray(tovec,&toarray);

    484:   /* no longer need the redundant vector and VecScatter to it */
    485:   VecScatterDestroy(&vscat);
    486:   VecDestroy(&tovec);

    488:   /* Create final scatter that goes directly from globalvec to localvec */
    489:   /* this is the one to be used in the application code */
    490:   PetscMalloc(nlocal*sizeof(PetscInt),&indices);
    491:   VecGetArray(localvec,&localarray);
    492:   for (i=0; i<nlocal; i++) {
    493:     indices[i] = (PetscInt) (localarray[i]);
    494:   }
    495:   VecRestoreArray(localvec,&localarray);
    496:   ISCreateBlock(PETSC_COMM_WORLD,2,nlocal,indices,PETSC_COPY_VALUES,&is);
    497:   PetscFree(indices);

    499:   VecCreateSeq(PETSC_COMM_SELF,2*nlocal,&fa->l);
    500:   VecCreateMPI(PETSC_COMM_WORLD,2*fromnglobal,PETSC_DETERMINE,&fa->g);

    502:   VecScatterCreate(fa->g,is,fa->l,NULL,&fa->vscat);
    503:   ISDestroy(&is);

    505:   VecDestroy(&globalvec);
    506:   VecDestroy(&localvec);
    507:   if (fa->comm[0]) {
    508:     DMRestoreLocalVector(da1,&vl1);
    509:     DMDestroy(&da1);
    510:   }
    511:   if (fa->comm[1]) {
    512:     DMRestoreLocalVector(da2,&vl2);
    513:     DMDestroy(&da2);
    514:   }
    515:   if (fa->comm[2]) {
    516:     DMRestoreLocalVector(da3,&vl3);
    517:     DMDestroy(&da3);
    518:   }
    519:   *infa = fa;
    520:   return(0);
    521: }

    523: /* Crude graphics to test that the ghost points are properly updated */
    524: #include <petscdraw.h>

    526: typedef struct {
    527:   PetscInt    m[3],n[3];
    528:   PetscScalar *xy[3];
    529: } ZoomCtx;

    531: PetscErrorCode DrawPatch(PetscDraw draw,void *ctx)
    532: {
    533:   ZoomCtx        *zctx = (ZoomCtx*)ctx;
    535:   PetscInt       m,n,i,j,id,k;
    536:   PetscReal      x1,x2,x3,x4,y_1,y2,y3,y4;
    537:   PetscScalar    *xy;

    540:   for (k=0; k<3; k++) {
    541:     m  = zctx->m[k];
    542:     n  = zctx->n[k];
    543:     xy = zctx->xy[k];

    545:     for (j=0; j<n-1; j++) {
    546:       for (i=0; i<m-1; i++) {
    547:         id   = i+j*m;    x1 = xy[2*id];y_1 = xy[2*id+1];
    548:         id   = i+j*m+1;  x2 = xy[2*id];y2  = xy[2*id+1];
    549:         id   = i+j*m+1+m;x3 = xy[2*id];y3  = xy[2*id+1];
    550:         id   = i+j*m+m;  x4 = xy[2*id];y4  = xy[2*id+1];
    551:         PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);
    552:         PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);
    553:         PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);
    554:         PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);
    555:       }
    556:     }
    557:   }
    558:   PetscDrawFlush(draw);
    559:   return(0);
    560: }

    562: PetscErrorCode DrawFA(FA fa,Vec v)
    563: {
    565:   PetscScalar    *va;
    566:   ZoomCtx        zctx;
    567:   PetscReal      xmint = 10000.0,xmaxt = -10000.0,ymint = 100000.0,ymaxt = -10000.0;
    568:   PetscReal      xmin,xmax,ymin,ymax;
    569:   PetscInt       i,vn,ln,j;

    572:   VecGetArray(v,&va);
    573:   VecGetSize(v,&vn);
    574:   VecGetSize(fa->l,&ln);
    575:   for (j=0; j<3; j++) {
    576:     if (vn == ln) {
    577:       zctx.xy[j] = va + 2*fa->offl[j];
    578:       zctx.m[j]  = fa->ml[j];
    579:       zctx.n[j]  = fa->nl[j];
    580:     } else {
    581:       zctx.xy[j] = va + 2*fa->offg[j];
    582:       zctx.m[j]  = fa->mg[j];
    583:       zctx.n[j]  = fa->ng[j];
    584:     }
    585:     for (i=0; i<zctx.m[j]*zctx.n[j]; i++) {
    586:       if (zctx.xy[j][2*i] > xmax) xmax = zctx.xy[j][2*i];
    587:       if (zctx.xy[j][2*i] < xmin) xmin = zctx.xy[j][2*i];
    588:       if (zctx.xy[j][2*i+1] > ymax) ymax = zctx.xy[j][2*i+1];
    589:       if (zctx.xy[j][2*i+1] < ymin) ymin = zctx.xy[j][2*i+1];
    590:     }
    591:   }
    592:   MPI_Allreduce(&xmin,&xmint,1,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD);
    593:   MPI_Allreduce(&xmax,&xmaxt,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD);
    594:   MPI_Allreduce(&ymin,&ymint,1,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD);
    595:   MPI_Allreduce(&ymax,&ymaxt,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD);
    596:   xmin = xmint - .2*(xmaxt - xmint);
    597:   xmax = xmaxt + .2*(xmaxt - xmint);
    598:   ymin = ymint - .2*(ymaxt - ymint);
    599:   ymax = ymaxt + .2*(ymaxt - ymint);
    600: #if defined(PETSC_HAVE_X) || defined(PETSC_HAVE_OPENGL)
    601:   {
    602:     PetscDraw      draw;

    604:     PetscDrawCreate(PETSC_COMM_WORLD,0,"meshes",PETSC_DECIDE,PETSC_DECIDE,700,700,&draw);
    605:     PetscDrawSetFromOptions(draw);
    606:     PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
    607:     PetscDrawZoom(draw,DrawPatch,&zctx);
    608:     VecRestoreArray(v,&va);
    609:     PetscDrawDestroy(&draw);
    610:   }
    611: #endif
    612:   return(0);
    613: }

    615: /* crude mappings from rectangular arrays to the true geometry. These are ONLY for testing!
    616:    they will not be used the actual code */
    617: PetscErrorCode FAMapRegion3(FA fa,Vec g)
    618: {
    620:   PetscReal      R = 1.0,Rscale,Ascale;
    621:   PetscInt       i,k,x,y,m,n;
    622:   Field          **ga;

    625:   Rscale = R/(fa->r2-1);
    626:   Ascale = 2.0*PETSC_PI/(3.0*(fa->p1 - fa->p2 - 1));

    628:   FAGetGlobalCorners(fa,2,&x,&y,&m,&n);
    629:   FAGetGlobalArray(fa,g,2,&ga);
    630:   for (k=y; k<y+n; k++) {
    631:     for (i=x; i<x+m; i++) {
    632:       ga[k][i].X = (R + k*Rscale)*PetscCosScalar(1.*PETSC_PI/6. + i*Ascale);
    633:       ga[k][i].Y = (R + k*Rscale)*PetscSinScalar(1.*PETSC_PI/6. + i*Ascale) - 4.*R;
    634:     }
    635:   }
    636:   FARestoreGlobalArray(fa,g,2,&ga);
    637:   return(0);
    638: }

    640: PetscErrorCode FAMapRegion2(FA fa,Vec g)
    641: {
    643:   PetscReal      R = 1.0,Rscale,Ascale;
    644:   PetscInt       i,k,x,y,m,n;
    645:   Field          **ga;

    648:   Rscale = R/(fa->r2-1);
    649:   Ascale = 2.0*PETSC_PI/fa->p2;

    651:   FAGetGlobalCorners(fa,1,&x,&y,&m,&n);
    652:   FAGetGlobalArray(fa,g,1,&ga);
    653:   for (k=y; k<y+n; k++) {
    654:     for (i=x; i<x+m; i++) {
    655:       ga[k][i].X = (R + k*Rscale)*PetscCosScalar(i*Ascale - PETSC_PI/2.0);
    656:       ga[k][i].Y = (R + k*Rscale)*PetscSinScalar(i*Ascale - PETSC_PI/2.0);
    657:     }
    658:   }
    659:   FARestoreGlobalArray(fa,g,1,&ga);
    660:   return(0);
    661: }

    663: PetscErrorCode FAMapRegion1(FA fa,Vec g)
    664: {
    666:   PetscReal      R = 1.0,Rscale,Ascale1,Ascale3;
    667:   PetscInt       i,k,x,y,m,n;
    668:   Field          **ga;

    671:   Rscale  = R/(fa->r1-1);
    672:   Ascale1 = 2.0*PETSC_PI/fa->p2;
    673:   Ascale3 = 2.0*PETSC_PI/(3.0*(fa->p1 - fa->p2 - 1));

    675:   FAGetGlobalCorners(fa,0,&x,&y,&m,&n);
    676:   FAGetGlobalArray(fa,g,0,&ga);

    678:   /* This mapping is WRONG! Not sure how to do it so I've done a poor job of
    679:      it. You can see that the grid connections are correct. */
    680:   for (k=y; k<y+n; k++) {
    681:     for (i=x; i<x+m; i++) {
    682:       if (i < (fa->p1-fa->p2)/2) {
    683:         ga[k][i].X = (2.0*R + k*Rscale)*PetscCosScalar(i*Ascale3);
    684:         ga[k][i].Y = (2.0*R + k*Rscale)*PetscSinScalar(i*Ascale3) - 4.*R;
    685:       } else if (i > fa->p2 + (fa->p1 - fa->p2)/2) {
    686:         ga[k][i].X = (2.0*R + k*Rscale)*PetscCosScalar(PETSC_PI+i*Ascale3);
    687:         ga[k][i].Y = (2.0*R + k*Rscale)*PetscSinScalar(PETSC_PI+i*Ascale3) - 4.*R;
    688:       } else {
    689:         ga[k][i].X = (2.*R + k*Rscale)*PetscCosScalar((i-(fa->p1-fa->p2)/2)*Ascale1 - PETSC_PI/2.0);
    690:         ga[k][i].Y = (2.*R + k*Rscale)*PetscSinScalar((i-(fa->p1-fa->p2)/2)*Ascale1 - PETSC_PI/2.0);
    691:       }
    692:     }
    693:   }
    694:   FARestoreGlobalArray(fa,g,0,&ga);
    695:   return(0);
    696: }

    698: /* Simple test to check that the ghost points are properly updated */
    699: PetscErrorCode FATest(FA fa)
    700: {
    702:   Vec            l,g;
    703:   Field          **la;
    704:   PetscInt       x,y,m,n,j,i,k,p;
    705:   PetscMPIInt    rank;

    708:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

    710:   FAGetGlobalVector(fa,&g);
    711:   FAGetLocalVector(fa,&l);

    713:   /* fill up global vector of one region at a time with ITS logical coordinates, then update LOCAL
    714:      vector; print local vectors to confirm they are correctly filled */
    715:   for (j=0; j<3; j++) {
    716:     VecSet(g,0.0);
    717:     FAGetGlobalCorners(fa,j,&x,&y,&m,&n);
    718:     PetscPrintf(PETSC_COMM_WORLD,"\nFilling global region %d, showing local results \n",j+1);
    719:     FAGetGlobalArray(fa,g,j,&la);
    720:     for (k=y; k<y+n; k++) {
    721:       for (i=x; i<x+m; i++) {
    722:         la[k][i].X = i;
    723:         la[k][i].Y = k;
    724:       }
    725:     }
    726:     FARestoreGlobalArray(fa,g,j,&la);

    728:     FAGlobalToLocal(fa,g,l);
    729:     DrawFA(fa,g);
    730:     DrawFA(fa,l);

    732:     for (p=0; p<3; p++) {
    733:       FAGetLocalCorners(fa,p,&x,&y,&m,&n);
    734:       FAGetLocalArray(fa,l,p,&la);
    735:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"\n[%d] Local array for region %d \n",rank,p+1);
    736:       for (k=y+n-1; k>=y; k--) { /* print in reverse order to match diagram in paper */
    737:         for (i=x; i<x+m; i++) {
    738:           PetscSynchronizedPrintf(PETSC_COMM_WORLD,"(%G,%G) ",la[k][i].X,la[k][i].Y);
    739:         }
    740:         PetscSynchronizedPrintf(PETSC_COMM_WORLD,"\n");
    741:       }
    742:       FARestoreLocalArray(fa,l,p,&la);
    743:       PetscSynchronizedFlush(PETSC_COMM_WORLD);
    744:     }
    745:   }
    746:   VecDestroy(&g);
    747:   VecDestroy(&l);
    748:   return(0);
    749: }

    753: int main(int argc,char **argv)
    754: {
    755:   FA             fa;
    757:   Vec            g,l;

    759:   PetscInitialize(&argc,&argv,0,help);

    761:   FACreate(&fa);
    762:   /* FATest(fa);*/

    764:   FAGetGlobalVector(fa,&g);
    765:   FAGetLocalVector(fa,&l);

    767:   FAMapRegion1(fa,g);
    768:   FAMapRegion2(fa,g);
    769:   FAMapRegion3(fa,g);

    771:   FAGlobalToLocal(fa,g,l);
    772:   DrawFA(fa,g);
    773:   DrawFA(fa,l);

    775:   VecDestroy(&g);
    776:   VecDestroy(&l);
    777:   FADestroy(&fa);

    779:   PetscFinalize();
    780:   return 0;
    781: }