Changeset 5178

Timestamp:

Mar 26, 2008, 5:25:59 PM (17 years ago)

Author:

ole

Message:

Rainfall restricted to polygon

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (3 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -273 +273 @@
             #raise Exception(msg)
 
-        
+        # Test that file function can be called
         q = self.F(0, point_id=0)
 

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -1502 +1502 @@
                   the inflow region and then applied to update the
                   quantity in question. 
+
+    polygon: Specifies a polygon to restrict the rainfall.
     
     Examples
@@ -1523 +1525 @@
     # FIXME (OLE): Generalise to all quantities
 
+    
     def __init__(self,
                  rain=0.0,
-                 quantity_name='stage'):
+                 quantity_name='stage',
+                 polygon=None):
 
         self.rain = rain
         self.quantity_name = quantity_name
+        self.polygon = polygon
     
     def __call__(self, domain):
+        
+        # FIXME(Ole): Move this to top of file eventually
+        from anuga.utilities.polygon import inside_polygon        
 
         # Update rainfall
@@ -1539 +1547 @@
 
         # Now rain is a number
+        # Converting mm/s to m/s to apply in ANUGA
+        # 1mm of rainfall is equivalent to 1 litre /m2 
+        # Flow is expressed as m3/s converted to a stage height in (m)
+        
+        # Note 1m3 = 1x10^9mm3 (mls)
+        # or is that m3 to Litres ??? Check this how is it applied !!!
+        rain_fall = rain/1000
+
+        # Find indices for coordinates inside polygon if specified
+        indices = None
+        if self.polygon is not None:
+            points = domain.get_centroid_coordinates(absolute=True)
+            indices = inside_polygon(points, self.polygon)
+            
+
+        # Assign rainfall value to explicit update vector to be used as
+        # a forcing term
         quantity = domain.quantities[self.quantity_name].explicit_update
-        quantity[:] += rain/1000  # Converting mm/s to m/s to apply in ANUGA
-                # 1mm of rainfall is equivalent to 1 litre /m2 
-                # Flow is expressed as m3/s converted to a stage height in (m)
-                
-                # Note 1m3 = 1x10^9mm3 (mls)
-                # or is that m3 to Litres ??? Check this how is it applied !!!
+        if indices is None:
+            quantity[:] += rain_fall
+        else:
+            # Brute force assignment of restricted rainfall
+            for i in indices:
+                quantity[i] += rain_fall
+
+        
+
+
 
 

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -1878 +1878 @@
 
 
+
+
     def test_windfield_from_file(self):
         from anuga.config import rho_a, rho_w, eta_w
@@ -2160 +2162 @@
             msg = 'Should have raised exception'
             raise msg
+
+
+
+    def test_rainfall(self):
+        from math import pi, cos, sin
+
+        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, constant rainfall
+        domain.forcing_terms = []
+        domain.forcing_terms.append( Rainfall(rain=2.0) )
+
+        domain.compute_forcing_terms()
+        assert allclose(domain.quantities['stage'].explicit_update, 2.0/1000)
+
+
+        # FIXME: Do Time dependent rainfall
+
+
+
+    def test_rainfall_restricted_by_polygon(self):
+        from math import pi, cos, sin
+
+        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, constant rainfall restricted to a polygon enclosing triangle #1 (bce)
+        domain.forcing_terms = []
+        domain.forcing_terms.append( Rainfall(rain=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]]))
+
+        domain.compute_forcing_terms()
+        #print domain.quantities['stage'].explicit_update
+        
+        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/1000)
+        assert allclose(domain.quantities['stage'].explicit_update[0], 0)
+        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)        
+        
+        # FIXME: Do Time dependent rainfall with poly
+