Changeset 8865 for trunk/anuga_work/development

Timestamp:

May 12, 2013, 8:11:11 PM (12 years ago)

Author:

davies

Message:

Updates to balanced_dev

Location:

trunk/anuga_work/development/gareth

Files:

• experimental/balanced_dev/swb2_domain.py (modified) (7 diffs)
• experimental/balanced_dev/swb2_domain_ext.c (modified) (37 diffs)
• tests/dam_break/plotme.py (modified) (1 diff)
• tests/runup/runup.py (modified) (2 diffs)
• tests/runup_sinusoid/runup_sinusoid.py (modified) (5 diffs)
• tests/shallow_steep_slope/channel_SU_sparse.py (modified) (1 diff)

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

@@ -59 +59 @@
                             number_of_full_triangles = number_of_full_triangles)
        
+        #self.forcing_terms.append(manning_friction_implicit)
         #------------------------------------------------
         # set some defaults
         # Most of these override the options in config.py
         #------------------------------------------------
-        self.set_CFL(1.00)
+        self.set_CFL(1.0)
         self.set_use_kinematic_viscosity(False)
-        self.timestepping_method='rk2' #'euler'#'rk2'#'euler'#'rk2' 
+        self.timestepping_method='rk2'#'rk2'#'euler'#'rk2' 
         # The following allows storage of the negative depths associated with this method
         self.minimum_storable_height=-99999999999.0 
@@ -76 +77 @@
         # extrapolate_second_order_and_limit_by_edge) uses a constant value for
         # all the betas.  
-        self.beta_w=1.0
+        self.beta_w=0.0
         self.beta_w_dry=0.0
-        self.beta_uh=1.0
+        self.beta_uh=0.0
         self.beta_uh_dry=0.0
-        self.beta_vh=1.0
+        self.beta_vh=0.0
         self.beta_vh_dry=0.0
 
@@ -241 +242 @@
 
         # Shortcuts
-        if(domain.timestepping_method!='rk2'):
-            err_mess='ERROR: Timestepping method is ' + domain.timestepping_method +'. '+\
-                     'You need to use rk2 timestepping with this solver, ' +\
-                     ' because there is presently a hack in compute fluxes central. \n'+\
-                     ' The HACK: The timestep will only ' +\
-                     'be recomputed on every 2nd call to compute_fluxes_central, to support'+\
-                     ' correct treatment of wetting and drying'
-
-            raise Exception, err_mess 
+        #if(domain.timestepping_method!='rk2'):
+        #    err_mess='ERROR: Timestepping method is ' + domain.timestepping_method +'. '+\
+        #             'You need to use rk2 timestepping with this solver, ' +\
+        #             ' because there is presently a hack in compute fluxes central. \n'+\
+        #             ' The HACK: The timestep will only ' +\
+        #             'be recomputed on every 2nd call to compute_fluxes_central, to support'+\
+        #             ' correct treatment of wetting and drying'
+
+        #    raise Exception, err_mess 
+
+        if(domain.timestepping_method=='rk2'):
+          timestep_order=2
+        elif(domain.timestepping_method=='euler'):
+          timestep_order=1
+        #elif(domain.timestepping_method=='rk3'):
+        #  timestep_order=3
+        else:
+          err_mess='ERROR: domain.timestepping_method= ' + domain.timestepping_method +' not supported in this solver'
+          raise Exception, err_mess
+
+        #print 'timestep_order=', timestep_order
 
         Stage = domain.quantities['stage']
