Changeset 9014

Timestamp:

Nov 3, 2013, 6:33:53 PM (11 years ago)

Author:

davies

Message:

Updates to improve well balancing at wet-dry edge

Location:

trunk/anuga_work/development/gareth

Files:

• experimental/bal_and/swb2_domain.py (modified) (1 diff)
• experimental/bal_and/swb2_domain_ext.c (modified) (13 diffs)
• tests/runup_sinusoid/runup_sinusoid.py (modified) (1 diff)
• tests/shallow_steep_slope/channel_SU_2plot.py (modified) (4 diffs)
• tests/shallow_steep_slope/channel_SU_sparse.py (modified) (2 diffs)

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

@@ -76 +76 @@
 
         self.beta_w=1.0
-        self.beta_w_dry=1.0
+        self.beta_w_dry=0.0
         self.beta_uh=1.0
         self.beta_uh_dry=0.0

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

@@ -68 +68 @@
     if (*h < epsilon) {
       //*h = max(0.0,*h);  // Could have been negative
-      //*h = 0.0;  // Could have been negative
       u = 0.0;
     } else {
@@ -146 +145 @@
 
 
-  //z=max(zl,zr);
-  //z = 0.5*(z_left + z_right); // Average elevation values. 
-                              // Even though this will nominally allow 
-                              // for discontinuities in the elevation data, 
-                              // there is currently no numerical support for
-                              // this so results may be strange near 
-                              // jumps in the bed.
-
   // Compute speeds in x-direction
   w_left = q_left_rotated[0];          
   uh_left=q_left_rotated[1];
   vh_left=q_left_rotated[2];
-  if(h_left>1.0e-03){
+  if(hle>1.0e-03){
     u_left = uh_left/(hle+0.0e-03) ; //max(h_left, 1.0e-06);
     uh_left=h_left*u_left;
@@ -167 +158 @@
     vh_left=0.;
   }
-  // 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);
@@ -174 +165 @@
   uh_right = q_right_rotated[1];
   vh_right = q_right_rotated[2];
-  if(h_right>1.0e-03){
+  if(hre>1.0e-03){
     u_right = uh_right/(hre+0.0e-03);//max(h_right, 1.0e-06);
     uh_right=h_right*u_right;
@@ -180 +171 @@
   }else{
     u_right=0.;
-    uh_right=0.;//h_right*u_right;
+    uh_right=0.;
     vh_right=0.;
   }
-  // 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);
 
   
-
-  // Limit y-momentum if necessary 
-  // Leaving this out, improves speed significantly (Ole 27/5/2009)
-  // All validation tests pass, so do we really need it anymore? 
-  //_compute_speed(&vh_left, &h_left, 
-  //    epsilon, h0, limiting_threshold);
-  //_compute_speed(&vh_right, &h_right, 
-  //     epsilon, h0, limiting_threshold);
-
   // Maximal and minimal wave speeds
   soundspeed_left  = sqrt(g*h_left);
@@ -441 +422 @@
             // Audusse magic 
             z_half=max(zl,zr);
+
             h_left= max(hle+zl-z_half,0.);
             h_right= max(hre+zr-z_half,0.);
+
             //if (fabs(zl-zr)>1.0e-10) {
             //    report_python_error(AT,"Discontinuous Elevation");
@@ -467 +450 @@
             // Don't allow an outward advective flux if the cell centroid stage
             // is < the edge value
-            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)| hre<H0) && edgeflux[0] < 0.){
-                edgeflux[0] = 0.;
-                edgeflux[1] = 0.;
-                edgeflux[2] = 0.;
-                //max_speed=0.;
-                //pressure_flux=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)| hre<H0) && edgeflux[0] < 0.){
+            //    edgeflux[0] = 0.;
+            //    edgeflux[1] = 0.;
+            //    edgeflux[2] = 0.;
+            //    //max_speed=0.;
+            //    //pressure_flux=0.;
+            //}
 
 
@@ -641 +624 @@
             // GD HACK
             // Compute bed slope term
-            //if(hc > -h0){
+            //if(hc > H0*0.5){
                 xmom_explicit_update[k] -= normals[ki2]*pressuregrad_store[ki];
                 ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_store[ki];
@@ -932 +915 @@
       // Check if cell k is wet
       //if(stage_centroid_values[k] > elevation_centroid_values[k]){
-      if(stage_centroid_values[k] > elevation_centroid_values[k] + 1.*minimum_allowed_height+epsilon){
+      if(stage_centroid_values[k] > elevation_centroid_values[k] + 1.0*minimum_allowed_height){
       //if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
           cell_wetness_scale[k] = 1.;  
@@ -1124 +1107 @@
           //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]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
+          if(k2<0 || (cell_wetness_scale[k2]==-10.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]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
+          if(k0 < 0 || (cell_wetness_scale[k0]==-10.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]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
+          //if(k1 < 0 || (cell_wetness_scale[k1]==0.0 && stage_centroid_values[k]<elevation_edge_values[3*k+1])){
+          if(k1 < 0 || (cell_wetness_scale[k1]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
               k1 = k ;
           }
@@ -1179 +1163 @@
            
           //// Treat triangles with zero or 1 wet neighbours (area2 <=0.) 
-          if ((area2 <= 0.) )// || (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
+          if ((area2 <= 0.) )//|| (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
           {
 
@@ -1214 +1198 @@
           h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
           h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
+          //h0=0.;
+          //h1=0.;
+          //h2=0.;
+
+          //if(k!=k0) h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
+          //if(k!=k1) h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
+          //if(k!=k2) 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);
@@ -1230 +1222 @@
           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= min(2.0*max(hc,0.)/max(hmax,minimum_allowed_height), 1.0);
+          //hfactor=0.;
+          //if(hmin==hc) hfactor=0.;
+          //hfactor= min(2.0*max(hmin,0.0)/max(hc,minimum_allowed_height+bedmax-bedmin), 
+          //             min(2.0*max(hc,0.)/max(hmax,minimum_allowed_height+bedmax-bedmin), 1.0));
 
           //hfactor=1.0;
-          //hfactor=max(hmin,0.)/(max(hmin,0.)+10.*minimum_allowed_height);
+          //hfactor=min( max(hmin,0.)/(max(hmin,0.)+10.*minimum_allowed_height), 1.0);
           //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;

trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

@@ -45 +45 @@
 bumpiness=50. # Higher = shorter wavelength oscillations in topography
 tstep=0.2
-lasttime=90.
+lasttime=150.
 
 #domain.minimum_allowed_height=domain.minimum_allowed_height*scale_me # Seems needed to make the algorithms behave

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

@@ -9 +9 @@
 import scipy
 from matplotlib import pyplot as pyplot
-from anuga.utilities import plot_utils as util
+#from anuga.utilities import plot_utils as util
+from bal_and import plot_utils as util
 #--------------
 # Get variables
@@ -29 +30 @@
 mann=0.03 # Manning's coef
 bedslope=-0.1
-fluxin=20./100. #The momentum flux at the upstream boundary ( = discharge / width)
+fluxin=0.5/100. #The momentum flux at the upstream boundary ( = discharge / width)
 
 #---------------------------------------------
@@ -46 +47 @@
 pyplot.ion()
 
-line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,( (p2.stage[:,v]-p2.elev[:,v]).max(),(p2.stage[:,v]-p2.elev[v]).min() ) )
+line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,( (p2.stage[:,v]-p2.elev[1,v]).max(),(p2.stage[:,v]-p2.elev[1,v]).min() ) )
 line.set_label('Numerical')
 pyplot.plot( (0,100),(dana,dana), 'r',label='Analytical' )
+pyplot.yscale('log')
 pyplot.legend()
 #pyplot.plot(x[v],elev[v],'green')
 for i in range(p2.xmom.shape[0]):
     line.set_xdata(p2.x[v])
-    line.set_ydata(p2.stage[i,v]-p2.elev[v])
+    line.set_ydata(p2.stage[i,v]-p2.elev[1,v])
     pyplot.draw()
     pyplot.title('Flow depth: should be converging to steady uniform flow ' )
@@ -97 +99 @@
 pyplot.plot(p2.y[v],p2.xvel[100,v], 'o', label='Numerical')
 pyplot.plot((0,100),(uana,uana),label='Analytical')
-pyplot.ylim([1.0,3.0])
+#pyplot.ylim([1.0,3.0])
 pyplot.xlabel('Distance along the line x==50 (across the slope) (m)')
 pyplot.ylabel('Velocity m/s')

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):
@@ -54 +54 @@
 # Setup boundary conditions
 #------------------------------------------------------------------------------
-#Bi = anuga.Dirichlet_boundary([0.06309625, 0.00, 0.00]) # Inflow for steady uniform flow with a depth integrated velocity of 0.1
+#Bi = anuga.Dirichlet_boundary([0.06309625, Qin/100., 0.00]) # Inflow for steady uniform flow with a depth integrated velocity of 0.1
 
 Br = anuga.Reflective_boundary(domain) # Solid reflective wall