Actual source code: pipeInterface.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  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: {

 20:   PetscNew(pipe);
 21:   return(0);
 22: }

 24: /*
 25:    PipeDestroy - Destroy Pipe object.

 27:    Input Parameters:
 28:    pipe - Reference to pipe intended to be destroyed.
 29: */
 30: PetscErrorCode PipeDestroy(Pipe *pipe)
 31: {

 35:   if (!*pipe) return(0);

 37:   PipeDestroyJacobian(*pipe);
 38:   VecDestroy(&(*pipe)->x);
 39:   DMDestroy(&(*pipe)->da);
 40:   return(0);
 41: }

 43: /*
 44:    PipeSetParameters - Set parameters for Pipe context

 46:    Input Parameter:
 47: +  pipe - PIPE object
 48: .  length -
 49: .  nnodes -
 50: .  D -
 51: .  a -
 52: -  fric -
 53: */
 54: PetscErrorCode PipeSetParameters(Pipe pipe,PetscReal length,PetscReal D,PetscReal a,PetscReal fric)
 55: {
 57:   pipe->length = length;
 58:   pipe->D      = D;
 59:   pipe->a      = a;
 60:   pipe->fric   = fric;
 61:   return(0);
 62: }

 64: /*
 65:     PipeSetUp - Set up pipe based on set parameters.
 66: */
 67: PetscErrorCode PipeSetUp(Pipe pipe)
 68: {
 69:   DMDALocalInfo  info;

 73:   DMDACreate1d(PETSC_COMM_SELF, DM_BOUNDARY_GHOSTED, pipe->nnodes, 2, 1, NULL, &pipe->da);
 74:   DMSetFromOptions(pipe->da);
 75:   DMSetUp(pipe->da);
 76:   DMDASetFieldName(pipe->da, 0, "Q");
 77:   DMDASetFieldName(pipe->da, 1, "H");
 78:   DMDASetUniformCoordinates(pipe->da, 0, pipe->length, 0, 0, 0, 0);
 79:   DMCreateGlobalVector(pipe->da, &(pipe->x));

 81:   DMDAGetLocalInfo(pipe->da, &info);

 83:   pipe->rad = pipe->D / 2;
 84:   pipe->A   = PETSC_PI*pipe->rad*pipe->rad;
 85:   pipe->R   = pipe->fric / (2*pipe->D*pipe->A);
 86:   return(0);
 87: }

 89: /*
 90:     PipeCreateJacobian - Create Jacobian matrix structures for a Pipe.

 92:     Collective on Pipe

 94:     Input Parameter:
 95: +   pipe - the Pipe object
 96: -   Jin - array of three constructed Jacobian matrices to be reused. Set NULL if it is not available

 98:     Output Parameter:
 99: .   J  - array of three empty Jacobian matrices

101:     Level: beginner
102: */
103: PetscErrorCode PipeCreateJacobian(Pipe pipe,Mat *Jin,Mat *J[])
104: {
106:   Mat            *Jpipe;
107:   PetscInt       M,rows[2],cols[2],*nz;
108:   PetscScalar    *aa;

111:   if (Jin) {
112:     *J = Jin;
113:     pipe->jacobian = Jin;
114:     PetscObjectReference((PetscObject)(Jin[0]));
115:     return(0);
116:   }
117:   PetscMalloc1(3,&Jpipe);

119:   /* Jacobian for this pipe */
120:   DMSetMatrixStructureOnly(pipe->da,PETSC_TRUE);
121:   DMCreateMatrix(pipe->da,&Jpipe[0]);
122:   DMSetMatrixStructureOnly(pipe->da,PETSC_FALSE);

124:   /* Jacobian for upstream vertex */
125:   MatGetSize(Jpipe[0],&M,NULL);
126:   PetscCalloc2(M,&nz,4,&aa);

128:   MatCreate(PETSC_COMM_SELF,&Jpipe[1]);
129:   MatSetSizes(Jpipe[1],PETSC_DECIDE,PETSC_DECIDE,M,2);
130:   MatSetFromOptions(Jpipe[1]);
131:   MatSetOption(Jpipe[1],MAT_STRUCTURE_ONLY,PETSC_TRUE);
132:   nz[0] = 2; nz[1] = 2;
133:   rows[0] = 0; rows[1] = 1;
134:   cols[0] = 0; cols[1] = 1;
135:   MatSeqAIJSetPreallocation(Jpipe[1],0,nz);
136:   MatSetValues(Jpipe[1],2,rows,2,cols,aa,INSERT_VALUES);
137:   MatAssemblyBegin(Jpipe[1],MAT_FINAL_ASSEMBLY);
138:   MatAssemblyEnd(Jpipe[1],MAT_FINAL_ASSEMBLY);

140:   /* Jacobian for downstream vertex */
141:   MatCreate(PETSC_COMM_SELF,&Jpipe[2]);
142:   MatSetSizes(Jpipe[2],PETSC_DECIDE,PETSC_DECIDE,M,2);
143:   MatSetFromOptions(Jpipe[2]);
144:   MatSetOption(Jpipe[2],MAT_STRUCTURE_ONLY,PETSC_TRUE);
145:   nz[0] = 0; nz[1] = 0; nz[M-2] = 2; nz[M-1] = 2;
146:   rows[0] = M - 2; rows[1] = M - 1;
147:   MatSeqAIJSetPreallocation(Jpipe[2],0,nz);
148:   MatSetValues(Jpipe[2],2,rows,2,cols,aa,INSERT_VALUES);
149:   MatAssemblyBegin(Jpipe[2],MAT_FINAL_ASSEMBLY);
150:   MatAssemblyEnd(Jpipe[2],MAT_FINAL_ASSEMBLY);

152:   PetscFree2(nz,aa);

154:   *J = Jpipe;
155:   pipe->jacobian = Jpipe;
156:   return(0);
157: }

