Index: /trunk/anuga_core/source/anuga/operators/boundary_flux_integral_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/boundary_flux_integral_operator.py	(revision 9257)
+++ /trunk/anuga_core/source/anuga/operators/boundary_flux_integral_operator.py	(revision 9257)
@@ -0,0 +1,79 @@
+"""
+
+Integrate the fluxes through the boundaries
+
+Relies on boundary_flux_sum being computed in
+compute_fluxes.
+
+
+
+"""
+
+__author__="gareth"
+__date__ ="$10/07/2014 $"
+
+
+
+import numpy as num
+from anuga.operators.base_operator import Operator
+
+class boundary_flux_integral_operator(Operator):
+    """
+    Simple operator to collect the integral of the boundary fluxes during a run
+    """
+
+    def __init__(self,
+                 domain,
+                 description = None,
+                 label = None,
+                 logging = False,
+                 verbose = False):
+
+
+        Operator.__init__(self, domain, description, label, logging, verbose)
+
+        #------------------------------------------
+        # Setup a quantity to store the boundary flux integral 
+        #------------------------------------------
+        self.boundary_flux_integral=num.array([0.])
+
+        # Alias for domain
+        self.domain=domain
+        
+
+    def __call__(self):
+        """
+        Calculate max_stage at each timestep
+        """
+
+        dt=self.domain.timestep
+        ts_method=self.domain.timestepping_method        
+        
+        if(ts_method=='euler'): 
+            self.boundary_flux_integral = self.boundary_flux_integral + dt*self.domain.boundary_flux_sum[0]
+        elif(ts_method=='rk2'):
+            self.boundary_flux_integral = self.boundary_flux_integral + 0.5*dt*self.domain.boundary_flux_sum[0:2].sum()
+        else:
+            raise Exception, 'Cannot compute boundary flux integral with this timestepping method'
+     
+        # Zero the boundary_flux_sum 
+        self.domain.boundary_flux_sum=self.domain.boundary_flux_sum*0.
+
+    def parallel_safe(self):
+        """Operator is applied independently on each parallel domain
+
+        """
+        return True
+
+    def statistics(self):
+
+        message = self.label + ': Boundary_flux_integral operator'
+        return message
+
+
+    def timestepping_statistics(self):
+        from anuga import indent
+
+        message  = indent + self.label + ': Integrating the boundary flux'
+        return message
+
Index: /trunk/anuga_core/source/anuga/operators/test_boundary_flux_integral_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/test_boundary_flux_integral_operator.py	(revision 9257)
+++ /trunk/anuga_core/source/anuga/operators/test_boundary_flux_integral_operator.py	(revision 9257)
@@ -0,0 +1,86 @@
+import unittest
+import anuga
+import numpy
+
+boundaryPolygon=[ [0., 0.], [0., 100.], [100.0, 100.0], [100.0, 0.0]]
+
+verbose=False
+
+class Test_boundary_flux_integral_operator(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+    def create_domain(self, flowalg):
+        # Riverwall = list of lists, each with a set of x,y,z (and optional QFactor) values
+
+        # Make the domain
+        anuga.create_mesh_from_regions(boundaryPolygon, 
+                                 boundary_tags={'left': [0],
+                                                'top': [1],
+                                                'right': [2],
+                                                'bottom': [3]},
+                                   maximum_triangle_area = 200.,
+                                   minimum_triangle_angle = 28.0,
+                                   filename = 'test_boundaryfluxintegral.msh',
+                                   use_cache=False,
+                                   verbose=verbose)
+
+        domain=anuga.create_domain_from_file('test_boundaryfluxintegral.msh')
+
+        # 05/05/2014 -- riverwalls only work with DE0 and DE1
+        domain.set_flow_algorithm(flowalg)
+        domain.set_name('test_boundaryfluxintegral')
+
+        domain.set_store_vertices_uniquely()
+       
+        def topography(x,y):
+            return -x/150. 
+
+        # NOTE: Setting quantities at centroids is important for exactness of tests
+        domain.set_quantity('elevation',topography,location='centroids')     
+        domain.set_quantity('friction',0.03)             
+        domain.set_quantity('stage', topography,location='centroids')            
+       
+        # Boundary conditions
+        Br=anuga.Reflective_boundary(domain)
+        Bd=anuga.Dirichlet_boundary([0., 0., 0.])
+        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
+
+        return domain
+
+    def test_boundary_flux_operator_DE0(self):
+        """
+            
+        """
+        
+        domain=self.create_domain('DE0')
+        for t in domain.evolve(yieldstep=0.1,finaltime=1.0):
+            pass
+        # The domain was initially dry
+        vol=domain.get_water_volume()
+        boundaryFluxInt=domain.get_boundary_flux_integral()
+
+        assert(numpy.allclose(vol,boundaryFluxInt))
+    
+    def test_boundary_flux_operator_DE1(self):
+        """
+            
+        """
+        
+        domain=self.create_domain('DE1')
+        for t in domain.evolve(yieldstep=0.1,finaltime=1.0):
+            pass
+        # The domain was initially dry
+        vol=domain.get_water_volume()
+        boundaryFluxInt=domain.get_boundary_flux_integral()
+
+        assert(numpy.allclose(vol,boundaryFluxInt))
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_boundary_flux_integral_operator, 'test')
+    runner = unittest.TextTestRunner(verbosity=1)
+    runner.run(suite)
+
Index: /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 9256)
+++ /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 9257)
@@ -250,4 +250,13 @@
         self.riverwallData=anuga.structures.riverwall.RiverWall(self)
 
