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