Read in a series of bathymetry files and update the cell values

DanGiles · Apr 13, 2021 · 47222a9 · 47222a9
1 parent c636962
commit 47222a9
Show file tree

Hide file tree

Showing 9 changed files with 51 additions and 70 deletions.
diff --git a/sp/EvolveValuesRK2_1.h b/sp/EvolveValuesRK2_1.h
@@ -7,18 +7,7 @@ inline void EvolveValuesRK2_1(const float *dT, const float *Lw_n, //OP_RW //temp
   out[2] = Lw_n[2] * *dT + in[2];
   out[3] = in[3]-in[0];
 
-  //call to ToPhysicalVariables inlined
   float TruncatedH = out[0] < EPS ? EPS : out[0];
   out[0] = TruncatedH;
-  //out[1] = out[1] / TruncatedH;
-  //out[2] = out[2] / TruncatedH;
   out[3] += TruncatedH;
-  /*if (out[0] <= EPS){
-   out[0] = EPS;
-   out[1] = 0.0f;
-   out[2] = 0.0f;
-   out[3] += EPS;
-  } else {
-    out[3] += out[0];
-  }*/  
 }
diff --git a/sp/EvolveValuesRK2_2.h b/sp/EvolveValuesRK2_2.h
@@ -10,15 +10,5 @@ inline void EvolveValuesRK2_2(const float *dT,const float *Lw_1, //OP_RW, discar
   out[3] = values[3]-values[0];
   float TruncatedH = out[0] < EPS ? EPS : out[0];
   out[0] = TruncatedH;
-  //out[1] = out[1] / TruncatedH;
-  //out[2] = out[2] / TruncatedH;
   out[3] += TruncatedH; 
-  /*if(out[0] <= EPS){
-    out[0] = EPS;
-    out[1] = 0.0f;
-    out[2] = 0.0f;
-    out[3] += EPS;
-  } else {
-    out[3] += out[0];
-  }*/
 }
