Actual source code: pipeImpls.c

  1: #include "pipe.h"

  3: /* Initial Function for PIPE       */
  4: /*-------------------------------- */
  5: /*
  6:      Q(x) = Q0 (constant)
  7:      H(x) = H0 - (R/gA) Q0*|Q0|* x
  8:  */
  9: /* ----------------------------------- */
 10: PetscErrorCode PipeComputeSteadyState(Pipe pipe, PetscScalar Q0, PetscScalar H0)
 11: {
 12:   DM           cda;
 13:   PipeField   *x;
 14:   PetscInt     i, start, n;
 15:   Vec          local;
 16:   PetscScalar *coords, c = pipe->R / (GRAV * pipe->A);

 18:   DMGetCoordinateDM(pipe->da, &cda);
 19:   DMGetCoordinatesLocal(pipe->da, &local);
 20:   DMDAVecGetArray(pipe->da, pipe->x, &x);
 21:   DMDAVecGetArrayRead(cda, local, &coords);
 22:   DMDAGetCorners(pipe->da, &start, 0, 0, &n, 0, 0);

 24:   for (i = start; i < start + n; i++) {
 25:     x[i].q = Q0;
 26:     x[i].h = H0 - c * Q0 * PetscAbsScalar(Q0) * coords[i];
 27:   }

 29:   DMDAVecRestoreArray(pipe->da, pipe->x, &x);
 30:   DMDAVecRestoreArrayRead(cda, local, &coords);
 31:   return 0;
 32: }

 34: /* Function evaluations for PIPE    */
 35: /*-------------------------------- */
 36: /* consider using a one-sided higher order fd derivative at boundary. */
 37: static inline PetscScalar dqdx(PipeField *x, PetscInt i, PetscInt ilast, PetscReal dx)
 38: {
 39:   if (i == 0) {
 40:     return (x[i + 1].q - x[i].q) / dx;
 41:   } else if (i == ilast) {
 42:     return (x[i].q - x[i - 1].q) / dx;
 43:   } else {
 44:     return (x[i + 1].q - x[i - 1].q) / (2 * dx);
 45:   }
 46: }

 48: static inline PetscScalar dhdx(PipeField *x, PetscInt i, PetscInt ilast, PetscReal dx)
 49: {
 50:   if (i == 0) {
 51:     return (x[i + 1].h - x[i].h) / dx;
 52:   } else if (i == ilast) {
 53:     return (x[i].h - x[i - 1].h) / dx;
 54:   } else {
 55:     return (x[i + 1].h - x[i - 1].h) / (2 * dx);
 56:   }
 57: }

 59: PetscErrorCode PipeIFunctionLocal_Lax(DMDALocalInfo *info, PetscReal ptime, PipeField *x, PipeField *xdot, PetscScalar *f, Pipe pipe)
 60: {
 61:   PetscInt    i, start, n, ilast;
 62:   PetscReal   a = pipe->a, A = pipe->A, R = pipe->R, c = a * a / (GRAV * A);
 63:   PetscReal   dx = pipe->length / (info->mx - 1), dt = pipe->dt;
 64:   PetscScalar qavg, xold_i, ha, hb, qa, qb;
 65:   PipeField  *xold = pipe->xold;

 67:   DMDAGetCorners(pipe->da, &start, 0, 0, &n, 0, 0);

 69:   /* interior and boundary */
 70:   ilast = start + n - 1;
 71:   for (i = start + 1; i < start + n - 1; i++) {
 72:     qavg = (xold[i + 1].q + xold[i - 1].q) / 2.0;
 73:     qa   = xold[i - 1].q;
 74:     qb   = xold[i + 1].q;
 75:     ha   = xold[i - 1].h;
 76:     hb   = xold[i + 1].h;

 78:     /* xdot[i].q = (x[i].q - old_i)/dt */
 79:     xold_i             = 0.5 * (qa + qb);
 80:     f[2 * (i - 1) + 2] = (x[i].q - xold_i) + dt * (GRAV * pipe->A * dhdx(xold, i, ilast, dx) + pipe->R * qavg * PetscAbsScalar(qavg));

 82:     /* xdot[i].h = (x[i].h - xold_i)/dt */
 83:     xold_i             = 0.5 * (ha + hb);
 84:     f[2 * (i - 1) + 3] = (x[i].h - xold_i) + dt * c * dqdx(xold, i, ilast, dx);
 85:   }

 87:   /* Characteristic equations */
 88:   f[start + 1] = x[start].q - xold[start + 1].q - ((GRAV * A) / a) * (x[start].h - xold[start + 1].h) + dt * R * xold[start + 1].q * PetscAbsScalar(xold[start + 1].q);
 89:   f[2 * ilast] = x[ilast].q - xold[ilast - 1].q + ((GRAV * A) / a) * (x[ilast].h - xold[ilast - 1].h) + dt * R * xold[ilast - 1].q * PetscAbsScalar(xold[ilast - 1].q);
 90:   return 0;
 91: }