+        # Keep track of the fluxes through the boundaries
+        # Only works for DE algorithms at present
+        max_time_substeps=3 # Maximum number of substeps supported by any timestepping method
+        # boundary_flux_sum holds boundary fluxes on each sub-step
+        self.boundary_flux_sum=num.array([0.]*max_time_substeps) 
+        from anuga.operators.boundary_flux_integral_operator import boundary_flux_integral_operator
+        self.boundary_flux_integral=boundary_flux_integral_operator(self)
+        # Make an integer counting how many times we call compute_fluxes_central -- so we know which substep we are on
+        self.call=1 
 
 
@@ -480,12 +489,4 @@
         self.maximum_allowed_speed=0.0
 
-        ## FIXME: Should implement tracking of boundary fluxes
-        ## Keep track of the fluxes through the boundaries
-        self.boundary_flux_integral=num.ndarray(1)
-        self.boundary_flux_integral[0]=0. 
-        self.boundary_flux_sum=num.ndarray(1)
-        self.boundary_flux_sum[0]=0.
-
-        self.call=1 # Integer counting how many times we call compute_fluxes_central
 
         if self.processor == 0 and self.verbose:
@@ -547,13 +548,4 @@
         self.maximum_allowed_speed=0.0
 
-        ## FIXME: Should implement tracking of boundary fluxes
-        ## Keep track of the fluxes through the boundaries
-        self.boundary_flux_integral=num.ndarray(1)
-        self.boundary_flux_integral[0]=0. 
-        self.boundary_flux_sum=num.ndarray(1)
-        self.boundary_flux_sum[0]=0.
-
-        self.call=1 # Integer counting how many times we call compute_fluxes_central
-
         if self.processor == 0 and self.verbose:
             print '##########################################################################'
@@ -615,13 +607,4 @@
         self.maximum_allowed_speed=0.0
 
-        ## FIXME: Should implement tracking of boundary fluxes
-        ## Keep track of the fluxes through the boundaries
-        self.boundary_flux_integral=num.ndarray(1)
-        self.boundary_flux_integral[0]=0. 
-        self.boundary_flux_sum=num.ndarray(1)
-        self.boundary_flux_sum[0]=0.
-
-        self.call=1 # Integer counting how many times we call compute_fluxes_central
-
         if self.processor == 0 and self.verbose:
             print '##########################################################################'
@@ -682,13 +665,4 @@
 
         self.maximum_allowed_speed=0.0
-
-        ## FIXME: Should implement tracking of boundary fluxes
-        ## Keep track of the fluxes through the boundaries
-        self.boundary_flux_integral=num.ndarray(1)
-        self.boundary_flux_integral[0]=0. 
-        self.boundary_flux_sum=num.ndarray(1)
-        self.boundary_flux_sum[0]=0.
-
-        self.call=1 # Integer counting how many times we call compute_fluxes_central
 
         if self.processor == 0 and self.verbose:
