Zmin and reconstruction of H+Z

Output of velocities at gauge points

sp/EvolveValuesRK2_1.h

@@ -1,14 +1,24 @@
-inline void EvolveValuesRK2_1(const float *dT, const float *inConservative, //OP_RW //temp
+inline void EvolveValuesRK2_1(const float *dT, const float *Lw_n, //OP_RW //temp
             const float *in, //OP_READ
             float *out) //OP_WRITE
 {
-  out[0] = inConservative[0] * *dT + in[0];
-  out[1] = inConservative[1] * *dT + in[1]*in[0];
-  out[2] = inConservative[2] * *dT + in[2]*in[0];
-  out[3] = in[3];
-
+  out[0] = Lw_n[0] * *dT + in[0];
+  out[1] = Lw_n[1] * *dT + in[1];
+  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[1] = out[1] / TruncatedH;
-  out[2] = out[2] / TruncatedH;
-
+  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];
+  }*/  
 }

sp/EvolveValuesRK2_2.h

@@ -5,10 +5,20 @@ inline void EvolveValuesRK2_2(const float *dT,const float *Lw_1, //OP_RW, discar
 
 {
   out[0] = 0.5*(Lw_1[0] * *dT + w_1[0] + values[0]);
-  out[1] = 0.5*(Lw_1[1] * *dT + w_1[1]*w_1[0] + values[1]*values[0]);
-  out[2] = 0.5*(Lw_1[2] * *dT + w_1[2]*w_1[0] + values[2]*values[0]);
-  out[3] = values[3];
+  out[1] = 0.5*(Lw_1[1] * *dT + w_1[1] + values[1]);
+  out[2] = 0.5*(Lw_1[2] * *dT + w_1[2] + values[2]);
+  out[3] = values[3]-values[0];
   float TruncatedH = out[0] < EPS ? EPS : out[0];
-  out[1] = out[1] / TruncatedH;
-  out[2] = out[2] / TruncatedH;
+  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];
+  }*/
 }

sp/Friction_manning.h

@@ -3,16 +3,22 @@ inline void Friction_manning(const float *dT,const float *M_n, //OP_RW, discard
 
 {
   float F;
-  float speed = sqrt(values[1]*values[1] + values[2]*values[2]);
+  float TruncatedH = values[0] < EPS ? EPS : values[0];
+  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;
-  float TruncatedH = values[0] < 10.f*EPS ? EPS : values[0];
   F = F/(pow(TruncatedH,4.0f/3.0f));
-	if (values[0] <= EPS){
-     values[0] = EPS;
+  //float Fx = F*TruncatedH*values[1];
+  //float Fy = F*TruncatedH*values[2];
+  // Update Momentum 
+  if (values[0] <= EPS){
+     values[0] = TruncatedH;
      values[1] = 0.0f;
      values[2] = 0.0f;
      values[3] = values[3];
-  } else {
+  } else {  
   values[0] = TruncatedH;
   values[1] = values[1] / (1.0f + F * *dT);
   values[2] = values[2] / (1.0f + F * *dT);

sp/Makefile

@@ -115,7 +115,7 @@ endif
 
 
 cuda/volna_kernels_cu.o:	cuda/volna_kernels.cu \
-	toConserved.h EvolveValuesRK2_1.h EvolveValuesRK2_2.h EvolveValuesRK3_1.h EvolveValuesRK3_2.h  EvolveValuesRK3_3.h EvolveValuesRK3_4.h applyConst.h getMaxElevation.h getTotalVol.h \
+	EvolveValuesRK2_1.h EvolveValuesRK2_2.h EvolveValuesRK3_1.h EvolveValuesRK3_2.h  EvolveValuesRK3_3.h EvolveValuesRK3_4.h applyConst.h getMaxElevation.h getTotalVol.h \
 	initBathymetry_formula.h initBathymetry_update.h initBore_select.h initEta_formula.h initGaussianLandslide.h \
 	initU_formula.h initV_formula.h computeGradient.h computeFluxes.h limiter.h SpaceDiscretization.h NumericalFluxes.h simulation_1.h \
 	values_operation2.h cuda/applyConst_kernel.cu cuda/EvolveValuesRK2_1_kernel.cu cuda/EvolveValuesRK2_2_kernel.cu  cuda/EvolveValuesRK3_1_kernel.cu\

sp/NumericalFluxes.h

@@ -1,6 +1,6 @@
 inline void NumericalFluxes(const float *maxEdgeEigenvalues0,
-          const float *maxEdgeEigenvalues1, 
-          const float *maxEdgeEigenvalues2, 
+          const float *maxEdgeEigenvalues1,
+          const float *maxEdgeEigenvalues2,
           const float *EdgeVolumes0,
           const float *EdgeVolumes1,
           const float *EdgeVolumes2,

sp/SpaceDiscretization.h

@@ -1,6 +1,5 @@
 inline void SpaceDiscretization(float *left, //OP_INC
               float *right, //OP_INC
-              const float *cellLeft, const float *cellRight,
               const float *edgeFluxes, //OP_READ
               const float *bathySource, //OP_READ
               const float *edgeNormals, const int *isRightBoundary,