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 { 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(0); 33 } 34 35 /* Function evalutions for PIPE */ 36 /*-------------------------------- */ 37 /* consider using a one-sided higher order fd derivative at boundary. */ 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 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 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; qb = xold[i+1].q; 76 ha = xold[i-1].h; hb = xold[i+1].h; 77 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)); 81 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 } 86 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 PetscFunctionReturn(0); 91 } 92