Actual source code: pipes.c
1: static char help[] = "This example demonstrates DMNetwork. It is used for testing parallel generation of dmnetwork, then redistribute. \n\\n";
2: /*
3: Example: mpiexec -n <np> ./pipes -ts_max_steps 10
4: */
6: #include "wash.h"
8: /*
9: WashNetworkDistribute - proc[0] distributes sequential wash object
10: Input Parameters:
11: . comm - MPI communicator
12: . wash - wash context with all network data in proc[0]
14: Output Parameter:
15: . wash - wash context with nedge, nvertex and edgelist distributed
17: Note: The routine is used for testing parallel generation of dmnetwork, then redistribute.
18: */
19: PetscErrorCode WashNetworkDistribute(MPI_Comm comm,Wash wash)
20: {
21: PetscMPIInt rank,size,tag=0;
22: PetscInt i,e,v,numEdges,numVertices,nedges,*eowners=NULL,estart,eend,*vtype=NULL,nvertices;
23: PetscInt *edgelist = wash->edgelist,*nvtx=NULL,*vtxDone=NULL;
25: MPI_Comm_size(comm,&size);
26: if (size == 1) return 0;
28: MPI_Comm_rank(comm,&rank);
29: numEdges = wash->nedge;
30: numVertices = wash->nvertex;
32: /* (1) all processes get global and local number of edges */
33: MPI_Bcast(&numEdges,1,MPIU_INT,0,comm);
34: nedges = numEdges/size; /* local nedges */
35: if (rank == 0) {
36: nedges += numEdges - size*(numEdges/size);
37: }
38: wash->Nedge = numEdges;
39: wash->nedge = nedges;
40: /* PetscPrintf(PETSC_COMM_SELF,"[%d] nedges %d, numEdges %d\n",rank,nedges,numEdges); */
42: PetscCalloc3(size+1,&eowners,size,&nvtx,numVertices,&vtxDone);
43: MPI_Allgather(&nedges,1,MPIU_INT,eowners+1,1,MPIU_INT,PETSC_COMM_WORLD);
44: eowners[0] = 0;
45: for (i=2; i<=size; i++) {
46: eowners[i] += eowners[i-1];
47: }
49: estart = eowners[rank];
50: eend = eowners[rank+1];
51: /* PetscPrintf(PETSC_COMM_SELF,"[%d] own lists row %d - %d\n",rank,estart,eend); */
53: /* (2) distribute row block edgelist to all processors */
54: if (rank == 0) {
55: vtype = wash->vtype;
56: for (i=1; i<size; i++) {
57: /* proc[0] sends edgelist to proc[i] */
58: MPI_Send(edgelist+2*eowners[i],2*(eowners[i+1]-eowners[i]),MPIU_INT,i,tag,comm);
60: /* proc[0] sends vtype to proc[i] */
61: MPI_Send(vtype+2*eowners[i],2*(eowners[i+1]-eowners[i]),MPIU_INT,i,tag,comm);
62: }
63: } else {
64: MPI_Status status;
65: PetscMalloc1(2*(eend-estart),&vtype);
66: PetscMalloc1(2*(eend-estart),&edgelist);
68: MPI_Recv(edgelist,2*(eend-estart),MPIU_INT,0,tag,comm,&status);
69: MPI_Recv(vtype,2*(eend-estart),MPIU_INT,0,tag,comm,&status);
70: }
72: wash->edgelist = edgelist;
74: /* (3) all processes get global and local number of vertices, without ghost vertices */
75: if (rank == 0) {
76: for (i=0; i<size; i++) {
77: for (e=eowners[i]; e<eowners[i+1]; e++) {
78: v = edgelist[2*e];
79: if (!vtxDone[v]) {
80: nvtx[i]++; vtxDone[v] = 1;
81: }
82: v = edgelist[2*e+1];
83: if (!vtxDone[v]) {
84: nvtx[i]++; vtxDone[v] = 1;
85: }
86: }
87: }
88: }
89: MPI_Bcast(&numVertices,1,MPIU_INT,0,PETSC_COMM_WORLD);
90: MPI_Scatter(nvtx,1,MPIU_INT,&nvertices,1,MPIU_INT,0,PETSC_COMM_WORLD);
91: PetscFree3(eowners,nvtx,vtxDone);
93: wash->Nvertex = numVertices;
94: wash->nvertex = nvertices;
95: wash->vtype = vtype;
96: return 0;
97: }
99: PetscErrorCode WASHIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void* ctx)
100: {
101: Wash wash=(Wash)ctx;
102: DM networkdm;
103: Vec localX,localXdot,localF, localXold;
104: const PetscInt *cone;
105: PetscInt vfrom,vto,offsetfrom,offsetto,varoffset;
106: PetscInt v,vStart,vEnd,e,eStart,eEnd;
107: PetscInt nend,type;
108: PetscBool ghost;
109: PetscScalar *farr,*juncf, *pipef;
110: PetscReal dt;
111: Pipe pipe;
112: PipeField *pipex,*pipexdot,*juncx;
113: Junction junction;
114: DMDALocalInfo info;
115: const PetscScalar *xarr,*xdotarr, *xoldarr;
117: localX = wash->localX;
118: localXdot = wash->localXdot;
120: TSGetSolution(ts,&localXold);
121: TSGetDM(ts,&networkdm);
122: TSGetTimeStep(ts,&dt);
123: DMGetLocalVector(networkdm,&localF);
125: /* Set F and localF as zero */
126: VecSet(F,0.0);
127: VecSet(localF,0.0);
129: /* Update ghost values of locaX and locaXdot */
130: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
131: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
133: DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot);
134: DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot);
136: VecGetArrayRead(localX,&xarr);
137: VecGetArrayRead(localXdot,&xdotarr);
138: VecGetArrayRead(localXold,&xoldarr);
139: VecGetArray(localF,&farr);
141: /* junction->type == JUNCTION:
142: juncf[0] = -qJ + sum(qin); juncf[1] = qJ - sum(qout)
143: junction->type == RESERVOIR (upper stream):
144: juncf[0] = -hJ + H0; juncf[1] = qJ - sum(qout)
145: junction->type == VALVE (down stream):
146: juncf[0] = -qJ + sum(qin); juncf[1] = qJ
147: */
148: /* Vertex/junction initialization */
149: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
150: for (v=vStart; v<vEnd; v++) {
151: DMNetworkIsGhostVertex(networkdm,v,&ghost);
152: if (ghost) continue;
154: DMNetworkGetComponent(networkdm,v,0,&type,(void**)&junction,NULL);
155: DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&varoffset);
156: juncx = (PipeField*)(xarr + varoffset);
157: juncf = (PetscScalar*)(farr + varoffset);
159: juncf[0] = -juncx[0].q;
160: juncf[1] = juncx[0].q;
162: if (junction->type == RESERVOIR) { /* upstream reservoir */
163: juncf[0] = juncx[0].h - wash->H0;
164: }
165: }
167: /* Edge/pipe */
168: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
169: for (e=eStart; e<eEnd; e++) {
170: DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe,NULL);
171: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&varoffset);
172: pipex = (PipeField*)(xarr + varoffset);
173: pipexdot = (PipeField*)(xdotarr + varoffset);
174: pipef = (PetscScalar*)(farr + varoffset);
176: /* Get some data into the pipe structure: note, some of these operations
177: * might be redundant. Will it consume too much time? */
178: pipe->dt = dt;
179: pipe->xold = (PipeField*)(xoldarr + varoffset);
181: /* Evaluate F over this edge/pipe: pipef[1], ...,pipef[2*nend] */
182: DMDAGetLocalInfo(pipe->da,&info);
183: PipeIFunctionLocal_Lax(&info,t,pipex,pipexdot,pipef,pipe);
185: /* Get boundary values from connected vertices */
186: DMNetworkGetConnectedVertices(networkdm,e,&cone);
187: vfrom = cone[0]; /* local ordering */
188: vto = cone[1];
189: DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
190: DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);
192: /* Evaluate upstream boundary */
193: DMNetworkGetComponent(networkdm,vfrom,0,&type,(void**)&junction,NULL);
195: juncx = (PipeField*)(xarr + offsetfrom);
196: juncf = (PetscScalar*)(farr + offsetfrom);
198: pipef[0] = pipex[0].h - juncx[0].h;
199: juncf[1] -= pipex[0].q;
201: /* Evaluate downstream boundary */
202: DMNetworkGetComponent(networkdm,vto,0,&type,(void**)&junction,NULL);
204: juncx = (PipeField*)(xarr + offsetto);
205: juncf = (PetscScalar*)(farr + offsetto);
206: nend = pipe->nnodes - 1;
208: pipef[2*nend + 1] = pipex[nend].h - juncx[0].h;
209: juncf[0] += pipex[nend].q;
210: }
212: VecRestoreArrayRead(localX,&xarr);
213: VecRestoreArrayRead(localXdot,&xdotarr);
214: VecRestoreArray(localF,&farr);
216: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
217: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
218: DMRestoreLocalVector(networkdm,&localF);
219: /*
220: PetscPrintf(PETSC_COMM_WORLD("F:\n");
221: VecView(F,PETSC_VIEWER_STDOUT_WORLD);
222: */
223: return 0;
224: }
226: PetscErrorCode WASHSetInitialSolution(DM networkdm,Vec X,Wash wash)
227: {
228: PetscInt k,nx,vkey,vfrom,vto,offsetfrom,offsetto;
229: PetscInt type,varoffset;
230: PetscInt e,eStart,eEnd;
231: Vec localX;
232: PetscScalar *xarr;
233: Pipe pipe;
234: Junction junction;
235: const PetscInt *cone;
236: const PetscScalar *xarray;
238: VecSet(X,0.0);
239: DMGetLocalVector(networkdm,&localX);
240: VecGetArray(localX,&xarr);
242: /* Edge */
243: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
244: for (e=eStart; e<eEnd; e++) {
245: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&varoffset);
246: DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe,NULL);
248: /* set initial values for this pipe */
249: PipeComputeSteadyState(pipe,wash->Q0,wash->H0);
250: VecGetSize(pipe->x,&nx);
252: VecGetArrayRead(pipe->x,&xarray);
253: /* copy pipe->x to xarray */
254: for (k=0; k<nx; k++) {
255: (xarr+varoffset)[k] = xarray[k];
256: }
258: /* set boundary values into vfrom and vto */
259: DMNetworkGetConnectedVertices(networkdm,e,&cone);
260: vfrom = cone[0]; /* local ordering */
261: vto = cone[1];
262: DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
263: DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);
265: /* if vform is a head vertex: */
266: DMNetworkGetComponent(networkdm,vfrom,0,&vkey,(void**)&junction,NULL);
267: if (junction->type == RESERVOIR) {
268: (xarr+offsetfrom)[1] = wash->H0; /* 1st H */
269: }
271: /* if vto is an end vertex: */
272: DMNetworkGetComponent(networkdm,vto,0,&vkey,(void**)&junction,NULL);
273: if (junction->type == VALVE) {
274: (xarr+offsetto)[0] = wash->QL; /* last Q */
275: }
276: VecRestoreArrayRead(pipe->x,&xarray);
277: }
279: VecRestoreArray(localX,&xarr);
280: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
281: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
282: DMRestoreLocalVector(networkdm,&localX);
284: #if 0
285: PetscInt N;
286: VecGetSize(X,&N);
287: PetscPrintf(PETSC_COMM_WORLD,"initial solution %d:\n",N);
288: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
289: #endif
290: return 0;
291: }
293: PetscErrorCode TSDMNetworkMonitor(TS ts, PetscInt step, PetscReal t, Vec x, void *context)
294: {
295: DMNetworkMonitor monitor;
297: monitor = (DMNetworkMonitor)context;
298: DMNetworkMonitorView(monitor,x);
299: return 0;
300: }
302: PetscErrorCode PipesView(DM networkdm,PetscInt KeyPipe,Vec X)
303: {
304: PetscInt i,numkeys=1,*blocksize,*numselectedvariable,**selectedvariables,n;
305: IS isfrom_q,isfrom_h,isfrom;
306: Vec Xto;
307: VecScatter ctx;
308: MPI_Comm comm;
310: PetscObjectGetComm((PetscObject)networkdm,&comm);
312: /* 1. Create isfrom_q for q-variable of pipes */
313: PetscMalloc3(numkeys,&blocksize,numkeys,&numselectedvariable,numkeys,&selectedvariables);
314: for (i=0; i<numkeys; i++) {
315: blocksize[i] = 2;
316: numselectedvariable[i] = 1;
317: PetscMalloc1(numselectedvariable[i],&selectedvariables[i]);
318: selectedvariables[i][0] = 0; /* q-variable */
319: }
320: DMNetworkCreateIS(networkdm,numkeys,&KeyPipe,blocksize,numselectedvariable,selectedvariables,&isfrom_q);
322: /* 2. Create Xto and isto */
323: ISGetLocalSize(isfrom_q, &n);
324: VecCreate(comm,&Xto);
325: VecSetSizes(Xto,n,PETSC_DECIDE);
326: VecSetFromOptions(Xto);
327: VecSet(Xto,0.0);
329: /* 3. Create scatter */
330: VecScatterCreate(X,isfrom_q,Xto,NULL,&ctx);
332: /* 4. Scatter to Xq */
333: VecScatterBegin(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
334: VecScatterEnd(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
335: VecScatterDestroy(&ctx);
336: ISDestroy(&isfrom_q);
338: PetscPrintf(PETSC_COMM_WORLD,"Xq:\n");
339: VecView(Xto,PETSC_VIEWER_STDOUT_WORLD);
341: /* 5. Create isfrom_h for h-variable of pipes; Create scatter; Scatter to Xh */
342: for (i=0; i<numkeys; i++) {
343: selectedvariables[i][0] = 1; /* h-variable */
344: }
345: DMNetworkCreateIS(networkdm,numkeys,&KeyPipe,blocksize,numselectedvariable,selectedvariables,&isfrom_h);
347: VecScatterCreate(X,isfrom_h,Xto,NULL,&ctx);
348: VecScatterBegin(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
349: VecScatterEnd(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
350: VecScatterDestroy(&ctx);
351: ISDestroy(&isfrom_h);
353: PetscPrintf(PETSC_COMM_WORLD,"Xh:\n");
354: VecView(Xto,PETSC_VIEWER_STDOUT_WORLD);
355: VecDestroy(&Xto);
357: /* 6. Create isfrom for all pipe variables; Create scatter; Scatter to Xpipes */
358: for (i=0; i<numkeys; i++) {
359: blocksize[i] = -1; /* select all the variables of the i-th component */
360: }
361: DMNetworkCreateIS(networkdm,numkeys,&KeyPipe,blocksize,NULL,NULL,&isfrom);
362: ISDestroy(&isfrom);
363: DMNetworkCreateIS(networkdm,numkeys,&KeyPipe,NULL,NULL,NULL,&isfrom);
365: ISGetLocalSize(isfrom, &n);
366: VecCreate(comm,&Xto);
367: VecSetSizes(Xto,n,PETSC_DECIDE);
368: VecSetFromOptions(Xto);
369: VecSet(Xto,0.0);
371: VecScatterCreate(X,isfrom,Xto,NULL,&ctx);
372: VecScatterBegin(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
373: VecScatterEnd(ctx,X,Xto,INSERT_VALUES,SCATTER_FORWARD);
374: VecScatterDestroy(&ctx);
375: ISDestroy(&isfrom);
377: PetscPrintf(PETSC_COMM_WORLD,"Xpipes:\n");
378: VecView(Xto,PETSC_VIEWER_STDOUT_WORLD);
380: /* 7. Free spaces */
381: for (i=0; i<numkeys; i++) {
382: PetscFree(selectedvariables[i]);
383: }
384: PetscFree3(blocksize,numselectedvariable,selectedvariables);
385: VecDestroy(&Xto);
386: return 0;
387: }
389: PetscErrorCode ISJunctionsView(DM networkdm,PetscInt KeyJunc)
390: {
391: PetscInt numkeys=1;
392: IS isfrom;
393: MPI_Comm comm;
394: PetscMPIInt rank;
396: PetscObjectGetComm((PetscObject)networkdm,&comm);
397: MPI_Comm_rank(comm,&rank);
399: /* Create a global isfrom for all junction variables */
400: DMNetworkCreateIS(networkdm,numkeys,&KeyJunc,NULL,NULL,NULL,&isfrom);
401: PetscPrintf(PETSC_COMM_WORLD,"ISJunctions:\n");
402: ISView(isfrom,PETSC_VIEWER_STDOUT_WORLD);
403: ISDestroy(&isfrom);
405: /* Create a local isfrom for all junction variables */
406: DMNetworkCreateLocalIS(networkdm,numkeys,&KeyJunc,NULL,NULL,NULL,&isfrom);
407: if (!rank) {
408: PetscPrintf(PETSC_COMM_SELF,"[%d] ISLocalJunctions:\n",rank);
409: ISView(isfrom,PETSC_VIEWER_STDOUT_SELF);
410: }
411: ISDestroy(&isfrom);
412: return 0;
413: }
415: PetscErrorCode WashNetworkCleanUp(Wash wash)
416: {
417: PetscMPIInt rank;
419: MPI_Comm_rank(wash->comm,&rank);
420: PetscFree(wash->edgelist);
421: PetscFree(wash->vtype);
422: if (rank == 0) {
423: PetscFree2(wash->junction,wash->pipe);
424: }
425: return 0;
426: }
428: PetscErrorCode WashNetworkCreate(MPI_Comm comm,PetscInt pipesCase,Wash *wash_ptr)
429: {
430: PetscInt npipes;
431: PetscMPIInt rank;
432: Wash wash=NULL;
433: PetscInt i,numVertices,numEdges,*vtype;
434: PetscInt *edgelist;
435: Junction junctions=NULL;
436: Pipe pipes=NULL;
437: PetscBool washdist=PETSC_TRUE;
439: MPI_Comm_rank(comm,&rank);
441: PetscCalloc1(1,&wash);
442: wash->comm = comm;
443: *wash_ptr = wash;
444: wash->Q0 = 0.477432; /* RESERVOIR */
445: wash->H0 = 150.0;
446: wash->HL = 143.488; /* VALVE */
447: wash->QL = 0.0;
448: wash->nnodes_loc = 0;
450: numVertices = 0;
451: numEdges = 0;
452: edgelist = NULL;
454: /* proc[0] creates a sequential wash and edgelist */
455: PetscPrintf(PETSC_COMM_WORLD,"Setup pipesCase %D\n",pipesCase);
457: /* Set global number of pipes, edges, and junctions */
458: /*-------------------------------------------------*/
459: switch (pipesCase) {
460: case 0:
461: /* pipeCase 0: */
462: /* =================================================
463: (RESERVOIR) v0 --E0--> v1--E1--> v2 --E2-->v3 (VALVE)
464: ==================================================== */
465: npipes = 3;
466: PetscOptionsGetInt(NULL,NULL, "-npipes", &npipes, NULL);
467: wash->nedge = npipes;
468: wash->nvertex = npipes + 1;
470: /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
471: numVertices = 0;
472: numEdges = 0;
473: edgelist = NULL;
474: if (rank == 0) {
475: numVertices = wash->nvertex;
476: numEdges = wash->nedge;
478: PetscCalloc1(2*numEdges,&edgelist);
479: for (i=0; i<numEdges; i++) {
480: edgelist[2*i] = i; edgelist[2*i+1] = i+1;
481: }
483: /* Add network components */
484: /*------------------------*/
485: PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
487: /* vertex */
488: for (i=0; i<numVertices; i++) {
489: junctions[i].id = i;
490: junctions[i].type = JUNCTION;
491: }
493: junctions[0].type = RESERVOIR;
494: junctions[numVertices-1].type = VALVE;
495: }
496: break;
497: case 1:
498: /* pipeCase 1: */
499: /* ==========================
500: v2 (VALVE)
501: ^
502: |
503: E2
504: |
505: v0 --E0--> v3--E1--> v1
506: (RESERVOIR) (RESERVOIR)
507: ============================= */
508: npipes = 3;
509: wash->nedge = npipes;
510: wash->nvertex = npipes + 1;
512: /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
513: if (rank == 0) {
514: numVertices = wash->nvertex;
515: numEdges = wash->nedge;
517: PetscCalloc1(2*numEdges,&edgelist);
518: edgelist[0] = 0; edgelist[1] = 3; /* edge[0] */
519: edgelist[2] = 3; edgelist[3] = 1; /* edge[1] */
520: edgelist[4] = 3; edgelist[5] = 2; /* edge[2] */
522: /* Add network components */
523: /*------------------------*/
524: PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
525: /* vertex */
526: for (i=0; i<numVertices; i++) {
527: junctions[i].id = i;
528: junctions[i].type = JUNCTION;
529: }
531: junctions[0].type = RESERVOIR;
532: junctions[1].type = VALVE;
533: junctions[2].type = VALVE;
534: }
535: break;
536: case 2:
537: /* pipeCase 2: */
538: /* ==========================
539: (RESERVOIR) v2--> E2
540: |
541: v0 --E0--> v3--E1--> v1
542: (RESERVOIR) (VALVE)
543: ============================= */
545: /* Set application parameters -- to be used in function evalutions */
546: npipes = 3;
547: wash->nedge = npipes;
548: wash->nvertex = npipes + 1;
550: /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
551: if (rank == 0) {
552: numVertices = wash->nvertex;
553: numEdges = wash->nedge;
555: PetscCalloc1(2*numEdges,&edgelist);
556: edgelist[0] = 0; edgelist[1] = 3; /* edge[0] */
557: edgelist[2] = 3; edgelist[3] = 1; /* edge[1] */
558: edgelist[4] = 2; edgelist[5] = 3; /* edge[2] */
560: /* Add network components */
561: /*------------------------*/
562: PetscCalloc2(numVertices,&junctions,numEdges,&pipes);
563: /* vertex */
564: for (i=0; i<numVertices; i++) {
565: junctions[i].id = i;
566: junctions[i].type = JUNCTION;
567: }
569: junctions[0].type = RESERVOIR;
570: junctions[1].type = VALVE;
571: junctions[2].type = RESERVOIR;
572: }
573: break;
574: default:
575: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"not done yet");
576: }
578: /* set edge global id */
579: for (i=0; i<numEdges; i++) pipes[i].id = i;
581: if (rank == 0) { /* set vtype for proc[0] */
582: PetscInt v;
583: PetscMalloc1(2*numEdges,&vtype);
584: for (i=0; i<2*numEdges; i++) {
585: v = edgelist[i];
586: vtype[i] = junctions[v].type;
587: }
588: wash->vtype = vtype;
589: }
591: *wash_ptr = wash;
592: wash->nedge = numEdges;
593: wash->nvertex = numVertices;
594: wash->edgelist = edgelist;
595: wash->junction = junctions;
596: wash->pipe = pipes;
598: /* Distribute edgelist to other processors */
599: PetscOptionsGetBool(NULL,NULL,"-wash_distribute",&washdist,NULL);
600: if (washdist) {
601: /*
602: PetscPrintf(PETSC_COMM_WORLD," Distribute sequential wash ...\n");
603: */
604: WashNetworkDistribute(comm,wash);
605: }
606: return 0;
607: }
609: /* ------------------------------------------------------- */
610: int main(int argc,char ** argv)
611: {
612: Wash wash;
613: Junction junctions,junction;
614: Pipe pipe,pipes;
615: PetscInt KeyPipe,KeyJunction,*edgelist = NULL,*vtype = NULL;
616: PetscInt i,e,v,eStart,eEnd,vStart,vEnd,key,vkey,type;
617: PetscInt steps=1,nedges,nnodes=6;
618: const PetscInt *cone;
619: DM networkdm;
620: PetscMPIInt size,rank;
621: PetscReal ftime;
622: Vec X;
623: TS ts;
624: TSConvergedReason reason;
625: PetscBool viewpipes,viewjuncs,monipipes=PETSC_FALSE,userJac=PETSC_TRUE,viewdm=PETSC_FALSE,viewX=PETSC_FALSE;
626: PetscInt pipesCase=0;
627: DMNetworkMonitor monitor;
628: MPI_Comm comm;
630: PetscInitialize(&argc,&argv,"pOption",help);
632: /* Read runtime options */
633: PetscOptionsGetInt(NULL,NULL, "-case", &pipesCase, NULL);
634: PetscOptionsGetBool(NULL,NULL,"-user_Jac",&userJac,NULL);
635: PetscOptionsGetBool(NULL,NULL,"-pipe_monitor",&monipipes,NULL);
636: PetscOptionsGetBool(NULL,NULL,"-viewdm",&viewdm,NULL);
637: PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);
638: PetscOptionsGetInt(NULL,NULL, "-npipenodes", &nnodes, NULL);
640: /* Create networkdm */
641: /*------------------*/
642: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
643: PetscObjectGetComm((PetscObject)networkdm,&comm);
644: MPI_Comm_rank(comm,&rank);
645: MPI_Comm_size(comm,&size);
647: if (size == 1 && monipipes) {
648: DMNetworkMonitorCreate(networkdm,&monitor);
649: }
651: /* Register the components in the network */
652: DMNetworkRegisterComponent(networkdm,"junctionstruct",sizeof(struct _p_Junction),&KeyJunction);
653: DMNetworkRegisterComponent(networkdm,"pipestruct",sizeof(struct _p_Pipe),&KeyPipe);
655: /* Create a distributed wash network (user-specific) */
656: WashNetworkCreate(comm,pipesCase,&wash);
657: nedges = wash->nedge;
658: edgelist = wash->edgelist;
659: vtype = wash->vtype;
660: junctions = wash->junction;
661: pipes = wash->pipe;
663: /* Set up the network layout */
664: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
665: DMNetworkAddSubnetwork(networkdm,NULL,nedges,edgelist,NULL);
667: DMNetworkLayoutSetUp(networkdm);
669: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
670: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
671: /* PetscPrintf(PETSC_COMM_SELF,"[%d] eStart/End: %d - %d; vStart/End: %d - %d\n",rank,eStart,eEnd,vStart,vEnd); */
673: if (rank) { /* junctions[] and pipes[] for proc[0] are allocated in WashNetworkCreate() */
674: /* vEnd - vStart = nvertices + number of ghost vertices! */
675: PetscCalloc2(vEnd - vStart,&junctions,nedges,&pipes);
676: }
678: /* Add Pipe component and number of variables to all local edges */
679: for (e = eStart; e < eEnd; e++) {
680: pipes[e-eStart].nnodes = nnodes;
681: DMNetworkAddComponent(networkdm,e,KeyPipe,&pipes[e-eStart],2*pipes[e-eStart].nnodes);
683: if (size == 1 && monipipes) { /* Add monitor -- show Q_{pipes[e-eStart].id}? */
684: pipes[e-eStart].length = 600.0;
685: DMNetworkMonitorAdd(monitor, "Pipe Q", e, pipes[e-eStart].nnodes, 0, 2, 0.0,pipes[e-eStart].length, -0.8, 0.8, PETSC_TRUE);
686: DMNetworkMonitorAdd(monitor, "Pipe H", e, pipes[e-eStart].nnodes, 1, 2, 0.0,pipes[e-eStart].length, -400.0, 800.0, PETSC_TRUE);
687: }
688: }
690: /* Add Junction component and number of variables to all local vertices */
691: for (v = vStart; v < vEnd; v++) {
692: DMNetworkAddComponent(networkdm,v,KeyJunction,&junctions[v-vStart],2);
693: }
695: if (size > 1) { /* must be called before DMSetUp()???. Other partitioners do not work yet??? -- cause crash in proc[0]! */
696: DM plexdm;
697: PetscPartitioner part;
698: DMNetworkGetPlex(networkdm,&plexdm);
699: DMPlexGetPartitioner(plexdm, &part);
700: PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);
701: PetscOptionsSetValue(NULL,"-dm_plex_csr_alg","mat"); /* for parmetis */
702: }
704: /* Set up DM for use */
705: DMSetUp(networkdm);
706: if (viewdm) {
707: PetscPrintf(PETSC_COMM_WORLD,"\nOriginal networkdm, DMView:\n");
708: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
709: }
711: /* Set user physical parameters to the components */
712: for (e = eStart; e < eEnd; e++) {
713: DMNetworkGetConnectedVertices(networkdm,e,&cone);
714: /* vfrom */
715: DMNetworkGetComponent(networkdm,cone[0],0,&vkey,(void**)&junction,NULL);
716: junction->type = (VertexType)vtype[2*e];
718: /* vto */
719: DMNetworkGetComponent(networkdm,cone[1],0,&vkey,(void**)&junction,NULL);
720: junction->type = (VertexType)vtype[2*e+1];
721: }
723: WashNetworkCleanUp(wash);
725: /* Network partitioning and distribution of data */
726: DMNetworkDistribute(&networkdm,0);
727: if (viewdm) {
728: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
729: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
730: }
732: /* create vectors */
733: DMCreateGlobalVector(networkdm,&X);
734: DMCreateLocalVector(networkdm,&wash->localX);
735: DMCreateLocalVector(networkdm,&wash->localXdot);
737: /* PipeSetUp -- each process only sets its own pipes */
738: /*---------------------------------------------------*/
739: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
741: userJac = PETSC_TRUE;
742: DMNetworkHasJacobian(networkdm,userJac,userJac);
743: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
744: for (e=eStart; e<eEnd; e++) { /* each edge has only one component, pipe */
745: DMNetworkGetComponent(networkdm,e,0,&type,(void**)&pipe,NULL);
747: wash->nnodes_loc += pipe->nnodes; /* local total number of nodes, will be used by PipesView() */
748: PetscCall(PipeSetParameters(pipe,
749: 600.0, /* length */
750: 0.5, /* diameter */
751: 1200.0, /* a */
752: 0.018)); /* friction */
753: PipeSetUp(pipe);
755: if (userJac) {
756: /* Create Jacobian matrix structures for a Pipe */
757: Mat *J;
758: PipeCreateJacobian(pipe,NULL,&J);
759: DMNetworkEdgeSetMatrix(networkdm,e,J);
760: }
761: }
763: if (userJac) {
764: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
765: for (v=vStart; v<vEnd; v++) {
766: Mat *J;
767: JunctionCreateJacobian(networkdm,v,NULL,&J);
768: DMNetworkVertexSetMatrix(networkdm,v,J);
770: DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction,NULL);
771: junction->jacobian = J;
772: }
773: }
775: /* Setup solver */
776: /*--------------------------------------------------------*/
777: TSCreate(PETSC_COMM_WORLD,&ts);
779: TSSetDM(ts,(DM)networkdm);
780: TSSetIFunction(ts,NULL,WASHIFunction,wash);
782: TSSetMaxSteps(ts,steps);
783: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
784: TSSetTimeStep(ts,0.1);
785: TSSetType(ts,TSBEULER);
786: if (size == 1 && monipipes) TSMonitorSet(ts, TSDMNetworkMonitor, monitor, NULL);
787: TSSetFromOptions(ts);
789: WASHSetInitialSolution(networkdm,X,wash);
791: TSSolve(ts,X);
793: TSGetSolveTime(ts,&ftime);
794: TSGetStepNumber(ts,&steps);
795: TSGetConvergedReason(ts,&reason);
796: PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
797: if (viewX) {
798: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
799: }
801: viewpipes = PETSC_FALSE;
802: PetscOptionsGetBool(NULL,NULL, "-Jac_view", &viewpipes,NULL);
803: if (viewpipes) {
804: SNES snes;
805: Mat Jac;
806: TSGetSNES(ts,&snes);
807: SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
808: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
809: }
811: /* View solutions */
812: /* -------------- */
813: viewpipes = PETSC_FALSE;
814: PetscOptionsGetBool(NULL,NULL, "-pipe_view", &viewpipes,NULL);
815: if (viewpipes) {
816: PipesView(networkdm,KeyPipe,X);
817: }
819: /* Test IS */
820: viewjuncs = PETSC_FALSE;
821: PetscOptionsGetBool(NULL,NULL, "-isjunc_view", &viewjuncs,NULL);
822: if (viewjuncs) {
823: ISJunctionsView(networkdm,KeyJunction);
824: }
826: /* Free spaces */
827: /* ----------- */
828: TSDestroy(&ts);
829: VecDestroy(&X);
830: VecDestroy(&wash->localX);
831: VecDestroy(&wash->localXdot);
833: /* Destroy objects from each pipe that are created in PipeSetUp() */
834: DMNetworkGetEdgeRange(networkdm,&eStart, &eEnd);
835: for (i = eStart; i < eEnd; i++) {
836: DMNetworkGetComponent(networkdm,i,0,&key,(void**)&pipe,NULL);
837: PipeDestroy(&pipe);
838: }
839: if (userJac) {
840: for (v=vStart; v<vEnd; v++) {
841: DMNetworkGetComponent(networkdm,v,0,&vkey,(void**)&junction,NULL);
842: JunctionDestroyJacobian(networkdm,v,junction);
843: }
844: }
846: if (size == 1 && monipipes) {
847: DMNetworkMonitorDestroy(&monitor);
848: }
849: DMDestroy(&networkdm);
850: PetscFree(wash);
852: if (rank) {
853: PetscFree2(junctions,pipes);
854: }
855: PetscFinalize();
856: return 0;
857: }
859: /*TEST
861: build:
862: depends: pipeInterface.c pipeImpls.c
863: requires: mumps
865: test:
866: args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -options_left no -viewX
867: localrunfiles: pOption
868: output_file: output/pipes_1.out
870: test:
871: suffix: 2
872: nsize: 2
873: args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
874: localrunfiles: pOption
875: output_file: output/pipes_2.out
877: test:
878: suffix: 3
879: nsize: 2
880: args: -ts_monitor -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
881: localrunfiles: pOption
882: output_file: output/pipes_3.out
884: test:
885: suffix: 4
886: args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -options_left no -viewX
887: localrunfiles: pOption
888: output_file: output/pipes_4.out
890: test:
891: suffix: 5
892: nsize: 3
893: args: -ts_monitor -case 2 -ts_max_steps 10 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
894: localrunfiles: pOption
895: output_file: output/pipes_5.out
897: test:
898: suffix: 6
899: nsize: 2
900: args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX
901: localrunfiles: pOption
902: output_file: output/pipes_6.out
904: test:
905: suffix: 7
906: nsize: 2
907: args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX
908: localrunfiles: pOption
909: output_file: output/pipes_7.out
911: test:
912: suffix: 8
913: nsize: 2
914: requires: parmetis
915: args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type parmetis -options_left no -wash_distribute 1
916: localrunfiles: pOption
917: output_file: output/pipes_8.out
919: test:
920: suffix: 9
921: nsize: 2
922: args: -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -pipe_view
923: localrunfiles: pOption
924: output_file: output/pipes_9.out
926: test:
927: suffix: 10
928: nsize: 2
929: args: -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -isjunc_view
930: localrunfiles: pOption
931: output_file: output/pipes_10.out
933: TEST*/