Actual source code: pipes.c

  1: static char help[] = "This example demonstrates DMNetwork. It is used for testing parallel generation of dmnetwork, then redistribute. \n\\n";
  2: /*
  3:   Example: mpiexec -n <np> ./pipes -ts_max_steps 10
  4: */

  6: #include "wash.h"

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

 14:    Output Parameter:
 15: .  wash - wash context with nedge, nvertex and edgelist distributed

 17:    Note: The routine is used for testing parallel generation of dmnetwork, then redistribute.
 18: */
 19: PetscErrorCode WashNetworkDistribute(MPI_Comm comm,Wash wash)
 20: {
 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;

 25:   MPI_Comm_size(comm,&size);
 26:   if (size == 1) return 0;

 28:   MPI_Comm_rank(comm,&rank);
 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 == 0) {
 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 == 0) {
 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 == 0) {
 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: {
101:   Wash           wash=(Wash)ctx;
102:   DM             networkdm;
103:   Vec            localX,localXdot,localF, localXold;
104:   const PetscInt *cone;
105:   PetscInt       vfrom,vto,offsetfrom,offsetto,varoffset;
106:   PetscInt       v,vStart,vEnd,e,eStart,eEnd;
107:   PetscInt       nend,type;
108:   PetscBool      ghost;
109:   PetscScalar    *farr,*juncf, *pipef;
110:   PetscReal      dt;
111:   Pipe           pipe;
112:   PipeField      *pipex,*pipexdot,*juncx;
113:   Junction       junction;
114:   DMDALocalInfo  info;
115:   const PetscScalar *xarr,*xdotarr, *xoldarr;

117:   localX    = wash->localX;
118:   localXdot = wash->localXdot;

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

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

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

133:   DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot);
134:   DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot);

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

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

154:     DMNetworkGetComponent(networkdm,v,0,&type,(void**)&junction,NULL);
155:     DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&varoffset);
156:     juncx      = (PipeField*)(xarr + varoffset);
157:     juncf      = (PetscScalar*)(farr + varoffset);

159:     juncf[0] = -juncx[0].q;
160:     juncf[1] =  juncx[0].q;

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

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

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

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

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

192:     /* Evaluate upstream boundary */
193:     DMNetworkGetComponent(networkdm,vfrom,0,&type,(void**)&junction,NULL);
195:     juncx = (PipeField*)(xarr + offsetfrom);
196:     juncf = (PetscScalar*)(farr + offsetfrom);

198:     pipef[0] = pipex[0].h - juncx[0].h;
199:     juncf[1] -= pipex[0].q;

201:     /* Evaluate downstream boundary */
202:     DMNetworkGetComponent(networkdm,vto,0,&type,(void**)&junction,NULL);
204:     juncx = (PipeField*)(xarr + offsetto);
205:     juncf = (PetscScalar*)(farr + offsetto);
206:     nend  = pipe->nnodes - 1;

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

212:   VecRestoreArrayRead(localX,&xarr);
213:   VecRestoreArrayRead(localXdot,&xdotarr);
214:   VecRestoreArray(localF,&farr);

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

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

238:   VecSet(X,0.0);
239:   DMGetLocalVector(networkdm,&localX);
240:   VecGetArray(localX,&xarr);

242:   /* Edge */
243:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
244:   for (e=eStart; e<eEnd; e++) {
245:     DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&varoffset);
246:     DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe,NULL);

248:     /* set initial values for this pipe */
249:     PipeComputeSteadyState(pipe,wash->Q0,wash->H0);
250:     VecGetSize(pipe->x,&nx);

252:     VecGetArrayRead(pipe->x,&xarray);
253:     /* copy pipe->x to xarray */
254:     for (k=0; k<nx; k++) {
255:       (xarr+varoffset)[k] = xarray[k];
256:     }

258:     /* set boundary values into vfrom and vto */
259:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
260:     vfrom = cone[0]; /* local ordering */
261:     vto   = cone[1];
262:     DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
263:     DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);

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

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

279:   VecRestoreArray(localX,&xarr);
280:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
281:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
282:   DMRestoreLocalVector(networkdm,&localX);

