Changeset 7562

Timestamp:

Nov 19, 2009, 5:23:52 PM (15 years ago)

Author:

steve

Message:

Updating the balanced and parallel code

Location:

anuga_core/source

Files:

• anuga/abstract_2d_finite_volumes/domain.py (modified) (4 diffs)
• anuga/abstract_2d_finite_volumes/test_domain.py (modified) (2 diffs)
• anuga/abstract_2d_finite_volumes/test_util.py (modified) (1 diff)
• anuga/culvert_flows/test_culvert_class.py (modified) (3 diffs)
• anuga/shallow_water/data_manager.py (modified) (4 diffs)
• anuga/shallow_water/shallow_water_domain.py (modified) (3 diffs)
• anuga/shallow_water/test_read_sww.py (added)
• anuga/shallow_water/test_shallow_water_domain.py (modified) (7 diffs)
• anuga/shallow_water_balanced/test_swb_distribute.py (modified) (6 diffs)
• anuga/shallow_water_balanced/test_swb_forcing_terms.py (modified) (2 diffs)
• anuga/test_all.py (modified) (1 diff)
• anuga_parallel/interface.py (modified) (1 diff)
• anuga_parallel/test_parallel_distribute_domain.py (modified) (2 diffs)
• anuga_parallel/test_parallel_distribute_mesh.py (modified) (2 diffs)
• anuga_parallel/test_parallel_sw_flow.py (modified) (3 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -1138 +1138 @@
         # Input checks
         if quantities is None:
-            quantities = self.conserved_quantities
+            quantities = self.evolved_quantities
         elif type(quantities) == types.StringType:
             quantities = [quantities] #Turn it into a list
@@ -1752 +1752 @@
     ##
     # @brief Mapping between conserved quantites and evolved quantities
-    # @param q_cons array of conserved quantity values
-    # @param q_evolved array of current evolved quantity values
-    def  conserved_to_evolved(self, q_cons, q_evolved):
+    # @param Input: q_cons array of conserved quantity values
+    # @param Input: q_evol array of current evolved quantity values
+    # @note  Output: Updated q_evol array
+    def  conserved_values_to_evolved_values(self, q_cons, q_evol):
         """Needs to be overridden by Domain subclass
         """
 
-        if len(q_cons) == len(q_evolved):
-            q_evolved[:] = q_cons
+        if len(q_cons) == len(q_evol):
+            q_evol[:] = q_cons
         else:
-            msg = 'Method conserved_to_evolved must be overridden by Domain subclass'
+            msg = 'Method conserved_values_to_evolved_values must be overridden by Domain subclass'
             raise Exception, msg
+
+        return q_evol
     
     ##
@@ -1781 +1784 @@
                 log.critical('WARNING: Ignored boundary segment (None)')
             else:
-                q_cons = B.evaluate(vol_id, edge_id)
-
-                if len(q_cons) == len(self.evolved_quantities):
+                q_bdry = B.evaluate(vol_id, edge_id)
+
+                if len(q_bdry) == len(self.evolved_quantities):
                     # conserved and evolved quantities are the same
-                    q_evol = q_cons
-                elif len(q_cons) == len(self.conserved_quantities):
+                    q_evol = q_bdry
+                elif len(q_bdry) == len(self.conserved_quantities):
                     # boundary just returns conserved quantities
                     # Need to calculate all the evolved quantities
@@ -1793 +1796 @@
                     q_evol = self.get_evolved_quantities(vol_id, edge = edge_id)
 
-                    self.conserved_to_evolved(q_cons, q_evol)
+                    q_evol = self.conserved_values_to_evolved_values(q_bdry, q_evol)
                 else:
                     msg = 'Boundary must return array of either conserved or evolved quantities'

anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

@@ -606 +606 @@
             pass
         else:
-            msg = 'Should have caught the lack of conserved_to_evolved member function'
+            msg = 'Should have caught the lack of conserved_values_to_evolved_values member function'
             raise Exception, msg
 
 
-        def  conserved_to_evolved(q_cons, q_evol):
+        def  conserved_values_to_evolved_values(q_cons, q_evol):
 
             q_evol[0:3] = q_cons
@@ -616 +616 @@
             q_evol[4] = q_cons[2]/q_cons[0]
 
-        domain.conserved_to_evolved = conserved_to_evolved
+            return q_evol
+
+        domain.conserved_values_to_evolved_values = conserved_values_to_evolved_values
 
         domain.update_boundary()

anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

@@ -140 +140 @@
         domain1.store = True
         domain1.set_datadir('.')
-        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
+        sww_file = 'spatio_temporal_boundary_source_%d' %(id(self))
+        domain1.set_name(sww_file)
 
         #Bed-slope, friction and IC at vertices (and interpolated elsewhere)

anuga_core/source/anuga/culvert_flows/test_culvert_class.py

@@ -124 +124 @@
         
 
+    os.remove('Test_culvert_shallow.sww')
 
 class Test_Culvert(unittest.TestCase):
@@ -243 +244 @@
         assert min_delta_w < 0
         assert max_delta_w > 10        
-        
+
+
+        os.remove('Test_culvert.sww')
 
     def test_that_culvert_dry_bed_rating_does_not_produce_flow(self):
@@ -346 +349 @@
 
     
-
+        os.remove('Test_culvert_dry.sww')
         
     def test_that_culvert_flows_conserves_volume(self):

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -644 +644 @@
             fid.close()
 
+
+##
+# @brief Class to open an sww file so that domain can be populated with quantity values 
+class Read_sww:
+
+    def __init__(self, source):
+        """The source parameter is assumed to be a NetCDF sww file.
+        """
+
+        self.source = source
+
+        self.frame_number = 0
+
+        fin = NetCDFFile(self.source, 'r')
+
+        self.time = num.array(fin.variables['time'], num.float)
+        self.last_frame_number = self.time.shape[0] - 1
+
+        self.frames = num.arange(self.last_frame_number+1)
+
+        fin.close()
+        
+        self.read_mesh()
+
+        self.quantities = {}
+
+        self.read_quantities()
+
+
+    def read_mesh(self):
+        fin = NetCDFFile(self.source, 'r')
+
+        self.vertices = num.array(fin.variables['volumes'], num.int)
+        
+        self.x = x = num.array(fin.variables['x'], num.float)
+        self.y = y = num.array(fin.variables['y'], num.float)
+
+        assert len(self.x) == len(self.y)
+        
+        self.xmin = num.min(x)
+        self.xmax = num.max(x)
+        self.ymin = num.min(y)
+        self.ymax = num.max(y)
+
+
+
+        fin.close()
+        
+    def read_quantities(self, frame_number=0):
+
+        assert frame_number >= 0 and frame_number <= self.last_frame_number
+
+        self.frame_number = frame_number
+
+        M = len(self.x)/3
+        
+        fin = NetCDFFile(self.source, 'r')
+        
+        for q in filter(lambda n:n != 'x' and n != 'y' and n != 'z' and n != 'time' and n != 'volumes' and \
+                        '_range' not in n, \
+                        fin.variables.keys()):
+            if len(fin.variables[q].shape) == 1: # Not a time-varying quantity
+                self.quantities[q] = num.ravel(num.array(fin.variables[q], num.float)).reshape(M,3)
+            else: # Time-varying, get the current timestep data
+                self.quantities[q] = num.array(fin.variables[q][self.frame_number], num.float).reshape(M,3)
+        fin.close()
+        return self.quantities
+
+    def get_bounds(self):
+        return [self.xmin, self.xmax, self.ymin, self.ymax]
+
+    def get_last_frame_number(self):
+        return self.last_frame_number
+
+    def get_time(self):
+        return self.time[self.frame_number]
 
 
@@ -2530 +2606 @@
             verbose = False,
             origin = None):
