Changeset 9059

Timestamp:

Mar 9, 2014, 9:11:36 PM (11 years ago)

Author:

steve

Message:

changed the order of calls in anuga.init.py to hopefully call the correct Boyd_box_operator in parallel

Location:

trunk/anuga_core/source/anuga

Files:

• __init__.py (modified) (5 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (3 diffs)
• abstract_2d_finite_volumes/test_quantity.py (modified) (1 diff)
• abstract_2d_finite_volumes/test_region.py (modified) (1 diff)
• geospatial_data/geospatial_data.py (modified) (1 diff)
• operators/erosion_operators.py (modified) (4 diffs)
• shallow_water/shallow_water_domain.py (modified) (5 diffs)
• shallow_water/swDE1_domain_ext.c (modified) (6 diffs)
• structures/test_boyd_pipe_operator.py (modified) (1 diff)

trunk/anuga_core/source/anuga/__init__.py

@@ -31 +31 @@
 from anuga.__metadata__ import __version__, __date__, __author__
 
+#--------------------------------
+# Important basic classes
+#--------------------------------
 from anuga.shallow_water.shallow_water_domain import Domain
 from anuga.abstract_2d_finite_volumes.quantity import Quantity
 from anuga.abstract_2d_finite_volumes.region import Region
-
+from anuga.operators.base_operator import Operator
+from anuga.structures.structure_operator import Structure_operator
 
 from anuga.abstract_2d_finite_volumes.util import file_function, \
@@ -63 +67 @@
 from anuga.geometry.polygon import read_polygon
 from anuga.caching import cache
+from os.path import join
+from anuga.config import indent
+
+
+
+#----------------------------
+# Parallel api (needs operators)
+#----------------------------
+from anuga_parallel.parallel_api import distribute
+from anuga_parallel.parallel_api import myid, numprocs, get_processor_name
+from anuga_parallel.parallel_api import send, receive
+from anuga_parallel.parallel_api import pypar_available, barrier, finalize
+
 
 #-----------------------------
@@ -135 +152 @@
 from anuga.shallow_water.sww_interrogate import get_flow_through_cross_section
     
-
 #---------------------------
 # Operators
 #---------------------------
-from anuga.operators.base_operator import Operator
 from anuga.operators.kinematic_viscosity_operator import Kinematic_viscosity_operator
 
@@ -145 +160 @@
 from anuga.operators.set_friction_operators import Depth_friction_operator 
 
+
 #---------------------------
 # Structure Operators
 #---------------------------
-from anuga.structures.structure_operator import Structure_operator
-from anuga.structures.boyd_box_operator import Boyd_box_operator
-from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
-from anuga.structures.weir_orifice_trapezoid_operator import Weir_orifice_trapezoid_operator
-from anuga.structures.inlet_operator import Inlet_operator
+
+
+if pypar_available:
+    from anuga_parallel.parallel_operator_factory import Inlet_operator
+    from anuga_parallel.parallel_operator_factory import Boyd_box_operator
+    from anuga_parallel.parallel_operator_factory import Boyd_pipe_operator
+    from anuga_parallel.parallel_operator_factory import Weir_orifice_trapezoid_operator
+else:
+    from anuga.structures.boyd_box_operator import Boyd_box_operator
+    from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
+    from anuga.structures.weir_orifice_trapezoid_operator import Weir_orifice_trapezoid_operator
+    from anuga.structures.inlet_operator import Inlet_operator
+
+#----------------------------
+# Parallel distribute
+#----------------------------
+
+
+#----------------------------
+# 
+#Added by Petar Milevski 10/09/2013
+#import time, os
+
+from anuga.utilities.model_tools import get_polygon_from_single_file
+from anuga.utilities.model_tools import get_polygons_from_Mid_Mif
+from anuga.utilities.model_tools import get_polygon_list_from_files
+from anuga.utilities.model_tools import get_polygon_dictionary
+from anuga.utilities.model_tools import get_polygon_value_list
+from anuga.utilities.model_tools import read_polygon_dir
+from anuga.utilities.model_tools import read_hole_dir_multi_files_with_single_poly
+from anuga.utilities.model_tools import read_multi_poly_file
+from anuga.utilities.model_tools import read_hole_dir_single_file_with_multi_poly
+from anuga.utilities.model_tools import read_multi_poly_file_value
+from anuga.utilities.model_tools import Create_culvert_bridge_Operator
+
 
 #---------------------------
@@ -374 +420 @@
 
 
-#Added by Petar Milevski 10/09/2013
-import time, os
-from os.path import join
-from anuga.config import indent
-
-from anuga.utilities.model_tools import get_polygon_from_single_file
-from anuga.utilities.model_tools import get_polygons_from_Mid_Mif
-from anuga.utilities.model_tools import get_polygon_list_from_files
-from anuga.utilities.model_tools import get_polygon_dictionary
-from anuga.utilities.model_tools import get_polygon_value_list
-from anuga.utilities.model_tools import read_polygon_dir
-from anuga.utilities.model_tools import read_hole_dir_multi_files_with_single_poly
-from anuga.utilities.model_tools import read_multi_poly_file
-from anuga.utilities.model_tools import read_hole_dir_single_file_with_multi_poly
-from anuga.utilities.model_tools import read_multi_poly_file_value
-
-from anuga.utilities.model_tools import Create_culvert_bridge_Operator
-
-from anuga_parallel.parallel_api import distribute
-from anuga_parallel.parallel_api import myid, numprocs, get_processor_name
-from anuga_parallel.parallel_api import send, receive
-from anuga_parallel.parallel_api import pypar_available, barrier, finalize
-
-if pypar_available:
-    #from anuga_parallel.parallel_meshes import parallel_rectangle
-    #from anuga_parallel.parallel_shallow_water import Parallel_domain as Parallel_shallow_water_domain
-    #from anuga_parallel.parallel_advection     import Parallel_domain as Parallel_advection_domain
-    from anuga_parallel.parallel_operator_factory import Inlet_operator, Boyd_box_operator, Boyd_pipe_operator, Weir_orifice_trapezoid_operator
-
-
-from anuga.structures.weir_orifice_trapezoid_operator import Weir_orifice_trapezoid_operator
-
+
+

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -335 +335 @@
 
 
+    def save_data_to_dem(self,filename=None):
+
+
+        #FIXME SR: Should add code to deal with parallel
+
+        ids = self.domain.tri_full_flag == 1
+        c_v = self.centroid_values[ids].reshape((-1,1))
+        c_x = self.domain.centroid_coordinates[ids,0].reshape((-1,1))
+        c_y = self.domain.centroid_coordinates[ids,1].reshape((-1,1))
+
+        import numpy
+
+        c_xyv = numpy.hstack((c_x, c_y, c_v))
+
+        if filename is None:
+            filename= self.name
+
+        if self.domain.parallel:
+            fullfilename = filename+'_centroid_data_P%g_%g.csv'% \
+                    ( self.domain.numproc, self.domain.processor)
+        else:
+            fullfilename = filename+'_centroid_data.csv'
+
+        numpy.savetxt(fullfilename, c_xyv, delimiter=',', fmt =  ['%.15e', '%.15e', '%.15e' ])
+
+
+        if self.domain.parallel:
+            import pypar
+            pypar.barrier()
+
+            # On processor 0 catenate the files
+            if self.domain.processor == 0:
+                import shutil
+                import os
+                destination = open(filename+'_centroid_data.csv','wb')
+                np = self.domain.numproc
+                files = [ filename+'_centroid_data'+"_P"+str(np)+"_"+str(v)+".csv" for v in range(np)]
+                for file in files:
+                    shutil.copyfileobj(open(file,'rb'), destination)
+                destination.close()
+                for file in files:
+                    os.remove(file)
+
 
 
@@ -357 +400 @@
             plt.show()
         
+
+    def save_to_array(self,
+                cellsize=10,
+                NODATA_value=-9999.0,
+                easting_min=None,
+                easting_max=None,
+                northing_min=None,
+                northing_max=None,
+                origin=None,
+                verbose=False):
+        
+        """Interpolate quantity to an array
+        """
+    
+
+        verbose = True
+        
+        from anuga.geometry.polygon import inside_polygon, outside_polygon
+        from anuga.abstract_2d_finite_volumes.util import \
+             apply_expression_to_dictionary
+    
+    
+    
+        #Get extent and reference
+        
+        domain = self.domain
+
+        volumes = domain.triangles
+        
+
+        x,y,_,v= self.get_vertex_values(xy=True)
+        
+        
+        vertex_coordinates = domain.vertex_coordinates
+
+        # store the connectivity data
+        #points = num.concatenate((x[:,num.newaxis],y[:,num.newaxis]), axis=1)
+#         self.writer.store_triangulation(fid,
+#                                         points,
+#                                         V.astype(num.float32),
+#                                         points_georeference=\
+#                                         domain.geo_reference)
+        
+        false_easting = 500000
+        false_northing = 10000000
+    
+    
+    
+        geo_ref = self.domain.geo_reference
+        
+        xllcorner = geo_ref.get_xllcorner()
+        yllcorner = geo_ref.get_yllcorner()
+        
+        if verbose: 
+            print 
+            print xllcorner
+            print yllcorner
+            print x
+            print y
+            print vertex_coordinates[:,0]
+            print vertex_coordinates[:,1]
+            
+        
+        
+        # Create grid and update xll/yll corner and x,y
+        # Relative extent
+        if easting_min is None:
+            xmin = min(x)
+        else:
+            xmin = easting_min - xllcorner
+    
+        if easting_max is None:
+            xmax = max(x)
+        else:
+            xmax = easting_max - xllcorner
+    
+        if northing_min is None:
+            ymin = min(y)
+        else:
+            ymin = northing_min - yllcorner
+    
+        if northing_max is None:
+            ymax = max(y)
+        else:
+            ymax = northing_max - yllcorner
+    
+    
+        """
+        msg = 'xmax must be greater than or equal to xmin.\n'
+        msg += 'I got xmin = %f, xmax = %f' %(xmin, xmax)
+        assert xmax >= xmin, msg
+    
+        msg = 'ymax must be greater than or equal to xmin.\n'
+        msg += 'I got ymin = %f, ymax = %f' %(ymin, ymax)
+        assert ymax >= ymin, msg
+    
+        if verbose: log.critical('Creating grid')
+        ncols = int((xmax-xmin)/cellsize) + 1
+        nrows = int((ymax-ymin)/cellsize) + 1
+    
+        # New absolute reference and coordinates
+        newxllcorner = xmin + xllcorner
+        newyllcorner = ymin + yllcorner
+    
+        x = x + xllcorner - newxllcorner
+        y = y + yllcorner - newyllcorner
+    
+    
+        grid_values = num.zeros( (nrows*ncols, ), num.float)
+        #print '---',grid_values.shape
+    
+        num_tri =  len(volumes)
+        norms = num.zeros(6*num_tri, num.float)
+    
+        #Use fasr method to calc grid values
+        from calc_grid_values_ext import calc_grid_values
+    
+        calc_grid_values(nrows, ncols, cellsize, NODATA_value,
+                         x,y, norms, volumes, result, grid_values)
+    
+    
+        fid.close()
+        
+        #print outside_indices
+    
+        if verbose:
+            log.critical('Interpolated values are in [%f, %f]'
+                         % (num.min(grid_values), num.max(grid_values)))
+    
+    
+    
+        return x,y, grid_values.reshape(nrows,ncols)[::-1,:]
+
+        """
+
+
+
+
+
+
 
 
@@ -1121 +1304 @@
         Parameters
         """
-
-
-
-
-#         msg = 'Filename must be a text string'
-#         assert isinstance(filename, basestring), msg
-#         
-#         msg = 'Extension should be .grd or asc'
-#         assert os.path.splitext(filename)[1] in ['.grd', '.asc'], msg
-#         
-# 
-#         msg = "set_values_from_grd_file is only defined for location='vertices' or 'centroids'"
-#         assert location in ['vertices', 'centroids'], msg
-# 
-#             
-# 
-#             
-#     
-#         root = filename[:-4]
-#     
-#     
-#         #Read DEM data
-#         datafile = open(filename)
-#     
-#         if verbose: log.critical('Reading data from %s' % (filename))
-#     
-#         lines = datafile.readlines()
-#         datafile.close()
-#     
-#         if verbose: log.critical('Got %d lines' % len(lines))
-#     
-# 
-#         ncols = int(lines[0].split()[1].strip())
-#         nrows = int(lines[1].split()[1].strip())
-# 
-#     
-#         # Do cellsize (line 4) before line 2 and 3
-#         cellsize = float(lines[4].split()[1].strip())
-#     
-#         # Checks suggested by Joaquim Luis
-#         # Our internal representation of xllcorner
-#         # and yllcorner is non-standard.
-#         xref = lines[2].split()
-#         if xref[0].strip() == 'xllcorner':
-#             xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
-#         elif xref[0].strip() == 'xllcenter':
-#             xllcorner = float(xref[1].strip())
-#         else:
-#             msg = 'Unknown keyword: %s' % xref[0].strip()
-#             raise Exception, msg
-#     
-#         yref = lines[3].split()
-#         if yref[0].strip() == 'yllcorner':
-#             yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
-#         elif yref[0].strip() == 'yllcenter':
-#             yllcorner = float(yref[1].strip())
-#         else:
-#             msg = 'Unknown keyword: %s' % yref[0].strip()
-#             raise Exception, msg
-#     
-#         NODATA_value = int(float(lines[5].split()[1].strip()))
-#     
-#         assert len(lines) == nrows + 6
-#     
-#   
-#         #Store data
-#         import numpy
-#     
-#         datafile = open(filename)
-#         Z = numpy.loadtxt(datafile, skiprows=6)
-#         datafile.close()
-#         
-#         #print Z.shape
-#         #print Z
-#         
-#         # For raster data we need to a flip and transpose
-#         Z = numpy.flipud(Z)
-# 
-#         # Transpose z to have y coordinates along the first axis and x coordinates
-#         # along the second axis
-#         Z = Z.transpose()
-#     
-#         x = num.linspace(xllcorner, xllcorner+cellsize*(ncols-1), ncols)
-#         y = num.linspace(yllcorner, yllcorner+cellsize*(nrows-1), nrows)
         
         

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -247 +247 @@
                                                    [2.2, 2.8, 4.0],
                                                    [2.5, -0.5, -2.0]])
-
-
-        #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
-        #                                       [5., 5., 5.],
-        #                                       [4.5, 4.5, 0.],
-        #                                       [3.0, -1.5, -1.5]])
+        
+    def test_save_to_array(self):
+        quantity = Quantity(self.mesh4,
+                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+        
+        
+        assert num.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
+
+
+        quantity.save_to_array()
+        
+        
+
+
+
 
     def test_get_extrema_1(self):

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py

@@ -67 +67 @@
         domain = Domain(points, vertices, boundary)
         
-        region = Region(domain, polygon=[[0.0,0.0], [0.5,0.0], [0.5,0.5]])
+        poly = [[0.0,0.0], [0.5,0.0], [0.5,0.5]]
+        
+        print poly
+        region = Region(domain, polygon=poly)
         
         expected_indices = [1]

trunk/anuga_core/source/anuga/geospatial_data/geospatial_data.py

@@ -701 +701 @@
             self.verbose_block_size = (self.last_row + 10)/10
             self.block_number = 0
-            self.number_of_blocks = self.number_of_points/self.max_read_lines
+            self.number_of_blocks = int(self.number_of_points/self.max_read_lines)
             # This computes the number of full blocks. The last block may be
             # smaller and won't be included in this estimate.
 
             if self.verbose is True:
-                log.critical('Geospatial_data: Reading %d points (in ~%d blocks) from file %s. '
+                log.critical('Geospatial_data: Reading %d points (in %d block(s)) from file %s. '
                              % (self.number_of_points, self.number_of_blocks+1,
                                 self.file_name))

trunk/anuga_core/source/anuga/operators/erosion_operators.py

@@ -104 +104 @@
             # Update all three vertices for each cell
             #--------------------------------------
-            self.elev_v[:] = self.elev_v + 0.0
+            self.elev_c[:] = self.elev_c + 0.0
 
         else:
@@ -114 +114 @@
             ind = self.indices
             m = num.sqrt(self.xmom_c[ind]**2 + self.ymom_c[ind]**2)
-            m = num.vstack((m,m,m)).T  # Stack up m to apply to vertices
-            m = num.where(m>self.threshold, m, 0.0)
-
-            de = m*dt
-            self.elev_v[ind] = num.maximum(self.elev_v[ind] - de, self.base)
+            
+            if self.domain.flow_algorithm == 'DE1':
+                m = num.where(m>self.threshold, m, 0.0)
+    
+                de = m*dt
+                height = self.stage_c[ind] - self.elev_c[ind]
+                self.elev_c[ind] = num.maximum(self.elev_v[ind] - de, self.base)
+                self.stage_c[ind] = self.elev_c[ind] + height
+            else:
+                m = num.vstack((m,m,m)).T  # Stack up m to apply to vertices
+                m = num.where(m>self.threshold, m, 0.0)
+    
+                de = m*dt
+                self.elev_v[ind] = num.maximum(self.elev_v[ind] - de, self.base)
 
             self.max_change = num.max(de)
@@ -141 +150 @@
 
         #------------------------------------------
-        # Apply changes to elevation vertex values
+        # Apply changes to elevation values
         # via the update_quantites routine
         #------------------------------------------
@@ -151 +160 @@
         # Cleanup elevation and stage quantity values
         #-----------------------------------------
+        self.clean_up_elevation_stage()
+        
+        
+        
+        
+    def clean_up_elevation_stage(self):
+        
+        #----------------------------------------------
+        # Don't need to clean up if using discontinuous
+        # elevation
+        #----------------------------------------------
+        if self.domain.flow_algorithm == 'DE1':
+            return 
+        
+        
+        #-----------------------------------------------
+        # Clean up to conserve volume and make elevation 
+        # continuous
+        #-----------------------------------------------
         if self.indices is None:
 
-            #--------------------------------------
-            # Update all three vertices for each cell
-            #--------------------------------------
-            self.elev_v[:] = self.elev_v + 0.0
 
             #--------------------------------------

trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -369 +369 @@
         #self.forcing_terms.append(manning_friction_explicit)
         #self.forcing_terms.remove(manning_friction_implicit)
-
-#        print '##########################################################################'
-#        print '#'
-#        print "#"
-#        print "# Here are some tips on using the 'tsunami' solver"
-#        print "#"
-#        print "# 1) When plotting outputs, I strongly suggest you examine centroid values, not vertex values"
-#        print "# , as the latter can be completely misleading near strong gradients in the flow. "
-#        print "# There is a plot_util.py script in anuga_core/utilities/ which might help you extract"
-#        print "# quantities at centroid values from sww files."
-#        print "# Note that to accuractely compute centroid values from sww files, the files need to store "
-#        print "# vertices uniquely. This makes for large sww files (3x), but is the price to pay for the right answer"
-#        print "# (unless we alter IO to allow centroids to be written to sww files, which would then affect"
-#        print "# ANUGA viewer as well -- I expect this would be lots of work)"
-#        print "#"
-#        print "# 2) In field scale applications (where the depth is typically > 1m), I suggest you set"
-#        print "# domain.minimum_allowed_height=0.01 (the default is 1.0e-3). "
-#        print "#"
-#        print "# 3) This solver is not expected to perform well in problems with very"
-#        print "# shallow water flowing down steep slopes (such that the stage_centroid_value "
-#        print "# is less than the maximum bed_edge_value on a given triangle). However, analytical tests"
-#        print "# suggest it can do typical wetting/drying situations very well (parabolic oscillations test case) "
-#        print "#"
-#        print "# 4) This solver allows the stage_centroid_value to drop to slightly below the minimum bed_vertex_value"
-#        print "# on it's triangle. In other ANUGA versions (e.g. 1.2.1), the limit would be the"
-#        print "# bed_centroid_value. This means that triangles store slightly more water than they are"
-#        print "# typically interpreted to store, which might have significance in some applications."
-#        print "#"
-#        print "# You will probably be able to tell this is causing you problems by convergence testing"
-#        print "#"
-#        print '# 5) Note that many options in config.py have been overridden by the solver -- we have '
-#        print '# deliberately attempted to get the solver to perform well with consistent values of '
-#        print '# these parameters -- so I advise against changing them unless you at least check that '
-#        print '# it really does improve things'
-#        print '#'
-#        print '##########################################################################'
+        if self.processor == 0 and self.verbose:
+            print '##########################################################################'
+            print '#'
+            print "#"
+            print "# Here are some tips on using the 'tsunami' solver"
+            print "#"
+            print "# 1) When plotting outputs, I strongly suggest you examine centroid values, not vertex values"
+            print "# , as the latter can be completely misleading near strong gradients in the flow. "
+            print "# There is a plot_util.py script in anuga_core/utilities/ which might help you extract"
+            print "# quantities at centroid values from sww files."
+            print "# Note that to accuractely compute centroid values from sww files, the files need to store "
+            print "# vertices uniquely. This makes for large sww files (3x), but is the price to pay for the right answer"
+            print "# (unless we alter IO to allow centroids to be written to sww files, which would then affect"
+            print "# ANUGA viewer as well -- I expect this would be lots of work)"
+            print "#"
+            print "# 2) In field scale applications (where the depth is typically > 1m), I suggest you set"
+            print "# domain.minimum_allowed_height=0.01 (the default is 1.0e-3). "
+            print "#"
+            print "# 3) This solver is not expected to perform well in problems with very"
+            print "# shallow water flowing down steep slopes (such that the stage_centroid_value "
+            print "# is less than the maximum bed_edge_value on a given triangle). However, analytical tests"
+            print "# suggest it can do typical wetting/drying situations very well (parabolic oscillations test case) "
+            print "#"
+            print "# 4) This solver allows the stage_centroid_value to drop to slightly below the minimum bed_vertex_value"
+            print "# on it's triangle. In other ANUGA versions (e.g. 1.2.1), the limit would be the"
+            print "# bed_centroid_value. This means that triangles store slightly more water than they are"
+            print "# typically interpreted to store, which might have significance in some applications."
+            print "#"
+            print "# You will probably be able to tell this is causing you problems by convergence testing"
+            print "#"
+            print '# 5) Note that many options in config.py have been overridden by the solver -- we have '
+            print '# deliberately attempted to get the solver to perform well with consistent values of '
+            print '# these parameters -- so I advise against changing them unless you at least check that '
+            print '# it really does improve things'
+            print '#'
+            print '##########################################################################'
 
     def set_DE1_defaults(self):
@@ -461 +461 @@
         self.call=1 # Integer counting how many times we call compute_fluxes_central
 
-        print '##########################################################################'
-        print '#'
-        print '# Using discontinuous elevation solver DE1 '
-        print '#'
-        print '# Mostly designed for rk2 timestepping'
-        print '#'
-        print '# Make sure you use centroid values when reporting on important output quantities'
-        print '#'
-        print '# NOTE: anuga-viewer sometimes shows very shallow regions in this solver as '
-        print '# erratically jumping from "totally dry" to "barely wet". In all my checks of centroid'
-        print '# values at such locations, I never found that it was really occurring. It may be'
-        print '# a problem with the viewer rather than this algorithm. '
-        print '##########################################################################'
+        if self.processor == 0 and self.verbose:
+            print '##########################################################################'
+            print '#'
+            print '# Using discontinuous elevation solver DE1 '
+            print '#'
+            print '# Mostly designed for rk2 timestepping'
+            print '#'
+            print '# Make sure you use centroid values when reporting on important output quantities'
+            print '#'
+            print '# NOTE: anuga-viewer sometimes shows very shallow regions in this solver as '
+            print '# erratically jumping from "totally dry" to "barely wet". In all my checks of centroid'
+            print '# values at such locations, I never found that it was really occurring. It may be'
+            print '# a problem with the viewer rather than this algorithm. '
+            print '##########################################################################'
 
 
@@ -615 +616 @@
            2.5
            tsunami
+           DE1
         """
 
@@ -1549 +1551 @@
                 self.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas)
 
