Actual source code: pipeImpls.c
petsc-3.9.3 2018-07-02
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: {
13: DM cda;
14: PipeField *x;
15: PetscInt i,start,n;
16: Vec local;
17: PetscScalar *coords,c=pipe->R/(GRAV*pipe->A);
20: DMGetCoordinateDM(pipe->da, &cda);
21: DMGetCoordinatesLocal(pipe->da, &local);
22: DMDAVecGetArray(pipe->da, pipe->x, &x);
23: DMDAVecGetArrayRead(cda, local, &coords);
24: DMDAGetCorners(pipe->da, &start, 0, 0, &n, 0, 0);
25:
26: for (i = start; i < start + n; i++) {
27: x[i].q = Q0;
28: x[i].h = H0 - c * Q0 * PetscAbsScalar(Q0) * coords[i];
29: }
31: DMDAVecRestoreArray(pipe->da, pipe->x, &x);
32: DMDAVecRestoreArrayRead(cda, local, &coords);
33: return(0);
34: }
36: /* Function evalutions for PIPE */
37: /*-------------------------------- */
38: /* consider using a one-sided higher order fd derivative at boundary. */
39: PETSC_STATIC_INLINE PetscScalar dqdx(PipeField *x,PetscInt i,PetscInt ilast,PetscReal dx)
40: {
41: if (i == 0) {
42: return (x[i+1].q - x[i].q) / dx;
43: } else if (i == ilast) {
44: return (x[i].q - x[i-1].q) / dx;
45: } else {
46: return (x[i+1].q - x[i-1].q) / (2*dx);
47: }
48: }
50: PETSC_STATIC_INLINE PetscScalar dhdx(PipeField *x,PetscInt i,PetscInt ilast,PetscReal dx)
51: {
52: if (i == 0) {
53: return (x[i+1].h - x[i].h) / dx;
54: } else if (i == ilast) {
55: return (x[i].h - x[i-1].h) / dx;
56: } else {
57: return (x[i+1].h - x[i-1].h) / (2*dx);
58: }
59: }
61: PetscErrorCode PipeIFunctionLocal(DMDALocalInfo *info,PetscReal ptime,PipeField *x,PipeField *xdot,PipeField *f,Pipe pipe)
62: {
64: PetscInt i, start, n, ilast;
65: PetscReal c = (pipe->a * pipe->a) / (GRAV * pipe->A);
66: PetscReal dx = pipe->length / (info->mx-1);
67: PetscScalar qavg;
70: DMDAGetCorners(pipe->da, &start, 0, 0, &n, 0, 0);
71:
72: /* interior and boundary */
73: ilast = start + n -1;
74: for (i = start; i < start + n; i++) {
75: if (i == start || i == ilast) {
76: qavg = x[i].q;
77: } else {
78: qavg = (x[i+1].q + x[i-1].q)/2.0; /* ok for single pipe with DM_BOUNDARY_GHOSTED, but mem corrupt for pipes! */
79: }
80: f[i].q = xdot[i].q + GRAV * pipe->A * dhdx(x, i, ilast, dx) + pipe->R * qavg * PetscAbsScalar(qavg);
81: f[i].h = xdot[i].h + c * dqdx(x, i, ilast, dx);
82: }
84: /* up-stream boundary */
85: if (info->xs == 0) {
86: if (pipe->boundary.Q0 == PIPE_CHARACTERISTIC) {
87: f[0].h = x[0].h - pipe->boundary.H0;
88: } else {
89: f[0].q = x[0].q - pipe->boundary.Q0;
90: }
91: }
92:
93: /* down-stream boundary */
94: if (start + n == info->mx) {
95: if (pipe->boundary.HL == PIPE_CHARACTERISTIC) {
96: f[info->mx-1].q = x[info->mx-1].q - pipe->boundary.QL;
97: } else {
98: f[info->mx-1].h = x[info->mx-1].h - pipe->boundary.HL;
99: }
100: }
101: return(0);
102: }