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: }