Actual source code: pipes1.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      1: static char help[] = "This example demonstrates the use of DMNetwork \n\\n";
      2: /*
      3:   Example: mpiexec -n <np> ./pipes1 -ts_max_steps 10
      4: */

      6: #include "wash.h"
      7:  #include <petscdmplex.h>

      9: /*
     10:   WashNetworkDistribute - proc[0] distributes sequential wash object
     11:    Input Parameters:
     12: .  comm - MPI communicator
     13: .  wash - wash context with all network data in proc[0]

     15:    Output Parameter:
     16: .  wash - wash context with nedge, nvertex and edgelist distributed
     17: */
     18: PetscErrorCode WashNetworkDistribute(MPI_Comm comm,Wash wash)
     19: {
     21:   PetscMPIInt    rank,size,tag=0;
     22:   PetscInt       i,e,v,numEdges,numVertices,nedges,*eowners=NULL,estart,eend,*vtype=NULL,nvertices;
     23:   PetscInt       *edgelist = wash->edgelist,*nvtx=NULL,*vtxDone=NULL;

     26:   MPI_Comm_rank(comm,&rank);
     27:   MPI_Comm_size(comm,&size);

     29:   numEdges    = wash->nedge;
     30:   numVertices = wash->nvertex;

     32:   /* (1) all processes get global and local number of edges */
     33:   MPI_Bcast(&numEdges,1,MPIU_INT,0,comm);
     34:   nedges = numEdges/size; /* local nedges */
     35:   if (!rank) {
     36:     nedges += numEdges - size*(numEdges/size);
     37:   }
     38:   wash->Nedge = numEdges;
     39:   wash->nedge = nedges;
     40:   /* PetscPrintf(PETSC_COMM_SELF,"[%d] nedges %d, numEdges %d\n",rank,nedges,numEdges); */

     42:   PetscCalloc3(size+1,&eowners,size,&nvtx,numVertices,&vtxDone);
     43:   MPI_Allgather(&nedges,1,MPIU_INT,eowners+1,1,MPIU_INT,PETSC_COMM_WORLD);
     44:   eowners[0] = 0;
     45:   for (i=2; i<=size; i++) {
     46:     eowners[i] += eowners[i-1];
     47:   }

     49:   estart = eowners[rank];
     50:   eend   = eowners[rank+1];
     51:   /* PetscPrintf(PETSC_COMM_SELF,"[%d] own lists row %d - %d\n",rank,estart,eend); */

     53:   /* (2) distribute row block edgelist to all processors */
     54:   if (!rank) {
     55:     vtype = wash->vtype;
     56:     for (i=1; i<size; i++) {
     57:       /* proc[0] sends edgelist to proc[i] */
     58:       MPI_Send(edgelist+2*eowners[i],2*(eowners[i+1]-eowners[i]),MPIU_INT,i,tag,comm);

     60:       /* proc[0] sends vtype to proc[i] */
     61:       MPI_Send(vtype+2*eowners[i],2*(eowners[i+1]-eowners[i]),MPIU_INT,i,tag,comm);
     62:     }
     63:   } else {
     64:     MPI_Status      status;
     65:     PetscMalloc1(2*(eend-estart),&vtype);
     66:     PetscMalloc1(2*(eend-estart),&edgelist);

     68:     MPI_Recv(edgelist,2*(eend-estart),MPIU_INT,0,tag,comm,&status);
     69:     MPI_Recv(vtype,2*(eend-estart),MPIU_INT,0,tag,comm,&status);
     70:   }

     72:   wash->edgelist = edgelist;

     74:   /* (3) all processes get global and local number of vertices, without ghost vertices */
     75:   if (!rank) {
     76:     for (i=0; i<size; i++) {
     77:       for (e=eowners[i]; e<eowners[i+1]; e++) {
     78:         v = edgelist[2*e];
     79:         if (!vtxDone[v]) {
     80:           nvtx[i]++; vtxDone[v] = 1;
     81:         }
     82:         v = edgelist[2*e+1];
     83:         if (!vtxDone[v]) {
     84:           nvtx[i]++; vtxDone[v] = 1;
     85:         }
     86:       }
     87:     }
     88:   }
     89:   MPI_Bcast(&numVertices,1,MPIU_INT,0,PETSC_COMM_WORLD);
     90:   MPI_Scatter(nvtx,1,MPIU_INT,&nvertices,1,MPIU_INT,0,PETSC_COMM_WORLD);
     91:   PetscFree3(eowners,nvtx,vtxDone);

     93:   wash->Nvertex = numVertices;
     94:   wash->nvertex = nvertices;
     95:   wash->vtype   = vtype;
     96:   return(0);
     97: }

     99: PetscErrorCode WASHIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void* ctx)
    100: {
    102:   Wash           wash=(Wash)ctx;
    103:   DM             networkdm;
    104:   Vec            localX,localXdot,localF, localXold;
    105:   const PetscInt *cone;
    106:   PetscInt       vfrom,vto,offsetfrom,offsetto,varoffset;
    107:   PetscInt       v,vStart,vEnd,e,eStart,eEnd;
    108:   PetscInt       nend,type;
    109:   PetscBool      ghost;
    110:   PetscScalar    *farr,*juncf, *pipef;
    111:   PetscReal      dt;
    112:   Pipe           pipe;
    113:   PipeField      *pipex,*pipexdot,*juncx;
    114:   Junction       junction;
    115:   DMDALocalInfo  info;
    116:   const PetscScalar *xarr,*xdotarr, *xoldarr;

    119:   localX    = wash->localX;
    120:   localXdot = wash->localXdot;

    122:   TSGetSolution(ts,&localXold);
    123:   TSGetDM(ts,&networkdm);
    124:   TSGetTimeStep(ts,&dt);
    125:   DMGetLocalVector(networkdm,&localF);

    127:   /* Set F and localF as zero */
    128:   VecSet(F,0.0);
    129:   VecSet(localF,0.0);

    131:   /* Update ghost values of locaX and locaXdot */
    132:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
    133:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

    135:   DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot);
    136:   DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot);

    138:   VecGetArrayRead(localX,&xarr);
    139:   VecGetArrayRead(localXdot,&xdotarr);
    140:   VecGetArrayRead(localXold,&xoldarr);
    141:   VecGetArray(localF,&farr);

    143:    /* junction->type == JUNCTION:
    144:            juncf[0] = -qJ + sum(qin); juncf[1] = qJ - sum(qout)
    145:        junction->type == RESERVOIR (upper stream):
    146:            juncf[0] = -hJ + H0; juncf[1] = qJ - sum(qout)
    147:        junction->type == VALVE (down stream):
    148:            juncf[0] =  -qJ + sum(qin); juncf[1] = qJ
    149:   */
    150:   /* Vertex/junction initialization */
    151:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
    152:   for (v=vStart; v<vEnd; v++) {
    153:     DMNetworkIsGhostVertex(networkdm,v,&ghost);
    154:     if (ghost) continue;

    156:     DMNetworkGetComponent(networkdm,v,0,&type,(void**)&junction);
    157:     DMNetworkGetVariableOffset(networkdm,v,&varoffset);
    158:     juncx      = (PipeField*)(xarr + varoffset);
    159:     juncf      = (PetscScalar*)(farr + varoffset);

    161:     juncf[0] = -juncx[0].q;
    162:     juncf[1] =  juncx[0].q;

    164:     if (junction->type == RESERVOIR) { /* upstream reservoir */
    165:       juncf[0] = juncx[0].h - wash->H0;
    166:     }
    167:   }

    169:   /* Edge/pipe */
    170:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    171:   for (e=eStart; e<eEnd; e++) {
    172:     DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe);
    173:     DMNetworkGetVariableOffset(networkdm,e,&varoffset);
    174:     pipex    = (PipeField*)(xarr + varoffset);
    175:     pipexdot = (PipeField*)(xdotarr + varoffset);
    176:     pipef    = (PetscScalar*)(farr + varoffset);

    178:     /* Get some data into the pipe structure: note, some of these operations
    179:      * might be redundant. Will it consume too much time? */
    180:     pipe->dt   = dt;
    181:     pipe->xold = (PipeField*)(xoldarr + varoffset);

    183:     /* Evaluate F over this edge/pipe: pipef[1], ...,pipef[2*nend] */
    184:     DMDAGetLocalInfo(pipe->da,&info);
    185:     PipeIFunctionLocal_Lax(&info,t,pipex,pipexdot,pipef,pipe);

    187:     /* Get boundary values from connected vertices */
    188:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
    189:     vfrom = cone[0]; /* local ordering */
    190:     vto   = cone[1];
    191:     DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
    192:     DMNetworkGetVariableOffset(networkdm,vto,&offsetto);

    194:     /* Evaluate upstream boundary */
    195:     DMNetworkGetComponent(networkdm,vfrom,0,&type,(void**)&junction);
    196:     if (junction->type != JUNCTION && junction->type != RESERVOIR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"junction type is not supported");
    197:     juncx = (PipeField*)(xarr + offsetfrom);
    198:     juncf = (PetscScalar*)(farr + offsetfrom);

    200:     pipef[0] = pipex[0].h - juncx[0].h;
    201:     juncf[1] -= pipex[0].q;

    203:     /* Evaluate downstream boundary */
    204:     DMNetworkGetComponent(networkdm,vto,0,&type,(void**)&junction);
    205:     if (junction->type != JUNCTION && junction->type != VALVE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"junction type is not supported");
    206:     juncx = (PipeField*)(xarr + offsetto);
    207:     juncf = (PetscScalar*)(farr + offsetto);
    208:     nend  = pipe->nnodes - 1;

    210:     pipef[2*nend + 1] = pipex[nend].h - juncx[0].h;
    211:     juncf[0] += pipex[nend].q;
    212:   }

    214:   VecRestoreArrayRead(localX,&xarr);
    215:   VecRestoreArrayRead(localXdot,&xdotarr);
    216:   VecRestoreArray(localF,&farr);

    218:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
    219:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
    220:   DMRestoreLocalVector(networkdm,&localF);
    221:   /*
    222:    PetscPrintf(PETSC_COMM_WORLD("F:\n");
    223:    VecView(F,PETSC_VIEWER_STDOUT_WORLD);
    224:    */
    225:   return(0);
    226: }

    228: PetscErrorCode WASHSetInitialSolution(DM networkdm,Vec X,Wash wash)
    229: {
    231:   PetscInt       k,nx,vkey,vfrom,vto,offsetfrom,offsetto;
    232:   PetscInt       type,varoffset;
    233:   PetscInt       e,eStart,eEnd;
    234:   Vec            localX;
    235:   PetscScalar    *xarr;
    236:   Pipe           pipe;
    237:   Junction       junction;
    238:   const PetscInt *cone;
    239:   const PetscScalar *xarray;

    242:   VecSet(X,0.0);
    243:   DMGetLocalVector(networkdm,&localX);
    244:   VecGetArray(localX,&xarr);

    246:   /* Edge */
    247:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    248:   for (e=eStart; e<eEnd; e++) {
    249:     DMNetworkGetVariableOffset(networkdm,e,&varoffset);
    250:     DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe);

    252:     /* set initial values for this pipe */
    253:     PipeComputeSteadyState(pipe,wash->Q0,wash->H0);
    254:     VecGetSize(pipe->x,&nx);

    256:     VecGetArrayRead(pipe->x,&xarray);
    257:     /* copy pipe->x to xarray */
    258:     for (k=0; k<nx; k++) {
    259:       (xarr+varoffset)[k] = xarray[k];
    260:     }

    262:     /* set boundary values into vfrom and vto */
    263:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
    264:     vfrom = cone[0]; /* local ordering */
    265:     vto   = cone[1];
    266:     DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
    267:     DMNetworkGetVariableOffset(networkdm,vto,&offsetto);

    269:     /* if vform is a head vertex: */
    270:     DMNetworkGetComponent(networkdm,vfrom,0,&vkey,(void**)&junction);
    271:     if (junction->type == RESERVOIR) {
    272:       (xarr+offsetfrom)[1] = wash->H0; /* 1st H */
    273:     }

    275:     /* if vto is an end vertex: */
    276:     DMNetworkGetComponent(networkdm,vto,0,&vkey,(void**)&junction);
    277:     if (junction->type == VALVE) {
    278:       (xarr+offsetto)[0] = wash->QL; /* last Q */
    279:     }
    280:     VecRestoreArrayRead(pipe->x,&xarray);
    281:   }

    283:   VecRestoreArray(localX,&xarr);
    284:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
    285:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
    286:   DMRestoreLocalVector(networkdm,&localX);

    288: #if 0
    289:   PetscInt N;
    290:   VecGetSize(X,&N);
    291:   PetscPrintf(PETSC_COMM_WORLD,"initial solution %d:\n",N);
    292:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
    293: #endif
    294:   return(0);
    295: }

    297: PetscErrorCode TSDMNetworkMonitor(TS ts, PetscInt step, PetscReal t, Vec x, void *context)
    298: {
    299:   PetscErrorCode     ierr;
    300:   DMNetworkMonitor   monitor;

    303:   monitor = (DMNetworkMonitor)context;
    304:   DMNetworkMonitorView(monitor,x);
    305:   return(0);
    306: }

    308: PetscErrorCode PipesView(Vec X,DM networkdm,Wash wash)
    309: {
    310:   PetscErrorCode       ierr;
    311:   Pipe                 pipe;
    312:   PetscInt             key,Start,End;
    313:   PetscMPIInt          rank;
    314:   PetscInt             nx,nnodes,nidx,*idx1,*idx2,*idx1_h,*idx2_h,idx_start,i,k,k1,xstart,j1;
    315:   Vec                  Xq,Xh,localX;
    316:   IS                   is1_q,is2_q,is1_h,is2_h;
    317:   VecScatter           ctx_q,ctx_h;

    320:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

    322:   /* get num of local and global total nnodes */
    323:   nidx = wash->nnodes_loc;
    324:   MPIU_Allreduce(&nidx,&nx,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);

    326:   VecCreate(PETSC_COMM_WORLD,&Xq);
    327:   if (rank == 0) { /* all entries of Xq are in proc[0] */
    328:     VecSetSizes(Xq,nx,PETSC_DECIDE);
    329:   } else {
    330:     VecSetSizes(Xq,0,PETSC_DECIDE);
    331:   }
    332:   VecSetFromOptions(Xq);
    333:   VecSet(Xq,0.0);
    334:   VecDuplicate(Xq,&Xh);

    336:   DMGetLocalVector(networkdm,&localX);

    338:   /* set idx1 and idx2 */
    339:   PetscCalloc4(nidx,&idx1,nidx,&idx2,nidx,&idx1_h,nidx,&idx2_h);

    341:   DMNetworkGetEdgeRange(networkdm,&Start, &End);

    343:   VecGetOwnershipRange(X,&xstart,NULL);
    344:   k1 = 0;
    345:   j1 = 0;
    346:   for (i = Start; i < End; i++) {
    347:     DMNetworkGetComponent(networkdm,i,0,&key,(void**)&pipe);
    348:     nnodes = pipe->nnodes;
    349:     idx_start = pipe->id*nnodes;
    350:     for (k=0; k<nnodes; k++) {
    351:       idx1[k1] = xstart + j1*2*nnodes + 2*k;
    352:       idx2[k1] = idx_start + k;

    354:       idx1_h[k1] = xstart + j1*2*nnodes + 2*k + 1;
    355:       idx2_h[k1] = idx_start + k;
    356:       k1++;
    357:     }
    358:     j1++;
    359:   }

    361:   ISCreateGeneral(PETSC_COMM_SELF,nidx,idx1,PETSC_COPY_VALUES,&is1_q);
    362:   ISCreateGeneral(PETSC_COMM_SELF,nidx,idx2,PETSC_COPY_VALUES,&is2_q);
    363:   VecScatterCreate(X,is1_q,Xq,is2_q,&ctx_q);
    364:   VecScatterBegin(ctx_q,X,Xq,INSERT_VALUES,SCATTER_FORWARD);
    365:   VecScatterEnd(ctx_q,X,Xq,INSERT_VALUES,SCATTER_FORWARD);

    367:   ISCreateGeneral(PETSC_COMM_SELF,nidx,idx1_h,PETSC_COPY_VALUES,&is1_h);
    368:   ISCreateGeneral(PETSC_COMM_SELF,nidx,idx2_h,PETSC_COPY_VALUES,&is2_h);
    369:   VecScatterCreate(X,is1_h,Xh,is2_h,&ctx_h);
    370:   VecScatterBegin(ctx_h,X,Xh,INSERT_VALUES,SCATTER_FORWARD);
    371:   VecScatterEnd(ctx_h,X,Xh,INSERT_VALUES,SCATTER_FORWARD);

    373:   PetscPrintf(PETSC_COMM_WORLD,"Xq: \n");
    374:   VecView(Xq,PETSC_VIEWER_STDOUT_WORLD);
    375:   PetscPrintf(PETSC_COMM_WORLD,"Xh: \n");
    376:   VecView(Xh,PETSC_VIEWER_STDOUT_WORLD);

    378:   VecScatterDestroy(&ctx_q);
    379:   PetscFree4(idx1,idx2,idx1_h,idx2_h);
    380:   ISDestroy(&is1_q);
    381:   ISDestroy(&is2_q);

    383:   VecScatterDestroy(&ctx_h);
    384:   ISDestroy(&is1_h);
    385:   ISDestroy(&is2_h);

    387:   VecDestroy(&Xq);
    388:   VecDestroy(&Xh);
    389:   DMRestoreLocalVector(networkdm,&localX);
    390:   return(0);
    391: }

    393: PetscErrorCode WashNetworkCleanUp(Wash wash)
    394: {
    396:   PetscMPIInt    rank;

    399:   MPI_Comm_rank(wash->comm,&rank);
    400:   PetscFree(wash->edgelist);
    401:   PetscFree(wash->vtype);
    402:   if (!rank) {
    403:     PetscFree2(wash->junction,wash->pipe);
    404:   }
    405:   return(0);
    406: }

    408: PetscErrorCode WashNetworkCreate(MPI_Comm comm,PetscInt pipesCase,Wash *wash_ptr)
    409: {
    411:   PetscInt       npipes;
    412:   PetscMPIInt    rank;
    413:   Wash           wash=NULL;
    414:   PetscInt       i,numVertices,numEdges,*vtype;
    415:   PetscInt       *edgelist;
    416:   Junction       junctions=NULL;
    417:   Pipe           pipes=NULL;
    418:   PetscBool      washdist=PETSC_TRUE;

    421:   MPI_Comm_rank(comm,&rank);

    423:   PetscCalloc1(1,&wash);
    424:   wash->comm = comm;
    425:   *wash_ptr  = wash;
    426:   wash->Q0   = 0.477432; /* RESERVOIR */
    427:   wash->H0   = 150.0;
    428:   wash->HL   = 143.488;  /* VALVE */
    429:   wash->QL   = 0.0;
    430:   wash->nnodes_loc = 0;

    432:   numVertices = 0;
    433:   numEdges    = 0;
    434:   edgelist    = NULL;

    436:   /* proc[0] creates a sequential wash and edgelist */
    437:   PetscPrintf(PETSC_COMM_WORLD,"Setup pipesCase %D\n",pipesCase);

    439:   /* Set global number of pipes, edges, and junctions */
    440:   /*-------------------------------------------------*/
    441:   switch (pipesCase) {
    442:   case 0:
    443:     /* pipeCase 0: */
    444:     /* =================================================
    445:     (RESERVOIR) v0 --E0--> v1--E1--> v2 --E2-->v3 (VALVE)
    446:     ====================================================  */
    447:     npipes = 3;
    448:     PetscOptionsGetInt(NULL,NULL, "-npipes", &npipes, NULL);
    449:     wash->nedge   = npipes;
    450:     wash->nvertex = npipes + 1;

    452:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
    453:     numVertices = 0;
    454:     numEdges    = 0;
    455:     edgelist    = NULL;
    456:     if (!rank) {
    457:       numVertices = wash->nvertex;
    458:       numEdges    = wash->nedge;

    460:       PetscCalloc1(2*numEdges,&edgelist);
    461:       for (i=0; i<numEdges; i++) {
    462:         edgelist[2*i] = i; edgelist[2*i+1] = i+1;
    463:       }

    465:       /* Add network components */
    466:       /*------------------------*/
    467:       PetscCalloc2(numVertices,&junctions,numEdges,&pipes);

    469:       /* vertex */
    470:       for (i=0; i<numVertices; i++) {
    471:         junctions[i].id = i;
    472:         junctions[i].type = JUNCTION;
    473:       }

    475:       junctions[0].type             = RESERVOIR;
    476:       junctions[numVertices-1].type = VALVE;
    477:     }
    478:     break;
    479:   case 1:
    480:     /* pipeCase 1: */
    481:     /* ==========================
    482:                 v2 (VALVE)
    483:                 ^
    484:                 |
    485:                E2
    486:                 |
    487:     v0 --E0--> v3--E1--> v1
    488:   (RESERVOIR)            (RESERVOIR)
    489:     =============================  */
    490:     npipes = 3;
    491:     wash->nedge   = npipes;
    492:     wash->nvertex = npipes + 1;

    494:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
    495:     if (!rank) {
    496:       numVertices = wash->nvertex;
    497:       numEdges    = wash->nedge;

    499:       PetscCalloc1(2*numEdges,&edgelist);
    500:       edgelist[0] = 0; edgelist[1] = 3;  /* edge[0] */
    501:       edgelist[2] = 3; edgelist[3] = 1;  /* edge[1] */
    502:       edgelist[4] = 3; edgelist[5] = 2;  /* edge[2] */

    504:       /* Add network components */
    505:       /*------------------------*/
    506:       PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
    507:       /* vertex */
    508:       for (i=0; i<numVertices; i++) {
    509:         junctions[i].id   = i;
    510:         junctions[i].type = JUNCTION;
    511:       }

    513:       junctions[0].type = RESERVOIR;
    514:       junctions[1].type = VALVE;
    515:       junctions[2].type = VALVE;
    516:     }
    517:     break;
    518:   case 2:
    519:     /* pipeCase 2: */
    520:     /* ==========================
    521:     (RESERVOIR)  v2--> E2
    522:                        |
    523:             v0 --E0--> v3--E1--> v1
    524:     (RESERVOIR)               (VALVE)
    525:     =============================  */

    527:     /* Set application parameters -- to be used in function evalutions */
    528:     npipes = 3;
    529:     wash->nedge   = npipes;
    530:     wash->nvertex = npipes + 1;

    532:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
    533:     if (!rank) {
    534:       numVertices = wash->nvertex;
    535:       numEdges    = wash->nedge;

    537:       PetscCalloc1(2*numEdges,&edgelist);
    538:       edgelist[0] = 0; edgelist[1] = 3;  /* edge[0] */
    539:       edgelist[2] = 3; edgelist[3] = 1;  /* edge[1] */
    540:       edgelist[4] = 2; edgelist[5] = 3;  /* edge[2] */

    542:       /* Add network components */
    543:       /*------------------------*/
    544:       PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
    545:       /* vertex */
    546:       for (i=0; i<numVertices; i++) {
    547:         junctions[i].id = i;
    548:         junctions[i].type = JUNCTION;
    549:       }

    551:       junctions[0].type = RESERVOIR;
    552:       junctions[1].type = VALVE;
    553:       junctions[2].type = RESERVOIR;
    554:     }
    555:     break;
    556:   default:
    557:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"not done yet");
    558:   }

    560:   /* set edge global id */
    561:   for (i=0; i<numEdges; i++) pipes[i].id = i;

    563:   if (!rank) { /* set vtype for proc[0] */
    564:     PetscInt v;
    565:     PetscMalloc1(2*numEdges,&vtype);
    566:     for (i=0; i<2*numEdges; i++) {
    567:       v        = edgelist[i];
    568:       vtype[i] = junctions[v].type;
    569:     }
    570:     wash->vtype = vtype;
    571:   }

    573:   *wash_ptr      = wash;
    574:   wash->nedge    = numEdges;
    575:   wash->nvertex  = numVertices;
    576:   wash->edgelist = edgelist;
    577:   wash->junction = junctions;
    578:   wash->pipe     = pipes;

    580:   /* Distribute edgelist to other processors */
    581:   PetscOptionsGetBool(NULL,NULL,"-wash_distribute",&washdist,NULL);
    582:   if (washdist) {
    583:     /*
    584:      PetscPrintf(PETSC_COMM_WORLD," Distribute sequential wash ...\n");
    585:      */
    586:     WashNetworkDistribute(comm,wash);
    587:   }
    588:   return(0);
    589: }

    591: /* ------------------------------------------------------- */
    592: int main(int argc,char ** argv)
    593: {
    594:   PetscErrorCode    ierr;
    595:   Wash              wash;
    596:   Junction          junctions,junction;
    597:   Pipe              pipe,pipes;
    598:   PetscInt          KeyPipe,KeyJunction;
    599:   PetscInt          *edgelist = NULL,*edgelists[1],*vtype = NULL;
    600:   PetscInt          i,e,v,eStart,eEnd,vStart,vEnd,key;
    601:   PetscInt          vkey,type;
    602:   const PetscInt    *cone;
    603:   DM                networkdm;
    604:   PetscMPIInt       size,rank;
    605:   PetscReal         ftime;
    606:   Vec               X;
    607:   TS                ts;
    608:   PetscInt          steps=1;
    609:   TSConvergedReason reason;
    610:   PetscBool         viewpipes,monipipes=PETSC_FALSE,userJac=PETSC_TRUE,viewdm=PETSC_FALSE,viewX=PETSC_FALSE;
    611:   PetscBool         test=PETSC_FALSE;
    612:   PetscInt          pipesCase=0;
    613:   DMNetworkMonitor  monitor;
    614:   MPI_Comm          comm;

    616:   PetscInt          nedges,nvertices; /* local num of edges and vertices */
    617:   PetscInt          nnodes = 6,nv,ne;
    618:   const PetscInt    *vtx,*edge;

    620:   PetscInitialize(&argc,&argv,"pOption",help);if (ierr) return ierr;

    622:   /* Read runtime options */
    623:   PetscOptionsGetInt(NULL,NULL, "-case", &pipesCase, NULL);
    624:   PetscOptionsGetBool(NULL,NULL,"-user_Jac",&userJac,NULL);
    625:   PetscOptionsGetBool(NULL,NULL,"-pipe_monitor",&monipipes,NULL);
    626:   PetscOptionsGetBool(NULL,NULL,"-viewdm",&viewdm,NULL);
    627:   PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);
    628:   PetscOptionsGetInt(NULL,NULL, "-npipenodes", &nnodes, NULL);
    629:   PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);

    631:   /* Create networkdm */
    632:   /*------------------*/
    633:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
    634:   PetscObjectGetComm((PetscObject)networkdm,&comm);
    635:   MPI_Comm_rank(comm,&rank);
    636:   MPI_Comm_size(comm,&size);

    638:   if (size == 1 && monipipes) {
    639:     DMNetworkMonitorCreate(networkdm,&monitor);
    640:   }

    642:   /* Register the components in the network */
    643:   DMNetworkRegisterComponent(networkdm,"junctionstruct",sizeof(struct _p_Junction),&KeyJunction);
    644:   DMNetworkRegisterComponent(networkdm,"pipestruct",sizeof(struct _p_Pipe),&KeyPipe);

    646:   /* Create a distributed wash network (user-specific) */
    647:   WashNetworkCreate(comm,pipesCase,&wash);
    648:   nedges      = wash->nedge;
    649:   nvertices   = wash->nvertex; /* local num of vertices, excluding ghosts */
    650:   edgelist    = wash->edgelist;
    651:   vtype       = wash->vtype;
    652:   junctions   = wash->junction;
    653:   pipes       = wash->pipe;

    655: #if 0
    656:   PetscPrintf(PETSC_COMM_SELF,"[%d] after WashNetworkDistribute...ne %d, nv %d\n",rank,nedges,nvertices);
    657:   if (rank == 0) {
    658:     for (e=0; e<nedges; e++) PetscPrintf(PETSC_COMM_SELF," edge %d --> %d\n",edgelist[2*e],edgelist[2*e+1]);
    659:   }
    660: #endif

    662:   /* Set up the network layout */
    663:   DMNetworkSetSizes(networkdm,1,&nvertices,&nedges,0,NULL);

    665:   /* Add local edge connectivity */
    666:   edgelists[0] = edgelist;
    667:   DMNetworkSetEdgeList(networkdm,edgelists,NULL);
    668:   DMNetworkLayoutSetUp(networkdm);

    670:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    671:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
    672:   /* PetscPrintf(PETSC_COMM_SELF,"[%d] eStart/End: %d - %d; vStart/End: %d - %d\n",rank,eStart,eEnd,vStart,vEnd); */

    674:   /* Test DMNetworkGetSubnetworkInfo() */
    675:   if (test) {
    676:     DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edge);
    677:     if (ne != eEnd - eStart || nv != vEnd - vStart) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"ne %D or nv %D is incorrect",ne,nv);
    678:   }

    680:   if (rank) { /* junctions[] and pipes[] for proc[0] are allocated in WashNetworkCreate() */
    681:     /* vEnd - vStart = nvertices + num of ghost vertices! */
    682:     PetscCalloc2(vEnd - vStart,&junctions,nedges,&pipes);
    683:   }

    685:   /* Add Pipe component to all local edges */
    686:   for (e = eStart; e < eEnd; e++) {
    687:     if (test) {
    688:       if (e != edge[e]) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"e %D != edge %D from DMNetworkGetSubnetworkInfo()",e,edge[e]);
    689:     }

    691:     pipes[e-eStart].nnodes = nnodes;
    692:     DMNetworkAddComponent(networkdm,e,KeyPipe,&pipes[e-eStart]);

    694:     /* Add number of variables to each edge */
    695:     DMNetworkAddNumVariables(networkdm,e,2*pipes[e-eStart].nnodes);

    697:     if (size == 1 && monipipes) { /* Add monitor -- show Q_{pipes[e-eStart].id}? */
    698:       pipes[e-eStart].length = 600.0;
    699:       DMNetworkMonitorAdd(monitor, "Pipe Q", e, pipes[e-eStart].nnodes, 0, 2, 0.0,pipes[e-eStart].length, -0.8, 0.8, PETSC_TRUE);
    700:       DMNetworkMonitorAdd(monitor, "Pipe H", e, pipes[e-eStart].nnodes, 1, 2, 0.0,pipes[e-eStart].length, -400.0, 800.0, PETSC_TRUE);
    701:     }
    702:   }

    704:   /* Add Junction component to all local vertices, including ghost vertices! */
    705:   for (v = vStart; v < vEnd; v++) {
    706:     if (test) {
    707:       if (v != vtx[v-vStart]) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"v %D != vtx %D from DMNetworkGetSubnetworkInfo()",v,vtx[v-vStart]);
    708:     }

    710:     DMNetworkAddComponent(networkdm,v,KeyJunction,&junctions[v-vStart]);

    712:     /* Add number of variables to vertex */
    713:     DMNetworkAddNumVariables(networkdm,v,2);
    714:   }

    716:   if (size > 1) {  /* must be called before DMSetUp()???. Other partitioners do not work yet??? -- cause crash in proc[0]! */
    717:     DM               plexdm;
    718:     PetscPartitioner part;
    719:     DMNetworkGetPlex(networkdm,&plexdm);
    720:     DMPlexGetPartitioner(plexdm, &part);
    721:     PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);
    722:   }

    724:   /* Set up DM for use */
    725:   DMSetUp(networkdm);
    726:   if (viewdm) {
    727:     PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");
    728:     DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
    729:   }

    731:   /* Set user physical parameters to the components */
    732:   for (e = eStart; e < eEnd; e++) {
    733:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
    734:     /* vfrom */
    735:     DMNetworkGetComponent(networkdm,cone[0],0,&vkey,(void**)&junction);
    736:     junction->type = (VertexType)vtype[2*e];

    738:     /* vto */
    739:     DMNetworkGetComponent(networkdm,cone[1],0,&vkey,(void**)&junction);
    740:     junction->type = (VertexType)vtype[2*e+1];
    741:   }

    743:   WashNetworkCleanUp(wash);

    745:   /* Network partitioning and distribution of data */
    746:   DMNetworkDistribute(&networkdm,0);
    747:   if (viewdm) {
    748:     PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
    749:     DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
    750:   }

    752:   /* create vectors */
    753:   DMCreateGlobalVector(networkdm,&X);
    754:   DMCreateLocalVector(networkdm,&wash->localX);
    755:   DMCreateLocalVector(networkdm,&wash->localXdot);

    757:   /* PipeSetUp -- each process only sets its own pipes */
    758:   /*---------------------------------------------------*/
    759:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

    761:   userJac = PETSC_TRUE;
    762:   DMNetworkHasJacobian(networkdm,userJac,userJac);
    763:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    764:   for (e=eStart; e<eEnd; e++) { /* each edge has only one component, pipe */
    765:     DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe);

    767:     wash->nnodes_loc += pipe->nnodes; /* local total num of nodes, will be used by PipesView() */
    768:     PipeSetParameters(pipe,
    769:                              600.0,          /* length */
    770:                              0.5,            /* diameter */
    771:                              1200.0,         /* a */
    772:                              0.018);    /* friction */
    773:     PipeSetUp(pipe);

    775:     if (userJac) {
    776:       /* Create Jacobian matrix structures for a Pipe */
    777:       Mat            *J;
    778:       PipeCreateJacobian(pipe,NULL,&J);
    779:       DMNetworkEdgeSetMatrix(networkdm,e,J);
    780:     }
    781:   }

    783:   if (userJac) {
    784:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
    785:     for (v=vStart; v<vEnd; v++) {
    786:       Mat            *J;
    787:       JunctionCreateJacobian(networkdm,v,NULL,&J);
    788:       DMNetworkVertexSetMatrix(networkdm,v,J);

    790:       DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction);
    791:       junction->jacobian = J;
    792:     }
    793:   }

    795:   /* Test DMNetworkGetSubnetworkInfo() */
    796:   if (test) {
    797:     DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    798:     DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edge);
    799:     if (ne != eEnd - eStart || nv != vEnd - vStart) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"ne %D or nv %D is incorrect",ne,nv);

    801:     for (e = eStart; e < eEnd; e++) {
    802:       if (e != edge[e]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"e %D != edge %D from DMNetworkGetSubnetworkInfo()",e,edge[e]);
    803:     }
    804:     for (v = vStart; v < vEnd; v++) {
    805:       if (v != vtx[v-vStart]) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"v %D != vtx %D from DMNetworkGetSubnetworkInfo()",v,vtx[v-vStart]);
    806:     }
    807:   }

    809:   /* Setup solver                                           */
    810:   /*--------------------------------------------------------*/
    811:   TSCreate(PETSC_COMM_WORLD,&ts);

    813:   TSSetDM(ts,(DM)networkdm);
    814:   TSSetIFunction(ts,NULL,WASHIFunction,wash);

    816:   TSSetMaxSteps(ts,steps);
    817:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
    818:   TSSetTimeStep(ts,0.1);
    819:   TSSetType(ts,TSBEULER);
    820:   if (size == 1 && monipipes) {
    821:     TSMonitorSet(ts, TSDMNetworkMonitor, monitor, NULL);
    822:   }
    823:   TSSetFromOptions(ts);

    825:   WASHSetInitialSolution(networkdm,X,wash);

    827:   TSSolve(ts,X);

    829:   TSGetSolveTime(ts,&ftime);
    830:   TSGetStepNumber(ts,&steps);
    831:   TSGetConvergedReason(ts,&reason);
    832:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
    833:   if (viewX) {
    834:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
    835:   }

    837:   viewpipes = PETSC_FALSE;
    838:   PetscOptionsGetBool(NULL,NULL, "-Jac_view", &viewpipes,NULL);
    839:   if (viewpipes) {
    840:     SNES snes;
    841:     Mat  Jac;
    842:     TSGetSNES(ts,&snes);
    843:     SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
    844:     MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
    845:   }

    847:   /* View solution q and h */
    848:   /* --------------------- */
    849:   viewpipes = PETSC_FALSE;
    850:   PetscOptionsGetBool(NULL,NULL, "-pipe_view", &viewpipes,NULL);
    851:   if (viewpipes) {
    852:     PipesView(X,networkdm,wash);
    853:   }

    855:   /* Free spaces */
    856:   /* ----------- */
    857:   TSDestroy(&ts);
    858:   VecDestroy(&X);
    859:   VecDestroy(&wash->localX);
    860:   VecDestroy(&wash->localXdot);

    862:   /* Destroy objects from each pipe that are created in PipeSetUp() */
    863:   DMNetworkGetEdgeRange(networkdm,&eStart, &eEnd);
    864:   for (i = eStart; i < eEnd; i++) {
    865:     DMNetworkGetComponent(networkdm,i,0,&key,(void**)&pipe);
    866:     PipeDestroy(&pipe);
    867:   }
    868:   if (userJac) {
    869:     for (v=vStart; v<vEnd; v++) {
    870:       DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction);
    871:       JunctionDestroyJacobian(networkdm,v,junction);
    872:     }
    873:   }

    875:   if (size == 1 && monipipes) {
    876:     DMNetworkMonitorDestroy(&monitor);
    877:   }
    878:   DMDestroy(&networkdm);
    879:   PetscFree(wash);

    881:   if (rank) {
    882:     PetscFree2(junctions,pipes);
    883:   }
    884:   PetscFinalize();
    885:   return ierr;
    886: }

    888: /*TEST

    890:    build:
    891:      depends: pipeInterface.c pipeImpls.c

    893:    test:
    894:       args: -ts_monitor -case 1 -ts_max_steps 1 -options_left no -viewX -test
    895:       localrunfiles: pOption
    896:       output_file: output/pipes1_1.out

    898:    test:
    899:       suffix: 2
    900:       nsize: 2
    901:       requires: mumps
    902:       args: -ts_monitor -case 1 -ts_max_steps 1 -petscpartitioner_type simple -options_left no -viewX -test
    903:       localrunfiles: pOption
    904:       output_file: output/pipes1_2.out

    906:    test:
    907:       suffix: 3
    908:       nsize: 2
    909:       requires: mumps
    910:       args: -ts_monitor -case 0 -ts_max_steps 1 -petscpartitioner_type simple -options_left no -viewX -test
    911:       localrunfiles: pOption
    912:       output_file: output/pipes1_3.out

    914:    test:
    915:       suffix: 4
    916:       args: -ts_monitor -case 2 -ts_max_steps 1 -options_left no -viewX -test
    917:       localrunfiles: pOption
    918:       output_file: output/pipes1_4.out

    920:    test:
    921:       suffix: 5
    922:       nsize: 3
    923:       requires: mumps
    924:       args: -ts_monitor -case 2 -ts_max_steps 10 -petscpartitioner_type simple -options_left no -viewX -test
    925:       localrunfiles: pOption
    926:       output_file: output/pipes1_5.out

    928:    test:
    929:       suffix: 6
    930:       nsize: 2
    931:       requires: mumps
    932:       args: -ts_monitor -case 1 -ts_max_steps 1 -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX -test
    933:       localrunfiles: pOption
    934:       output_file: output/pipes1_6.out

    936:    test:
    937:       suffix: 7
    938:       nsize: 2
    939:       requires: mumps
    940:       args: -ts_monitor -case 2 -ts_max_steps 1 -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX -test
    941:       localrunfiles: pOption
    942:       output_file: output/pipes1_7.out

    944: TEST*/