Actual source code: pipes1.c
petsc-3.12.0 2019-09-29
1: static char help[] = "This example demonstrates the use of DMNetwork \n\\n";
2: /*
3: Example: mpiexec -n <np> ./pipes1 -ts_max_steps 10
4: */
6: #include "wash.h"
7: #include <petscdmplex.h>
9: /*
10: WashNetworkDistribute - proc[0] distributes sequential wash object
11: Input Parameters:
12: . comm - MPI communicator
13: . wash - wash context with all network data in proc[0]
15: Output Parameter:
16: . wash - wash context with nedge, nvertex and edgelist distributed
17: */
18: PetscErrorCode WashNetworkDistribute(MPI_Comm comm,Wash wash)
19: {
21: PetscMPIInt rank,size,tag=0;
22: PetscInt i,e,v,numEdges,numVertices,nedges,*eowners=NULL,estart,eend,*vtype=NULL,nvertices;
23: PetscInt *edgelist = wash->edgelist,*nvtx=NULL,*vtxDone=NULL;
26: MPI_Comm_rank(comm,&rank);
27: MPI_Comm_size(comm,&size);
29: numEdges = wash->nedge;
30: numVertices = wash->nvertex;
32: /* (1) all processes get global and local number of edges */
33: MPI_Bcast(&numEdges,1,MPIU_INT,0,comm);
34: nedges = numEdges/size; /* local nedges */
35: if (!rank) {
36: nedges += numEdges - size*(numEdges/size);
37: }
38: wash->Nedge = numEdges;
39: wash->nedge = nedges;
40: /* PetscPrintf(PETSC_COMM_SELF,"[%d] nedges %d, numEdges %d\n",rank,nedges,numEdges); */
42: PetscCalloc3(size+1,&eowners,size,&nvtx,numVertices,&vtxDone);
43: MPI_Allgather(&nedges,1,MPIU_INT,eowners+1,1,MPIU_INT,PETSC_COMM_WORLD);
44: eowners[0] = 0;
45: for (i=2; i<=size; i++) {
46: eowners[i] += eowners[i-1];
47: }
49: estart = eowners[rank];
50: eend = eowners[rank+1];
51: /* PetscPrintf(PETSC_COMM_SELF,"[%d] own lists row %d - %d\n",rank,estart,eend); */
53: /* (2) distribute row block edgelist to all processors */
54: if (!rank) {
55: vtype = wash->vtype;
56: for (i=1; i<size; i++) {
57: /* proc[0] sends edgelist to proc[i] */
58: MPI_Send(edgelist+2*eowners[i],2*(eowners[i+1]-eowners[i]),MPIU_INT,i,tag,comm);
60: /* proc[0] sends vtype to proc[i] */
61: MPI_Send(vtype+2*eowners[i],2*(eowners[i+1]-eowners[i]),MPIU_INT,i,tag,comm);
62: }
63: } else {
64: MPI_Status status;
65: PetscMalloc1(2*(eend-estart),&vtype);
66: PetscMalloc1(2*(eend-estart),&edgelist);
68: MPI_Recv(edgelist,2*(eend-estart),MPIU_INT,0,tag,comm,&status);
69: MPI_Recv(vtype,2*(eend-estart),MPIU_INT,0,tag,comm,&status);
70: }
72: wash->edgelist = edgelist;
74: /* (3) all processes get global and local number of vertices, without ghost vertices */
75: if (!rank) {
76: for (i=0; i<size; i++) {
77: for (e=eowners[i]; e<eowners[i+1]; e++) {
78: v = edgelist[2*e];
79: if (!vtxDone[v]) {
80: nvtx[i]++; vtxDone[v] = 1;
81: }
82: v = edgelist[2*e+1];
83: if (!vtxDone[v]) {
84: nvtx[i]++; vtxDone[v] = 1;
85: }
86: }
87: }
88: }
89: MPI_Bcast(&numVertices,1,MPIU_INT,0,PETSC_COMM_WORLD);
90: MPI_Scatter(nvtx,1,MPIU_INT,&nvertices,1,MPIU_INT,0,PETSC_COMM_WORLD);
91: PetscFree3(eowners,nvtx,vtxDone);
93: wash->Nvertex = numVertices;
94: wash->nvertex = nvertices;
95: wash->vtype = vtype;
96: return(0);
97: }
99: PetscErrorCode WASHIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void* ctx)
100: {
102: Wash wash=(Wash)ctx;
103: DM networkdm;
104: Vec localX,localXdot,localF, localXold;
105: const PetscInt *cone;
106: PetscInt vfrom,vto,offsetfrom,offsetto,varoffset;
107: PetscInt v,vStart,vEnd,e,eStart,eEnd;
108: PetscInt nend,type;
109: PetscBool ghost;
110: PetscScalar *farr,*juncf, *pipef;
111: PetscReal dt;
112: Pipe pipe;
113: PipeField *pipex,*pipexdot,*juncx;
114: Junction junction;
115: DMDALocalInfo info;
116: const PetscScalar *xarr,*xdotarr, *xoldarr;
119: localX = wash->localX;
120: localXdot = wash->localXdot;
122: TSGetSolution(ts,&localXold);
123: TSGetDM(ts,&networkdm);
124: TSGetTimeStep(ts,&dt);
125: DMGetLocalVector(networkdm,&localF);
127: /* Set F and localF as zero */
128: VecSet(F,0.0);
129: VecSet(localF,0.0);
131: /* Update ghost values of locaX and locaXdot */
132: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
133: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
135: DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot);
136: DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot);
138: VecGetArrayRead(localX,&xarr);
139: VecGetArrayRead(localXdot,&xdotarr);
140: VecGetArrayRead(localXold,&xoldarr);
141: VecGetArray(localF,&farr);
143: /* junction->type == JUNCTION:
144: juncf[0] = -qJ + sum(qin); juncf[1] = qJ - sum(qout)
145: junction->type == RESERVOIR (upper stream):
146: juncf[0] = -hJ + H0; juncf[1] = qJ - sum(qout)
147: junction->type == VALVE (down stream):
148: juncf[0] = -qJ + sum(qin); juncf[1] = qJ
149: */
150: /* Vertex/junction initialization */
151: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
152: for (v=vStart; v<vEnd; v++) {
153: DMNetworkIsGhostVertex(networkdm,v,&ghost);
154: if (ghost) continue;
156: DMNetworkGetComponent(networkdm,v,0,&type,(void**)&junction);
157: DMNetworkGetVariableOffset(networkdm,v,&varoffset);
158: juncx = (PipeField*)(xarr + varoffset);
159: juncf = (PetscScalar*)(farr + varoffset);
161: juncf[0] = -juncx[0].q;
162: juncf[1] = juncx[0].q;
164: if (junction->type == RESERVOIR) { /* upstream reservoir */
165: juncf[0] = juncx[0].h - wash->H0;
166: }
167: }
169: /* Edge/pipe */
170: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
171: for (e=eStart; e<eEnd; e++) {
172: DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe);
173: DMNetworkGetVariableOffset(networkdm,e,&varoffset);
174: pipex = (PipeField*)(xarr + varoffset);
175: pipexdot = (PipeField*)(xdotarr + varoffset);
176: pipef = (PetscScalar*)(farr + varoffset);
178: /* Get some data into the pipe structure: note, some of these operations
179: * might be redundant. Will it consume too much time? */
180: pipe->dt = dt;
181: pipe->xold = (PipeField*)(xoldarr + varoffset);
183: /* Evaluate F over this edge/pipe: pipef[1], ...,pipef[2*nend] */
184: DMDAGetLocalInfo(pipe->da,&info);
185: PipeIFunctionLocal_Lax(&info,t,pipex,pipexdot,pipef,pipe);
187: /* Get boundary values from connected vertices */
188: DMNetworkGetConnectedVertices(networkdm,e,&cone);
189: vfrom = cone[0]; /* local ordering */
190: vto = cone[1];
191: DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
192: DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
194: /* Evaluate upstream boundary */
195: DMNetworkGetComponent(networkdm,vfrom,0,&type,(void**)&junction);
196: if (junction->type != JUNCTION && junction->type != RESERVOIR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"junction type is not supported");
197: juncx = (PipeField*)(xarr + offsetfrom);
198: juncf = (PetscScalar*)(farr + offsetfrom);
200: pipef[0] = pipex[0].h - juncx[0].h;
201: juncf[1] -= pipex[0].q;
203: /* Evaluate downstream boundary */
204: DMNetworkGetComponent(networkdm,vto,0,&type,(void**)&junction);
205: if (junction->type != JUNCTION && junction->type != VALVE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"junction type is not supported");
206: juncx = (PipeField*)(xarr + offsetto);
207: juncf = (PetscScalar*)(farr + offsetto);
208: nend = pipe->nnodes - 1;
210: pipef[2*nend + 1] = pipex[nend].h - juncx[0].h;
211: juncf[0] += pipex[nend].q;
212: }
214: VecRestoreArrayRead(localX,&xarr);
215: VecRestoreArrayRead(localXdot,&xdotarr);
216: VecRestoreArray(localF,&farr);
218: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
219: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
220: DMRestoreLocalVector(networkdm,&localF);
221: /*
222: PetscPrintf(PETSC_COMM_WORLD("F:\n");
223: VecView(F,PETSC_VIEWER_STDOUT_WORLD);
224: */
225: return(0);
226: }
228: PetscErrorCode WASHSetInitialSolution(DM networkdm,Vec X,Wash wash)
229: {
231: PetscInt k,nx,vkey,vfrom,vto,offsetfrom,offsetto;
232: PetscInt type,varoffset;
233: PetscInt e,eStart,eEnd;
234: Vec localX;
235: PetscScalar *xarr;
236: Pipe pipe;
237: Junction junction;
238: const PetscInt *cone;
239: const PetscScalar *xarray;
242: VecSet(X,0.0);
243: DMGetLocalVector(networkdm,&localX);
244: VecGetArray(localX,&xarr);
246: /* Edge */
247: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
248: for (e=eStart; e<eEnd; e++) {
249: DMNetworkGetVariableOffset(networkdm,e,&varoffset);
250: DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe);
252: /* set initial values for this pipe */
253: PipeComputeSteadyState(pipe,wash->Q0,wash->H0);
254: VecGetSize(pipe->x,&nx);
256: VecGetArrayRead(pipe->x,&xarray);
257: /* copy pipe->x to xarray */
258: for (k=0; k<nx; k++) {
259: (xarr+varoffset)[k] = xarray[k];
260: }
262: /* set boundary values into vfrom and vto */
263: DMNetworkGetConnectedVertices(networkdm,e,&cone);
264: vfrom = cone[0]; /* local ordering */
265: vto = cone[1];
266: DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
267: DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
269: /* if vform is a head vertex: */
270: DMNetworkGetComponent(networkdm,vfrom,0,&vkey,(void**)&junction);
271: if (junction->type == RESERVOIR) {
272: (xarr+offsetfrom)[1] = wash->H0; /* 1st H */
273: }
275: /* if vto is an end vertex: */
276: DMNetworkGetComponent(networkdm,vto,0,&vkey,(void**)&junction);
277: if (junction->type == VALVE) {
278: (xarr+offsetto)[0] = wash->QL; /* last Q */
279: }
280: VecRestoreArrayRead(pipe->x,&xarray);
281: }
283: VecRestoreArray(localX,&xarr);
284: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
285: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
286: DMRestoreLocalVector(networkdm,&localX);
288: #if 0
289: PetscInt N;
290: VecGetSize(X,&N);
291: PetscPrintf(PETSC_COMM_WORLD,"initial solution %d:\n",N);
292: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
293: #endif
294: return(0);
295: }
297: PetscErrorCode TSDMNetworkMonitor(TS ts, PetscInt step, PetscReal t, Vec x, void *context)
298: {
299: PetscErrorCode ierr;
300: DMNetworkMonitor monitor;
303: monitor = (DMNetworkMonitor)context;
304: DMNetworkMonitorView(monitor,x);
305: return(0);
306: }
308: PetscErrorCode PipesView(Vec X,DM networkdm,Wash wash)
309: {
310: PetscErrorCode ierr;
311: Pipe pipe;
312: PetscInt key,Start,End;
313: PetscMPIInt rank;
314: PetscInt nx,nnodes,nidx,*idx1,*idx2,*idx1_h,*idx2_h,idx_start,i,k,k1,xstart,j1;
315: Vec Xq,Xh,localX;
316: IS is1_q,is2_q,is1_h,is2_h;
317: VecScatter ctx_q,ctx_h;
320: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
322: /* get num of local and global total nnodes */
323: nidx = wash->nnodes_loc;
324: MPIU_Allreduce(&nidx,&nx,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
326: VecCreate(PETSC_COMM_WORLD,&Xq);
327: if (rank == 0) { /* all entries of Xq are in proc[0] */
328: VecSetSizes(Xq,nx,PETSC_DECIDE);
329: } else {
330: VecSetSizes(Xq,0,PETSC_DECIDE);
331: }
332: VecSetFromOptions(Xq);
333: VecSet(Xq,0.0);
334: VecDuplicate(Xq,&Xh);
336: DMGetLocalVector(networkdm,&localX);
338: /* set idx1 and idx2 */
339: PetscCalloc4(nidx,&idx1,nidx,&idx2,nidx,&idx1_h,nidx,&idx2_h);
341: DMNetworkGetEdgeRange(networkdm,&Start, &End);
343: VecGetOwnershipRange(X,&xstart,NULL);
344: k1 = 0;
345: j1 = 0;
346: for (i = Start; i < End; i++) {
347: DMNetworkGetComponent(networkdm,i,0,&key,(void**)&pipe);
348: nnodes = pipe->nnodes;
349: idx_start = pipe->id*nnodes;
350: for (k=0; k<nnodes; k++) {
351: idx1[k1] = xstart + j1*2*nnodes + 2*k;
352: idx2[k1] = idx_start + k;
354: idx1_h[k1] = xstart + j1*2*nnodes + 2*k + 1;
355: idx2_h[k1] = idx_start + k;
356: k1++;
357: }
358: j1++;
359: }
361: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx1,PETSC_COPY_VALUES,&is1_q);
362: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx2,PETSC_COPY_VALUES,&is2_q);
363: VecScatterCreate(X,is1_q,Xq,is2_q,&ctx_q);
364: VecScatterBegin(ctx_q,X,Xq,INSERT_VALUES,SCATTER_FORWARD);
365: VecScatterEnd(ctx_q,X,Xq,INSERT_VALUES,SCATTER_FORWARD);
367: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx1_h,PETSC_COPY_VALUES,&is1_h);
368: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx2_h,PETSC_COPY_VALUES,&is2_h);
369: VecScatterCreate(X,is1_h,Xh,is2_h,&ctx_h);
370: VecScatterBegin(ctx_h,X,Xh,INSERT_VALUES,SCATTER_FORWARD);
371: VecScatterEnd(ctx_h,X,Xh,INSERT_VALUES,SCATTER_FORWARD);
373: PetscPrintf(PETSC_COMM_WORLD,"Xq: \n");
374: VecView(Xq,PETSC_VIEWER_STDOUT_WORLD);
375: PetscPrintf(PETSC_COMM_WORLD,"Xh: \n");
376: VecView(Xh,PETSC_VIEWER_STDOUT_WORLD);
378: VecScatterDestroy(&ctx_q);
379: PetscFree4(idx1,idx2,idx1_h,idx2_h);
380: ISDestroy(&is1_q);
381: ISDestroy(&is2_q);
383: VecScatterDestroy(&ctx_h);
384: ISDestroy(&is1_h);
385: ISDestroy(&is2_h);
387: VecDestroy(&Xq);
388: VecDestroy(&Xh);
389: DMRestoreLocalVector(networkdm,&localX);
390: return(0);
391: }
393: PetscErrorCode WashNetworkCleanUp(Wash wash)
394: {
396: PetscMPIInt rank;
399: MPI_Comm_rank(wash->comm,&rank);
400: PetscFree(wash->edgelist);
401: PetscFree(wash->vtype);
402: if (!rank) {
403: PetscFree2(wash->junction,wash->pipe);
404: }
405: return(0);
406: }
408: PetscErrorCode WashNetworkCreate(MPI_Comm comm,PetscInt pipesCase,Wash *wash_ptr)
409: {
411: PetscInt npipes;
412: PetscMPIInt rank;
413: Wash wash=NULL;
414: PetscInt i,numVertices,numEdges,*vtype;
415: PetscInt *edgelist;
416: Junction junctions=NULL;
417: Pipe pipes=NULL;
418: PetscBool washdist=PETSC_TRUE;
421: MPI_Comm_rank(comm,&rank);
423: PetscCalloc1(1,&wash);
424: wash->comm = comm;
425: *wash_ptr = wash;
426: wash->Q0 = 0.477432; /* RESERVOIR */
427: wash->H0 = 150.0;
428: wash->HL = 143.488; /* VALVE */
429: wash->QL = 0.0;
430: wash->nnodes_loc = 0;
432: numVertices = 0;
433: numEdges = 0;
434: edgelist = NULL;
436: /* proc[0] creates a sequential wash and edgelist */
437: PetscPrintf(PETSC_COMM_WORLD,"Setup pipesCase %D\n",pipesCase);
439: /* Set global number of pipes, edges, and junctions */
440: /*-------------------------------------------------*/
441: switch (pipesCase) {
442: case 0:
443: /* pipeCase 0: */
444: /* =================================================
445: (RESERVOIR) v0 --E0--> v1--E1--> v2 --E2-->v3 (VALVE)
446: ==================================================== */
447: npipes = 3;
448: PetscOptionsGetInt(NULL,NULL, "-npipes", &npipes, NULL);
449: wash->nedge = npipes;
450: wash->nvertex = npipes + 1;
452: /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
453: numVertices = 0;
454: numEdges = 0;
455: edgelist = NULL;
456: if (!rank) {
457: numVertices = wash->nvertex;
458: numEdges = wash->nedge;
460: PetscCalloc1(2*numEdges,&edgelist);
461: for (i=0; i<numEdges; i++) {
462: edgelist[2*i] = i; edgelist[2*i+1] = i+1;
463: }
465: /* Add network components */
466: /*------------------------*/
467: PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
469: /* vertex */
470: for (i=0; i<numVertices; i++) {
471: junctions[i].id = i;
472: junctions[i].type = JUNCTION;
473: }
475: junctions[0].type = RESERVOIR;
476: junctions[numVertices-1].type = VALVE;
477: }
478: break;
479: case 1:
480: /* pipeCase 1: */
481: /* ==========================
482: v2 (VALVE)
483: ^
484: |
485: E2
486: |
487: v0 --E0--> v3--E1--> v1
488: (RESERVOIR) (RESERVOIR)
489: ============================= */
490: npipes = 3;
491: wash->nedge = npipes;
492: wash->nvertex = npipes + 1;
494: /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
495: if (!rank) {
496: numVertices = wash->nvertex;
497: numEdges = wash->nedge;
499: PetscCalloc1(2*numEdges,&edgelist);
500: edgelist[0] = 0; edgelist[1] = 3; /* edge[0] */
501: edgelist[2] = 3; edgelist[3] = 1; /* edge[1] */
502: edgelist[4] = 3; edgelist[5] = 2; /* edge[2] */
504: /* Add network components */
505: /*------------------------*/
506: PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
507: /* vertex */
508: for (i=0; i<numVertices; i++) {
509: junctions[i].id = i;
510: junctions[i].type = JUNCTION;
511: }
513: junctions[0].type = RESERVOIR;
514: junctions[1].type = VALVE;
515: junctions[2].type = VALVE;
516: }
517: break;
518: case 2:
519: /* pipeCase 2: */
520: /* ==========================
521: (RESERVOIR) v2--> E2
522: |
523: v0 --E0--> v3--E1--> v1
524: (RESERVOIR) (VALVE)
525: ============================= */
527: /* Set application parameters -- to be used in function evalutions */
528: npipes = 3;
529: wash->nedge = npipes;
530: wash->nvertex = npipes + 1;
532: /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
533: if (!rank) {
534: numVertices = wash->nvertex;
535: numEdges = wash->nedge;
537: PetscCalloc1(2*numEdges,&edgelist);
538: edgelist[0] = 0; edgelist[1] = 3; /* edge[0] */
539: edgelist[2] = 3; edgelist[3] = 1; /* edge[1] */
540: edgelist[4] = 2; edgelist[5] = 3; /* edge[2] */
542: /* Add network components */
543: /*------------------------*/
544: PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
545: /* vertex */
546: for (i=0; i<numVertices; i++) {
547: junctions[i].id = i;
548: junctions[i].type = JUNCTION;
549: }
551: junctions[0].type = RESERVOIR;
552: junctions[1].type = VALVE;
553: junctions[2].type = RESERVOIR;
554: }
555: break;
556: default:
557: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"not done yet");
558: }
560: /* set edge global id */
561: for (i=0; i<numEdges; i++) pipes[i].id = i;
563: if (!rank) { /* set vtype for proc[0] */
564: PetscInt v;
565: PetscMalloc1(2*numEdges,&vtype);
566: for (i=0; i<2*numEdges; i++) {
567: v = edgelist[i];
568: vtype[i] = junctions[v].type;
569: }
570: wash->vtype = vtype;
571: }
573: *wash_ptr = wash;
574: wash->nedge = numEdges;
575: wash->nvertex = numVertices;
576: wash->edgelist = edgelist;
577: wash->junction = junctions;
578: wash->pipe = pipes;
580: /* Distribute edgelist to other processors */
581: PetscOptionsGetBool(NULL,NULL,"-wash_distribute",&washdist,NULL);
582: if (washdist) {
583: /*
584: PetscPrintf(PETSC_COMM_WORLD," Distribute sequential wash ...\n");
585: */
586: WashNetworkDistribute(comm,wash);
587: }
588: return(0);
589: }
591: /* ------------------------------------------------------- */
592: int main(int argc,char ** argv)
593: {
594: PetscErrorCode ierr;
595: Wash wash;
596: Junction junctions,junction;
597: Pipe pipe,pipes;
598: PetscInt KeyPipe,KeyJunction;
599: PetscInt *edgelist = NULL,*edgelists[1],*vtype = NULL;
600: PetscInt i,e,v,eStart,eEnd,vStart,vEnd,key;
601: PetscInt vkey,type;
602: const PetscInt *cone;
603: DM networkdm;
604: PetscMPIInt size,rank;
605: PetscReal ftime;
606: Vec X;
607: TS ts;
608: PetscInt steps=1;
609: TSConvergedReason reason;
610: PetscBool viewpipes,monipipes=PETSC_FALSE,userJac=PETSC_TRUE,viewdm=PETSC_FALSE,viewX=PETSC_FALSE;
611: PetscBool test=PETSC_FALSE;
612: PetscInt pipesCase=0;
613: DMNetworkMonitor monitor;
614: MPI_Comm comm;
616: PetscInt nedges,nvertices; /* local num of edges and vertices */
617: PetscInt nnodes = 6,nv,ne;
618: const PetscInt *vtx,*edge;
620: PetscInitialize(&argc,&argv,"pOption",help);if (ierr) return ierr;
622: /* Read runtime options */
623: PetscOptionsGetInt(NULL,NULL, "-case", &pipesCase, NULL);
624: PetscOptionsGetBool(NULL,NULL,"-user_Jac",&userJac,NULL);
625: PetscOptionsGetBool(NULL,NULL,"-pipe_monitor",&monipipes,NULL);
626: PetscOptionsGetBool(NULL,NULL,"-viewdm",&viewdm,NULL);
627: PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);
628: PetscOptionsGetInt(NULL,NULL, "-npipenodes", &nnodes, NULL);
629: PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);
631: /* Create networkdm */
632: /*------------------*/
633: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
634: PetscObjectGetComm((PetscObject)networkdm,&comm);
635: MPI_Comm_rank(comm,&rank);
636: MPI_Comm_size(comm,&size);
638: if (size == 1 && monipipes) {
639: DMNetworkMonitorCreate(networkdm,&monitor);
640: }
642: /* Register the components in the network */
643: DMNetworkRegisterComponent(networkdm,"junctionstruct",sizeof(struct _p_Junction),&KeyJunction);
644: DMNetworkRegisterComponent(networkdm,"pipestruct",sizeof(struct _p_Pipe),&KeyPipe);
646: /* Create a distributed wash network (user-specific) */
647: WashNetworkCreate(comm,pipesCase,&wash);
648: nedges = wash->nedge;
649: nvertices = wash->nvertex; /* local num of vertices, excluding ghosts */
650: edgelist = wash->edgelist;
651: vtype = wash->vtype;
652: junctions = wash->junction;
653: pipes = wash->pipe;
655: #if 0
656: PetscPrintf(PETSC_COMM_SELF,"[%d] after WashNetworkDistribute...ne %d, nv %d\n",rank,nedges,nvertices);
657: if (rank == 0) {
658: for (e=0; e<nedges; e++) PetscPrintf(PETSC_COMM_SELF," edge %d --> %d\n",edgelist[2*e],edgelist[2*e+1]);
659: }
660: #endif
662: /* Set up the network layout */
663: DMNetworkSetSizes(networkdm,1,&nvertices,&nedges,0,NULL);
665: /* Add local edge connectivity */
666: edgelists[0] = edgelist;
667: DMNetworkSetEdgeList(networkdm,edgelists,NULL);
668: DMNetworkLayoutSetUp(networkdm);
670: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
671: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
672: /* PetscPrintf(PETSC_COMM_SELF,"[%d] eStart/End: %d - %d; vStart/End: %d - %d\n",rank,eStart,eEnd,vStart,vEnd); */
674: /* Test DMNetworkGetSubnetworkInfo() */
675: if (test) {
676: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edge);
677: if (ne != eEnd - eStart || nv != vEnd - vStart) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"ne %D or nv %D is incorrect",ne,nv);
678: }
680: if (rank) { /* junctions[] and pipes[] for proc[0] are allocated in WashNetworkCreate() */
681: /* vEnd - vStart = nvertices + num of ghost vertices! */
682: PetscCalloc2(vEnd - vStart,&junctions,nedges,&pipes);
683: }
685: /* Add Pipe component to all local edges */
686: for (e = eStart; e < eEnd; e++) {
687: if (test) {
688: if (e != edge[e]) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"e %D != edge %D from DMNetworkGetSubnetworkInfo()",e,edge[e]);
689: }
691: pipes[e-eStart].nnodes = nnodes;
692: DMNetworkAddComponent(networkdm,e,KeyPipe,&pipes[e-eStart]);
694: /* Add number of variables to each edge */
695: DMNetworkAddNumVariables(networkdm,e,2*pipes[e-eStart].nnodes);
697: if (size == 1 && monipipes) { /* Add monitor -- show Q_{pipes[e-eStart].id}? */
698: pipes[e-eStart].length = 600.0;
699: DMNetworkMonitorAdd(monitor, "Pipe Q", e, pipes[e-eStart].nnodes, 0, 2, 0.0,pipes[e-eStart].length, -0.8, 0.8, PETSC_TRUE);
700: DMNetworkMonitorAdd(monitor, "Pipe H", e, pipes[e-eStart].nnodes, 1, 2, 0.0,pipes[e-eStart].length, -400.0, 800.0, PETSC_TRUE);
701: }
702: }
704: /* Add Junction component to all local vertices, including ghost vertices! */
705: for (v = vStart; v < vEnd; v++) {
706: if (test) {
707: if (v != vtx[v-vStart]) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"v %D != vtx %D from DMNetworkGetSubnetworkInfo()",v,vtx[v-vStart]);
708: }
710: DMNetworkAddComponent(networkdm,v,KeyJunction,&junctions[v-vStart]);
712: /* Add number of variables to vertex */
713: DMNetworkAddNumVariables(networkdm,v,2);
714: }
716: if (size > 1) { /* must be called before DMSetUp()???. Other partitioners do not work yet??? -- cause crash in proc[0]! */
717: DM plexdm;
718: PetscPartitioner part;
719: DMNetworkGetPlex(networkdm,&plexdm);
720: DMPlexGetPartitioner(plexdm, &part);
721: PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);
722: }
724: /* Set up DM for use */
725: DMSetUp(networkdm);
726: if (viewdm) {
727: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");
728: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
729: }
731: /* Set user physical parameters to the components */
732: for (e = eStart; e < eEnd; e++) {
733: DMNetworkGetConnectedVertices(networkdm,e,&cone);
734: /* vfrom */
735: DMNetworkGetComponent(networkdm,cone[0],0,&vkey,(void**)&junction);
736: junction->type = (VertexType)vtype[2*e];
738: /* vto */
739: DMNetworkGetComponent(networkdm,cone[1],0,&vkey,(void**)&junction);
740: junction->type = (VertexType)vtype[2*e+1];
741: }
743: WashNetworkCleanUp(wash);
745: /* Network partitioning and distribution of data */
746: DMNetworkDistribute(&networkdm,0);
747: if (viewdm) {
748: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
749: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
750: }
752: /* create vectors */
753: DMCreateGlobalVector(networkdm,&X);
754: DMCreateLocalVector(networkdm,&wash->localX);
755: DMCreateLocalVector(networkdm,&wash->localXdot);
757: /* PipeSetUp -- each process only sets its own pipes */
758: /*---------------------------------------------------*/
759: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
761: userJac = PETSC_TRUE;
762: DMNetworkHasJacobian(networkdm,userJac,userJac);
763: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
764: for (e=eStart; e<eEnd; e++) { /* each edge has only one component, pipe */
765: DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe);
767: wash->nnodes_loc += pipe->nnodes; /* local total num of nodes, will be used by PipesView() */
768: PipeSetParameters(pipe,
769: 600.0, /* length */
770: 0.5, /* diameter */
771: 1200.0, /* a */
772: 0.018); /* friction */
773: PipeSetUp(pipe);
775: if (userJac) {
776: /* Create Jacobian matrix structures for a Pipe */
777: Mat *J;
778: PipeCreateJacobian(pipe,NULL,&J);
779: DMNetworkEdgeSetMatrix(networkdm,e,J);
780: }
781: }
783: if (userJac) {
784: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
785: for (v=vStart; v<vEnd; v++) {
786: Mat *J;
787: JunctionCreateJacobian(networkdm,v,NULL,&J);
788: DMNetworkVertexSetMatrix(networkdm,v,J);
790: DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction);
791: junction->jacobian = J;
792: }
793: }
795: /* Test DMNetworkGetSubnetworkInfo() */
796: if (test) {
797: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
798: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edge);
799: if (ne != eEnd - eStart || nv != vEnd - vStart) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"ne %D or nv %D is incorrect",ne,nv);
801: for (e = eStart; e < eEnd; e++) {
802: if (e != edge[e]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"e %D != edge %D from DMNetworkGetSubnetworkInfo()",e,edge[e]);
803: }
804: for (v = vStart; v < vEnd; v++) {
805: if (v != vtx[v-vStart]) SETERRQ2(PetscObjectComm((PetscObject)networkdm),PETSC_ERR_ARG_WRONG,"v %D != vtx %D from DMNetworkGetSubnetworkInfo()",v,vtx[v-vStart]);
806: }
807: }
809: /* Setup solver */
810: /*--------------------------------------------------------*/
811: TSCreate(PETSC_COMM_WORLD,&ts);
813: TSSetDM(ts,(DM)networkdm);
814: TSSetIFunction(ts,NULL,WASHIFunction,wash);
816: TSSetMaxSteps(ts,steps);
817: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
818: TSSetTimeStep(ts,0.1);
819: TSSetType(ts,TSBEULER);
820: if (size == 1 && monipipes) {
821: TSMonitorSet(ts, TSDMNetworkMonitor, monitor, NULL);
822: }
823: TSSetFromOptions(ts);
825: WASHSetInitialSolution(networkdm,X,wash);
827: TSSolve(ts,X);
829: TSGetSolveTime(ts,&ftime);
830: TSGetStepNumber(ts,&steps);
831: TSGetConvergedReason(ts,&reason);
832: PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
833: if (viewX) {
834: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
835: }
837: viewpipes = PETSC_FALSE;
838: PetscOptionsGetBool(NULL,NULL, "-Jac_view", &viewpipes,NULL);
839: if (viewpipes) {
840: SNES snes;
841: Mat Jac;
842: TSGetSNES(ts,&snes);
843: SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
844: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
845: }
847: /* View solution q and h */
848: /* --------------------- */
849: viewpipes = PETSC_FALSE;
850: PetscOptionsGetBool(NULL,NULL, "-pipe_view", &viewpipes,NULL);
851: if (viewpipes) {
852: PipesView(X,networkdm,wash);
853: }
855: /* Free spaces */
856: /* ----------- */
857: TSDestroy(&ts);
858: VecDestroy(&X);
859: VecDestroy(&wash->localX);
860: VecDestroy(&wash->localXdot);
862: /* Destroy objects from each pipe that are created in PipeSetUp() */
863: DMNetworkGetEdgeRange(networkdm,&eStart, &eEnd);
864: for (i = eStart; i < eEnd; i++) {
865: DMNetworkGetComponent(networkdm,i,0,&key,(void**)&pipe);
866: PipeDestroy(&pipe);
867: }
868: if (userJac) {
869: for (v=vStart; v<vEnd; v++) {
870: DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction);
871: JunctionDestroyJacobian(networkdm,v,junction);
872: }
873: }
875: if (size == 1 && monipipes) {
876: DMNetworkMonitorDestroy(&monitor);
877: }
878: DMDestroy(&networkdm);
879: PetscFree(wash);
881: if (rank) {
882: PetscFree2(junctions,pipes);
883: }
884: PetscFinalize();
885: return ierr;
886: }
888: /*TEST
890: build:
891: depends: pipeInterface.c pipeImpls.c
893: test:
894: args: -ts_monitor -case 1 -ts_max_steps 1 -options_left no -viewX -test
895: localrunfiles: pOption
896: output_file: output/pipes1_1.out
898: test:
899: suffix: 2
900: nsize: 2
901: requires: mumps
902: args: -ts_monitor -case 1 -ts_max_steps 1 -petscpartitioner_type simple -options_left no -viewX -test
903: localrunfiles: pOption
904: output_file: output/pipes1_2.out
906: test:
907: suffix: 3
908: nsize: 2
909: requires: mumps
910: args: -ts_monitor -case 0 -ts_max_steps 1 -petscpartitioner_type simple -options_left no -viewX -test
911: localrunfiles: pOption
912: output_file: output/pipes1_3.out
914: test:
915: suffix: 4
916: args: -ts_monitor -case 2 -ts_max_steps 1 -options_left no -viewX -test
917: localrunfiles: pOption
918: output_file: output/pipes1_4.out
920: test:
921: suffix: 5
922: nsize: 3
923: requires: mumps
924: args: -ts_monitor -case 2 -ts_max_steps 10 -petscpartitioner_type simple -options_left no -viewX -test
925: localrunfiles: pOption
926: output_file: output/pipes1_5.out
928: test:
929: suffix: 6
930: nsize: 2
931: requires: mumps
932: args: -ts_monitor -case 1 -ts_max_steps 1 -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX -test
933: localrunfiles: pOption
934: output_file: output/pipes1_6.out
936: test:
937: suffix: 7
938: nsize: 2
939: requires: mumps
940: args: -ts_monitor -case 2 -ts_max_steps 1 -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX -test
941: localrunfiles: pOption
942: output_file: output/pipes1_7.out
944: TEST*/