Changeset 5178


Ignore:
Timestamp:
Mar 26, 2008, 5:25:59 PM (15 years ago)
Author:
ole
Message:

Rainfall restricted to polygon

Location:
anuga_core/source/anuga
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r4820 r5178  
    273273            #raise Exception(msg)
    274274
    275        
     275        # Test that file function can be called
    276276        q = self.F(0, point_id=0)
    277277
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r5176 r5178  
    15021502                  the inflow region and then applied to update the
    15031503                  quantity in question.
     1504
     1505    polygon: Specifies a polygon to restrict the rainfall.
    15041506   
    15051507    Examples
     
    15231525    # FIXME (OLE): Generalise to all quantities
    15241526
     1527   
    15251528    def __init__(self,
    15261529                 rain=0.0,
    1527                  quantity_name='stage'):
     1530                 quantity_name='stage',
     1531                 polygon=None):
    15281532
    15291533        self.rain = rain
    15301534        self.quantity_name = quantity_name
     1535        self.polygon = polygon
    15311536   
    15321537    def __call__(self, domain):
     1538       
     1539        # FIXME(Ole): Move this to top of file eventually
     1540        from anuga.utilities.polygon import inside_polygon       
    15331541
    15341542        # Update rainfall
     
    15391547
    15401548        # Now rain is a number
     1549        # Converting mm/s to m/s to apply in ANUGA
     1550        # 1mm of rainfall is equivalent to 1 litre /m2
     1551        # Flow is expressed as m3/s converted to a stage height in (m)
     1552       
     1553        # Note 1m3 = 1x10^9mm3 (mls)
     1554        # or is that m3 to Litres ??? Check this how is it applied !!!
     1555        rain_fall = rain/1000
     1556
     1557        # Find indices for coordinates inside polygon if specified
     1558        indices = None
     1559        if self.polygon is not None:
     1560            points = domain.get_centroid_coordinates(absolute=True)
     1561            indices = inside_polygon(points, self.polygon)
     1562           
     1563
     1564        # Assign rainfall value to explicit update vector to be used as
     1565        # a forcing term
    15411566        quantity = domain.quantities[self.quantity_name].explicit_update
    1542         quantity[:] += rain/1000  # Converting mm/s to m/s to apply in ANUGA
    1543                 # 1mm of rainfall is equivalent to 1 litre /m2
    1544                 # Flow is expressed as m3/s converted to a stage height in (m)
    1545                
    1546                 # Note 1m3 = 1x10^9mm3 (mls)
    1547                 # or is that m3 to Litres ??? Check this how is it applied !!!
     1567        if indices is None:
     1568            quantity[:] += rain_fall
     1569        else:
     1570            # Brute force assignment of restricted rainfall
     1571            for i in indices:
     1572                quantity[i] += rain_fall
     1573
     1574       
     1575
     1576
    15481577
    15491578
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r5175 r5178  
    18781878
    18791879
     1880
     1881
    18801882    def test_windfield_from_file(self):
    18811883        from anuga.config import rho_a, rho_w, eta_w
     
    21602162            msg = 'Should have raised exception'
    21612163            raise msg
     2164
     2165
     2166
     2167    def test_rainfall(self):
     2168        from math import pi, cos, sin
     2169
     2170        a = [0.0, 0.0]
     2171        b = [0.0, 2.0]
     2172        c = [2.0, 0.0]
     2173        d = [0.0, 4.0]
     2174        e = [2.0, 2.0]
     2175        f = [4.0, 0.0]
     2176
     2177        points = [a, b, c, d, e, f]
     2178        #bac, bce, ecf, dbe
     2179        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2180
     2181
     2182        domain = Domain(points, vertices)
     2183
     2184        #Flat surface with 1m of water
     2185        domain.set_quantity('elevation', 0)
     2186        domain.set_quantity('stage', 1.0)
     2187        domain.set_quantity('friction', 0)
     2188
     2189        Br = Reflective_boundary(domain)
     2190        domain.set_boundary({'exterior': Br})
     2191
     2192        # Setup only one forcing term, constant rainfall
     2193        domain.forcing_terms = []
     2194        domain.forcing_terms.append( Rainfall(rain=2.0) )
     2195
     2196        domain.compute_forcing_terms()
     2197        assert allclose(domain.quantities['stage'].explicit_update, 2.0/1000)
     2198
     2199
     2200        # FIXME: Do Time dependent rainfall
     2201
     2202
     2203
     2204    def test_rainfall_restricted_by_polygon(self):
     2205        from math import pi, cos, sin
     2206
     2207        a = [0.0, 0.0]
     2208        b = [0.0, 2.0]
     2209        c = [2.0, 0.0]
     2210        d = [0.0, 4.0]
     2211        e = [2.0, 2.0]
     2212        f = [4.0, 0.0]
     2213
     2214        points = [a, b, c, d, e, f]
     2215        #bac, bce, ecf, dbe
     2216        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2217
     2218
     2219        domain = Domain(points, vertices)
     2220
     2221        #Flat surface with 1m of water
     2222        domain.set_quantity('elevation', 0)
     2223        domain.set_quantity('stage', 1.0)
     2224        domain.set_quantity('friction', 0)
     2225
     2226        Br = Reflective_boundary(domain)
     2227        domain.set_boundary({'exterior': Br})
     2228
     2229        # Setup only one forcing term, constant rainfall restricted to a polygon enclosing triangle #1 (bce)
     2230        domain.forcing_terms = []
     2231        domain.forcing_terms.append( Rainfall(rain=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]]))
     2232
     2233        domain.compute_forcing_terms()
     2234        #print domain.quantities['stage'].explicit_update
     2235       
     2236        assert allclose(domain.quantities['stage'].explicit_update[1], 2.0/1000)
     2237        assert allclose(domain.quantities['stage'].explicit_update[0], 0)
     2238        assert allclose(domain.quantities['stage'].explicit_update[2:], 0)       
     2239       
     2240        # FIXME: Do Time dependent rainfall with poly
     2241
    21622242
    21632243
Note: See TracChangeset for help on using the changeset viewer.