Changeset 9010

Timestamp:

Oct 30, 2013, 5:34:25 PM (12 years ago)

Author:

davies

Message:

Reduction of velocity spikes in Audusse version?

Location:

trunk/anuga_work/development/gareth

Files:

• experimental/bal_and/swb2_domain.py (modified) (2 diffs)
• experimental/bal_and/swb2_domain_ext.c (modified) (18 diffs)
• tests/shallow_steep_slope/channel_SU_sparse.py (modified) (1 diff)

trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain.py

@@ -64 +64 @@
         # Most of these override the options in config.py
         #------------------------------------------------
-        self.set_CFL(1.0)
+        self.set_CFL(1.00)
         self.set_use_kinematic_viscosity(False)
         self.timestepping_method='rk2'#'rk3'#'euler'#'rk2' 
@@ -76 +76 @@
 
         self.beta_w=1.0
-        self.beta_w_dry=0.0
+        self.beta_w_dry=1.0
         self.beta_uh=1.0
         self.beta_uh_dry=0.0

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

@@ -67 +67 @@
     // Apply limiting of speeds according to the ANUGA manual
     if (*h < epsilon) {
-      *h = max(0.0,*h);  // Could have been negative
+      //*h = max(0.0,*h);  // Could have been negative
       //*h = 0.0;  // Could have been negative
       u = 0.0;
@@ -167 +167 @@
   //  vh_left=0.;
   //}
-  u_left = _compute_speed(&uh_left, &h_left, 
+  // NOTE: using hle instead of h_left here seems to help prevent velocity spikes
+  u_left = _compute_speed(&uh_left, &hle, 
               epsilon, h0, limiting_threshold);
 
@@ -182 +183 @@
   //  vh_right=0.;
   //}
-  u_right = _compute_speed(&uh_right, &h_right, 
+  // NOTE: using hre instead of h_right here seems to help prevent velocity spikes
+  u_right = _compute_speed(&uh_right, &hre, 
                 epsilon, h0, limiting_threshold);
 
@@ -284 +286 @@
       //}else{
           // Standard smoothing term
-      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
+      edgeflux[i] += 1.0*(s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
       //}
 
@@ -339 +341 @@
     // Local variables
     double max_speed, length, inv_area, zl, zr;
-    double h0 = H0;//H0*H0;//H0*0.1;//H0*H0;//H0*H0; // This ensures a good balance when h approaches H0.
+    double h0 = H0*H0;//H0*H0;//H0*0.1;//H0*H0;//H0*H0; // This ensures a good balance when h approaches H0.
     double h_left, h_right, z_half ;  // For andusse scheme
 
@@ -465 +467 @@
             // Don't allow an outward advective flux if the cell centroid stage
             // is < the edge value
-            if(stage_centroid_values[k] < z_half && edgeflux[0] > 0.){
+            if(((stage_centroid_values[k] < z_half)| hle<H0) && edgeflux[0] > 0.){
                 edgeflux[0] = 0.;
                 edgeflux[1] = 0.;
                 edgeflux[2] = 0.;
+                //max_speed=0.;
                 //pressure_flux=0.;
             }
-            
-            if(stage_centroid_values[n] < z_half && edgeflux[0] < 0.){
+            //
+            if(((stage_centroid_values[n] < z_half)| hre<H0) && edgeflux[0] < 0.){
                 edgeflux[0] = 0.;
                 edgeflux[1] = 0.;
                 edgeflux[2] = 0.;
+                //max_speed=0.;
                 //pressure_flux=0.;
             }
@@ -604 +608 @@
      }
 
+    //printf("%e \n", edgeflux_store[3*30*3]);
+
     // Now add up stage, xmom, ymom explicit updates
     for(k=0; k<number_of_elements; k++){
@@ -619 +625 @@
             // GD HACK
             // Option to limit advective fluxes
-            //if(hc > -h0){
+            if(hc > H0){
                 stage_explicit_update[k] += edgeflux_store[ki3+0];
-                // Cut these terms for shallow flows, and protect stationary states from round-off error
                 xmom_explicit_update[k] += edgeflux_store[ki3+1];
                 ymom_explicit_update[k] += edgeflux_store[ki3+2];
-            //}else{
-            //    stage_explicit_update[k] += edgeflux_store[ki3+0];
-            //}
+            }else{
+                stage_explicit_update[k] += edgeflux_store[ki3+0];
+            }
 
 
@@ -927 +932 @@
       // Check if cell k is wet
       //if(stage_centroid_values[k] > elevation_centroid_values[k]){
-      if(stage_centroid_values[k] > elevation_centroid_values[k] + 2.*minimum_allowed_height+epsilon){
+      if(stage_centroid_values[k] > elevation_centroid_values[k] + 1.*minimum_allowed_height+epsilon){
       //if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
           cell_wetness_scale[k] = 1.;  
       }
+
+      min_elevation_edgevalue[k] = min(elevation_edge_values[3*k], 
+                                       min(elevation_edge_values[3*k+1],
+                                           elevation_edge_values[3*k+2]));
+      max_elevation_edgevalue[k] = max(elevation_edge_values[3*k], 
+                                       max(elevation_edge_values[3*k+1],
+                                           elevation_edge_values[3*k+2]));
   }
 
@@ -944 +956 @@
 
     
-    // Try some Froude-number limiting in shallow depths
-    dpth=max(stage_centroid_values[k] - elevation_centroid_values[k], 0.0);
-
-    if(dpth<0.1){
-        // momnorm = momentum^2
-        momnorm=(xmom_centroid_values[k]*xmom_centroid_values[k]+
-                 ymom_centroid_values[k]*ymom_centroid_values[k]);
-        
-        // vel^2 < constant*g*dpth [-- then multiply both sides by dpth^2]
-        if(momnorm > 4*9.81*dpth*dpth*dpth){
-            // Down-scale momentum so that Froude number < constant
-            tmp=sqrt((4*9.81*dpth*dpth*dpth)/momnorm);
-            xmom_centroid_values[k] *=tmp;
-            ymom_centroid_values[k] *=tmp;
-        }
-    }
+    ////// Try some Froude-number limiting in shallow depths
+    //dpth=max(stage_centroid_values[k] - elevation_centroid_values[k], 0.0);
+    ////
+    //if(dpth< max(max_elevation_edgevalue[k]-min_elevation_edgevalue[k],10*minimum_allowed_height)){
+    //    // momnorm = momentum^2
+    //    momnorm=(xmom_centroid_values[k]*xmom_centroid_values[k]+
+    //             ymom_centroid_values[k]*ymom_centroid_values[k]);
+    //    
+    //    // vel^2 < constant*g*dpth [-- then multiply both sides by dpth^2]
+    //    if(momnorm > 1*9.81*dpth*dpth*dpth){
+    //        // Down-scale momentum so that Froude number < constant
+    //        tmp=sqrt((1*9.81*dpth*dpth*dpth)/momnorm);
+    //        xmom_centroid_values[k] *=tmp;
+    //        ymom_centroid_values[k] *=tmp;
+    //    }
+    //}
   }
 
@@ -994 +1006 @@
           ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
 
-          min_elevation_edgevalue[k] = min(elevation_edge_values[3*k], 
-                                           min(elevation_edge_values[3*k+1],
-                                               elevation_edge_values[3*k+2]));
-          max_elevation_edgevalue[k] = max(elevation_edge_values[3*k], 
-                                           max(elevation_edge_values[3*k+1],
-                                               elevation_edge_values[3*k+2]));
 
           // SET HEIGHT IN PLACE
@@ -1033 +1039 @@
           stage_edge_values[k3+1] = stage_centroid_values[k];
           stage_edge_values[k3+2] = stage_centroid_values[k];
+    
+          //xmom_centroid_values[k] = 0.;
+          //ymom_centroid_values[k] = 0.;
+          
           xmom_edge_values[k3]    = xmom_centroid_values[k];
           xmom_edge_values[k3+1]  = xmom_centroid_values[k];
@@ -1114 +1124 @@
           //if(stage_centroid_values[k2]<=elevation_centroid_values[k2]+minimum_allowed_height+epsilon){
           // FIXME: Remove cell_wetness_scale if you don't need it
-          if(k2<0 || (cell_wetness_scale[k2]==-10.0 && stage_centroid_values[k]<stage_centroid_values[k2]) ){
+          if(k2<0 || (cell_wetness_scale[k2]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
               k2 = k ;
           }
           //if(stage_centroid_values[k0]<=max_elevation_edgevalue[k0]+epsilon){
           //if(stage_centroid_values[k0]<=elevation_centroid_values[k0]+minimum_allowed_height+epsilon){
-          if(k0 < 0 || (cell_wetness_scale[k0]==-10.0 && stage_centroid_values[k]<stage_centroid_values[k0])){
+          if(k0 < 0 || (cell_wetness_scale[k0]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
               k0 = k ;
           }
           //if(stage_centroid_values[k1]<=max_elevation_edgevalue[k1]+epsilon){
           //if(stage_centroid_values[k1]<=elevation_centroid_values[k1]+minimum_allowed_height+epsilon){
-          if(k1 < 0 || (cell_wetness_scale[k1]==-10.0 && stage_centroid_values[k]<stage_centroid_values[k1])){
+          if(k1 < 0 || (cell_wetness_scale[k1]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
               k1 = k ;
           }
@@ -1177 +1187 @@
               //stage_edge_values[k3+2] = stage_centroid_values[k];
 
-              dk=max(stage_centroid_values[k]-elevation_centroid_values[k],0.);
+              dk=height_centroid_values[k]; //max(stage_centroid_values[k]-elevation_centroid_values[k],0.);
               height_edge_values[k3] = dk;
               height_edge_values[k3+1] = dk;
@@ -1205 +1215 @@
           h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
           hmin = min(min(h0, min(h1, h2)), hc);
+          hmax = max(max(h0, max(h1, h2)), hc);
           //// GD FIXME: No longer needed? 
           //hfactor = 0.0;
@@ -1214 +1225 @@
             //hfactor=hmin/(hmin + 0.004);
           //}
-          //hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height*10.)), 1.0);
-          hfactor= min(2.0*max(hmin,0.0)/max(hc,minimum_allowed_height), 1.0);
+          //hfactor= min(2.0*max(hmin,0.0)/max(hc,max(1.0*(max_elevation_edgevalue[k]-min_elevation_edgevalue[k]),minimum_allowed_height*1.)), 1.0);
+
+          // Look for strong changes in cell wetness as an indicator of near-wet-dry
+          hfactor= min(2.0*max(hmin,0.0)/max(hc,minimum_allowed_height), 
+                       min(2.0*max(hc,0.)/max(hmax,minimum_allowed_height), 1.0));
+
+          //hfactor=1.0;
+          //hfactor=max(hmin,0.)/(max(hmin,0.)+10.*minimum_allowed_height);
+          //if(hc<minimum_allowed_height*10.) hfactor=0.;
           //if(cell_wetness_scale[k0]==0 | cell_wetness_scale[k1]==0 | cell_wetness_scale[k2]==0) hfactor=0;
+          //hfactor=1.0;
           
           //-----------------------------------
@@ -1250 +1269 @@
           beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
           //beta_tmp = beta_w_dry*0. + (beta_w - beta_w_dry*0.) * hfactor;
+          //beta_tmp=1.0;
         
           
@@ -1291 +1311 @@
         
           beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
-        
+          //beta_tmp=beta_w; 
           
           // Limit the gradient

trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py

@@ -36 +36 @@
 #------------------------------------------------------------------------------
 def topography(x, y):
-        return -x/10.  + 1.*(numpy.sin(x/10.) +abs(y-50.)/10.) -0.*(x>80.) # linear bed slope
+        return -x/10.  #+ 1.*(numpy.sin(x/10.) +abs(y-50.)/10.) -0.*(x>80.) # linear bed slope
 
 def stagetopo(x,y):