1 #include "pipe.h" 2 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 DM cda; 12 PipeField *x; 13 PetscInt i, start, n; 14 Vec local; 15 PetscScalar *coords, c = pipe->R / (GRAV * pipe->A); 16 17 PetscFunctionBegin; 18 PetscCall(DMGetCoordinateDM(pipe->da, &cda)); 19 PetscCall(DMGetCoordinatesLocal(pipe->da, &local)); 20 PetscCall(DMDAVecGetArray(pipe->da, pipe->x, &x)); 21 PetscCall(DMDAVecGetArrayRead(cda, local, &coords)); 22 PetscCall(DMDAGetCorners(pipe->da, &start, 0, 0, &n, 0, 0)); 23 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 } 28 29 PetscCall(DMDAVecRestoreArray(pipe->da, pipe->x, &x)); 30 PetscCall(DMDAVecRestoreArrayRead(cda, local, &coords)); 31 PetscFunctionReturn(0); 32 } 33 34 /* Function evalutions 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 if (i == 0) { 39 return (x[i + 1].q - x[i].q) / dx; 40 } else if (i == ilast) { 41 return (x[i].q - x[i - 1].q) / dx; 42 } else { 43 return (x[i + 1].q - x[i - 1].q) / (2 * dx); 44 } 45 } 46 47 static inline PetscScalar dhdx(PipeField *x, PetscInt i, PetscInt ilast, PetscReal dx) { 48 if (i == 0) { 49 return (x[i + 1].h - x[i].h) / dx; 50 } else if (i == ilast) { 51 return (x[i].h - x[i - 1].h) / dx; 52 } else { 53 return (x[i + 1].h - x[i - 1].h) / (2 * dx); 54 } 55 } 56 57 PetscErrorCode PipeIFunctionLocal_Lax(DMDALocalInfo *info, PetscReal ptime, PipeField *x, PipeField *xdot, PetscScalar *f, Pipe pipe) { 58 PetscInt i, start, n, ilast; 59 PetscReal a = pipe->a, A = pipe->A, R = pipe->R, c = a * a / (GRAV * A); 60 PetscReal dx = pipe->length / (info->mx - 1), dt = pipe->dt; 61 PetscScalar qavg, xold_i, ha, hb, qa, qb; 62 PipeField *xold = pipe->xold; 63 64 PetscFunctionBegin; 65 PetscCall(DMDAGetCorners(pipe->da, &start, 0, 0, &n, 0, 0)); 66 67 /* interior and boundary */ 68 ilast = start + n - 1; 69 for (i = start + 1; i < start + n - 1; i++) { 70 qavg = (xold[i + 1].q + xold[i - 1].q) / 2.0; 71 qa = xold[i - 1].q; 72 qb = xold[i + 1].q; 73 ha = xold[i - 1].h; 74 hb = xold[i + 1].h; 75 76 /* xdot[i].q = (x[i].q - old_i)/dt */ 77 xold_i = 0.5 * (qa + qb); 78 f[2 * (i - 1) + 2] = (x[i].q - xold_i) + dt * (GRAV * pipe->A * dhdx(xold, i, ilast, dx) + pipe->R * qavg * PetscAbsScalar(qavg)); 79 80 /* xdot[i].h = (x[i].h - xold_i)/dt */ 81 xold_i = 0.5 * (ha + hb); 82 f[2 * (i - 1) + 3] = (x[i].h - xold_i) + dt * c * dqdx(xold, i, ilast, dx); 83 } 84 85 /* Characteristic equations */ 86 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); 87 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); 88 PetscFunctionReturn(0); 89 } 90