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 PetscErrorCode ierr; 13 DM cda; 14 PipeField *x; 15 PetscInt i,start,n; 16 Vec local; 17 PetscScalar *coords,c=pipe->R/(GRAV*pipe->A); 18 19 PetscFunctionBegin; 20 ierr = DMGetCoordinateDM(pipe->da, &cda);CHKERRQ(ierr); 21 ierr = DMGetCoordinatesLocal(pipe->da, &local);CHKERRQ(ierr); 22 ierr = DMDAVecGetArray(pipe->da, pipe->x, &x);CHKERRQ(ierr); 23 ierr = DMDAVecGetArrayRead(cda, local, &coords);CHKERRQ(ierr); 24 ierr = DMDAGetCorners(pipe->da, &start, 0, 0, &n, 0, 0);CHKERRQ(ierr); 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 } 30 31 ierr = DMDAVecRestoreArray(pipe->da, pipe->x, &x);CHKERRQ(ierr); 32 ierr = DMDAVecRestoreArrayRead(cda, local, &coords);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 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 } 49 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 } 60 61 PetscErrorCode PipeIFunctionLocal_Lax(DMDALocalInfo *info,PetscReal ptime,PipeField *x,PipeField *xdot,PetscScalar *f,Pipe pipe) 62 { 63 PetscErrorCode ierr; 64 PetscInt i,start,n,ilast; 65 PetscReal a=pipe->a,A=pipe->A,R=pipe->R,c=a*a/(GRAV*A); 66 PetscReal dx=pipe->length/(info->mx-1),dt=pipe->dt; 67 PetscScalar qavg,xold_i,ha,hb,qa,qb; 68 PipeField *xold=pipe->xold; 69 70 PetscFunctionBegin; 71 ierr = DMDAGetCorners(pipe->da, &start, 0, 0, &n, 0, 0);CHKERRQ(ierr); 72 73 /* interior and boundary */ 74 ilast = start + n - 1; 75 for (i = start + 1; i < start + n - 1; i++) { 76 qavg = (xold[i+1].q + xold[i-1].q)/2.0; 77 qa = xold[i-1].q; qb = xold[i+1].q; 78 ha = xold[i-1].h; 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(0); 93 } 94