+    
     log.critical('sww2asc will soon be obsoleted - please use sww2dem')
     sww2dem(basename_in,
@@ -2540 +2617 @@
             verbose = verbose,
             origin = origin,
-        datum = 'WGS84',
-        format = 'asc')
+            datum = 'WGS84',
+            format = 'asc')
 
 
@@ -6379 +6456 @@
 
 
+
+
 ##
 # @brief A class to write STS files.

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -102 +102 @@
 from anuga.geospatial_data.geospatial_data import ensure_geospatial
 
-from anuga.config import minimum_storable_height
-from anuga.config import minimum_allowed_height, maximum_allowed_speed
-from anuga.config import g, epsilon, beta_w, beta_w_dry,\
-     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
-from anuga.config import alpha_balance
-from anuga.config import optimise_dry_cells
-from anuga.config import optimised_gradient_limiter
-from anuga.config import use_edge_limiter
-from anuga.config import use_centroid_velocities
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 
@@ -199 +190 @@
                                 number_of_full_triangles=number_of_full_triangles)
 
+        self.set_defaults()
+
+ 
+        self.forcing_terms.append(manning_friction_implicit)
+        self.forcing_terms.append(gravity)
+
+        # Stored output
+        self.store = True
+        self.set_store_vertices_uniquely(False)
+
+        self.quantities_to_be_stored = {'elevation': 1, 
+                                        'stage': 2, 
+                                        'xmomentum': 2, 
+                                        'ymomentum': 2}
+
+
+
+
+    ##
+    # @brief Set default values, usually read in from a config file
+    # @param flag
+    def set_defaults(self):
+        """Set the default values in this routine. That way we can inherit class
+        and just over redefine the defaults for the new class
+        """
+
+        from anuga.config import minimum_storable_height
+        from anuga.config import minimum_allowed_height, maximum_allowed_speed
+        from anuga.config import g, epsilon, beta_w, beta_w_dry,\
+             beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
+        from anuga.config import alpha_balance
+        from anuga.config import optimise_dry_cells
+        from anuga.config import optimised_gradient_limiter
+        from anuga.config import use_edge_limiter
+        from anuga.config import use_centroid_velocities
+
+
+
         self.set_minimum_allowed_height(minimum_allowed_height)
