Actual source code: ex2.c

      2: static char help[] = "Reads a a simple unstructured grid from a file. Partitions it,\n\
      3: and distributes the grid data accordingly\n\n";

      5: /*T
      6:    Concepts: Mat^partitioning a matrix;
      7:    Processors: n
      8: T*/

     10: /*
     11:     Updates of this example MAY be found at
     12:        http://www.mcs.anl.gov/petsc/src/dm/ao/examples/tutorials/ex2.c

     14:     This is a very basic, even crude, example of managing an unstructured
     15:     grid in parallel.

     17:     This particular code is for a Galerkin-style finite element method.

     19:     After the calls below, each processor will have
     20:       1) a list of elements it "owns"; for each "owned" element it will have the global
     21:          numbering of the three vertices; stored in gdata->ele;
     22:       2) a list of vertices it "owns". For each owned it will have the x and y
     23:          coordinates; stored in gdata->vert

     25:     It will not have
     26:       1) list of ghost elements (since they are not needed for traditional
     27:          Galerkin style finite element methods). For various finite volume methods
     28:          you may need the ghost element lists, these may be generated using the
     29:          element neighbor information given in the file database.

     31:     To compute the local element stiffness matrix and load vector, each processor
     32:     will need the vertex coordinates for all of the vertices on the locally
     33:     "owned" elements.  This data could be obtained by doing the appropriate vector
     34:     scatter on the data stored in gdata->vert; we haven't had time to demonstrate this.

     36:     Clearly writing a complete parallel unstructured grid code with PETSc is still
     37:     a good deal of work and requires a lot of application level coding.  BUT, at least
     38:     PETSc can manage all of the nonlinear and linear solvers (including matrix assembly
     39:     etc.), which allows the programmer to concentrate his or her efforts on managing
     40:     the unstructured grid. The PETSc team is developing additional library objects
     41:     to help manage parallel unstructured grid computations.  Unfortunately we have
     42:     not had time to complete these yet, so the application programmer still must
     43:     manage much of the parallel grid manipulation as indicated below.

     45: */

     47: /*
     48:   Include "petscmat.h" so that we can use matrices.
     49:   automatically includes:
     50:      petsc.h       - base PETSc routines   petscvec.h    - vectors
     51:      petscsys.h    - system routines       petscmat.h    - matrices
     52:      petscis.h     - index sets            petscviewer.h - viewers

     54:   Include "petscao.h" allows use of the AO (application ordering) commands,
     55:   used below for renumbering the vertex numbers after the partitioning.

     57:   Include "petscbt.h" for managing logical bit arrays that are used to
     58:   conserve space. Note that the code does use order N bit arrays on each
     59:   processor so is theoretically not scalable, but even with 64 million
     60:   vertices it will only need temporarily 8 megabytes of memory for the
     61:   bit array so one can still do very large problems with this approach,
     62:   since the bit arrays are freed before the vectors and matrices are
     63:   created.
     64: */
     65:  #include petscmat.h
     66:  #include petscao.h
     67:  #include petscbt.h

     69: /*
     70:     This is the user-defined grid data context
     71: */
     72: typedef struct {
     73:   PetscInt    n_vert,n_ele;
     74:   PetscInt    mlocal_vert,mlocal_ele;
     75:   PetscInt    *ele;
     76:   PetscScalar *vert;
     77:   PetscInt    *ia,*ja;
     78:   IS     isnewproc;
     79:   PetscInt    *localvert,nlocal; /* used to stash temporarily old global vertex number of new vertex */
     80: } GridData;

     82: /*
     83:    Variables on all processors:
     84:       n_vert          - total number of vertices
     85:       mlocal_vert     - number of vertices on this processor
     86:       vert            - x,y coordinates of local vertices

     88:       n_ele           - total number of elements
     89:       mlocal_ele      - number of elements on this processor
     90:       ele             - vertices of elements on this processor

     92:       ia, ja          - adjacency graph of elements (for partitioning)
     93:
     94:    Variables on processor 0 during data reading from file:
     95:       mmlocal_vert[i] - number of vertices on each processor
     96:       tmpvert         - x,y coordinates of vertices on any processor (as read in)

     98:       mmlocal_ele[i]  - number of elements on each processor

    100:       tmpia, tmpja    - adjacency graph of elements for other processors

    102:    Notes:
    103:    The code below has a great deal of IO (print statements). This is to allow one to track
    104:    the renumbering and movement of data among processors. In an actual
    105:    production run, IO of this type would be deactivated.

    107:    To use the ParMETIS partitioner run with the option -mat_partitioning_type parmetis
    108:    otherwise it defaults to the initial element partitioning induced when the data
    109:    is read in.

    111:    To understand the parallel performance of this type of code, it is important
    112:    to profile the time spent in various events in the code; running with the
    113:    option -log_summary will indicate how much time is spent in the routines
    114:    below.   Of course, for very small problems, such as the sample grid used here,
    115:    the profiling results are meaningless.
    116: */

    127: int main(int argc,char **args)
    128: {
    130:   PetscEvent     READ_EVENT,PARTITION_ELEMENT_EVENT,MOVE_ELEMENT_EVENT;
    131:   PetscEvent     PARTITION_VERTEX_EVENT,MOVE_VERTEX_EVENT;
    132:   GridData       gdata;

    134:   PetscInitialize(&argc,&args,(char *)0,help);

    136:   PetscLogEventRegister(&READ_EVENT,             "Read Data",0);
    137:   PetscLogEventRegister(&PARTITION_ELEMENT_EVENT,"Partition elemen",0);
    138:   PetscLogEventRegister(&MOVE_ELEMENT_EVENT,     "Move elements",0);
    139:   PetscLogEventRegister(&PARTITION_VERTEX_EVENT, "Partition vertic",0);
    140:   PetscLogEventRegister(&MOVE_VERTEX_EVENT,      "Move vertices",0);

    142:   PetscLogEventBegin(READ_EVENT,0,0,0,0);
    143:   DataRead(&gdata);
    144:   PetscLogEventEnd(READ_EVENT,0,0,0,0);
    145:   PetscLogEventBegin(PARTITION_ELEMENT_EVENT,0,0,0,0);
    146:   DataPartitionElements(&gdata);
    147:   PetscLogEventEnd(PARTITION_ELEMENT_EVENT,0,0,0,0);
    148:   PetscLogEventBegin(MOVE_ELEMENT_EVENT,0,0,0,0);
    149:   DataMoveElements(&gdata);
    150:   PetscLogEventEnd(MOVE_ELEMENT_EVENT,0,0,0,0);
    151:   PetscLogEventBegin(PARTITION_VERTEX_EVENT,0,0,0,0);
    152:   DataPartitionVertices(&gdata);
    153:   PetscLogEventEnd(PARTITION_VERTEX_EVENT,0,0,0,0);
    154:   PetscLogEventBegin(MOVE_VERTEX_EVENT,0,0,0,0);
    155:   DataMoveVertices(&gdata);
    156:   PetscLogEventEnd(MOVE_VERTEX_EVENT,0,0,0,0);
    157:   DataDestroy(&gdata);

    159:   PetscFinalize();
    160:   return 0;
    161: }

    166: /*
    167:      Reads in the grid data from a file; each processor is naively
    168:   assigned a continuous chunk of vertex and element data. Later the data
    169:   will be partitioned and moved to the appropriate processor.
    170: */
    171: PetscErrorCode DataRead(GridData *gdata)
    172: {
    173:   PetscMPIInt    rank,size;
    174:   PetscInt       n_vert,*mmlocal_vert,mlocal_vert,i,*ia,*ja,cnt,j;
    175:   PetscInt       mlocal_ele,*mmlocal_ele,*ele,*tmpele,n_ele,net,a1,a2,a3;
    176:   PetscInt       *iatmp,*jatmp;
    178:   char           msg[128];
    179:   PetscScalar    *vert,*tmpvert;
    180:   MPI_Status     status;

    183:   /*
    184:      Processor 0 opens the file, reads in chunks of data and sends a portion off to
    185:    each other processor.

    187:      Note: For a truely scalable IO portion of the code, one would store
    188:    the grid data in a binary file and use MPI-IO commands to have each
    189:    processor read in the parts that it needs. However in most circumstances
    190:    involving up to a say a million nodes and 100 processors this approach
    191:    here is fine.
    192:   */
    193:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
    194:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

    196:   if (!rank) {
    197:     FILE *fd;
    198:     fd = fopen("usgdata","r"); if (!fd) SETERRQ(1,"Cannot open grid file");

    200:     /* read in number of vertices */
    201:     fgets(msg,128,fd);
    202:     printf("File msg:%s",msg);
    203:     fscanf(fd,"Number Vertices = %d\n",&n_vert);
    204:     printf("Number of grid vertices %d\n",n_vert);

    206:     /* broadcast number of vertices to all processors */
    207:     MPI_Bcast(&n_vert,1,MPI_INT,0,PETSC_COMM_WORLD);
    208:     mlocal_vert  = n_vert/size + ((n_vert % size) > 0);

    210:     /*
    211:       allocate enough room for the first processor to keep track of how many
    212:       vertices are assigned to each processor. Splitting vertices equally amoung
    213:       all processors.
    214:     */
    215:     PetscMalloc(size*sizeof(PetscInt),&mmlocal_vert);
    216:     for (i=0; i<size; i++) {
    217:       mmlocal_vert[i] = n_vert/size + ((n_vert % size) > i);
    218:       printf("Processor %d assigned %d vertices\n",i,mmlocal_vert[i]);
    219:     }

    221:     /*
    222:        Read in vertices assigned to first processor
    223:     */
    224:     PetscMalloc(2*mmlocal_vert[0]*sizeof(PetscScalar),&vert);
    225:     printf("Vertices assigned to processor 0\n");
    226:     for (i=0; i<mlocal_vert; i++) {
    227:       fscanf(fd,"%d %lf %lf\n",&cnt,vert+2*i,vert+2*i+1);
    228:       printf("%d %g %g\n",cnt,PetscRealPart(vert[2*i]),PetscRealPart(vert[2*i+1]));
    229:     }

    231:     /*
    232:        Read in vertices for all the other processors
    233:     */
    234:     PetscMalloc(2*mmlocal_vert[0]*sizeof(PetscScalar),&tmpvert);
    235:     for (j=1; j<size; j++) {
    236:       printf("Vertices assigned to processor %d\n",j);
    237:       for (i=0; i<mmlocal_vert[j]; i++) {
    238:         fscanf(fd,"%d %lf %lf\n",&cnt,tmpvert+2*i,tmpvert+2*i+1);
    239:         printf("%d %g %g\n",cnt,tmpvert[2*i],tmpvert[2*i+1]);
    240:       }
    241:       MPI_Send(tmpvert,2*mmlocal_vert[j],MPIU_SCALAR,j,0,PETSC_COMM_WORLD);
    242:     }
    243:     PetscFree(tmpvert);
    244:     PetscFree(mmlocal_vert);

    246:     fscanf(fd,"Number Elements = %d\n",&n_ele);
    247:     printf("Number of grid elements %d\n",n_ele);

    249:     /*
    250:        Broadcast number of elements to all processors
    251:     */
    252:     MPI_Bcast(&n_ele,1,MPI_INT,0,PETSC_COMM_WORLD);
    253:     mlocal_ele  = n_ele/size + ((n_ele % size) > 0);

    255:     /*
    256:       Allocate enough room for the first processor to keep track of how many
    257:       elements are assigned to each processor.
    258:     */
    259:     PetscMalloc(size*sizeof(PetscInt),&mmlocal_ele);
    260:     for (i=0; i<size; i++) {
    261:       mmlocal_ele[i] = n_ele/size + ((n_ele % size) > i);
    262:       printf("Processor %d assigned %d elements\n",i,mmlocal_ele[i]);
    263:     }
    264:
    265:     /*
    266:         read in element information for the first processor
    267:     */
    268:     PetscMalloc(3*mmlocal_ele[0]*sizeof(PetscInt),&ele);
    269:     printf("Elements assigned to processor 0\n");
    270:     for (i=0; i<mlocal_ele; i++) {
    271:       fscanf(fd,"%d %d %d %d\n",&cnt,ele+3*i,ele+3*i+1,ele+3*i+2);
    272:       printf("%d %d %d %d\n",cnt,ele[3*i],ele[3*i+1],ele[3*i+2]);
    273:     }

    275:     /*
    276:        Read in elements for all the other processors
    277:     */
    278:     PetscMalloc(3*mmlocal_ele[0]*sizeof(PetscInt),&tmpele);
    279:     for (j=1; j<size; j++) {
    280:       printf("Elements assigned to processor %d\n",j);
    281:       for (i=0; i<mmlocal_ele[j]; i++) {
    282:         fscanf(fd,"%d %d %d %d\n",&cnt,tmpele+3*i,tmpele+3*i+1,tmpele+3*i+2);
    283:         printf("%d %d %d %d\n",cnt,tmpele[3*i],tmpele[3*i+1],tmpele[3*i+2]);
    284:       }
    285:       MPI_Send(tmpele,3*mmlocal_ele[j],MPI_INT,j,0,PETSC_COMM_WORLD);
    286:     }
    287:     PetscFree(tmpele);

    289:     /*
    290:          Read in element neighbors for processor 0
    291:          We don't know how many spaces in ja[] to allocate so we allocate
    292:        3*the number of local elements, this is the maximum it could be
    293:     */
    294:     PetscMalloc((mlocal_ele+1)*sizeof(PetscInt),&ia);
    295:     PetscMalloc((3*mlocal_ele+1)*sizeof(PetscInt),&ja);
    296:     net   = 0;
    297:     ia[0] = 0;
    298:     printf("Element neighbors on processor 0\n");
    299:     fgets(msg,128,fd);
    300:     for (i=0; i<mlocal_ele; i++) {
    301:       fscanf(fd,"%d %d %d %d\n",&cnt,&a1,&a2,&a3);
    302:       printf("%d %d %d %d\n",cnt,a1,a2,a3);
    303:       if (a1 >= 0) {ja[net++] = a1;}
    304:       if (a2 >= 0) {ja[net++] = a2;}
    305:       if (a3 >= 0) {ja[net++] = a3;}
    306:       ia[i+1] = net;
    307:     }

    309:     printf("ia values for processor 0\n");
    310:     for (i=0; i<mlocal_ele+1; i++) {
    311:       printf("%d ",ia[i]);
    312:     }
    313:     printf("\n");
    314:     printf("ja values for processor 0\n");
    315:     for (i=0; i<ia[mlocal_ele]; i++) {
    316:       printf("%d ",ja[i]);
    317:     }
    318:     printf("\n");

    320:     /*
    321:        Read in element neighbor information for all other processors
    322:     */
    323:     PetscMalloc((mlocal_ele+1)*sizeof(PetscInt),&iatmp);
    324:     PetscMalloc((3*mlocal_ele+1)*sizeof(PetscInt),&jatmp);
    325:     for (j=1; j<size; j++) {
    326:       net   = 0;
    327:       iatmp[0] = 0;
    328:       printf("Element neighbors on processor %d\n",j);
    329:       for (i=0; i<mmlocal_ele[j]; i++) {
    330:         fscanf(fd,"%d %d %d %d\n",&cnt,&a1,&a2,&a3);
    331:         printf("%d %d %d %d\n",cnt,a1,a2,a3);
    332:         if (a1 >= 0) {jatmp[net++] = a1;}
    333:         if (a2 >= 0) {jatmp[net++] = a2;}
    334:         if (a3 >= 0) {jatmp[net++] = a3;}
    335:         iatmp[i+1] = net;
    336:       }

    338:       printf("ia values for processor %d\n",j);
    339:       for (i=0; i<mmlocal_ele[j]+1; i++) {
    340:         printf("%d ",iatmp[i]);
    341:       }
    342:       printf("\n");
    343:       printf("ja values for processor %d\n",j);
    344:       for (i=0; i<iatmp[mmlocal_ele[j]]; i++) {
    345:         printf("%d ",jatmp[i]);
    346:       }
    347:       printf("\n");

    349:       /* send graph off to appropriate processor */
    350:       MPI_Send(iatmp,mmlocal_ele[j]+1,MPI_INT,j,0,PETSC_COMM_WORLD);
    351:       MPI_Send(jatmp,iatmp[mmlocal_ele[j]],MPI_INT,j,0,PETSC_COMM_WORLD);
    352:     }
    353:     PetscFree(iatmp);
    354:     PetscFree(jatmp);
    355:     PetscFree(mmlocal_ele);

    357:     fclose(fd);
    358:   } else {
    359:     /*
    360:         We are not the zeroth processor so we do not open the file
    361:       rather we wait for processor 0 to send us our data.
    362:     */

    364:     /* receive total number of vertices */
    365:     MPI_Bcast(&n_vert,1,MPI_INT,0,PETSC_COMM_WORLD);
    366:     mlocal_vert = n_vert/size + ((n_vert % size) > rank);

    368:     /* receive vertices */
    369:     PetscMalloc(2*(mlocal_vert+1)*sizeof(PetscScalar),&vert);
    370:     MPI_Recv(vert,2*mlocal_vert,MPIU_SCALAR,0,0,PETSC_COMM_WORLD,&status);

    372:     /* receive total number of elements */
    373:     MPI_Bcast(&n_ele,1,MPI_INT,0,PETSC_COMM_WORLD);
    374:     mlocal_ele = n_ele/size + ((n_ele % size) > rank);

    376:     /* receive elements */
    377:     PetscMalloc(3*(mlocal_ele+1)*sizeof(PetscInt),&ele);
    378:     MPI_Recv(ele,3*mlocal_ele,MPI_INT,0,0,PETSC_COMM_WORLD,&status);

    380:     /* receive element adjacency graph */
    381:     PetscMalloc((mlocal_ele+1)*sizeof(PetscInt),&ia);
    382:     MPI_Recv(ia,mlocal_ele+1,MPI_INT,0,0,PETSC_COMM_WORLD,&status);

    384:     PetscMalloc((ia[mlocal_ele]+1)*sizeof(PetscInt),&ja);
    385:     MPI_Recv(ja,ia[mlocal_ele],MPI_INT,0,0,PETSC_COMM_WORLD,&status);
    386:   }

    388:   gdata->n_vert      = n_vert;
    389:   gdata->n_ele       = n_ele;
    390:   gdata->mlocal_vert = mlocal_vert;
    391:   gdata->mlocal_ele  = mlocal_ele;
    392:   gdata->ele         = ele;
    393:   gdata->vert        = vert;

    395:   gdata->ia          = ia;
    396:   gdata->ja          = ja;

    398:   return(0);
    399: }

    404: /*
    405:          Given the grid data spread across the processors, determines a
    406:    new partitioning of the CELLS (elements) to reduce the number of cut edges between
    407:    cells (elements).
    408: */
    409: PetscErrorCode DataPartitionElements(GridData *gdata)
    410: {
    411:   Mat             Adj;                /* adjacency matrix */
    412:   PetscInt        *ia,*ja;
    413:   PetscInt        mlocal_ele,n_ele;
    414:   PetscErrorCode  ierr;
    415:   MatPartitioning part;
    416:   IS              isnewproc;

    419:   n_ele       = gdata->n_ele;
    420:   mlocal_ele  = gdata->mlocal_ele;

    422:   ia          = gdata->ia;
    423:   ja          = gdata->ja;

    425:   /*
    426:       Create the adjacency graph matrix
    427:   */
    428:   MatCreateMPIAdj(PETSC_COMM_WORLD,mlocal_ele,n_ele,ia,ja,PETSC_NULL,&Adj);

    430:   /*
    431:       Create the partioning object
    432:   */
    433:   MatPartitioningCreate(PETSC_COMM_WORLD,&part);
    434:   MatPartitioningSetAdjacency(part,Adj);
    435:   MatPartitioningSetFromOptions(part);
    436:   MatPartitioningApply(part,&isnewproc);
    437:   MatPartitioningDestroy(part);

    439:   /*
    440:        isnewproc - indicates for each local element the new processor it is assigned to
    441:   */
    442:   PetscPrintf(PETSC_COMM_WORLD,"New processor assignment for each element\n");
    443:   ISView(isnewproc,PETSC_VIEWER_STDOUT_WORLD);
    444:   gdata->isnewproc = isnewproc;

    446:   /*
    447:       Free the adjacency graph data structures
    448:   */
    449:   MatDestroy(Adj);

    452:   return(0);
    453: }

    457: /*
    458:       Moves the grid element data to be on the correct processor for the new
    459:    element partitioning.
    460: */
    461: PetscErrorCode DataMoveElements(GridData *gdata)
    462: {
    464:   PetscInt       *counts,i,*idx;
    465:   PetscMPIInt    rank,size;
    466:   Vec            vele,veleold;
    467:   PetscScalar    *array;
    468:   IS             isscat,isnum;
    469:   VecScatter     vecscat;

    473:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
    474:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

    476:   /*
    477:       Determine how many elements are assigned to each processor
    478:   */
    479:   PetscMalloc(size*sizeof(PetscInt),&counts);
    480:   ISPartitioningCount(gdata->isnewproc,counts);

    482:   /*
    483:      Create a vector to contain the newly ordered element information

    485:      Note: we use vectors to communicate this data since we can use the powerful
    486:      VecScatter routines to get the data to the correct location. This is a little
    487:      wasteful since the vectors hold double precision numbers instead of integers,
    488:      but since this is just a setup phase in the entire numerical computation that
    489:      is only called once it is not a measureable performance bottleneck.
    490:   */
    491:   VecCreate(PETSC_COMM_WORLD,&vele);
    492:   VecSetSizes(vele,3*counts[rank],PETSC_DECIDE);
    493:   VecSetFromOptions(vele);

    495:   /*
    496:       Create an index set from the isnewproc index set to indicate the mapping TO
    497:   */
    498:   ISPartitioningToNumbering(gdata->isnewproc,&isnum);
    499:   ISDestroy(gdata->isnewproc);
    500:   /*
    501:       There are three data items per cell (element), the integer vertex numbers of its three
    502:     coordinates (we convert to double to use the scatter) (one can think
    503:     of the vectors of having a block size of 3, then there is one index in idx[] for each element)
    504:   */
    505:   ISGetIndices(isnum,&idx);
    506:   for (i=0; i<gdata->mlocal_ele; i++) {
    507:     idx[i] *= 3;
    508:   }
    509:   ISCreateBlock(PETSC_COMM_WORLD,3,gdata->mlocal_ele,idx,&isscat);
    510:   ISRestoreIndices(isnum,&idx);
    511:   ISDestroy(isnum);

    513:   /*
    514:      Create a vector to contain the original vertex information for each element
    515:   */
    516:   VecCreateSeq(PETSC_COMM_SELF,3*gdata->mlocal_ele,&veleold);
    517:   VecGetArray(veleold,&array);
    518:   for (i=0; i<3*gdata->mlocal_ele; i++) {
    519:     array[i] = gdata->ele[i];
    520:   }
    521:   VecRestoreArray(veleold,&array);
    522:
    523:   /*
    524:      Scatter the element vertex information (still in the original vertex ordering) to the correct processor
    525:   */
    526:   VecScatterCreate(veleold,PETSC_NULL,vele,isscat,&vecscat);
    527:   ISDestroy(isscat);
    528:   VecScatterBegin(veleold,vele,INSERT_VALUES,SCATTER_FORWARD,vecscat);
    529:   VecScatterEnd(veleold,vele,INSERT_VALUES,SCATTER_FORWARD,vecscat);
    530:   VecScatterDestroy(vecscat);
    531:   VecDestroy(veleold);

    533:   /*
    534:      Put the element vertex data into a new allocation of the gdata->ele
    535:   */
    536:   PetscFree(gdata->ele);
    537:   gdata->mlocal_ele = counts[rank];
    538:   PetscFree(counts);
    539:   PetscMalloc(3*gdata->mlocal_ele*sizeof(PetscInt),&gdata->ele);
    540:   VecGetArray(vele,&array);
    541:   for (i=0; i<3*gdata->mlocal_ele; i++) {
    542:     gdata->ele[i] = (int)PetscRealPart(array[i]);
    543:   }
    544:   VecRestoreArray(vele,&array);
    545:   VecDestroy(vele);

    547:   PetscPrintf(PETSC_COMM_WORLD,"Old vertex numbering in new element ordering\n");
    548:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor %d\n",rank);
    549:   for (i=0; i<gdata->mlocal_ele; i++) {
    550:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%d %d %d %d\n",i,gdata->ele[3*i],gdata->ele[3*i+1],
    551:                             gdata->ele[3*i+2]);
    552:   }
    553:   PetscSynchronizedFlush(PETSC_COMM_WORLD);

    555:   return(0);
    556: }

    560: /*
    561:          Given the newly partitioned cells (elements), this routine partitions the
    562:      vertices.

    564:      The code is not completely scalable since it requires
    565:      1) O(n_vert) bits per processor memory
    566:      2) uses O(size) stages of communication; each of size O(n_vert) bits
    567:      3) it is sequential (each processor marks it vertices ONLY after all processors
    568:         to the left have marked theirs.
    569:      4) the amount of work on the last processor is O(n_vert)

    571:      The algorithm also does not take advantage of vertices that are "interior" to a
    572:      processors elements (that is; is not contained in any element on another processor).
    573:      A better algorithm would first have all processors determine "interior" vertices and
    574:      make sure they are retained on that processor before listing "boundary" vertices.

    576:      The algorithm is:
    577:      a) each processor waits for a message from the left containing mask of all marked vertices
    578:      b) it loops over all local elements, generating a list of vertices it will
    579:         claim (not claiming ones that have already been marked in the bit-array)
    580:         it claims at most n_vert/size vertices
    581:      c) it sends to the right the mask

    583:      This is a quick-and-dirty implementation; it should work fine for many problems,
    584:      but will need to be replaced once profiling shows that it takes a large amount of
    585:      time. An advantage is it requires no searching or sorting.
    586:
    587: */
    588: PetscErrorCode DataPartitionVertices(GridData *gdata)
    589: {
    590:   PetscInt       n_vert = gdata->n_vert,*localvert;
    591:   PetscInt       mlocal_ele = gdata->mlocal_ele,*ele = gdata->ele,i,j,nlocal = 0,nmax;
    592:   PetscMPIInt    rank,size;
    594:   PetscBT        mask;
    595:   MPI_Status     status;

    598:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    599:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

    601:   /*
    602:       Allocated space to store bit-array indicting vertices marked
    603:   */
    604:   PetscBTCreate(n_vert,mask);

    606:   /*
    607:      All processors except last can have a maximum of n_vert/size vertices assigned
    608:      (because last processor needs to handle any leftovers)
    609:   */
    610:   nmax = n_vert/size;
    611:   if (rank == size-1) {
    612:     nmax = n_vert;
    613:   }

    615:   /*
    616:      Receive list of marked vertices from left
    617:   */
    618:   if (rank) {
    619:     MPI_Recv(mask,PetscBTLength(n_vert),MPI_CHAR,rank-1,0,PETSC_COMM_WORLD,&status);
    620:   }

    622:   if (rank == size-1) {
    623:     /* last processor gets all the rest */
    624:     for (i=0; i<n_vert; i++) {
    625:       if (!PetscBTLookup(mask,i)) {
    626:         nlocal++;
    627:       }
    628:     }
    629:     nmax = nlocal;
    630:   }

    632:   /*
    633:      Now we know how many are local, allocated enough space for them and mark them
    634:   */
    635:   PetscMalloc((nmax+1)*sizeof(PetscInt),&localvert);

    637:   /* generate local list and fill in mask */
    638:   nlocal = 0;
    639:   if (rank < size-1) {
    640:     /* count my vertices */
    641:     for (i=0; i<mlocal_ele; i++) {
    642:       for (j=0; j<3; j++) {
    643:         if (!PetscBTLookupSet(mask,ele[3*i+j])) {
    644:           localvert[nlocal++] = ele[3*i+j];
    645:           if (nlocal >= nmax) goto foundenough2;
    646:         }
    647:       }
    648:     }
    649:     foundenough2:;
    650:   } else {
    651:     /* last processor gets all the rest */
    652:     for (i=0; i<n_vert; i++) {
    653:       if (!PetscBTLookup(mask,i)) {
    654:         localvert[nlocal++] = i;
    655:       }
    656:     }
    657:   }
    658:   /*
    659:       Send bit mask on to next processor
    660:   */
    661:   if (rank < size-1) {
    662:     MPI_Send(mask,PetscBTLength(n_vert),MPI_CHAR,rank+1,0,PETSC_COMM_WORLD);
    663:   }
    664:   PetscBTDestroy(mask);

    666:   gdata->localvert = localvert;
    667:   gdata->nlocal    = nlocal;

    669:   /* print lists of owned vertices */
    670:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Number vertices assigned %d\n",rank,nlocal);
    671:   PetscSynchronizedFlush(PETSC_COMM_WORLD);
    672:   PetscIntView(nlocal,localvert,PETSC_VIEWER_STDOUT_WORLD);

    674:   return(0);
    675: }

    679: /*
    680:      Given the partitioning of the vertices; renumbers the element vertex lists for the
    681:      new vertex numbering and moves the vertex coordinate values to the correct processor
    682: */
    683: PetscErrorCode DataMoveVertices(GridData *gdata)
    684: {
    685:   AO             ao;
    687:   PetscMPIInt    rank;
    688:   PetscInt       i;
    689:   Vec            vert,overt;
    690:   VecScatter     vecscat;
    691:   IS             isscat;
    692:   PetscScalar    *avert;

    695:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

    697:   /* ---------------------------------------------------------------------
    698:       Create a global reodering of the vertex numbers
    699:   */
    700:   AOCreateBasic(PETSC_COMM_WORLD,gdata->nlocal,gdata->localvert,PETSC_NULL,&ao);

    702:   /*
    703:      Change the element vertex information to the new vertex numbering
    704:   */
    705:   AOApplicationToPetsc(ao,3*gdata->mlocal_ele,gdata->ele);
    706:   PetscPrintf(PETSC_COMM_WORLD,"New vertex numbering in new element ordering\n");
    707:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor %d\n",rank);
    708:   for (i=0; i<gdata->mlocal_ele; i++) {
    709:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%d %d %d %d\n",i,gdata->ele[3*i],gdata->ele[3*i+1],
    710:                             gdata->ele[3*i+2]);
    711:   }
    712:   PetscSynchronizedFlush(PETSC_COMM_WORLD);

    714:   /*
    715:      Destroy the AO that is no longer needed
    716:   */
    717:   AODestroy(ao);

    719:   /* --------------------------------------------------------------------
    720:       Finally ship the vertex coordinate information to its owning process
    721:       note, we do this in a way very similar to what was done for the element info
    722:   */
    723:   /* create a vector to contain the newly ordered vertex information */
    724:   VecCreateSeq(PETSC_COMM_SELF,2*gdata->nlocal,&vert);

    726:   /* create a vector to contain the old ordered vertex information */
    727:   VecCreateMPIWithArray(PETSC_COMM_WORLD,2*gdata->mlocal_vert,PETSC_DECIDE,gdata->vert,
    728:                                &overt);

    730:   /*
    731:       There are two data items per vertex, the x and y coordinates (i.e. one can think
    732:     of the vectors of having a block size of 2 and there is one index in localvert[] for each block)
    733:   */
    734:   for (i=0; i<gdata->nlocal; i++) gdata->localvert[i] *= 2;
    735:   ISCreateBlock(PETSC_COMM_WORLD,2,gdata->nlocal,gdata->localvert,&isscat);
    736:   PetscFree(gdata->localvert);

    738:   /*
    739:       Scatter the element vertex information to the correct processor
    740:   */
    741:   VecScatterCreate(overt,isscat,vert,PETSC_NULL,&vecscat);
    742:   ISDestroy(isscat);
    743:   VecScatterBegin(overt,vert,INSERT_VALUES,SCATTER_FORWARD,vecscat);
    744:   VecScatterEnd(overt,vert,INSERT_VALUES,SCATTER_FORWARD,vecscat);
    745:   VecScatterDestroy(vecscat);

    747:   VecDestroy(overt);
    748:   PetscFree(gdata->vert);
    749:
    750:   /*
    751:         Put resulting vertex information into gdata->vert array
    752:   */
    753:   PetscMalloc(2*gdata->nlocal*sizeof(PetscScalar),&gdata->vert);
    754:   VecGetArray(vert,&avert);
    755:   PetscMemcpy(gdata->vert,avert,2*gdata->nlocal*sizeof(PetscScalar));
    756:   VecRestoreArray(vert,&avert);
    757:   gdata->mlocal_vert = gdata->nlocal;
    758:   VecDestroy(vert);

    760:   PetscPrintf(PETSC_COMM_WORLD,"Vertex coordinates in new numbering\n");
    761:   for (i=0; i<gdata->mlocal_vert; i++) {
    762:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"(%g,%g)\n",gdata->vert[2*i],gdata->vert[2*i+1]);
    763:   }
    764:   PetscSynchronizedFlush(PETSC_COMM_WORLD);

    766:   return(0);
    767: }

    772: PetscErrorCode DataDestroy(GridData *gdata)
    773: {

    777:   PetscFree(gdata->ele);
    778:   PetscFree(gdata->vert);
    779:   return(0);
    780: }