1 #include "wash.h"
2
3 /* Subroutines for Pipe */
4 /* -------------------------------------------------------*/
5
6 /*
7 PipeCreate - Create Pipe object.
8
9 Input Parameters:
10 comm - MPI communicator
11
12 Output Parameter:
13 . pipe - location to put the PIPE context
14 */
PipeCreate(MPI_Comm comm,Pipe * pipe)15 PetscErrorCode PipeCreate(MPI_Comm comm, Pipe *pipe)
16 {
17 PetscFunctionBegin;
18 PetscCall(PetscNew(pipe));
19 PetscFunctionReturn(PETSC_SUCCESS);
20 }
21
22 /*
23 PipeDestroy - Destroy Pipe object.
24
25 Input Parameters:
26 pipe - Reference to pipe intended to be destroyed.
27 */
PipeDestroy(Pipe * pipe)28 PetscErrorCode PipeDestroy(Pipe *pipe)
29 {
30 PetscFunctionBegin;
31 if (!*pipe) PetscFunctionReturn(PETSC_SUCCESS);
32
33 PetscCall(PipeDestroyJacobian(*pipe));
34 PetscCall(VecDestroy(&(*pipe)->x));
35 PetscCall(DMDestroy(&(*pipe)->da));
36 PetscFunctionReturn(PETSC_SUCCESS);
37 }
38
39 /*
40 PipeSetParameters - Set parameters for Pipe context
41
42 Input Parameter:
43 + pipe - PIPE object
44 . length -
45 . nnodes -
46 . D -
47 . a -
48 - fric -
49 */
PipeSetParameters(Pipe pipe,PetscReal length,PetscReal D,PetscReal a,PetscReal fric)50 PetscErrorCode PipeSetParameters(Pipe pipe, PetscReal length, PetscReal D, PetscReal a, PetscReal fric)
51 {
52 PetscFunctionBegin;
53 pipe->length = length;
54 pipe->D = D;
55 pipe->a = a;
56 pipe->fric = fric;
57 PetscFunctionReturn(PETSC_SUCCESS);
58 }
59
60 /*
61 PipeSetUp - Set up pipe based on set parameters.
62 */
PipeSetUp(Pipe pipe)63 PetscErrorCode PipeSetUp(Pipe pipe)
64 {
65 DMDALocalInfo info;
66
67 PetscFunctionBegin;
68 PetscCall(DMDACreate1d(PETSC_COMM_SELF, DM_BOUNDARY_GHOSTED, pipe->nnodes, 2, 1, NULL, &pipe->da));
69 PetscCall(DMSetFromOptions(pipe->da));
70 PetscCall(DMSetUp(pipe->da));
71 PetscCall(DMDASetFieldName(pipe->da, 0, "Q"));
72 PetscCall(DMDASetFieldName(pipe->da, 1, "H"));
73 PetscCall(DMDASetUniformCoordinates(pipe->da, 0, pipe->length, 0, 0, 0, 0));
74 PetscCall(DMCreateGlobalVector(pipe->da, &pipe->x));
75
76 PetscCall(DMDAGetLocalInfo(pipe->da, &info));
77
78 pipe->rad = pipe->D / 2;
79 pipe->A = PETSC_PI * pipe->rad * pipe->rad;
80 pipe->R = pipe->fric / (2 * pipe->D * pipe->A);
81 PetscFunctionReturn(PETSC_SUCCESS);
82 }
83
84 /*
85 PipeCreateJacobian - Create Jacobian matrix nonzero structures for a Pipe.
86
87 Collective
88
89 Input Parameter:
90 + pipe - the Pipe object
91 - Jin - array of three constructed Jacobian matrices to be reused. Set NULL if it is not available
92
93 Output Parameter:
94 . J - array of three empty Jacobian matrices
95
96 Level: beginner
97 */
PipeCreateJacobian(Pipe pipe,Mat * Jin,Mat * J[])98 PetscErrorCode PipeCreateJacobian(Pipe pipe, Mat *Jin, Mat *J[])
99 {
100 Mat *Jpipe;
101 PetscInt M, rows[2], cols[2], *nz;
102 PetscScalar *aa;
103
104 PetscFunctionBegin;
105 if (Jin) {
106 *J = Jin;
107 pipe->jacobian = Jin;
108 PetscCall(PetscObjectReference((PetscObject)Jin[0]));
109 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111 PetscCall(PetscMalloc1(3, &Jpipe));
112
113 /* Jacobian for this pipe */
114 PetscCall(DMSetMatrixStructureOnly(pipe->da, PETSC_TRUE));
115 PetscCall(DMCreateMatrix(pipe->da, &Jpipe[0]));
116 PetscCall(DMSetMatrixStructureOnly(pipe->da, PETSC_FALSE));
117
118 /* Jacobian for upstream vertex */
119 PetscCall(MatGetSize(Jpipe[0], &M, NULL));
120 PetscCall(PetscCalloc2(M, &nz, 4, &aa));
121
122 PetscCall(MatCreate(PETSC_COMM_SELF, &Jpipe[1]));
123 PetscCall(MatSetSizes(Jpipe[1], PETSC_DECIDE, PETSC_DECIDE, M, 2));
124 PetscCall(MatSetFromOptions(Jpipe[1]));
125 PetscCall(MatSetOption(Jpipe[1], MAT_STRUCTURE_ONLY, PETSC_TRUE));
126 nz[0] = 2;
127 nz[1] = 2;
128 rows[0] = 0;
129 rows[1] = 1;
130 cols[0] = 0;
131 cols[1] = 1;
132 PetscCall(MatSeqAIJSetPreallocation(Jpipe[1], 0, nz));
133 PetscCall(MatSetValues(Jpipe[1], 2, rows, 2, cols, aa, INSERT_VALUES));
134 PetscCall(MatAssemblyBegin(Jpipe[1], MAT_FINAL_ASSEMBLY));
135 PetscCall(MatAssemblyEnd(Jpipe[1], MAT_FINAL_ASSEMBLY));
136
137 /* Jacobian for downstream vertex */
138 PetscCall(MatCreate(PETSC_COMM_SELF, &Jpipe[2]));
139 PetscCall(MatSetSizes(Jpipe[2], PETSC_DECIDE, PETSC_DECIDE, M, 2));
140 PetscCall(MatSetFromOptions(Jpipe[2]));
141 PetscCall(MatSetOption(Jpipe[2], MAT_STRUCTURE_ONLY, PETSC_TRUE));
142 nz[0] = 0;
143 nz[1] = 0;
144 nz[M - 2] = 2;
145 nz[M - 1] = 2;
146 rows[0] = M - 2;
147 rows[1] = M - 1;
148 PetscCall(MatSeqAIJSetPreallocation(Jpipe[2], 0, nz));
149 PetscCall(MatSetValues(Jpipe[2], 2, rows, 2, cols, aa, INSERT_VALUES));
150 PetscCall(MatAssemblyBegin(Jpipe[2], MAT_FINAL_ASSEMBLY));
151 PetscCall(MatAssemblyEnd(Jpipe[2], MAT_FINAL_ASSEMBLY));
152
153 PetscCall(PetscFree2(nz, aa));
154
155 *J = Jpipe;
156 pipe->jacobian = Jpipe;
157 PetscFunctionReturn(PETSC_SUCCESS);
158 }
159
PipeDestroyJacobian(Pipe pipe)160 PetscErrorCode PipeDestroyJacobian(Pipe pipe)
161 {
162 Mat *Jpipe = pipe->jacobian;
163 PetscInt i;
164
165 PetscFunctionBegin;
166 if (Jpipe) {
167 for (i = 0; i < 3; i++) PetscCall(MatDestroy(&Jpipe[i]));
168 }
169 PetscCall(PetscFree(Jpipe));
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
173 /*
174 JunctionCreateJacobian - Create Jacobian matrices for a vertex.
175
176 Collective
177
178 Input Parameter:
179 + dm - the DMNetwork object
180 . v - vertex point
181 - Jin - Jacobian patterns created by JunctionCreateJacobianSample() for reuse
182
183 Output Parameter:
184 . J - array of Jacobian matrices (see dmnetworkimpl.h)
185
186 Level: beginner
187 */
JunctionCreateJacobian(DM dm,PetscInt v,Mat * Jin,Mat * J[])188 PetscErrorCode JunctionCreateJacobian(DM dm, PetscInt v, Mat *Jin, Mat *J[])
189 {
190 Mat *Jv;
191 PetscInt nedges, e, i, M, N, *rows, *cols;
192 PetscBool isSelf;
193 const PetscInt *edges, *cone;
194 PetscScalar *zeros;
195
196 PetscFunctionBegin;
197 /* Get array size of Jv */
198 PetscCall(DMNetworkGetSupportingEdges(dm, v, &nedges, &edges));
199 PetscCheck(nedges > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "%" PetscInt_FMT " vertex, nedges %" PetscInt_FMT, v, nedges);
200
201 /* two Jacobians for each connected edge: J(v,e) and J(v,vc); adding J(v,v), total 2*nedges+1 Jacobians */
202 PetscCall(PetscCalloc1(2 * nedges + 1, &Jv));
203
204 /* Create dense zero block for this vertex: J[0] = Jacobian(v,v) */
205 PetscCall(DMNetworkGetComponent(dm, v, -1, NULL, NULL, &M));
206 PetscCheck(M == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "%" PetscInt_FMT " != 2", M);
207 PetscCall(PetscMalloc3(M, &rows, M, &cols, M * M, &zeros));
208 PetscCall(PetscArrayzero(zeros, M * M));
209 for (i = 0; i < M; i++) rows[i] = i;
210
211 for (e = 0; e < nedges; e++) {
212 /* create Jv[2*e+1] = Jacobian(v,e), e: supporting edge */
213 PetscCall(DMNetworkGetConnectedVertices(dm, edges[e], &cone));
214 isSelf = (v == cone[0]) ? PETSC_TRUE : PETSC_FALSE;
215
216 if (Jin) {
217 if (isSelf) {
218 Jv[2 * e + 1] = Jin[0];
219 } else {
220 Jv[2 * e + 1] = Jin[1];
221 }
222 Jv[2 * e + 2] = Jin[2];
223 PetscCall(PetscObjectReference((PetscObject)Jv[2 * e + 1]));
224 PetscCall(PetscObjectReference((PetscObject)Jv[2 * e + 2]));
225 } else {
226 /* create J(v,e) */
227 PetscCall(MatCreate(PETSC_COMM_SELF, &Jv[2 * e + 1]));
228 PetscCall(DMNetworkGetComponent(dm, edges[e], -1, NULL, NULL, &N));
229 PetscCall(MatSetSizes(Jv[2 * e + 1], PETSC_DECIDE, PETSC_DECIDE, M, N));
230 PetscCall(MatSetFromOptions(Jv[2 * e + 1]));
231 PetscCall(MatSetOption(Jv[2 * e + 1], MAT_STRUCTURE_ONLY, PETSC_TRUE));
232 PetscCall(MatSeqAIJSetPreallocation(Jv[2 * e + 1], 2, NULL));
233 if (N) {
234 if (isSelf) { /* coupling at upstream */
235 for (i = 0; i < 2; i++) cols[i] = i;
236 } else { /* coupling at downstream */
237 cols[0] = N - 2;
238 cols[1] = N - 1;
239 }
240 PetscCall(MatSetValues(Jv[2 * e + 1], 2, rows, 2, cols, zeros, INSERT_VALUES));
241 }
242 PetscCall(MatAssemblyBegin(Jv[2 * e + 1], MAT_FINAL_ASSEMBLY));
243 PetscCall(MatAssemblyEnd(Jv[2 * e + 1], MAT_FINAL_ASSEMBLY));
244
245 /* create Jv[2*e+2] = Jacobian(v,vc), vc: connected vertex.
246 In WashNetwork, v and vc are not connected, thus Jacobian(v,vc) is empty */
247 PetscCall(MatCreate(PETSC_COMM_SELF, &Jv[2 * e + 2]));
248 PetscCall(MatSetSizes(Jv[2 * e + 2], PETSC_DECIDE, PETSC_DECIDE, M, M)); /* empty matrix, sizes can be arbitrary */
249 PetscCall(MatSetFromOptions(Jv[2 * e + 2]));
250 PetscCall(MatSetOption(Jv[2 * e + 2], MAT_STRUCTURE_ONLY, PETSC_TRUE));
251 PetscCall(MatSeqAIJSetPreallocation(Jv[2 * e + 2], 1, NULL));
252 PetscCall(MatAssemblyBegin(Jv[2 * e + 2], MAT_FINAL_ASSEMBLY));
253 PetscCall(MatAssemblyEnd(Jv[2 * e + 2], MAT_FINAL_ASSEMBLY));
254 }
255 }
256 PetscCall(PetscFree3(rows, cols, zeros));
257
258 *J = Jv;
259 PetscFunctionReturn(PETSC_SUCCESS);
260 }
261
JunctionDestroyJacobian(DM dm,PetscInt v,Junction junc)262 PetscErrorCode JunctionDestroyJacobian(DM dm, PetscInt v, Junction junc)
263 {
264 Mat *Jv = junc->jacobian;
265 const PetscInt *edges;
266 PetscInt nedges, e;
267
268 PetscFunctionBegin;
269 if (!Jv) PetscFunctionReturn(PETSC_SUCCESS);
270
271 PetscCall(DMNetworkGetSupportingEdges(dm, v, &nedges, &edges));
272 for (e = 0; e < nedges; e++) {
273 PetscCall(MatDestroy(&Jv[2 * e + 1]));
274 PetscCall(MatDestroy(&Jv[2 * e + 2]));
275 }
276 PetscCall(PetscFree(Jv));
277 PetscFunctionReturn(PETSC_SUCCESS);
278 }
279