Actual source code: pipes1.c

petsc-3.11.4 2019-09-28
Report Typos and Errors
  1:      static char help[] = "This example demonstrates the use of DMNetwork \n\\n";

  3: /*
  4:   Example: mpiexec -n <np> ./pipes1 -ts_max_steps 10
  5: */

  7: #include "wash.h"
  8:  #include <petscdmnetwork.h>

 10: PetscErrorCode WASHIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void* ctx)
 11: {
 12:   PetscErrorCode    ierr;
 13:   Wash              wash=(Wash)ctx;
 14:   DM                networkdm;
 15:   Vec               localX,localXdot,localF;
 16:   const PetscInt    *cone;
 17:   PetscInt          vfrom,vto,offsetfrom,offsetto,type,varoffset;
 18:   PetscInt          v,vStart,vEnd,e,eStart,eEnd;
 19:   PetscBool         ghost;
 20:   PetscScalar       *farr,*vf,*juncx,*juncf;
 21:   Pipe              pipe;
 22:   PipeField         *pipex,*pipexdot,*pipef;
 23:   DMDALocalInfo     info;
 24:   Junction          junction;
 25:   MPI_Comm          comm;
 26:   PetscMPIInt       rank,size;
 27:   const PetscScalar *xarr,*xdotarr;


 31:   PetscObjectGetComm((PetscObject)ts,&comm);
 32:   MPI_Comm_rank(comm,&rank);
 33:   MPI_Comm_size(comm,&size);

 35:   VecSet(F,0.0);

 37:   localX    = wash->localX;
 38:   localXdot = wash->localXdot;

 40:   TSGetDM(ts,&networkdm);
 41:   DMGetLocalVector(networkdm,&localF);
 42:   VecSet(localF,0.0);

 44:   /* update ghost values of locaX and locaXdot */
 45:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
 46:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

 48:   DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot);
 49:   DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot);

 51:   VecGetArrayRead(localX,&xarr);
 52:   VecGetArrayRead(localXdot,&xdotarr);
 53:   VecGetArray(localF,&farr);

 55:   /* Initialize localF = localX at non-ghost vertices */
 56:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
 57:   for (v=vStart; v<vEnd; v++) {
 58:     DMNetworkIsGhostVertex(networkdm,v,&ghost);
 59:     if (!ghost) {
 60:       DMNetworkGetVariableOffset(networkdm,v,&varoffset);
 61:       juncx  = (PetscScalar*)(xarr+varoffset);
 62:       juncf  = (PetscScalar*)(farr+varoffset);
 63:       juncf[0] = juncx[0];
 64:       juncf[1] = juncx[1];
 65:     }
 66:   }

 68:   /* Edge */
 69:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
 70:   for (e=eStart; e<eEnd; e++) {
 71:     DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe);
 72:     DMNetworkGetVariableOffset(networkdm,e,&varoffset);
 73:     pipex    = (PipeField*)(xarr + varoffset);
 74:     pipexdot = (PipeField*)(xdotarr + varoffset);
 75:     pipef    = (PipeField*)(farr + varoffset);

 77:     /* Get boundary values H0 and QL from connected vertices */
 78:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
 79:     vfrom = cone[0]; /* local ordering */
 80:     vto   = cone[1];
 81:     DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
 82:     DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
 83:     if (pipe->boundary.Q0 == PIPE_CHARACTERISTIC) {
 84:       pipe->boundary.H0 = (xarr+offsetfrom)[1]; /* h_from */
 85:     } else {
 86:       pipe->boundary.Q0 = (xarr+offsetfrom)[0]; /* q_from */
 87:     }
 88:     if (pipe->boundary.HL == PIPE_CHARACTERISTIC) {
 89:       pipe->boundary.QL = (xarr+offsetto)[0];   /* q_to */
 90:     } else {
 91:       pipe->boundary.HL = (xarr+offsetto)[1];   /* h_to */
 92:     }

 94:     /* Evaluate PipeIFunctionLocal() */
 95:     DMDAGetLocalInfo(pipe->da,&info);
 96:     PipeIFunctionLocal(&info, t, pipex, pipexdot, pipef, pipe);

 98:     /* Set F at vfrom */
 99:     vf = (PetscScalar*)(farr+offsetfrom);
100:     if (pipe->boundary.Q0 == PIPE_CHARACTERISTIC) {
101:       vf[0] -= pipex[0].q; /* q_vfrom - q[0] */
102:     } else {
103:       vf[1] -= pipex[0].h; /* h_vfrom - h[0] */
104:     }

106:     /* Set F at vto */
107:     vf = (PetscScalar*)(farr+offsetto);
108:     if (pipe->boundary.HL == PIPE_CHARACTERISTIC) {
109:       vf[1] -= pipex[pipe->nnodes-1].h; /* h_vto - h[last] */
110:     } else {
111:       vf[0] -= pipex[pipe->nnodes-1].q; /* q_vto - q[last] */
112:     }
113:   }

115:   /* Set F at boundary vertices */
116:   for (v=vStart; v<vEnd; v++) {
117:     DMNetworkGetComponent(networkdm,v,0,&type,(void**)&junction);
118:     DMNetworkGetVariableOffset(networkdm,v,&varoffset);
119:     juncf = (PetscScalar *)(farr + varoffset);
120:     if (junction->isEnd == -1) {
121:       juncf[1] -= wash->H0;
122:       } else if (junction->isEnd == 1) {
123:       juncf[0] -= wash->QL;
124:     }
125:   }