@@ -260 +273 @@
         flux_timestep = compute_fluxes_ext(timestep,
                                            domain.epsilon,
-                                           domain.H0,
+                                           domain.minimum_allowed_height,
                                            domain.g,
                                            domain.boundary_flux_sum,
@@ -284 +297 @@
                                            domain.max_speed,
                                            int(domain.optimise_dry_cells),
+                                           timestep_order,
                                            Stage.centroid_values,
                                            Xmom.centroid_values,
@@ -293 +307 @@
         #pdb.set_trace()
         #print 'flux timestep: ', flux_timestep, domain.timestep
-        if(flux_timestep == float(sys.maxint)):
-            domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
-                                              domain.boundary_flux_sum[0]*domain.timestep*0.5
-            #print 'dbfi ', domain.boundary_flux_integral, domain.boundary_flux_sum
-            domain.boundary_flux_sum[0]=0.
+        if(domain.timestepping_method=='rk2'):
+          if(flux_timestep == float(sys.maxint)):
+              domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
+                                                domain.boundary_flux_sum[0]*domain.timestep*0.5
+              #print 'dbfi ', domain.boundary_flux_integral, domain.boundary_flux_sum
+              domain.boundary_flux_sum[0]=0.
 
         domain.flux_timestep = flux_timestep
@@ -590 +605 @@
 
 
+#def manning_friction_implicit(domain):
+#    """Apply (Manning) friction to water momentum
+#    Wrapper for c version
+#    """
+#
+#    from shallow_water_ext import manning_friction_flat
+#    from shallow_water_ext import manning_friction_sloped
+#
+#    xmom = domain.quantities['xmomentum']
+#    ymom = domain.quantities['ymomentum']
+#
+#    x = domain.get_vertex_coordinates()
+#    
+#    w = domain.quantities['stage'].centroid_values
+#    z = domain.quantities['elevation'].vertex_values
+#
+#    uh = xmom.centroid_values
+#    vh = ymom.centroid_values
+#    eta = domain.quantities['friction'].centroid_values
+#
+#    xmom_update = xmom.semi_implicit_update
+#    ymom_update = ymom.semi_implicit_update
+#
+#    eps = domain.epsilon
+#    g = domain.g
+#
+#    if domain.use_sloped_mannings:
+#        manning_friction_sloped(g, eps, x, w, uh, vh, z, eta, xmom_update, \
+#                                ymom_update)
+#    else:
+#        manning_friction_flat(g, eps, w, uh, vh, z, eta, xmom_update, \
+#                                #ymom_update)

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

@@ -94 +94 @@
                            double g,
                            double *edgeflux, double *max_speed, 
-                           double *pressure_flux) 
+                           double *pressure_flux, double hc, 
+                           double hc_n, 
+                           double speed_max_last) 
 {
 
@@ -141 +143 @@
   h_left = max(w_left - z,0.);
   uh_left = q_left_rotated[1];
-  //if(h_left>1.e-06){
-  //  u_left = uh_left/max(h_left, 1.0e-06);
-  //}else{
-  //  u_left=0.;
-  //}
-  u_left = _compute_speed(&uh_left, &h_left, 
-                          epsilon, h0, limiting_threshold);
+  if(h_left>h0){
+    u_left = uh_left/max(h_left, 1.0e-06);
+  }else{
+    u_left=0.;
+  }
+  //u_left = _compute_speed(&uh_left, &h_left, 
+  //                      epsilon, h0, limiting_threshold);
 
   w_right = q_right_rotated[0];
   h_right = max(w_right - z, 0.);
   uh_right = q_right_rotated[1];
-  u_right = _compute_speed(&uh_right, &h_right, 
-                           epsilon, h0, limiting_threshold);
-  //if(h_right>1.0e-06){
-  //  u_right = uh_right/max(h_right, 1.0e-06);
-  //}else{
-  //  u_right=0.;
-  //}
+  //u_right = _compute_speed(&uh_right, &h_right, 
+  //                       epsilon, h0, limiting_threshold);
+  if(h_right>h0){
+    u_right = uh_right/max(h_right, 1.0e-06);
+  }else{
+    u_right=0.;
+  }
 
   // Momentum in y-direction
@@ -198 +200 @@
   }
 
+  if( hc < h0){
+    s_max = 0.0;
+  }
+
+
   s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
   if (s_min > 0.0)
@@ -204 +211 @@
   }
   
+  if( hc_n < h0){
+    s_min = 0.0;
+  }
+  
   // Flux formulas
   flux_left[0] = u_left*h_left;
-  flux_left[1] = u_left*uh_left ; //+ 0.5*g*h_left*h_left;
-  flux_left[2] = u_left*vh_left;
+  //if(hc>h0){
+    flux_left[1] = u_left*uh_left; //+ 0.5*g*h_left*h_left;
+    flux_left[2] = u_left*vh_left;
+  //}else{
+  //  flux_left[1] = 0.;//u_left*uh_left; //+ 0.5*g*h_left*h_left;
+  //  flux_left[2] = 0.;//u_left*vh_left;
+  //}  
 
   flux_right[0] = u_right*h_right;
