- Timestamp:
- May 19, 2010, 10:13:06 AM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7731 r7733 83 83 import numpy as num 84 84 85 from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints86 85 from 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 87 from anuga.shallow_water.forcing import Cross_section 100 88 from anuga.pmesh.mesh_interface import create_mesh_from_regions 101 89 from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric … … 1432 1420 1433 1421 1434 ################################################################################1435 # Experimental auxiliary functions1436 ################################################################################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 vectors1447 and that it returns an array or a list of same length1448 as x and y1449 2: a scalar1450 """1451 1452 if callable(f):1453 N = 31454 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 semantics1461 raise Exception, msg1462 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, msg1470 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_info1476 try:1477 result_len = len(q)1478 except:1479 msg = '%s must return vector' % func_msg1480 self.fail(msg)1481 msg = '%s must return vector of length %d' % (func_msg, N)1482 assert result_len == N, msg1483 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, msg1490 1491 return f1492 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 of1498 wind speed [m/s] and wind direction [degrees]1499 """1500 1501 ##1502 # @brief Create an instance of Wind_stress.1503 # @param *args1504 # @param **kwargs1505 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. the1514 vector (1,0) has zero degrees.1515 We may need to convert from 'compass' degrees later on and also1516 map from True north to grid north.1517 1518 Arguments can also be Python functions of t,x,y as in1519 1520 def speed(t,x,y):1521 ...1522 return s1523 1524 def angle(t,x,y):1525 ...1526 return phi1527 1528 where x and y are vectors.1529 1530 and then pass the functions in1531 1532 W = Wind_stress(speed, angle)1533 1534 The instantiated object W can be appended to the list of1535 forcing_terms as in1536 1537 Alternatively, one vector valued function for (speed, angle)1538 can be applied, providing both quantities simultaneously.1539 As in1540 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_w1546 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 arguments1557 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_w1567 1568 ##1569 # @brief 'execute' this class instance.1570 # @param domain1571 def __call__(self, domain):1572 """Evaluate windfield based on values found in domain"""1573 1574 from math import pi, cos, sin, sqrt1575 1576 xmom_update = domain.quantities['xmomentum'].explicit_update1577 ymom_update = domain.quantities['ymomentum'].explicit_update1578 1579 N = len(domain) # number_of_triangles1580 t = domain.time1581 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 scalar1587 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.s1591 raise msg1592 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 scalar1598 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.phi1603 raise msg1604 1605 assign_windfield_values(xmom_update, ymom_update,1606 s_vec, phi_vec, self.const)1607 1608 1609 ##1610 # @brief Assign wind field values1611 # @param xmom_update1612 # @param ymom_update1613 # @param s_vec1614 # @param phi_vec1615 # @param const1616 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, sqrt1623 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 radians1630 phi = phi*pi/1801631 1632 # Compute velocity vector (u, v)1633 u = s*cos(phi)1634 v = s*sin(phi)1635 1636 # Compute wind stress1637 S = const * sqrt(u**2 + v**2)1638 xmom_update[k] += S*u1639 ymom_update[k] += S*v1640 1641 1642 ##1643 # @brief A class for a general explicit forcing term.1644 class General_forcing:1645 """General explicit forcing term for update of quantity1646 1647 This is used by Inflow and Rainfall for instance1648 1649 1650 General_forcing(quantity_name, rate, center, radius, polygon)1651 1652 domain: ANUGA computational domain1653 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 a1658 function of time. Positive values indicate increases,1659 negative values indicate decreases.1660 Rate can be None at initialisation but must be specified1661 before forcing term is applied (i.e. simulation has started).1662 1663 center [m]: Coordinates at center of flow point1664 radius [m]: Size of circular area1665 polygon: Arbitrary polygon1666 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 t1668 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 use1675 """1676 1677 1678 # FIXME (AnyOne) : Add various methods to allow spatial variations1679 1680 ##1681 # @brief Create an instance of this forcing term.1682 # @param domain1683 # @param quantity_name1684 # @param rate1685 # @param center1686 # @param radius1687 # @param polygon1688 # @param default_rate1689 # @param verbose1690 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, sin1701 1702 if center is None:1703 msg = 'I got radius but no center.'1704 assert radius is None, msg1705 1706 if radius is None:1707 msg += 'I got center but no radius.'1708 assert center is None, msg1709 1710 self.domain = domain1711 self.quantity_name = quantity_name1712 self.rate = rate1713 self.center = ensure_numeric(center)1714 self.radius = radius1715 self.polygon = polygon1716 self.verbose = verbose1717 self.value = 0.0 # Can be used to remember value at1718 # previous timestep in order to obtain rate1719 1720 # Get boundary (in absolute coordinates)1721 bounding_polygon = domain.get_boundary_polygon()1722 1723 # Update area if applicable1724 if center is not None and radius is not None:1725 assert len(center) == 21726 msg = 'Polygon cannot be specified when center and radius are'1727 assert polygon is None, msg1728 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), msg1733 1734 # Check that circle periphery lies within the mesh.1735 N = 1001736 periphery_points = []1737 for i in range(N):1738 theta = 2*pi*i/1001739 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), msg1749 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), msg1756 1757 # Pointer to update vector1758 self.update = domain.quantities[self.quantity_name].explicit_update1759 1760 # Determine indices in flow area1761 N = len(domain)1762 points = domain.get_centroid_coordinates(absolute=True)1763 1764 # Calculate indices in exchange area for this forcing term1765 self.exchange_indices = None1766 if self.center is not None and self.radius is not None:1767 # Inlet is circular1768 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,:] # Centroid1773 1774 c = self.center1775 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 polygon1780 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_region1789 raise Exception, msg1790 1791 # Compute exchange area as the sum of areas of triangles identified1792 # by circle or polygon1793 self.exchange_area = 0.01794 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_area1800 assert self.exchange_area > 0.01801 1802 1803 1804 1805 # Check and store default_rate1806 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 or1809 type(default_rate) in [IntType, FloatType] or1810 callable(default_rate)), msg1811 1812 if default_rate is not None:1813 # If it is a constant, make it a function1814 if not callable(default_rate):1815 tmp = default_rate1816 default_rate = lambda t: tmp1817 1818 # Check that default_rate is a function of one argument1819 try:1820 default_rate(0.0)1821 except:1822 raise Exception, msg1823 1824 self.default_rate = default_rate1825 self.default_rate_invoked = False # Flag1826 1827 ##1828 # @brief Execute this instance.1829 # @param domain1830 def __call__(self, domain):1831 """Apply inflow function at time specified in domain, update stage"""1832 1833 # Call virtual method allowing local modifications1834 t = domain.get_time()1835 try:1836 rate = self.update_rate(t)1837 except Modeltime_too_early, e:1838 raise Modeltime_too_early, e1839 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, msg1845 else:1846 # Pass control to default rate function1847 rate = self.default_rate(t)1848 1849 if self.default_rate_invoked is False:1850 # Issue warning the first time1851 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 with1858 # Python's ability to print warnings only once.1859 # See http://docs.python.org/lib/warning-filter.html1860 self.default_rate_invoked = True1861 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, msg1866 1867 # Now rate is a number1868 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[:] += rate1874 else:1875 # Brute force assignment of restricted rate1876 for k in self.exchange_indices:1877 self.update[k] += rate1878 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 an1885 overriding version in descendant1886 """1887 1888 if callable(self.rate):1889 rate = self.rate(t)1890 else:1891 rate = self.rate1892 1893 return rate1894 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 opening1902 1903 Optionally a quantity name can be specified if values from another1904 quantity is sought1905 """1906 1907 if quantity_name is None:1908 quantity_name = self.quantity_name1909 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 opening1921 1922 Optionally a quantity name can be specified if values from another1923 quantity is sought1924 """1925 1926 if quantity_name is None:1927 quantity_name = self.quantity_name1928 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 Capability1946 (This module came from copying and amending the Inflow Code)1947 1948 Rainfall(rain)1949 1950 domain1951 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 input1954 here is as mm/(timeStep) so 10mm in 5minutes becomes1955 10/(5x60) = 0.0333mm/s.1956 1957 This parameter can be either a constant or a1958 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 of1962 the inflow region and then applied to update the1963 stage quantity.1964 1965 polygon: Specifies a polygon to restrict the rainfall.1966 1967 Examples1968 How to put them in a run File...1969 1970 #------------------------------------------------------------------------1971 # Setup specialised forcing terms1972 #------------------------------------------------------------------------1973 # This is the new element implemented by Ole and Rudy to allow direct1974 # input of Rainfall in mm/s1975 1976 catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms'))1977 # Note need path to File in String.1978 # Else assumed in same directory1979 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 center1988 # @param radius1989 # @param polygon Polygon to restrict rainfall.1990 # @param default_rate1991 # @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.02004 else:2005 rain = rate/1000.02006 2007 if default_rate is not None:2008 if callable(default_rate):2009 default_rain = lambda t: default_rate(t)/1000.02010 else:2011 default_rain = default_rate/1000.02012 else:2013 default_rain = None2014 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 domain2039 rate [m^3/s]: Total flow rate over the specified area.2040 This parameter can be either a constant or a2041 function of time. Positive values indicate inflow,2042 negative values indicate outflow.2043 The specified flow will be divided by the area of2044 the inflow region and then applied to update stage.2045 center [m]: Coordinates at center of flow point2046 radius [m]: Size of circular area2047 polygon: Arbitrary polygon.2048 2049 Either center, radius or polygon must be specified2050 2051 Examples2052 2053 # Constant drain at 0.003 m^3/s.2054 # The outflow area is 0.07**2*pi=0.0154 m^22055 # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s2056 #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^22062 # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s2063 # over the specified area2064 Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))2065 2066 2067 #------------------------------------------------------------------------2068 # Setup specialised forcing terms2069 #------------------------------------------------------------------------2070 # This is the new element implemented by Ole to allow direct input2071 # of Inflow in m^3/s2072 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 center2084 # @param radius2085 # @param polygon Polygon to restrict rainfall.2086 # @param default_rate2087 # @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 available2097 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 an2112 overriding version in descendant2113 2114 This one converts m^3/s to m/s which can be added directly2115 to 'stage' in ANUGA2116 """2117 2118 if callable(self.rate):2119 _rate = self.rate(t)/self.exchange_area2120 else:2121 _rate = self.rate/self.exchange_area2122 2123 return _rate2124 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 from2131 which you can then calculate flow and energy through cross section2132 2133 2134 Cross_section(domain, polyline)2135 2136 domain:2137 polyline: Representation of desired cross section - it may contain2138 multiple sections allowing for complex shapes. Assume2139 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 section2148 # @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 = domain2155 self.polyline = polyline2156 self.verbose = verbose2157 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 midpoints2164 self.midpoints = segment_midpoints(self.segments)2165 2166 # Make midpoints Geospatial instances2167 self.midpoints = ensure_geospatial(self.midpoints, self.domain.geo_reference)2168 2169 ##2170 # @brief set verbose mode2171 def set_verbose(self,verbose=True):2172 """Set verbose mode true or flase2173 """2174 2175 self.verbose=verbose2176 2177 ##2178 # @brief calculate current flow through cross section2179 def get_flow_through_cross_section(self):2180 """ Output: Total flow [m^3/s] across cross section.2181 """2182 2183 # Get interpolated values2184 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 segment2193 total_flow = 02194 for i in range(len(uh)):2195 # Inner product of momentum vector with segment normal [m^2/s]2196 normal = self.segments[i].normal2197 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].length2201 2202 # Accumulate2203 total_flow += segment_flow2204 2205 return total_flow2206 2207 2208 ##2209 # @brief calculate current energy flow through cross section2210 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 by2217 the polyline and averaged weighted by segment lengths.2218 2219 The typical usage of this function would be to get average energy of2220 flow in a channel, and the polyline would then be a cross section2221 perpendicular to the flow.2222 2223 #FIXME (Ole) - need name for this energy reflecting that its dimension2224 is [m].2225 """2226 2227 from anuga.config import g, epsilon, velocity_protection as h02228 2229 # Get interpolated values2230 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 # Depth2242 2243 # Compute total length of polyline for use with weighted averages2244 total_line_length = 0.02245 for segment in self.segments:2246 total_line_length += segment.length2247 2248 # Compute and sum flows across each segment2249 average_energy = 0.02250 for i in range(len(w)):2251 # Average velocity across this segment2252 if h[i] > epsilon:2253 # Use protection against degenerate velocities2254 u = uh[i]/(h[i] + h0/h[i])2255 v = vh[i]/(h[i] + h0/h[i])2256 else:2257 u = v = 0.02258 2259 speed_squared = u*u + v*v2260 kinetic_energy = 0.5*speed_squared/g2261 2262 if kind == 'specific':2263 segment_energy = h[i] + kinetic_energy2264 elif kind == 'total':2265 segment_energy = w[i] + kinetic_energy2266 else:2267 msg = 'Energy kind must be either "specific" or "total".'2268 msg += ' I got %s' %kind2269 2270 # Add to weighted average2271 weigth = self.segments[i].length/total_line_length2272 average_energy += segment_energy*weigth2273 2274 return average_energy2275 2276 2277 1422 2278 1423 ################################################################################
Note: See TracChangeset
for help on using the changeset viewer.