127:   VecRestoreArrayRead(localX,&xarr);
128:   VecRestoreArrayRead(localXdot,&xdotarr);
129:   VecRestoreArray(localF,&farr);

131:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
132:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
133:   DMRestoreLocalVector(networkdm,&localF);
134:   /* VecView(F,PETSC_VIEWER_STDOUT_WORLD); */
135:   return(0);
136: }

138: PetscErrorCode WASHSetInitialSolution(DM networkdm,Vec X,Wash wash)
139: {
141:   PetscInt       k,nx,vkey,vfrom,vto,offsetfrom,offsetto;
142:   PetscInt       type,varoffset;
143:   PetscInt       e,eStart,eEnd;
144:   Vec            localX;
145:   PetscScalar    *xarr;
146:   Pipe           pipe;
147:   Junction       junction;
148:   const PetscInt *cone;
149:   const PetscScalar *xarray;

152:   VecSet(X,0.0);
153:   DMGetLocalVector(networkdm,&localX);
154:   VecGetArray(localX,&xarr);

156:   /* Edge */
157:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
158:   for (e=eStart; e<eEnd; e++) {
159:     DMNetworkGetVariableOffset(networkdm,e,&varoffset);
160:     DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe);

162:     /* get from and to vertices */
163:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
164:     vfrom = cone[0]; /* local ordering */
165:     vto   = cone[1];

167:     DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
168:     DMNetworkGetVariableOffset(networkdm,vto,&offsetto);

170:     /* set initial values for this pipe */
171:     /* Q0=0.477432; H0=150.0; needs to be updated by its succeeding pipe. Use SNESSolve()? */
172:     PipeComputeSteadyState(pipe, 0.477432, wash->H0);
173:     VecGetSize(pipe->x,&nx);

175:     VecGetArrayRead(pipe->x,&xarray);
176:     /* copy pipe->x to xarray */
177:     for (k=0; k<nx; k++) {
178:       (xarr+varoffset)[k] = xarray[k];
179:     }

181:     /* set boundary values into vfrom and vto */
182:     if (pipe->boundary.Q0 == PIPE_CHARACTERISTIC) {
183:       (xarr+offsetfrom)[0] += xarray[0];    /* Q0 -> vfrom[0] */
184:     } else {
185:       (xarr+offsetfrom)[1] += xarray[1];    /* H0 -> vfrom[1] */
186:     }

188:     if (pipe->boundary.HL == PIPE_CHARACTERISTIC) {
189:       (xarr+offsetto)[1]   += xarray[nx-1]; /* HL -> vto[1]   */
190:     } else {
191:       (xarr+offsetto)[0]   += xarray[nx-2]; /* QL -> vto[0]   */
192:     }

194:     /* if vform is a head vertex: */
195:     DMNetworkGetComponent(networkdm,vfrom,0,&vkey,(void**)&junction);
196:     if (junction->isEnd == -1) { /* head junction */
197:       (xarr+offsetfrom)[0] = 0.0;      /* 1st Q -- not used */
198:       (xarr+offsetfrom)[1] = wash->H0; /* 1st H */
199:     }

201:     /* if vto is an end vertex: */
202:     DMNetworkGetComponent(networkdm,vto,0,&vkey,(void**)&junction);
203:     if (junction->isEnd == 1) { /* end junction */
204:       (xarr+offsetto)[0] = wash->QL; /* last Q */
205:       (xarr+offsetto)[1] = 0.0;      /* last H -- not used */
206:     }
207:     VecRestoreArrayRead(pipe->x,&xarray);
208:   }

210:   VecRestoreArray(localX,&xarr);
211:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
212:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
213:   DMRestoreLocalVector(networkdm,&localX);

215: #if 0
216:   PetscInt N;
217:   VecGetSize(X,&N);
218:   printf("initial solution %d:\n",N);
219:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
220: #endif
221:   return(0);
222: }

224: PetscErrorCode TSDMNetworkMonitor(TS ts, PetscInt step, PetscReal t, Vec x, void *context)
225: {
226:   PetscErrorCode     ierr;
227:   DMNetworkMonitor   monitor;

230:   monitor = (DMNetworkMonitor)context;
231:   DMNetworkMonitorView(monitor,x);
232:   return(0);
233: }