@@ -1337,4 +1311,35 @@
         return water_volume
 
+    def get_boundary_flux_integral(self):
+        """
+            Compute the boundary flux integral.
+            Should work in parallel
+        """
+    
+        from anuga import myid, numprocs, send, receive, barrier
+
+        if not self.compute_fluxes_method=='DE':
+            msg='Boundary flux integral only supported for DE fluxes '+\
+                '(because computation of boundary_flux_sum is only implemented there)'
+            raise Exception, msg
+            
+        flux_integral = self.boundary_flux_integral.boundary_flux_integral
+            
+        if myid == 0:
+            for i in range(1,numprocs):
+                remote_flux_integral = receive(i)
+                flux_integral = flux_integral + remote_flux_integral
+        else:
+            send(flux_integral,0)
+        
+        barrier()
+    
+        if myid == 0:
+            for i in range(1,numprocs):
+                send(flux_integral,i)
+        else:
+            flux_integral = receive(0)
+        
+        return flux_integral 
 
     def get_flow_through_cross_section(self, polyline, verbose=False):
@@ -1587,29 +1592,4 @@
                                                self.riverwallData.hydraulic_properties)
 
-            ## FIXME: This won't work in parallel since the timestep has not been updated to include other processors. 
-            ## Update the boundary flux integral
-            ## FIXME GD: This should be moved to its own routine -- something like this
-            ##        could be included in all solvers.
-            #if(self.timestepping_method=='rk2'):
-            #    if(self.call%2==1):
-            #      self.boundary_flux_integral[0]= self.boundary_flux_integral[0] +\
-            #                                        self.boundary_flux_sum[0]*self.timestep*0.5
-            #      #print 'dbfi ', self.boundary_flux_integral, self.boundary_flux_sum
-            #      self.boundary_flux_sum[0]=0.
-            #elif(self.timestepping_method=='euler'):
-            #    # FIXME GD: Unsupported at present.
-            #    self.boundary_flux_integral=0.
-            #    # This presently doesn't work -- this section of code may need to go elsewhere
-            #    #self.boundary_flux_integral[0]= self.boundary_flux_integral[0] +\
-            #    #                                  self.boundary_flux_sum[0]*self.timestep
-            #    #self.boundary_flux_sum[0]=0.
-            #elif(self.timestepping_method=='rk3'):
-            #    # FIXME GD: Unsupported at present.
-            #    self.boundary_flux_integral=0.
-            #    # FIXME GD: This needs to be implemented
-            #else:
-            #    mess='ERROR: domain.timestepping_method', self.timestepping_method,' method not recognised'
-            #    raise Exception, mess
-
             self.flux_timestep = flux_timestep
 
@@ -2338,4 +2318,6 @@
         These calculations are only approximate since they don't use the
         flux calculation used in evolve
