Changeset 8772

Timestamp:

Mar 20, 2013, 12:01:46 AM (12 years ago)

Author:

davies

Message:

Updates to gd experimental

Location:

trunk/anuga_work/development/gareth

Files:

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

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

@@ -2 +2 @@
 from anuga.shallow_water.shallow_water_domain import *
 from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain
-
+import numpy
 
 ##############################################################################
@@ -77 +77 @@
         # all the betas.  
         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=1.0
+        self.beta_uh_dry=0.0
         self.beta_vh=1.0
-        self.beta_vh_dry=1.0
+        self.beta_vh_dry=0.0
 
         #self.optimise_dry_cells=True
@@ -104 +104 @@
         #self.forcing_terms.append(manning_friction_explicit)
         #self.forcing_terms.remove(manning_friction_implicit)
+
+        # Keep track of the fluxes through the boundaries
+        self.boundary_flux_integral=numpy.ndarray(1)
+        self.boundary_flux_integral[0]=0. 
+        self.boundary_flux_sum=numpy.ndarray(1)
+        self.boundary_flux_sum[0]=0.
 
         print '##########################################################################'
@@ -256 +262 @@
                                            domain.H0,
                                            domain.g,
+                                           domain.boundary_flux_sum,
                                            domain.neighbours,
                                            domain.neighbour_edges,
@@ -285 +292 @@
         #import pdb
         #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.
 
         domain.flux_timestep = flux_timestep
+
 
 

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

@@ -139 +139 @@
   // Compute speeds in x-direction
   w_left = q_left_rotated[0];          
-  h_left = w_left - z;
+  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);
+                          epsilon, h0, limiting_threshold);
 
   w_right = q_right_rotated[0];
-  h_right = w_right - z;
+  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);
+                           epsilon, h0, limiting_threshold);
+  //if(h_right>1.0e-06){
+  //  u_right = uh_right/max(h_right, 1.0e-06);
+  //}else{
+  //  u_right=0.;
+  //}
 
   // Momentum in y-direction
@@ -157 +167 @@
   // 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);
+  //_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
@@ -219 +229 @@
     {
       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]);
+      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
       edgeflux[i] *= inverse_denominator;
     }
@@ -242 +252 @@
         double H0,
         double g,
+        double* boundary_flux_sum,
         long* neighbours,
         long* neighbour_edges,
@@ -272 +283 @@
     double h0 = H0*H0; // This ensures a good balance when h approaches H0.
 
-    double limiting_threshold = 10 * H0; // Avoid applying limiter below this
+    double limiting_threshold = H0; //10 * H0; // Avoid applying limiter below this
     // threshold for performance reasons.
     // See ANUGA manual under flux limiting
-    int k, i, m, n,j;
+    int k, i, m, n,j, ii;
     int ki, nm = 0, ki2,ki3, nm3; // Index shorthands
 
@@ -281 +292 @@
     double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
     double stage_edges[3];//Work array
-    double bedslope_work;
+    double bedslope_work, local_timestep;
     int neighbours_wet[3];//Work array
     int useint;
-    double stage_edge_lim, outgoing_flux_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2;
+    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
 
@@ -300 +311 @@
     memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double));
     memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
+
 
     // Compute minimum bed edge value on each triangle
@@ -424 +436 @@
             // every 2nd call.  This works for rk2 timestepping only, and is
             // needed for correct wetting and drying treatment