235: PetscErrorCode PipesView(Vec X,DM networkdm,Wash wash)
236: {
237:   PetscErrorCode       ierr;
238:   Pipe                 pipe;
239:   PetscInt             key,Start,End;
240:   PetscMPIInt          rank;
241:   PetscInt             nx,nnodes,nidx,*idx1,*idx2,*idx1_h,*idx2_h,idx_start,i,k,k1,xstart,j1;
242:   Vec                  Xq,Xh,localX;
243:   IS                   is1_q,is2_q,is1_h,is2_h;
244:   VecScatter           ctx_q,ctx_h;

247:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

253:   VecCreate(PETSC_COMM_WORLD,&Xq);
254:   if (rank == 0) { /* all entries of Xq are in proc[0] */
255:     VecSetSizes(Xq,nx,PETSC_DECIDE);
256:   } else {
257:     VecSetSizes(Xq,0,PETSC_DECIDE);
258:   }
259:   VecSetFromOptions(Xq);
260:   VecSet(Xq,0.0);
261:   VecDuplicate(Xq,&Xh);

263:   DMGetLocalVector(networkdm,&localX);

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

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

270:   VecGetOwnershipRange(X,&xstart,NULL);
271:   k1 = 0;
272:   j1 = 0;
273:   for (i = Start; i < End; i++) {
274:     DMNetworkGetComponent(networkdm,i,0,&key,(void**)&pipe);
275:     nnodes = pipe->nnodes;
276:     idx_start = pipe->id*nnodes;
277:     for (k=0; k<nnodes; k++) {
278:       idx1[k1] = xstart + j1*2*nnodes + 2*k;
279:       idx2[k1] = idx_start + k;

281:       idx1_h[k1] = xstart + j1*2*nnodes + 2*k + 1;
282:       idx2_h[k1] = idx_start + k;
283:       k1++;
284:     }
285:     j1++;
286:   }

288:   ISCreateGeneral(PETSC_COMM_SELF,nidx,idx1,PETSC_COPY_VALUES,&is1_q);
289:   ISCreateGeneral(PETSC_COMM_SELF,nidx,idx2,PETSC_COPY_VALUES,&is2_q);
290:   VecScatterCreate(X,is1_q,Xq,is2_q,&ctx_q);
291:   VecScatterBegin(ctx_q,X,Xq,INSERT_VALUES,SCATTER_FORWARD);
292:   VecScatterEnd(ctx_q,X,Xq,INSERT_VALUES,SCATTER_FORWARD);

294:   ISCreateGeneral(PETSC_COMM_SELF,nidx,idx1_h,PETSC_COPY_VALUES,&is1_h);
295:   ISCreateGeneral(PETSC_COMM_SELF,nidx,idx2_h,PETSC_COPY_VALUES,&is2_h);
296:   VecScatterCreate(X,is1_h,Xh,is2_h,&ctx_h);
297:   VecScatterBegin(ctx_h,X,Xh,INSERT_VALUES,SCATTER_FORWARD);
298:   VecScatterEnd(ctx_h,X,Xh,INSERT_VALUES,SCATTER_FORWARD);

300:   if (!rank) printf("Xq: \n");
301:   VecView(Xq,PETSC_VIEWER_STDOUT_WORLD);
302:   if (!rank) printf("Xh: \n");
303:   VecView(Xh,PETSC_VIEWER_STDOUT_WORLD);

305:   VecScatterDestroy(&ctx_q);
306:   PetscFree4(idx1,idx2,idx1_h,idx2_h);
307:   ISDestroy(&is1_q);
308:   ISDestroy(&is2_q);

310:   VecScatterDestroy(&ctx_h);
311:   ISDestroy(&is1_h);
312:   ISDestroy(&is2_h);

314:   VecDestroy(&Xq);
315:   VecDestroy(&Xh);
316:   DMRestoreLocalVector(networkdm,&localX);
317:   return(0);
318: }

320: PetscErrorCode WashNetworkCleanUp(Wash wash,PetscInt *edgelist)
321: {
323:   PetscMPIInt    rank;

326:   MPI_Comm_rank(wash->comm,&rank);
327:   if (!rank) {
328:     PetscFree(edgelist);
329:     PetscFree2(wash->junction,wash->pipe);
330:   }
331:   return(0);
332: }