-
         self.maximum_allowed_speed = maximum_allowed_speed
+
         self.g = g
         self.beta_w = beta_w
@@ -214 +243 @@
         self.optimise_dry_cells = optimise_dry_cells
 
-        self.use_new_mannings = False
-        self.forcing_terms.append(manning_friction_implicit)
-        self.forcing_terms.append(gravity)
-
-        # Stored output
-        self.store = True
-        self.set_store_vertices_uniquely(False)
+        
+        self.set_new_mannings_function(False)
+
         self.minimum_storable_height = minimum_storable_height
-        self.quantities_to_be_stored = {'elevation': 1, 
-                                        'stage': 2, 
-                                        'xmomentum': 2, 
-                                        'ymomentum': 2}
-
-        # Limiters
+
+         # Limiters
         self.use_edge_limiter = use_edge_limiter
         self.optimised_gradient_limiter = optimised_gradient_limiter
-        self.use_centroid_velocities = use_centroid_velocities
-
+        self.use_centroid_velocities = use_centroid_velocities       
 
     ##

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -6677 +6677 @@
                 assert num.allclose(q, ref_flow)
 
+        os.remove('Inflow_flowline_test.sww') 
 
     def test_volume_conservation_inflow(self):
@@ -6776 +6777 @@
 
 
+        os.remove('Inflow_volume_test.sww')
+
+
+        
     def test_volume_conservation_rain(self):
         """test_volume_conservation
@@ -6881 +6886 @@
             ref_volume += ys * delta_V
 
+        os.remove('Rain_volume_test.sww')
 
     def Xtest_rain_conservation_and_runoff(self):
@@ -7084 +7090 @@
         # Check that quantities have been stored correctly    
         from Scientific.IO.NetCDF import NetCDFFile
-        fid = NetCDFFile(domain.get_name() + '.sww')
+        sww_file = domain.get_name() + '.sww'
+        fid = NetCDFFile(sww_file)
 
         x = fid.variables['x'][:]
@@ -7091 +7098 @@
         elevation = fid.variables['elevation'][:]
         fid.close()
+
+        os.remove(sww_file)
+        
                    
         assert len(stage.shape) == 2
@@ -7261 +7271 @@
                        % (q, ref_flow))
                 assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
+
+        os.remove('inflow_flowline_test.sww')
 
         
@@ -7552 +7564 @@
                                % (domain_depth, normal_depth))
 
+        os.remove('Inflow_flowline_test.sww')
+
 #################################################################################
 

anuga_core/source/anuga/shallow_water_balanced/test_swb_distribute.py

@@ -147 +147 @@
 
         # First order
-        domain._order_ = 1
-        domain.distribute_to_vertices_and_edges()
+        domain.set_default_order(1)
+        domain.distribute_to_vertices_and_edges()
+        
         assert num.allclose(L[1], val1)
 
         # Second order
-        domain._order_ = 2
-        domain.beta_w = 0.9
-        domain.beta_w_dry = 0.9
-        domain.beta_uh = 0.9
-        domain.beta_uh_dry = 0.9
-        domain.beta_vh = 0.9
-        domain.beta_vh_dry = 0.9
-        domain.distribute_to_vertices_and_edges()
-        assert num.allclose(L[1], [2.2, 4.9, 4.9])
+        domain.set_default_order(2)
+        domain.distribute_to_vertices_and_edges()
+
+        assert num.allclose(L[1], [1.00000000e-03,   5.99950000e+00,  5.99950000e+00])
+
+        assert num.allclose(val1, num.sum(L[1])/3)
 
     def test_distribute_away_from_bed(self):
@@ -192 +190 @@
         assert num.allclose(b[1], 0.0)
 
-        domain._order_ = 1
-        domain.distribute_to_vertices_and_edges()
+        domain.set_default_order(1)
+        domain.distribute_to_vertices_and_edges()
+        
         assert num.allclose(L[1], 1.77777778)
 
-        domain._order_ = 2
-        domain.beta_w = 0.9
-        domain.beta_w_dry = 0.9
-        domain.beta_uh = 0.9
-        domain.beta_uh_dry = 0.9
-        domain.beta_vh = 0.9
-        domain.beta_vh_dry = 0.9
-        domain.distribute_to_vertices_and_edges()
-        assert num.allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
+        domain.set_default_order(2)
+        domain.distribute_to_vertices_and_edges()
+
+        assert num.allclose(L[1], [1.00000000e-03,   2.66616667e+00,   2.66616667e+00])
+
+        assert num.allclose(1.77777778, num.sum(L[1])/3)
 
     def test_distribute_away_from_bed1(self):
@@ -234 +230 @@
         assert num.allclose(b[1], 3.33333333)
 
-        domain._order_ = 1
-        domain.distribute_to_vertices_and_edges()
+        domain.set_default_order(1)
+        domain.distribute_to_vertices_and_edges()
+        
         assert num.allclose(L[1], 4.9382716)
 
-        domain._order_ = 2
-        domain.beta_w = 0.9
-        domain.beta_w_dry = 0.9
-        domain.beta_uh = 0.9
-        domain.beta_uh_dry = 0.9
-        domain.beta_vh = 0.9
-        domain.beta_vh_dry = 0.9
-        domain.distribute_to_vertices_and_edges()
-        assert num.allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
+        domain.set_default_order(2)
+        domain.distribute_to_vertices_and_edges()
+    
+        assert num.allclose(L[1], [ 1.00000000e-03,   6.88207932e+00,   7.93173549e+00])
+
 
     def test_distribute_near_bed(self):
@@ -283 +276 @@
                                 domain.quantities['stage'].centroid_values[i])
 
-        domain._order_ = 1
-
-        domain.tight_slope_limiters = 0
-        domain.distribute_to_vertices_and_edges()
-        assert num.allclose(L[1], [0.1, 20.1, 20.1])
+        domain.set_default_order(1)
+        domain.distribute_to_vertices_and_edges()
+
+        
+        assert num.allclose(L[1], [0.298,  20.001,  20.001])
         for i in range(len(L)):
             assert num.allclose(volumes[i], num.sum(L[i])/3)
 
-        # Allow triangle to be flatter (closer to bed)
-        domain.tight_slope_limiters = 1
-
-        domain.distribute_to_vertices_and_edges()
-        assert num.allclose(L[1], [0.298, 20.001, 20.001])
-        for i in range(len(L)):
-            assert num.allclose(volumes[i], num.sum(L[i])/3)
-
-        domain._order_ = 2
-
-        domain.tight_slope_limiters = 0
-        domain.distribute_to_vertices_and_edges()
-        assert num.allclose(L[1], [0.1, 20.1, 20.1])
-        for i in range(len(L)):
-            assert num.allclose(volumes[i], num.sum(L[i])/3)
-
-        # Allow triangle to be flatter (closer to bed)
-        domain.tight_slope_limiters = 1
-
-        domain.distribute_to_vertices_and_edges()
-        assert num.allclose(L[1], [0.298, 20.001, 20.001])
+        domain.set_default_order(2)
+        domain.distribute_to_vertices_and_edges()
+
+        
+        assert num.allclose(L[1], [  0.1,  20.1,  20.1])
         for i in range(len(L)):
             assert num.allclose(volumes[i], num.sum(L[i])/3)
@@ -378 +355 @@
 
         domain.distribute_to_vertices_and_edges()
-        assert num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\
-               num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\
-               num.allclose(L[1], [4.19351461, 16.10548539, 20.001]) # old limiters
+
+
+
+
+        assert num.allclose(L[1], [4.001, 16.15377472, 20.14522528],rtol=1.0e-3)
         
         for i in range(len(L)):
@@ -440 +419 @@
         domain.distribute_to_vertices_and_edges()
 
-        assert num.allclose(L[1,:], [-0.00825735775384,
-                                     -0.00801881482869,
-                                     0.0272445025825])
-        assert num.allclose(X[1,:], [0.0143507718962,
-                                     0.0142502147066,
-                                     0.00931268339717])
-        assert num.allclose(Y[1,:], [-0.000117062180693,
-                                     7.94434448109e-005,
-                                     -0.000151505429018])
+
+        assert num.allclose(L[1,:],[-0.01434766, -0.01292565, 0.03824164])
+
+        assert num.allclose(X[1,:],[ 0.01649334,  0.016267,    0.00515332])
+
+        assert num.allclose(Y[1,:],[-0.00027927,  0.00050728, -0.00041714])
 
 

anuga_core/source/anuga/shallow_water_balanced/test_swb_forcing_terms.py

@@ -24 +24 @@
 
 
-
-
-class Test_swb_clean(unittest.TestCase):
+# Variable windfield implemented using functions
+def speed(t, x, y):
+    """Large speeds halfway between center and edges
+
+    Low speeds at center and edges
+    """
+
+    from math import exp, cos, pi
+
+    x = num.array(x)
+    y = num.array(y)
+
+    N = len(x)
+    s = 0*x  #New array
+
+    for k in range(N):
+        r = num.sqrt(x[k]**2 + y[k]**2)
+        factor = exp(-(r-0.15)**2)
+        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
+
+    return s
+
+def scalar_func(t, x, y):
+    """Function that returns a scalar.
+
+    Used to test error message when numeric array is expected
+    """
+
+    return 17.7
+
+def scalar_func_list(t, x, y):
+    """Function that returns a scalar.
+
+    Used to test error message when numeric array is expected
+    """
+
+    return [17.7]
+
+
+def angle(t, x, y):
+    """Rotating field
+    """
+    from math import atan, pi
+
+    x = num.array(x)
+    y = num.array(y)
+
+    N = len(x)
+    a = 0 * x    # New array
+
+    for k in range(N):
+        r = num.sqrt(x[k]**2 + y[k]**2)
+
+        angle = atan(y[k]/x[k])
+
+        if x[k] < 0:
+            angle += pi
+
+        # Take normal direction
+        angle -= pi/2
+
+        # Ensure positive radians
+        if angle < 0:
+            angle += 2*pi
+
+        a[k] = angle/pi*180
+
+    return a
+
+
+###############################################################################
+
+class Test_swb_forcing_terms(unittest.TestCase):
     def setUp(self):
         pass
@@ -1723 +1793 @@
 
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_swb_clean, 'test')
+    suite = unittest.makeSuite(Test_swb_forcing_terms, 'test')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)

anuga_core/source/anuga/test_all.py

@@ -23 +23 @@
 
 # Directories that should not be searched for test files.
-exclude_dirs = ['pypar_dist',                        # Special requirements
-                '.svn',                              # subversion
+exclude_dirs = ['pypar_dist',    # Special requirements
+                '.svn',          # subversion
                 'props', 'wcprops', 'prop-base', 'text-base', 'tmp']
 

anuga_core/source/anuga_parallel/interface.py

@@ -8 +8 @@
 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
+from anuga_parallel.parallel_api import pypar_available, barrier, finalize
 
 from anuga_parallel.parallel_meshes import parallel_rectangle

anuga_core/source/anuga_parallel/test_parallel_distribute_domain.py

@@ -35 +35 @@
 
 
-from anuga_parallel.interface import distribute, myid, numprocs
+from anuga_parallel.interface import distribute, myid, numprocs, finalize
 
 
@@ -222 +222 @@
 
 
+
+    finalize()

anuga_core/source/anuga_parallel/test_parallel_distribute_mesh.py

@@ -15 +15 @@
 from anuga_parallel.distribute_mesh import extract_hostmesh, rec_submesh, send_submesh
 
-from anuga_parallel.interface import myid, numprocs, barrier
+from anuga_parallel.interface import myid, numprocs, barrier, finalize
 
 import numpy as num
@@ -302 +302 @@
         distibute_three_processors()
 
-
+    finalize()
         
 

anuga_core/source/anuga_parallel/test_parallel_sw_flow.py

@@ -27 +27 @@
 
 
-from anuga_parallel.interface import distribute, myid, numprocs, send, receive, barrier
+from anuga_parallel.interface import distribute, myid, numprocs, send, receive, barrier, finalize
 
 #--------------------------------------------------------------------------
@@ -37 +37 @@
 N = 25
 M = 25
-verbose = False
+verbose = True 
 
 #---------------------------------
@@ -212 +212 @@
         evolution_test(parallel=True, G=G, seq_interpolation_points = interpolation_points)
         
-
-
-
+        finalize()
+
+