-            if ((tri_full_flag[k] == 1) && (call%2==0)) {
+            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 (max_speed > epsilon) {
                     // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
@@ -449 +464 @@
     // Limit edgefluxes, for mass conservation near wet/dry cells
     for(k=0; k< number_of_elements; k++){
-        // Only check cells that are **near dry**
-        // Found I could get problems by checking all cells
-        if(stage_centroid_values[k] <= max_bed_edgevalue[k] + H0){
             hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.);
             // Loop over every edge
@@ -457 +469 @@
                 if(i==0){
                     // Add up the outgoing flux through the cell
-                    outgoing_flux_edges=0.0;
+                    outgoing_mass_edges=0.0;
                     for(useint=0; useint<3; useint++){
                         if(edgeflux_store[3*(3*k+useint)]< 0.){
-                            //outgoing_flux_edges+=1.0;
-                            outgoing_flux_edges+=(edgeflux_store[3*(3*k+useint)]*timestep);
+                            //outgoing_mass_edges+=1.0;
+                            outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)]*timestep);
                         }
                     }
@@ -473 +485 @@
                 // Idea: The cell will not go dry if:
                 // total_outgoing_flux <= cell volume = Area_triangle*hc
-                if((edgeflux_store[ki3]< 0.0) && (outgoing_flux_edges<0.)){
+                if((edgeflux_store[ki3]< 0.0) && (outgoing_mass_edges<0.)){
                     
                     // This bound could be improved (e.g. we could actually sum
@@ -480 +492 @@
                     // about subsequent changes to the + edgeflux caused by
                     // constraints associated with neighbouring triangles.
-                    tmp = (areas[k]*hc)/(fabs(outgoing_flux_edges)+1.0e-10) ;
+                    tmp = (areas[k]*hc)/(fabs(outgoing_mass_edges)+1.0e-10) ;
                     if(tmp< 1.0){
-                        // Idea: we should be close to conserving mass if this holds, 
-                        // --- consider CFL condition. 
-                        tmp2 = min(hc/((stage_edge_values[ki]-bed_edge_values[ki])+1.0e-10), 1.0);
-                        //tmp2=(xmom_centroid_values[k]*xmom_centroid_values[k]+ymom_centroid_values[k]*ymom_centroid_values[k]);
-                        //tmp2/=(xmom_edge_values[ki]*xmom_edge_values[ki]+ ymom_edge_values[ki]*ymom_edge_values[ki]+1.0e-10);
-                        //tmp2 = sqrt(tmp2);
-                        //tmp2 = min(tmp2, 1.0);
-                        edgeflux_store[ki3+0]*=tmp2;
-                        edgeflux_store[ki3+1]*=tmp2;
-                        edgeflux_store[ki3+2]*=tmp2;
+                        edgeflux_store[ki3+0]*=tmp;
+                        edgeflux_store[ki3+1]*=tmp;
+                        edgeflux_store[ki3+2]*=tmp;
 
                         // Compute neighbour edge index
@@ -498 +503 @@
                             nm = 3*n + neighbour_edges[ki];
                             nm3 = nm*3;
-                            edgeflux_store[nm3+0]*=tmp2;
-                            edgeflux_store[nm3+1]*=tmp2;
-                            edgeflux_store[nm3+2]*=tmp2;
+                            edgeflux_store[nm3+0]*=tmp;
+                            edgeflux_store[nm3+1]*=tmp;
+                            edgeflux_store[nm3+2]*=tmp;
                         }
                     }
                 }
-            }
         }
      }
-
 
     // Now add up stage, xmom, ymom explicit updates
@@ -520 +523 @@
             ki2=ki*2;
             ki3 = ki*3;
+            n=neighbours[ki];
 
             // GD HACK
@@ -533 +537 @@
 
 
+            // If neighbour is a boundary condition, add the flux to the boundary_flux_integral
+            if(n<0){ 
+                boundary_flux_sum[0] += edgeflux_store[ki3];
+                //printf("%f, %f, %d \n", boundary_flux_sum[0], timestep, ki3);
+            }
+
             // GD HACK
             // Compute bed slope term
@@ -552 +562 @@
         ymom_explicit_update[k] *= inv_area;
     }  // end cell k
+
+    //if(call%2==0) timestep=local_timestep;
 
     // Look for 'dry' edges, and set momentum (& component of momentum_update)
@@ -790 +802 @@
   double dqv[3], qmin, qmax, hmin, hmax, bedmax, 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;
+  double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vect_norm, l1, l2, a_bs, b_bs, area_e, inv_area_e;
   
   double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store, *min_elevation_edgevalue, *max_elevation_edgevalue;
@@ -835 +847 @@
       cell_wetness_scale[k] = 0.;
       // Check if cell k is wet