334: PetscErrorCode WashNetworkCreate(MPI_Comm comm,PetscInt pipesCase,Wash *wash_ptr,PetscInt **elist)
335: {
337:   PetscInt       nnodes,npipes;
338:   PetscMPIInt    rank;
339:   Wash           wash;
340:   PetscInt       i,numVertices,numEdges;
341:   PetscInt       *edgelist;
342:   Junction       junctions=NULL;
343:   Pipe           pipes=NULL;

346:   MPI_Comm_rank(comm,&rank);

348:   PetscCalloc1(1,&wash);
349:   wash->comm = comm;
350:   *wash_ptr  = wash;
351:   wash->Q0   = 0.477432; /* copied from initial soluiton */
352:   wash->H0   = 150.0;
353:   wash->HL   = 143.488; /* copied from initial soluiton */
354:   wash->nnodes_loc = 0;

356:   numVertices = 0;
357:   numEdges    = 0;
358:   edgelist    = NULL;

360:   if (!rank) {
361:     PetscPrintf(PETSC_COMM_SELF,"Setup pipesCase %D\n",pipesCase);
362:   }
363:   nnodes = 6;
364:   PetscOptionsGetInt(NULL,NULL, "-npipenodes", &nnodes, NULL);

366:   /* Set global number of pipes, edges, and junctions */
367:   /*-------------------------------------------------*/
368:   switch (pipesCase) {
369:   case 0:
370:     /* pipeCase 0: */
371:     /* =============================
372:     v0 --E0--> v1--E1--> v2 --E2-->v3
373:     ================================  */
374:     npipes = 3;
375:     PetscOptionsGetInt(NULL,NULL, "-npipes", &npipes, NULL);
376:     wash->nedge   = npipes;
377:     wash->nvertex = npipes + 1;

379:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
380:     numVertices = 0;
381:     numEdges    = 0;
382:     edgelist    = NULL;
383:     if (!rank) {
384:       numVertices = wash->nvertex;
385:       numEdges    = wash->nedge;

387:       PetscCalloc1(2*numEdges,&edgelist);
388:       for (i=0; i<numEdges; i++) {
389:         edgelist[2*i] = i; edgelist[2*i+1] = i+1;
390:       }

392:       /* Add network components */
393:       /*------------------------*/
394:       PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
395:       /* vertex */
396:       for (i=0; i<numVertices; i++) {
397:         junctions[i].id = i;
398:         junctions[i].isEnd = 0;
399:         junctions[i].nedges_in = 1; junctions[i].nedges_out = 1;

401:         /* Set GPS data */
402:         junctions[i].latitude  = 0.0;
403:         junctions[i].longitude = 0.0;
404:       }
405:       junctions[0].isEnd                  = -1;
406:       junctions[0].nedges_in              =  0;
407:       junctions[numVertices-1].isEnd      =  1;
408:       junctions[numVertices-1].nedges_out =  0;

410:       /* edge and pipe */
411:       for (i=0; i<numEdges; i++) {
412:         pipes[i].id   = i;
413:         pipes[i].nnodes = nnodes;
414:       }
415:     }
416:     break;
417:   case 1:
418:     /* pipeCase 1: */
419:     /* ==========================
420:                 v2
421:                 ^
422:                 |
423:                E2
424:                 |
425:     v0 --E0--> v3--E1--> v1
426:     =============================  */
427:     npipes = 3;
428:     wash->nedge   = npipes;
429:     wash->nvertex = npipes + 1;

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

436:       PetscCalloc1(2*numEdges,&edgelist);
437:       edgelist[0] = 0; edgelist[1] = 3;  /* edge[0] */
438:       edgelist[2] = 3; edgelist[3] = 1;  /* edge[1] */
439:       edgelist[4] = 3; edgelist[5] = 2;  /* edge[2] */

441:       /* Add network components */
442:       /*------------------------*/
443:       PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
444:       /* vertex */
445:       for (i=0; i<numVertices; i++) {
446:         junctions[i].id = i;

448:         /* Set GPS data */
449:         junctions[i].latitude  = 0.0;
450:         junctions[i].longitude = 0.0;
451:       }
452:       junctions[0].isEnd = -1; junctions[0].nedges_in = 0; junctions[0].nedges_out = 1;
453:       junctions[1].isEnd =  1; junctions[1].nedges_in = 1; junctions[1].nedges_out = 0;
454:       junctions[2].isEnd =  1; junctions[2].nedges_in = 1; junctions[2].nedges_out = 0;
455:       junctions[3].isEnd =  0; junctions[3].nedges_in = 1; junctions[3].nedges_out = 2;

457:       /* edge and pipe */
458:       for (i=0; i<numEdges; i++) {
459:         pipes[i].id     = i;
460:         pipes[i].nnodes = nnodes;
461:       }
462:     }
463:     break;
464:   case 2:
465:     /* pipeCase 2: */
466:     /* ==========================
467:          v2--> E2
468:                 |
469:     v0 --E0--> v3--E1--> v1
470:     =============================  */

472:     /* Set application parameters -- to be used in function evalutions */
473:     npipes = 3;
474:     wash->nedge   = npipes;
475:     wash->nvertex = npipes + 1;

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

482:       PetscCalloc1(2*numEdges,&edgelist);
483:       edgelist[0] = 0; edgelist[1] = 3;  /* edge[0] */
484:       edgelist[2] = 3; edgelist[3] = 1;  /* edge[1] */
485:       edgelist[4] = 2; edgelist[5] = 3;  /* edge[2] */

487:       /* Add network components */
488:       /*------------------------*/
489:       PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
490:       /* vertex */
491:       for (i=0; i<numVertices; i++) {
492:         junctions[i].id = i;

494:         /* Set GPS data */
495:         junctions[i].latitude  = 0.0;
496:         junctions[i].longitude = 0.0;
497:       }
498:       junctions[0].isEnd = -1; junctions[0].nedges_in = 0; junctions[0].nedges_out = 1;
499:       junctions[1].isEnd =  1; junctions[1].nedges_in = 1; junctions[1].nedges_out = 0;
500:       junctions[2].isEnd = -1; junctions[2].nedges_in = 0; junctions[2].nedges_out = 1;
501:       junctions[3].isEnd =  0; junctions[3].nedges_in = 2; junctions[3].nedges_out = 1;

503:       /* edge and pipe */
504:       for (i=0; i<numEdges; i++) {
505:         pipes[i].id     = i;
506:         pipes[i].nnodes = nnodes;
507:       }
508:     }
509:     break;
510:   default:
511:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"not done yet");
512:   }

514:   *wash_ptr      = wash;
515:   wash->nedge    = numEdges;
516:   wash->nvertex  = numVertices;
517:   *elist         = edgelist;
518:   wash->junction = junctions;
519:   wash->pipe     = pipes;
520:   return(0);
521: }

