diff --git a/sp/EvolveValuesRK2_1.h b/sp/EvolveValuesRK2_1.h index 73a23ad..4a39883 100644 --- 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 index 5969bf4..e766a9a 100644 --- 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 index 7884976..62517b1 100644 --- 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 index 4aeccdd..c7e9159 100644 --- 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,11 +95,12 @@ 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. @@ -104,10 +108,10 @@ inline void computeFluxes(const float *cellLeft, const float *cellRight, 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 index 88ff486..5f19400 100644 --- 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 index 6dadae8..37c3aec 100644 --- 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 index d4b05a0..bcddb70 100644 --- 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 index bde1309..5a9ef73 100644 --- 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 index ac18991..06898f5 100644 --- 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