-      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] + minimum_allowed_height+epsilon){
       //if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
           cell_wetness_scale[k] = 1.;  
@@ -918 +931 @@
           dyv1 = yv1 - y; 
           dyv2 = yv2 - y;
+
+
+          // 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;
+
+          //printf("slopes: %f, %f \n", a_bs, b_bs);
         }
 
@@ -995 +1022 @@
           area2 = dy2*dx1 - dy1*dx2; 
            
-          // Treat triangles with zero or 1 wet neighbours (area2 <=0.) 
+          //// Treat triangles with zero or 1 wet neighbours (area2 <=0.) 
           if ((area2 <= 0.))// || (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
           {
-            //if((cell_wetness_scale[k]==1.0) | (k==k0)&&(k==k1)&&(k==k2)){
-            //if(count_wet_neighbours[k] > 5.){
+
+              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;
+            
+                // 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
                 stage_edge_values[k3]   = stage_centroid_values[k];
                 stage_edge_values[k3+1] = stage_centroid_values[k];
                 stage_edge_values[k3+2] = stage_centroid_values[k];
+              }
             //}else{
             //}else{
@@ -1116 +1173 @@
           
           // Calculate heights of neighbouring cells
-          //hc = stage_centroid_values[k]  - elevation_centroid_values[k];
-          //h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
-          //h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
-          //h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
-          //hmin = min(min(h0, min(h1, h2)), hc);
+          hc = stage_centroid_values[k]  - elevation_centroid_values[k];
+          h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
+          h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
+          h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
+          hmin = min(min(h0, min(h1, h2)), hc);
           //// GD FIXME: No longer needed? 
           //hfactor = 0.0;
@@ -1130 +1187 @@
             //hfactor=hmin/(hmin + 0.004);
           //}
+          hfactor= min(2.0*max(hmin,0)/max(hc,1.0e-06), 1.0);
           
           //-----------------------------------
@@ -1143 +1201 @@
           dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
           dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
-          
+         
+          //if( (area2>0.0) ){ 
+
           inv_area2 = 1.0/area2;
-          // Calculate the gradient of stage on the auxiliary triangle
-          a = dy2*dq1 - dy1*dq2;
-          a *= inv_area2;
-          b = dx1*dq2 - dx2*dq1;
-          b *= inv_area2;
-          // d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth
-          //hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0);
-          //tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 
-          //           ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03;
-          //tmp /=max(hc*hc*hc, 1.0e-09);
-          //a = -xmom_centroid_values[k]*tmp;
-          //b = -ymom_centroid_values[k]*tmp;
+          //if(count_wet_neighbours[k] > 1){
+          //if(1==1){
+          //    // Calculate the gradient of stage on the auxiliary triangle
+              a = dy2*dq1 - dy1*dq2;
+              a *= inv_area2;
+              b = dx1*dq2 - dx2*dq1;
+              b *= inv_area2;
+          //}else{
+          ////    //// 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 /=max(pow(hc,10.0/3.0), 1.0e-30);
+          //    }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;
+          //}
           // Calculate provisional jumps in stage from the centroid 
           // of triangle k to its vertices, to be limited
@@ -1175 +1249 @@
           //if((count_wet_neighbours[k]>0)&&(cell_wetness_scale[k]=1.0)){//&&(k_wetdry==0)){
           //if(stage_centroid_values[k] > max_elevation_edgevalue[k]){
+          //if( (area2>0) ){ 
+          //if(count_wet_neighbours[k]>0){
               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];
@@ -1194 +1272 @@
           //}
           //}else{
-          // d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth
+          // 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);
           //    tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 
@@ -1265 +1343 @@
           //if((k!=k0)&(k!=k1)&(k!=k2))
           //if(stagemin>=bedmax){
-          if((count_wet_neighbours[k]>0)){ // && 
+          //if((count_wet_neighbours[k]>0)){ // && 
             // (cell_wetness_scale[k]==1.0)){
           //if(stage_centroid_values[k]>(minimum_allowed_height+max_elevation_edgevalue[k])){
             limit_gradient(dqv, qmin, qmax, beta_tmp);
-          }else{
-            dqv[0]=0.;
-            dqv[1]=0.;
-            dqv[2]=0.;
-          //  limit_gradient(dqv, qmin, qmax, 0.);
-          } 
+          //}else{
+          //  dqv[0]=0.;
+          //  dqv[1]=0.;
+          //  dqv[2]=0.;
+          ////  limit_gradient(dqv, qmin, qmax, 0.);
+          //} 
           //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
           
