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