284: #if 0
285:   PetscInt N;
286:   VecGetSize(X,&N);
287:   PetscPrintf(PETSC_COMM_WORLD,"initial solution %d:\n",N);
288:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
289: #endif
290:   return 0;
291: }

293: PetscErrorCode TSDMNetworkMonitor(TS ts, PetscInt step, PetscReal t, Vec x, void *context)
294: {
295:   DMNetworkMonitor   monitor;

297:   monitor = (DMNetworkMonitor)context;
298:   DMNetworkMonitorView(monitor,x);
299:   return 0;
300: }

302: PetscErrorCode PipesView(DM networkdm,PetscInt KeyPipe,Vec X)
303: {
304:   PetscInt       i,numkeys=1,*blocksize,*numselectedvariable,**selectedvariables,n;
305:   IS             isfrom_q,isfrom_h,isfrom;
306:   Vec            Xto;
307:   VecScatter     ctx;
308:   MPI_Comm       comm;

310:   PetscObjectGetComm((PetscObject)networkdm,&comm);

312:   /* 1. Create isfrom_q for q-variable of pipes */
313:   PetscMalloc3(numkeys,&blocksize,numkeys,&numselectedvariable,numkeys,&selectedvariables);
314:   for (i=0; i<numkeys; i++) {
315:     blocksize[i]           = 2;
316:     numselectedvariable[i] = 1;
317:     PetscMalloc1(numselectedvariable[i],&selectedvariables[i]);
318:     selectedvariables[i][0] = 0; /* q-variable */
319:   }
320:   DMNetworkCreateIS(networkdm,numkeys,&KeyPipe,blocksize,numselectedvariable,selectedvariables,&isfrom_q);

322:   /* 2. Create Xto and isto */
323:   ISGetLocalSize(isfrom_q, &n);
324:   VecCreate(comm,&Xto);
325:   VecSetSizes(Xto,n,PETSC_DECIDE);
326:   VecSetFromOptions(Xto);
327:   VecSet(Xto,0.0);

329:   /* 3. Create scatter */
330:   VecScatterCreate(X,isfrom_q,Xto,NULL,&ctx);

332:   /* 4. Scatter to Xq */
333:   VecScatterBegin(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
334:   VecScatterEnd(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
335:   VecScatterDestroy(&ctx);
336:   ISDestroy(&isfrom_q);

338:   PetscPrintf(PETSC_COMM_WORLD,"Xq:\n");
339:   VecView(Xto,PETSC_VIEWER_STDOUT_WORLD);

341:   /* 5. Create isfrom_h for h-variable of pipes; Create scatter; Scatter to Xh */
342:   for (i=0; i<numkeys; i++) {
343:     selectedvariables[i][0] = 1; /* h-variable */
344:   }
345:   DMNetworkCreateIS(networkdm,numkeys,&KeyPipe,blocksize,numselectedvariable,selectedvariables,&isfrom_h);

347:   VecScatterCreate(X,isfrom_h,Xto,NULL,&ctx);
348:   VecScatterBegin(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
349:   VecScatterEnd(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
350:   VecScatterDestroy(&ctx);
351:   ISDestroy(&isfrom_h);

353:   PetscPrintf(PETSC_COMM_WORLD,"Xh:\n");
354:   VecView(Xto,PETSC_VIEWER_STDOUT_WORLD);
355:   VecDestroy(&Xto);

357:   /* 6. Create isfrom for all pipe variables; Create scatter; Scatter to Xpipes */
358:   for (i=0; i<numkeys; i++) {
359:     blocksize[i] = -1; /* select all the variables of the i-th component */
360:   }
361:   DMNetworkCreateIS(networkdm,numkeys,&KeyPipe,blocksize,NULL,NULL,&isfrom);
362:   ISDestroy(&isfrom);
363:   DMNetworkCreateIS(networkdm,numkeys,&KeyPipe,NULL,NULL,NULL,&isfrom);

365:   ISGetLocalSize(isfrom, &n);
366:   VecCreate(comm,&Xto);
367:   VecSetSizes(Xto,n,PETSC_DECIDE);
368:   VecSetFromOptions(Xto);
369:   VecSet(Xto,0.0);

371:   VecScatterCreate(X,isfrom,Xto,NULL,&ctx);
372:   VecScatterBegin(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
373:   VecScatterEnd(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
374:   VecScatterDestroy(&ctx);
375:   ISDestroy(&isfrom);

377:   PetscPrintf(PETSC_COMM_WORLD,"Xpipes:\n");
378:   VecView(Xto,PETSC_VIEWER_STDOUT_WORLD);

380:   /* 7. Free spaces */
381:   for (i=0; i<numkeys; i++) {
382:     PetscFree(selectedvariables[i]);
383:   }
384:   PetscFree3(blocksize,numselectedvariable,selectedvariables);
385:   VecDestroy(&Xto);
386:   return 0;
387: }

389: PetscErrorCode ISJunctionsView(DM networkdm,PetscInt KeyJunc)
390: {
391:   PetscInt       numkeys=1;
392:   IS             isfrom;
393:   MPI_Comm       comm;
394:   PetscMPIInt    rank;

396:   PetscObjectGetComm((PetscObject)networkdm,&comm);
397:   MPI_Comm_rank(comm,&rank);

399:   /* Create a global isfrom for all junction variables */
400:   DMNetworkCreateIS(networkdm,numkeys,&KeyJunc,NULL,NULL,NULL,&isfrom);
401:   PetscPrintf(PETSC_COMM_WORLD,"ISJunctions:\n");
402:   ISView(isfrom,PETSC_VIEWER_STDOUT_WORLD);
403:   ISDestroy(&isfrom);

405:   /* Create a local isfrom for all junction variables */
406:   DMNetworkCreateLocalIS(networkdm,numkeys,&KeyJunc,NULL,NULL,NULL,&isfrom);
407:   if (!rank) {
408:     PetscPrintf(PETSC_COMM_SELF,"[%d] ISLocalJunctions:\n",rank);
409:     ISView(isfrom,PETSC_VIEWER_STDOUT_SELF);
410:   }
411:   ISDestroy(&isfrom);
412:   return 0;
413: }

415: PetscErrorCode WashNetworkCleanUp(Wash wash)
416: {
417:   PetscMPIInt    rank;

419:   MPI_Comm_rank(wash->comm,&rank);
420:   PetscFree(wash->edgelist);
421:   PetscFree(wash->vtype);
422:   if (rank == 0) {
423:     PetscFree2(wash->junction,wash->pipe);
424:   }
425:   return 0;
426: }

428: PetscErrorCode WashNetworkCreate(MPI_Comm comm,PetscInt pipesCase,Wash *wash_ptr)
429: {
430:   PetscInt       npipes;
431:   PetscMPIInt    rank;
432:   Wash           wash=NULL;
433:   PetscInt       i,numVertices,numEdges,*vtype;
434:   PetscInt       *edgelist;
435:   Junction       junctions=NULL;
436:   Pipe           pipes=NULL;
437:   PetscBool      washdist=PETSC_TRUE;

439:   MPI_Comm_rank(comm,&rank);

441:   PetscCalloc1(1,&wash);
442:   wash->comm = comm;
443:   *wash_ptr  = wash;
444:   wash->Q0   = 0.477432; /* RESERVOIR */
445:   wash->H0   = 150.0;
446:   wash->HL   = 143.488;  /* VALVE */
447:   wash->QL   = 0.0;
448:   wash->nnodes_loc = 0;

450:   numVertices = 0;
451:   numEdges    = 0;
452:   edgelist    = NULL;

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

457:   /* Set global number of pipes, edges, and junctions */
458:   /*-------------------------------------------------*/
459:   switch (pipesCase) {
460:   case 0:
461:     /* pipeCase 0: */
462:     /* =================================================
463:     (RESERVOIR) v0 --E0--> v1--E1--> v2 --E2-->v3 (VALVE)
464:     ====================================================  */
465:     npipes = 3;
466:     PetscOptionsGetInt(NULL,NULL, "-npipes", &npipes, NULL);
467:     wash->nedge   = npipes;
468:     wash->nvertex = npipes + 1;

470:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
471:     numVertices = 0;
472:     numEdges    = 0;
473:     edgelist    = NULL;
474:     if (rank == 0) {
475:       numVertices = wash->nvertex;
476:       numEdges    = wash->nedge;

478:       PetscCalloc1(2*numEdges,&edgelist);
479:       for (i=0; i<numEdges; i++) {
480:         edgelist[2*i] = i; edgelist[2*i+1] = i+1;
481:       }

483:       /* Add network components */
484:       /*------------------------*/
485:       PetscCalloc2(numVertices,&junctions,numEdges,&pipes);

487:       /* vertex */
488:       for (i=0; i<numVertices; i++) {
489:         junctions[i].id = i;
490:         junctions[i].type = JUNCTION;
491:       }

493:       junctions[0].type             = RESERVOIR;
494:       junctions[numVertices-1].type = VALVE;
495:     }
496:     break;
497:   case 1:
498:     /* pipeCase 1: */
499:     /* ==========================
500:                 v2 (VALVE)
501:                 ^
502:                 |
503:                E2
504:                 |
505:     v0 --E0--> v3--E1--> v1
506:   (RESERVOIR)            (RESERVOIR)
507:     =============================  */
508:     npipes = 3;
509:     wash->nedge   = npipes;
510:     wash->nvertex = npipes + 1;

512:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
513:     if (rank == 0) {
514:       numVertices = wash->nvertex;
515:       numEdges    = wash->nedge;

517:       PetscCalloc1(2*numEdges,&edgelist);
518:       edgelist[0] = 0; edgelist[1] = 3;  /* edge[0] */
519:       edgelist[2] = 3; edgelist[3] = 1;  /* edge[1] */
520:       edgelist[4] = 3; edgelist[5] = 2;  /* edge[2] */

522:       /* Add network components */
523:       /*------------------------*/
524:       PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
525:       /* vertex */
526:       for (i=0; i<numVertices; i++) {
527:         junctions[i].id   = i;
528:         junctions[i].type = JUNCTION;
529:       }

531:       junctions[0].type = RESERVOIR;
532:       junctions[1].type = VALVE;
533:       junctions[2].type = VALVE;
534:     }
535:     break;
536:   case 2:
537:     /* pipeCase 2: */
538:     /* ==========================
539:     (RESERVOIR)  v2--> E2
540:                        |
541:             v0 --E0--> v3--E1--> v1
542:     (RESERVOIR)               (VALVE)
543:     =============================  */

545:     /* Set application parameters -- to be used in function evalutions */
546:     npipes = 3;
547:     wash->nedge   = npipes;
548:     wash->nvertex = npipes + 1;

550:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
551:     if (rank == 0) {
552:       numVertices = wash->nvertex;
553:       numEdges    = wash->nedge;

555:       PetscCalloc1(2*numEdges,&edgelist);
556:       edgelist[0] = 0; edgelist[1] = 3;  /* edge[0] */
557:       edgelist[2] = 3; edgelist[3] = 1;  /* edge[1] */
558:       edgelist[4] = 2; edgelist[5] = 3;  /* edge[2] */

560:       /* Add network components */
561:       /*------------------------*/
562:       PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
563:       /* vertex */
564:       for (i=0; i<numVertices; i++) {
565:         junctions[i].id = i;
566:         junctions[i].type = JUNCTION;
567:       }

569:       junctions[0].type = RESERVOIR;
570:       junctions[1].type = VALVE;
571:       junctions[2].type = RESERVOIR;
572:     }
573:     break;
574:   default:
575:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"not done yet");
576:   }

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

581:   if (rank == 0) { /* set vtype for proc[0] */
582:     PetscInt v;
583:     PetscMalloc1(2*numEdges,&vtype);
584:     for (i=0; i<2*numEdges; i++) {
585:       v        = edgelist[i];
586:       vtype[i] = junctions[v].type;
587:     }
588:     wash->vtype = vtype;
589:   }

591:   *wash_ptr      = wash;
592:   wash->nedge    = numEdges;
593:   wash->nvertex  = numVertices;
594:   wash->edgelist = edgelist;
595:   wash->junction = junctions;
596:   wash->pipe     = pipes;

598:   /* Distribute edgelist to other processors */
599:   PetscOptionsGetBool(NULL,NULL,"-wash_distribute",&washdist,NULL);
600:   if (washdist) {
601:     /*
602:      PetscPrintf(PETSC_COMM_WORLD," Distribute sequential wash ...\n");
603:      */
604:     WashNetworkDistribute(comm,wash);
605:   }
606:   return 0;
607: }

609: /* ------------------------------------------------------- */
610: int main(int argc,char ** argv)
611: {
612:   Wash              wash;
613:   Junction          junctions,junction;
614:   Pipe              pipe,pipes;
615:   PetscInt          KeyPipe,KeyJunction,*edgelist = NULL,*vtype = NULL;
616:   PetscInt          i,e,v,eStart,eEnd,vStart,vEnd,key,vkey,type;
617:   PetscInt          steps=1,nedges,nnodes=6;
618:   const PetscInt    *cone;
619:   DM                networkdm;
620:   PetscMPIInt       size,rank;
621:   PetscReal         ftime;
622:   Vec               X;
623:   TS                ts;
624:   TSConvergedReason reason;
625:   PetscBool         viewpipes,viewjuncs,monipipes=PETSC_FALSE,userJac=PETSC_TRUE,viewdm=PETSC_FALSE,viewX=PETSC_FALSE;
626:   PetscInt          pipesCase=0;
627:   DMNetworkMonitor  monitor;
628:   MPI_Comm          comm;

630:   PetscInitialize(&argc,&argv,"pOption",help);

632:   /* Read runtime options */
633:   PetscOptionsGetInt(NULL,NULL, "-case", &pipesCase, NULL);
634:   PetscOptionsGetBool(NULL,NULL,"-user_Jac",&userJac,NULL);
635:   PetscOptionsGetBool(NULL,NULL,"-pipe_monitor",&monipipes,NULL);
636:   PetscOptionsGetBool(NULL,NULL,"-viewdm",&viewdm,NULL);
637:   PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);
638:   PetscOptionsGetInt(NULL,NULL, "-npipenodes", &nnodes, NULL);

640:   /* Create networkdm */
641:   /*------------------*/
642:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
643:   PetscObjectGetComm((PetscObject)networkdm,&comm);
644:   MPI_Comm_rank(comm,&rank);
645:   MPI_Comm_size(comm,&size);

647:   if (size == 1 && monipipes) {
648:     DMNetworkMonitorCreate(networkdm,&monitor);
649:   }

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

655:   /* Create a distributed wash network (user-specific) */
656:   WashNetworkCreate(comm,pipesCase,&wash);
657:   nedges      = wash->nedge;
658:   edgelist    = wash->edgelist;
659:   vtype       = wash->vtype;
660:   junctions   = wash->junction;
661:   pipes       = wash->pipe;

663:   /* Set up the network layout */
664:   DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
665:   DMNetworkAddSubnetwork(networkdm,NULL,nedges,edgelist,NULL);

667:   DMNetworkLayoutSetUp(networkdm);

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

673:   if (rank) { /* junctions[] and pipes[] for proc[0] are allocated in WashNetworkCreate() */
674:     /* vEnd - vStart = nvertices + number of ghost vertices! */
675:     PetscCalloc2(vEnd - vStart,&junctions,nedges,&pipes);
676:   }

678:   /* Add Pipe component and number of variables to all local edges */
679:   for (e = eStart; e < eEnd; e++) {
680:     pipes[e-eStart].nnodes = nnodes;
681:     DMNetworkAddComponent(networkdm,e,KeyPipe,&pipes[e-eStart],2*pipes[e-eStart].nnodes);

683:     if (size == 1 && monipipes) { /* Add monitor -- show Q_{pipes[e-eStart].id}? */
684:       pipes[e-eStart].length = 600.0;
685:       DMNetworkMonitorAdd(monitor, "Pipe Q", e, pipes[e-eStart].nnodes, 0, 2, 0.0,pipes[e-eStart].length, -0.8, 0.8, PETSC_TRUE);
686:       DMNetworkMonitorAdd(monitor, "Pipe H", e, pipes[e-eStart].nnodes, 1, 2, 0.0,pipes[e-eStart].length, -400.0, 800.0, PETSC_TRUE);
687:     }
688:   }

690:   /* Add Junction component and number of variables to all local vertices */
691:   for (v = vStart; v < vEnd; v++) {
692:     DMNetworkAddComponent(networkdm,v,KeyJunction,&junctions[v-vStart],2);
693:   }

695:   if (size > 1) {  /* must be called before DMSetUp()???. Other partitioners do not work yet??? -- cause crash in proc[0]! */
696:     DM               plexdm;
697:     PetscPartitioner part;
698:     DMNetworkGetPlex(networkdm,&plexdm);
699:     DMPlexGetPartitioner(plexdm, &part);
700:     PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);
701:     PetscOptionsSetValue(NULL,"-dm_plex_csr_alg","mat"); /* for parmetis */
702:   }

704:   /* Set up DM for use */
705:   DMSetUp(networkdm);
706:   if (viewdm) {
707:     PetscPrintf(PETSC_COMM_WORLD,"\nOriginal networkdm, DMView:\n");
708:     DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
709:   }

711:   /* Set user physical parameters to the components */
712:   for (e = eStart; e < eEnd; e++) {
713:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
714:     /* vfrom */
715:     DMNetworkGetComponent(networkdm,cone[0],0,&vkey,(void**)&junction,NULL);
716:     junction->type = (VertexType)vtype[2*e];

718:     /* vto */
719:     DMNetworkGetComponent(networkdm,cone[1],0,&vkey,(void**)&junction,NULL);
720:     junction->type = (VertexType)vtype[2*e+1];
721:   }

723:   WashNetworkCleanUp(wash);

725:   /* Network partitioning and distribution of data */
726:   DMNetworkDistribute(&networkdm,0);
727:   if (viewdm) {
728:     PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
729:     DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
730:   }

732:   /* create vectors */
733:   DMCreateGlobalVector(networkdm,&X);
734:   DMCreateLocalVector(networkdm,&wash->localX);
735:   DMCreateLocalVector(networkdm,&wash->localXdot);

737:   /* PipeSetUp -- each process only sets its own pipes */
738:   /*---------------------------------------------------*/
739:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

741:   userJac = PETSC_TRUE;
742:   DMNetworkHasJacobian(networkdm,userJac,userJac);
743:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
744:   for (e=eStart; e<eEnd; e++) { /* each edge has only one component, pipe */
745:     DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe,NULL);

747:     wash->nnodes_loc += pipe->nnodes; /* local total number of nodes, will be used by PipesView() */
748:     PetscCall(PipeSetParameters(pipe,
749:                               600.0,   /* length   */
750:                               0.5,     /* diameter */
751:                               1200.0,  /* a        */
752:                               0.018)); /* friction */
753:     PipeSetUp(pipe);

755:     if (userJac) {
756:       /* Create Jacobian matrix structures for a Pipe */
757:       Mat            *J;
758:       PipeCreateJacobian(pipe,NULL,&J);
759:       DMNetworkEdgeSetMatrix(networkdm,e,J);
760:     }
761:   }

763:   if (userJac) {
764:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
765:     for (v=vStart; v<vEnd; v++) {
766:       Mat            *J;
767:       JunctionCreateJacobian(networkdm,v,NULL,&J);
768:       DMNetworkVertexSetMatrix(networkdm,v,J);

770:       DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction,NULL);
771:       junction->jacobian = J;
772:     }
773:   }

775:   /* Setup solver                                           */
776:   /*--------------------------------------------------------*/
777:   TSCreate(PETSC_COMM_WORLD,&ts);

779:   TSSetDM(ts,(DM)networkdm);
780:   TSSetIFunction(ts,NULL,WASHIFunction,wash);

782:   TSSetMaxSteps(ts,steps);
783:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
784:   TSSetTimeStep(ts,0.1);
785:   TSSetType(ts,TSBEULER);
786:   if (size == 1 && monipipes) TSMonitorSet(ts, TSDMNetworkMonitor, monitor, NULL);
787:   TSSetFromOptions(ts);

789:   WASHSetInitialSolution(networkdm,X,wash);

791:   TSSolve(ts,X);

793:   TSGetSolveTime(ts,&ftime);
794:   TSGetStepNumber(ts,&steps);
795:   TSGetConvergedReason(ts,&reason);
796:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
797:   if (viewX) {
798:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
799:   }

801:   viewpipes = PETSC_FALSE;
802:   PetscOptionsGetBool(NULL,NULL, "-Jac_view", &viewpipes,NULL);
803:   if (viewpipes) {
804:     SNES snes;
805:     Mat  Jac;
806:     TSGetSNES(ts,&snes);
807:     SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
808:     MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
809:   }

811:   /* View solutions */
812:   /* -------------- */
813:   viewpipes = PETSC_FALSE;
814:   PetscOptionsGetBool(NULL,NULL, "-pipe_view", &viewpipes,NULL);
815:   if (viewpipes) {
816:     PipesView(networkdm,KeyPipe,X);
817:   }

819:   /* Test IS */
820:   viewjuncs = PETSC_FALSE;
821:   PetscOptionsGetBool(NULL,NULL, "-isjunc_view", &viewjuncs,NULL);
822:   if (viewjuncs) {
823:     ISJunctionsView(networkdm,KeyJunction);
824:   }

826:   /* Free spaces */
827:   /* ----------- */
828:   TSDestroy(&ts);
829:   VecDestroy(&X);
830:   VecDestroy(&wash->localX);
831:   VecDestroy(&wash->localXdot);

833:   /* Destroy objects from each pipe that are created in PipeSetUp() */
834:   DMNetworkGetEdgeRange(networkdm,&eStart, &eEnd);
835:   for (i = eStart; i < eEnd; i++) {
836:     DMNetworkGetComponent(networkdm,i,0,&key,(void**)&pipe,NULL);
837:     PipeDestroy(&pipe);
838:   }
839:   if (userJac) {
840:     for (v=vStart; v<vEnd; v++) {
841:       DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction,NULL);
842:       JunctionDestroyJacobian(networkdm,v,junction);
843:     }
844:   }

846:   if (size == 1 && monipipes) {
847:     DMNetworkMonitorDestroy(&monitor);
848:   }
849:   DMDestroy(&networkdm);
850:   PetscFree(wash);

852:   if (rank) {
853:     PetscFree2(junctions,pipes);
854:   }
855:   PetscFinalize();
856:   return 0;
857: }

