Changeset 9032

Timestamp:

Dec 3, 2013, 10:06:36 AM (11 years ago)

Author:

davies

Message:

Editing hfactor treatment

File:

• trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain_ext.c (modified) (6 diffs)

trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain_ext.c

@@ -181 +181 @@
   soundspeed_left  = sqrt(g*h_left);
   soundspeed_right = sqrt(g*h_right);  
+  //soundspeed_left  = sqrt(g*hle);
+  //soundspeed_right = sqrt(g*hre);  
   
   s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
@@ -235 +237 @@
       // Standard smoothing term
       // edgeflux[i] += 1.0*(s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
-      // Smoothing by stage can cause high velocities / slow draining for nearly dry cells
+      // Smoothing by stage alone can cause high velocities / slow draining for nearly dry cells
       if(i==0) edgeflux[i] += (s_max*s_min)*(max(q_right_rotated[i],ze) - max(q_left_rotated[i],ze));
       if(i==1) edgeflux[i] += (s_max*s_min)*(uh_right - uh_left);
@@ -753 +755 @@
   double dqv[3], qmin, qmax, hmin, hmax, bedmax,bedmin, stagemin;
   double hc, h0, h1, h2, beta_tmp, hfactor, xtmp, ytmp, weight, tmp;
-  double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vel_norm, l1, l2, a_bs, b_bs, area_e, inv_area_e;
+  double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vel_norm, l1, l2, a_tmp, b_tmp, c_tmp,d_tmp;
   
   double *xmom_centroid_store, *ymom_centroid_store; // , *min_elevation_edgevalue, *max_elevation_edgevalue;
@@ -1033 +1035 @@
       hmax = max(max(h0, max(h1, h2)), hc);
 
-      // Look for strong changes in cell depth, or shallow cell depths, as an indicator of near-wet-dry
-      // Reduce beta linearly from 1-0 between depth ratio (hmin/hc) of 0.3-0.1
-      hfactor= max(0., min(5.0*max(hmin,0.0)/max(hc,1.0e-06)-0.5, 
-                           min(5.0*max(hc,0.)/max(hmax,1.0e-06)-0.5, 1.0))
+      // Look for strong changes in cell depth as an indicator of near-wet-dry
+      // Reduce beta linearly from 1-0 between depth ratio (hmin/hc) of [a_tmp , b_tmp]
+      // NOTE: If we have a more 'second order' treatment in near dry areas (e.g. with b_tmp being negative), then
+      //       the water tends to dry more rapidly (which is in agreement with analytical results),
+      //       but is also more 'artefacty' in important cases (tendency for high velocities, etc).
+      //       
+      //hfactor=1.0;
+      a_tmp=0.2; // Highest depth ratio with hfactor=1
+      b_tmp=0.01; // Highest depth ratio with hfactor=0
+      c_tmp=1.0/(a_tmp-b_tmp); 
+      d_tmp= 1.0-(c_tmp*a_tmp);
+      hfactor= max(0., min(c_tmp*max(hmin,0.0)/max(hc,1.0e-06)+d_tmp, 
+                           min(c_tmp*max(hc,0.)/max(hmax,1.0e-06)+d_tmp, 1.0))
                   );
+      //printf("%e, %e, \n", c_tmp, d_tmp);
       //hfactor= max(0., min(5.0*max(hmin,0.0)/max(hmax,1.0e-06)-0.5, 1.0)
       //            );
       //hfactor=1.0;
-      hfactor=min( max(hmin,0.)/(max(hmin,0.)+1.*minimum_allowed_height), hfactor);
+      
+      //hfactor=min( 1.2*max(hmin,0.)/(max(hmin,0.)+1.*minimum_allowed_height), hfactor);
 
       //-----------------------------------
@@ -1089 +1102 @@
       // height
       //-----------------------------------
+      //hfactor=1.0;
+      //hfactor= max(0., min(5.0*max(hmin,0.0)/max(hc,1.0e-06)-0.5, 
+      //                     min(5.0*max(hc,0.)/max(hmax,1.0e-06)-0.5, 1.0))
+      //            );
        
       // Calculate the difference between vertex 0 of the auxiliary 
@@ -1134 +1151 @@
       //                      1.0));
       //hfactor=min( max(hmin,0.)/(max(hmin,0.)+10.*minimum_allowed_height), hfactor);
+      //hfactor=1.0;
       
       //-----------------------------------