523: extern PetscErrorCode PipeCreateJacobian(Pipe,Mat*,Mat*[]);
524: extern PetscErrorCode PipeDestroyJacobian(Pipe);
525: extern PetscErrorCode JunctionCreateJacobian(DM,PetscInt,Mat*,Mat*[]);
526: extern PetscErrorCode JunctionDestroyJacobian(DM,PetscInt,Junction);
527: /* ------------------------------------------------------- */
528: int main(int argc,char ** argv)
529: {
530:   PetscErrorCode    ierr;
531:   Wash              wash;
532:   Junction          junctions,junction;
533:   Pipe              pipe,pipes;
534:   PetscInt          numEdges,numVertices,KeyPipe,KeyJunction;
535:   PetscInt          *edgelist = NULL,*edgelists[1];
536:   PetscInt          i,e,v,eStart,eEnd,vStart,vEnd,key,frombType,tobType;
537:   PetscInt          vfrom,vto,vkey,type,varoffset;
538:   PetscInt          from_nedge_in,from_nedge_out,to_nedge_in;
539:   const PetscInt    *cone;
540:   DM                networkdm;
541:   PetscMPIInt       size,rank;
542:   PetscReal         ftime = 20.0;
543:   Vec               X;
544:   TS                ts;
545:   PetscInt          steps;
546:   TSConvergedReason reason;
547:   PetscBool         viewpipes,monipipes=PETSC_FALSE,userJac=PETSC_TRUE;
548:   PetscInt          pipesCase;
549:   DMNetworkMonitor  monitor;

551:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
552:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
553:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

555:   /* Create and setup network */
556:   /*--------------------------*/
557:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
558:   PetscOptionsGetBool(NULL,NULL,"-user_Jac",&userJac,NULL);
559:   PetscOptionsGetBool(NULL,NULL,"-pipe_monitor",&monipipes,NULL);
560:   if (size == 1 && monipipes) {
561:     DMNetworkMonitorCreate(networkdm,&monitor);
562:   }
563:   /* Register the components in the network */
564:   DMNetworkRegisterComponent(networkdm,"junctionstruct",sizeof(struct _p_Junction),&KeyJunction);
565:   DMNetworkRegisterComponent(networkdm,"pipestruct",sizeof(struct _p_Pipe),&KeyPipe);

567:   /* Set global number of pipes, edges, and vertices */
568:   pipesCase = 2;
569:   PetscOptionsGetInt(NULL,NULL, "-case", &pipesCase, NULL);

571:   WashNetworkCreate(PETSC_COMM_WORLD,pipesCase,&wash,&edgelist);
572:   numEdges    = wash->nedge;
573:   numVertices = wash->nvertex;
574:   junctions    = wash->junction;
575:   pipes       = wash->pipe;

577:   /* Set number of vertices and edges */
578:   DMNetworkSetSizes(networkdm,1,0,&numVertices,&numEdges,NULL,NULL);
579:   /* Add edge connectivity */
580:   edgelists[0] = edgelist;
581:   DMNetworkSetEdgeList(networkdm,edgelists,NULL);
582:   /* Set up the network layout */
583:   DMNetworkLayoutSetUp(networkdm);

585:   /* Add EDGEDATA component to all edges -- currently networkdm is a sequential network */
586:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
587:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

589:   for (e = eStart; e < eEnd; e++) {
590:     /* Add Pipe component to all edges -- create pipe here */
591:     DMNetworkAddComponent(networkdm,e,KeyPipe,&pipes[e-eStart]);

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

596:     if (size == 1 && monipipes) { /* Add monitor -- show Q_{pipes[e-eStart].id}? */
597:       DMNetworkMonitorAdd(monitor, "Pipe Q", e, pipes[e-eStart].nnodes, 0, 2, -0.8, 0.8, PETSC_TRUE);
598:       DMNetworkMonitorAdd(monitor, "Pipe H", e, pipes[e-eStart].nnodes, 1, 2, -400.0, 800.0, PETSC_TRUE);
599:     }
600:   }

602:   /* Add Junction component to all vertices */
603:   for (v = vStart; v < vEnd; v++) {
604:     DMNetworkAddComponent(networkdm,v,KeyJunction,&junctions[v-vStart]);

606:     /* Add number of variables to vertex */
607:     DMNetworkAddNumVariables(networkdm,v,2);
608:   }

610:   /* Set up DM for use */
611:   DMSetUp(networkdm);
612:   WashNetworkCleanUp(wash,edgelist);

614:   /* Network partitioning and distribution of data */
615:   DMNetworkDistribute(&networkdm,0);

617:   /* PipeSetUp -- each process only sets its own pipes */
618:   /*---------------------------------------------------*/
619:   DMNetworkHasJacobian(networkdm,userJac,userJac);
620:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
621:   for (e=eStart; e<eEnd; e++) { /* each edge has only one component, pipe */
622:     DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe);
623:     DMNetworkGetVariableOffset(networkdm,e,&varoffset);

625:     /* Setup conntected vertices */
626:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
627:     vfrom = cone[0]; /* local ordering */
628:     vto   = cone[1];

630:     /* vfrom */
631:     DMNetworkGetComponent(networkdm,vfrom,0,&vkey,(void**)&junction);
632:     from_nedge_in  = junction->nedges_in;
633:     from_nedge_out = junction->nedges_out;

635:     /* vto */
636:     DMNetworkGetComponent(networkdm,vto,0,&vkey,(void**)&junction);
637:     to_nedge_in = junction->nedges_in;

639:     pipe->comm = PETSC_COMM_SELF; /* must be set here, otherwise crashes in my mac??? */
640:     wash->nnodes_loc += pipe->nnodes; /* local total num of nodes, will be used by PipesView() */
641:     PipeSetParameters(pipe,
642:                              600.0,          /* length */
643:                              pipe->nnodes,   /* nnodes -- rm from PipeSetParameters */
644:                              0.5,            /* diameter */
645:                              1200.0,         /* a */
646:                              0.018);    /* friction */

648:     /* set boundary conditions for this pipe */
649:     if (from_nedge_in <= 1 && from_nedge_out > 0) {
650:       frombType = 0;
651:     } else {
652:       frombType = 1;
653:     }

655:     if (to_nedge_in == 1) {
656:       tobType = 0;
657:     } else {
658:       tobType = 1;
659:     }

661:     if (frombType == 0) {
662:       pipe->boundary.Q0 = PIPE_CHARACTERISTIC; /* will be obtained from characteristic */
663:       pipe->boundary.H0 = wash->H0;
664:     } else {
665:       pipe->boundary.Q0 = wash->Q0;
666:       pipe->boundary.H0 = PIPE_CHARACTERISTIC; /* will be obtained from characteristic */
667:     }
668:     if (tobType == 0) {
669:       pipe->boundary.QL = wash->QL;
670:       pipe->boundary.HL = PIPE_CHARACTERISTIC; /* will be obtained from characteristic */
671:     } else {
672:       pipe->boundary.QL = PIPE_CHARACTERISTIC; /* will be obtained from characteristic */
673:       pipe->boundary.HL = wash->HL;
674:     }

676:     PipeSetUp(pipe);

678:     if (userJac) {
679:       /* Create Jacobian matrix structures for a Pipe */
680:       Mat            *J;
681:       PipeCreateJacobian(pipe,NULL,&J);
682:       DMNetworkEdgeSetMatrix(networkdm,e,J);
683:     }
684:   }

