Actual source code: pipes.c

      1: static char help[] = "This example demonstrates DMNetwork. It is used for testing parallel generation of dmnetwork, then redistribute. \n\\n";
      2: /*
      3:   Example: mpiexec -n <np> ./pipes -ts_max_steps 10
      4: */

      6: #include "wash.h"

      8: /*
      9:   WashNetworkDistribute - proc[0] distributes sequential wash object
     10:    Input Parameters:
     11: .  comm - MPI communicator
     12: .  wash - wash context with all network data in proc[0]

     14:    Output Parameter:
     15: .  wash - wash context with nedge, nvertex and edgelist distributed

     17:    Note: The routine is used for testing parallel generation of dmnetwork, then redistribute.
     18: */
     19: PetscErrorCode WashNetworkDistribute(MPI_Comm comm, Wash wash)
     20: {
     21:   PetscMPIInt rank, size, tag = 0;
     22:   PetscInt    i, e, v, numEdges, numVertices, nedges, *eowners = NULL, estart, eend, *vtype = NULL, nvertices;
     23:   PetscInt   *edgelist = wash->edgelist, *nvtx = NULL, *vtxDone = NULL;

     25:   PetscFunctionBegin;
     26:   PetscCallMPI(MPI_Comm_size(comm, &size));
     27:   if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);

     29:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
     30:   numEdges    = wash->nedge;
     31:   numVertices = wash->nvertex;

     33:   /* (1) all processes get global and local number of edges */
     34:   PetscCallMPI(MPI_Bcast(&numEdges, 1, MPIU_INT, 0, comm));
     35:   nedges = numEdges / size; /* local nedges */
     36:   if (rank == 0) nedges += numEdges - size * (numEdges / size);
     37:   wash->Nedge = numEdges;
     38:   wash->nedge = nedges;
     39:   /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d] nedges %d, numEdges %d\n",rank,nedges,numEdges)); */

     41:   PetscCall(PetscCalloc3(size + 1, &eowners, size, &nvtx, numVertices, &vtxDone));
     42:   PetscCallMPI(MPI_Allgather(&nedges, 1, MPIU_INT, eowners + 1, 1, MPIU_INT, PETSC_COMM_WORLD));
     43:   eowners[0] = 0;
     44:   for (i = 2; i <= size; i++) eowners[i] += eowners[i - 1];

     46:   estart = eowners[rank];
     47:   eend   = eowners[rank + 1];
     48:   /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d] own lists row %d - %d\n",rank,estart,eend)); */

     50:   /* (2) distribute row block edgelist to all processors */
     51:   if (rank == 0) {
     52:     vtype = wash->vtype;
     53:     for (i = 1; i < size; i++) {
     54:       /* proc[0] sends edgelist to proc[i] */
     55:       PetscCallMPI(MPI_Send(edgelist + 2 * eowners[i], 2 * (eowners[i + 1] - eowners[i]), MPIU_INT, i, tag, comm));

     57:       /* proc[0] sends vtype to proc[i] */
     58:       PetscCallMPI(MPI_Send(vtype + 2 * eowners[i], 2 * (eowners[i + 1] - eowners[i]), MPIU_INT, i, tag, comm));
     59:     }
     60:   } else {
     61:     MPI_Status status;
     62:     PetscCall(PetscMalloc1(2 * (eend - estart), &vtype));
     63:     PetscCall(PetscMalloc1(2 * (eend - estart), &edgelist));

     65:     PetscCallMPI(MPI_Recv(edgelist, 2 * (eend - estart), MPIU_INT, 0, tag, comm, &status));
     66:     PetscCallMPI(MPI_Recv(vtype, 2 * (eend - estart), MPIU_INT, 0, tag, comm, &status));
     67:   }

     69:   wash->edgelist = edgelist;

     71:   /* (3) all processes get global and local number of vertices, without ghost vertices */
     72:   if (rank == 0) {
     73:     for (i = 0; i < size; i++) {
     74:       for (e = eowners[i]; e < eowners[i + 1]; e++) {
     75:         v = edgelist[2 * e];
     76:         if (!vtxDone[v]) {
     77:           nvtx[i]++;
     78:           vtxDone[v] = 1;
     79:         }
     80:         v = edgelist[2 * e + 1];
     81:         if (!vtxDone[v]) {
     82:           nvtx[i]++;
     83:           vtxDone[v] = 1;
     84:         }
     85:       }
     86:     }
     87:   }
     88:   PetscCallMPI(MPI_Bcast(&numVertices, 1, MPIU_INT, 0, PETSC_COMM_WORLD));
     89:   PetscCallMPI(MPI_Scatter(nvtx, 1, MPIU_INT, &nvertices, 1, MPIU_INT, 0, PETSC_COMM_WORLD));
     90:   PetscCall(PetscFree3(eowners, nvtx, vtxDone));

     92:   wash->Nvertex = numVertices;
     93:   wash->nvertex = nvertices;
     94:   wash->vtype   = vtype;
     95:   PetscFunctionReturn(PETSC_SUCCESS);
     96: }

     98: PetscErrorCode WASHIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
     99: {
    100:   Wash               wash = (Wash)ctx;
    101:   DM                 networkdm;
    102:   Vec                localX, localXdot, localF, localXold;
    103:   const PetscInt    *cone;
    104:   PetscInt           vfrom, vto, offsetfrom, offsetto, varoffset;
    105:   PetscInt           v, vStart, vEnd, e, eStart, eEnd;
    106:   PetscInt           nend, type;
    107:   PetscBool          ghost;
    108:   PetscScalar       *farr, *juncf, *pipef;
    109:   PetscReal          dt;
    110:   Pipe               pipe;
    111:   PipeField         *pipex, *pipexdot, *juncx;
    112:   Junction           junction;
    113:   DMDALocalInfo      info;
    114:   const PetscScalar *xarr, *xdotarr, *xoldarr;

    116:   PetscFunctionBegin;
    117:   localX    = wash->localX;
    118:   localXdot = wash->localXdot;

    120:   PetscCall(TSGetSolution(ts, &localXold));
    121:   PetscCall(TSGetDM(ts, &networkdm));
    122:   PetscCall(TSGetTimeStep(ts, &dt));
    123:   PetscCall(DMGetLocalVector(networkdm, &localF));

    125:   /* Set F and localF as zero */
    126:   PetscCall(VecSet(F, 0.0));
    127:   PetscCall(VecSet(localF, 0.0));

    129:   /* Update ghost values of locaX and locaXdot */
    130:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
    131:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

    133:   PetscCall(DMGlobalToLocalBegin(networkdm, Xdot, INSERT_VALUES, localXdot));
    134:   PetscCall(DMGlobalToLocalEnd(networkdm, Xdot, INSERT_VALUES, localXdot));

    136:   PetscCall(VecGetArrayRead(localX, &xarr));
    137:   PetscCall(VecGetArrayRead(localXdot, &xdotarr));
    138:   PetscCall(VecGetArrayRead(localXold, &xoldarr));
    139:   PetscCall(VecGetArray(localF, &farr));

    141:   /* junction->type == JUNCTION:
    142:            juncf[0] = -qJ + sum(qin); juncf[1] = qJ - sum(qout)
    143:        junction->type == RESERVOIR (upper stream):
    144:            juncf[0] = -hJ + H0; juncf[1] = qJ - sum(qout)
    145:        junction->type == VALVE (down stream):
    146:            juncf[0] =  -qJ + sum(qin); juncf[1] = qJ
    147:   */
    148:   /* Vertex/junction initialization */
    149:   PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
    150:   for (v = vStart; v < vEnd; v++) {
    151:     PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghost));
    152:     if (ghost) continue;

    154:     PetscCall(DMNetworkGetComponent(networkdm, v, 0, &type, (void **)&junction, NULL));
    155:     PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &varoffset));
    156:     juncx = (PipeField *)(xarr + varoffset);
    157:     juncf = (PetscScalar *)(farr + varoffset);

    159:     juncf[0] = -juncx[0].q;
    160:     juncf[1] = juncx[0].q;

    162:     if (junction->type == RESERVOIR) { /* upstream reservoir */
    163:       juncf[0] = juncx[0].h - wash->H0;
    164:     }
    165:   }

    167:   /* Edge/pipe */
    168:   PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
    169:   for (e = eStart; e < eEnd; e++) {
    170:     PetscCall(DMNetworkGetComponent(networkdm, e, 0, &type, (void **)&pipe, NULL));
    171:     PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &varoffset));
    172:     pipex    = (PipeField *)(xarr + varoffset);
    173:     pipexdot = (PipeField *)(xdotarr + varoffset);
    174:     pipef    = (PetscScalar *)(farr + varoffset);

    176:     /* Get some data into the pipe structure: note, some of these operations
    177:      * might be redundant. Will it consume too much time? */
    178:     pipe->dt   = dt;
    179:     pipe->xold = (PipeField *)(xoldarr + varoffset);

    181:     /* Evaluate F over this edge/pipe: pipef[1], ...,pipef[2*nend] */
    182:     PetscCall(DMDAGetLocalInfo(pipe->da, &info));
    183:     PetscCall(PipeIFunctionLocal_Lax(&info, t, pipex, pipexdot, pipef, pipe));

    185:     /* Get boundary values from connected vertices */
    186:     PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
    187:     vfrom = cone[0]; /* local ordering */
    188:     vto   = cone[1];
    189:     PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom));
    190:     PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto));

    192:     /* Evaluate upstream boundary */
    193:     PetscCall(DMNetworkGetComponent(networkdm, vfrom, 0, &type, (void **)&junction, NULL));
    194:     PetscCheck(junction->type == JUNCTION || junction->type == RESERVOIR, PETSC_COMM_SELF, PETSC_ERR_SUP, "junction type is not supported");
    195:     juncx = (PipeField *)(xarr + offsetfrom);
    196:     juncf = (PetscScalar *)(farr + offsetfrom);

    198:     pipef[0] = pipex[0].h - juncx[0].h;
    199:     juncf[1] -= pipex[0].q;

    201:     /* Evaluate downstream boundary */
    202:     PetscCall(DMNetworkGetComponent(networkdm, vto, 0, &type, (void **)&junction, NULL));
    203:     PetscCheck(junction->type == JUNCTION || junction->type == VALVE, PETSC_COMM_SELF, PETSC_ERR_SUP, "junction type is not supported");
    204:     juncx = (PipeField *)(xarr + offsetto);
    205:     juncf = (PetscScalar *)(farr + offsetto);
    206:     nend  = pipe->nnodes - 1;

    208:     pipef[2 * nend + 1] = pipex[nend].h - juncx[0].h;
    209:     juncf[0] += pipex[nend].q;
    210:   }

    212:   PetscCall(VecRestoreArrayRead(localX, &xarr));
    213:   PetscCall(VecRestoreArrayRead(localXdot, &xdotarr));
    214:   PetscCall(VecRestoreArray(localF, &farr));

    216:   PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
    217:   PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
    218:   PetscCall(DMRestoreLocalVector(networkdm, &localF));
    219:   /*
    220:    PetscCall(PetscPrintf(PETSC_COMM_WORLD("F:\n"));
    221:    PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD));
    222:    */
    223:   PetscFunctionReturn(PETSC_SUCCESS);
    224: }

    226: PetscErrorCode WASHSetInitialSolution(DM networkdm, Vec X, Wash wash)
    227: {
    228:   PetscInt           k, nx, vkey, vfrom, vto, offsetfrom, offsetto;
    229:   PetscInt           type, varoffset;
    230:   PetscInt           e, eStart, eEnd;
    231:   Vec                localX;
    232:   PetscScalar       *xarr;
    233:   Pipe               pipe;
    234:   Junction           junction;
    235:   const PetscInt    *cone;
    236:   const PetscScalar *xarray;

    238:   PetscFunctionBegin;
    239:   PetscCall(VecSet(X, 0.0));
    240:   PetscCall(DMGetLocalVector(networkdm, &localX));
    241:   PetscCall(VecGetArray(localX, &xarr));

    243:   /* Edge */
    244:   PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
    245:   for (e = eStart; e < eEnd; e++) {
    246:     PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &varoffset));
    247:     PetscCall(DMNetworkGetComponent(networkdm, e, 0, &type, (void **)&pipe, NULL));

    249:     /* set initial values for this pipe */
    250:     PetscCall(PipeComputeSteadyState(pipe, wash->Q0, wash->H0));
    251:     PetscCall(VecGetSize(pipe->x, &nx));

    253:     PetscCall(VecGetArrayRead(pipe->x, &xarray));
    254:     /* copy pipe->x to xarray */
    255:     for (k = 0; k < nx; k++) (xarr + varoffset)[k] = xarray[k];

    257:     /* set boundary values into vfrom and vto */
    258:     PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
    259:     vfrom = cone[0]; /* local ordering */
    260:     vto   = cone[1];
    261:     PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom));
    262:     PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto));

    264:     /* if vform is a head vertex: */
    265:     PetscCall(DMNetworkGetComponent(networkdm, vfrom, 0, &vkey, (void **)&junction, NULL));
    266:     if (junction->type == RESERVOIR) (xarr + offsetfrom)[1] = wash->H0; /* 1st H */

    268:     /* if vto is an end vertex: */
    269:     PetscCall(DMNetworkGetComponent(networkdm, vto, 0, &vkey, (void **)&junction, NULL));
    270:     if (junction->type == VALVE) (xarr + offsetto)[0] = wash->QL; /* last Q */
    271:     PetscCall(VecRestoreArrayRead(pipe->x, &xarray));
    272:   }

    274:   PetscCall(VecRestoreArray(localX, &xarr));
    275:   PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
    276:   PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
    277:   PetscCall(DMRestoreLocalVector(networkdm, &localX));

    279: #if 0
    280:   PetscInt N;
    281:   PetscCall(VecGetSize(X,&N));
    282:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"initial solution %d:\n",N));
    283:   PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD));
    284: #endif
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }

    288: PetscErrorCode TSDMNetworkMonitor(TS ts, PetscInt step, PetscReal t, Vec x, void *context)
    289: {
    290:   DMNetworkMonitor monitor;

    292:   PetscFunctionBegin;
    293:   monitor = (DMNetworkMonitor)context;
    294:   PetscCall(DMNetworkMonitorView(monitor, x));
    295:   PetscFunctionReturn(PETSC_SUCCESS);
    296: }

    298: PetscErrorCode PipesView(DM networkdm, PetscInt KeyPipe, Vec X)
    299: {
    300:   PetscInt   i, numkeys = 1, *blocksize, *numselectedvariable, **selectedvariables, n;
    301:   IS         isfrom_q, isfrom_h, isfrom;
    302:   Vec        Xto;
    303:   VecScatter ctx;
    304:   MPI_Comm   comm;

    306:   PetscFunctionBegin;
    307:   PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));

    309:   /* 1. Create isfrom_q for q-variable of pipes */
    310:   PetscCall(PetscMalloc3(numkeys, &blocksize, numkeys, &numselectedvariable, numkeys, &selectedvariables));
    311:   for (i = 0; i < numkeys; i++) {
    312:     blocksize[i]           = 2;
    313:     numselectedvariable[i] = 1;
    314:     PetscCall(PetscMalloc1(numselectedvariable[i], &selectedvariables[i]));
    315:     selectedvariables[i][0] = 0; /* q-variable */
    316:   }
    317:   PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyPipe, blocksize, numselectedvariable, selectedvariables, &isfrom_q));

    319:   /* 2. Create Xto and isto */
    320:   PetscCall(ISGetLocalSize(isfrom_q, &n));
    321:   PetscCall(VecCreate(comm, &Xto));
    322:   PetscCall(VecSetSizes(Xto, n, PETSC_DECIDE));
    323:   PetscCall(VecSetFromOptions(Xto));
    324:   PetscCall(VecSet(Xto, 0.0));

    326:   /* 3. Create scatter */
    327:   PetscCall(VecScatterCreate(X, isfrom_q, Xto, NULL, &ctx));

    329:   /* 4. Scatter to Xq */
    330:   PetscCall(VecScatterBegin(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
    331:   PetscCall(VecScatterEnd(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
    332:   PetscCall(VecScatterDestroy(&ctx));
    333:   PetscCall(ISDestroy(&isfrom_q));

    335:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Xq:\n"));
    336:   PetscCall(VecView(Xto, PETSC_VIEWER_STDOUT_WORLD));

    338:   /* 5. Create isfrom_h for h-variable of pipes; Create scatter; Scatter to Xh */
    339:   for (i = 0; i < numkeys; i++) selectedvariables[i][0] = 1; /* h-variable */
    340:   PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyPipe, blocksize, numselectedvariable, selectedvariables, &isfrom_h));

    342:   PetscCall(VecScatterCreate(X, isfrom_h, Xto, NULL, &ctx));
    343:   PetscCall(VecScatterBegin(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
    344:   PetscCall(VecScatterEnd(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
    345:   PetscCall(VecScatterDestroy(&ctx));
    346:   PetscCall(ISDestroy(&isfrom_h));

    348:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Xh:\n"));
    349:   PetscCall(VecView(Xto, PETSC_VIEWER_STDOUT_WORLD));
    350:   PetscCall(VecDestroy(&Xto));

    352:   /* 6. Create isfrom for all pipe variables; Create scatter; Scatter to Xpipes */
    353:   for (i = 0; i < numkeys; i++) blocksize[i] = -1; /* select all the variables of the i-th component */
    354:   PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyPipe, blocksize, NULL, NULL, &isfrom));
    355:   PetscCall(ISDestroy(&isfrom));
    356:   PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyPipe, NULL, NULL, NULL, &isfrom));

    358:   PetscCall(ISGetLocalSize(isfrom, &n));
    359:   PetscCall(VecCreate(comm, &Xto));
    360:   PetscCall(VecSetSizes(Xto, n, PETSC_DECIDE));
    361:   PetscCall(VecSetFromOptions(Xto));
    362:   PetscCall(VecSet(Xto, 0.0));

    364:   PetscCall(VecScatterCreate(X, isfrom, Xto, NULL, &ctx));
    365:   PetscCall(VecScatterBegin(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
    366:   PetscCall(VecScatterEnd(ctx, X, Xto, INSERT_VALUES, SCATTER_FORWARD));
    367:   PetscCall(VecScatterDestroy(&ctx));
    368:   PetscCall(ISDestroy(&isfrom));

    370:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Xpipes:\n"));
    371:   PetscCall(VecView(Xto, PETSC_VIEWER_STDOUT_WORLD));

    373:   /* 7. Free spaces */
    374:   for (i = 0; i < numkeys; i++) PetscCall(PetscFree(selectedvariables[i]));
    375:   PetscCall(PetscFree3(blocksize, numselectedvariable, selectedvariables));
    376:   PetscCall(VecDestroy(&Xto));
    377:   PetscFunctionReturn(PETSC_SUCCESS);
    378: }

    380: PetscErrorCode ISJunctionsView(DM networkdm, PetscInt KeyJunc)
    381: {
    382:   PetscInt    numkeys = 1;
    383:   IS          isfrom;
    384:   MPI_Comm    comm;
    385:   PetscMPIInt rank;

    387:   PetscFunctionBegin;
    388:   PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
    389:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    391:   /* Create a global isfrom for all junction variables */
    392:   PetscCall(DMNetworkCreateIS(networkdm, numkeys, &KeyJunc, NULL, NULL, NULL, &isfrom));
    393:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ISJunctions:\n"));
    394:   PetscCall(ISView(isfrom, PETSC_VIEWER_STDOUT_WORLD));
    395:   PetscCall(ISDestroy(&isfrom));

    397:   /* Create a local isfrom for all junction variables */
    398:   PetscCall(DMNetworkCreateLocalIS(networkdm, numkeys, &KeyJunc, NULL, NULL, NULL, &isfrom));
    399:   if (rank == 0) {
    400:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] ISLocalJunctions:\n", rank));
    401:     PetscCall(ISView(isfrom, PETSC_VIEWER_STDOUT_SELF));
    402:   }
    403:   PetscCall(ISDestroy(&isfrom));
    404:   PetscFunctionReturn(PETSC_SUCCESS);
    405: }

    407: PetscErrorCode WashNetworkCleanUp(Wash wash)
    408: {
    409:   PetscMPIInt rank;

    411:   PetscFunctionBegin;
    412:   PetscCallMPI(MPI_Comm_rank(wash->comm, &rank));
    413:   PetscCall(PetscFree(wash->edgelist));
    414:   PetscCall(PetscFree(wash->vtype));
    415:   if (rank == 0) PetscCall(PetscFree2(wash->junction, wash->pipe));
    416:   PetscFunctionReturn(PETSC_SUCCESS);
    417: }

    419: PetscErrorCode WashNetworkCreate(MPI_Comm comm, PetscInt pipesCase, Wash *wash_ptr)
    420: {
    421:   PetscInt    npipes;
    422:   PetscMPIInt rank;
    423:   Wash        wash = NULL;
    424:   PetscInt    i, numVertices, numEdges, *vtype;
    425:   PetscInt   *edgelist;
    426:   Junction    junctions = NULL;
    427:   Pipe        pipes     = NULL;
    428:   PetscBool   washdist  = PETSC_TRUE;

    430:   PetscFunctionBegin;
    431:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    433:   PetscCall(PetscCalloc1(1, &wash));
    434:   wash->comm       = comm;
    435:   *wash_ptr        = wash;
    436:   wash->Q0         = 0.477432; /* RESERVOIR */
    437:   wash->H0         = 150.0;
    438:   wash->HL         = 143.488; /* VALVE */
    439:   wash->QL         = 0.0;
    440:   wash->nnodes_loc = 0;

    442:   numVertices = 0;
    443:   numEdges    = 0;
    444:   edgelist    = NULL;

    446:   /* proc[0] creates a sequential wash and edgelist */
    447:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setup pipesCase %" PetscInt_FMT "\n", pipesCase));

    449:   /* Set global number of pipes, edges, and junctions */
    450:   /*-------------------------------------------------*/
    451:   switch (pipesCase) {
    452:   case 0:
    453:     /* pipeCase 0: */
    454:     /* =================================================
    455:     (RESERVOIR) v0 --E0--> v1--E1--> v2 --E2-->v3 (VALVE)
    456:     ====================================================  */
    457:     npipes = 3;
    458:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-npipes", &npipes, NULL));
    459:     wash->nedge   = npipes;
    460:     wash->nvertex = npipes + 1;

    462:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
    463:     numVertices = 0;
    464:     numEdges    = 0;
    465:     edgelist    = NULL;
    466:     if (rank == 0) {
    467:       numVertices = wash->nvertex;
    468:       numEdges    = wash->nedge;

    470:       PetscCall(PetscCalloc1(2 * numEdges, &edgelist));
    471:       for (i = 0; i < numEdges; i++) {
    472:         edgelist[2 * i]     = i;
    473:         edgelist[2 * i + 1] = i + 1;
    474:       }

    476:       /* Add network components */
    477:       /*------------------------*/
    478:       PetscCall(PetscCalloc2(numVertices, &junctions, numEdges, &pipes));

    480:       /* vertex */
    481:       for (i = 0; i < numVertices; i++) {
    482:         junctions[i].id   = i;
    483:         junctions[i].type = JUNCTION;
    484:       }

    486:       junctions[0].type               = RESERVOIR;
    487:       junctions[numVertices - 1].type = VALVE;
    488:     }
    489:     break;
    490:   case 1:
    491:     /* pipeCase 1: */
    492:     /* ==========================
    493:                 v2 (VALVE)
    494:                 ^
    495:                 |
    496:                E2
    497:                 |
    498:     v0 --E0--> v3--E1--> v1
    499:   (RESERVOIR)            (RESERVOIR)
    500:     =============================  */
    501:     npipes        = 3;
    502:     wash->nedge   = npipes;
    503:     wash->nvertex = npipes + 1;

    505:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
    506:     if (rank == 0) {
    507:       numVertices = wash->nvertex;
    508:       numEdges    = wash->nedge;

    510:       PetscCall(PetscCalloc1(2 * numEdges, &edgelist));
    511:       edgelist[0] = 0;
    512:       edgelist[1] = 3; /* edge[0] */
    513:       edgelist[2] = 3;
    514:       edgelist[3] = 1; /* edge[1] */
    515:       edgelist[4] = 3;
    516:       edgelist[5] = 2; /* edge[2] */

    518:       /* Add network components */
    519:       /*------------------------*/
    520:       PetscCall(PetscCalloc2(numVertices, &junctions, numEdges, &pipes));
    521:       /* vertex */
    522:       for (i = 0; i < numVertices; i++) {
    523:         junctions[i].id   = i;
    524:         junctions[i].type = JUNCTION;
    525:       }

    527:       junctions[0].type = RESERVOIR;
    528:       junctions[1].type = VALVE;
    529:       junctions[2].type = VALVE;
    530:     }
    531:     break;
    532:   case 2:
    533:     /* pipeCase 2: */
    534:     /* ==========================
    535:     (RESERVOIR)  v2--> E2
    536:                        |
    537:             v0 --E0--> v3--E1--> v1
    538:     (RESERVOIR)               (VALVE)
    539:     =============================  */

    541:     /* Set application parameters -- to be used in function evaluations */
    542:     npipes        = 3;
    543:     wash->nedge   = npipes;
    544:     wash->nvertex = npipes + 1;

    546:     /* Set local edges and vertices -- proc[0] sets entire network, then distributes */
    547:     if (rank == 0) {
    548:       numVertices = wash->nvertex;
    549:       numEdges    = wash->nedge;

    551:       PetscCall(PetscCalloc1(2 * numEdges, &edgelist));
    552:       edgelist[0] = 0;
    553:       edgelist[1] = 3; /* edge[0] */
    554:       edgelist[2] = 3;
    555:       edgelist[3] = 1; /* edge[1] */
    556:       edgelist[4] = 2;
    557:       edgelist[5] = 3; /* edge[2] */

    559:       /* Add network components */
    560:       /*------------------------*/
    561:       PetscCall(PetscCalloc2(numVertices, &junctions, numEdges, &pipes));
    562:       /* vertex */
    563:       for (i = 0; i < numVertices; i++) {
    564:         junctions[i].id   = i;
    565:         junctions[i].type = JUNCTION;
    566:       }

    568:       junctions[0].type = RESERVOIR;
    569:       junctions[1].type = VALVE;
    570:       junctions[2].type = RESERVOIR;
    571:     }
    572:     break;
    573:   default:
    574:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "not done yet");
    575:   }

    577:   /* set edge global id */
    578:   for (i = 0; i < numEdges; i++) pipes[i].id = i;

    580:   if (rank == 0) { /* set vtype for proc[0] */
    581:     PetscInt v;
    582:     PetscCall(PetscMalloc1(2 * numEdges, &vtype));
    583:     for (i = 0; i < 2 * numEdges; i++) {
    584:       v        = edgelist[i];
    585:       vtype[i] = junctions[v].type;
    586:     }
    587:     wash->vtype = vtype;
    588:   }

    590:   *wash_ptr      = wash;
    591:   wash->nedge    = numEdges;
    592:   wash->nvertex  = numVertices;
    593:   wash->edgelist = edgelist;
    594:   wash->junction = junctions;
    595:   wash->pipe     = pipes;

    597:   /* Distribute edgelist to other processors */
    598:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-wash_distribute", &washdist, NULL));
    599:   if (washdist) {
    600:     /*
    601:      PetscCall(PetscPrintf(PETSC_COMM_WORLD," Distribute sequential wash ...\n"));
    602:      */
    603:     PetscCall(WashNetworkDistribute(comm, wash));
    604:   }
    605:   PetscFunctionReturn(PETSC_SUCCESS);
    606: }

    608: /* ------------------------------------------------------- */
    609: int main(int argc, char **argv)
    610: {
    611:   Wash              wash;
    612:   Junction          junctions, junction;
    613:   Pipe              pipe, pipes;
    614:   PetscInt          KeyPipe, KeyJunction, *edgelist = NULL, *vtype = NULL;
    615:   PetscInt          i, e, v, eStart, eEnd, vStart, vEnd, key, vkey, type;
    616:   PetscInt          steps = 1, nedges, nnodes = 6;
    617:   const PetscInt   *cone;
    618:   DM                networkdm;
    619:   PetscMPIInt       size, rank;
    620:   PetscReal         ftime;
    621:   Vec               X;
    622:   TS                ts;
    623:   TSConvergedReason reason;
    624:   PetscBool         viewpipes, viewjuncs, monipipes = PETSC_FALSE, userJac = PETSC_TRUE, viewdm = PETSC_FALSE, viewX = PETSC_FALSE;
    625:   PetscInt          pipesCase = 0;
    626:   DMNetworkMonitor  monitor;
    627:   MPI_Comm          comm;

    629:   PetscFunctionBeginUser;
    630:   PetscCall(PetscInitialize(&argc, &argv, "pOption", help));

    632:   /* Read runtime options */
    633:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-case", &pipesCase, NULL));
    634:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-user_Jac", &userJac, NULL));
    635:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-pipe_monitor", &monipipes, NULL));
    636:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewdm", &viewdm, NULL));
    637:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewX", &viewX, NULL));
    638:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-npipenodes", &nnodes, NULL));

    640:   /* Create networkdm */
    641:   /*------------------*/
    642:   PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
    643:   PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
    644:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    645:   PetscCallMPI(MPI_Comm_size(comm, &size));

    647:   if (size == 1 && monipipes) PetscCall(DMNetworkMonitorCreate(networkdm, &monitor));

    649:   /* Register the components in the network */
    650:   PetscCall(DMNetworkRegisterComponent(networkdm, "junctionstruct", sizeof(struct _p_Junction), &KeyJunction));
    651:   PetscCall(DMNetworkRegisterComponent(networkdm, "pipestruct", sizeof(struct _p_Pipe), &KeyPipe));

    653:   /* Create a distributed wash network (user-specific) */
    654:   PetscCall(WashNetworkCreate(comm, pipesCase, &wash));
    655:   nedges    = wash->nedge;
    656:   edgelist  = wash->edgelist;
    657:   vtype     = wash->vtype;
    658:   junctions = wash->junction;
    659:   pipes     = wash->pipe;

    661:   /* Set up the network layout */
    662:   PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1));
    663:   PetscCall(DMNetworkAddSubnetwork(networkdm, NULL, nedges, edgelist, NULL));

    665:   PetscCall(DMNetworkLayoutSetUp(networkdm));

    667:   PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
    668:   PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
    669:   /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d] eStart/End: %d - %d; vStart/End: %d - %d\n",rank,eStart,eEnd,vStart,vEnd)); */

    671:   if (rank) { /* junctions[] and pipes[] for proc[0] are allocated in WashNetworkCreate() */
    672:     /* vEnd - vStart = nvertices + number of ghost vertices! */
    673:     PetscCall(PetscCalloc2(vEnd - vStart, &junctions, nedges, &pipes));
    674:   }

    676:   /* Add Pipe component and number of variables to all local edges */
    677:   for (e = eStart; e < eEnd; e++) {
    678:     pipes[e - eStart].nnodes = nnodes;
    679:     PetscCall(DMNetworkAddComponent(networkdm, e, KeyPipe, &pipes[e - eStart], 2 * pipes[e - eStart].nnodes));

    681:     if (size == 1 && monipipes) { /* Add monitor -- show Q_{pipes[e-eStart].id}? */
    682:       pipes[e - eStart].length = 600.0;
    683:       PetscCall(DMNetworkMonitorAdd(monitor, "Pipe Q", e, pipes[e - eStart].nnodes, 0, 2, 0.0, pipes[e - eStart].length, -0.8, 0.8, PETSC_TRUE));
    684:       PetscCall(DMNetworkMonitorAdd(monitor, "Pipe H", e, pipes[e - eStart].nnodes, 1, 2, 0.0, pipes[e - eStart].length, -400.0, 800.0, PETSC_TRUE));
    685:     }
    686:   }

    688:   /* Add Junction component and number of variables to all local vertices */
    689:   for (v = vStart; v < vEnd; v++) PetscCall(DMNetworkAddComponent(networkdm, v, KeyJunction, &junctions[v - vStart], 2));

    691:   if (size > 1) { /* must be called before DMSetUp()???. Other partitioners do not work yet??? -- cause crash in proc[0]! */
    692:     DM               plexdm;
    693:     PetscPartitioner part;
    694:     PetscCall(DMNetworkGetPlex(networkdm, &plexdm));
    695:     PetscCall(DMPlexGetPartitioner(plexdm, &part));
    696:     PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSIMPLE));
    697:     PetscCall(PetscOptionsSetValue(NULL, "-dm_plex_csr_alg", "mat")); /* for parmetis */
    698:   }

    700:   /* Set up DM for use */
    701:   PetscCall(DMSetUp(networkdm));
    702:   if (viewdm) {
    703:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nOriginal networkdm, DMView:\n"));
    704:     PetscCall(DMView(networkdm, PETSC_VIEWER_STDOUT_WORLD));
    705:   }

    707:   /* Set user physical parameters to the components */
    708:   for (e = eStart; e < eEnd; e++) {
    709:     PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
    710:     /* vfrom */
    711:     PetscCall(DMNetworkGetComponent(networkdm, cone[0], 0, &vkey, (void **)&junction, NULL));
    712:     junction->type = (VertexType)vtype[2 * e];

    714:     /* vto */
    715:     PetscCall(DMNetworkGetComponent(networkdm, cone[1], 0, &vkey, (void **)&junction, NULL));
    716:     junction->type = (VertexType)vtype[2 * e + 1];
    717:   }

    719:   PetscCall(WashNetworkCleanUp(wash));

    721:   /* Network partitioning and distribution of data */
    722:   PetscCall(DMNetworkDistribute(&networkdm, 0));
    723:   if (viewdm) {
    724:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nAfter DMNetworkDistribute, DMView:\n"));
    725:     PetscCall(DMView(networkdm, PETSC_VIEWER_STDOUT_WORLD));
    726:   }

    728:   /* create vectors */
    729:   PetscCall(DMCreateGlobalVector(networkdm, &X));
    730:   PetscCall(DMCreateLocalVector(networkdm, &wash->localX));
    731:   PetscCall(DMCreateLocalVector(networkdm, &wash->localXdot));

    733:   /* PipeSetUp -- each process only sets its own pipes */
    734:   /*---------------------------------------------------*/
    735:   PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));

    737:   userJac = PETSC_TRUE;
    738:   PetscCall(DMNetworkHasJacobian(networkdm, userJac, userJac));
    739:   PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
    740:   for (e = eStart; e < eEnd; e++) { /* each edge has only one component, pipe */
    741:     PetscCall(DMNetworkGetComponent(networkdm, e, 0, &type, (void **)&pipe, NULL));

    743:     wash->nnodes_loc += pipe->nnodes;        /* local total number of nodes, will be used by PipesView() */
    744:     PetscCall(PipeSetParameters(pipe, 600.0, /* length   */
    745:                                 0.5,         /* diameter */
    746:                                 1200.0,      /* a        */
    747:                                 0.018));     /* friction */
    748:     PetscCall(PipeSetUp(pipe));

    750:     if (userJac) {
    751:       /* Create Jacobian matrix nonzero structures for a Pipe */
    752:       Mat *J;
    753:       PetscCall(PipeCreateJacobian(pipe, NULL, &J));
    754:       PetscCall(DMNetworkEdgeSetMatrix(networkdm, e, J));
    755:     }
    756:   }

    758:   if (userJac) {
    759:     PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
    760:     for (v = vStart; v < vEnd; v++) {
    761:       Mat *J;
    762:       PetscCall(JunctionCreateJacobian(networkdm, v, NULL, &J));
    763:       PetscCall(DMNetworkVertexSetMatrix(networkdm, v, J));

    765:       PetscCall(DMNetworkGetComponent(networkdm, v, 0, &vkey, (void **)&junction, NULL));
    766:       junction->jacobian = J;
    767:     }
    768:   }

    770:   /* Setup solver                                           */
    771:   /*--------------------------------------------------------*/
    772:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));

    774:   PetscCall(TSSetDM(ts, (DM)networkdm));
    775:   PetscCall(TSSetIFunction(ts, NULL, WASHIFunction, wash));

    777:   PetscCall(TSSetMaxSteps(ts, steps));
    778:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    779:   PetscCall(TSSetTimeStep(ts, 0.1));
    780:   PetscCall(TSSetType(ts, TSBEULER));
    781:   if (size == 1 && monipipes) PetscCall(TSMonitorSet(ts, TSDMNetworkMonitor, monitor, NULL));
    782:   PetscCall(TSSetFromOptions(ts));

    784:   PetscCall(WASHSetInitialSolution(networkdm, X, wash));

    786:   PetscCall(TSSolve(ts, X));

    788:   PetscCall(TSGetSolveTime(ts, &ftime));
    789:   PetscCall(TSGetStepNumber(ts, &steps));
    790:   PetscCall(TSGetConvergedReason(ts, &reason));
    791:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));
    792:   if (viewX) PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));

    794:   viewpipes = PETSC_FALSE;
    795:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-Jac_view", &viewpipes, NULL));
    796:   if (viewpipes) {
    797:     SNES snes;
    798:     Mat  Jac;
    799:     PetscCall(TSGetSNES(ts, &snes));
    800:     PetscCall(SNESGetJacobian(snes, &Jac, NULL, NULL, NULL));
    801:     PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
    802:   }

    804:   /* View solutions */
    805:   /* -------------- */
    806:   viewpipes = PETSC_FALSE;
    807:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-pipe_view", &viewpipes, NULL));
    808:   if (viewpipes) PetscCall(PipesView(networkdm, KeyPipe, X));

    810:   /* Test IS */
    811:   viewjuncs = PETSC_FALSE;
    812:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-isjunc_view", &viewjuncs, NULL));
    813:   if (viewjuncs) PetscCall(ISJunctionsView(networkdm, KeyJunction));

    815:   /* Free spaces */
    816:   /* ----------- */
    817:   PetscCall(TSDestroy(&ts));
    818:   PetscCall(VecDestroy(&X));
    819:   PetscCall(VecDestroy(&wash->localX));
    820:   PetscCall(VecDestroy(&wash->localXdot));

    822:   /* Destroy objects from each pipe that are created in PipeSetUp() */
    823:   PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
    824:   for (i = eStart; i < eEnd; i++) {
    825:     PetscCall(DMNetworkGetComponent(networkdm, i, 0, &key, (void **)&pipe, NULL));
    826:     PetscCall(PipeDestroy(&pipe));
    827:   }
    828:   if (userJac) {
    829:     for (v = vStart; v < vEnd; v++) {
    830:       PetscCall(DMNetworkGetComponent(networkdm, v, 0, &vkey, (void **)&junction, NULL));
    831:       PetscCall(JunctionDestroyJacobian(networkdm, v, junction));
    832:     }
    833:   }

    835:   if (size == 1 && monipipes) PetscCall(DMNetworkMonitorDestroy(&monitor));
    836:   PetscCall(DMDestroy(&networkdm));
    837:   PetscCall(PetscFree(wash));

    839:   if (rank) PetscCall(PetscFree2(junctions, pipes));
    840:   PetscCall(PetscFinalize());
    841:   return 0;
    842: }

    844: /*TEST

    846:    build:
    847:      depends: pipeInterface.c pipeImpls.c
    848:      requires: mumps

    850:    test:
    851:       args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -options_left no -viewX
    852:       localrunfiles: pOption
    853:       output_file: output/pipes_1.out

    855:    test:
    856:       suffix: 2
    857:       nsize: 2
    858:       args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
    859:       localrunfiles: pOption
    860:       output_file: output/pipes_2.out

    862:    test:
    863:       suffix: 3
    864:       nsize: 2
    865:       args: -ts_monitor -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
    866:       localrunfiles: pOption
    867:       output_file: output/pipes_3.out

    869:    test:
    870:       suffix: 4
    871:       args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -options_left no -viewX
    872:       localrunfiles: pOption
    873:       output_file: output/pipes_4.out

    875:    test:
    876:       suffix: 5
    877:       nsize: 3
    878:       args: -ts_monitor -case 2 -ts_max_steps 10 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -viewX
    879:       localrunfiles: pOption
    880:       output_file: output/pipes_5.out

    882:    test:
    883:       suffix: 6
    884:       nsize: 2
    885:       args: -ts_monitor -case 1 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX
    886:       localrunfiles: pOption
    887:       output_file: output/pipes_6.out

    889:    test:
    890:       suffix: 7
    891:       nsize: 2
    892:       args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -viewX
    893:       localrunfiles: pOption
    894:       output_file: output/pipes_7.out

    896:    test:
    897:       suffix: 8
    898:       nsize: 2
    899:       requires: parmetis
    900:       args: -ts_monitor -case 2 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type parmetis -options_left no -wash_distribute 1
    901:       localrunfiles: pOption
    902:       output_file: output/pipes_8.out

    904:    test:
    905:       suffix: 9
    906:       nsize: 2
    907:       args: -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -pipe_view
    908:       localrunfiles: pOption
    909:       output_file: output/pipes_9.out

    911:    test:
    912:       suffix: 10
    913:       nsize: 2
    914:       args: -case 0 -ts_max_steps 1 -pc_factor_mat_solver_type mumps -petscpartitioner_type simple -options_left no -wash_distribute 0 -isjunc_view
    915:       localrunfiles: pOption
    916:       output_file: output/pipes_10.out

    918: TEST*/