@@ -1318 +1396 @@
           // Limit the gradient
           //if(stagemin>=bedmax){
-          if((count_wet_neighbours[k]>0)){//&&
+          //if((count_wet_neighbours[k]>0)){//&&
              //(cell_wetness_scale[k]==1.0)){
             //(stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
           //if(stage_centroid_values[k]>(minimum_allowed_height+max_elevation_edgevalue[k])){
             limit_gradient(dqv, qmin, qmax, beta_tmp);
-          }else{
-            dqv[0]=0.;
-            dqv[1]=0.;
-            dqv[2]=0.;
-          } 
+          //}else{
+          //  dqv[0]=0.;
+          //  dqv[1]=0.;
+          //  dqv[2]=0.;
+          //} 
           
           for (i=0;i<3;i++)
@@ -1644 +1722 @@
 
 
-  PyArrayObject *neighbours, *neighbour_edges,
+  PyArrayObject *boundary_flux_sum, *neighbours, *neighbour_edges,
     *normals, *edgelengths, *radii, *areas,
     *tri_full_flag,
@@ -1670 +1748 @@
     
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOiOOOOO",
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOiOOOOO",
             &timestep,
             &epsilon,
             &H0,
             &g,
+            &boundary_flux_sum,
             &neighbours,
             &neighbour_edges,
@@ -1739 +1818 @@
                      H0,
                      g,
+                     (double*) boundary_flux_sum -> data,
                      (long*) neighbours -> data,
                      (long*) neighbour_edges -> data,

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_20130219_202649/dam_break.sww', 0.001)
+p_dev = util.get_output('dam_break_20130319_235611/dam_break.sww', 0.001)
 p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
 

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

@@ -38 +38 @@
 
 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)
+    #stg=-0.2 # Stage
     #topo=topography(x,y) #Bed
     return stg #*(stg>topo) + topo*(stg<=topo)
@@ -46 +46 @@
 domain.get_quantity('elevation').smooth_vertex_values() # Steve's fix -- without this, substantial artificial velcities are generated everywhere in the domain. With this fix, there are artificial velocities near the coast, but not elsewhere.
 
-domain.set_quantity('friction',0.10)             # Constant friction
+domain.set_quantity('friction',0.03)             # Constant friction
 
 domain.set_quantity('stage', stagefun)              # Constant negative initial stage

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

@@ -80 +80 @@
 for t in domain.evolve(yieldstep=0.200,finaltime=40.00):
     print domain.timestepping_statistics()
+    print domain.boundary_flux_integral
     xx = domain.quantities['xmomentum'].centroid_values
     yy = domain.quantities['ymomentum'].centroid_values
@@ -89 +90 @@
     print 'Peak velocity is: ', vv.max(), vv.argmax()
     print 'Volume is', sum(dd_raw*domain.areas)    
+    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

@@ -2 +2 @@
 Water flowing down a channel
 """
-import sys
 import numpy
 #---------------
@@ -21 +20 @@
 from anuga import rectangular_cross as rectangular_cross
 from anuga.structures.inlet_operator import Inlet_operator
-#from anuga.shallow_water.shallow_water_domain import Domain as Domain
+from anuga.shallow_water.shallow_water_domain import Domain as Domain
 from balanced_dev import *
 from balanced_dev import Domain as Domain
@@ -37 +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):
@@ -76 +75 @@
 for t in domain.evolve(yieldstep=4.0, finaltime=400.0):
     print domain.timestepping_statistics()
+    print domain.boundary_flux_integral
+    print (domain.areas*(domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values)).sum()
     #print domain.quantities['stage'].centroid_values[myindex] - domain.quantities['elevation'].centroid_values[myindex]
     s3 = domain.get_flow_through_cross_section([[30., 0.0], [30., 100.]])