686:   if (userJac) {
687:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
688:     for (v=vStart; v<vEnd; v++) {
689:       Mat            *J;
690:       JunctionCreateJacobian(networkdm,v,NULL,&J);
691:       DMNetworkVertexSetMatrix(networkdm,v,J);

693:       DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction);
694:       junction->jacobian = J;
695:     }
696:   }

698:   /* create vectors */
699:   DMCreateGlobalVector(networkdm,&X);
700:   DMCreateLocalVector(networkdm,&wash->localX);
701:   DMCreateLocalVector(networkdm,&wash->localXdot);

703:   /* Setup solver                                           */
704:   /*--------------------------------------------------------*/
705:   TSCreate(PETSC_COMM_WORLD,&ts);

707:   TSSetDM(ts,(DM)networkdm);
708:   TSSetIFunction(ts,NULL,WASHIFunction,wash);

710:   TSSetMaxTime(ts,ftime);
711:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
712:   TSSetTimeStep(ts,0.1);
713:   TSSetType(ts,TSBEULER);
714:   if (size == 1 && monipipes) {
715:     TSMonitorSet(ts, TSDMNetworkMonitor, monitor, NULL);
716:   }
717:   TSSetFromOptions(ts);

719:   WASHSetInitialSolution(networkdm,X,wash);

721:   TSSolve(ts,X);

723:   TSGetSolveTime(ts,&ftime);
724:   TSGetStepNumber(ts,&steps);
725:   TSGetConvergedReason(ts,&reason);
726:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
727:   /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */

729:   viewpipes = PETSC_FALSE;
730:   PetscOptionsGetBool(NULL,NULL, "-Jac_view", &viewpipes,NULL);
731:   if (viewpipes) {
732:     SNES snes;
733:     Mat  Jac;
734:     TSGetSNES(ts,&snes);
735:     SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
736:     MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
737:   }

739:   /* View solution q and h */
740:   /* --------------------- */
741:   viewpipes = PETSC_FALSE;
742:   PetscOptionsGetBool(NULL,NULL, "-pipe_view", &viewpipes,NULL);
743:   if (viewpipes) {
744:     PipesView(X,networkdm,wash);
745:   }

747:   /* Free spaces */
748:   /* ----------- */
749:   TSDestroy(&ts);
750:   VecDestroy(&X);
751:   VecDestroy(&wash->localX);
752:   VecDestroy(&wash->localXdot);

754:   /* Destroy objects from each pipe that are created in PipeSetUp() */
755:   DMNetworkGetEdgeRange(networkdm,&eStart, &eEnd);
756:   for (i = eStart; i < eEnd; i++) {
757:     DMNetworkGetComponent(networkdm,i,0,&key,(void**)&pipe);
758:     DMDestroy(&(pipe->da));
759:     VecDestroy(&pipe->x);
760:     PipeDestroyJacobian(pipe);
761:   }
762:   if (userJac) {
763:     for (v=vStart; v<vEnd; v++) {
764:       DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction);
765:       JunctionDestroyJacobian(networkdm,v,junction);
766:     }
767:   }

769:   if (size == 1 && monipipes) {
770:     DMNetworkMonitorDestroy(&monitor);
771:   }
772:   DMDestroy(&networkdm);
773:   PetscFree(wash);
774:   PetscFinalize();
775:   return ierr;
776: }

