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 on Pipe
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; nz[1] = 2;
122: rows[0] = 0; rows[1] = 1;
123: cols[0] = 0; cols[1] = 1;
124: MatSeqAIJSetPreallocation(Jpipe[1],0,nz);
125: MatSetValues(Jpipe[1],2,rows,2,cols,aa,INSERT_VALUES);
126: MatAssemblyBegin(Jpipe[1],MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(Jpipe[1],MAT_FINAL_ASSEMBLY);
129: /* Jacobian for downstream vertex */
130: MatCreate(PETSC_COMM_SELF,&Jpipe[2]);
131: MatSetSizes(Jpipe[2],PETSC_DECIDE,PETSC_DECIDE,M,2);
132: MatSetFromOptions(Jpipe[2]);
133: MatSetOption(Jpipe[2],MAT_STRUCTURE_ONLY,PETSC_TRUE);
134: nz[0] = 0; nz[1] = 0; nz[M-2] = 2; nz[M-1] = 2;
135: rows[0] = M - 2; rows[1] = M - 1;
136: MatSeqAIJSetPreallocation(Jpipe[2],0,nz);
137: MatSetValues(Jpipe[2],2,rows,2,cols,aa,INSERT_VALUES);
138: MatAssemblyBegin(Jpipe[2],MAT_FINAL_ASSEMBLY);
139: MatAssemblyEnd(Jpipe[2],MAT_FINAL_ASSEMBLY);
141: PetscFree2(nz,aa);
143: *J = Jpipe;
144: pipe->jacobian = Jpipe;
145: return 0;
146: }
148: PetscErrorCode PipeDestroyJacobian(Pipe pipe)
149: {
150: Mat *Jpipe = pipe->jacobian;
151: PetscInt i;
153: if (Jpipe) {
154: for (i=0; i<3; i++) {
155: MatDestroy(&Jpipe[i]);
156: }
157: }
158: PetscFree(Jpipe);
159: return 0;
160: }
162: /*
163: JunctionCreateJacobian - Create Jacobian matrices for a vertex.
165: Collective on Pipe
167: Input Parameter:
168: + dm - the DMNetwork object
169: . v - vertex point
170: - Jin - Jacobian patterns created by JunctionCreateJacobianSample() for reuse
172: Output Parameter:
173: . J - array of Jacobian matrices (see dmnetworkimpl.h)
175: Level: beginner
176: */
177: PetscErrorCode JunctionCreateJacobian(DM dm,PetscInt v,Mat *Jin,Mat *J[])
178: {
179: Mat *Jv;
180: PetscInt nedges,e,i,M,N,*rows,*cols;
181: PetscBool isSelf;
182: const PetscInt *edges,*cone;
183: PetscScalar *zeros;
185: /* Get array size of Jv */
186: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
189: /* two Jacobians for each connected edge: J(v,e) and J(v,vc); adding J(v,v), total 2*nedges+1 Jacobians */
190: PetscCalloc1(2*nedges+1,&Jv);
192: /* Create dense zero block for this vertex: J[0] = Jacobian(v,v) */
193: DMNetworkGetComponent(dm,v,-1,NULL,NULL,&M);
195: PetscMalloc3(M,&rows,M,&cols,M*M,&zeros);
196: PetscArrayzero(zeros,M*M);
197: for (i=0; i<M; i++) rows[i] = i;
199: for (e=0; e<nedges; e++) {
200: /* create Jv[2*e+1] = Jacobian(v,e), e: supporting edge */
201: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
202: isSelf = (v == cone[0]) ? PETSC_TRUE:PETSC_FALSE;
204: if (Jin) {
205: if (isSelf) {
206: Jv[2*e+1] = Jin[0];
207: } else {
208: Jv[2*e+1] = Jin[1];
209: }
210: Jv[2*e+2] = Jin[2];
211: PetscObjectReference((PetscObject)(Jv[2*e+1]));
212: PetscObjectReference((PetscObject)(Jv[2*e+2]));
213: } else {
214: /* create J(v,e) */
215: MatCreate(PETSC_COMM_SELF,&Jv[2*e+1]);
216: DMNetworkGetComponent(dm,edges[e],-1,NULL,NULL,&N);
217: MatSetSizes(Jv[2*e+1],PETSC_DECIDE,PETSC_DECIDE,M,N);
218: MatSetFromOptions(Jv[2*e+1]);
219: MatSetOption(Jv[2*e+1],MAT_STRUCTURE_ONLY,PETSC_TRUE);
220: MatSeqAIJSetPreallocation(Jv[2*e+1],2,NULL);
221: if (N) {
222: if (isSelf) { /* coupling at upstream */
223: for (i=0; i<2; i++) cols[i] = i;
224: } else { /* coupling at downstream */
225: cols[0] = N-2; cols[1] = N-1;
226: }
227: MatSetValues(Jv[2*e+1],2,rows,2,cols,zeros,INSERT_VALUES);
228: }
229: MatAssemblyBegin(Jv[2*e+1],MAT_FINAL_ASSEMBLY);
230: MatAssemblyEnd(Jv[2*e+1],MAT_FINAL_ASSEMBLY);
232: /* create Jv[2*e+2] = Jacobian(v,vc), vc: connected vertex.
233: In WashNetwork, v and vc are not connected, thus Jacobian(v,vc) is empty */
234: MatCreate(PETSC_COMM_SELF,&Jv[2*e+2]);
235: MatSetSizes(Jv[2*e+2],PETSC_DECIDE,PETSC_DECIDE,M,M); /* empty matrix, sizes can be arbitrary */
236: MatSetFromOptions(Jv[2*e+2]);
237: MatSetOption(Jv[2*e+2],MAT_STRUCTURE_ONLY,PETSC_TRUE);
238: MatSeqAIJSetPreallocation(Jv[2*e+2],1,NULL);
239: MatAssemblyBegin(Jv[2*e+2],MAT_FINAL_ASSEMBLY);
240: MatAssemblyEnd(Jv[2*e+2],MAT_FINAL_ASSEMBLY);
241: }
242: }
243: PetscFree3(rows,cols,zeros);
245: *J = Jv;
246: return 0;
247: }
249: PetscErrorCode JunctionDestroyJacobian(DM dm,PetscInt v,Junction junc)
250: {
251: Mat *Jv=junc->jacobian;
252: const PetscInt *edges;
253: PetscInt nedges,e;
255: if (!Jv) return 0;
257: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
258: for (e=0; e<nedges; e++) {
259: MatDestroy(&Jv[2*e+1]);
260: MatDestroy(&Jv[2*e+2]);
261: }
262: PetscFree(Jv);
263: return 0;
264: }