159: PetscErrorCode PipeDestroyJacobian(Pipe pipe)
160: {
162:   Mat            *Jpipe = pipe->jacobian;
163:   PetscInt       i;

166:   if (Jpipe) {
167:     for (i=0; i<3; i++) {
168:       MatDestroy(&Jpipe[i]);
169:     }
170:   }
171:   PetscFree(Jpipe);
172:   return(0);
173: }

175: /*
176:     JunctionCreateJacobian - Create Jacobian matrices for a vertex.

178:     Collective on Pipe

180:     Input Parameter:
181: +   dm - the DMNetwork object
182: .   v - vertex point
183: -   Jin - Jacobian patterns created by JunctionCreateJacobianSample() for reuse

185:     Output Parameter:
186: .   J  - array of Jacobian matrices (see dmnetworkimpl.h)

188:     Level: beginner
189: */
190: PetscErrorCode JunctionCreateJacobian(DM dm,PetscInt v,Mat *Jin,Mat *J[])
191: {
193:   Mat            *Jv;
194:   PetscInt       nedges,e,i,M,N,*rows,*cols;
195:   PetscBool      isSelf;
196:   const PetscInt *edges,*cone;
197:   PetscScalar    *zeros;

200:   /* Get arrary size of Jv */
201:   DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
202:   if (nedges <= 0) SETERRQ2(PETSC_COMM_SELF,1,"%d vertex, nedges %d\n",v,nedges);

204:   /* two Jacobians for each connected edge: J(v,e) and J(v,vc); adding J(v,v), total 2*nedges+1 Jacobians */
205:   PetscCalloc1(2*nedges+1,&Jv);

207:   /* Create dense zero block for this vertex: J[0] = Jacobian(v,v) */
208:   DMNetworkGetNumVariables(dm,v,&M);
209:   if (M !=2) SETERRQ1(PETSC_COMM_SELF,1,"M != 2",M);
210:   PetscMalloc3(M,&rows,M,&cols,M*M,&zeros);
211:   PetscArrayzero(zeros,M*M);
212:   for (i=0; i<M; i++) rows[i] = i;

214:   for (e=0; e<nedges; e++) {
215:     /* create Jv[2*e+1] = Jacobian(v,e), e: supporting edge */
216:     DMNetworkGetConnectedVertices(dm,edges[e],&cone);
217:     isSelf = (v == cone[0]) ? PETSC_TRUE:PETSC_FALSE;

219:     if (Jin) {
220:       if (isSelf) {
221:         Jv[2*e+1] = Jin[0];
222:       } else {
223:         Jv[2*e+1] = Jin[1];
224:       }
225:       Jv[2*e+2] = Jin[2];
226:       PetscObjectReference((PetscObject)(Jv[2*e+1]));
227:       PetscObjectReference((PetscObject)(Jv[2*e+2]));
228:     } else {
229:       /* create J(v,e) */
230:       MatCreate(PETSC_COMM_SELF,&Jv[2*e+1]);
231:       DMNetworkGetNumVariables(dm,edges[e],&N);
232:       MatSetSizes(Jv[2*e+1],PETSC_DECIDE,PETSC_DECIDE,M,N);
233:       MatSetFromOptions(Jv[2*e+1]);
234:       MatSetOption(Jv[2*e+1],MAT_STRUCTURE_ONLY,PETSC_TRUE);
235:       MatSeqAIJSetPreallocation(Jv[2*e+1],2,NULL);
236:       if (N) {
237:         if (isSelf) { /* coupling at upstream */
238:           for (i=0; i<2; i++) cols[i] = i;
239:         } else { /* coupling at downstream */
240:           cols[0] = N-2; cols[1] = N-1;
241:         }
242:         MatSetValues(Jv[2*e+1],2,rows,2,cols,zeros,INSERT_VALUES);
243:       }
244:       MatAssemblyBegin(Jv[2*e+1],MAT_FINAL_ASSEMBLY);
245:       MatAssemblyEnd(Jv[2*e+1],MAT_FINAL_ASSEMBLY);

247:       /* create Jv[2*e+2] = Jacobian(v,vc), vc: connected vertex.
248:        In WashNetwork, v and vc are not connected, thus Jacobian(v,vc) is empty */
249:       MatCreate(PETSC_COMM_SELF,&Jv[2*e+2]);
250:       MatSetSizes(Jv[2*e+2],PETSC_DECIDE,PETSC_DECIDE,M,M); /* empty matrix, sizes can be arbitrary */
251:       MatSetFromOptions(Jv[2*e+2]);
252:       MatSetOption(Jv[2*e+2],MAT_STRUCTURE_ONLY,PETSC_TRUE);
253:       MatSeqAIJSetPreallocation(Jv[2*e+2],1,NULL);
254:       MatAssemblyBegin(Jv[2*e+2],MAT_FINAL_ASSEMBLY);
255:       MatAssemblyEnd(Jv[2*e+2],MAT_FINAL_ASSEMBLY);
256:     }
257:   }
258:   PetscFree3(rows,cols,zeros);

260:   *J = Jv;
261:   return(0);
262: }

264: PetscErrorCode JunctionDestroyJacobian(DM dm,PetscInt v,Junction junc)
265: {
267:   Mat            *Jv=junc->jacobian;
268:   const PetscInt *edges;
269:   PetscInt       nedges,e;

272:   if (!Jv) return(0);

274:   DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
275:   for (e=0; e<nedges; e++) {
276:     MatDestroy(&Jv[2*e+1]);
277:     MatDestroy(&Jv[2*e+2]);
278:   }
279:   PetscFree(Jv);
280:   return(0);
281: }