-  flux_right[1] = u_right*uh_right ; //+ 0.5*g*h_right*h_right;
-  flux_right[2] = u_right*vh_right;
+  //if(hc_n>h0){
+    flux_right[1] = u_right*uh_right ; //+ 0.5*g*h_right*h_right;
+    flux_right[2] = u_right*vh_right;
+  //}else{
+  //  flux_right[1] = 0.; //u_right*uh_right ; //+ 0.5*g*h_right*h_right;
+  //  flux_right[2] = 0.; //u_right*vh_right;
+  //}
 
   // Flux computation
@@ -222 +243 @@
     *max_speed = 0.0;
     *pressure_flux = 0.0;
+    //*pressure_flux = 0.5*g*0.5*(h_left*h_left+h_right*h_right);
   } 
   else 
   {
+    // Maximal wavespeed
+    *max_speed = max(s_max, -s_min);
+    
     inverse_denominator = 1.0/denom;
     for (i = 0; i < 3; i++) 
     {
       edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
-      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
+      // Standard smoothing term
+      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
+      //GD HACK: Here, we supress the 'smoothing' part of the flux -- like scaling diffusion by local wave speed
+      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])*
+                      (min(*max_speed/(max(speed_max_last,1e-06)),1.0));//(min(hc/h_left,min(hc_n/h_right,1.0)));
       edgeflux[i] *= inverse_denominator;
     }
@@ -235 +264 @@
     *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
     
-    // Maximal wavespeed
-    *max_speed = max(fabs(s_max), fabs(s_min));
 
     // Rotate back
@@ -274 +301 @@
         double* max_speed_array,
         int optimise_dry_cells, 
+        int timestep_order,
         double* stage_centroid_values,
         double* xmom_centroid_values,
@@ -281 +309 @@
     // Local variables
     double max_speed, length, inv_area, zl, zr;
-    double h0 = H0*H0; // This ensures a good balance when h approaches H0.
-
-    double limiting_threshold = H0; //10 * H0; // Avoid applying limiter below this
+    double h0 = H0*2.0;//H0*H0;//H0*H0; // This ensures a good balance when h approaches H0.
+
+    double limiting_threshold = 10 * H0 ;//H0; //10 * H0; // Avoid applying limiter below this
     // threshold for performance reasons.
     // See ANUGA manual under flux limiting
@@ -297 +325 @@
     double stage_edge_lim, outgoing_mass_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2;
     static long call = 1; // Static local variable flagging already computed flux
+    double speed_max_last, vol;
 
     double *max_bed_edgevalue, *min_bed_edgevalue, *edgeflux_store, *pressuregrad_store;
@@ -314 +343 @@
     local_timestep=timestep;
 
+    //printf("timestep_order %lf \n", timestep_order*1.0);
     // Compute minimum bed edge value on each triangle
+    speed_max_last=0.0;
     for (k = 0; k < number_of_elements; k++){
         max_bed_edgevalue[k] = max(bed_edge_values[3*k], 
                                    max(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
+        speed_max_last=max(speed_max_last, max_speed_array[k]);
+
+        //if(stage_centroid_values[k]<bed_centroid_values[k]){
+        //    printf("Stage < bed");
+        //}
     }
 
+    //printf("SML: %.12lf, %.12lf, %.12lf \n", speed_max_last, fabs(-speed_max_last), max(speed_max_last*2.0, 0.5*speed_max_last)/2.0);
+
+    //printf("SPM: %lf \n", speed_max_last);
 
     // For all triangles
@@ -375 +414 @@
             //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
             // unless the local centroid value is smaller 
-            if(n>=0){
-                if((hc<=H0)&&(hc_n>H0)){
-                    ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
-                }
-                if((hc_n<=H0)&&(hc>H0)){
-                    qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
-                }
-
-            }else{
-                // Treat the boundary case
-                // ??
-            }
+            //if(n>=0){
+            //    if((hc<=H0)&&(hc_n>H0)){
+            //        ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
+            //    }
+            //    if((hc_n<=H0)&&(hc>H0)){
+            //        qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
+            //    }
+
+            //}else{
+            //    // Treat the boundary case
+            //    // ??
+            //}
 
             // Edge flux computation (triangle k, edge i)
@@ -392 +431 @@
                     normals[ki2], normals[ki2 + 1],
                     epsilon, h0, limiting_threshold, g,
-                    edgeflux, &max_speed, &pressure_flux);
+                    edgeflux, &max_speed, &pressure_flux, hc, hc_n, 
+                    speed_max_last);
 
            
@@ -409 +449 @@
             // bedslope_work contains all gravity related terms -- weighting of
             // conservative and non-conservative versions.