+
+        See get_boundary_flux_integral for an exact computation
         """
 		
Index: /trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c	(revision 9256)
+++ /trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c	(revision 9257)
@@ -415,4 +415,5 @@
       // 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==0) edgeflux[i] += (s_max*s_min)*(h_right - h_left);
       if(i==1) edgeflux[i] += (s_max*s_min)*(uh_right - uh_left);
       if(i==2) edgeflux[i] += (s_max*s_min)*(vh_right - vh_left);
@@ -584,5 +585,5 @@
     double bedslope_work, local_timestep;
     int neighbours_wet[3];//Work array
-    int RiverWall_count;
+    int RiverWall_count, substep_count;
     double hle, hre, zc, zc_n, Qfactor, s1, s2, h1, h2; 
     double stage_edge_lim, outgoing_mass_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2;
@@ -606,4 +607,5 @@
     local_timestep=timestep;
     RiverWall_count=0;
+    substep_count=(call-2)%timestep_fluxcalls;
 
 
@@ -790,5 +792,5 @@
             edgeflux_store[ki3 + 2 ] = -edgeflux[2];
 
-            // bedslope_work contains all gravity related terms -- weighting of
+            // bedslope_work contains all gravity related terms
             bedslope_work=length*(-g*0.5*(h_left*h_left - hle*hle -(hle+hc)*(zl-zc))+pressure_flux);
 
@@ -836,5 +838,5 @@
         } // End edge i (and neighbour n)
         // Keep track of maximal speeds
-        if((call-2)%timestep_fluxcalls==0) max_speed_array[k] = speed_max_last; //max_speed;
+        if(substep_count==0) max_speed_array[k] = speed_max_last; //max_speed;
 
 
@@ -843,4 +845,5 @@
     //// GD HACK 
     //// Limit edgefluxes, for mass conservation near wet/dry cells
+    //// This doesn't seem to be needed anymore
     //for(k=0; k< number_of_elements; k++){
     //    //continue;
@@ -925,7 +928,8 @@
             // boundary_flux_integral
             if( (n<0 & tri_full_flag[k]==1) | ( n>=0 && (tri_full_flag[k]==1 & tri_full_flag[n]==0)) ){ 
-                boundary_flux_sum[0] += edgeflux_store[ki3];
+                // boundary_flux_sum is an array with length = timestep_fluxcalls
+                // For each sub-step, we put the boundary flux sum in.
+                boundary_flux_sum[substep_count] += edgeflux_store[ki3];
             }
-
             // GD HACK
             // Compute bed slope term
@@ -950,6 +954,6 @@
 
     // Hack to ensure we only update the timestep on the first call within each rk2/rk3 step
-    if((call-2)%timestep_fluxcalls==0) timestep=local_timestep; 
-
+    if(substep_count==0) timestep=local_timestep; 
+            
     free(edgeflux_store);
     free(pressuregrad_store);
Index: /trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py
===================================================================
--- /trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 9256)
+++ /trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py	(revision 9257)
@@ -9,13 +9,6 @@
 
 from math import sin, pi, exp
-#from anuga.shallow_water.shallow_water_domain import Domain as Domain
-#from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
-#path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
-#from swb2_domain import *
-#from balanced_basic import *
-#from balanced_dev import *
 from anuga import Domain as Domain
-#from bal_and import *
-#from anuga_tsunami import *
+from anuga.operators.boundary_flux_integral_operator import boundary_flux_integral_operator
 
 #---------
@@ -27,11 +20,7 @@
 domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
 #domain.set_store_centroids(True)
-#domain.set_timestepping_method('euler')
-#domain.set_flow_algorithm('DE1')
-#domain.set_flow_algorithm('tsunami')
 domain.set_flow_algorithm('DE0')
 domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 
          'ymomentum': 2, 'elevation': 2})
-#domain.minimum_storable_height=1.0e-03
 domain.set_store_vertices_uniquely()
 #------------------
@@ -70,9 +59,4 @@
 
 print '#-#', domain.minimum_allowed_height
-#domain.minimum_allowed_height=0.001
-#def frict_change(x,y):
-#	return 0.2*(x>0.5)+0.1*(x<=0.5)
-#
-#domain.set_quantity('friction',frict_change)
 
 domain.set_quantity('stage', stagefun)              # Constant negative initial stage
@@ -108,5 +92,4 @@
 #Evolve the system through time
 #------------------------------
-
 for t in domain.evolve(yieldstep=tstep,finaltime=lasttime):
 #for t in domain.evolve(yieldstep=0.2,finaltime=60.):
@@ -121,7 +104,11 @@
     vv = vv*(dd>1.0e-03)
     print 'Peak velocity is: ', vv.max(), vv.argmax(), dd[vv.argmax()]
-    print 'Volume is', sum(dd_raw*domain.areas)    
+    vol=sum(dd_raw*domain.areas)
+    print 'Volume is', vol    
+    flux_integral=domain.get_boundary_flux_integral()
+    print 'Flux integral is: ', flux_integral
+    print 'Their difference is:', vol-flux_integral
     #print 'Volume less flux int', sum(dd_raw*domain.areas) - domain.boundary_flux_integral
-
+    
 
 vel_final_inds=(vv>1.0e-01).nonzero()
