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 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: {
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; qb = xold[i+1].q;
74: ha = xold[i-1].h; hb = xold[i+1].h;
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));
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: }
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: return 0;
89: }