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