-            print 'Cumulative mass protection: python'+str(mass_error)+'m^3 '
+            if mass_error > 0.0 and self.verbose :
+                print 'Cumulative mass protection: '+str(mass_error)+' m^3 '
 
         elif self.flow_algorithm == 'DE1':
@@ -1566 +1569 @@
             yc = self.centroid_coordinates[:,1] 
 
-            protect(self.minimum_allowed_height, self.maximum_allowed_speed,
+            mass_error = protect(self.minimum_allowed_height, self.maximum_allowed_speed,
                     self.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas, xc, yc)
+            
+            if mass_error > 0.0 and self.verbose :
+                print 'Cumulative mass protection: '+str(mass_error)+' m^3 '
             
         else:

trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c

@@ -596 +596 @@
 
 // Protect against the water elevation falling below the triangle bed
-int _protect(int N,
+double  _protect(int N,
          double minimum_allowed_height,
          double maximum_allowed_speed,
@@ -613 +613 @@
   double hc, bmin, bmax;
   double u, v, reduced_speed;
-  static double mass_error=0.; 
+  double mass_error=0.;
   // This acts like minimum_allowed height, but scales with the vertical
   // distance between the bed_centroid_value and the max bed_edge_value of
@@ -651 +651 @@
     }
 
-  if(mass_added==1){
-     printf("Cumulative mass protection: %f m^3 \n", mass_error);
-  }
-  return 0;
+  //if(mass_added==1){
+  //  printf("Cumulative mass protection: %f m^3 \n", mass_error);
+  //}
+  return mass_error;
 }
 
@@ -2031 +2031 @@
 
   int N;
+  double mass_error;
   double minimum_allowed_height, maximum_allowed_speed, epsilon;
 
@@ -2045 +2046 @@
   N = wc -> dimensions[0];
 
-  _protect(N,
+  mass_error = _protect(N,
        minimum_allowed_height,
        maximum_allowed_speed,
@@ -2059 +2060 @@
        (double*) yc -> data );
 
-  return Py_BuildValue("");
+  return Py_BuildValue("d", mass_error);
 }
 

trunk/anuga_core/source/anuga/structures/test_boyd_pipe_operator.py

@@ -12 +12 @@
 import numpy
 
-verbose = True
+verbose = False
 #diameter = width