778: /*
779:     PipeCreateJacobian - Create Jacobian matrix structures for a Pipe.

781:     Collective on Pipe

783:     Input Parameter:
784: +   pipe - the Pipe object
785: -   Jin - array of three constructed Jacobian matrices to be reused. Set NULL if it is not available

787:     Output Parameter:
788: .   J  - array of three empty Jacobian matrices

790:     Level: beginner
791: */
792: PetscErrorCode PipeCreateJacobian(Pipe pipe,Mat *Jin,Mat *J[])
793: {
795:   Mat            *Jpipe;
796:   PetscInt       i,M,rows[2],cols[2],*nz;
797:   PetscScalar    *aa;

800:   if (Jin) {
801:     *J = Jin;
802:     pipe->jacobian = Jin;
803:     PetscObjectReference((PetscObject)(Jin[0]));
804:     return(0);
805:   }
806:   PetscMalloc1(3,&Jpipe);

808:   /* Jacobian for this pipe */
809:   DMSetMatrixStructureOnly(pipe->da,PETSC_TRUE);
810:   DMCreateMatrix(pipe->da,&Jpipe[0]);
811:   DMSetMatrixStructureOnly(pipe->da,PETSC_FALSE);

813:   /* Jacobian for upstream vertex */
814:   MatGetSize(Jpipe[0],&M,NULL);
815:   PetscCalloc2(M,&nz,4,&aa);
816:   for (i=0; i<4; i++) aa[i] = 0.0;

818:   MatCreate(PETSC_COMM_SELF,&Jpipe[1]);
819:   MatSetSizes(Jpipe[1],PETSC_DECIDE,PETSC_DECIDE,M,2);
820:   MatSetFromOptions(Jpipe[1]);
821:   MatSetOption(Jpipe[1],MAT_STRUCTURE_ONLY,PETSC_TRUE);
822:   nz[0] = 2; nz[1] = 2;
823:   rows[0] = 0; rows[1] = 1;
824:   cols[0] = 0; cols[1] = 1;
825:   MatSeqAIJSetPreallocation(Jpipe[1],0,nz);
826:   MatSetValues(Jpipe[1],2,rows,2,cols,aa,INSERT_VALUES);
827:   MatAssemblyBegin(Jpipe[1],MAT_FINAL_ASSEMBLY);
828:   MatAssemblyEnd(Jpipe[1],MAT_FINAL_ASSEMBLY);

830:   /* Jacobian for downstream vertex */
831:   MatCreate(PETSC_COMM_SELF,&Jpipe[2]);
832:   MatSetSizes(Jpipe[2],PETSC_DECIDE,PETSC_DECIDE,M,2);
833:   MatSetFromOptions(Jpipe[2]);
834:   MatSetOption(Jpipe[2],MAT_STRUCTURE_ONLY,PETSC_TRUE);
835:   nz[0] = 0; nz[1] = 0; nz[M-2] = 2; nz[M-1] = 2;
836:   rows[0] = M - 2; rows[1] = M - 1;
837:   MatSeqAIJSetPreallocation(Jpipe[2],0,nz);
838:   MatSetValues(Jpipe[2],2,rows,2,cols,aa,INSERT_VALUES);
839:   MatAssemblyBegin(Jpipe[2],MAT_FINAL_ASSEMBLY);
840:   MatAssemblyEnd(Jpipe[2],MAT_FINAL_ASSEMBLY);

842:   PetscFree2(nz,aa);

844:   *J = Jpipe;
845:   pipe->jacobian = Jpipe;
846:   return(0);
847: }

849: PetscErrorCode PipeDestroyJacobian(Pipe pipe)
850: {
852:   Mat            *Jpipe = pipe->jacobian;
853:   PetscInt       i;

856:   if (Jpipe) {
857:     for (i=0; i<3; i++) {
858:       MatDestroy(&Jpipe[i]);
859:     }
860:   }
861:   PetscFree(Jpipe);
862:   return(0);
863: }