-            bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux);
-            //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 
-            pressuregrad_store[ki] = bedslope_work;
-
+            if(hc> -h0){
+              bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux);
+              //bedslope_work = length*(g*(hc*stage_centroid_values[k]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux);
+              //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 
+              pressuregrad_store[ki] = bedslope_work;
+            }else{
+              bedslope_work = length*(g*(hc*ql[0]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux);
+              //bedslope_work = length*(g*(hc*stage_centroid_values[k]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux);
+              //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 
+              pressuregrad_store[ki] = bedslope_work;
+            }
             
             already_computed_flux[ki] = call; // #k Done
@@ -424 +471 @@
                 edgeflux_store[nm3 + 2 ] = edgeflux[2];
 
-                bedslope_work = length*(g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux);
-                //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0];  
-                pressuregrad_store[nm] = bedslope_work;
-
+                if(hc_n> -h0){
+                  bedslope_work = length*(g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux);
+                  //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux);
+                  //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0];  
+                  pressuregrad_store[nm] = bedslope_work;
+                }else{
+                  bedslope_work = length*(g*(hc_n*qr[0]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux);
+                  //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux);
+                  //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0];  
+                  pressuregrad_store[nm] = bedslope_work;
+                }
 
                 already_computed_flux[nm] = call; // #n Done
@@ -439 +493 @@
             if ((tri_full_flag[k] == 1) ) {
 
-                if(call%2==1) max_speed = max_speed_array[k]; // HACK to Ensure that local timestep is the same as the last timestep 
+                if(call%timestep_order!=0) max_speed = max_speed_array[k]; // HACK to Ensure that local timestep is the same as the last timestep 
 
                 if (max_speed > epsilon) {
@@ -465 +519 @@
     // Limit edgefluxes, for mass conservation near wet/dry cells
     for(k=0; k< number_of_elements; k++){
+            //continue;
             hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.);
             // Loop over every edge
@@ -474 +529 @@
                         if(edgeflux_store[3*(3*k+useint)]< 0.){
                             //outgoing_mass_edges+=1.0;
-                            outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)]*local_timestep);
+                            outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)]);
                         }
                     }
+                    outgoing_mass_edges*=local_timestep;
                 }
 
@@ -486 +542 @@
                 // Idea: The cell will not go dry if:
                 // total_outgoing_flux <= cell volume = Area_triangle*hc
