Changeset 5873

Timestamp:

Oct 28, 2008, 5:19:47 PM (16 years ago)

Author:

ole

Message:

Implemented default_rate to allow ANUGA to continue beyond data used by forcing functions such as Rainfall
This should address ticket:265

Location:

anuga_core/source/anuga/shallow_water

Files:

• shallow_water_domain.py (modified) (11 diffs)
• test_shallow_water_domain.py (modified) (3 diffs)

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -109 +109 @@
 from anuga.config import use_centroid_velocities
 
+from anuga.fit_interpolate.interpolate import Modeltime_too_late, Modeltime_too_early
 
 from anuga.utilities.polygon import inside_polygon, polygon_area, is_inside_polygon
 
+from types import IntType, FloatType
 
 
@@ -1643 +1645 @@
     radius [m]: Size of circular area
     polygon:    Arbitrary polygon.
+    default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
+                  Admissible types: None, constant number or function of t
 
 
@@ -1660 +1664 @@
                  center=None, radius=None,
                  polygon=None,
+                 default_rate=None,
                  verbose=False):
                      
@@ -1768 +1773 @@
                 msg += 'specified region: %s' %inlet_region
                 raise Exception, msg
-
+            
+        # Check and store default_rate
+        msg = 'Keyword argument default_rate must be either None ' 
+        msg += 'or a function of time.\n I got %s' %(str(default_rate)) 
+        assert default_rate is None or \
+               type(default_rate) in [IntType, FloatType] or \
+               callable(default_rate), msg
+        
+        if default_rate is not None:
+
+            # If it is a constant, make it a function
+            if not callable(default_rate):
+                tmp = default_rate
+                default_rate = lambda t: tmp
+
+            
+            # Check that default_rate is a function of one argument
+            try:
+                default_rate(0.0)
+            except:
+                raise Exception, msg
+
+        self.default_rate = default_rate
+        
 
     def __call__(self, domain):
@@ -1775 +1803 @@
 
         # Call virtual method allowing local modifications
-        rate = self.update_rate(domain.get_time())
+
+        t = domain.get_time()
+        try:
+            rate = self.update_rate(t)
+        except (Modeltime_too_late, Modeltime_too_early), e: 
+            if self.default_rate is None: 
+                raise Exception, e # Reraise exception
+            else:
+                # FIXME: Issue a warning first time this happens (See changeset:5657)
+                rate = self.default_rate(t)
+                
+
         if rate is None:
             msg = 'Attribute rate must be specified in General_forcing'
@@ -1879 +1918 @@
                  center=None, radius=None,
                  polygon=None,
+                 default_rate=None,                 
                  verbose=False):
 
@@ -1885 +1925 @@
             rain = lambda t: rate(t)/1000.0
         else:
-            rain = rate/1000.0            
+            rain = rate/1000.0
+
+        if default_rate is not None:    
+            if callable(default_rate):
+                default_rain = lambda t: default_rate(t)/1000.0
+            else:
+                default_rain = default_rate/1000.0
+        else:
+            default_rain = None
+
+
+            
             
         General_forcing.__init__(self,
@@ -1893 +1944 @@
                                  center=center, radius=radius,
                                  polygon=polygon,
+                                 default_rate=default_rain,                                 
                                  verbose=verbose)
 
@@ -1955 +2007 @@
                  center=None, radius=None,
                  polygon=None,
+                 default_rate=None,                                                  
                  verbose=False):                 
 
@@ -1965 +2018 @@
                                  center=center, radius=radius,
                                  polygon=polygon,
+                                 default_rate=default_rate,
                                  verbose=verbose)
 
@@ -1975 +2029 @@
         """
 
-        
-        
         if callable(self.rate):
             _rate = self.rate(t)/self.exchange_area

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -2401 +2401 @@
 
 
-        # FIXME: Do Time dependent rainfall
-
-
 
     def test_rainfall_restricted_by_polygon(self):
@@ -2445 +2442 @@
         assert allclose(domain.quantities['stage'].explicit_update[2:], 0)        
         
-        # FIXME: Do Time dependent rainfall with poly
+
+
+    def test_time_dependent_rainfall_restricted_by_polygon(self):
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+
+        domain = Domain(points, vertices)
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+        # Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce)
+        domain.forcing_terms = []
+        R = Rainfall(domain, rate=lambda t: 3*t + 7, polygon = [[1,1], [2,1], [2,2], [1,2]])
+
+        assert allclose(R.exchange_area, 1)
+        
+        domain.forcing_terms.append(R)
+
+
+        domain.time = 10.
+
+        domain.compute_forcing_terms()
+        #print domain.quantities['stage'].explicit_update
+        
+        assert allclose(domain.quantities['stage'].explicit_update[1], (3*domain.time+7)/1000)
+        assert allclose(domain.quantities['stage'].explicit_update[0], 0)
+        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)        
+        
+
+
+    def test_time_dependent_rainfall_restricted_by_polygon_with_default(self):
+        """test_time_dependent_rainfall_restricted_by_polygon_with_default
+
+        Test that default rainfall can be used when given rate runs out of data.
+        """
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+
+        domain = Domain(points, vertices)
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+        # Setup only one forcing term, time dependent rainfall that expires at t==20
+        from anuga.fit_interpolate.interpolate import Modeltime_too_late
+        def main_rate(t):
+            if t > 20:
+                msg = 'Model time exceeded.'
+                raise Modeltime_too_late, msg
+            else:
+                return 3*t + 7
+        
+        domain.forcing_terms = []
+        R = Rainfall(domain, rate=main_rate, polygon = [[1,1], [2,1], [2,2], [1,2]],
+                     default_rate=5.0)
+
+        assert allclose(R.exchange_area, 1)
+        
+        domain.forcing_terms.append(R)
+
+
+        domain.time = 10.
+
+        domain.compute_forcing_terms()
+        #print domain.quantities['stage'].explicit_update
+        
+        assert allclose(domain.quantities['stage'].explicit_update[1], (3*domain.time+7)/1000)
+        assert allclose(domain.quantities['stage'].explicit_update[0], 0)
+        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)        
+
+
+        domain.time = 100.
+        domain.quantities['stage'].explicit_update[:] = 0.0  # Reset
+        domain.compute_forcing_terms()
+        #print domain.quantities['stage'].explicit_update
+        
+        assert allclose(domain.quantities['stage'].explicit_update[1], 5.0/1000) # Default value
+        assert allclose(domain.quantities['stage'].explicit_update[0], 0)
+        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)        
+        
 
 
@@ -6013 +6119 @@
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')
-    #suite = unittest.makeSuite(Test_Shallow_Water,'test_temp')    
+    #suite = unittest.makeSuite(Test_Shallow_Water,'test_time_dependent_rainfall_restricted_by_polygon')