865: /*
866:     JunctionCreateJacobian - Create Jacobian matrices for a vertex.

868:     Collective on Pipe

870:     Input Parameter:
871: +   dm - the DMNetwork object
872: .   v - vertex point
873: -   Jin - Jacobian patterns created by JunctionCreateJacobianSample() for reuse

875:     Output Parameter:
876: .   J  - array of Jacobian matrices (see dmnetworkimpl.h)

878:     Level: beginner
879: */
880: PetscErrorCode JunctionCreateJacobian(DM dm,PetscInt v,Mat *Jin,Mat *J[])
881: {
883:   Mat            *Jv;
884:   PetscInt       nedges,e,i,M,N,*rows,*cols;
885:   PetscBool      isSelf;
886:   const PetscInt *edges,*cone;
887:   PetscScalar    *zeros;

890:   /* Get arrary size of Jv */
891:   DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
892:   if (nedges <= 0) SETERRQ2(PETSC_COMM_SELF,1,"%d vertex, nedges %d\n",v,nedges);

894:   /* two Jacobians for each connected edge: J(v,e) and J(v,vc); adding J(v,v), total 2*nedges+1 Jacobians */
895:   PetscCalloc1(2*nedges+1,&Jv);

897:   /* Create dense zero block for this vertex: J[0] = Jacobian(v,v) */
898:   DMNetworkGetNumVariables(dm,v,&M);
899:   if (M !=2) SETERRQ1(PETSC_COMM_SELF,1,"M != 2",M);
900:   PetscMalloc3(M,&rows,M,&cols,M*M,&zeros);
901:   PetscMemzero(zeros,M*M*sizeof(PetscScalar));
902:   for (i=0; i<M; i++) rows[i] = i;

904:   for (e=0; e<nedges; e++) {
905:     /* create Jv[2*e+1] = Jacobian(v,e), e: supporting edge */
906:     DMNetworkGetConnectedVertices(dm,edges[e],&cone);
907:     isSelf = (v == cone[0]) ? PETSC_TRUE:PETSC_FALSE;

909:     if (Jin) {
910:       if (isSelf) {
911:         Jv[2*e+1] = Jin[0];
912:       } else {
913:         Jv[2*e+1] = Jin[1];
914:       }
915:       Jv[2*e+2] = Jin[2];
916:       PetscObjectReference((PetscObject)(Jv[2*e+1]));
917:       PetscObjectReference((PetscObject)(Jv[2*e+2]));
918:     } else {
919:       /* create J(v,e) */
920:       MatCreate(PETSC_COMM_SELF,&Jv[2*e+1]);
921:       DMNetworkGetNumVariables(dm,edges[e],&N);
922:       MatSetSizes(Jv[2*e+1],PETSC_DECIDE,PETSC_DECIDE,M,N);
923:       MatSetFromOptions(Jv[2*e+1]);
924:       MatSetOption(Jv[2*e+1],MAT_STRUCTURE_ONLY,PETSC_TRUE);
925:       MatSeqAIJSetPreallocation(Jv[2*e+1],2,NULL);
926:       if (N) {
927:         if (isSelf) { /* coupling at upstream */
928:           for (i=0; i<2; i++) cols[i] = i;
929:         } else { /* coupling at downstream */
930:           cols[0] = N-2; cols[1] = N-1;
931:         }
932:         MatSetValues(Jv[2*e+1],2,rows,2,cols,zeros,INSERT_VALUES);
933:       }
934:       MatAssemblyBegin(Jv[2*e+1],MAT_FINAL_ASSEMBLY);
935:       MatAssemblyEnd(Jv[2*e+1],MAT_FINAL_ASSEMBLY);

937:       /* create Jv[2*e+2] = Jacobian(v,vc), vc: connected vertex.
938:        In WashNetwork, v and vc are not connected, thus Jacobian(v,vc) is empty */
939:       MatCreate(PETSC_COMM_SELF,&Jv[2*e+2]);
940:       MatSetSizes(Jv[2*e+2],PETSC_DECIDE,PETSC_DECIDE,M,M); /* empty matrix, sizes can be arbitrary */
941:       MatSetFromOptions(Jv[2*e+2]);
942:       MatSetOption(Jv[2*e+2],MAT_STRUCTURE_ONLY,PETSC_TRUE);
943:       MatSeqAIJSetPreallocation(Jv[2*e+2],1,NULL);
944:       MatAssemblyBegin(Jv[2*e+2],MAT_FINAL_ASSEMBLY);
945:       MatAssemblyEnd(Jv[2*e+2],MAT_FINAL_ASSEMBLY);
946:     }
947:   }
948:   PetscFree3(rows,cols,zeros);

950:   *J = Jv;
951:   return(0);
952: }

954: PetscErrorCode JunctionDestroyJacobian(DM dm,PetscInt v,Junction junc)
955: {
957:   Mat            *Jv=junc->jacobian;
958:   const PetscInt *edges;
959:   PetscInt       nedges,e;

962:   if (!Jv) return(0);

964:   DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
965:   for (e=0; e<nedges; e++) {
966:     MatDestroy(&Jv[2*e+1]);
967:     MatDestroy(&Jv[2*e+2]);
968:   }
969:   PetscFree(Jv);
970:   return(0);
971: }

973: /*TEST

975:    build:
976:      depends: pipeInterface.c pipeImpls.c

978:    test:
979:       args: -ts_monitor -case 1 -ts_max_steps 1

981:    test:
982:       suffix: 2
983:       nsize: 2
984:       args: -ts_monitor -case 1 -ts_max_steps 1 -petscpartitioner_type simple
985:       output_file: output/pipes1_1.out

987:    test:
988:       suffix: 3
989:       nsize: 4
990:       args: -ts_monitor -case 1 -ts_max_steps 1 -petscpartitioner_type simple
991:       output_file: output/pipes1_1.out

993:    test:
994:       suffix: 4
995:       args: -ts_monitor -case 0 -ts_max_steps 1

997:    test:
998:       suffix: 5
999:       nsize: 2
1000:       args: -ts_monitor -case 0 -ts_max_steps 1 -petscpartitioner_type simple
1001:       output_file: output/pipes1_4.out

1003:    test:
1004:       suffix: 6
1005:       nsize: 4
1006:       args: -ts_monitor -case 0 -ts_max_steps 1 -petscpartitioner_type simple
1007:       output_file: output/pipes1_4.out

1009:    test:
1010:       suffix: 7
1011:       args: -ts_monitor -case 2 -ts_max_steps 1

1013:    test:
1014:       suffix: 8
1015:       nsize: 2
1016:       args: -ts_monitor -case 2 -ts_max_steps 1 -petscpartitioner_type simple
1017:       output_file: output/pipes1_7.out

1019:    test:
1020:       suffix: 9
1021:       nsize: 4
1022:       args: -ts_monitor -case 2 -ts_max_steps 1 -petscpartitioner_type simple
1023:       output_file: output/pipes1_7.out

1025: TEST*/