-                if((edgeflux_store[ki3]< 0.0) && (outgoing_mass_edges<0.)){
+                vol=areas[k]*hc;
+                if((edgeflux_store[ki3]< 0.0) && (-outgoing_mass_edges> vol)){
                     
                     // This bound could be improved (e.g. we could actually sum
@@ -493 +550 @@
                     // about subsequent changes to the + edgeflux caused by
                     // constraints associated with neighbouring triangles.
-                    tmp = (areas[k]*hc)/(fabs(outgoing_mass_edges)+1.0e-10) ;
+                    tmp = vol/(-(outgoing_mass_edges)) ;
                     if(tmp< 1.0){
                         edgeflux_store[ki3+0]*=tmp;
@@ -528 +585 @@
             // 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
@@ -546 +603 @@
             // GD HACK
             // Compute bed slope term
-            if(hc>H0){
+            if(hc > h0){
                 xmom_explicit_update[k] -= normals[ki2]*pressuregrad_store[ki];
                 ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_store[ki];
@@ -564 +621 @@
     }  // end cell k
 
-    if(call%2==0) timestep=local_timestep;
+    if(call%timestep_order==0) timestep=local_timestep;
 
     // Look for 'dry' edges, and set momentum (& component of momentum_update)
@@ -668 +725 @@
 
       //if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
-      if (hc < minimum_allowed_height) {
+      if (hc < minimum_allowed_height*2.0) {
             // Set momentum to zero and ensure h is non negative
-            xmomc[k] = 0.0;
-            ymomc[k] = 0.0;
+            //xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
+            //ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
+            xmomc[k] = xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
+            ymomc[k] = ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
 
         if (hc <= 0.0){
@@ -678 +737 @@
              //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
              //bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
-             bmin=zc[k]-minimum_allowed_height;
+             bmin=zc[k];//-1.0e-03;
              // Minimum allowed stage = bmin
 
@@ -687 +746 @@
 
                  wc[k] = max(wc[k], bmin); 
-                 printf(" mass_add, %f, %f, %f,%f,%f,%d\n", xc[k], yc[k], wc[k],zc[k],wc[k]-zc[k], k) ;
+                 //printf(" mass_add, %f, %f, %f,%f,%f,%d\n", xc[k], yc[k], wc[k],zc[k],wc[k]-zc[k], k) ;
              
 
@@ -801 +860 @@
   double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
   double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
-  double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
+  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, vect_norm, l1, l2, a_bs, b_bs, area_e, inv_area_e;
@@ -848 +907 @@
       cell_wetness_scale[k] = 0.;
       // Check if cell k is wet
-      if(stage_centroid_values[k] > elevation_centroid_values[k]){
-      //if(stage_centroid_values[k] > elevation_centroid_values[k] + minimum_allowed_height+epsilon){
+      //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] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
           cell_wetness_scale[k] = 1.;  
@@ -935 +994 @@
 
           // Compute bed slope as a reference
-          area_e = (yv2-yv0)*(xv1-xv0) - (yv1-yv0)*(xv2-xv0);
-          inv_area_e=1.0/area_e;
-          a_bs = (yv2 - yv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) -
-                 (yv1 - yv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);
-          a_bs *= inv_area_e;
-
-          b_bs = -(xv2 - xv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) +
-                 (xv1 - xv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);
-          b_bs *= inv_area_e;
+          //area_e = (yv2-yv0)*(xv1-xv0) - (yv1-yv0)*(xv2-xv0);
+          //inv_area_e=1.0/area_e;
+          //a_bs = (yv2 - yv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) -
+          //       (yv1 - yv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);
+          //a_bs *= inv_area_e;
+
+          //b_bs = -(xv2 - xv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) +
+          //       (xv1 - xv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);
+          //b_bs *= inv_area_e;
 
           //printf("slopes: %f, %f \n", a_bs, b_bs);
@@ -989 +1048 @@
           // Take note if the max neighbour bed elevation is greater than the min
           // neighbour stage -- suggests a 'steep' bed relative to the flow
-          //bedmax = max(elevation_centroid_values[k], 
-          //             max(elevation_centroid_values[k0],
-          //                 max(elevation_centroid_values[k1], 
-          //                     elevation_centroid_values[k2])));
+          bedmax = max(elevation_centroid_values[k], 
+                       max(elevation_centroid_values[k0],
+                           max(elevation_centroid_values[k1], 
+                               elevation_centroid_values[k2])));
+          bedmin = min(elevation_centroid_values[k], 
+                       min(elevation_centroid_values[k0],
+                           min(elevation_centroid_values[k1], 
+                               elevation_centroid_values[k2])));
           //stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]), 
           //               min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),
@@ -1028 +1091 @@
 
               if(0==1){
-                // Friction slope type extrapolation -- emulating steady flow
-                // d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth**(4/3)
-                hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0);
-                if(hc>1.0e-06){
-                  tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 
-                             ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03;
-                  tmp /=pow(hc,10.0/3.0);
-                }else{
-                  tmp=0.0;
-                }
-                a = -xmom_centroid_values[k]*tmp;
-                b = -ymom_centroid_values[k]*tmp;
+                //// Friction slope type extrapolation -- emulating steady flow
+                //// d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth**(4/3)
+                //hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0);
+                //if(hc>1.0e-06){
+                //  tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 
+                //             ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03;
+                //  tmp /=pow(hc,10.0/3.0);
+                //}else{
+                //  tmp=0.0;
+                //}
+                //a = -xmom_centroid_values[k]*tmp;
+                //b = -ymom_centroid_values[k]*tmp;
             
-                // Limit gradient to be between the bed slope, and zero
-                if(a*a_bs<0.) a=0.;
-                if(fabs(a)>fabs(a_bs)) a=a_bs;
-                if(b*b_bs<0.) b=0.;
-                if(fabs(b)>fabs(b_bs)) b=b_bs;
-
-
-                dqv[0] = a*dxv0 + b*dyv0;
-                dqv[1] = a*dxv1 + b*dyv1;
-                dqv[2] = a*dxv2 + b*dyv2;
-
-                stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
-                stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
-                stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
+                //// Limit gradient to be between the bed slope, and zero
+                //if(a*a_bs<0.) a=0.;
+                //if(fabs(a)>fabs(a_bs)) a=a_bs;
+                //if(b*b_bs<0.) b=0.;
+                //if(fabs(b)>fabs(b_bs)) b=b_bs;
+
+
+                //dqv[0] = a*dxv0 + b*dyv0;
+                //dqv[1] = a*dxv1 + b*dyv1;
+                //dqv[2] = a*dxv2 + b*dyv2;
+
+                //stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
+                //stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
+                //stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
               }else{
                 // Constant stage extrapolation
@@ -1185 +1248 @@
           ////if (hc>0.0)
           //{
-            hfactor = 1.0 ;//hmin/(hmin + 0.004);
+          //  hfactor = 1.0 ;//hmin/(hmin + 0.004);
             //hfactor=hmin/(hmin + 0.004);
           //}
-          hfactor= min(2.0*max(hmin,0)/max(hc,1.0e-06), 1.0);
+          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,max(0.5*(bedmax-bedmin),minimum_allowed_height)), 1.0);
           
           //-----------------------------------
@@ -1252 +1316 @@
           //if( (area2>0) ){ 
           //if(count_wet_neighbours[k]>0){
-              limit_gradient(dqv, qmin, qmax, beta_tmp);
+          limit_gradient(dqv, qmin, qmax, beta_tmp);
           //}
           //}
-              stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
-              stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
-              stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
+          stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
+          stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
+          stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
           // IDEA: extrapolate assuming slope is that of steady flow??
           // GD HACK - FIXME
@@ -1746 +1810 @@
     
   double timestep, epsilon, H0, g;
-  int optimise_dry_cells;
+  int optimise_dry_cells, timestep_order;
     
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOiOOOOO",
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOiiOOOOO",
             &timestep,
             &epsilon,
@@ -1774 +1838 @@
             &max_speed_array,
             &optimise_dry_cells,
+            &timestep_order,
             &stage_centroid_values,
             &xmom_centroid_values,
@@ -1841 +1906 @@
                      (double*) max_speed_array -> data,
                      optimise_dry_cells, 
+                     timestep_order,
                      (double*) stage_centroid_values -> data, 
                      (double*) xmom_centroid_values -> data, 
@@ -1871 +1937 @@
 
   
-  h0 = H0*H0; // This ensures a good balance when h approaches H0.
-              // But evidence suggests that h0 can be as little as
-              // epsilon!
-              
+  h0 = H0;            
   limiting_threshold = 10*H0; // Avoid applying limiter below this
                               // threshold for performance reasons.
@@ -1890 +1953 @@
                          (double*) edgeflux -> data, 
                          &max_speed,
-             &pressure_flux);
+             &pressure_flux,
+                         ((double*) normal -> data)[0],
+                         ((double*) normal -> data)[1], 
+                         ((double*) normal -> data)[1] 
+             );
   
   return Py_BuildValue("d", max_speed);  

trunk/anuga_work/development/gareth/tests/dam_break/plotme.py

@@ -9 +9 @@
 p2_st=util.get_centroids(p_st)
 
-p_dev = util.get_output('dam_break_20130320_001049/dam_break.sww', 0.001)
+p_dev = util.get_output('dam_break_20130504_145918/dam_break.sww', 0.001)
 p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
 

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

@@ -34 +34 @@
 #------------------
 
+scale_me=100.0
+
 def topography(x,y):
-        return -x/2 #+0.05*numpy.sin((x+y)*50.0) #+0.1*(numpy.random.rand(len(x)) -0.5)       # Linear bed slope + small random perturbation 
+        return -x/2*scale_me #+0.05*numpy.sin((x+y)*50.0) #+0.1*(numpy.random.rand(len(x)) -0.5)       # Linear bed slope + small random perturbation 
 
 def stagefun(x,y):
-    stg=-0.2*(x<0.5) -0.1*(x>=0.5)
-    #stg=-0.2 # Stage
+    #stg=(-0.2*(x<0.5) -0.1*(x>=0.5))*scale_me
+    stg=-0.2*scale_me # Stage
     #topo=topography(x,y) #Bed
     return stg #*(stg>topo) + topo*(stg<=topo)
@@ -55 +57 @@
 Br=anuga.Reflective_boundary(domain)            # Solid reflective wall
 Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
-Bd=anuga.Dirichlet_boundary([-0.2,0.,0.])       # Constant boundary values -- not used in this example
+Bd=anuga.Dirichlet_boundary([-0.2*scale_me,0.,0.])       # Constant boundary values -- not used in this example
 #Bw=anuga.Time_boundary(domain=domain,
 #       f=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,0.0,0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.

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

@@ -24 +24 @@
 domain=Domain(points,vertices,boundary)    # Create Domain
 domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
-
+#domain.set_timestepping_method('euler')
 #------------------
 # Define topography
 #------------------
 scale_me=1.0
+boundary_ws=-0.1
+init_ws=-0.2
+bumpiness=10. # Higher = shorter wavelength oscillations in topography
 #domain.minimum_allowed_height=domain.minimum_allowed_height*scale_me # Seems needed to make the algorithms behave
 
 def topography(x,y):
-        return (-x/2.0 +0.05*numpy.sin((x+y)*50.0))*scale_me
+        return (-x/2.0 +0.05*numpy.sin((x+y)*bumpiness))*scale_me
 
 def stagefun(x,y):
-    stge=-0.2*scale_me #+0.1*(x>0.9) 
+    stge=init_ws*scale_me# +0.1*(x>0.9)*scale_me 
     #topo=topography(x,y) 
     return stge#*(stge>topo) + (topo)*(stge<=topo)
@@ -41 +44 @@
 domain.set_quantity('elevation',topography)     # Use function for elevation
 domain.get_quantity('elevation').smooth_vertex_values() 
-domain.set_quantity('friction',0.03)             # Constant friction
+domain.set_quantity('friction',0.00)             # Constant friction
 
 #def frict_change(x,y):
@@ -64 +67 @@
 Bd=anuga.Dirichlet_boundary([-0.1*scale_me,0.,0.])       # Constant boundary values -- not used in this example
 def waveform(t):
-    return -0.1 #(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1
+    return boundary_ws*scale_me
+    #return -0.2*scale_me #-0.1 #(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1
 Bt2=anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain,waveform)
 #Bw=anuga.Time_boundary(domain=domain,
@@ -78 +82 @@
 #------------------------------
 
-for t in domain.evolve(yieldstep=0.200,finaltime=40.00):
+for t in domain.evolve(yieldstep=0.2,finaltime=120.00):
     print domain.timestepping_statistics()
-    print domain.boundary_flux_integral
+    #print domain.boundary_flux_integral
     xx = domain.quantities['xmomentum'].centroid_values
     yy = domain.quantities['ymomentum'].centroid_values
@@ -88 +92 @@
     vv = ( (xx/dd)**2 + (yy/dd)**2)**0.5
     vv = vv*(dd>1.0e-03)
-    print 'Peak velocity is: ', vv.max(), vv.argmax()
+    print 'Peak velocity is: ', vv.max(), vv.argmax(), dd[vv.argmax()]
     print 'Volume is', sum(dd_raw*domain.areas)    
-    print 'Volume less flux int', sum(dd_raw*domain.areas) - domain.boundary_flux_integral
+    #print 'Volume less flux int', sum(dd_raw*domain.areas) - domain.boundary_flux_integral
 
 

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):