859: /*TEST

861:    build:
862:      depends: pipeInterface.c pipeImpls.c
863:      requires: mumps

865:    test:
866:       args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -options_left no -viewX
867:       localrunfiles: pOption
868:       output_file: output/pipes_1.out

870:    test:
871:       suffix: 2
872:       nsize: 2
873:       args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
874:       localrunfiles: pOption
875:       output_file: output/pipes_2.out

877:    test:
878:       suffix: 3
879:       nsize: 2
880:       args: -ts_monitor -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
881:       localrunfiles: pOption
882:       output_file: output/pipes_3.out

884:    test:
885:       suffix: 4
886:       args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -options_left no -viewX
887:       localrunfiles: pOption
888:       output_file: output/pipes_4.out

890:    test:
891:       suffix: 5
892:       nsize: 3
893:       args: -ts_monitor -case 2 -ts_max_steps 10 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
894:       localrunfiles: pOption
895:       output_file: output/pipes_5.out

897:    test:
898:       suffix: 6
899:       nsize: 2
900:       args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX
901:       localrunfiles: pOption
902:       output_file: output/pipes_6.out

904:    test:
905:       suffix: 7
906:       nsize: 2
907:       args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX
908:       localrunfiles: pOption
909:       output_file: output/pipes_7.out

911:    test:
912:       suffix: 8
913:       nsize: 2
914:       requires: parmetis
915:       args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type parmetis -options_left no -wash_distribute 1
916:       localrunfiles: pOption
917:       output_file: output/pipes_8.out

919:    test:
920:       suffix: 9
921:       nsize: 2
922:       args: -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -pipe_view
923:       localrunfiles: pOption
924:       output_file: output/pipes_9.out

926:    test:
927:       suffix: 10
928:       nsize: 2
929:       args: -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -isjunc_view
930:       localrunfiles: pOption
931:       output_file: output/pipes_10.out

933: TEST*/