Ignore:
Timestamp:
May 19, 2010, 10:13:06 AM (13 years ago)
Author:
hudson
Message:

Fixed unit test failures.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r7731 r7733  
    8383import numpy as num
    8484
    85 from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
    8685from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain
    87 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    88      import Boundary
    89 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    90      import File_boundary
    91 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    92      import Dirichlet_boundary
    93 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    94      import Time_boundary
    95 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    96      import Transmissive_boundary
    97 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    98      import AWI_boundary
    99 
     86
     87from anuga.shallow_water.forcing import Cross_section
    10088from anuga.pmesh.mesh_interface import create_mesh_from_regions
    10189from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
     
    14321420
    14331421
    1434 ################################################################################
    1435 # Experimental auxiliary functions
    1436 ################################################################################
    1437 
    1438 ##
    1439 # @brief Check forcefield parameter.
    1440 # @param f Object to check.
    1441 # @note 'f' may be a callable object or a scalar value.
    1442 def check_forcefield(f):
    1443     """Check that force object is as expected.
    1444    
    1445     Check that f is either:
    1446     1: a callable object f(t,x,y), where x and y are vectors
    1447        and that it returns an array or a list of same length
    1448        as x and y
    1449     2: a scalar
    1450     """
    1451 
    1452     if callable(f):
    1453         N = 3
    1454         x = num.ones(3, num.float)
    1455         y = num.ones(3, num.float)
    1456         try:
    1457             q = f(1.0, x=x, y=y)
    1458         except Exception, e:
    1459             msg = 'Function %s could not be executed:\n%s' %(f, e)
    1460             # FIXME: Reconsider this semantics
    1461             raise Exception, msg
    1462 
    1463         try:
    1464             q = num.array(q, num.float)
    1465         except:
    1466             msg = ('Return value from vector function %s could not '
    1467                    'be converted into a numeric array of floats.\nSpecified '
    1468                    'function should return either list or array.' % f)
    1469             raise Exception, msg
    1470 
    1471         # Is this really what we want?
    1472         # info is "(func name, filename, defining line)"
    1473         func_info = (f.func_name, f.func_code.co_filename,
    1474                      f.func_code.co_firstlineno)
    1475         func_msg = 'Function %s (defined in %s, line %d)' % func_info
    1476         try:
    1477             result_len = len(q)
    1478         except:
    1479             msg = '%s must return vector' % func_msg
    1480             self.fail(msg)
    1481         msg = '%s must return vector of length %d' % (func_msg, N)
    1482         assert result_len == N, msg
    1483     else:
    1484         try:
    1485             f = float(f)
    1486         except:
    1487             msg = ('Force field %s must be a scalar value coercible to float.'
    1488                    % str(f))
    1489             raise Exception, msg
    1490 
    1491     return f
    1492 
    1493 
    1494 ##
    1495 # Class to apply a wind stress to a domain.
    1496 class Wind_stress:
    1497     """Apply wind stress to water momentum in terms of
    1498     wind speed [m/s] and wind direction [degrees]
    1499     """
    1500 
    1501     ##
    1502     # @brief Create an instance of Wind_stress.
    1503     # @param *args
    1504     # @param **kwargs
    1505     def __init__(self, *args, **kwargs):
    1506         """Initialise windfield from wind speed s [m/s]
    1507         and wind direction phi [degrees]
    1508 
    1509         Inputs v and phi can be either scalars or Python functions, e.g.
    1510 
    1511         W = Wind_stress(10, 178)
    1512 
    1513         #FIXME - 'normal' degrees are assumed for now, i.e. the
    1514         vector (1,0) has zero degrees.
    1515         We may need to convert from 'compass' degrees later on and also
    1516         map from True north to grid north.
    1517 
    1518         Arguments can also be Python functions of t,x,y as in
    1519 
    1520         def speed(t,x,y):
    1521             ...
    1522             return s
    1523 
    1524         def angle(t,x,y):
    1525             ...
    1526             return phi
    1527 
    1528         where x and y are vectors.
    1529 
    1530         and then pass the functions in
    1531 
    1532         W = Wind_stress(speed, angle)
    1533 
    1534         The instantiated object W can be appended to the list of
    1535         forcing_terms as in
    1536 
    1537         Alternatively, one vector valued function for (speed, angle)
    1538         can be applied, providing both quantities simultaneously.
    1539         As in
    1540         W = Wind_stress(F), where returns (speed, angle) for each t.
    1541 
    1542         domain.forcing_terms.append(W)
    1543         """
    1544 
    1545         from anuga.config import rho_a, rho_w, eta_w
    1546 
    1547         if len(args) == 2:
    1548             s = args[0]
    1549             phi = args[1]
    1550         elif len(args) == 1:
    1551             # Assume vector function returning (s, phi)(t,x,y)
    1552             vector_function = args[0]
    1553             s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
    1554             phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
    1555         else:
    1556            # Assume info is in 2 keyword arguments
    1557            if len(kwargs) == 2:
    1558                s = kwargs['s']
    1559                phi = kwargs['phi']
    1560            else:
    1561                raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
    1562 
    1563         self.speed = check_forcefield(s)
    1564         self.phi = check_forcefield(phi)
    1565 
    1566         self.const = eta_w*rho_a/rho_w
    1567 
    1568     ##
    1569     # @brief 'execute' this class instance.
    1570     # @param domain
    1571     def __call__(self, domain):
    1572         """Evaluate windfield based on values found in domain"""
    1573 
    1574         from math import pi, cos, sin, sqrt
    1575 
    1576         xmom_update = domain.quantities['xmomentum'].explicit_update
    1577         ymom_update = domain.quantities['ymomentum'].explicit_update
    1578 
    1579         N = len(domain)    # number_of_triangles
    1580         t = domain.time
    1581 
    1582         if callable(self.speed):
    1583             xc = domain.get_centroid_coordinates()
    1584             s_vec = self.speed(t, xc[:,0], xc[:,1])
    1585         else:
    1586             # Assume s is a scalar
    1587             try:
    1588                 s_vec = self.speed * num.ones(N, num.float)
    1589             except:
    1590                 msg = 'Speed must be either callable or a scalar: %s' %self.s
    1591                 raise msg
    1592 
    1593         if callable(self.phi):
    1594             xc = domain.get_centroid_coordinates()
    1595             phi_vec = self.phi(t, xc[:,0], xc[:,1])
    1596         else:
    1597             # Assume phi is a scalar
    1598 
    1599             try:
    1600                 phi_vec = self.phi * num.ones(N, num.float)
    1601             except:
    1602                 msg = 'Angle must be either callable or a scalar: %s' %self.phi
    1603                 raise msg
    1604 
    1605         assign_windfield_values(xmom_update, ymom_update,
    1606                                 s_vec, phi_vec, self.const)
    1607 
    1608 
    1609 ##
    1610 # @brief Assign wind field values
    1611 # @param xmom_update
    1612 # @param ymom_update
    1613 # @param s_vec
    1614 # @param phi_vec
    1615 # @param const
    1616 def assign_windfield_values(xmom_update, ymom_update,
    1617                             s_vec, phi_vec, const):
    1618     """Python version of assigning wind field to update vectors.
    1619     A C version also exists (for speed)
    1620     """
    1621 
    1622     from math import pi, cos, sin, sqrt
    1623 
    1624     N = len(s_vec)
    1625     for k in range(N):
    1626         s = s_vec[k]
    1627         phi = phi_vec[k]
    1628 
    1629         # Convert to radians
    1630         phi = phi*pi/180
    1631 
    1632         # Compute velocity vector (u, v)
    1633         u = s*cos(phi)
    1634         v = s*sin(phi)
    1635 
    1636         # Compute wind stress
    1637         S = const * sqrt(u**2 + v**2)
    1638         xmom_update[k] += S*u
    1639         ymom_update[k] += S*v
    1640 
    1641 
    1642 ##
    1643 # @brief A class for a general explicit forcing term.
    1644 class General_forcing:
    1645     """General explicit forcing term for update of quantity
    1646 
    1647     This is used by Inflow and Rainfall for instance
    1648 
    1649 
    1650     General_forcing(quantity_name, rate, center, radius, polygon)
    1651 
    1652     domain:     ANUGA computational domain
    1653     quantity_name: Name of quantity to update.
    1654                    It must be a known conserved quantity.
    1655 
    1656     rate [?/s]: Total rate of change over the specified area.
    1657                 This parameter can be either a constant or a
    1658                 function of time. Positive values indicate increases,
    1659                 negative values indicate decreases.
    1660                 Rate can be None at initialisation but must be specified
    1661                 before forcing term is applied (i.e. simulation has started).
    1662 
    1663     center [m]: Coordinates at center of flow point
    1664     radius [m]: Size of circular area
    1665     polygon:    Arbitrary polygon
    1666     default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
    1667                   Admissible types: None, constant number or function of t
    1668 
    1669 
    1670     Either center, radius or polygon can be specified but not both.
    1671     If neither are specified the entire domain gets updated.
    1672     All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.
    1673 
    1674     Inflow or Rainfall for examples of use
    1675     """
    1676 
    1677 
    1678     # FIXME (AnyOne) : Add various methods to allow spatial variations
    1679 
    1680     ##
    1681     # @brief Create an instance of this forcing term.
    1682     # @param domain
    1683     # @param quantity_name
    1684     # @param rate
    1685     # @param center
    1686     # @param radius
    1687     # @param polygon
    1688     # @param default_rate
    1689     # @param verbose
    1690     def __init__(self,
    1691                  domain,
    1692                  quantity_name,
    1693                  rate=0.0,
    1694                  center=None,
    1695                  radius=None,
    1696                  polygon=None,
    1697                  default_rate=None,
    1698                  verbose=False):
    1699 
    1700         from math import pi, cos, sin
    1701 
    1702         if center is None:
    1703             msg = 'I got radius but no center.'
    1704             assert radius is None, msg
    1705 
    1706         if radius is None:
    1707             msg += 'I got center but no radius.'
    1708             assert center is None, msg
    1709 
    1710         self.domain = domain
    1711         self.quantity_name = quantity_name
    1712         self.rate = rate
    1713         self.center = ensure_numeric(center)
    1714         self.radius = radius
    1715         self.polygon = polygon
    1716         self.verbose = verbose
    1717         self.value = 0.0    # Can be used to remember value at
    1718                             # previous timestep in order to obtain rate
    1719 
    1720         # Get boundary (in absolute coordinates)
    1721         bounding_polygon = domain.get_boundary_polygon()
    1722 
    1723         # Update area if applicable
    1724         if center is not None and radius is not None:
    1725             assert len(center) == 2
    1726             msg = 'Polygon cannot be specified when center and radius are'
    1727             assert polygon is None, msg
    1728 
    1729             # Check that circle center lies within the mesh.
    1730             msg = 'Center %s specified for forcing term did not' % str(center)
    1731             msg += 'fall within the domain boundary.'
    1732             assert is_inside_polygon(center, bounding_polygon), msg
    1733 
    1734             # Check that circle periphery lies within the mesh.
    1735             N = 100
    1736             periphery_points = []
    1737             for i in range(N):
    1738                 theta = 2*pi*i/100
    1739 
    1740                 x = center[0] + radius*cos(theta)
    1741                 y = center[1] + radius*sin(theta)
    1742 
    1743                 periphery_points.append([x,y])
    1744 
    1745             for point in periphery_points:
    1746                 msg = 'Point %s on periphery for forcing term' % str(point)
    1747                 msg += ' did not fall within the domain boundary.'
    1748                 assert is_inside_polygon(point, bounding_polygon), msg
    1749 
    1750         if polygon is not None:
    1751             # Check that polygon lies within the mesh.
    1752             for point in self.polygon:
    1753                 msg = 'Point %s in polygon for forcing term' % str(point)
    1754                 msg += ' did not fall within the domain boundary.'
    1755                 assert is_inside_polygon(point, bounding_polygon), msg
    1756 
    1757         # Pointer to update vector
    1758         self.update = domain.quantities[self.quantity_name].explicit_update
    1759 
    1760         # Determine indices in flow area
    1761         N = len(domain)
    1762         points = domain.get_centroid_coordinates(absolute=True)
    1763 
    1764         # Calculate indices in exchange area for this forcing term
    1765         self.exchange_indices = None
    1766         if self.center is not None and self.radius is not None:
    1767             # Inlet is circular
    1768             inlet_region = 'center=%s, radius=%s' % (self.center, self.radius)
    1769 
    1770             self.exchange_indices = []
    1771             for k in range(N):
    1772                 x, y = points[k,:]    # Centroid
    1773 
    1774                 c = self.center
    1775                 if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
    1776                     self.exchange_indices.append(k)
    1777 
    1778         if self.polygon is not None:
    1779             # Inlet is polygon
    1780             inlet_region = 'polygon=%s' % (self.polygon)
    1781             self.exchange_indices = inside_polygon(points, self.polygon)
    1782 
    1783         if self.exchange_indices is None:
    1784             self.exchange_area = polygon_area(bounding_polygon)
    1785         else:   
    1786             if len(self.exchange_indices) == 0:
    1787                 msg = 'No triangles have been identified in '
    1788                 msg += 'specified region: %s' % inlet_region
    1789                 raise Exception, msg
    1790 
    1791             # Compute exchange area as the sum of areas of triangles identified
    1792             # by circle or polygon
    1793             self.exchange_area = 0.0
    1794             for i in self.exchange_indices:
    1795                 self.exchange_area += domain.areas[i]
    1796            
    1797 
    1798         msg = 'Exchange area in forcing term'
    1799         msg += ' has area = %f' %self.exchange_area
    1800         assert self.exchange_area > 0.0           
    1801            
    1802                
    1803 
    1804            
    1805         # Check and store default_rate
    1806         msg = ('Keyword argument default_rate must be either None '
    1807                'or a function of time.\nI got %s.' % str(default_rate))
    1808         assert (default_rate is None or
    1809                 type(default_rate) in [IntType, FloatType] or
    1810                 callable(default_rate)), msg
    1811 
    1812         if default_rate is not None:
    1813             # If it is a constant, make it a function
    1814             if not callable(default_rate):
    1815                 tmp = default_rate
    1816                 default_rate = lambda t: tmp
    1817 
    1818             # Check that default_rate is a function of one argument
    1819             try:
    1820                 default_rate(0.0)
    1821             except:
    1822                 raise Exception, msg
    1823 
    1824         self.default_rate = default_rate
    1825         self.default_rate_invoked = False    # Flag
    1826 
    1827     ##
    1828     # @brief Execute this instance.
    1829     # @param domain
    1830     def __call__(self, domain):
    1831         """Apply inflow function at time specified in domain, update stage"""
    1832 
    1833         # Call virtual method allowing local modifications
    1834         t = domain.get_time()
    1835         try:
    1836             rate = self.update_rate(t)
    1837         except Modeltime_too_early, e:
    1838             raise Modeltime_too_early, e
    1839         except Modeltime_too_late, e:
    1840             if self.default_rate is None:
    1841                 msg = '%s: ANUGA is trying to run longer than specified data.\n' %str(e)
    1842                 msg += 'You can specify keyword argument default_rate in the '
    1843                 msg += 'forcing function to tell it what to do in the absence of time data.'
    1844                 raise Modeltime_too_late, msg   
    1845             else:
    1846                 # Pass control to default rate function
    1847                 rate = self.default_rate(t)
    1848 
    1849                 if self.default_rate_invoked is False:
    1850                     # Issue warning the first time
    1851                     msg = ('%s\n'
    1852                            'Instead I will use the default rate: %s\n'
    1853                            'Note: Further warnings will be supressed'
    1854                            % (str(e), str(self.default_rate)))
    1855                     warn(msg)
    1856 
    1857                     # FIXME (Ole): Replace this crude flag with
    1858                     # Python's ability to print warnings only once.
    1859                     # See http://docs.python.org/lib/warning-filter.html
    1860                     self.default_rate_invoked = True
    1861 
    1862         if rate is None:
    1863             msg = ('Attribute rate must be specified in General_forcing '
    1864                    'or its descendants before attempting to call it')
    1865             raise Exception, msg
    1866 
    1867         # Now rate is a number
    1868         if self.verbose is True:
    1869             log.critical('Rate of %s at time = %.2f = %f'
    1870                          % (self.quantity_name, domain.get_time(), rate))
    1871 
    1872         if self.exchange_indices is None:
    1873             self.update[:] += rate
    1874         else:
    1875             # Brute force assignment of restricted rate
    1876             for k in self.exchange_indices:
    1877                 self.update[k] += rate
    1878 
    1879     ##
    1880     # @brief Update the internal rate.
    1881     # @param t A callable or scalar used to set the rate.
    1882     # @return The new rate.
    1883     def update_rate(self, t):
    1884         """Virtual method allowing local modifications by writing an
    1885         overriding version in descendant
    1886         """
    1887 
    1888         if callable(self.rate):
    1889             rate = self.rate(t)
    1890         else:
    1891             rate = self.rate
    1892 
    1893         return rate
    1894 
    1895     ##
    1896     # @brief Get values for the specified quantity.
    1897     # @param quantity_name Name of the quantity of interest.
    1898     # @return The value(s) of the quantity.
    1899     # @note If 'quantity_name' is None, use self.quantity_name.
    1900     def get_quantity_values(self, quantity_name=None):
    1901         """Return values for specified quantity restricted to opening
    1902 
    1903         Optionally a quantity name can be specified if values from another
    1904         quantity is sought
    1905         """
    1906 
    1907         if quantity_name is None:
    1908             quantity_name = self.quantity_name
    1909 
    1910         q = self.domain.quantities[quantity_name]
    1911         return q.get_values(location='centroids',
    1912                             indices=self.exchange_indices)
    1913 
    1914     ##
    1915     # @brief Set value for the specified quantity.
    1916     # @param val The value object used to set value.
    1917     # @param quantity_name Name of the quantity of interest.
    1918     # @note If 'quantity_name' is None, use self.quantity_name.
    1919     def set_quantity_values(self, val, quantity_name=None):
    1920         """Set values for specified quantity restricted to opening
    1921 
    1922         Optionally a quantity name can be specified if values from another
    1923         quantity is sought
    1924         """
    1925 
    1926         if quantity_name is None:
    1927             quantity_name = self.quantity_name
    1928 
    1929         q = self.domain.quantities[self.quantity_name]
    1930         q.set_values(val,
    1931                      location='centroids',
    1932                      indices=self.exchange_indices)
    1933 
    1934 
    1935 ##
    1936 # @brief A class for rainfall forcing function.
    1937 # @note Inherits from General_forcing.
    1938 class Rainfall(General_forcing):
    1939     """Class Rainfall - general 'rain over entire domain' forcing term.
    1940 
    1941     Used for implementing Rainfall over the entire domain.
    1942 
    1943         Current Limited to only One Gauge..
    1944 
    1945         Need to add Spatial Varying Capability
    1946         (This module came from copying and amending the Inflow Code)
    1947 
    1948     Rainfall(rain)
    1949 
    1950     domain
    1951     rain [mm/s]:  Total rain rate over the specified domain.
    1952                   NOTE: Raingauge Data needs to reflect the time step.
    1953                   IE: if Gauge is mm read at a time step, then the input
    1954                   here is as mm/(timeStep) so 10mm in 5minutes becomes
    1955                   10/(5x60) = 0.0333mm/s.
    1956 
    1957                   This parameter can be either a constant or a
    1958                   function of time. Positive values indicate inflow,
    1959                   negative values indicate outflow.
    1960                   (and be used for Infiltration - Write Seperate Module)
    1961                   The specified flow will be divided by the area of
    1962                   the inflow region and then applied to update the
    1963                   stage quantity.
    1964 
    1965     polygon: Specifies a polygon to restrict the rainfall.
    1966 
    1967     Examples
    1968     How to put them in a run File...
    1969 
    1970     #------------------------------------------------------------------------
    1971     # Setup specialised forcing terms
    1972     #------------------------------------------------------------------------
    1973     # This is the new element implemented by Ole and Rudy to allow direct
    1974     # input of Rainfall in mm/s
    1975 
    1976     catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms'))
    1977                         # Note need path to File in String.
    1978                         # Else assumed in same directory
    1979 
    1980     domain.forcing_terms.append(catchmentrainfall)
    1981     """
    1982 
    1983     ##
    1984     # @brief Create an instance of the class.
    1985     # @param domain Domain of interest.
    1986     # @param rate Total rain rate over the specified domain (mm/s).
    1987     # @param center
    1988     # @param radius
    1989     # @param polygon Polygon  to restrict rainfall.
    1990     # @param default_rate
    1991     # @param verbose True if this instance is to be verbose.
    1992     def __init__(self,
    1993                  domain,
    1994                  rate=0.0,
    1995                  center=None,
    1996                  radius=None,
    1997                  polygon=None,
    1998                  default_rate=None,
    1999                  verbose=False):
    2000 
    2001         # Converting mm/s to m/s to apply in ANUGA)
    2002         if callable(rate):
    2003             rain = lambda t: rate(t)/1000.0
    2004         else:
    2005             rain = rate/1000.0
    2006 
    2007         if default_rate is not None:
    2008             if callable(default_rate):
    2009                 default_rain = lambda t: default_rate(t)/1000.0
    2010             else:
    2011                 default_rain = default_rate/1000.0
    2012         else:
    2013             default_rain = None
    2014 
    2015            
    2016            
    2017         General_forcing.__init__(self,
    2018                                  domain,
    2019                                  'stage',
    2020                                  rate=rain,
    2021                                  center=center,
    2022                                  radius=radius,
    2023                                  polygon=polygon,
    2024                                  default_rate=default_rain,
    2025                                  verbose=verbose)
    2026 
    2027 
    2028 ##
    2029 # @brief A class for inflow (rain and drain) forcing function.
    2030 # @note Inherits from General_forcing.
    2031 class Inflow(General_forcing):
    2032     """Class Inflow - general 'rain and drain' forcing term.
    2033 
    2034     Useful for implementing flows in and out of the domain.
    2035 
    2036     Inflow(flow, center, radius, polygon)
    2037 
    2038     domain
    2039     rate [m^3/s]: Total flow rate over the specified area.
    2040                   This parameter can be either a constant or a
    2041                   function of time. Positive values indicate inflow,
    2042                   negative values indicate outflow.
    2043                   The specified flow will be divided by the area of
    2044                   the inflow region and then applied to update stage.
    2045     center [m]: Coordinates at center of flow point
    2046     radius [m]: Size of circular area
    2047     polygon:    Arbitrary polygon.
    2048 
    2049     Either center, radius or polygon must be specified
    2050 
    2051     Examples
    2052 
    2053     # Constant drain at 0.003 m^3/s.
    2054     # The outflow area is 0.07**2*pi=0.0154 m^2
    2055     # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
    2056     #
    2057     Inflow((0.7, 0.4), 0.07, -0.003)
    2058 
    2059 
    2060     # Tap turning up to a maximum inflow of 0.0142 m^3/s.
    2061     # The inflow area is 0.03**2*pi = 0.00283 m^2
    2062     # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s
    2063     # over the specified area
    2064     Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
    2065 
    2066 
    2067     #------------------------------------------------------------------------
    2068     # Setup specialised forcing terms
    2069     #------------------------------------------------------------------------
    2070     # This is the new element implemented by Ole to allow direct input
    2071     # of Inflow in m^3/s
    2072 
    2073     hydrograph = Inflow(center=(320, 300), radius=10,
    2074                         rate=file_function('Q/QPMF_Rot_Sub13.tms'))
    2075 
    2076     domain.forcing_terms.append(hydrograph)
    2077     """
    2078 
    2079     ##
    2080     # @brief Create an instance of the class.
    2081     # @param domain Domain of interest.
    2082     # @param rate Total rain rate over the specified domain (mm/s).
    2083     # @param center
    2084     # @param radius
    2085     # @param polygon Polygon  to restrict rainfall.
    2086     # @param default_rate
    2087     # @param verbose True if this instance is to be verbose.
    2088     def __init__(self,
    2089                  domain,
    2090                  rate=0.0,
    2091                  center=None,
    2092                  radius=None,
    2093                  polygon=None,
    2094                  default_rate=None,
    2095                  verbose=False):
    2096         # Create object first to make area is available
    2097         General_forcing.__init__(self,
    2098                                  domain,
    2099                                  'stage',
    2100                                  rate=rate,
    2101                                  center=center,
    2102                                  radius=radius,
    2103                                  polygon=polygon,
    2104                                  default_rate=default_rate,
    2105                                  verbose=verbose)
    2106 
    2107     ##
    2108     # @brief Update the instance rate.
    2109     # @param t New rate object.
    2110     def update_rate(self, t):
    2111         """Virtual method allowing local modifications by writing an
    2112         overriding version in descendant
    2113 
    2114         This one converts m^3/s to m/s which can be added directly
    2115         to 'stage' in ANUGA
    2116         """
    2117 
    2118         if callable(self.rate):
    2119             _rate = self.rate(t)/self.exchange_area
    2120         else:
    2121             _rate = self.rate/self.exchange_area
    2122 
    2123         return _rate
    2124 
    2125 
    2126 ##
    2127 # @brief A class for creating cross sections.
    2128 # @note Inherits from General_forcing.
    2129 class Cross_section:
    2130     """Class Cross_section - a class to setup a cross section from
    2131     which you can then calculate flow and energy through cross section
    2132 
    2133 
    2134     Cross_section(domain, polyline)
    2135 
    2136     domain:
    2137     polyline: Representation of desired cross section - it may contain
    2138               multiple sections allowing for complex shapes. Assume
    2139               absolute UTM coordinates.
    2140               Format [[x0, y0], [x1, y1], ...]
    2141     verbose:
    2142     """
    2143 
    2144     ##
    2145     # @brief Create an instance of the class.
    2146     # @param domain Domain of interest.
    2147     # @param polyline Polyline defining cross section
    2148     # @param verbose True if this instance is to be verbose.
    2149     def __init__(self,
    2150                  domain,
    2151                  polyline=None,
    2152                  verbose=False):
    2153        
    2154         self.domain = domain
    2155         self.polyline = polyline
    2156         self.verbose = verbose
    2157        
    2158         # Find all intersections and associated triangles.
    2159         self.segments = self.domain.get_intersecting_segments(self.polyline,
    2160                                                               use_cache=True,
    2161                                                               verbose=self.verbose)
    2162        
    2163         # Get midpoints
    2164         self.midpoints = segment_midpoints(self.segments)
    2165 
    2166         # Make midpoints Geospatial instances
    2167         self.midpoints = ensure_geospatial(self.midpoints, self.domain.geo_reference)
    2168 
    2169     ##
    2170     # @brief set verbose mode
    2171     def set_verbose(self,verbose=True):
    2172         """Set verbose mode true or flase
    2173         """
    2174 
    2175         self.verbose=verbose
    2176 
    2177     ##
    2178     # @brief calculate current flow through cross section
    2179     def get_flow_through_cross_section(self):
    2180         """ Output: Total flow [m^3/s] across cross section.
    2181         """
    2182 
    2183         # Get interpolated values
    2184         xmomentum = self.domain.get_quantity('xmomentum')
    2185         ymomentum = self.domain.get_quantity('ymomentum')
    2186 
    2187         uh = xmomentum.get_values(interpolation_points=self.midpoints,
    2188                                   use_cache=True)
    2189         vh = ymomentum.get_values(interpolation_points=self.midpoints,
    2190                                   use_cache=True)
    2191 
    2192         # Compute and sum flows across each segment
    2193         total_flow = 0
    2194         for i in range(len(uh)):
    2195             # Inner product of momentum vector with segment normal [m^2/s]
    2196             normal = self.segments[i].normal
    2197             normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]
    2198 
    2199             # Flow across this segment [m^3/s]
    2200             segment_flow = normal_momentum*self.segments[i].length
    2201 
    2202             # Accumulate
    2203             total_flow += segment_flow
    2204 
    2205         return total_flow
    2206  
    2207 
    2208     ##
    2209     # @brief calculate current energy flow through cross section
    2210     def get_energy_through_cross_section(self, kind='total'):
    2211         """Obtain average energy head [m] across specified cross section.
    2212 
    2213         Output:
    2214             E: Average energy [m] across given segments for all stored times.
    2215 
    2216         The average velocity is computed for each triangle intersected by
    2217         the polyline and averaged weighted by segment lengths.
    2218 
    2219         The typical usage of this function would be to get average energy of
    2220         flow in a channel, and the polyline would then be a cross section
    2221         perpendicular to the flow.
    2222 
    2223         #FIXME (Ole) - need name for this energy reflecting that its dimension
    2224         is [m].
    2225         """
    2226 
    2227         from anuga.config import g, epsilon, velocity_protection as h0
    2228        
    2229         # Get interpolated values
    2230         stage = self.domain.get_quantity('stage')
    2231         elevation = self.domain.get_quantity('elevation')
    2232         xmomentum = self.domain.get_quantity('xmomentum')
    2233         ymomentum = self.domain.get_quantity('ymomentum')
    2234 
    2235         w = stage.get_values(interpolation_points=self.midpoints, use_cache=True)
    2236         z = elevation.get_values(interpolation_points=self.midpoints, use_cache=True)
    2237         uh = xmomentum.get_values(interpolation_points=self.midpoints,
    2238                                   use_cache=True)
    2239         vh = ymomentum.get_values(interpolation_points=self.midpoints,
    2240                                   use_cache=True)
    2241         h = w-z                # Depth
    2242 
    2243         # Compute total length of polyline for use with weighted averages
    2244         total_line_length = 0.0
    2245         for segment in self.segments:
    2246             total_line_length += segment.length
    2247 
    2248         # Compute and sum flows across each segment
    2249         average_energy = 0.0
    2250         for i in range(len(w)):
    2251             # Average velocity across this segment
    2252             if h[i] > epsilon:
    2253                 # Use protection against degenerate velocities
    2254                 u = uh[i]/(h[i] + h0/h[i])
    2255                 v = vh[i]/(h[i] + h0/h[i])
    2256             else:
    2257                 u = v = 0.0
    2258 
    2259             speed_squared = u*u + v*v
    2260             kinetic_energy = 0.5*speed_squared/g
    2261 
    2262             if kind == 'specific':
    2263                 segment_energy = h[i] + kinetic_energy
    2264             elif kind == 'total':
    2265                 segment_energy = w[i] + kinetic_energy
    2266             else:
    2267                 msg = 'Energy kind must be either "specific" or "total".'
    2268                 msg += ' I got %s' %kind
    2269 
    2270             # Add to weighted average
    2271             weigth = self.segments[i].length/total_line_length
    2272             average_energy += segment_energy*weigth
    2273 
    2274         return average_energy
    2275 
    2276 
    22771422
    22781423################################################################################
Note: See TracChangeset for help on using the changeset viewer.