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: PetscFunctionBegin;
26: PetscCallMPI(MPI_Comm_size(comm, &size));
27: if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);
29: PetscCallMPI(MPI_Comm_rank(comm, &rank));
30: numEdges = wash->nedge;
31: numVertices = wash->nvertex;
33: /* (1) all processes get global and local number of edges */
34: PetscCallMPI(MPI_Bcast(&numEdges, 1, MPIU_INT, 0, comm));
35: nedges = numEdges / size; /* local nedges */
36: if (rank == 0) nedges += numEdges - size * (numEdges / size);
37: wash->Nedge = numEdges;
38: wash->nedge = nedges;
39: /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d] nedges %d, numEdges %d\n",rank,nedges,numEdges)); */
41: PetscCall(PetscCalloc3(size + 1, &eowners, size, &nvtx, numVertices, &vtxDone));
42: PetscCallMPI(MPI_Allgather(&nedges, 1, MPIU_INT, eowners + 1, 1, MPIU_INT, PETSC_COMM_WORLD));
43: eowners[0] = 0;
44: for (i = 2; i <= size; i++) eowners[i] += eowners[i - 1];
46: estart = eowners[rank];
47: eend = eowners[rank + 1];
48: /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d] own lists row %d - %d\n",rank,estart,eend)); */
50: /* (2) distribute row block edgelist to all processors */
51: if (rank == 0) {
52: vtype = wash->vtype;
53: for (i = 1; i < size; i++) {
54: /* proc[0] sends edgelist to proc[i] */
55: PetscCallMPI(MPI_Send(edgelist + 2 * eowners[i], 2 * (eowners[i + 1] - eowners[i]), MPIU_INT, i, tag, comm));
57: /* proc[0] sends vtype to proc[i] */
58: PetscCallMPI(MPI_Send(vtype + 2 * eowners[i], 2 * (eowners[i + 1] - eowners[i]), MPIU_INT, i, tag, comm));
59: }
60: } else {
61: MPI_Status status;
62: PetscCall(PetscMalloc1(2 * (eend - estart), &vtype));
63: PetscCall(PetscMalloc1(2 * (eend - estart), &edgelist));
65: PetscCallMPI(MPI_Recv(edgelist, 2 * (eend - estart), MPIU_INT, 0, tag, comm, &status));
66: PetscCallMPI(MPI_Recv(vtype, 2 * (eend - estart), MPIU_INT, 0, tag, comm, &status));
67: }
69: wash->edgelist = edgelist;
71: /* (3) all processes get global and local number of vertices, without ghost vertices */
72: if (rank == 0) {
73: for (i = 0; i < size; i++) {
74: for (e = eowners[i]; e < eowners[i + 1]; e++) {
75: v = edgelist[2 * e];
76: if (!vtxDone[v]) {
77: nvtx[i]++;
78: vtxDone[v] = 1;
79: }
80: v = edgelist[2 * e + 1];
81: if (!vtxDone[v]) {
82: nvtx[i]++;
83: vtxDone[v] = 1;
84: }
85: }
86: }
87: }
88: PetscCallMPI(MPI_Bcast(&numVertices, 1, MPIU_INT, 0, PETSC_COMM_WORLD));
89: PetscCallMPI(MPI_Scatter(nvtx, 1, MPIU_INT, &nvertices, 1, MPIU_INT, 0, PETSC_COMM_WORLD));
90: PetscCall(PetscFree3(eowners, nvtx, vtxDone));
92: wash->Nvertex = numVertices;
93: wash->nvertex = nvertices;
94: wash->vtype = vtype;
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: PetscErrorCode WASHIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx)
99: {
100: Wash wash = (Wash)ctx;
101: DM networkdm;
102: Vec localX, localXdot, localF, localXold;
103: const PetscInt *cone;
104: PetscInt vfrom, vto, offsetfrom, offsetto, varoffset;
105: PetscInt v, vStart, vEnd, e, eStart, eEnd;
106: PetscInt nend, type;
107: PetscBool ghost;
108: PetscScalar *farr, *juncf, *pipef;
109: PetscReal dt;
110: Pipe pipe;
111: PipeField *pipex, *pipexdot, *juncx;
112: Junction junction;
113: DMDALocalInfo info;
114: const PetscScalar *xarr, *xdotarr, *xoldarr;
116: PetscFunctionBegin;
117: localX = wash->localX;
118: localXdot = wash->localXdot;
120: PetscCall(TSGetSolution(ts, &localXold));
121: PetscCall(TSGetDM(ts, &networkdm));
122: PetscCall(TSGetTimeStep(ts, &dt));
123: PetscCall(DMGetLocalVector(networkdm, &localF));
125: /* Set F and localF as zero */
126: PetscCall(VecSet(F, 0.0));
127: PetscCall(VecSet(localF, 0.0));
129: /* Update ghost values of locaX and locaXdot */
130: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
131: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
133: PetscCall(DMGlobalToLocalBegin(networkdm, Xdot, INSERT_VALUES, localXdot));
134: PetscCall(DMGlobalToLocalEnd(networkdm, Xdot, INSERT_VALUES, localXdot));
136: PetscCall(VecGetArrayRead(localX, &xarr));
137: PetscCall(VecGetArrayRead(localXdot, &xdotarr));
138: PetscCall(VecGetArrayRead(localXold, &xoldarr));
139: PetscCall(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: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
150: for (v = vStart; v < vEnd; v++) {
151: PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghost));
152: if (ghost) continue;
154: PetscCall(DMNetworkGetComponent(networkdm, v, 0, &type, (void **)&junction, NULL));
155: PetscCall(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: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
169: for (e = eStart; e < eEnd; e++) {
170: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &type, (void **)&pipe, NULL));
171: PetscCall(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: PetscCall(DMDAGetLocalInfo(pipe->da, &info));
183: PetscCall(PipeIFunctionLocal_Lax(&info, t, pipex, pipexdot, pipef, pipe));
185: /* Get boundary values from connected vertices */
186: PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
187: vfrom = cone[0]; /* local ordering */
188: vto = cone[1];
189: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom));
190: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto));
192: /* Evaluate upstream boundary */
193: PetscCall(DMNetworkGetComponent(networkdm, vfrom, 0, &type, (void **)&junction, NULL));
194: PetscCheck(junction->type == JUNCTION || junction->type == RESERVOIR, PETSC_COMM_SELF, PETSC_ERR_SUP, "junction type is not supported");
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: PetscCall(DMNetworkGetComponent(networkdm, vto, 0, &type, (void **)&junction, NULL));
203: PetscCheck(junction->type == JUNCTION || junction->type == VALVE, PETSC_COMM_SELF, PETSC_ERR_SUP, "junction type is not supported");
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: PetscCall(VecRestoreArrayRead(localX, &xarr));
213: PetscCall(VecRestoreArrayRead(localXdot, &xdotarr));
214: PetscCall(VecRestoreArray(localF, &farr));
216: PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
217: PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
218: PetscCall(DMRestoreLocalVector(networkdm, &localF));
219: /*
220: PetscCall(PetscPrintf(PETSC_COMM_WORLD("F:\n"));
221: PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD));
222: */
223: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
239: PetscCall(VecSet(X, 0.0));
240: PetscCall(DMGetLocalVector(networkdm, &localX));
241: PetscCall(VecGetArray(localX, &xarr));
243: /* Edge */
244: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
245: for (e = eStart; e < eEnd; e++) {
246: PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &varoffset));
247: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &type, (void **)&pipe, NULL));
249: /* set initial values for this pipe */
250: PetscCall(PipeComputeSteadyState(pipe, wash->Q0, wash->H0));
251: PetscCall(VecGetSize(pipe->x, &nx));
253: PetscCall(VecGetArrayRead(pipe->x, &xarray));
254: /* copy pipe->x to xarray */
255: for (k = 0; k < nx; k++) (xarr + varoffset)[k] = xarray[k];
257: /* set boundary values into vfrom and vto */
258: PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
259: vfrom = cone[0]; /* local ordering */
260: vto = cone[1];
261: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom));
262: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto));
264: /* if vform is a head vertex: */
265: PetscCall(DMNetworkGetComponent(networkdm, vfrom, 0, &vkey, (void **)&junction, NULL));
266: if (junction->type == RESERVOIR) (xarr + offsetfrom)[1] = wash->H0; /* 1st H */
268: /* if vto is an end vertex: */
269: PetscCall(DMNetworkGetComponent(networkdm, vto, 0, &vkey, (void **)&junction, NULL));
270: if (junction->type == VALVE) (xarr + offsetto)[0] = wash->QL; /* last Q */
271: PetscCall(VecRestoreArrayRead(pipe->x, &xarray));
272: }
274: PetscCall(VecRestoreArray(localX, &xarr));
275: PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
276: PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
277: PetscCall(DMRestoreLocalVector(networkdm, &localX));
279: #if 0
280: PetscInt N;
281: PetscCall(VecGetSize(X,&N));
282: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"initial solution %d:\n",N));
283: PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD));
284: #endif
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: PetscErrorCode TSDMNetworkMonitor(TS ts, PetscInt step, PetscReal t, Vec x, void *context)
289: {
290: DMNetworkMonitor monitor;
292: PetscFunctionBegin;
293: monitor = (DMNetworkMonitor)context;
294: PetscCall(DMNetworkMonitorView(monitor, x));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: PetscErrorCode PipesView(DM networkdm, PetscInt KeyPipe, Vec X)
299: {
300: PetscInt i, numkeys = 1, *blocksize, *numselectedvariable, **selectedvariables, n;
301: IS isfrom_q, isfrom_h, isfrom;
302: Vec Xto;
303: VecScatter ctx;
304: MPI_Comm comm;
306: PetscFunctionBegin;
307: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
309: /* 1. Create isfrom_q for q-variable of pipes */
310: PetscCall(PetscMalloc3(numkeys, &blocksize, numkeys, &numselectedvariable, numkeys, &selectedvariables));
311: for (i = 0; i < numkeys; i++) {
312: blocksize[i] = 2;
313: numselectedvariable[i] = 1;
314: PetscCall(PetscMalloc1(numselectedvariable[i], &selectedvariables[i]));
315: selectedvariables[i][0] = 0; /* q-variable */
316: }
317: PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyPipe, blocksize, numselectedvariable, selectedvariables, &isfrom_q));
319: /* 2. Create Xto and isto */
320: PetscCall(ISGetLocalSize(isfrom_q, &n));
321: PetscCall(VecCreate(comm, &Xto));
322: PetscCall(VecSetSizes(Xto, n, PETSC_DECIDE));
323: PetscCall(VecSetFromOptions(Xto));
325: /* 3. Create scatter */
326: PetscCall(VecScatterCreate(X, isfrom_q, Xto, NULL, &ctx));
328: /* 4. Scatter to Xq */
329: PetscCall(VecScatterBegin(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
330: PetscCall(VecScatterEnd(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
331: PetscCall(VecScatterDestroy(&ctx));
332: PetscCall(ISDestroy(&isfrom_q));
334: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Xq:\n"));
335: PetscCall(VecView(Xto, PETSC_VIEWER_STDOUT_WORLD));
337: /* 5. Create isfrom_h for h-variable of pipes; Create scatter; Scatter to Xh */
338: for (i = 0; i < numkeys; i++) selectedvariables[i][0] = 1; /* h-variable */
339: PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyPipe, blocksize, numselectedvariable, selectedvariables, &isfrom_h));
341: PetscCall(VecScatterCreate(X, isfrom_h, Xto, NULL, &ctx));
342: PetscCall(VecScatterBegin(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
343: PetscCall(VecScatterEnd(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
344: PetscCall(VecScatterDestroy(&ctx));
345: PetscCall(ISDestroy(&isfrom_h));
347: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Xh:\n"));
348: PetscCall(VecView(Xto, PETSC_VIEWER_STDOUT_WORLD));
349: PetscCall(VecDestroy(&Xto));
351: /* 6. Create isfrom for all pipe variables; Create scatter; Scatter to Xpipes */
352: for (i = 0; i < numkeys; i++) blocksize[i] = -1; /* select all the variables of the i-th component */
353: PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyPipe, blocksize, NULL, NULL, &isfrom));
354: PetscCall(ISDestroy(&isfrom));
355: PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyPipe, NULL, NULL, NULL, &isfrom));
357: PetscCall(ISGetLocalSize(isfrom, &n));
358: PetscCall(VecCreate(comm, &Xto));
359: PetscCall(VecSetSizes(Xto, n, PETSC_DECIDE));
360: PetscCall(VecSetFromOptions(Xto));
362: PetscCall(VecScatterCreate(X, isfrom, Xto, NULL, &ctx));
363: PetscCall(VecScatterBegin(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
364: PetscCall(VecScatterEnd(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
365: PetscCall(VecScatterDestroy(&ctx));
366: PetscCall(ISDestroy(&isfrom));
368: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Xpipes:\n"));
369: PetscCall(VecView(Xto, PETSC_VIEWER_STDOUT_WORLD));
371: /* 7. Free spaces */
372: for (i = 0; i < numkeys; i++) PetscCall(PetscFree(selectedvariables[i]));
373: PetscCall(PetscFree3(blocksize, numselectedvariable, selectedvariables));
374: PetscCall(VecDestroy(&Xto));
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: PetscErrorCode ISJunctionsView(DM networkdm, PetscInt KeyJunc)
379: {
380: PetscInt numkeys = 1;
381: IS isfrom;
382: MPI_Comm comm;
383: PetscMPIInt rank;
385: PetscFunctionBegin;
386: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
387: PetscCallMPI(MPI_Comm_rank(comm, &rank));
389: /* Create a global isfrom for all junction variables */
390: PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyJunc, NULL, NULL, NULL, &isfrom));
391: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ISJunctions:\n"));
392: PetscCall(ISView(isfrom, PETSC_VIEWER_STDOUT_WORLD));
393: PetscCall(ISDestroy(&isfrom));
395: /* Create a local isfrom for all junction variables */
396: PetscCall(DMNetworkCreateLocalIS(networkdm, numkeys, &KeyJunc, NULL, NULL, NULL, &isfrom));
397: if (rank == 0) {
398: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] ISLocalJunctions:\n", rank));
399: PetscCall(ISView(isfrom, PETSC_VIEWER_STDOUT_SELF));
400: }
401: PetscCall(ISDestroy(&isfrom));
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: PetscErrorCode WashNetworkCleanUp(Wash wash)
406: {
407: PetscMPIInt rank;
409: PetscFunctionBegin;
410: PetscCallMPI(MPI_Comm_rank(wash->comm, &rank));
411: PetscCall(PetscFree(wash->edgelist));
412: PetscCall(PetscFree(wash->vtype));
413: if (rank == 0) PetscCall(PetscFree2(wash->junction, wash->pipe));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: PetscErrorCode WashNetworkCreate(MPI_Comm comm, PetscInt pipesCase, Wash *wash_ptr)
418: {
419: PetscInt npipes;
420: PetscMPIInt rank;
421: Wash wash = NULL;
422: PetscInt i, numVertices, numEdges, *vtype;
423: PetscInt *edgelist;
424: Junction junctions = NULL;
425: Pipe pipes = NULL;
426: PetscBool washdist = PETSC_TRUE;
428: PetscFunctionBegin;
429: PetscCallMPI(MPI_Comm_rank(comm, &rank));
431: PetscCall(PetscCalloc1(1, &wash));
432: wash->comm = comm;
433: *wash_ptr = wash;
434: wash->Q0 = 0.477432; /* RESERVOIR */
435: wash->H0 = 150.0;
436: wash->HL = 143.488; /* VALVE */
437: wash->QL = 0.0;
438: wash->nnodes_loc = 0;
440: numVertices = 0;
441: numEdges = 0;
442: edgelist = NULL;
444: /* proc[0] creates a sequential wash and edgelist */
445: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setup pipesCase %" PetscInt_FMT "\n", pipesCase));
447: /* Set global number of pipes, edges, and junctions */
448: /*-------------------------------------------------*/
449: switch (pipesCase) {
450: case 0:
451: /* pipeCase 0: */
452: /* =================================================
453: (RESERVOIR) v0 --E0--> v1--E1--> v2 --E2-->v3 (VALVE)
454: ==================================================== */
455: npipes = 3;
456: PetscCall(PetscOptionsGetInt(NULL, NULL, "-npipes", &npipes, NULL));
457: wash->nedge = npipes;
458: wash->nvertex = npipes + 1;
460: /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
461: numVertices = 0;
462: numEdges = 0;
463: edgelist = NULL;
464: if (rank == 0) {
465: numVertices = wash->nvertex;
466: numEdges = wash->nedge;
468: PetscCall(PetscCalloc1(2 * numEdges, &edgelist));
469: for (i = 0; i < numEdges; i++) {
470: edgelist[2 * i] = i;
471: edgelist[2 * i + 1] = i + 1;
472: }
474: /* Add network components */
475: /*------------------------*/
476: PetscCall(PetscCalloc2(numVertices, &junctions, numEdges, &pipes));
478: /* vertex */
479: for (i = 0; i < numVertices; i++) {
480: junctions[i].id = i;
481: junctions[i].type = JUNCTION;
482: }
484: junctions[0].type = RESERVOIR;
485: junctions[numVertices - 1].type = VALVE;
486: }
487: break;
488: case 1:
489: /* pipeCase 1: */
490: /* ==========================
491: v2 (VALVE)
492: ^
493: |
494: E2
495: |
496: v0 --E0--> v3--E1--> v1
497: (RESERVOIR) (RESERVOIR)
498: ============================= */
499: npipes = 3;
500: wash->nedge = npipes;
501: wash->nvertex = npipes + 1;
503: /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
504: if (rank == 0) {
505: numVertices = wash->nvertex;
506: numEdges = wash->nedge;
508: PetscCall(PetscCalloc1(2 * numEdges, &edgelist));
509: edgelist[0] = 0;
510: edgelist[1] = 3; /* edge[0] */
511: edgelist[2] = 3;
512: edgelist[3] = 1; /* edge[1] */
513: edgelist[4] = 3;
514: edgelist[5] = 2; /* edge[2] */
516: /* Add network components */
517: /*------------------------*/
518: PetscCall(PetscCalloc2(numVertices, &junctions, numEdges, &pipes));
519: /* vertex */
520: for (i = 0; i < numVertices; i++) {
521: junctions[i].id = i;
522: junctions[i].type = JUNCTION;
523: }
525: junctions[0].type = RESERVOIR;
526: junctions[1].type = VALVE;
527: junctions[2].type = VALVE;
528: }
529: break;
530: case 2:
531: /* pipeCase 2: */
532: /* ==========================
533: (RESERVOIR) v2--> E2
534: |
535: v0 --E0--> v3--E1--> v1
536: (RESERVOIR) (VALVE)
537: ============================= */
539: /* Set application parameters -- to be used in function evaluations */
540: npipes = 3;
541: wash->nedge = npipes;
542: wash->nvertex = npipes + 1;
544: /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
545: if (rank == 0) {
546: numVertices = wash->nvertex;
547: numEdges = wash->nedge;
549: PetscCall(PetscCalloc1(2 * numEdges, &edgelist));
550: edgelist[0] = 0;
551: edgelist[1] = 3; /* edge[0] */
552: edgelist[2] = 3;
553: edgelist[3] = 1; /* edge[1] */
554: edgelist[4] = 2;
555: edgelist[5] = 3; /* edge[2] */
557: /* Add network components */
558: /*------------------------*/
559: PetscCall(PetscCalloc2(numVertices, &junctions, numEdges, &pipes));
560: /* vertex */
561: for (i = 0; i < numVertices; i++) {
562: junctions[i].id = i;
563: junctions[i].type = JUNCTION;
564: }
566: junctions[0].type = RESERVOIR;
567: junctions[1].type = VALVE;
568: junctions[2].type = RESERVOIR;
569: }
570: break;
571: default:
572: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "not done yet");
573: }
575: /* set edge global id */
576: for (i = 0; i < numEdges; i++) pipes[i].id = i;
578: if (rank == 0) { /* set vtype for proc[0] */
579: PetscInt v;
580: PetscCall(PetscMalloc1(2 * numEdges, &vtype));
581: for (i = 0; i < 2 * numEdges; i++) {
582: v = edgelist[i];
583: vtype[i] = junctions[v].type;
584: }
585: wash->vtype = vtype;
586: }
588: *wash_ptr = wash;
589: wash->nedge = numEdges;
590: wash->nvertex = numVertices;
591: wash->edgelist = edgelist;
592: wash->junction = junctions;
593: wash->pipe = pipes;
595: /* Distribute edgelist to other processors */
596: PetscCall(PetscOptionsGetBool(NULL, NULL, "-wash_distribute", &washdist, NULL));
597: if (washdist) {
598: /*
599: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Distribute sequential wash ...\n"));
600: */
601: PetscCall(WashNetworkDistribute(comm, wash));
602: }
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: /* ------------------------------------------------------- */
607: int main(int argc, char **argv)
608: {
609: Wash wash;
610: Junction junctions, junction;
611: Pipe pipe, pipes;
612: PetscInt KeyPipe, KeyJunction, *edgelist = NULL, *vtype = NULL;
613: PetscInt i, e, v, eStart, eEnd, vStart, vEnd, key, vkey, type;
614: PetscInt steps = 1, nedges, nnodes = 6;
615: const PetscInt *cone;
616: DM networkdm;
617: PetscMPIInt size, rank;
618: PetscReal ftime;
619: Vec X;
620: TS ts;
621: TSConvergedReason reason;
622: PetscBool viewpipes, viewjuncs, monipipes = PETSC_FALSE, userJac = PETSC_TRUE, viewdm = PETSC_FALSE, viewX = PETSC_FALSE;
623: PetscInt pipesCase = 0;
624: DMNetworkMonitor monitor;
625: MPI_Comm comm;
627: PetscFunctionBeginUser;
628: PetscCall(PetscInitialize(&argc, &argv, "pOption", help));
630: /* Read runtime options */
631: PetscCall(PetscOptionsGetInt(NULL, NULL, "-case", &pipesCase, NULL));
632: PetscCall(PetscOptionsGetBool(NULL, NULL, "-user_Jac", &userJac, NULL));
633: PetscCall(PetscOptionsGetBool(NULL, NULL, "-pipe_monitor", &monipipes, NULL));
634: PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewdm", &viewdm, NULL));
635: PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewX", &viewX, NULL));
636: PetscCall(PetscOptionsGetInt(NULL, NULL, "-npipenodes", &nnodes, NULL));
638: /* Create networkdm */
639: /*------------------*/
640: PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
641: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
642: PetscCallMPI(MPI_Comm_rank(comm, &rank));
643: PetscCallMPI(MPI_Comm_size(comm, &size));
645: if (size == 1 && monipipes) PetscCall(DMNetworkMonitorCreate(networkdm, &monitor));
647: /* Register the components in the network */
648: PetscCall(DMNetworkRegisterComponent(networkdm, "junctionstruct", sizeof(struct _p_Junction), &KeyJunction));
649: PetscCall(DMNetworkRegisterComponent(networkdm, "pipestruct", sizeof(struct _p_Pipe), &KeyPipe));
651: /* Create a distributed wash network (user-specific) */
652: PetscCall(WashNetworkCreate(comm, pipesCase, &wash));
653: nedges = wash->nedge;
654: edgelist = wash->edgelist;
655: vtype = wash->vtype;
656: junctions = wash->junction;
657: pipes = wash->pipe;
659: /* Set up the network layout */
660: PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1));
661: PetscCall(DMNetworkAddSubnetwork(networkdm, NULL, nedges, edgelist, NULL));
663: PetscCall(DMNetworkLayoutSetUp(networkdm));
665: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
666: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
667: /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d] eStart/End: %d - %d; vStart/End: %d - %d\n",rank,eStart,eEnd,vStart,vEnd)); */
669: if (rank) { /* junctions[] and pipes[] for proc[0] are allocated in WashNetworkCreate() */
670: /* vEnd - vStart = nvertices + number of ghost vertices! */
671: PetscCall(PetscCalloc2(vEnd - vStart, &junctions, nedges, &pipes));
672: }
674: /* Add Pipe component and number of variables to all local edges */
675: for (e = eStart; e < eEnd; e++) {
676: pipes[e - eStart].nnodes = nnodes;
677: PetscCall(DMNetworkAddComponent(networkdm, e, KeyPipe, &pipes[e - eStart], 2 * pipes[e - eStart].nnodes));
679: if (size == 1 && monipipes) { /* Add monitor -- show Q_{pipes[e-eStart].id}? */
680: pipes[e - eStart].length = 600.0;
681: PetscCall(DMNetworkMonitorAdd(monitor, "Pipe Q", e, pipes[e - eStart].nnodes, 0, 2, 0.0, pipes[e - eStart].length, -0.8, 0.8, PETSC_TRUE));
682: PetscCall(DMNetworkMonitorAdd(monitor, "Pipe H", e, pipes[e - eStart].nnodes, 1, 2, 0.0, pipes[e - eStart].length, -400.0, 800.0, PETSC_TRUE));
683: }
684: }
686: /* Add Junction component and number of variables to all local vertices */
687: for (v = vStart; v < vEnd; v++) PetscCall(DMNetworkAddComponent(networkdm, v, KeyJunction, &junctions[v - vStart], 2));
689: if (size > 1) { /* must be called before DMSetUp()???. Other partitioners do not work yet??? -- cause crash in proc[0]! */
690: DM plexdm;
691: PetscPartitioner part;
692: PetscCall(DMNetworkGetPlex(networkdm, &plexdm));
693: PetscCall(DMPlexGetPartitioner(plexdm, &part));
694: PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE));
695: PetscCall(PetscOptionsSetValue(NULL, "-dm_plex_csr_alg", "mat")); /* for parmetis */
696: }
698: /* Set up DM for use */
699: PetscCall(DMSetUp(networkdm));
700: if (viewdm) {
701: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nOriginal networkdm, DMView:\n"));
702: PetscCall(DMView(networkdm, PETSC_VIEWER_STDOUT_WORLD));
703: }
705: /* Set user physical parameters to the components */
706: for (e = eStart; e < eEnd; e++) {
707: PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
708: /* vfrom */
709: PetscCall(DMNetworkGetComponent(networkdm, cone[0], 0, &vkey, (void **)&junction, NULL));
710: junction->type = (VertexType)vtype[2 * e];
712: /* vto */
713: PetscCall(DMNetworkGetComponent(networkdm, cone[1], 0, &vkey, (void **)&junction, NULL));
714: junction->type = (VertexType)vtype[2 * e + 1];
715: }
717: PetscCall(WashNetworkCleanUp(wash));
719: /* Network partitioning and distribution of data */
720: PetscCall(DMNetworkDistribute(&networkdm, 0));
721: if (viewdm) {
722: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nAfter DMNetworkDistribute, DMView:\n"));
723: PetscCall(DMView(networkdm, PETSC_VIEWER_STDOUT_WORLD));
724: }
726: /* create vectors */
727: PetscCall(DMCreateGlobalVector(networkdm, &X));
728: PetscCall(DMCreateLocalVector(networkdm, &wash->localX));
729: PetscCall(DMCreateLocalVector(networkdm, &wash->localXdot));
731: /* PipeSetUp -- each process only sets its own pipes */
732: /*---------------------------------------------------*/
733: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
735: userJac = PETSC_TRUE;
736: PetscCall(DMNetworkHasJacobian(networkdm, userJac, userJac));
737: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
738: for (e = eStart; e < eEnd; e++) { /* each edge has only one component, pipe */
739: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &type, (void **)&pipe, NULL));
741: wash->nnodes_loc += pipe->nnodes; /* local total number of nodes, will be used by PipesView() */
742: PetscCall(PipeSetParameters(pipe, 600.0, /* length */
743: 0.5, /* diameter */
744: 1200.0, /* a */
745: 0.018)); /* friction */
746: PetscCall(PipeSetUp(pipe));
748: if (userJac) {
749: /* Create Jacobian matrix nonzero structures for a Pipe */
750: Mat *J;
751: PetscCall(PipeCreateJacobian(pipe, NULL, &J));
752: PetscCall(DMNetworkEdgeSetMatrix(networkdm, e, J));
753: }
754: }
756: if (userJac) {
757: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
758: for (v = vStart; v < vEnd; v++) {
759: Mat *J;
760: PetscCall(JunctionCreateJacobian(networkdm, v, NULL, &J));
761: PetscCall(DMNetworkVertexSetMatrix(networkdm, v, J));
763: PetscCall(DMNetworkGetComponent(networkdm, v, 0, &vkey, (void **)&junction, NULL));
764: junction->jacobian = J;
765: }
766: }
768: /* Setup solver */
769: /*--------------------------------------------------------*/
770: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
772: PetscCall(TSSetDM(ts, (DM)networkdm));
773: PetscCall(TSSetIFunction(ts, NULL, WASHIFunction, wash));
775: PetscCall(TSSetMaxSteps(ts, steps));
776: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
777: PetscCall(TSSetTimeStep(ts, 0.1));
778: PetscCall(TSSetType(ts, TSBEULER));
779: if (size == 1 && monipipes) PetscCall(TSMonitorSet(ts, TSDMNetworkMonitor, monitor, NULL));
780: PetscCall(TSSetFromOptions(ts));
782: PetscCall(WASHSetInitialSolution(networkdm, X, wash));
784: PetscCall(TSSolve(ts, X));
786: PetscCall(TSGetSolveTime(ts, &ftime));
787: PetscCall(TSGetStepNumber(ts, &steps));
788: PetscCall(TSGetConvergedReason(ts, &reason));
789: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));
790: if (viewX) PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
792: viewpipes = PETSC_FALSE;
793: PetscCall(PetscOptionsGetBool(NULL, NULL, "-Jac_view", &viewpipes, NULL));
794: if (viewpipes) {
795: SNES snes;
796: Mat Jac;
797: PetscCall(TSGetSNES(ts, &snes));
798: PetscCall(SNESGetJacobian(snes, &Jac, NULL, NULL, NULL));
799: PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
800: }
802: /* View solutions */
803: /* -------------- */
804: viewpipes = PETSC_FALSE;
805: PetscCall(PetscOptionsGetBool(NULL, NULL, "-pipe_view", &viewpipes, NULL));
806: if (viewpipes) PetscCall(PipesView(networkdm, KeyPipe, X));
808: /* Test IS */
809: viewjuncs = PETSC_FALSE;
810: PetscCall(PetscOptionsGetBool(NULL, NULL, "-isjunc_view", &viewjuncs, NULL));
811: if (viewjuncs) PetscCall(ISJunctionsView(networkdm, KeyJunction));
813: /* Free spaces */
814: /* ----------- */
815: PetscCall(TSDestroy(&ts));
816: PetscCall(VecDestroy(&X));
817: PetscCall(VecDestroy(&wash->localX));
818: PetscCall(VecDestroy(&wash->localXdot));
820: /* Destroy objects from each pipe that are created in PipeSetUp() */
821: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
822: for (i = eStart; i < eEnd; i++) {
823: PetscCall(DMNetworkGetComponent(networkdm, i, 0, &key, (void **)&pipe, NULL));
824: PetscCall(PipeDestroy(&pipe));
825: }
826: if (userJac) {
827: for (v = vStart; v < vEnd; v++) {
828: PetscCall(DMNetworkGetComponent(networkdm, v, 0, &vkey, (void **)&junction, NULL));
829: PetscCall(JunctionDestroyJacobian(networkdm, v, junction));
830: }
831: }
833: if (size == 1 && monipipes) PetscCall(DMNetworkMonitorDestroy(&monitor));
834: PetscCall(DMDestroy(&networkdm));
835: PetscCall(PetscFree(wash));
837: if (rank) PetscCall(PetscFree2(junctions, pipes));
838: PetscCall(PetscFinalize());
839: return 0;
840: }
842: /*TEST
844: build:
845: depends: pipeInterface.c pipeImpls.c
846: #TODO: bugs in DMNetwork causing segfault with __float128
847: requires: mumps !__float128
849: test:
850: args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -options_left no -viewX
851: localrunfiles: pOption
852: output_file: output/pipes_1.out
854: test:
855: suffix: 2
856: nsize: 2
857: args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
858: localrunfiles: pOption
859: output_file: output/pipes_2.out
861: test:
862: suffix: 3
863: nsize: 2
864: args: -ts_monitor -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
865: localrunfiles: pOption
866: output_file: output/pipes_3.out
868: test:
869: suffix: 4
870: args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -options_left no -viewX
871: localrunfiles: pOption
872: output_file: output/pipes_4.out
874: test:
875: suffix: 5
876: nsize: 3
877: args: -ts_monitor -case 2 -ts_max_steps 10 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
878: localrunfiles: pOption
879: output_file: output/pipes_5.out
881: test:
882: suffix: 6
883: nsize: 2
884: 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
885: localrunfiles: pOption
886: output_file: output/pipes_6.out
888: test:
889: suffix: 7
890: nsize: 2
891: 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
892: localrunfiles: pOption
893: output_file: output/pipes_7.out
895: test:
896: suffix: 8
897: nsize: 2
898: requires: parmetis
899: args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type parmetis -options_left no -wash_distribute 1
900: localrunfiles: pOption
901: output_file: output/pipes_8.out
903: test:
904: suffix: 9
905: nsize: 2
906: args: -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -pipe_view
907: localrunfiles: pOption
908: output_file: output/pipes_9.out
910: test:
911: suffix: 10
912: nsize: 2
913: args: -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -isjunc_view
914: localrunfiles: pOption
915: output_file: output/pipes_10.out
917: TEST*/