Actual source code: pipes1.c
petsc-3.11.1 2019-04-12
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*/