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:   PetscFunctionBegin;
     18:   PetscCall(PetscNew(pipe));
     19:   PetscFunctionReturn(PETSC_SUCCESS);
     20: }

     22: /*
     23:    PipeDestroy - Destroy Pipe object.

     25:    Input Parameters:
     26:    pipe - Reference to pipe intended to be destroyed.
     27: */
     28: PetscErrorCode PipeDestroy(Pipe *pipe)
     29: {
     30:   PetscFunctionBegin;
     31:   if (!*pipe) PetscFunctionReturn(PETSC_SUCCESS);

     33:   PetscCall(PipeDestroyJacobian(*pipe));
     34:   PetscCall(VecDestroy(&(*pipe)->x));
     35:   PetscCall(DMDestroy(&(*pipe)->da));
     36:   PetscFunctionReturn(PETSC_SUCCESS);
     37: }

     39: /*
     40:    PipeSetParameters - Set parameters for Pipe context

     42:    Input Parameter:
     43: +  pipe - PIPE object
     44: .  length -
     45: .  nnodes -
     46: .  D -
     47: .  a -
     48: -  fric -
     49: */
     50: PetscErrorCode PipeSetParameters(Pipe pipe, PetscReal length, PetscReal D, PetscReal a, PetscReal fric)
     51: {
     52:   PetscFunctionBegin;
     53:   pipe->length = length;
     54:   pipe->D      = D;
     55:   pipe->a      = a;
     56:   pipe->fric   = fric;
     57:   PetscFunctionReturn(PETSC_SUCCESS);
     58: }

     60: /*
     61:     PipeSetUp - Set up pipe based on set parameters.
     62: */
     63: PetscErrorCode PipeSetUp(Pipe pipe)
     64: {
     65:   DMDALocalInfo info;

     67:   PetscFunctionBegin;
     68:   PetscCall(DMDACreate1d(PETSC_COMM_SELF, DM_BOUNDARY_GHOSTED, pipe->nnodes, 2, 1, NULL, &pipe->da));
     69:   PetscCall(DMSetFromOptions(pipe->da));
     70:   PetscCall(DMSetUp(pipe->da));
     71:   PetscCall(DMDASetFieldName(pipe->da, 0, "Q"));
     72:   PetscCall(DMDASetFieldName(pipe->da, 1, "H"));
     73:   PetscCall(DMDASetUniformCoordinates(pipe->da, 0, pipe->length, 0, 0, 0, 0));
     74:   PetscCall(DMCreateGlobalVector(pipe->da, &pipe->x));

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

     78:   pipe->rad = pipe->D / 2;
     79:   pipe->A   = PETSC_PI * pipe->rad * pipe->rad;
     80:   pipe->R   = pipe->fric / (2 * pipe->D * pipe->A);
     81:   PetscFunctionReturn(PETSC_SUCCESS);
     82: }

     84: /*
     85:     PipeCreateJacobian - Create Jacobian matrix nonzero structures for a Pipe.

     87:     Collective

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

     93:     Output Parameter:
     94: .   J  - array of three empty Jacobian matrices

     96:     Level: beginner
     97: */
     98: PetscErrorCode PipeCreateJacobian(Pipe pipe, Mat *Jin, Mat *J[])
     99: {
    100:   Mat         *Jpipe;
    101:   PetscInt     M, rows[2], cols[2], *nz;
    102:   PetscScalar *aa;

    104:   PetscFunctionBegin;
    105:   if (Jin) {
    106:     *J             = Jin;
    107:     pipe->jacobian = Jin;
    108:     PetscCall(PetscObjectReference((PetscObject)Jin[0]));
    109:     PetscFunctionReturn(PETSC_SUCCESS);
    110:   }
    111:   PetscCall(PetscMalloc1(3, &Jpipe));

    113:   /* Jacobian for this pipe */
    114:   PetscCall(DMSetMatrixStructureOnly(pipe->da, PETSC_TRUE));
    115:   PetscCall(DMCreateMatrix(pipe->da, &Jpipe[0]));
    116:   PetscCall(DMSetMatrixStructureOnly(pipe->da, PETSC_FALSE));

    118:   /* Jacobian for upstream vertex */
    119:   PetscCall(MatGetSize(Jpipe[0], &M, NULL));
    120:   PetscCall(PetscCalloc2(M, &nz, 4, &aa));

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

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

    153:   PetscCall(PetscFree2(nz, aa));

    155:   *J             = Jpipe;
    156:   pipe->jacobian = Jpipe;
    157:   PetscFunctionReturn(PETSC_SUCCESS);
    158: }

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

    165:   PetscFunctionBegin;
    166:   if (Jpipe) {
    167:     for (i = 0; i < 3; i++) PetscCall(MatDestroy(&Jpipe[i]));
    168:   }
    169:   PetscCall(PetscFree(Jpipe));
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: /*
    174:     JunctionCreateJacobian - Create Jacobian matrices for a vertex.

    176:     Collective

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

    183:     Output Parameter:
    184: .   J  - array of Jacobian matrices (see dmnetworkimpl.h)

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

    196:   PetscFunctionBegin;
    197:   /* Get array size of Jv */
    198:   PetscCall(DMNetworkGetSupportingEdges(dm, v, &nedges, &edges));
    199:   PetscCheck(nedges > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "%" PetscInt_FMT " vertex, nedges %" PetscInt_FMT, v, nedges);

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

    204:   /* Create dense zero block for this vertex: J[0] = Jacobian(v,v) */
    205:   PetscCall(DMNetworkGetComponent(dm, v, -1, NULL, NULL, &M));
    206:   PetscCheck(M == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "%" PetscInt_FMT " != 2", M);
    207:   PetscCall(PetscMalloc3(M, &rows, M, &cols, M * M, &zeros));
    208:   PetscCall(PetscArrayzero(zeros, M * M));
    209:   for (i = 0; i < M; i++) rows[i] = i;

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

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

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

    258:   *J = Jv;
    259:   PetscFunctionReturn(PETSC_SUCCESS);
    260: }

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

    268:   PetscFunctionBegin;
    269:   if (!Jv) PetscFunctionReturn(PETSC_SUCCESS);

    271:   PetscCall(DMNetworkGetSupportingEdges(dm, v, &nedges, &edges));
    272:   for (e = 0; e < nedges; e++) {
    273:     PetscCall(MatDestroy(&Jv[2 * e + 1]));
    274:     PetscCall(MatDestroy(&Jv[2 * e + 2]));
    275:   }
    276:   PetscCall(PetscFree(Jv));
    277:   PetscFunctionReturn(PETSC_SUCCESS);
    278: }