Actual source code: pipeInterface.c
1: #include "wash.h"
3: /* Subroutines for Pipe */
4: /* -------------------------------------------------------*/
6: /*
7: PipeCreate - Create Pipe object.
9: Input Parameters:
10: comm - MPI communicator
12: Output Parameter:
13: . pipe - location to put the PIPE context
14: */
15: PetscErrorCode PipeCreate(MPI_Comm comm, Pipe *pipe)
16: {
17: PetscFunctionBegin;
18: PetscCall(PetscNew(pipe));
19: PetscFunctionReturn(PETSC_SUCCESS);
20: }
22: /*
23: PipeDestroy - Destroy Pipe object.
25: Input Parameters:
26: pipe - Reference to pipe intended to be destroyed.
27: */
28: PetscErrorCode PipeDestroy(Pipe *pipe)
29: {
30: PetscFunctionBegin;
31: if (!*pipe) PetscFunctionReturn(PETSC_SUCCESS);
33: PetscCall(PipeDestroyJacobian(*pipe));
34: PetscCall(VecDestroy(&(*pipe)->x));
35: PetscCall(DMDestroy(&(*pipe)->da));
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: /*
40: PipeSetParameters - Set parameters for Pipe context
42: Input Parameter:
43: + pipe - PIPE object
44: . length -
45: . nnodes -
46: . D -
47: . a -
48: - fric -
49: */
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: }
60: /*
61: PipeSetUp - Set up pipe based on set parameters.
62: */
63: PetscErrorCode PipeSetUp(Pipe pipe)
64: {
65: DMDALocalInfo info;
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));
76: PetscCall(DMDAGetLocalInfo(pipe->da, &info));
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: }
84: /*
85: PipeCreateJacobian - Create Jacobian matrix structures for a Pipe.
87: Collective
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
93: Output Parameter:
94: . J - array of three empty Jacobian matrices
96: Level: beginner
97: */
98: PetscErrorCode PipeCreateJacobian(Pipe pipe, Mat *Jin, Mat *J[])
99: {
100: Mat *Jpipe;
101: PetscInt M, rows[2], cols[2], *nz;
102: PetscScalar *aa;
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));
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));
118: /* Jacobian for upstream vertex */
119: PetscCall(MatGetSize(Jpipe[0], &M, NULL));
120: PetscCall(PetscCalloc2(M, &nz, 4, &aa));
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));
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));
153: PetscCall(PetscFree2(nz, aa));
155: *J = Jpipe;
156: pipe->jacobian = Jpipe;
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: PetscErrorCode PipeDestroyJacobian(Pipe pipe)
161: {
162: Mat *Jpipe = pipe->jacobian;
163: PetscInt i;
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: }
173: /*
174: JunctionCreateJacobian - Create Jacobian matrices for a vertex.
176: Collective
178: Input Parameter:
179: + dm - the DMNetwork object
180: . v - vertex point
181: - Jin - Jacobian patterns created by JunctionCreateJacobianSample() for reuse
183: Output Parameter:
184: . J - array of Jacobian matrices (see dmnetworkimpl.h)
186: Level: beginner
187: */
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;
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);
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));
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;
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;
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));
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));
258: *J = Jv;
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: PetscErrorCode JunctionDestroyJacobian(DM dm, PetscInt v, Junction junc)
263: {
264: Mat *Jv = junc->jacobian;
265: const PetscInt *edges;
266: PetscInt nedges, e;
268: PetscFunctionBegin;
269: if (!Jv) PetscFunctionReturn(PETSC_SUCCESS);
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: }