Actual source code: waterfunctions.c
1: /* function subroutines used by water.c */
3: #include "water.h"
4: #include <petscdmnetwork.h>
6: PetscScalar Flow_Pipe(Pipe *pipe, PetscScalar hf, PetscScalar ht)
7: {
8: PetscScalar flow_pipe;
10: flow_pipe = PetscSign(hf - ht) * PetscPowScalar(PetscAbsScalar(hf - ht) / pipe->k, (1 / pipe->n));
11: return flow_pipe;
12: }
14: PetscScalar Flow_Pump(Pump *pump, PetscScalar hf, PetscScalar ht)
15: {
16: PetscScalar flow_pump;
17: flow_pump = PetscPowScalar((hf - ht + pump->h0) / pump->r, (1 / pump->n));
18: return flow_pump;
19: }
21: PetscErrorCode FormFunction_Water(DM networkdm, Vec localX, Vec localF, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
22: {
23: const PetscScalar *xarr;
24: const PetscInt *cone;
25: PetscScalar *farr, hf, ht, flow;
26: PetscInt i, key, vnode1, vnode2, offsetnode1, offsetnode2, offset, ncomp;
27: PetscBool ghostvtex;
28: VERTEX_Water vertex, vertexnode1, vertexnode2;
29: EDGE_Water edge;
30: Pipe *pipe;
31: Pump *pump;
32: Reservoir *reservoir;
33: Tank *tank;
35: /* Get arrays for the vectors */
36: VecGetArrayRead(localX, &xarr);
37: VecGetArray(localF, &farr);
39: for (i = 0; i < ne; i++) { /* for each edge */
40: /* Get the offset and the key for the component for edge number e[i] */
41: DMNetworkGetComponent(networkdm, edges[i], 0, &key, (void **)&edge, NULL);
43: /* Get the numbers for the vertices covering this edge */
44: DMNetworkGetConnectedVertices(networkdm, edges[i], &cone);
45: vnode1 = cone[0];
46: vnode2 = cone[1];
48: /* Get the components at the two vertices, their variable offsets */
49: DMNetworkGetNumComponents(networkdm, vnode1, &ncomp);
50: DMNetworkGetComponent(networkdm, vnode1, ncomp - 1, &key, (void **)&vertexnode1, NULL);
51: DMNetworkGetLocalVecOffset(networkdm, vnode1, ncomp - 1, &offsetnode1);
53: DMNetworkGetNumComponents(networkdm, vnode2, &ncomp);
54: DMNetworkGetComponent(networkdm, vnode2, ncomp - 1, &key, (void **)&vertexnode2, NULL);
55: DMNetworkGetLocalVecOffset(networkdm, vnode2, ncomp - 1, &offsetnode2);
57: /* Variables at node1 and node 2 */
58: hf = xarr[offsetnode1];
59: ht = xarr[offsetnode2];
61: flow = 0.0;
62: if (edge->type == EDGE_TYPE_PIPE) {
63: pipe = &edge->pipe;
64: flow = Flow_Pipe(pipe, hf, ht);
65: } else if (edge->type == EDGE_TYPE_PUMP) {
66: pump = &edge->pump;
67: flow = Flow_Pump(pump, hf, ht);
68: }
69: /* Convention: Node 1 has outgoing flow and Node 2 has incoming flow */
70: if (vertexnode1->type == VERTEX_TYPE_JUNCTION) farr[offsetnode1] -= flow;
71: if (vertexnode2->type == VERTEX_TYPE_JUNCTION) farr[offsetnode2] += flow;
72: }
74: /* Subtract Demand flows at the vertices */
75: for (i = 0; i < nv; i++) {
76: DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex);
77: if (ghostvtex) continue;
79: DMNetworkGetLocalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset);
80: DMNetworkGetNumComponents(networkdm, vtx[i], &ncomp);
81: DMNetworkGetComponent(networkdm, vtx[i], ncomp - 1, &key, (void **)&vertex, NULL);
83: if (vertex->type == VERTEX_TYPE_JUNCTION) {
84: farr[offset] -= vertex->junc.demand;
85: } else if (vertex->type == VERTEX_TYPE_RESERVOIR) {
86: reservoir = &vertex->res;
87: farr[offset] = xarr[offset] - reservoir->head;
88: } else if (vertex->type == VERTEX_TYPE_TANK) {
89: tank = &vertex->tank;
90: farr[offset] = xarr[offset] - (tank->elev + tank->initlvl);
91: }
92: }
94: VecRestoreArrayRead(localX, &xarr);
95: VecRestoreArray(localF, &farr);
96: return 0;
97: }
99: PetscErrorCode WaterFormFunction(SNES snes, Vec X, Vec F, void *user)
100: {
101: DM networkdm;
102: Vec localX, localF;
103: const PetscInt *v, *e;
104: PetscInt nv, ne;
106: /* Get the DM attached with the SNES */
107: SNESGetDM(snes, &networkdm);
109: /* Get local vertices and edges */
110: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &v, &e);
112: /* Get local vectors */
113: DMGetLocalVector(networkdm, &localX);
114: DMGetLocalVector(networkdm, &localF);
116: /* Scatter values from global to local vector */
117: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
118: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
120: /* Initialize residual */
121: VecSet(F, 0.0);
122: VecSet(localF, 0.0);
124: FormFunction_Water(networkdm, localX, localF, nv, ne, v, e, NULL);
126: DMRestoreLocalVector(networkdm, &localX);
127: DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F);
128: DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F);
130: DMRestoreLocalVector(networkdm, &localF);
131: return 0;
132: }
134: PetscErrorCode WaterSetInitialGuess(DM networkdm, Vec X)
135: {
136: PetscInt nv, ne;
137: const PetscInt *vtx, *edges;
138: Vec localX;
140: DMGetLocalVector(networkdm, &localX);
142: VecSet(X, 0.0);
143: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
144: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
146: /* Get subnetwork */
147: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
148: SetInitialGuess_Water(networkdm, localX, nv, ne, vtx, edges, NULL);
150: DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X);
151: DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X);
152: DMRestoreLocalVector(networkdm, &localX);
153: return 0;
154: }
156: PetscErrorCode GetListofEdges_Water(WATERDATA *water, PetscInt *edgelist)
157: {
158: PetscInt i, j, node1, node2;
159: Pipe *pipe;
160: Pump *pump;
161: PetscBool netview = PETSC_FALSE;
163: PetscOptionsHasName(NULL, NULL, "-water_view", &netview);
164: for (i = 0; i < water->nedge; i++) {
165: if (water->edge[i].type == EDGE_TYPE_PIPE) {
166: pipe = &water->edge[i].pipe;
167: node1 = pipe->node1;
168: node2 = pipe->node2;
169: if (netview) PetscPrintf(PETSC_COMM_SELF, "edge %" PetscInt_FMT ", pipe v[%" PetscInt_FMT "] -> v[%" PetscInt_FMT "]\n", i, node1, node2);
170: } else {
171: pump = &water->edge[i].pump;
172: node1 = pump->node1;
173: node2 = pump->node2;
174: if (netview) PetscPrintf(PETSC_COMM_SELF, "edge %" PetscInt_FMT ", pump v[%" PetscInt_FMT "] -> v[%" PetscInt_FMT "]\n", i, node1, node2);
175: }
177: for (j = 0; j < water->nvertex; j++) {
178: if (water->vertex[j].id == node1) {
179: edgelist[2 * i] = j;
180: break;
181: }
182: }
184: for (j = 0; j < water->nvertex; j++) {
185: if (water->vertex[j].id == node2) {
186: edgelist[2 * i + 1] = j;
187: break;
188: }
189: }
190: }
191: return 0;
192: }
194: PetscErrorCode SetInitialGuess_Water(DM networkdm, Vec localX, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
195: {
196: PetscInt i, offset, key;
197: PetscBool ghostvtex, sharedv;
198: VERTEX_Water vertex;
199: PetscScalar *xarr;
201: VecGetArray(localX, &xarr);
202: for (i = 0; i < nv; i++) {
203: DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex);
204: DMNetworkIsSharedVertex(networkdm, vtx[i], &sharedv);
205: if (ghostvtex || sharedv) continue;
207: DMNetworkGetComponent(networkdm, vtx[i], 0, &key, (void **)&vertex, NULL);
208: DMNetworkGetLocalVecOffset(networkdm, vtx[i], 0, &offset);
209: if (vertex->type == VERTEX_TYPE_JUNCTION) {
210: xarr[offset] = 100;
211: } else if (vertex->type == VERTEX_TYPE_RESERVOIR) {
212: xarr[offset] = vertex->res.head;
213: } else {
214: xarr[offset] = vertex->tank.initlvl + vertex->tank.elev;
215: }
216: }
217: VecRestoreArray(localX, &xarr);
218: return 0;
219: }