diff --git a/sp/Friction_manning.h b/sp/Friction_manning.h
@@ -7,7 +7,6 @@ inline void Friction_manning(const float *dT,const float *M_n, //OP_RW, discard
   float u = values[1]/TruncatedH;
   float v = values[2]/TruncatedH;
   float speed = sqrt(u*u + v*v);
-  //float speed = sqrt(values[1]*values[1] + values[2]*values[2]);
   F = g* (*M_n * *M_n) *speed;
   F = F/(pow(TruncatedH,4.0f/3.0f));
   //float Fx = F*TruncatedH*values[1];

diff --git a/sp/computeFluxes.h b/sp/computeFluxes.h
@@ -23,26 +23,29 @@ inline void computeFluxes(const float *cellLeft, const float *cellRight,
   dyl = (edgeCenters[1] - leftcellCenters[1]);
   dxr = (edgeCenters[0] - rightcellCenters[0]);
   dyr = (edgeCenters[1] - rightcellCenters[1]);
-
+  
   // Second order Reconstruction
   leftCellValues[0] += alphaleft[0] * ((dxl * leftGradient[0])+(dyl * leftGradient[1]));
-  //leftCellValues[1] += alphaleft[0] * ((dxl * leftGradient[2])+(dyl * leftGradient[3]));
-  //leftCellValues[2] += alphaleft[0] * ((dxl * leftGradient[4])+(dyl * leftGradient[5]));
-  leftCellValues[3] += alphaleft[0] * ((dxl * leftGradient[6])+(dyl * leftGradient[7]));
+  leftCellValues[1] += alphaleft[1] * ((dxl * leftGradient[2])+(dyl * leftGradient[3]));
+  leftCellValues[2] += alphaleft[2] * ((dxl * leftGradient[4])+(dyl * leftGradient[5]));
+  leftCellValues[3] += alphaleft[3] * ((dxl * leftGradient[6])+(dyl * leftGradient[7]));
+
   float TruncatedHL = leftCellValues[0] > EPS ? leftCellValues[0] : EPS;
   uL = leftCellValues[1]/TruncatedHL;
   vL = leftCellValues[2]/TruncatedHL;
   zL = cellLeft[3] - cellLeft[0];
-
+  
   if (!*isRightBoundary) {
     rightCellValues[0] = cellRight[0];
     rightCellValues[1] = cellRight[1];
     rightCellValues[2] = cellRight[2];
     rightCellValues[3] = cellRight[3];
+
+    // Second order Reconstruction
     rightCellValues[0] += alpharight[0] * ((dxr * rightGradient[0])+(dyr * rightGradient[1]));
-    //rightCellValues[1] += alpharight[0] * ((dxr * rightGradient[2])+(dyr * rightGradient[3]));
-    //rightCellValues[2] += alpharight[0] * ((dxr * rightGradient[4])+(dyr * rightGradient[5]));
-    rightCellValues[3] += alpharight[0] * ((dxr * rightGradient[6])+(dyr * rightGradient[7]));
+    rightCellValues[1] += alpharight[1] * ((dxr * rightGradient[2])+(dyr * rightGradient[3]));
+    rightCellValues[2] += alpharight[2] * ((dxr * rightGradient[4])+(dyr * rightGradient[5]));
+    rightCellValues[3] += alpharight[3] * ((dxr * rightGradient[6])+(dyr * rightGradient[7]));
     float TruncatedHR = rightCellValues[0] > EPS ? rightCellValues[0] : EPS;
     uR = rightCellValues[1]/TruncatedHR;
     vR = rightCellValues[2]/TruncatedHR;
@@ -54,30 +57,30 @@ inline void computeFluxes(const float *cellLeft, const float *cellRight,
     float outNormalVelocity = 0.0f;
     float outTangentVelocity = 0.0f;
     rightCellValues[3] = leftCellValues[3];
-    //Outflow
+    // Outflow
     rightCellValues[0] = leftCellValues[0];
+    outNormalVelocity = inNormalVelocity;
     outTangentVelocity = inTangentVelocity;
-    outNormalVelocity =  inNormalVelocity;
     uR = outNormalVelocity * nx - outTangentVelocity * ny;
     vR = outNormalVelocity * ny + outTangentVelocity * nx;
-   // // } else {
-   //     // Wall
-   //     rightCellValues[0] = leftCellValues[0];
-   //     outNormalVelocity =  -1.0f*inNormalVelocity;
-   //     outTangentVelocity = inTangentVelocity;
-   //     uR = outNormalVelocity * nx - outTangentVelocity * ny;
-   //     vR = outNormalVelocity * ny + outTangentVelocity * nx;
-   // // }
+
+    // Wall 
+    /*rightCellValues[0] = leftCellValues[0];
+    outNormalVelocity =  -1.0f*inNormalVelocity;
+    outTangentVelocity = inTangentVelocity;
+    uR = outNormalVelocity * nx - outTangentVelocity * ny;
+    vR = outNormalVelocity * ny + outTangentVelocity * nx;
+    */
 
     rightCellValues[1] = uR*rightCellValues[0];
     rightCellValues[2] = vR*rightCellValues[0];
   }
-  //zL = cellLeft[3] - cellLeft[0];
   zR = cellRight[3] - cellRight[0];
   rightCellValues[3] -= rightCellValues[0];
   leftCellValues[3] -= leftCellValues[0];
   // ------------------------------------------------------------------------------------
   // Audusse Reconstruction(2004) 1st order Source Discretization
+  // ------------------------------------------------------------------------------------
   InterfaceBathy = leftCellValues[3] > rightCellValues[3] ? leftCellValues[3] : rightCellValues[3];
   bathySource[0] =0.5f * g * (leftCellValues[0]*leftCellValues[0]);
   bathySource[1] =0.5f * g * (rightCellValues[0]*rightCellValues[0]);
@@ -92,22 +95,23 @@ inline void computeFluxes(const float *cellLeft, const float *cellRight,
   // Audusse Reconstruction(2005) 2nd order Centered term
   bathySource[2] = -.5f * g *(leftCellValues[0] + cellLeft[0])*(leftCellValues[3] - zL);
   bathySource[3] = -.5f * g *(rightCellValues[0] + cellRight[0])*(rightCellValues[3] - zR);
+
   bathySource[0] *= *edgeLength;
   bathySource[1] *= *edgeLength;
   bathySource[2] *= *edgeLength;
   bathySource[3] *= *edgeLength;
-
+  
   // ------------------------------------------------------------------------------------
   // HLL Riemann Solver
   // Estimation of the wave speeds at the interface.
   float cL = sqrt(g * hL);
   cL = cL > 0.0f ? cL : 0.0f;
   float cR = sqrt(g * hR);
   cR = cR > 0.0f ? cR : 0.0f;
-
+  
   float uLn = uL * edgeNormals[0] + vL * edgeNormals[1];
   float uRn = uR * edgeNormals[0] + vR * edgeNormals[1];
-
+ 
   float unStar = 0.5f * (uLn + uRn) + (cL-cR);
   float cStar = 0.5f * (cL + cR) - 0.25f* (uRn-uLn);
   float sL = (uLn - cL) < (unStar - cStar) ? (uLn - cL) : (unStar - cStar);
@@ -117,7 +121,7 @@ inline void computeFluxes(const float *cellLeft, const float *cellRight,
   float sStar;
   sStar = (sL*hR*(uRn - sR) - sR*hL*(uLn - sL))/
           (hR*(uRn - sR) - hL*(uLn - sL));
-
+  
   if ((leftCellValues[0] <= EPS) && (rightCellValues[0] > EPS)) {
       sL = uRn - 2.0f*cR;
       sR = uRn + cR;
@@ -137,13 +141,13 @@ inline void computeFluxes(const float *cellLeft, const float *cellRight,
   //inlined ProjectedPhysicalFluxes(leftCellValues, Normals, params, LeftFluxes);
   float HuDotN = (leftCellValues[1]) * edgeNormals[0] +
   (leftCellValues[2]) * edgeNormals[1];
-
+  
   LeftFluxes_H = HuDotN;
   LeftFluxes_U = HuDotN * uL;
   LeftFluxes_V = HuDotN * vL;
   // Normal Momentum flux term
   LeftFluxes_N = HuDotN * uLn;
-
+  
   LeftFluxes_U += (.5f * g * edgeNormals[0] ) * ( hL * hL );
   LeftFluxes_V += (.5f * g * edgeNormals[1] ) * ( hL * hL );
   LeftFluxes_N += (.5f * g ) * ( hL * hL );
@@ -153,7 +157,7 @@ inline void computeFluxes(const float *cellLeft, const float *cellRight,
   //inlined ProjectedPhysicalFluxes(rightCellValues, Normals, params, RightFluxes);
   HuDotN = (rightCellValues[1]) * edgeNormals[0] +
   (rightCellValues[2]) * edgeNormals[1];
-
+  
   RightFluxes_H =   HuDotN;
   RightFluxes_U =   HuDotN * uR;
   RightFluxes_V =   HuDotN * vR;
@@ -217,7 +221,7 @@ inline void computeFluxes(const float *cellLeft, const float *cellRight,
   ( t2 * RightFluxes_V ) +
   ( t3 * ( (rightCellValues[2]) -
           (leftCellValues[2]) ) );
-
+    
   // ------------------------------------------------------------------------
   //     HLLC Flux Solver (Huang et al. 2013)
   //     "Well-balanced finite volume scheme for shallow water flooding and drying
@@ -265,7 +269,7 @@ inline void computeFluxes(const float *cellLeft, const float *cellRight,
   out[0] *= *edgeLength;
   out[1] *= *edgeLength;
   out[2] *= *edgeLength;
-
+  
   float maximum = fabs(uLn + cL);
   maximum = maximum > fabs(uLn - cL) ? maximum : fabs(uLn - cL);
   maximum = maximum > fabs(uRn + cR) ? maximum : fabs(uRn + cR);

diff --git a/sp/computeGradient.h b/sp/computeGradient.h
@@ -9,7 +9,7 @@ inline void computeGradient(const float *center,
                             float *q, float *out) //OP_WRITE
 {
   // Least-Squares Gradient Reconstruction
-  //if(center[0]> EPS){
+  if(center[0]> EPS){
     float total, Rhs[8];
     float dh[3], dz[3],du[3], dv[3], weights[3];
     float Gram[2][2], inverse[2][2], delta[3][2];
@@ -109,10 +109,7 @@ inline void computeGradient(const float *center,
     out[6] = (inverse[0][0] * Rhs[6]) + (inverse[0][1] * Rhs[7]);
     out[7] = (inverse[1][0] * Rhs[6]) + (inverse[1][1] * Rhs[7]);
 
-    //if( (cellCenter[0] == nb3Center[0]) && (cellCenter[1] == nb3Center[1])){
-    // op_printf("out %g %g \n", out[0], out[1]);
-    //}
-  /*} else {
+  } else {
     // Gradients for the dry cells are set to zero.
     out[0] = 0.0f;
     out[1] = 0.0f;
@@ -122,7 +119,7 @@ inline void computeGradient(const float *center,
     out[5] = 0.0f;
     out[6] = 0.0f;
     out[7] = 0.0f;
- }*/
+ }
   // Computed the local max and min values for H,U,V,Z
   // q[0] - Hmin , q[1] - Hmax
   // q[2] - Umin , q[3] - Umax

diff --git a/sp/initBathymetry_update.h b/sp/initBathymetry_update.h
@@ -1,10 +1,13 @@
 inline void initBathymetry_update(float *values, const float *zmin, const int *firstTime) {
-  if (*firstTime)
-    values[0] -= values[3];
+    if (*firstTime){
+      //op_printf("values %g %g \n", values[0], values[3]);
+      values[0] -= values[3];
+      values[0] = values[0] < EPS ? EPS : values[0];
+      values[3] += fabs(*zmin) + values[0];
+      //op_printf("After values %g %g \n", values[0], values[3]);
+    } else {
     values[0] = values[0] < EPS ? EPS : values[0];
-    //values[1] *= values[0];
-    //values[2] *= values[0];
     values[3] += fabs(*zmin) + values[0];
-
-  values[0] = values[0] < EPS ? EPS : values[0];
+    }
+   //op_printf("Second Time values %g %g \n", values[0], values[3]);
 }
diff --git a/sp/limiter.h b/sp/limiter.h
@@ -1,14 +1,13 @@
 inline void limiter(const float *q, float *lim,
                     const float *value, const float *gradient,
                     const float *edgecenter1, const float *edgecenter2,
-                    const float *edgecenter3, 
+                    const float *edgecenter3,
                     const float *cellcenter)
 {
 
   float facevalue[3], dx[3], dy[3];
   int i, j;
   float max[3], edgealpha[3];
-
   dx[0] = (edgecenter1[0] - cellcenter[0]);
   dy[0] = (edgecenter1[1] - cellcenter[1]);
   dx[1] = (edgecenter2[0] - cellcenter[0]);
@@ -34,13 +33,13 @@ inline void limiter(const float *q, float *lim,
       edgealpha[i] = 1.0f;
      }
     max[i] = edgealpha[i] < 1.0f ? edgealpha[i] : 1.0f;
-   } 
+   }
    lim[j] = max[0] < max[1] ? max[0] : max[1];
    lim[j] = lim[j] < max[2] ? lim[j]: max[2];
   }
   //lim[0] = lim[0] < lim[1] ? lim[0]: lim[1];
   //lim[0] = lim[0] < lim[2] ? lim[0]: lim[2];
-  lim[0] = lim[0] < lim[3] ? lim[0]: lim[3];
+  //lim[0] = lim[0] < lim[3] ? lim[0]: lim[3];
   } else {
     lim[0] = 0.0f;
     lim[1] = 0.0f;

diff --git a/sp/volna.cpp b/sp/volna.cpp
@@ -388,7 +388,7 @@ int main(int argc, char **argv) {
 
 
       timestep=dT;
-      float Mn = 0.025f;
+      float Mn = 0.013f;
       op_par_loop(Friction_manning, "Friction_manning", cells,
           op_arg_gbl(&dT,1,"float", OP_READ),
           op_arg_gbl(&Mn,1,"float", OP_READ),

diff --git a/sp/volna_writeVTK.cpp b/sp/volna_writeVTK.cpp
@@ -219,19 +219,19 @@ void WriteMeshToVTKAscii(const char* filename, op_dat nodeCoords, int nnode, op_
   for ( i=0; i<ncell; ++i )
     fprintf(fp, "%10.20g\n", values_data[i*N_STATEVAR] + (values_data[i*N_STATEVAR+3] - values_data[i*N_STATEVAR] + *zmin));
   fprintf(fp, "\n");
-  
+
   fprintf(fp, "SCALARS U float 1\n"
               "LOOKUP_TABLE default\n");
   for ( i=0; i<ncell; ++i )
     fprintf(fp, "%10.20g\n", values_data[i*N_STATEVAR+1]/values_data[i*N_STATEVAR]);
   fprintf(fp, "\n");
-  
+
   fprintf(fp, "SCALARS V float 1\n"
               "LOOKUP_TABLE default\n");
   for ( i=0; i<ncell; ++i )
     fprintf(fp, "%10.20g\n", values_data[i*N_STATEVAR+2]/values_data[i*N_STATEVAR]);
   fprintf(fp, "\n");
-  
+
   fprintf(fp, "SCALARS Bathymetry float 1\n"
               "LOOKUP_TABLE default\n");
   for ( i=0; i<ncell; ++i ) {