Changeset 8126
 Timestamp:
 Mar 7, 2011, 5:09:14 PM (13 years ago)
 Location:
 trunk/anuga_core/source/anuga/shallow_water
 Files:

 1 deleted
 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/shallow_water/test_forcing.py
r8122 r8126 1917 1917 os.remove(field_sts_filename+'.sww') 1918 1918 1919 def test_gravity(self): 1920 #Assuming no friction 1921 1922 from anuga.config import g 1923 1924 a = [0.0, 0.0] 1925 b = [0.0, 2.0] 1926 c = [2.0, 0.0] 1927 d = [0.0, 4.0] 1928 e = [2.0, 2.0] 1929 f = [4.0, 0.0] 1930 1931 points = [a, b, c, d, e, f] 1932 # bac, bce, ecf, dbe 1933 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1934 1935 domain = Domain(points, vertices) 1936 1937 #Set up for a gradient of (3,0) at mid triangle (bce) 1938 def slope(x, y): 1939 return 3*x 1940 1941 h = 0.1 1942 def stage(x, y): 1943 return slope(x, y) + h 1944 1945 domain.set_quantity('elevation', slope) 1946 domain.set_quantity('stage', stage) 1947 1948 for name in domain.conserved_quantities: 1949 assert num.allclose(domain.quantities[name].explicit_update, 0) 1950 assert num.allclose(domain.quantities[name].semi_implicit_update, 0) 1951 1952 domain.compute_forcing_terms() 1953 1954 assert num.allclose(domain.quantities['stage'].explicit_update, 0) 1955 assert num.allclose(domain.quantities['xmomentum'].explicit_update, 1956 g*h*3) 1957 assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0) 1958 1959 1960 1961 def test_manning_friction_old(self): 1962 from anuga.config import g 1963 1964 a = [0.0, 0.0] 1965 b = [0.0, 2.0] 1966 c = [2.0, 0.0] 1967 d = [0.0, 4.0] 1968 e = [2.0, 2.0] 1969 f = [4.0, 0.0] 1970 1971 points = [a, b, c, d, e, f] 1972 # bac, bce, ecf, dbe 1973 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1974 1975 domain = Domain(points, vertices) 1976 1977 #Set up for a gradient of (3,0) at mid triangle (bce) 1978 def slope(x, y): 1979 return 3*x 1980 1981 h = 0.1 1982 def stage(x, y): 1983 return slope(x, y) + h 1984 1985 eta = 0.07 1986 domain.set_quantity('elevation', slope) 1987 domain.set_quantity('stage', stage) 1988 domain.set_quantity('friction', eta) 1989 1990 for name in domain.conserved_quantities: 1991 assert num.allclose(domain.quantities[name].explicit_update, 0) 1992 assert num.allclose(domain.quantities[name].semi_implicit_update, 0) 1993 1994 domain.compute_forcing_terms() 1995 1996 assert num.allclose(domain.quantities['stage'].explicit_update, 0) 1997 assert num.allclose(domain.quantities['xmomentum'].explicit_update, 1998 g*h*3) 1999 assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0) 2000 2001 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2002 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2003 0) 2004 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2005 0) 2006 2007 #Create some momentum for friction to work with 2008 domain.set_quantity('xmomentum', 1) 2009 S = g*eta**2 / h**(7.0/3) 2010 2011 domain.compute_forcing_terms() 2012 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2013 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2014 S) 2015 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2016 0) 2017 2018 #A more complex example 2019 domain.quantities['stage'].semi_implicit_update[:] = 0.0 2020 domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0 2021 domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0 2022 2023 domain.set_quantity('xmomentum', 3) 2024 domain.set_quantity('ymomentum', 4) 2025 # sqrt(3^2 +4^2) = 5 2026 2027 S = g*eta**2 / h**(7.0/3) * 5 2028 2029 domain.compute_forcing_terms() 2030 2031 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2032 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,3*S) 2033 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,4*S) 2034 2035 2036 def test_manning_friction_new(self): 2037 from anuga.config import g 2038 import math 2039 2040 a = [0.0, 0.0] 2041 b = [0.0, 2.0] 2042 c = [2.0, 0.0] 2043 d = [0.0, 4.0] 2044 e = [2.0, 2.0] 2045 f = [4.0, 0.0] 2046 2047 points = [a, b, c, d, e, f] 2048 # bac, bce, ecf, dbe 2049 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2050 2051 domain = Domain(points, vertices) 2052 2053 # Use the new function which takes into account the extra 2054 # wetted area due to slope of bed 2055 domain.set_sloped_mannings_function(True) 2056 2057 #Set up for a gradient of (3,0) at mid triangle (bce) 2058 def slope(x, y): 2059 return 3*x 2060 2061 h = 0.1 2062 def stage(x, y): 2063 return slope(x, y) + h 2064 2065 eta = 0.07 2066 domain.set_quantity('elevation', slope) 2067 domain.set_quantity('stage', stage) 2068 domain.set_quantity('friction', eta) 2069 2070 for name in domain.conserved_quantities: 2071 assert num.allclose(domain.quantities[name].explicit_update, 0) 2072 assert num.allclose(domain.quantities[name].semi_implicit_update, 0) 2073 2074 domain.compute_forcing_terms() 2075 2076 assert num.allclose(domain.quantities['stage'].explicit_update, 0) 2077 assert num.allclose(domain.quantities['xmomentum'].explicit_update, 2078 g*h*3) 2079 assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0) 2080 2081 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2082 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2083 0) 2084 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2085 0) 2086 2087 #Create some momentum for friction to work with 2088 domain.set_quantity('xmomentum', 1) 2089 S = g*eta**2 / h**(7.0/3) * math.sqrt(10) 2090 2091 domain.compute_forcing_terms() 2092 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2093 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2094 S) 2095 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2096 0) 2097 2098 #A more complex example 2099 domain.quantities['stage'].semi_implicit_update[:] = 0.0 2100 domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0 2101 domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0 2102 2103 domain.set_quantity('xmomentum', 3) 2104 domain.set_quantity('ymomentum', 4) 2105 2106 S = g*eta**2 *5 / h**(7.0/3) * math.sqrt(10.0) 2107 2108 domain.compute_forcing_terms() 2109 2110 #print 'S', S 2111 #print domain.quantities['xmomentum'].semi_implicit_update 2112 #print domain.quantities['ymomentum'].semi_implicit_update 2113 2114 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2115 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,3*S) 2116 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,4*S) 2117 2118 2119 2120 2121 2122 def test_inflow_using_circle(self): 2123 from math import pi, cos, sin 2124 2125 a = [0.0, 0.0] 2126 b = [0.0, 2.0] 2127 c = [2.0, 0.0] 2128 d = [0.0, 4.0] 2129 e = [2.0, 2.0] 2130 f = [4.0, 0.0] 2131 2132 points = [a, b, c, d, e, f] 2133 # bac, bce, ecf, dbe 2134 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2135 2136 domain = Domain(points, vertices) 2137 2138 # Flat surface with 1m of water 2139 domain.set_quantity('elevation', 0) 2140 domain.set_quantity('stage', 1.0) 2141 domain.set_quantity('friction', 0) 2142 2143 Br = Reflective_boundary(domain) 2144 domain.set_boundary({'exterior': Br}) 2145 2146 # Setup only one forcing term, constant inflow of 2 m^3/s 2147 # on a circle affecting triangles #0 and #1 (bac and bce) 2148 domain.forcing_terms = [] 2149 2150 I = Inflow(domain, rate=2.0, center=(1,1), radius=1) 2151 domain.forcing_terms.append(I) 2152 domain.compute_forcing_terms() 2153 2154 2155 A = I.exchange_area 2156 assert num.allclose(A, 4) # Two triangles 2157 2158 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A) 2159 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A) 2160 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2161 2162 2163 def test_inflow_using_circle_function(self): 2164 from math import pi, cos, sin 2165 2166 a = [0.0, 0.0] 2167 b = [0.0, 2.0] 2168 c = [2.0, 0.0] 2169 d = [0.0, 4.0] 2170 e = [2.0, 2.0] 2171 f = [4.0, 0.0] 2172 2173 points = [a, b, c, d, e, f] 2174 # bac, bce, ecf, dbe 2175 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2176 2177 domain = Domain(points, vertices) 2178 2179 # Flat surface with 1m of water 2180 domain.set_quantity('elevation', 0) 2181 domain.set_quantity('stage', 1.0) 2182 domain.set_quantity('friction', 0) 2183 2184 Br = Reflective_boundary(domain) 2185 domain.set_boundary({'exterior': Br}) 2186 2187 # Setup only one forcing term, time dependent inflow of 2 m^3/s 2188 # on a circle affecting triangles #0 and #1 (bac and bce) 2189 domain.forcing_terms = [] 2190 I = Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) 2191 domain.forcing_terms.append(I) 2192 2193 domain.compute_forcing_terms() 2194 2195 A = I.exchange_area 2196 assert num.allclose(A, 4) # Two triangles 2197 2198 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A) 2199 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A) 2200 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2201 2202 2203 2204 2205 def test_inflow_catch_too_few_triangles(self): 2206 """ 2207 Test that exception is thrown if no triangles are covered 2208 by the inflow area 2209 """ 2210 2211 from math import pi, cos, sin 2212 2213 a = [0.0, 0.0] 2214 b = [0.0, 2.0] 2215 c = [2.0, 0.0] 2216 d = [0.0, 4.0] 2217 e = [2.0, 2.0] 2218 f = [4.0, 0.0] 2219 2220 points = [a, b, c, d, e, f] 2221 # bac, bce, ecf, dbe 2222 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2223 2224 domain = Domain(points, vertices) 2225 2226 # Flat surface with 1m of water 2227 domain.set_quantity('elevation', 0) 2228 domain.set_quantity('stage', 1.0) 2229 domain.set_quantity('friction', 0) 2230 2231 Br = Reflective_boundary(domain) 2232 domain.set_boundary({'exterior': Br}) 2233 2234 # Setup only one forcing term, constant inflow of 2 m^3/s 2235 # on a circle affecting triangles #0 and #1 (bac and bce) 2236 try: 2237 Inflow(domain, rate=2.0, center=(1,1.1), radius=0.01) 2238 except: 2239 pass 2240 else: 2241 msg = 'Should have raised exception' 2242 raise Exception, msg 2243 1919 2244 1920 2245
Note: See TracChangeset
for help on using the changeset viewer.