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