Changeset 6451
- Timestamp:
- Mar 4, 2009, 3:04:42 PM (16 years ago)
- Location:
- branches/numpy/anuga/shallow_water
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/shallow_water/shallow_water_domain.py
r6441 r6451 1602 1602 msg = 'Function %s could not be executed:\n%s' %(f, e) 1603 1603 # FIXME: Reconsider this semantics 1604 raise msg1604 raise Exception, msg 1605 1605 1606 1606 try: 1607 1607 q = num.array(q, num.float) 1608 1608 except: 1609 msg = 'Return value from vector function %s could ' %f1610 msg += 'not be converted into a numeric array of floats.\n'1611 msg += 'Specified function should return either list or array.'1609 msg = ('Return value from vector function %s could not ' 1610 'be converted into a numeric array of floats.\nSpecified ' 1611 'function should return either list or array.' % f) 1612 1612 raise Exception, msg 1613 1613 … … 1616 1616 func_info = (f.func_name, f.func_code.co_filename, 1617 1617 f.func_code.co_firstlineno) 1618 msg = ('Function %s() must return vector (defined in %s, line %d)' 1619 % func_info) 1620 assert hasattr(q, 'len'), msg 1621 1622 msg = ('Return vector from function %s() must have same ' 1623 'length as input vectors\nq=%s' % (f.func_name, str(q))) 1624 assert len(q) == N, msg 1618 func_msg = 'Function %s (defined in %s, line %d)' % func_info 1619 try: 1620 result_len = len(q) 1621 except: 1622 msg = '%s must return vector' % func_msg 1623 self.fail(msg) 1624 msg = '%s must return vector of length %d' % (func_msg, N) 1625 assert result_len == N, msg 1625 1626 else: 1626 1627 try: -
branches/numpy/anuga/shallow_water/test_shallow_water_domain.py
r6441 r6451 7 7 from anuga.config import g, epsilon 8 8 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a 9 import numpy as num10 9 from anuga.utilities.numerical_tools import mean 11 10 from anuga.utilities.polygon import is_inside_polygon … … 14 13 from anuga.geospatial_data.geospatial_data import Geospatial_data 15 14 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 15 16 16 from shallow_water_domain import * 17 18 import numpy as num 17 19 18 20 # Get gateway to C implementation of flux function for direct testing 19 21 from shallow_water_ext import flux_function_central as flux_function 22 20 23 21 24 # For test_fitting_using_shallow_water_domain example 22 25 def linear_function(point): 23 26 point = num.array(point) 24 return point[:,0]+point[:,1] 27 return point[:,0] + point[:,1] 28 25 29 26 30 class Weir: 27 31 """Set a bathymetry for weir with a hole and a downstream gutter 32 28 33 x,y are assumed to be in the unit square 29 34 """ … … 33 38 34 39 def __call__(self, x, y): 35 36 40 N = len(x) 37 41 assert N == len(y) … … 39 43 z = num.zeros(N, num.float) 40 44 for i in range(N): 41 z[i] = -x[i] /2 #General slope42 43 # Flattish bit to the left45 z[i] = -x[i] / 2 # General slope 46 47 # Flattish bit to the left 44 48 if x[i] < 0.3: 45 49 z[i] = -x[i]/10 46 50 47 # Weir51 # Weir 48 52 if x[i] >= 0.3 and x[i] < 0.4: 49 53 z[i] = -x[i]+0.9 50 54 51 # Dip55 # Dip 52 56 x0 = 0.6 53 #depth = -1.354 57 depth = -1.0 55 #plateaux = -0.956 58 plateaux = -0.6 57 59 if y[i] < 0.7: 58 60 if x[i] > x0 and x[i] < 0.9: 59 61 z[i] = depth 60 61 #RHS plateaux 62 # RHS plateaux 62 63 if x[i] >= 0.9: 63 64 z[i] = plateaux 64 65 66 65 elif y[i] >= 0.7 and y[i] < 1.5: 67 # Restrict and deepen66 # Restrict and deepen 68 67 if x[i] >= x0 and x[i] < 0.8: 69 z[i] = depth-(y[i]/3-0.3) 70 #z[i] = depth-y[i]/5 71 #z[i] = depth 68 z[i] = depth - (y[i]/3 - 0.3) 72 69 elif x[i] >= 0.8: 73 # RHS plateaux70 # RHS plateaux 74 71 z[i] = plateaux 75 76 72 elif y[i] >= 1.5: 77 73 if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2: 78 # Widen up and stay at constant depth79 z[i] = depth -1.5/580 elif x[i] >= 0.8 + (y[i] -1.5)/1.2:81 # RHS plateaux74 # Widen up and stay at constant depth 75 z[i] = depth - 1.5/5 76 elif x[i] >= 0.8 + (y[i] - 1.5)/1.2: 77 # RHS plateaux 82 78 z[i] = plateaux 83 79 84 85 #Hole in weir (slightly higher than inflow condition) 80 # Hole in weir (slightly higher than inflow condition) 86 81 if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4: 87 z[i] = -x[i] +self.inflow_stage + 0.0288 89 # Channel behind weir82 z[i] = -x[i] + self.inflow_stage + 0.02 83 84 # Channel behind weir 90 85 x0 = 0.5 91 86 if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4: 92 z[i] = -x[i] +self.inflow_stage + 0.0287 z[i] = -x[i] + self.inflow_stage + 0.02 93 88 94 89 if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4: 95 # Flatten it out towards the end96 z[i] = -x0 +self.inflow_stage + 0.02 + (x0-x[i])/597 98 # Hole to the east99 x0 = 1.1 ; y0 = 0.35100 #if x[i] < -0.2 and y < 0.5:90 # Flatten it out towards the end 91 z[i] = -x0 + self.inflow_stage + 0.02 + (x0 - x[i])/5 92 93 # Hole to the east 94 x0 = 1.1 95 y0 = 0.35 101 96 if num.sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2: 102 z[i] = num.sqrt(( (x[i]-x0))**2 + ((y[i]-y0))**2)-1.0103 104 # Tiny channel draining hole97 z[i] = num.sqrt((x[i]-x0)**2 + (y[i]-y0)**2) - 1.0 98 99 # Tiny channel draining hole 105 100 if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6: 106 z[i] = -0.9 #North south101 z[i] = -0.9 # North south 107 102 108 103 if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65: 109 z[i] = -1.0 + (x[i]-0.9)/3 #East west 110 111 112 113 #Stuff not in use 114 115 #Upward slope at inlet to the north west 116 #if x[i] < 0.0: # and y[i] > 0.5: 104 z[i] = -1.0 + (x[i]-0.9)/3 # East west 105 106 # Stuff not in use 107 108 # Upward slope at inlet to the north west 109 # if x[i] < 0.0: # and y[i] > 0.5: 117 110 # #z[i] = -y[i]+0.5 #-x[i]/2 118 111 # z[i] = x[i]/4 - y[i]**2 + 0.5 119 112 120 # Hole to the west121 # x0 = -0.4; y0 = 0.35 # center122 # if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:113 # Hole to the west 114 # x0 = -0.4; y0 = 0.35 # center 115 # if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2: 123 116 # z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2 124 117 125 126 127 128 129 118 return z/2 119 130 120 131 121 class Weir_simple: 132 122 """Set a bathymetry for weir with a hole and a downstream gutter 123 133 124 x,y are assumed to be in the unit square 134 125 """ … … 138 129 139 130 def __call__(self, x, y): 140 141 131 N = len(x) 142 132 assert N == len(y) … … 144 134 z = num.zeros(N, num.float) 145 135 for i in range(N): 146 z[i] = -x[i] #General slope147 148 # Flat bit to the left136 z[i] = -x[i] # General slope 137 138 # Flat bit to the left 149 139 if x[i] < 0.3: 150 z[i] = -x[i]/10 #General slope151 152 # Weir140 z[i] = -x[i]/10 # General slope 141 142 # Weir 153 143 if x[i] > 0.3 and x[i] < 0.4: 154 z[i] = -x[i] +0.9155 156 # Dip144 z[i] = -x[i] + 0.9 145 146 # Dip 157 147 if x[i] > 0.6 and x[i] < 0.9: 158 z[i] = -x[i] -0.5 #-y[i]/5159 160 # Hole in weir (slightly higher than inflow condition)148 z[i] = -x[i] - 0.5 # -y[i]/5 149 150 # Hole in weir (slightly higher than inflow condition) 161 151 if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4: 162 z[i] = -x[i]+self.inflow_stage + 0.05 163 152 z[i] = -x[i] + self.inflow_stage + 0.05 164 153 165 154 return z/2 166 155 167 156 168 169 170 #Variable windfield implemented using functions 171 def speed(t,x,y): 157 # Variable windfield implemented using functions 158 def speed(t, x, y): 172 159 """Large speeds halfway between center and edges 160 173 161 Low speeds at center and edges 174 162 """ … … 180 168 181 169 N = len(x) 182 s = 0 *x #New array170 s = 0 * x # New array 183 171 184 172 for k in range(N): 185 186 173 r = num.sqrt(x[k]**2 + y[k]**2) 187 188 factor = exp( -(r-0.15)**2 ) 189 174 factor = exp(-(r-0.15)**2) 190 175 s[k] = 4000 * factor * (cos(t*2*pi/150) + 2) 191 176 192 177 return s 193 178 194 195 def scalar_func(t,x,y): 179 def scalar_func(t, x, y): 196 180 """Function that returns a scalar. 181 197 182 Used to test error message when numeric array is expected 198 183 """ … … 200 185 return 17.7 201 186 202 203 def scalar_func_list(t,x,y): 187 def scalar_func_list(t, x, y): 204 188 """Function that returns a scalar. 189 205 190 Used to test error message when numeric array is expected 206 191 """ … … 208 193 return [17.7] 209 194 210 211 def angle(t,x,y): 195 def angle(t, x, y): 212 196 """Rotating field 213 197 """ … … 218 202 219 203 N = len(x) 220 a = 0 *x #New array204 a = 0 * x # New array 221 205 222 206 for k in range(N): … … 226 210 227 211 if x[k] < 0: 228 angle +=pi #atan in ]-pi/2; pi/2[229 230 # Take normal direction212 angle += pi 213 214 # Take normal direction 231 215 angle -= pi/2 232 216 233 # Ensure positive radians217 # Ensure positive radians 234 218 if angle < 0: 235 219 angle += 2*pi … … 248 232 249 233 def test_rotate(self): 250 normal = num.array([0.0, -1.0])251 252 q = num.array([1.0, 2.0,3.0])234 normal = num.array([0.0, -1.0]) 235 236 q = num.array([1.0, 2.0, 3.0]) 253 237 254 238 r = rotate(q, normal, direction = 1) … … 260 244 assert num.allclose(w, q) 261 245 262 # Check error check246 # Check error check 263 247 try: 264 rotate(r, num.array([1, 1,1]))248 rotate(r, num.array([1, 1, 1])) 265 249 except: 266 250 pass 267 251 else: 268 raise 'Should have raised an exception' 269 252 raise Exception, 'Should have raised an exception' 270 253 271 254 # Individual flux tests 272 255 def test_flux_zero_case(self): 273 ql = num.zeros( 3, num.float)274 qr = num.zeros( 3, num.float)275 normal = num.zeros( 2, num.float)276 edgeflux = num.zeros( 3, num.float)256 ql = num.zeros(3, num.float) 257 qr = num.zeros(3, num.float) 258 normal = num.zeros(2, num.float) 259 edgeflux = num.zeros(3, num.float) 277 260 zl = zr = 0. 278 261 H0 = 0.0 279 280 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) 281 282 assert num.allclose(edgeflux, [0,0,0]) 262 263 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, 264 epsilon, g, H0) 265 266 assert num.allclose(edgeflux, [0, 0, 0]) 283 267 assert max_speed == 0. 284 268 … … 286 270 w = 2.0 287 271 288 normal = num.array([1., 0])272 normal = num.array([1., 0]) 289 273 ql = num.array([w, 0, 0]) 290 274 qr = num.array([w, 0, 0]) 291 edgeflux = num.zeros(3, num.float) 275 edgeflux = num.zeros(3, num.float) 292 276 zl = zr = 0. 293 277 h = w - (zl+zr)/2 294 278 H0 = 0.0 295 279 296 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) 280 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, 281 epsilon, g, H0) 282 297 283 assert num.allclose(edgeflux, [0., 0.5*g*h**2, 0.]) 298 284 assert max_speed == num.sqrt(g*h) … … 302 288 # w = 2.0 303 289 # 304 # normal = array([1., 0])290 # normal = array([1., 0]) 305 291 # ql = array([w, 0, 0]) 306 292 # qr = array([w, 0, 0]) … … 313 299 # assert max_speed == sqrt(g*h) 314 300 315 316 301 def test_flux1(self): 317 # Use data from previous version of abstract_2d_finite_volumes318 normal = num.array([1., 0])302 # Use data from previous version of abstract_2d_finite_volumes 303 normal = num.array([1., 0]) 319 304 ql = num.array([-0.2, 2, 3]) 320 305 qr = num.array([-0.2, 2, 3]) 321 306 zl = zr = -0.5 322 edgeflux = num.zeros(3, num.float) 307 edgeflux = num.zeros(3, num.float) 323 308 324 309 H0 = 0.0 325 310 326 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) 327 328 assert num.allclose(edgeflux, [2.,13.77433333, 20.]) 311 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, 312 epsilon, g, H0) 313 314 assert num.allclose(edgeflux, [2., 13.77433333, 20.]) 329 315 assert num.allclose(max_speed, 8.38130948661) 330 316 331 332 317 def test_flux2(self): 333 # Use data from previous version of abstract_2d_finite_volumes318 # Use data from previous version of abstract_2d_finite_volumes 334 319 normal = num.array([0., -1.]) 335 320 ql = num.array([-0.075, 2, 3]) … … 337 322 zl = zr = -0.375 338 323 339 edgeflux = num.zeros(3, num.float) 324 edgeflux = num.zeros(3, num.float) 340 325 H0 = 0.0 341 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) 342 343 assert num.allclose(edgeflux, [-3.,-20.0, -30.441]) 326 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, 327 epsilon, g, H0) 328 329 assert num.allclose(edgeflux, [-3., -20.0, -30.441]) 344 330 assert num.allclose(max_speed, 11.7146428199) 345 331 346 332 def test_flux3(self): 347 # Use data from previous version of abstract_2d_finite_volumes333 # Use data from previous version of abstract_2d_finite_volumes 348 334 normal = num.array([-sqrt(2)/2, sqrt(2)/2]) 349 335 ql = num.array([-0.075, 2, 3]) … … 351 337 zl = zr = -0.375 352 338 353 edgeflux = num.zeros(3, num.float) 339 edgeflux = num.zeros(3, num.float) 354 340 H0 = 0.0 355 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) 341 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, 342 epsilon, g, H0) 356 343 357 344 assert num.allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019]) … … 359 346 360 347 def test_flux4(self): 361 # Use data from previous version of abstract_2d_finite_volumes348 # Use data from previous version of abstract_2d_finite_volumes 362 349 normal = num.array([-sqrt(2)/2, sqrt(2)/2]) 363 350 ql = num.array([-0.34319278, 0.10254161, 0.07273855]) … … 365 352 zl = zr = -0.375 366 353 367 edgeflux = num.zeros(3, num.float) 354 edgeflux = num.zeros(3, num.float) 368 355 H0 = 0.0 369 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) 356 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, 357 epsilon, g, H0) 370 358 371 359 assert num.allclose(edgeflux, [-0.04072676, -0.07096636, -0.01604364]) 372 360 assert num.allclose(max_speed, 1.31414103233) 373 361 374 def test_flux_computation(self): 375 """test_flux_computation - test flux calculation (actual C implementation) 376 This one tests the constant case where only the pressure term contributes to each edge and cancels out 377 once the total flux has been summed up. 362 def test_flux_computation(self): 363 """test flux calculation (actual C implementation) 364 365 This one tests the constant case where only the pressure term 366 contributes to each edge and cancels out once the total flux has 367 been summed up. 378 368 """ 379 369 380 370 a = [0.0, 0.0] 381 371 b = [0.0, 2.0] 382 c = [2.0, 0.0]372 c = [2.0, 0.0] 383 373 d = [0.0, 4.0] 384 374 e = [2.0, 2.0] 385 f = [4.0, 0.0]375 f = [4.0, 0.0] 386 376 387 377 points = [a, b, c, d, e, f] 388 # bac, bce, ecf, dbe, daf, dae378 # bac, bce, ecf, dbe 389 379 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 390 380 … … 392 382 domain.check_integrity() 393 383 394 # The constant case 384 # The constant case 395 385 domain.set_quantity('elevation', -1) 396 domain.set_quantity('stage', 1) 397 386 domain.set_quantity('stage', 1) 387 398 388 domain.compute_fluxes() 399 assert num.allclose(domain.get_quantity('stage').explicit_update[1], 0)# Central triangle400 401 402 # The more general case 403 def surface(x, y):404 return -x/2 405 389 # Central triangle 390 assert num.allclose(domain.get_quantity('stage').explicit_update[1], 0) 391 392 # The more general case 393 def surface(x, y): 394 return -x/2 395 406 396 domain.set_quantity('elevation', -10) 407 domain.set_quantity('stage', surface) 408 domain.set_quantity('xmomentum', 1) 409 397 domain.set_quantity('stage', surface) 398 domain.set_quantity('xmomentum', 1) 399 410 400 domain.compute_fluxes() 411 401 412 402 #print domain.get_quantity('stage').explicit_update 413 403 # FIXME (Ole): TODO the general case 414 #assert allclose(domain.get_quantity('stage').explicit_update[1], ........??) 415 416 417 404 #assert allclose(domain.get_quantity('stage').explicit_update[1], ...??) 405 418 406 def test_sw_domain_simple(self): 419 407 a = [0.0, 0.0] 420 408 b = [0.0, 2.0] 421 c = [2.0, 0.0]409 c = [2.0, 0.0] 422 410 d = [0.0, 4.0] 423 411 e = [2.0, 2.0] 424 f = [4.0, 0.0]412 f = [4.0, 0.0] 425 413 426 414 points = [a, b, c, d, e, f] 427 #bac, bce, ecf, dbe, daf, dae 428 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 429 415 # bac, bce, ecf, dbe 416 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 430 417 431 418 #from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain … … 443 430 assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.) 444 431 445 446 432 def test_boundary_conditions(self): 447 448 433 a = [0.0, 0.0] 449 434 b = [0.0, 2.0] 450 c = [2.0, 0.0]435 c = [2.0, 0.0] 451 436 d = [0.0, 4.0] 452 437 e = [2.0, 2.0] 453 f = [4.0, 0.0]438 f = [4.0, 0.0] 454 439 455 440 points = [a, b, c, d, e, f] 456 #bac, bce, ecf, dbe 457 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 458 boundary = { (0, 0): 'Third', 459 (0, 2): 'First', 460 (2, 0): 'Second', 461 (2, 1): 'Second', 462 (3, 1): 'Second', 463 (3, 2): 'Third'} 464 441 # bac, bce, ecf, dbe 442 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 443 boundary = {(0, 0): 'Third', 444 (0, 2): 'First', 445 (2, 0): 'Second', 446 (2, 1): 'Second', 447 (3, 1): 'Second', 448 (3, 2): 'Third'} 465 449 466 450 domain = Domain(points, vertices, boundary) 467 451 domain.check_integrity() 468 452 469 470 domain.set_quantity('stage', [[1,2,3], [5,5,5], 471 [0,0,9], [-6, 3, 3]]) 472 473 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], 474 [3,3,3], [4, 4, 4]]) 453 domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]]) 454 455 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]]) 475 456 476 457 domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], 477 [30,30,30], [40, 40, 40]]) 478 479 480 D = Dirichlet_boundary([5,2,1]) 458 [30,30,30], [40,40,40]]) 459 460 D = Dirichlet_boundary([5, 2, 1]) 481 461 T = Transmissive_boundary(domain) 482 462 R = Reflective_boundary(domain) 483 domain.set_boundary( 463 domain.set_boundary({'First': D, 'Second': T, 'Third': R}) 484 464 485 465 domain.update_boundary() 486 466 487 # Stage467 # Stage 488 468 assert domain.quantities['stage'].boundary_values[0] == 2.5 489 assert domain.quantities['stage'].boundary_values[0] ==\ 490 domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5) 491 assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet 492 assert domain.quantities['stage'].boundary_values[2] ==\ 493 domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5) 494 assert domain.quantities['stage'].boundary_values[3] ==\ 495 domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5) 496 assert domain.quantities['stage'].boundary_values[4] ==\ 497 domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5) 498 assert domain.quantities['stage'].boundary_values[5] ==\ 499 domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5) 500 501 #Xmomentum 502 assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective 503 assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet 504 assert domain.quantities['xmomentum'].boundary_values[2] ==\ 505 domain.get_conserved_quantities(2, edge=0)[1] #Transmissive 506 assert domain.quantities['xmomentum'].boundary_values[3] ==\ 507 domain.get_conserved_quantities(2, edge=1)[1] #Transmissive 508 assert domain.quantities['xmomentum'].boundary_values[4] ==\ 509 domain.get_conserved_quantities(3, edge=1)[1] #Transmissive 510 assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective 511 512 #Ymomentum 513 assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective 514 assert domain.quantities['ymomentum'].boundary_values[1] == 1. #Dirichlet 515 assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive 516 assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive 517 assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive 518 assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective 519 469 # Reflective (2.5) 470 assert (domain.quantities['stage'].boundary_values[0] == 471 domain.get_conserved_quantities(0, edge=0)[0]) 472 # Dirichlet 473 assert domain.quantities['stage'].boundary_values[1] == 5. 474 # Transmissive (4.5) 475 assert (domain.quantities['stage'].boundary_values[2] == 476 domain.get_conserved_quantities(2, edge=0)[0]) 477 # Transmissive (4.5) 478 assert (domain.quantities['stage'].boundary_values[3] == 479 domain.get_conserved_quantities(2, edge=1)[0]) 480 # Transmissive (-1.5) 481 assert (domain.quantities['stage'].boundary_values[4] == 482 domain.get_conserved_quantities(3, edge=1)[0]) 483 # Reflective (-1.5) 484 assert (domain.quantities['stage'].boundary_values[5] == 485 domain.get_conserved_quantities(3, edge=2)[0]) 486 487 # Xmomentum 488 # Reflective 489 assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 490 # Dirichlet 491 assert domain.quantities['xmomentum'].boundary_values[1] == 2. 492 # Transmissive 493 assert (domain.quantities['xmomentum'].boundary_values[2] == 494 domain.get_conserved_quantities(2, edge=0)[1]) 495 # Transmissive 496 assert (domain.quantities['xmomentum'].boundary_values[3] == 497 domain.get_conserved_quantities(2, edge=1)[1]) 498 # Transmissive 499 assert (domain.quantities['xmomentum'].boundary_values[4] == 500 domain.get_conserved_quantities(3, edge=1)[1]) 501 # Reflective 502 assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 503 504 # Ymomentum 505 # Reflective 506 assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 507 # Dirichlet 508 assert domain.quantities['ymomentum'].boundary_values[1] == 1. 509 # Transmissive 510 assert domain.quantities['ymomentum'].boundary_values[2] == 30. 511 # Transmissive 512 assert domain.quantities['ymomentum'].boundary_values[3] == 30. 513 # Transmissive 514 assert domain.quantities['ymomentum'].boundary_values[4] == 40. 515 # Reflective 516 assert domain.quantities['ymomentum'].boundary_values[5] == 40. 520 517 521 518 def test_boundary_conditionsII(self): 522 523 519 a = [0.0, 0.0] 524 520 b = [0.0, 2.0] 525 c = [2.0, 0.0]521 c = [2.0, 0.0] 526 522 d = [0.0, 4.0] 527 523 e = [2.0, 2.0] 528 f = [4.0, 0.0]524 f = [4.0, 0.0] 529 525 530 526 points = [a, b, c, d, e, f] 531 #bac, bce, ecf, dbe 532 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 533 boundary = { (0, 0): 'Third', 534 (0, 2): 'First', 535 (2, 0): 'Second', 536 (2, 1): 'Second', 537 (3, 1): 'Second', 538 (3, 2): 'Third', 539 (0, 1): 'Internal'} 540 527 # bac, bce, ecf, dbe 528 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 529 boundary = {(0, 0): 'Third', 530 (0, 2): 'First', 531 (2, 0): 'Second', 532 (2, 1): 'Second', 533 (3, 1): 'Second', 534 (3, 2): 'Third', 535 (0, 1): 'Internal'} 541 536 542 537 domain = Domain(points, vertices, boundary) 543 538 domain.check_integrity() 544 539 545 546 domain.set_quantity('stage', [[1,2,3], [5,5,5], 547 [0,0,9], [-6, 3, 3]]) 548 549 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], 550 [3,3,3], [4, 4, 4]]) 540 domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]]) 541 542 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]]) 551 543 552 544 domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], 553 [30,30,30], [40, 40, 40]]) 554 555 556 D = Dirichlet_boundary([5,2,1]) 545 [30,30,30], [40,40,40]]) 546 547 D = Dirichlet_boundary([5, 2, 1]) 557 548 T = Transmissive_boundary(domain) 558 549 R = Reflective_boundary(domain) 559 domain.set_boundary( 560 550 domain.set_boundary({'First': D, 'Second': T, 551 'Third': R, 'Internal': None}) 561 552 562 553 domain.update_boundary() 563 554 domain.check_integrity() 564 555 565 566 567 568 556 def test_boundary_conditionsIII(self): 569 557 """test_boundary_conditionsIII 570 558 571 559 Test Transmissive_stage_zero_momentum_boundary 572 560 """ 573 561 574 562 a = [0.0, 0.0] 575 563 b = [0.0, 2.0] 576 c = [2.0, 0.0]564 c = [2.0, 0.0] 577 565 d = [0.0, 4.0] 578 566 e = [2.0, 2.0] 579 f = [4.0, 0.0]567 f = [4.0, 0.0] 580 568 581 569 points = [a, b, c, d, e, f] 582 #bac, bce, ecf, dbe 583 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 584 boundary = { (0, 0): 'Third', 585 (0, 2): 'First', 586 (2, 0): 'Second', 587 (2, 1): 'Second', 588 (3, 1): 'Second', 589 (3, 2): 'Third'} 590 570 # bac, bce, ecf, dbe 571 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 572 boundary = {(0, 0): 'Third', 573 (0, 2): 'First', 574 (2, 0): 'Second', 575 (2, 1): 'Second', 576 (3, 1): 'Second', 577 (3, 2): 'Third'} 591 578 592 579 domain = Domain(points, vertices, boundary) 593 580 domain.check_integrity() 594 581 595 596 domain.set_quantity('stage', [[1,2,3], [5,5,5], 597 [0,0,9], [-6, 3, 3]]) 598 599 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], 600 [3,3,3], [4, 4, 4]]) 582 domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]]) 583 584 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]]) 601 585 602 586 domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], 603 [30,30,30], [40, 40, 40]]) 604 605 606 D = Dirichlet_boundary([5,2,1]) 587 [30,30,30], [40,40,40]]) 588 589 D = Dirichlet_boundary([5, 2, 1]) 607 590 T = Transmissive_stage_zero_momentum_boundary(domain) 608 591 R = Reflective_boundary(domain) 609 domain.set_boundary( 592 domain.set_boundary({'First': D, 'Second': T, 'Third': R}) 610 593 611 594 domain.update_boundary() 612 595 613 596 # Stage 597 # Reflective (2.5) 614 598 assert domain.quantities['stage'].boundary_values[0] == 2.5 615 assert domain.quantities['stage'].boundary_values[0] ==\ 616 domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5) 617 assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet 618 assert domain.quantities['stage'].boundary_values[2] ==\ 619 domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5) 620 assert domain.quantities['stage'].boundary_values[3] ==\ 621 domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5) 622 assert domain.quantities['stage'].boundary_values[4] ==\ 623 domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5) 624 assert domain.quantities['stage'].boundary_values[5] ==\ 625 domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5) 599 assert (domain.quantities['stage'].boundary_values[0] == 600 domain.get_conserved_quantities(0, edge=0)[0]) 601 # Dirichlet 602 assert domain.quantities['stage'].boundary_values[1] == 5. 603 # Transmissive (4.5) 604 assert (domain.quantities['stage'].boundary_values[2] == 605 domain.get_conserved_quantities(2, edge=0)[0]) 606 # Transmissive (4.5) 607 assert (domain.quantities['stage'].boundary_values[3] == 608 domain.get_conserved_quantities(2, edge=1)[0]) 609 # Transmissive (-1.5) 610 assert (domain.quantities['stage'].boundary_values[4] == 611 domain.get_conserved_quantities(3, edge=1)[0]) 612 # Reflective (-1.5) 613 assert (domain.quantities['stage'].boundary_values[5] == 614 domain.get_conserved_quantities(3, edge=2)[0]) 626 615 627 616 # Xmomentum 628 assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective 629 assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet 630 assert domain.quantities['xmomentum'].boundary_values[2] == 0.0 631 assert domain.quantities['xmomentum'].boundary_values[3] == 0.0 632 assert domain.quantities['xmomentum'].boundary_values[4] == 0.0 633 assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective 617 # Reflective 618 assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 619 # Dirichlet 620 assert domain.quantities['xmomentum'].boundary_values[1] == 2. 621 assert domain.quantities['xmomentum'].boundary_values[2] == 0.0 622 assert domain.quantities['xmomentum'].boundary_values[3] == 0.0 623 assert domain.quantities['xmomentum'].boundary_values[4] == 0.0 624 # Reflective 625 assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 634 626 635 627 # Ymomentum 636 assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective637 assert domain.quantities['ymomentum'].boundary_values[ 1] == 1. #Dirichlet638 assert domain.quantities['ymomentum'].boundary_values[2] == 0.0639 assert domain.quantities['ymomentum'].boundary_values[ 3] == 0.0640 assert domain.quantities['ymomentum'].boundary_values[ 4] == 0.0641 assert domain.quantities['ymomentum'].boundary_values[ 5] == 40. #Reflective642 643 644 645 628 # Reflective 629 assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 630 # Dirichlet 631 assert domain.quantities['ymomentum'].boundary_values[1] == 1. 632 assert domain.quantities['ymomentum'].boundary_values[2] == 0.0 633 assert domain.quantities['ymomentum'].boundary_values[3] == 0.0 634 assert domain.quantities['ymomentum'].boundary_values[4] == 0.0 635 # Reflective 636 assert domain.quantities['ymomentum'].boundary_values[5] == 40. 637 646 638 def test_boundary_condition_time(self): 647 639 """test_boundary_condition_time 648 640 649 641 This tests that boundary conditions are evaluated 650 642 using the right time from domain. 651 652 643 """ 653 644 654 645 # Setup computational domain 655 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross656 646 from anuga.abstract_2d_finite_volumes.mesh_factory \ 647 import rectangular_cross 657 648 658 649 #-------------------------------------------------------------- … … 660 651 #-------------------------------------------------------------- 661 652 N = 5 662 points, vertices, boundary = rectangular_cross(N, N) 653 points, vertices, boundary = rectangular_cross(N, N) 663 654 domain = Domain(points, vertices, boundary) 664 655 … … 670 661 domain.set_quantity('stage', 0.0) 671 662 672 673 663 #-------------------------------------------------------------- 674 664 # Setup boundary conditions 675 665 #-------------------------------------------------------------- 676 Bt = Time_boundary(domain=domain, # Time dependent boundary677 678 666 # Time dependent boundary 667 Bt = Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0]) 668 679 669 Br = Reflective_boundary(domain) # Reflective wall 680 670 681 671 domain.set_boundary({'left': Bt, 'right': Br, 'top': Br, 'bottom': Br}) 682 672 683 673 for t in domain.evolve(yieldstep = 10, finaltime = 20.0): 684 674 q = Bt.evaluate() 685 686 # FIXME (Ole): This test would not have passed in 675 676 # FIXME (Ole): This test would not have passed in 687 677 # changeset:5846. 688 678 msg = 'Time boundary not evaluated correctly' 689 679 assert num.allclose(t, q[0]), msg 690 691 692 680 693 681 def test_compute_fluxes0(self): … … 697 685 a = [0.0, 0.0] 698 686 b = [0.0, 2.0] 699 c = [2.0, 0.0]687 c = [2.0, 0.0] 700 688 d = [0.0, 4.0] 701 689 e = [2.0, 2.0] 702 f = [4.0, 0.0]690 f = [4.0, 0.0] 703 691 704 692 points = [a, b, c, d, e, f] 705 # bac, bce, ecf,dbe706 vertices = [ 693 # bac, bce, ecf, dbe 694 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 707 695 708 696 domain = Domain(points, vertices) 709 domain.set_quantity('stage', [[2,2,2], [2,2,2], 710 [2,2,2], [2,2,2]]) 697 domain.set_quantity('stage', [[2,2,2], [2,2,2], [2,2,2], [2,2,2]]) 711 698 domain.check_integrity() 712 699 713 assert num.allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]]) 714 assert num.allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]]) 715 716 zl=zr=0. # Assume flat bed 717 718 edgeflux = num.zeros(3, num.float) 700 assert num.allclose(domain.neighbours, 701 [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]]) 702 assert num.allclose(domain.neighbour_edges, 703 [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]]) 704 705 zl = zr = 0. # Assume flat bed 706 707 edgeflux = num.zeros(3, num.float) 719 708 edgeflux0 = num.zeros(3, num.float) 720 709 edgeflux1 = num.zeros(3, num.float) 721 edgeflux2 = num.zeros(3, num.float) 722 H0 = 0.0 710 edgeflux2 = num.zeros(3, num.float) 711 H0 = 0.0 723 712 724 713 # Flux across right edge of volume 1 725 normal = domain.get_normal(1, 0)714 normal = domain.get_normal(1, 0) 726 715 ql = domain.get_conserved_quantities(vol_id=1, edge=0) 727 716 qr = domain.get_conserved_quantities(vol_id=2, edge=2) 728 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0) 717 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, 718 epsilon, g, H0) 729 719 730 720 # Check that flux seen from other triangles is inverse 731 tmp = qr; qr=ql; ql=tmp 732 normal = domain.get_normal(2,2) 733 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) 721 (ql, qr) = (qr, ql) 722 #tmp = qr 723 #qr = ql 724 #ql = tmp 725 normal = domain.get_normal(2, 2) 726 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, 727 epsilon, g, H0) 734 728 735 729 assert num.allclose(edgeflux0 + edgeflux, 0.) 736 730 737 731 # Flux across upper edge of volume 1 738 normal = domain.get_normal(1, 1)732 normal = domain.get_normal(1, 1) 739 733 ql = domain.get_conserved_quantities(vol_id=1, edge=1) 740 734 qr = domain.get_conserved_quantities(vol_id=3, edge=0) 741 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0) 735 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, 736 epsilon, g, H0) 742 737 743 738 # Check that flux seen from other triangles is inverse 744 tmp = qr; qr=ql; ql=tmp 745 normal = domain.get_normal(3,0) 746 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) 747 748 assert num.allclose(edgeflux1 + edgeflux, 0.) 749 739 (ql, qr) = (qr, ql) 740 #tmp = qr 741 #qr = ql 742 #ql = tmp 743 normal = domain.get_normal(3, 0) 744 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, 745 epsilon, g, H0) 746 747 assert num.allclose(edgeflux1 + edgeflux, 0.) 750 748 751 749 # Flux across lower left hypotenuse of volume 1 752 normal = domain.get_normal(1, 2)750 normal = domain.get_normal(1, 2) 753 751 ql = domain.get_conserved_quantities(vol_id=1, edge=2) 754 752 qr = domain.get_conserved_quantities(vol_id=0, edge=1) 755 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0) 753 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, 754 epsilon, g, H0) 756 755 757 756 # Check that flux seen from other triangles is inverse 758 tmp = qr; qr=ql; ql=tmp 759 normal = domain.get_normal(0,1) 760 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) 757 (ql, qr) = (qr, ql) 758 #tmp = qr 759 #qr=ql 760 #ql=tmp 761 normal = domain.get_normal(0, 1) 762 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, 763 epsilon, g, H0) 761 764 assert num.allclose(edgeflux2 + edgeflux, 0.) 762 763 765 764 766 # Scale by edgelengths, add up anc check that total flux is zero … … 767 769 e2 = domain.edgelengths[1, 2] 768 770 769 assert num.allclose(e0*edgeflux0 +e1*edgeflux1+e2*edgeflux2, 0.)771 assert num.allclose(e0*edgeflux0 + e1*edgeflux1 + e2*edgeflux2, 0.) 770 772 771 773 # Now check that compute_flux yields zeros as well … … 773 775 774 776 for name in ['stage', 'xmomentum', 'ymomentum']: 775 #print name, domain.quantities[name].explicit_update776 777 assert num.allclose(domain.quantities[name].explicit_update[1], 0) 777 778 779 778 780 779 def test_compute_fluxes1(self): 781 780 #Use values from previous version 782 783 781 a = [0.0, 0.0] 784 782 b = [0.0, 2.0] 785 c = [2.0, 0.0]783 c = [2.0, 0.0] 786 784 d = [0.0, 4.0] 787 785 e = [2.0, 2.0] 788 f = [4.0, 0.0]786 f = [4.0, 0.0] 789 787 790 788 points = [a, b, c, d, e, f] 791 # bac, bce, ecf,dbe792 vertices = [ 789 # bac, bce, ecf, dbe 790 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 793 791 794 792 domain = Domain(points, vertices) 795 val0 = 2. +2.0/3796 val1 = 4. +4.0/3797 val2 = 8. +2.0/3798 val3 = 2. +8.0/3793 val0 = 2. + 2.0/3 794 val1 = 4. + 4.0/3 795 val2 = 8. + 2.0/3 796 val3 = 2. + 8.0/3 799 797 800 798 domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1], … … 802 800 domain.check_integrity() 803 801 804 zl =zr=0. #Assume flat bed805 806 edgeflux = num.zeros(3, num.float) 802 zl = zr = 0. # Assume flat bed 803 804 edgeflux = num.zeros(3, num.float) 807 805 edgeflux0 = num.zeros(3, num.float) 808 806 edgeflux1 = num.zeros(3, num.float) 809 edgeflux2 = num.zeros(3, num.float) 810 H0 = 0.0 811 807 edgeflux2 = num.zeros(3, num.float) 808 H0 = 0.0 812 809 813 810 # Flux across right edge of volume 1 814 normal = domain.get_normal(1, 0) #Get normal 0 of triangle 1811 normal = domain.get_normal(1, 0) # Get normal 0 of triangle 1 815 812 assert num.allclose(normal, [1, 0]) 816 813 817 814 ql = domain.get_conserved_quantities(vol_id=1, edge=0) 818 815 assert num.allclose(ql, [val1, 0, 0]) 819 816 820 817 qr = domain.get_conserved_quantities(vol_id=2, edge=2) 821 818 assert num.allclose(qr, [val2, 0, 0]) 822 819 823 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0) 820 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, 821 epsilon, g, H0) 824 822 825 823 # Flux across edge in the east direction (as per normal vector) … … 827 825 assert num.allclose(max_speed, 9.21592824046) 828 826 829 830 827 #Flux across edge in the west direction (opposite sign for xmomentum) 831 normal_opposite = domain.get_normal(2, 2) #Get normal 2 of triangle 2828 normal_opposite = domain.get_normal(2, 2) # Get normal 2 of triangle 2 832 829 assert num.allclose(normal_opposite, [-1, 0]) 833 830 834 max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux, epsilon, g, H0)835 #flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr)831 max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux, 832 epsilon, g, H0) 836 833 assert num.allclose(edgeflux, [-15.3598804, -253.71111111, 0.]) 837 838 834 839 835 #Flux across upper edge of volume 1 840 normal = domain.get_normal(1, 1)836 normal = domain.get_normal(1, 1) 841 837 ql = domain.get_conserved_quantities(vol_id=1, edge=1) 842 838 qr = domain.get_conserved_quantities(vol_id=3, edge=0) 843 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0) 839 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, 840 epsilon, g, H0) 844 841 845 842 assert num.allclose(edgeflux1, [2.4098563, 0., 123.04444444]) … … 847 844 848 845 #Flux across lower left hypotenuse of volume 1 849 normal = domain.get_normal(1, 2)846 normal = domain.get_normal(1, 2) 850 847 ql = domain.get_conserved_quantities(vol_id=1, edge=2) 851 848 qr = domain.get_conserved_quantities(vol_id=0, edge=1) 852 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0) 849 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, 850 epsilon, g, H0) 853 851 854 852 assert num.allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738]) … … 860 858 e2 = domain.edgelengths[1, 2] 861 859 862 total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1] 860 total_flux = -(e0*edgeflux0 + 861 e1*edgeflux1 + 862 e2*edgeflux2) / domain.areas[1] 863 863 864 assert num.allclose(total_flux, [-0.68218178, -166.6, -35.93333333]) 864 865 865 866 866 domain.compute_fluxes() 867 867 868 #assert num.allclose(total_flux, domain.explicit_update[1,:])869 868 for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']): 870 869 assert num.allclose(total_flux[i], 871 870 domain.quantities[name].explicit_update[1]) 872 873 #assert allclose(domain.explicit_update, [874 # [0., -69.68888889, -69.68888889],875 # [-0.68218178, -166.6, -35.93333333],876 # [-111.77316251, 69.68888889, 0.],877 # [-35.68522449, 0., 69.68888889]])878 871 879 872 assert num.allclose(domain.quantities['stage'].explicit_update, … … 884 877 [-69.68888889, -35.93333333, 0., 69.68888889]) 885 878 886 887 #assert allclose(domain.quantities[name].explicit_update888 889 890 891 892 893 879 def test_compute_fluxes2(self): 894 880 #Random values, incl momentum 895 896 881 a = [0.0, 0.0] 897 882 b = [0.0, 2.0] 898 c = [2.0, 0.0]883 c = [2.0, 0.0] 899 884 d = [0.0, 4.0] 900 885 e = [2.0, 2.0] 901 f = [4.0, 0.0]886 f = [4.0, 0.0] 902 887 903 888 points = [a, b, c, d, e, f] 904 # bac, bce, ecf,dbe905 vertices = [ 889 # bac, bce, ecf, dbe 890 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 906 891 907 892 domain = Domain(points, vertices) 893 val0 = 2. + 2.0/3 894 val1 = 4. + 4.0/3 895 val2 = 8. + 2.0/3 896 val3 = 2. + 8.0/3 897 898 zl = zr = 0 # Assume flat zero bed 899 edgeflux = num.zeros(3, num.float) 900 edgeflux0 = num.zeros(3, num.float) 901 edgeflux1 = num.zeros(3, num.float) 902 edgeflux2 = num.zeros(3, num.float) 903 H0 = 0.0 904 905 domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default# 906 907 domain.set_quantity('stage', [[val0, val0-1, val0-2], 908 [val1, val1+1, val1], 909 [val2, val2-2, val2], 910 [val3-0.5, val3, val3]]) 911 912 domain.set_quantity('xmomentum', 913 [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]]) 914 915 domain.set_quantity('ymomentum', 916 [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]]) 917 918 domain.check_integrity() 919 920 # Flux across right edge of volume 1 921 normal = domain.get_normal(1, 0) 922 ql = domain.get_conserved_quantities(vol_id=1, edge=0) 923 qr = domain.get_conserved_quantities(vol_id=2, edge=2) 924 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, 925 epsilon, g, H0) 926 927 # Flux across upper edge of volume 1 928 normal = domain.get_normal(1, 1) 929 ql = domain.get_conserved_quantities(vol_id=1, edge=1) 930 qr = domain.get_conserved_quantities(vol_id=3, edge=0) 931 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, 932 epsilon, g, H0) 933 934 # Flux across lower left hypotenuse of volume 1 935 normal = domain.get_normal(1, 2) 936 ql = domain.get_conserved_quantities(vol_id=1, edge=2) 937 qr = domain.get_conserved_quantities(vol_id=0, edge=1) 938 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, 939 epsilon, g, H0) 940 941 # Scale, add up and check that compute_fluxes is correct for vol 1 942 e0 = domain.edgelengths[1, 0] 943 e1 = domain.edgelengths[1, 1] 944 e2 = domain.edgelengths[1, 2] 945 946 total_flux = -(e0*edgeflux0 + 947 e1*edgeflux1 + 948 e2*edgeflux2) / domain.areas[1] 949 950 domain.compute_fluxes() 951 952 for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']): 953 assert num.allclose(total_flux[i], 954 domain.quantities[name].explicit_update[1]) 955 956 # FIXME (Ole): Need test like this for fluxes in very shallow water. 957 def test_compute_fluxes3(self): 958 #Random values, incl momentum 959 a = [0.0, 0.0] 960 b = [0.0, 2.0] 961 c = [2.0, 0.0] 962 d = [0.0, 4.0] 963 e = [2.0, 2.0] 964 f = [4.0, 0.0] 965 966 points = [a, b, c, d, e, f] 967 # bac, bce, ecf, dbe 968 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 969 970 domain = Domain(points, vertices) 971 908 972 val0 = 2.+2.0/3 909 973 val1 = 4.+4.0/3 … … 911 975 val3 = 2.+8.0/3 912 976 913 zl=zr=0 #Assume flat zero bed 914 edgeflux = num.zeros(3, num.float) 977 zl = zr = -3.75 # Assume constant bed (must be less than stage) 978 domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default# 979 980 edgeflux = num.zeros(3, num.float) 915 981 edgeflux0 = num.zeros(3, num.float) 916 982 edgeflux1 = num.zeros(3, num.float) 917 edgeflux2 = num.zeros(3, num.float) 918 H0 = 0.0 919 920 921 domain.set_quantity('elevation', zl*num.ones( (4,3), num.int )) #array default# 922 983 edgeflux2 = num.zeros(3, num.float) 984 H0 = 0.0 923 985 924 986 domain.set_quantity('stage', [[val0, val0-1, val0-2], … … 928 990 929 991 domain.set_quantity('xmomentum', 930 [[1, 2, 3], [3, 4, 5], 931 [1, -1, 0], [0, -2, 2]]) 992 [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]]) 932 993 933 994 domain.set_quantity('ymomentum', 934 [[1, -1, 0], [0, -3, 2], 935 [0, 1, 0], [-1, 2, 2]]) 936 995 [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]]) 937 996 938 997 domain.check_integrity() 939 998 940 941 942 #Flux across right edge of volume 1 943 normal = domain.get_normal(1,0) 999 # Flux across right edge of volume 1 1000 normal = domain.get_normal(1, 0) 944 1001 ql = domain.get_conserved_quantities(vol_id=1, edge=0) 945 1002 qr = domain.get_conserved_quantities(vol_id=2, edge=2) 946 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0) 947 948 #Flux across upper edge of volume 1 949 normal = domain.get_normal(1,1) 1003 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, 1004 epsilon, g, H0) 1005 1006 # Flux across upper edge of volume 1 1007 normal = domain.get_normal(1, 1) 950 1008 ql = domain.get_conserved_quantities(vol_id=1, edge=1) 951 1009 qr = domain.get_conserved_quantities(vol_id=3, edge=0) 952 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0) 953 954 #Flux across lower left hypotenuse of volume 1 955 normal = domain.get_normal(1,2) 1010 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, 1011 epsilon, g, H0) 1012 1013 # Flux across lower left hypotenuse of volume 1 1014 normal = domain.get_normal(1, 2) 956 1015 ql = domain.get_conserved_quantities(vol_id=1, edge=2) 957 1016 qr = domain.get_conserved_quantities(vol_id=0, edge=1) 958 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0) 959 960 #Scale, add up and check that compute_fluxes is correct for vol 1 1017 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, 1018 epsilon, g, H0) 1019 1020 # Scale, add up and check that compute_fluxes is correct for vol 1 961 1021 e0 = domain.edgelengths[1, 0] 962 1022 e1 = domain.edgelengths[1, 1] 963 1023 e2 = domain.edgelengths[1, 2] 964 1024 965 total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1] 966 1025 total_flux = -(e0*edgeflux0 + 1026 e1*edgeflux1 + 1027 e2*edgeflux2) / domain.areas[1] 967 1028 968 1029 domain.compute_fluxes() 1030 969 1031 for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']): 970 1032 assert num.allclose(total_flux[i], 971 1033 domain.quantities[name].explicit_update[1]) 972 #assert allclose(total_flux, domain.explicit_update[1,:]) 973 974 975 # FIXME (Ole): Need test like this for fluxes in very shallow water. 976 def test_compute_fluxes3(self): 977 #Random values, incl momentum 1034 1035 def xtest_catching_negative_heights(self): 1036 #OBSOLETE 978 1037 979 1038 a = [0.0, 0.0] 980 1039 b = [0.0, 2.0] 981 c = [2.0, 0.0]1040 c = [2.0, 0.0] 982 1041 d = [0.0, 4.0] 983 1042 e = [2.0, 2.0] 984 f = [4.0, 0.0]1043 f = [4.0, 0.0] 985 1044 986 1045 points = [a, b, c, d, e, f] 987 # bac, bce, ecf,dbe988 vertices = [ 1046 # bac, bce, ecf, dbe 1047 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 989 1048 990 1049 domain = Domain(points, vertices) 991 val0 = 2.+2.0/3 992 val1 = 4.+4.0/3 993 val2 = 8.+2.0/3 994 val3 = 2.+8.0/3 995 996 zl=zr=-3.75 #Assume constant bed (must be less than stage) 997 domain.set_quantity('elevation', zl*num.ones( (4,3), num.int )) #array default# 998 999 1000 edgeflux = num.zeros(3, num.float) 1001 edgeflux0 = num.zeros(3, num.float) 1002 edgeflux1 = num.zeros(3, num.float) 1003 edgeflux2 = num.zeros(3, num.float) 1004 H0 = 0.0 1005 1006 1007 1008 domain.set_quantity('stage', [[val0, val0-1, val0-2], 1009 [val1, val1+1, val1], 1010 [val2, val2-2, val2], 1011 [val3-0.5, val3, val3]]) 1012 1013 domain.set_quantity('xmomentum', 1014 [[1, 2, 3], [3, 4, 5], 1015 [1, -1, 0], [0, -2, 2]]) 1016 1017 domain.set_quantity('ymomentum', 1018 [[1, -1, 0], [0, -3, 2], 1019 [0, 1, 0], [-1, 2, 2]]) 1020 1021 1022 domain.check_integrity() 1023 1024 1025 1026 #Flux across right edge of volume 1 1027 normal = domain.get_normal(1,0) 1028 ql = domain.get_conserved_quantities(vol_id=1, edge=0) 1029 qr = domain.get_conserved_quantities(vol_id=2, edge=2) 1030 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0) 1031 1032 #Flux across upper edge of volume 1 1033 normal = domain.get_normal(1,1) 1034 ql = domain.get_conserved_quantities(vol_id=1, edge=1) 1035 qr = domain.get_conserved_quantities(vol_id=3, edge=0) 1036 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0) 1037 1038 #Flux across lower left hypotenuse of volume 1 1039 normal = domain.get_normal(1,2) 1040 ql = domain.get_conserved_quantities(vol_id=1, edge=2) 1041 qr = domain.get_conserved_quantities(vol_id=0, edge=1) 1042 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0) 1043 1044 #Scale, add up and check that compute_fluxes is correct for vol 1 1045 e0 = domain.edgelengths[1, 0] 1046 e1 = domain.edgelengths[1, 1] 1047 e2 = domain.edgelengths[1, 2] 1048 1049 total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1] 1050 1051 domain.compute_fluxes() 1052 for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']): 1053 assert num.allclose(total_flux[i], 1054 domain.quantities[name].explicit_update[1]) 1055 1056 1057 1058 def xtest_catching_negative_heights(self): 1059 1060 #OBSOLETE 1061 1062 a = [0.0, 0.0] 1063 b = [0.0, 2.0] 1064 c = [2.0,0.0] 1065 d = [0.0, 4.0] 1066 e = [2.0, 2.0] 1067 f = [4.0,0.0] 1068 1069 points = [a, b, c, d, e, f] 1070 #bac, bce, ecf, dbe 1071 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1072 1073 domain = Domain(points, vertices) 1074 val0 = 2.+2.0/3 1075 val1 = 4.+4.0/3 1076 val2 = 8.+2.0/3 1077 val3 = 2.+8.0/3 1078 1079 zl=zr=4 #Too large 1080 domain.set_quantity('elevation', zl*num.ones( (4,3), num.int )) #array default# 1050 val0 = 2. + 2.0/3 1051 val1 = 4. + 4.0/3 1052 val2 = 8. + 2.0/3 1053 val3 = 2. + 8.0/3 1054 1055 zl = zr = 4 # Too large 1056 domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default# 1081 1057 domain.set_quantity('stage', [[val0, val0-1, val0-2], 1082 1058 [val1, val1+1, val1], … … 1090 1066 pass 1091 1067 1092 1093 1094 1068 def test_get_wet_elements(self): 1095 1096 1069 a = [0.0, 0.0] 1097 1070 b = [0.0, 2.0] 1098 c = [2.0, 0.0]1071 c = [2.0, 0.0] 1099 1072 d = [0.0, 4.0] 1100 1073 e = [2.0, 2.0] 1101 f = [4.0, 0.0]1074 f = [4.0, 0.0] 1102 1075 1103 1076 points = [a, b, c, d, e, f] 1104 # bac, bce, ecf,dbe1105 vertices = [ 1077 # bac, bce, ecf, dbe 1078 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1106 1079 1107 1080 domain = Domain(points, vertices) 1108 val0 = 2.+2.0/3 1109 val1 = 4.+4.0/3 1110 val2 = 8.+2.0/3 1111 val3 = 2.+8.0/3 1112 1113 zl=zr=5 1114 domain.set_quantity('elevation', zl*num.ones( (4,3), num.int )) #array default# 1081 1082 val0 = 2. + 2.0/3 1083 val1 = 4. + 4.0/3 1084 val2 = 8. + 2.0/3 1085 val3 = 2. + 8.0/3 1086 1087 zl = zr = 5 1088 domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default# 1115 1089 domain.set_quantity('stage', [[val0, val0-1, val0-2], 1116 1090 [val1, val1+1, val1], … … 1118 1092 [val3-0.5, val3, val3]]) 1119 1093 1120 1121 1122 #print domain.get_quantity('elevation').get_values(location='centroids')1123 #print domain.get_quantity('stage').get_values(location='centroids')1124 1094 domain.check_integrity() 1125 1095 1126 1096 indices = domain.get_wet_elements() 1127 assert num.allclose(indices, [1, 2])1128 1129 indices = domain.get_wet_elements(indices=[0, 1,3])1097 assert num.allclose(indices, [1, 2]) 1098 1099 indices = domain.get_wet_elements(indices=[0, 1, 3]) 1130 1100 assert num.allclose(indices, [1]) 1131 1132 1133 1101 1134 1102 def test_get_maximum_inundation_1(self): 1135 1136 1103 a = [0.0, 0.0] 1137 1104 b = [0.0, 2.0] 1138 c = [2.0, 0.0]1105 c = [2.0, 0.0] 1139 1106 d = [0.0, 4.0] 1140 1107 e = [2.0, 2.0] 1141 f = [4.0, 0.0]1108 f = [4.0, 0.0] 1142 1109 1143 1110 points = [a, b, c, d, e, f] 1144 # bac, bce, ecf,dbe1145 vertices = [ 1111 # bac, bce, ecf, dbe 1112 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1146 1113 1147 1114 domain = Domain(points, vertices) 1148 1115 1149 domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 61116 domain.set_quantity('elevation', lambda x, y: x+2*y) # 2 4 4 6 1150 1117 domain.set_quantity('stage', 3) 1151 1118 … … 1156 1123 1157 1124 q = domain.get_maximum_inundation_elevation() 1158 assert num.allclose(q, domain.get_quantity('elevation').get_values(location='centroids')[0]) 1125 assert num.allclose(q, domain.get_quantity('elevation').\ 1126 get_values(location='centroids')[0]) 1159 1127 1160 1128 x, y = domain.get_maximum_inundation_location() 1161 1129 assert num.allclose([x, y], domain.get_centroid_coordinates()[0]) 1162 1130 1163 1164 1131 def test_get_maximum_inundation_2(self): 1165 1132 """test_get_maximum_inundation_2(self) … … 1167 1134 Test multiple wet cells with same elevation 1168 1135 """ 1169 1136 1170 1137 a = [0.0, 0.0] 1171 1138 b = [0.0, 2.0] 1172 c = [2.0, 0.0]1139 c = [2.0, 0.0] 1173 1140 d = [0.0, 4.0] 1174 1141 e = [2.0, 2.0] 1175 f = [4.0, 0.0]1142 f = [4.0, 0.0] 1176 1143 1177 1144 points = [a, b, c, d, e, f] 1178 # bac, bce, ecf,dbe1179 vertices = [ 1145 # bac, bce, ecf, dbe 1146 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1180 1147 1181 1148 domain = Domain(points, vertices) 1182 1149 1183 domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 61150 domain.set_quantity('elevation', lambda x, y: x+2*y) # 2 4 4 6 1184 1151 domain.set_quantity('stage', 4.1) 1185 1152 … … 1187 1154 1188 1155 indices = domain.get_wet_elements() 1189 assert num.allclose(indices, [0, 1,2])1156 assert num.allclose(indices, [0, 1, 2]) 1190 1157 1191 1158 q = domain.get_maximum_inundation_elevation() 1192 assert num.allclose(q, 4) 1159 assert num.allclose(q, 4) 1193 1160 1194 1161 x, y = domain.get_maximum_inundation_location() 1195 assert num.allclose([x, y], domain.get_centroid_coordinates()[1]) 1196 1162 assert num.allclose([x, y], domain.get_centroid_coordinates()[1]) 1197 1163 1198 1164 def test_get_maximum_inundation_3(self): … … 1202 1168 """ 1203 1169 1204 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 1170 from anuga.abstract_2d_finite_volumes.mesh_factory \ 1171 import rectangular_cross 1205 1172 1206 1173 initial_runup_height = -0.4 1207 1174 final_runup_height = -0.3 1208 1209 1175 1210 1176 #-------------------------------------------------------------- … … 1212 1178 #-------------------------------------------------------------- 1213 1179 N = 5 1214 points, vertices, boundary = rectangular_cross(N, N) 1180 points, vertices, boundary = rectangular_cross(N, N) 1215 1181 domain = Domain(points, vertices, boundary) 1216 domain.set_maximum_allowed_speed(1.0) 1182 domain.set_maximum_allowed_speed(1.0) 1217 1183 1218 1184 #-------------------------------------------------------------- 1219 1185 # Setup initial conditions 1220 1186 #-------------------------------------------------------------- 1221 def topography(x, y):1187 def topography(x, y): 1222 1188 return -x/2 # linear bed slope 1223 1224 1225 domain.set_quantity('elevation', topography) # Use function for elevation1226 domain.set_quantity('friction', 0.) # Zero friction 1227 domain.set_quantity('stage', initial_runup_height)# Constant negative initial stage1228 1189 1190 # Use function for elevation 1191 domain.set_quantity('elevation', topography) 1192 domain.set_quantity('friction', 0.) # Zero friction 1193 # Constant negative initial stage 1194 domain.set_quantity('stage', initial_runup_height) 1229 1195 1230 1196 #-------------------------------------------------------------- 1231 1197 # Setup boundary conditions 1232 1198 #-------------------------------------------------------------- 1233 Br = Reflective_boundary(domain) # Reflective wall 1234 Bd = Dirichlet_boundary([final_runup_height, # Constant inflow 1235 0, 1236 0]) 1237 1238 # All reflective to begin with (still water) 1199 Br = Reflective_boundary(domain) # Reflective wall 1200 Bd = Dirichlet_boundary([final_runup_height, 0, 0]) # Constant inflow 1201 1202 # All reflective to begin with (still water) 1239 1203 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 1240 1241 1204 1242 1205 #-------------------------------------------------------------- … … 1245 1208 1246 1209 indices = domain.get_wet_elements() 1247 z = domain.get_quantity('elevation'). \1248 get_values(location='centroids',indices=indices)1249 assert num.alltrue(z <initial_runup_height)1210 z = domain.get_quantity('elevation').get_values(location='centroids', 1211 indices=indices) 1212 assert num.alltrue(z < initial_runup_height) 1250 1213 1251 1214 q = domain.get_maximum_inundation_elevation() 1252 assert num.allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy 1215 # First order accuracy 1216 assert num.allclose(q, initial_runup_height, rtol=1.0/N) 1253 1217 1254 1218 x, y = domain.get_maximum_inundation_location() 1255 1219 1256 qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]]) 1220 qref = domain.get_quantity('elevation').\ 1221 get_values(interpolation_points=[[x, y]]) 1257 1222 assert num.allclose(q, qref) 1258 1223 1259 1260 1224 wet_elements = domain.get_wet_elements() 1261 wet_elevations = domain.get_quantity('elevation').get_values(location='centroids', 1262 indices=wet_elements) 1263 assert num.alltrue(wet_elevations<initial_runup_height) 1264 assert num.allclose(wet_elevations, z) 1265 1266 1267 #print domain.get_quantity('elevation').get_maximum_value(indices=wet_elements) 1268 #print domain.get_quantity('elevation').get_maximum_location(indices=wet_elements) 1269 #print domain.get_quantity('elevation').get_maximum_index(indices=wet_elements) 1270 1271 1225 wet_elevations = domain.get_quantity('elevation').\ 1226 get_values(location='centroids', 1227 indices=wet_elements) 1228 assert num.alltrue(wet_elevations < initial_runup_height) 1229 assert num.allclose(wet_elevations, z) 1230 1272 1231 #-------------------------------------------------------------- 1273 1232 # Let triangles adjust … … 1276 1235 pass 1277 1236 1278 1279 1237 #-------------------------------------------------------------- 1280 1238 # Test inundation height again 1281 1239 #-------------------------------------------------------------- 1282 1283 1240 indices = domain.get_wet_elements() 1284 z = domain.get_quantity('elevation'). \1285 get_values(location='centroids',indices=indices)1286 1287 assert num.alltrue(z <initial_runup_height)1241 z = domain.get_quantity('elevation').get_values(location='centroids', 1242 indices=indices) 1243 1244 assert num.alltrue(z < initial_runup_height) 1288 1245 1289 1246 q = domain.get_maximum_inundation_elevation() 1290 assert num.allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy 1291 1247 # First order accuracy 1248 assert num.allclose(q, initial_runup_height, rtol=1.0/N) 1249 1292 1250 x, y = domain.get_maximum_inundation_location() 1293 qref = domain.get_quantity('elevation'). get_values(interpolation_points = [[x, y]])1294 assert num.allclose(q, qref)1295 1251 qref = domain.get_quantity('elevation').\ 1252 get_values(interpolation_points=[[x, y]]) 1253 assert num.allclose(q, qref) 1296 1254 1297 1255 #-------------------------------------------------------------- … … 1300 1258 domain.set_boundary({'right': Bd}) 1301 1259 1302 1303 1260 #-------------------------------------------------------------- 1304 1261 # Evolve system through time 1305 1262 #-------------------------------------------------------------- 1306 1263 for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0): 1307 #print domain.timestepping_statistics(track_speeds=True)1308 #domain.write_time()1309 1264 pass 1310 1265 1311 1266 #-------------------------------------------------------------- 1312 1267 # Test inundation height again 1313 1268 #-------------------------------------------------------------- 1314 1315 1269 indices = domain.get_wet_elements() 1316 1270 z = domain.get_quantity('elevation').\ 1317 get_values(location='centroids', indices=indices)1318 1319 assert num.alltrue(z <final_runup_height)1271 get_values(location='centroids', indices=indices) 1272 1273 assert num.alltrue(z < final_runup_height) 1320 1274 1321 1275 q = domain.get_maximum_inundation_elevation() 1322 assert num.allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy 1276 # First order accuracy 1277 assert num.allclose(q, final_runup_height, rtol=1.0/N) 1323 1278 1324 1279 x, y = domain.get_maximum_inundation_location() 1325 qref = domain.get_quantity('elevation'). get_values(interpolation_points = [[x, y]])1326 assert num.allclose(q, qref)1327 1280 qref = domain.get_quantity('elevation').\ 1281 get_values(interpolation_points=[[x, y]]) 1282 assert num.allclose(q, qref) 1328 1283 1329 1284 wet_elements = domain.get_wet_elements() 1330 wet_elevations = domain.get_quantity('elevation').get_values(location='centroids', 1331 indices=wet_elements) 1332 assert num.alltrue(wet_elevations<final_runup_height) 1333 assert num.allclose(wet_elevations, z) 1334 1335 1285 wet_elevations = domain.get_quantity('elevation').\ 1286 get_values(location='centroids', 1287 indices=wet_elements) 1288 assert num.alltrue(wet_elevations < final_runup_height) 1289 assert num.allclose(wet_elevations, z) 1336 1290 1337 1291 def test_get_maximum_inundation_from_sww(self): … … 1340 1294 Test of get_maximum_inundation_elevation() 1341 1295 and get_maximum_inundation_location() from data_manager.py 1342 1296 1343 1297 This is based on test_get_maximum_inundation_3(self) but works with the 1344 1298 stored results instead of with the internal data structure. … … 1347 1301 """ 1348 1302 1349 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 1303 from anuga.abstract_2d_finite_volumes.mesh_factory \ 1304 import rectangular_cross 1350 1305 from data_manager import get_maximum_inundation_elevation 1351 1306 from data_manager import get_maximum_inundation_location 1352 1307 from data_manager import get_maximum_inundation_data 1353 1354 1308 1355 1309 initial_runup_height = -0.4 1356 1310 final_runup_height = -0.3 1357 1358 1311 1359 1312 #-------------------------------------------------------------- … … 1361 1314 #-------------------------------------------------------------- 1362 1315 N = 10 1363 points, vertices, boundary = rectangular_cross(N, N) 1316 points, vertices, boundary = rectangular_cross(N, N) 1364 1317 domain = Domain(points, vertices, boundary) 1365 1318 domain.set_name('runup_test') 1366 1319 domain.set_maximum_allowed_speed(1.0) 1367 1320 1368 domain.tight_slope_limiters = 0 # FIXME: This works better with old limiters so far 1321 # FIXME: This works better with old limiters so far 1322 domain.tight_slope_limiters = 0 1369 1323 1370 1324 #-------------------------------------------------------------- 1371 1325 # Setup initial conditions 1372 1326 #-------------------------------------------------------------- 1373 def topography(x, y):1327 def topography(x, y): 1374 1328 return -x/2 # linear bed slope 1375 1376 1377 domain.set_quantity('elevation', topography) # Use function for elevation1378 domain.set_quantity('friction', 0.) # Zero friction 1379 domain.set_quantity('stage', initial_runup_height)# Constant negative initial stage1380 1329 1330 # Use function for elevation 1331 domain.set_quantity('elevation', topography) 1332 domain.set_quantity('friction', 0.) # Zero friction 1333 # Constant negative initial stage 1334 domain.set_quantity('stage', initial_runup_height) 1381 1335 1382 1336 #-------------------------------------------------------------- 1383 1337 # Setup boundary conditions 1384 1338 #-------------------------------------------------------------- 1385 Br = Reflective_boundary(domain) # Reflective wall 1386 Bd = Dirichlet_boundary([final_runup_height, # Constant inflow 1387 0, 1388 0]) 1389 1390 # All reflective to begin with (still water) 1339 Br = Reflective_boundary(domain) # Reflective wall 1340 Bd = Dirichlet_boundary([final_runup_height, 0, 0]) # Constant inflow 1341 1342 # All reflective to begin with (still water) 1391 1343 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 1392 1393 1344 1394 1345 #-------------------------------------------------------------- 1395 1346 # Test initial inundation height 1396 1347 #-------------------------------------------------------------- 1397 1398 1348 indices = domain.get_wet_elements() 1399 1349 z = domain.get_quantity('elevation').\ 1400 get_values(location='centroids', indices=indices)1401 assert num.alltrue(z <initial_runup_height)1350 get_values(location='centroids', indices=indices) 1351 assert num.alltrue(z < initial_runup_height) 1402 1352 1403 1353 q_ref = domain.get_maximum_inundation_elevation() 1404 assert num.allclose(q_ref, initial_runup_height, rtol = 1.0/N) # First order accuracy1405 1406 1354 # First order accuracy 1355 assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N) 1356 1407 1357 #-------------------------------------------------------------- 1408 1358 # Let triangles adjust … … 1411 1361 pass 1412 1362 1413 1414 1363 #-------------------------------------------------------------- 1415 1364 # Test inundation height again 1416 1365 #-------------------------------------------------------------- 1417 1418 1366 q_ref = domain.get_maximum_inundation_elevation() 1419 1367 q = get_maximum_inundation_elevation('runup_test.sww') 1420 msg = 'We got %f, should have been %f' % (q, q_ref)1368 msg = 'We got %f, should have been %f' % (q, q_ref) 1421 1369 assert num.allclose(q, q_ref, rtol=1.0/N), msg 1422 1370 1423 1424 1371 q = get_maximum_inundation_elevation('runup_test.sww') 1425 msg = 'We got %f, should have been %f' %(q, initial_runup_height) 1426 assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg 1427 1372 msg = 'We got %f, should have been %f' % (q, initial_runup_height) 1373 assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg 1428 1374 1429 1375 # Test error condition if time interval is out … … 1439 1385 # Check correct time interval 1440 1386 q, loc = get_maximum_inundation_data('runup_test.sww', 1441 time_interval=[0.0, 3.0]) 1442 msg = 'We got %f, should have been %f' % (q, initial_runup_height)1387 time_interval=[0.0, 3.0]) 1388 msg = 'We got %f, should have been %f' % (q, initial_runup_height) 1443 1389 assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg 1444 assert num.allclose(-loc[0]/2, q) # From topography formula 1445 1390 assert num.allclose(-loc[0]/2, q) # From topography formula 1446 1391 1447 1392 #-------------------------------------------------------------- … … 1450 1395 domain.set_boundary({'right': Bd}) 1451 1396 1452 1453 1397 #-------------------------------------------------------------- 1454 1398 # Evolve system through time … … 1456 1400 q_max = None 1457 1401 for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0, 1458 skip_initial_step = True): 1402 skip_initial_step = True): 1459 1403 q = domain.get_maximum_inundation_elevation() 1460 if q > q_max: q_max = q1461 1462 1404 if q > q_max: 1405 q_max = q 1406 1463 1407 #-------------------------------------------------------------- 1464 1408 # Test inundation height again 1465 1409 #-------------------------------------------------------------- 1466 1467 1410 indices = domain.get_wet_elements() 1468 1411 z = domain.get_quantity('elevation').\ 1469 get_values(location='centroids', indices=indices)1470 1471 assert num.alltrue(z <final_runup_height)1412 get_values(location='centroids', indices=indices) 1413 1414 assert num.alltrue(z < final_runup_height) 1472 1415 1473 1416 q = domain.get_maximum_inundation_elevation() 1474 assert num.allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy 1475 1476 q, loc = get_maximum_inundation_data('runup_test.sww', time_interval=[3.0, 3.0]) 1477 msg = 'We got %f, should have been %f' %(q, final_runup_height) 1417 # First order accuracy 1418 assert num.allclose(q, final_runup_height, rtol=1.0/N) 1419 1420 q, loc = get_maximum_inundation_data('runup_test.sww', 1421 time_interval=[3.0, 3.0]) 1422 msg = 'We got %f, should have been %f' % (q, final_runup_height) 1478 1423 assert num.allclose(q, final_runup_height, rtol=1.0/N), msg 1479 #print 'loc', loc, q 1480 assert num.allclose(-loc[0]/2, q) # From topography formula 1424 assert num.allclose(-loc[0]/2, q) # From topography formula 1481 1425 1482 1426 q = get_maximum_inundation_elevation('runup_test.sww') 1483 loc = get_maximum_inundation_location('runup_test.sww') 1484 msg = 'We got %f, should have been %f' % (q, q_max)1427 loc = get_maximum_inundation_location('runup_test.sww') 1428 msg = 'We got %f, should have been %f' % (q, q_max) 1485 1429 assert num.allclose(q, q_max, rtol=1.0/N), msg 1486 #print 'loc', loc, q 1487 assert num.allclose(-loc[0]/2, q) # From topography formula 1488 1489 1490 1491 q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[0, 3]) 1492 msg = 'We got %f, should have been %f' %(q, q_max) 1430 assert num.allclose(-loc[0]/2, q) # From topography formula 1431 1432 q = get_maximum_inundation_elevation('runup_test.sww', 1433 time_interval=[0, 3]) 1434 msg = 'We got %f, should have been %f' % (q, q_max) 1493 1435 assert num.allclose(q, q_max, rtol=1.0/N), msg 1494 1436 1495 1496 1437 # Check polygon mode 1497 polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]] # Runup region 1438 # Runup region 1439 polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]] 1498 1440 q = get_maximum_inundation_elevation('runup_test.sww', 1499 1441 polygon = polygon, 1500 1442 time_interval=[0, 3]) 1501 msg = 'We got %f, should have been %f' % (q, q_max)1443 msg = 'We got %f, should have been %f' % (q, q_max) 1502 1444 assert num.allclose(q, q_max, rtol=1.0/N), msg 1503 1445 1504 1505 polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]] # Offshore region1446 # Offshore region 1447 polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]] 1506 1448 q, loc = get_maximum_inundation_data('runup_test.sww', 1507 1449 polygon = polygon, 1508 1450 time_interval=[0, 3]) 1509 msg = 'We got %f, should have been %f' % (q, -0.475)1451 msg = 'We got %f, should have been %f' % (q, -0.475) 1510 1452 assert num.allclose(q, -0.475, rtol=1.0/N), msg 1511 1453 assert is_inside_polygon(loc, polygon) 1512 assert num.allclose(-loc[0]/2, q) # From topography formula1513 1514 1515 polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]] # Dry region1454 assert num.allclose(-loc[0]/2, q) # From topography formula 1455 1456 # Dry region 1457 polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]] 1516 1458 q, loc = get_maximum_inundation_data('runup_test.sww', 1517 1459 polygon = polygon, 1518 1460 time_interval=[0, 3]) 1519 msg = 'We got %s, should have been None' % (q)1461 msg = 'We got %s, should have been None' % (q) 1520 1462 assert q is None, msg 1521 msg = 'We got %s, should have been None' % (loc)1522 assert loc is None, msg 1463 msg = 'We got %s, should have been None' % (loc) 1464 assert loc is None, msg 1523 1465 1524 1466 # Check what happens if no time point is within interval 1525 1467 try: 1526 q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[2.75, 2.75]) 1468 q = get_maximum_inundation_elevation('runup_test.sww', 1469 time_interval=[2.75, 2.75]) 1527 1470 except AssertionError: 1528 1471 pass 1529 1472 else: 1530 1473 msg = 'Time interval should have raised an exception' 1531 raise msg1474 raise Exception, msg 1532 1475 1533 1476 # Cleanup … … 1537 1480 pass 1538 1481 #FIXME(Ole): Windows won't allow removal of this 1539 1540 1541 1482 1542 1483 def test_get_flow_through_cross_section_with_geo(self): 1543 1484 """test_get_flow_through_cross_section(self): … … 1545 1486 Test that the total flow through a cross section can be 1546 1487 correctly obtained at run-time from the ANUGA domain. 1547 1488 1548 1489 This test creates a flat bed with a known flow through it and tests 1549 1490 that the function correctly returns the expected flow. … … 1557 1498 q = u*h*w = 12 m^3/s 1558 1499 1559 1560 1500 This run tries it with georeferencing and with elevation = -1 1561 1562 1501 """ 1563 1502 1564 1503 import time, os 1565 1504 from Scientific.IO.NetCDF import NetCDFFile 1566 1567 # Setup1568 1505 from mesh_factory import rectangular 1569 1506 … … 1572 1509 length = 20 1573 1510 t_end = 1 1574 points, vertices, boundary = rectangular(length, width, 1575 length, width) 1511 points, vertices, boundary = rectangular(length, width, length, width) 1576 1512 1577 1513 # Create shallow water domain 1578 1514 domain = Domain(points, vertices, boundary, 1579 geo_reference=Geo_reference(56, 308500,6189000))1515 geo_reference=Geo_reference(56, 308500, 6189000)) 1580 1516 1581 1517 domain.default_order = 2 1582 1518 domain.set_quantities_to_be_stored(None) 1583 1584 1519 1585 1520 e = -1.0 … … 1589 1524 uh = u*h 1590 1525 1591 Br = Reflective_boundary(domain) # Side walls 1592 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet: 1593 1526 Br = Reflective_boundary(domain) # Side walls 1527 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet: 1594 1528 1595 1529 # Initial conditions … … 1597 1531 domain.set_quantity('stage', w) 1598 1532 domain.set_quantity('xmomentum', uh) 1599 domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br}) 1600 1601 1533 domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br}) 1534 1602 1535 # Interpolation points down the middle 1603 1536 I = [[0, width/2.], 1604 1537 [length/2., width/2.], 1605 1538 [length, width/2.]] 1606 interpolation_points = domain.geo_reference.get_absolute(I) 1607 1539 interpolation_points = domain.geo_reference.get_absolute(I) 1540 1608 1541 # Shortcuts to quantites 1609 stage = domain.get_quantity('stage') 1610 xmomentum = domain.get_quantity('xmomentum') 1611 ymomentum = domain.get_quantity('ymomentum') 1612 1613 for t in domain.evolve(yieldstep=0.1, finaltime =t_end):1542 stage = domain.get_quantity('stage') 1543 xmomentum = domain.get_quantity('xmomentum') 1544 ymomentum = domain.get_quantity('ymomentum') 1545 1546 for t in domain.evolve(yieldstep=0.1, finaltime=t_end): 1614 1547 # Check that quantities are they should be in the interior 1615 1616 w_t = stage.get_values(interpolation_points) 1548 w_t = stage.get_values(interpolation_points) 1617 1549 uh_t = xmomentum.get_values(interpolation_points) 1618 vh_t = ymomentum.get_values(interpolation_points) 1619 1550 vh_t = ymomentum.get_values(interpolation_points) 1551 1620 1552 assert num.allclose(w_t, w) 1621 assert num.allclose(uh_t, uh) 1622 assert num.allclose(vh_t, 0.0) 1623 1624 1553 assert num.allclose(uh_t, uh) 1554 assert num.allclose(vh_t, 0.0) 1555 1625 1556 # Check flows through the middle 1626 1557 for i in range(5): 1627 x = length/2. + i*0.23674563 # Arbitrary1558 x = length/2. + i*0.23674563 # Arbitrary 1628 1559 cross_section = [[x, 0], [x, width]] 1629 1560 1630 cross_section = domain.geo_reference.get_absolute(cross_section) 1561 cross_section = domain.geo_reference.get_absolute(cross_section) 1631 1562 Q = domain.get_flow_through_cross_section(cross_section, 1632 1563 verbose=False) … … 1634 1565 assert num.allclose(Q, uh*width) 1635 1566 1636 1637 1638 1567 def test_get_energy_through_cross_section_with_geo(self): 1639 1568 """test_get_energy_through_cross_section(self): … … 1641 1570 Test that the total and specific energy through a cross section can be 1642 1571 correctly obtained at run-time from the ANUGA domain. 1643 1572 1644 1573 This test creates a flat bed with a known flow through it and tests 1645 1574 that the function correctly returns the expected energies. … … 1653 1582 q = u*h*w = 12 m^3/s 1654 1583 1655 1656 1584 This run tries it with georeferencing and with elevation = -1 1657 1658 1585 """ 1659 1586 1660 1587 import time, os 1661 1588 from Scientific.IO.NetCDF import NetCDFFile 1662 1663 # Setup1664 1589 from mesh_factory import rectangular 1665 1590 … … 1668 1593 length = 20 1669 1594 t_end = 1 1670 points, vertices, boundary = rectangular(length, width, 1671 length, width) 1595 points, vertices, boundary = rectangular(length, width, length, width) 1672 1596 1673 1597 # Create shallow water domain 1674 1598 domain = Domain(points, vertices, boundary, 1675 geo_reference=Geo_reference(56, 308500,6189000))1599 geo_reference=Geo_reference(56, 308500, 6189000)) 1676 1600 1677 1601 domain.default_order = 2 1678 1602 domain.set_quantities_to_be_stored(None) 1679 1680 1603 1681 1604 e = -1.0 … … 1685 1608 uh = u*h 1686 1609 1687 Br = Reflective_boundary(domain) # Side walls 1688 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet: 1689 1610 Br = Reflective_boundary(domain) # Side walls 1611 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet: 1690 1612 1691 1613 # Initial conditions … … 1693 1615 domain.set_quantity('stage', w) 1694 1616 domain.set_quantity('xmomentum', uh) 1695 domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br}) 1696 1697 1617 domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br}) 1618 1698 1619 # Interpolation points down the middle 1699 1620 I = [[0, width/2.], 1700 1621 [length/2., width/2.], 1701 1622 [length, width/2.]] 1702 interpolation_points = domain.geo_reference.get_absolute(I) 1703 1623 interpolation_points = domain.geo_reference.get_absolute(I) 1624 1704 1625 # Shortcuts to quantites 1705 stage = domain.get_quantity('stage') 1706 xmomentum = domain.get_quantity('xmomentum') 1707 ymomentum = domain.get_quantity('ymomentum') 1708 1709 for t in domain.evolve(yieldstep=0.1, finaltime =t_end):1626 stage = domain.get_quantity('stage') 1627 xmomentum = domain.get_quantity('xmomentum') 1628 ymomentum = domain.get_quantity('ymomentum') 1629 1630 for t in domain.evolve(yieldstep=0.1, finaltime=t_end): 1710 1631 # Check that quantities are they should be in the interior 1711 1712 w_t = stage.get_values(interpolation_points) 1632 w_t = stage.get_values(interpolation_points) 1713 1633 uh_t = xmomentum.get_values(interpolation_points) 1714 vh_t = ymomentum.get_values(interpolation_points) 1715 1634 vh_t = ymomentum.get_values(interpolation_points) 1635 1716 1636 assert num.allclose(w_t, w) 1717 assert num.allclose(uh_t, uh) 1718 assert num.allclose(vh_t, 0.0) 1719 1720 1637 assert num.allclose(uh_t, uh) 1638 assert num.allclose(vh_t, 0.0) 1639 1721 1640 # Check energies through the middle 1722 1641 for i in range(5): 1723 x = length/2. + i*0.23674563 # Arbitrary1642 x = length/2. + i*0.23674563 # Arbitrary 1724 1643 cross_section = [[x, 0], [x, width]] 1725 1644 1726 cross_section = domain.geo_reference.get_absolute(cross_section) 1645 cross_section = domain.geo_reference.get_absolute(cross_section) 1727 1646 Es = domain.get_energy_through_cross_section(cross_section, 1728 1647 kind='specific', 1729 1648 verbose=False) 1730 1649 1731 1650 assert num.allclose(Es, h + 0.5*u*u/g) 1732 1651 1733 1652 Et = domain.get_energy_through_cross_section(cross_section, 1734 1653 kind='total', 1735 1654 verbose=False) 1736 assert num.allclose(Et, w + 0.5*u*u/g) 1737 1738 1739 1740 1655 assert num.allclose(Et, w + 0.5*u*u/g) 1741 1656 1742 1657 def test_another_runup_example(self): … … 1747 1662 """ 1748 1663 1749 #-----------------------------------------------------------------1750 # Import necessary modules1751 #-----------------------------------------------------------------1752 1753 1664 from anuga.pmesh.mesh_interface import create_mesh_from_regions 1754 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 1665 from anuga.abstract_2d_finite_volumes.mesh_factory \ 1666 import rectangular_cross 1755 1667 from anuga.shallow_water import Domain 1756 1668 from anuga.shallow_water import Reflective_boundary 1757 1669 from anuga.shallow_water import Dirichlet_boundary 1758 1670 1759 1760 1671 #----------------------------------------------------------------- 1761 1672 # Setup computational domain 1762 1673 #----------------------------------------------------------------- 1763 points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh1764 domain = Domain(points, vertices, boundary) # Create domain1674 points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh 1675 domain = Domain(points, vertices, boundary) # Create domain 1765 1676 domain.set_quantities_to_be_stored(None) 1766 1677 domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this 1767 1678 1768 1679 # FIXME (Ole): Need tests where this is commented out 1769 domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7) 1770 domain.H0 = 0 # Backwards compatibility (6/2/7) 1771 domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8) 1772 1680 domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7) 1681 domain.H0 = 0 # Backwards compatibility (6/2/7) 1682 domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8) 1773 1683 1774 1684 #----------------------------------------------------------------- 1775 1685 # Setup initial conditions 1776 1686 #----------------------------------------------------------------- 1777 1778 def topography(x,y): 1687 def topography(x, y): 1779 1688 return -x/2 # linear bed slope 1780 1689 1781 domain.set_quantity('elevation', topography) 1782 domain.set_quantity('friction', 0.0) 1783 domain.set_quantity('stage', expression='elevation') 1784 1690 domain.set_quantity('elevation', topography) 1691 domain.set_quantity('friction', 0.0) 1692 domain.set_quantity('stage', expression='elevation') 1785 1693 1786 1694 #---------------------------------------------------------------- 1787 1695 # Setup boundary conditions 1788 1696 #---------------------------------------------------------------- 1789 1790 Br = Reflective_boundary(domain) # Solid reflective wall 1791 Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values 1697 Br = Reflective_boundary(domain) # Solid reflective wall 1698 Bd = Dirichlet_boundary([-0.2, 0., 0.]) # Constant boundary values 1792 1699 domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br}) 1793 1794 1700 1795 1701 #---------------------------------------------------------------- 1796 1702 # Evolve system through time 1797 1703 #---------------------------------------------------------------- 1798 1799 1704 interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]] 1800 1705 gauge_values = [] … … 1803 1708 1804 1709 time = [] 1805 for t in domain.evolve(yieldstep = 0.1, finaltime =5.0):1710 for t in domain.evolve(yieldstep=0.1, finaltime=5.0): 1806 1711 # Record time series at known points 1807 1712 time.append(domain.get_time()) 1808 1713 1809 1714 stage = domain.get_quantity('stage') 1810 1715 w = stage.get_values(interpolation_points=interpolation_points) 1811 1716 1812 1717 for i, _ in enumerate(interpolation_points): 1813 1718 gauge_values[i].append(w[i]) 1814 1719 1815 1816 #print1817 #print time1818 #print1819 #for i, (x,y) in enumerate(interpolation_points):1820 # print i, gauge_values[i]1821 # print1822 1823 1720 #Reference (nautilus 26/6/2008) 1824 1825 G0 = [-0.20000000000000001, -0.20000000000000001, -0.19920600846161715, -0.19153647344085376, -0.19127622768281194, -0.1770671909675095, -0.16739412133181927, -0.16196038919122191, -0.15621633053131384, -0.15130021599977705, -0.13930978857215484, -0.19349274358263582, -0.19975307598803765, -0.19999897143103357, -0.1999999995532111, -0.19999999999949952, -0.19999999999949952, -0.19999999999949952, -0.19997270012494556, -0.19925805948554556, -0.19934513778450533, -0.19966484196394893, -0.1997352860102084, -0.19968260481750394, -0.19980280797303882, -0.19998804881822749, -0.19999999778075916, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167] 1826 1827 G1 = [-0.29999999999999993, -0.29999588068034899, -0.29250047332330331, -0.28335081844518584, -0.26142206997410805, -0.22656028856329835, -0.21224087216745585, -0.19934324109114465, -0.1889857939783175, -0.18146311603911383, -0.17401078727434263, -0.15419361061257214, -0.16225060576782063, -0.19010941396999181, -0.20901161407004412, -0.21670683975774699, -0.21771386270738891, -0.21481284465869752, -0.21063120869004387, -0.20669243364582401, -0.20320707386714859, -0.19984087691926442, -0.19725417448019505, -0.19633783049069981, -0.19650494599999785, -0.19708316838336942, -0.19779309449413818, -0.19853070294429562, -0.19917342167307153, -0.19964814677795845, -0.19991627610824922, -0.20013162970144974, -0.20029864969405509, -0.20036259676501131, -0.20030682824965193, -0.20016105135750167, -0.19997664501985918, -0.19980185871568762, -0.19966836175417696, -0.19958856744312226, -0.19955954696194517, -0.19956950051110917, -0.19960377086336181, -0.19964885299433241, -0.19969427478531132, -0.19973301547655564, -0.19976121574277764, -0.19977765285688653, -0.19978315117522441, -0.19977994634841739, -0.19977101394878494] 1828 1829 G2 = [-0.40000000000000002, -0.39077401254732241, -0.33350466136630474, -0.29771023004255281, -0.27605439066140897, -0.25986156218997497, -0.24502185018573647, -0.231792624329521, -0.21981564668803993, -0.20870707082936543, -0.19877739883776599, -0.18980922837977957, -0.17308011674005838, -0.16306400164013773, -0.17798470933304333, -0.1929554075869116, -0.20236705191987037, -0.20695767560655007, -0.20841025876092567, -0.20792102174869989, -0.20655350005579293, -0.20492002526259828, -0.20310627026780645, -0.20105983335287836, -0.19937394565794653, -0.19853917506699659, -0.19836389977624452, -0.19850305023602796, -0.19877764028836831, -0.19910928131034669, -0.19943705712418805, -0.19970344172958865, -0.19991076989870474, -0.20010020127747646, -0.20025937787100062, -0.20035087292905965, -0.20035829921463297, -0.20029606557316171, -0.20019606915365515, -0.20009096093399206, -0.20000371608204368, -0.19994495432920584, -0.19991535665176338, -0.19990981826533513, -0.19992106419898723, -0.19994189853516578, -0.19996624091229293, -0.19998946016985167, -0.20000842303470234, -0.20002144460718174, -0.20002815561337187] 1830 1831 G3 = [-0.45000000000000001, -0.37631169657400332, -0.33000044342859486, -0.30586045469008522, -0.28843572253009941, -0.27215308978603808, -0.25712951540331219, -0.2431608296216613, -0.23032023651386374, -0.2184546873456619, -0.20735123704254332, -0.19740397194806389, -0.1859829564064375, -0.16675980728362105, -0.16951575032846536, -0.1832860872609344, -0.19485758939241243, -0.20231368291811427, -0.20625610376074754, -0.20758116241495619, -0.20721445402086161, -0.20603406830353785, -0.20450262808396991, -0.2026769581185151, -0.2007401212066364, -0.19931160535777592, -0.19863606301128725, -0.19848511940572691, -0.19860091042948352, -0.19885490669377764, -0.19916542732701112, -0.19946678238611959, -0.19971209594104697, -0.19991912886512292, -0.2001058430788881, -0.20024959409472989, -0.20032160254609382, -0.20031583165752354, -0.20025051539293123, -0.2001556115816068, -0.20005952955420872, -0.1999814429561611, -0.19992977821558131, -0.19990457708664208, -0.19990104785490476, -0.19991257153954825, -0.19993258231880562, -0.19995548502882532, -0.19997700760919687, -0.19999429663503748, -0.20000588800248761] 1832 1721 G0 = [-0.20000000000000001, -0.20000000000000001, -0.19920600846161715, 1722 -0.19153647344085376, -0.19127622768281194, -0.1770671909675095, 1723 -0.16739412133181927, -0.16196038919122191, -0.15621633053131384, 1724 -0.15130021599977705, -0.13930978857215484, -0.19349274358263582, 1725 -0.19975307598803765, -0.19999897143103357, -0.1999999995532111, 1726 -0.19999999999949952, -0.19999999999949952, -0.19999999999949952, 1727 -0.19997270012494556, -0.19925805948554556, -0.19934513778450533, 1728 -0.19966484196394893, -0.1997352860102084, -0.19968260481750394, 1729 -0.19980280797303882, -0.19998804881822749, -0.19999999778075916, 1730 -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, 1731 -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, 1732 -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, 1733 -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, 1734 -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, 1735 -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, 1736 -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, 1737 -0.19999999999966167, -0.19999999999966167, -0.19999999999966167] 1738 1739 G1 = [-0.29999999999999993, -0.29999588068034899, -0.29250047332330331, 1740 -0.28335081844518584, -0.26142206997410805, -0.22656028856329835, 1741 -0.21224087216745585, -0.19934324109114465, -0.1889857939783175, 1742 -0.18146311603911383, -0.17401078727434263, -0.15419361061257214, 1743 -0.16225060576782063, -0.19010941396999181, -0.20901161407004412, 1744 -0.21670683975774699, -0.21771386270738891, -0.21481284465869752, 1745 -0.21063120869004387, -0.20669243364582401, -0.20320707386714859, 1746 -0.19984087691926442, -0.19725417448019505, -0.19633783049069981, 1747 -0.19650494599999785, -0.19708316838336942, -0.19779309449413818, 1748 -0.19853070294429562, -0.19917342167307153, -0.19964814677795845, 1749 -0.19991627610824922, -0.20013162970144974, -0.20029864969405509, 1750 -0.20036259676501131, -0.20030682824965193, -0.20016105135750167, 1751 -0.19997664501985918, -0.19980185871568762, -0.19966836175417696, 1752 -0.19958856744312226, -0.19955954696194517, -0.19956950051110917, 1753 -0.19960377086336181, -0.19964885299433241, -0.19969427478531132, 1754 -0.19973301547655564, -0.19976121574277764, -0.19977765285688653, 1755 -0.19978315117522441, -0.19977994634841739, -0.19977101394878494] 1756 1757 G2 = [-0.40000000000000002, -0.39077401254732241, -0.33350466136630474, 1758 -0.29771023004255281, -0.27605439066140897, -0.25986156218997497, 1759 -0.24502185018573647, -0.231792624329521, -0.21981564668803993, 1760 -0.20870707082936543, -0.19877739883776599, -0.18980922837977957, 1761 -0.17308011674005838, -0.16306400164013773, -0.17798470933304333, 1762 -0.1929554075869116, -0.20236705191987037, -0.20695767560655007, 1763 -0.20841025876092567, -0.20792102174869989, -0.20655350005579293, 1764 -0.20492002526259828, -0.20310627026780645, -0.20105983335287836, 1765 -0.19937394565794653, -0.19853917506699659, -0.19836389977624452, 1766 -0.19850305023602796, -0.19877764028836831, -0.19910928131034669, 1767 -0.19943705712418805, -0.19970344172958865, -0.19991076989870474, 1768 -0.20010020127747646, -0.20025937787100062, -0.20035087292905965, 1769 -0.20035829921463297, -0.20029606557316171, -0.20019606915365515, 1770 -0.20009096093399206, -0.20000371608204368, -0.19994495432920584, 1771 -0.19991535665176338, -0.19990981826533513, -0.19992106419898723, 1772 -0.19994189853516578, -0.19996624091229293, -0.19998946016985167, 1773 -0.20000842303470234, -0.20002144460718174, -0.20002815561337187] 1774 1775 G3 = [-0.45000000000000001, -0.37631169657400332, -0.33000044342859486, 1776 -0.30586045469008522, -0.28843572253009941, -0.27215308978603808, 1777 -0.25712951540331219, -0.2431608296216613, -0.23032023651386374, 1778 -0.2184546873456619, -0.20735123704254332, -0.19740397194806389, 1779 -0.1859829564064375, -0.16675980728362105, -0.16951575032846536, 1780 -0.1832860872609344, -0.19485758939241243, -0.20231368291811427, 1781 -0.20625610376074754, -0.20758116241495619, -0.20721445402086161, 1782 -0.20603406830353785, -0.20450262808396991, -0.2026769581185151, 1783 -0.2007401212066364, -0.19931160535777592, -0.19863606301128725, 1784 -0.19848511940572691, -0.19860091042948352, -0.19885490669377764, 1785 -0.19916542732701112, -0.19946678238611959, -0.19971209594104697, 1786 -0.19991912886512292, -0.2001058430788881, -0.20024959409472989, 1787 -0.20032160254609382, -0.20031583165752354, -0.20025051539293123, 1788 -0.2001556115816068, -0.20005952955420872, -0.1999814429561611, 1789 -0.19992977821558131, -0.19990457708664208, -0.19990104785490476, 1790 -0.19991257153954825, -0.19993258231880562, -0.19995548502882532, 1791 -0.19997700760919687, -0.19999429663503748, -0.20000588800248761] 1792 1833 1793 #FIXME (DSG):This is a hack so the anuga install, not precompiled 1834 1794 # works on DSG's win2000, python 2.3 … … 1841 1801 #if len(gauge_values[3]) == 52: gauge_values[3].pop() 1842 1802 1843 ## print len(G0), len(gauge_values[0])1844 ## print len(G1), len(gauge_values[1])1845 1846 #print gauge_values[3]1847 #print G0[:4]1848 #print array(gauge_values[0])-array(G0)1849 1850 1851 1803 assert num.allclose(gauge_values[0], G0) 1852 1804 assert num.allclose(gauge_values[1], G1) 1853 1805 assert num.allclose(gauge_values[2], G2) 1854 assert num.allclose(gauge_values[3], G3) 1855 1856 1857 1858 1859 1860 1806 assert num.allclose(gauge_values[3], G3) 1861 1807 1862 1808 ##################################################### … … 1864 1810 def test_flux_optimisation(self): 1865 1811 """test_flux_optimisation 1812 1866 1813 Test that fluxes are correctly computed using 1867 1814 dry and still cell exclusions … … 1879 1826 1880 1827 points = [a, b, c, d, e, f] 1881 # bac, bce, ecf,dbe1882 vertices = [ 1828 # bac, bce, ecf, dbe 1829 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1883 1830 1884 1831 domain = Domain(points, vertices) … … 1889 1836 1890 1837 h = 0.1 1891 def stage(x, y):1892 return slope(x, y)+h1838 def stage(x, y): 1839 return slope(x, y) + h 1893 1840 1894 1841 domain.set_quantity('elevation', slope) 1895 1842 domain.set_quantity('stage', stage) 1896 1843 1897 # Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?) 1898 domain.distribute_to_vertices_and_edges() 1844 # Allow slope limiters to work 1845 # (FIXME (Ole): Shouldn't this be automatic in ANUGA?) 1846 domain.distribute_to_vertices_and_edges() 1899 1847 1900 1848 initial_stage = copy.copy(domain.quantities['stage'].vertex_values) … … 1902 1850 domain.set_boundary({'exterior': Reflective_boundary(domain)}) 1903 1851 1904 1905 # Check that update arrays are initialised to zero 1852 # Check that update arrays are initialised to zero 1906 1853 assert num.allclose(domain.get_quantity('stage').explicit_update, 0) 1907 1854 assert num.allclose(domain.get_quantity('xmomentum').explicit_update, 0) 1908 assert num.allclose(domain.get_quantity('ymomentum').explicit_update, 0) 1909 1855 assert num.allclose(domain.get_quantity('ymomentum').explicit_update, 0) 1910 1856 1911 1857 # Get true values … … 1914 1860 stage_ref = copy.copy(domain.get_quantity('stage').explicit_update) 1915 1861 xmom_ref = copy.copy(domain.get_quantity('xmomentum').explicit_update) 1916 ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update) 1862 ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update) 1917 1863 1918 1864 # Try with flux optimisation … … 1920 1866 domain.compute_fluxes() 1921 1867 1922 assert num.allclose(stage_ref, domain.get_quantity('stage').explicit_update) 1923 assert num.allclose(xmom_ref, domain.get_quantity('xmomentum').explicit_update) 1924 assert num.allclose(ymom_ref, domain.get_quantity('ymomentum').explicit_update) 1925 1926 1927 1868 assert num.allclose(stage_ref, 1869 domain.get_quantity('stage').explicit_update) 1870 assert num.allclose(xmom_ref, 1871 domain.get_quantity('xmomentum').explicit_update) 1872 assert num.allclose(ymom_ref, 1873 domain.get_quantity('ymomentum').explicit_update) 1874 1928 1875 def test_initial_condition(self): 1929 1876 """test_initial_condition 1877 1930 1878 Test that initial condition is output at time == 0 and that 1931 1879 computed values change as system evolves … … 1943 1891 1944 1892 points = [a, b, c, d, e, f] 1945 # bac, bce, ecf,dbe1946 vertices = [ 1893 # bac, bce, ecf, dbe 1894 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1947 1895 1948 1896 domain = Domain(points, vertices) … … 1953 1901 1954 1902 h = 0.1 1955 def stage(x, y):1956 return slope(x, y)+h1903 def stage(x, y): 1904 return slope(x, y) + h 1957 1905 1958 1906 domain.set_quantity('elevation', slope) 1959 1907 domain.set_quantity('stage', stage) 1960 1908 1961 # Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?) 1962 domain.distribute_to_vertices_and_edges() 1909 # Allow slope limiters to work 1910 # (FIXME (Ole): Shouldn't this be automatic in ANUGA?) 1911 domain.distribute_to_vertices_and_edges() 1963 1912 1964 1913 initial_stage = copy.copy(domain.quantities['stage'].vertex_values) … … 1967 1916 1968 1917 domain.optimise_dry_cells = True 1918 1969 1919 #Evolution 1970 for t in domain.evolve(yieldstep = 0.5, finaltime =2.0):1920 for t in domain.evolve(yieldstep=0.5, finaltime=2.0): 1971 1921 stage = domain.quantities['stage'].vertex_values 1972 1922 … … 1976 1926 assert not num.allclose(stage, initial_stage) 1977 1927 1978 1979 1928 os.remove(domain.get_name() + '.sww') 1980 1929 1981 1982 1983 1930 ##################################################### 1931 1984 1932 def test_gravity(self): 1985 1933 #Assuming no friction … … 1995 1943 1996 1944 points = [a, b, c, d, e, f] 1997 # bac, bce, ecf,dbe1998 vertices = [ 1945 # bac, bce, ecf, dbe 1946 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1999 1947 2000 1948 domain = Domain(points, vertices) … … 2005 1953 2006 1954 h = 0.1 2007 def stage(x, y):2008 return slope(x, y)+h1955 def stage(x, y): 1956 return slope(x, y) + h 2009 1957 2010 1958 domain.set_quantity('elevation', slope) … … 2018 1966 2019 1967 assert num.allclose(domain.quantities['stage'].explicit_update, 0) 2020 assert num.allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3) 1968 assert num.allclose(domain.quantities['xmomentum'].explicit_update, 1969 -g*h*3) 2021 1970 assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0) 2022 2023 1971 2024 1972 def test_manning_friction(self): … … 2033 1981 2034 1982 points = [a, b, c, d, e, f] 2035 # bac, bce, ecf,dbe2036 vertices = [ 1983 # bac, bce, ecf, dbe 1984 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2037 1985 2038 1986 domain = Domain(points, vertices) … … 2043 1991 2044 1992 h = 0.1 2045 def stage(x, y):2046 return slope(x, y)+h1993 def stage(x, y): 1994 return slope(x, y) + h 2047 1995 2048 1996 eta = 0.07 … … 2058 2006 2059 2007 assert num.allclose(domain.quantities['stage'].explicit_update, 0) 2060 assert num.allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3) 2008 assert num.allclose(domain.quantities['xmomentum'].explicit_update, 2009 -g*h*3) 2061 2010 assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0) 2062 2011 2063 2012 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2064 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 0) 2065 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 0) 2013 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2014 0) 2015 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2016 0) 2066 2017 2067 2018 #Create some momentum for friction to work with 2068 2019 domain.set_quantity('xmomentum', 1) 2069 S = -g *eta**2 / h**(7.0/3)2020 S = -g*eta**2 / h**(7.0/3) 2070 2021 2071 2022 domain.compute_forcing_terms() 2072 2023 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2073 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, S) 2074 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 0) 2024 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2025 S) 2026 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2027 0) 2075 2028 2076 2029 #A more complex example … … 2082 2035 domain.set_quantity('ymomentum', 4) 2083 2036 2084 S = -g * eta**2 * 5 / h**(7.0/3) 2085 2037 S = -g*eta**2*5 / h**(7.0/3) 2086 2038 2087 2039 domain.compute_forcing_terms() 2088 2040 2089 2041 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2090 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S) 2091 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S) 2042 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2043 3*S) 2044 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2045 4*S) 2092 2046 2093 2047 def test_constant_wind_stress(self): … … 2103 2057 2104 2058 points = [a, b, c, d, e, f] 2105 #bac, bce, ecf, dbe 2106 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2107 2059 # bac, bce, ecf, dbe 2060 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2108 2061 2109 2062 domain = Domain(points, vertices) … … 2121 2074 phi = 135 2122 2075 domain.forcing_terms = [] 2123 domain.forcing_terms.append( Wind_stress(s, phi))2076 domain.forcing_terms.append(Wind_stress(s, phi)) 2124 2077 2125 2078 domain.compute_forcing_terms() 2126 2079 2127 2128 const = eta_w*rho_a/rho_w 2080 const = eta_w*rho_a / rho_w 2129 2081 2130 2082 #Convert to radians 2131 phi = phi*pi /1802083 phi = phi*pi / 180 2132 2084 2133 2085 #Compute velocity vector (u, v) … … 2141 2093 assert num.allclose(domain.quantities['xmomentum'].explicit_update, S*u) 2142 2094 assert num.allclose(domain.quantities['ymomentum'].explicit_update, S*v) 2143 2144 2095 2145 2096 def test_variable_wind_stress(self): … … 2155 2106 2156 2107 points = [a, b, c, d, e, f] 2157 # bac, bce, ecf,dbe2158 vertices = [ 2108 # bac, bce, ecf, dbe 2109 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2159 2110 2160 2111 domain = Domain(points, vertices) … … 2168 2119 domain.set_boundary({'exterior': Br}) 2169 2120 2170 2171 domain.time = 5.54 #Take a random time (not zero) 2121 domain.time = 5.54 # Take a random time (not zero) 2172 2122 2173 2123 #Setup only one forcing term, constant wind stress … … 2175 2125 phi = 135 2176 2126 domain.forcing_terms = [] 2177 domain.forcing_terms.append( Wind_stress(s = speed, phi = angle))2127 domain.forcing_terms.append(Wind_stress(s=speed, phi=angle)) 2178 2128 2179 2129 domain.compute_forcing_terms() 2180 2130 2181 2131 #Compute reference solution 2182 const = eta_w*rho_a /rho_w2183 2184 N = len(domain) # number_of_triangles2132 const = eta_w*rho_a / rho_w 2133 2134 N = len(domain) # number_of_triangles 2185 2135 2186 2136 xc = domain.get_centroid_coordinates() … … 2192 2142 phi_vec = angle(t,x,y) 2193 2143 2194 2195 2144 for k in range(N): 2196 # Convert to radians2197 phi = phi_vec[k]*pi /1802145 # Convert to radians 2146 phi = phi_vec[k]*pi / 180 2198 2147 s = s_vec[k] 2199 2148 2200 # Compute velocity vector (u, v)2149 # Compute velocity vector (u, v) 2201 2150 u = s*cos(phi) 2202 2151 v = s*sin(phi) 2203 2152 2204 # Compute wind stress2153 # Compute wind stress 2205 2154 S = const * num.sqrt(u**2 + v**2) 2206 2155 2207 assert num.allclose(domain.quantities['stage'].explicit_update[k], 0)2208 assert num.allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)2209 assert num.allclose(domain.quantities[' ymomentum'].explicit_update[k], S*v)2210 2211 2212 2213 2214 2156 assert num.allclose(domain.quantities['stage'].explicit_update[k], 2157 0) 2158 assert num.allclose(domain.quantities['xmomentum'].\ 2159 explicit_update[k], 2160 S*u) 2161 assert num.allclose(domain.quantities['ymomentum'].\ 2162 explicit_update[k], 2163 S*v) 2215 2164 2216 2165 def test_windfield_from_file(self): 2166 import time 2217 2167 from anuga.config import rho_a, rho_w, eta_w 2218 2168 from math import pi, cos, sin 2219 2169 from anuga.config import time_format 2220 2170 from anuga.abstract_2d_finite_volumes.util import file_function 2221 import time2222 2223 2171 2224 2172 a = [0.0, 0.0] … … 2230 2178 2231 2179 points = [a, b, c, d, e, f] 2232 # bac, bce, ecf,dbe2233 vertices = [ 2180 # bac, bce, ecf, dbe 2181 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2234 2182 2235 2183 domain = Domain(points, vertices) 2236 2184 2237 # Flat surface with 1m of water2185 # Flat surface with 1m of water 2238 2186 domain.set_quantity('elevation', 0) 2239 2187 domain.set_quantity('stage', 1.0) … … 2243 2191 domain.set_boundary({'exterior': Br}) 2244 2192 2245 2246 domain.time = 7 #Take a time that is represented in file (not zero) 2247 2248 #Write wind stress file (ensure that domain.time is covered) 2249 #Take x=1 and y=0 2193 domain.time = 7 # Take a time that is represented in file (not zero) 2194 2195 # Write wind stress file (ensure that domain.time is covered) 2196 # Take x=1 and y=0 2250 2197 filename = 'test_windstress_from_file' 2251 2198 start = time.mktime(time.strptime('2000', '%Y')) 2252 2199 fid = open(filename + '.txt', 'w') 2253 dt = 1 #One second interval2200 dt = 1 # One second interval 2254 2201 t = 0.0 2255 2202 while t <= 10.0: 2256 2203 t_string = time.strftime(time_format, time.gmtime(t+start)) 2257 2204 2258 fid.write('%s, %f %f\n' %(t_string, 2259 speed(t,[1],[0])[0], 2260 angle(t,[1],[0])[0])) 2205 fid.write('%s, %f %f\n' % 2206 (t_string, speed(t,[1],[0])[0], angle(t,[1],[0])[0])) 2261 2207 t += dt 2262 2208 2263 2209 fid.close() 2264 2210 2265 2266 #Convert ASCII file to NetCDF (Which is what we really like!)2267 from data_manager import timefile2netcdf 2211 # Convert ASCII file to NetCDF (Which is what we really like!) 2212 from data_manager import timefile2netcdf 2213 2268 2214 timefile2netcdf(filename) 2269 2215 os.remove(filename + '.txt') 2270 2216 2271 2272 #Setup wind stress 2273 F = file_function(filename + '.tms', quantities = ['Attribute0', 2274 'Attribute1']) 2217 # Setup wind stress 2218 F = file_function(filename + '.tms', 2219 quantities=['Attribute0', 'Attribute1']) 2275 2220 os.remove(filename + '.tms') 2276 2277 2278 #print 'F(5)', F(5) 2279 2280 #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3)) 2281 2282 #print dir(F) 2283 #print F.T 2284 #print F.precomputed_values 2285 # 2286 #F = file_function(filename + '.txt') 2287 # 2288 #print dir(F) 2289 #print F.T 2290 #print F.Q 2291 2221 2292 2222 W = Wind_stress(F) 2293 2223 … … 2297 2227 domain.compute_forcing_terms() 2298 2228 2299 # Compute reference solution2300 const = eta_w*rho_a /rho_w2301 2302 N = len(domain) # number_of_triangles2229 # Compute reference solution 2230 const = eta_w*rho_a / rho_w 2231 2232 N = len(domain) # number_of_triangles 2303 2233 2304 2234 t = domain.time 2305 2235 2306 s = speed(t,[1],[0])[0] 2307 phi = angle(t,[1],[0])[0] 2308 2309 #Convert to radians 2310 phi = phi*pi/180 2311 2312 2313 #Compute velocity vector (u, v) 2236 s = speed(t, [1], [0])[0] 2237 phi = angle(t, [1], [0])[0] 2238 2239 # Convert to radians 2240 phi = phi*pi / 180 2241 2242 # Compute velocity vector (u, v) 2314 2243 u = s*cos(phi) 2315 2244 v = s*sin(phi) 2316 2245 2317 # Compute wind stress2246 # Compute wind stress 2318 2247 S = const * num.sqrt(u**2 + v**2) 2319 2248 2320 2249 for k in range(N): 2321 assert num.allclose(domain.quantities['stage'].explicit_update[k], 0) 2322 assert num.allclose(domain.quantities['xmomentum'].explicit_update[k], S*u) 2323 assert num.allclose(domain.quantities['ymomentum'].explicit_update[k], S*v) 2324 2250 assert num.allclose(domain.quantities['stage'].explicit_update[k], 2251 0) 2252 assert num.allclose(domain.quantities['xmomentum'].\ 2253 explicit_update[k], 2254 S*u) 2255 assert num.allclose(domain.quantities['ymomentum'].\ 2256 explicit_update[k], 2257 S*v) 2325 2258 2326 2259 def test_windfield_from_file_seconds(self): 2260 import time 2327 2261 from anuga.config import rho_a, rho_w, eta_w 2328 2262 from math import pi, cos, sin 2329 2263 from anuga.config import time_format 2330 2264 from anuga.abstract_2d_finite_volumes.util import file_function 2331 import time2332 2333 2265 2334 2266 a = [0.0, 0.0] … … 2340 2272 2341 2273 points = [a, b, c, d, e, f] 2342 # bac, bce, ecf,dbe2343 vertices = [ 2274 # bac, bce, ecf, dbe 2275 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2344 2276 2345 2277 domain = Domain(points, vertices) 2346 2278 2347 # Flat surface with 1m of water2279 # Flat surface with 1m of water 2348 2280 domain.set_quantity('elevation', 0) 2349 2281 domain.set_quantity('stage', 1.0) … … 2353 2285 domain.set_boundary({'exterior': Br}) 2354 2286 2355 2356 domain.time = 7 #Take a time that is represented in file (not zero) 2357 2358 #Write wind stress file (ensure that domain.time is covered) 2359 #Take x=1 and y=0 2287 domain.time = 7 # Take a time that is represented in file (not zero) 2288 2289 # Write wind stress file (ensure that domain.time is covered) 2290 # Take x=1 and y=0 2360 2291 filename = 'test_windstress_from_file' 2361 2292 start = time.mktime(time.strptime('2000', '%Y')) 2362 2293 fid = open(filename + '.txt', 'w') 2363 dt = 0.5 #1 #Onesecond interval2294 dt = 0.5 # Half second interval 2364 2295 t = 0.0 2365 2296 while t <= 10.0: 2366 fid.write('%s, %f %f\n' %(str(t), 2367 speed(t,[1],[0])[0], 2368 angle(t,[1],[0])[0])) 2297 fid.write('%s, %f %f\n' 2298 % (str(t), speed(t, [1], [0])[0], angle(t, [1], [0])[0])) 2369 2299 t += dt 2370 2300 2371 2301 fid.close() 2372 2302 2373 2374 #Convert ASCII file to NetCDF (Which is what we really like!)2375 from data_manager import timefile2netcdf 2303 # Convert ASCII file to NetCDF (Which is what we really like!) 2304 from data_manager import timefile2netcdf 2305 2376 2306 timefile2netcdf(filename, time_as_seconds=True) 2377 2307 os.remove(filename + '.txt') 2378 2308 2379 2380 #Setup wind stress 2381 F = file_function(filename + '.tms', quantities = ['Attribute0', 2382 'Attribute1']) 2309 # Setup wind stress 2310 F = file_function(filename + '.tms', 2311 quantities=['Attribute0', 'Attribute1']) 2383 2312 os.remove(filename + '.tms') 2384 2385 2386 #print 'F(5)', F(5) 2387 2388 #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3)) 2389 2390 #print dir(F) 2391 #print F.T 2392 #print F.precomputed_values 2393 # 2394 #F = file_function(filename + '.txt') 2395 # 2396 #print dir(F) 2397 #print F.T 2398 #print F.Q 2399 2313 2400 2314 W = Wind_stress(F) 2401 2315 … … 2405 2319 domain.compute_forcing_terms() 2406 2320 2407 # Compute reference solution2408 const = eta_w*rho_a /rho_w2409 2410 N = len(domain) # number_of_triangles2321 # Compute reference solution 2322 const = eta_w*rho_a / rho_w 2323 2324 N = len(domain) # number_of_triangles 2411 2325 2412 2326 t = domain.time 2413 2327 2414 s = speed(t,[1],[0])[0] 2415 phi = angle(t,[1],[0])[0] 2416 2417 #Convert to radians 2418 phi = phi*pi/180 2419 2420 2421 #Compute velocity vector (u, v) 2328 s = speed(t, [1], [0])[0] 2329 phi = angle(t, [1], [0])[0] 2330 2331 # Convert to radians 2332 phi = phi*pi / 180 2333 2334 # Compute velocity vector (u, v) 2422 2335 u = s*cos(phi) 2423 2336 v = s*sin(phi) 2424 2337 2425 # Compute wind stress2338 # Compute wind stress 2426 2339 S = const * num.sqrt(u**2 + v**2) 2427 2340 2428 2341 for k in range(N): 2429 assert num.allclose(domain.quantities['stage'].explicit_update[k], 0) 2430 assert num.allclose(domain.quantities['xmomentum'].explicit_update[k], S*u) 2431 assert num.allclose(domain.quantities['ymomentum'].explicit_update[k], S*v) 2432 2433 2434 2342 assert num.allclose(domain.quantities['stage'].explicit_update[k], 2343 0) 2344 assert num.allclose(domain.quantities['xmomentum'].\ 2345 explicit_update[k], 2346 S*u) 2347 assert num.allclose(domain.quantities['ymomentum'].\ 2348 explicit_update[k], 2349 S*v) 2435 2350 2436 2351 def test_wind_stress_error_condition(self): … … 2439 2354 """ 2440 2355 2356 from math import pi, cos, sin 2441 2357 from anuga.config import rho_a, rho_w, eta_w 2442 from math import pi, cos, sin2443 2358 2444 2359 a = [0.0, 0.0] … … 2450 2365 2451 2366 points = [a, b, c, d, e, f] 2452 # bac, bce, ecf,dbe2453 vertices = [ 2367 # bac, bce, ecf, dbe 2368 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2454 2369 2455 2370 domain = Domain(points, vertices) 2456 2371 2457 # Flat surface with 1m of water2372 # Flat surface with 1m of water 2458 2373 domain.set_quantity('elevation', 0) 2459 2374 domain.set_quantity('stage', 1.0) … … 2463 2378 domain.set_boundary({'exterior': Br}) 2464 2379 2465 2466 domain.time = 5.54 #Take a random time (not zero) 2467 2468 #Setup only one forcing term, bad func 2380 domain.time = 5.54 # Take a random time (not zero) 2381 2382 # Setup only one forcing term, bad func 2469 2383 domain.forcing_terms = [] 2470 2384 2471 2385 try: 2472 domain.forcing_terms.append(Wind_stress(s =scalar_func_list,2473 phi =angle))2386 domain.forcing_terms.append(Wind_stress(s=scalar_func_list, 2387 phi=angle)) 2474 2388 except AssertionError: 2475 2389 pass 2476 2390 else: 2477 2391 msg = 'Should have raised exception' 2478 raise msg 2479 2392 raise Exception, msg 2480 2393 2481 2394 try: 2482 domain.forcing_terms.append(Wind_stress(s = speed, 2483 phi = scalar_func)) 2395 domain.forcing_terms.append(Wind_stress(s=speed, phi=scalar_func)) 2484 2396 except Exception: 2485 2397 pass 2486 2398 else: 2487 2399 msg = 'Should have raised exception' 2488 raise msg2400 raise Exception, msg 2489 2401 2490 2402 try: 2491 domain.forcing_terms.append(Wind_stress(s = speed, 2492 phi = 'xx')) 2403 domain.forcing_terms.append(Wind_stress(s=speed, phi='xx')) 2493 2404 except: 2494 2405 pass 2495 2406 else: 2496 2407 msg = 'Should have raised exception' 2497 raise msg 2498 2499 2408 raise Exception, msg 2500 2409 2501 2410 def test_rainfall(self): … … 2510 2419 2511 2420 points = [a, b, c, d, e, f] 2512 #bac, bce, ecf, dbe 2513 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2514 2421 # bac, bce, ecf, dbe 2422 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2515 2423 2516 2424 domain = Domain(points, vertices) 2517 2425 2518 # Flat surface with 1m of water2426 # Flat surface with 1m of water 2519 2427 domain.set_quantity('elevation', 0) 2520 2428 domain.set_quantity('stage', 1.0) … … 2526 2434 # Setup only one forcing term, constant rainfall 2527 2435 domain.forcing_terms = [] 2528 domain.forcing_terms.append( Rainfall(domain, rate=2.0))2436 domain.forcing_terms.append(Rainfall(domain, rate=2.0)) 2529 2437 2530 2438 domain.compute_forcing_terms() 2531 assert num.allclose(domain.quantities['stage'].explicit_update, 2.0/1000) 2532 2533 2439 assert num.allclose(domain.quantities['stage'].explicit_update, 2440 2.0/1000) 2534 2441 2535 2442 def test_rainfall_restricted_by_polygon(self): … … 2544 2451 2545 2452 points = [a, b, c, d, e, f] 2546 #bac, bce, ecf, dbe 2547 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2548 2453 # bac, bce, ecf, dbe 2454 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2549 2455 2550 2456 domain = Domain(points, vertices) 2551 2457 2552 # Flat surface with 1m of water2458 # Flat surface with 1m of water 2553 2459 domain.set_quantity('elevation', 0) 2554 2460 domain.set_quantity('stage', 1.0) … … 2558 2464 domain.set_boundary({'exterior': Br}) 2559 2465 2560 # Setup only one forcing term, constant rainfall restricted to a polygon enclosing triangle #1 (bce) 2466 # Setup only one forcing term, constant rainfall 2467 # restricted to a polygon enclosing triangle #1 (bce) 2561 2468 domain.forcing_terms = [] 2562 R = Rainfall(domain, 2563 rate=2.0, 2564 polygon = [[1,1], [2,1], [2,2], [1,2]]) 2469 R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]]) 2565 2470 2566 2471 assert num.allclose(R.exchange_area, 1) 2567 2472 2568 2473 domain.forcing_terms.append(R) 2569 2474 2570 2475 domain.compute_forcing_terms() 2571 #print domain.quantities['stage'].explicit_update 2572 2476 2573 2477 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2574 2478 2.0/1000) 2575 2479 assert num.allclose(domain.quantities['stage'].explicit_update[0], 0) 2576 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2577 2578 2480 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2579 2481 2580 2482 def test_time_dependent_rainfall_restricted_by_polygon(self): 2581 2582 2483 a = [0.0, 0.0] 2583 2484 b = [0.0, 2.0] … … 2588 2489 2589 2490 points = [a, b, c, d, e, f] 2590 #bac, bce, ecf, dbe 2591 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2592 2491 # bac, bce, ecf, dbe 2492 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2593 2493 2594 2494 domain = Domain(points, vertices) 2595 2495 2596 # Flat surface with 1m of water2496 # Flat surface with 1m of water 2597 2497 domain.set_quantity('elevation', 0) 2598 2498 domain.set_quantity('stage', 1.0) … … 2602 2502 domain.set_boundary({'exterior': Br}) 2603 2503 2604 # Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce) 2504 # Setup only one forcing term, time dependent rainfall 2505 # restricted to a polygon enclosing triangle #1 (bce) 2605 2506 domain.forcing_terms = [] 2606 R = Rainfall(domain, 2607 rate=lambda t: 3*t + 7, 2608 polygon = [[1,1], [2,1], [2,2], [1,2]]) 2507 R = Rainfall(domain, rate=lambda t: 3*t + 7, 2508 polygon=[[1,1], [2,1], [2,2], [1,2]]) 2609 2509 2610 2510 assert num.allclose(R.exchange_area, 1) 2611 2511 2612 2512 domain.forcing_terms.append(R) 2613 2513 2614 2615 2514 domain.time = 10. 2616 2515 2617 2516 domain.compute_forcing_terms() 2618 #print domain.quantities['stage'].explicit_update 2619 2517 2620 2518 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2621 (3*domain.time +7)/1000)2519 (3*domain.time + 7)/1000) 2622 2520 assert num.allclose(domain.quantities['stage'].explicit_update[0], 0) 2623 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2624 2625 2626 2521 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2627 2522 2628 2523 def test_time_dependent_rainfall_using_starttime(self): 2629 2630 rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float) 2524 rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float) 2631 2525 2632 2526 a = [0.0, 0.0] … … 2638 2532 2639 2533 points = [a, b, c, d, e, f] 2640 #bac, bce, ecf, dbe 2641 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2642 2534 # bac, bce, ecf, dbe 2535 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2643 2536 2644 2537 domain = Domain(points, vertices) 2645 2538 2646 # Flat surface with 1m of water2539 # Flat surface with 1m of water 2647 2540 domain.set_quantity('elevation', 0) 2648 2541 domain.set_quantity('stage', 1.0) … … 2652 2545 domain.set_boundary({'exterior': Br}) 2653 2546 2654 # Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce) 2547 # Setup only one forcing term, time dependent rainfall 2548 # restricted to a polygon enclosing triangle #1 (bce) 2655 2549 domain.forcing_terms = [] 2656 R = Rainfall(domain, 2657 rate=lambda t: 3*t + 7, 2658 polygon=rainfall_poly) 2550 R = Rainfall(domain, rate=lambda t: 3*t + 7, 2551 polygon=rainfall_poly) 2659 2552 2660 2553 assert num.allclose(R.exchange_area, 1) 2661 2554 2662 2555 domain.forcing_terms.append(R) 2663 2556 … … 2669 2562 2670 2563 domain.compute_forcing_terms() 2671 #print domain.quantities['stage'].explicit_update 2672 2673 #print domain.get_time() 2564 2674 2565 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2675 (3*domain.get_time() +7)/1000)2566 (3*domain.get_time() + 7)/1000) 2676 2567 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2677 (3*(domain.time + domain.starttime) +7)/1000)2568 (3*(domain.time + domain.starttime) + 7)/1000) 2678 2569 2679 2570 # Using internal time her should fail 2680 2571 assert not num.allclose(domain.quantities['stage'].explicit_update[1], 2681 (3*domain.time +7)/1000)2572 (3*domain.time + 7)/1000) 2682 2573 2683 2574 assert num.allclose(domain.quantities['stage'].explicit_update[0], 0) 2684 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2685 2686 2687 2688 2575 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2576 2689 2577 def test_time_dependent_rainfall_using_georef(self): 2690 2578 """test_time_dependent_rainfall_using_georef 2691 2579 2692 2580 This will also test the General forcing term using georef 2693 2581 """ 2694 2695 #Mesh in zone 56 (absolute coords) 2696 2582 2583 # Mesh in zone 56 (absolute coords) 2697 2584 x0 = 314036.58727982 2698 2585 y0 = 6224951.2960092 2699 2586 2700 2701 2587 rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float) 2702 2588 rainfall_poly += [x0, y0] … … 2710 2596 2711 2597 points = [a, b, c, d, e, f] 2712 #bac, bce, ecf, dbe 2713 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2714 2598 # bac, bce, ecf, dbe 2599 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2715 2600 2716 2601 domain = Domain(points, vertices, 2717 geo_reference =Geo_reference(56, x0, y0))2718 2719 # Flat surface with 1m of water2602 geo_reference=Geo_reference(56, x0, y0)) 2603 2604 # Flat surface with 1m of water 2720 2605 domain.set_quantity('elevation', 0) 2721 2606 domain.set_quantity('stage', 1.0) … … 2725 2610 domain.set_boundary({'exterior': Br}) 2726 2611 2727 # Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce) 2612 # Setup only one forcing term, time dependent rainfall 2613 # restricted to a polygon enclosing triangle #1 (bce) 2728 2614 domain.forcing_terms = [] 2729 R = Rainfall(domain, 2730 rate=lambda t: 3*t + 7, 2731 polygon=rainfall_poly) 2615 R = Rainfall(domain, rate=lambda t: 3*t + 7, polygon=rainfall_poly) 2732 2616 2733 2617 assert num.allclose(R.exchange_area, 1) 2734 2618 2735 2619 domain.forcing_terms.append(R) 2736 2620 … … 2742 2626 2743 2627 domain.compute_forcing_terms() 2744 #print domain.quantities['stage'].explicit_update 2745 2746 #print domain.get_time() 2628 2747 2629 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2748 (3*domain.get_time() +7)/1000)2630 (3*domain.get_time() + 7)/1000) 2749 2631 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2750 (3*(domain.time + domain.starttime) +7)/1000)2632 (3*(domain.time + domain.starttime) + 7)/1000) 2751 2633 2752 2634 # Using internal time her should fail 2753 2635 assert not num.allclose(domain.quantities['stage'].explicit_update[1], 2754 (3*domain.time +7)/1000)2636 (3*domain.time + 7)/1000) 2755 2637 2756 2638 assert num.allclose(domain.quantities['stage'].explicit_update[0], 0) 2757 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2758 2759 2760 2761 2762 2639 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2763 2640 2764 2641 def test_time_dependent_rainfall_restricted_by_polygon_with_default(self): 2765 """test_time_dependent_rainfall_restricted_by_polygon_with_default 2766 2642 """ 2767 2643 Test that default rainfall can be used when given rate runs out of data. 2768 2644 """ 2645 2769 2646 a = [0.0, 0.0] 2770 2647 b = [0.0, 2.0] … … 2775 2652 2776 2653 points = [a, b, c, d, e, f] 2777 #bac, bce, ecf, dbe 2778 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2779 2654 # bac, bce, ecf, dbe 2655 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2780 2656 2781 2657 domain = Domain(points, vertices) 2782 2658 2783 # Flat surface with 1m of water2659 # Flat surface with 1m of water 2784 2660 domain.set_quantity('elevation', 0) 2785 2661 domain.set_quantity('stage', 1.0) … … 2789 2665 domain.set_boundary({'exterior': Br}) 2790 2666 2791 # Setup only one forcing term, time dependent rainfall that expires at t==20 2667 # Setup only one forcing term, time dependent rainfall 2668 # that expires at t==20 2792 2669 from anuga.fit_interpolate.interpolate import Modeltime_too_late 2670 2793 2671 def main_rate(t): 2794 2672 if t > 20: … … 2797 2675 else: 2798 2676 return 3*t + 7 2799 2677 2800 2678 domain.forcing_terms = [] 2801 R = Rainfall(domain, 2802 rate=main_rate, 2803 polygon = [[1,1], [2,1], [2,2], [1,2]], 2804 default_rate=5.0) 2679 R = Rainfall(domain, rate=main_rate, 2680 polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0) 2805 2681 2806 2682 assert num.allclose(R.exchange_area, 1) 2807 2683 2808 2684 domain.forcing_terms.append(R) 2809 2685 2810 2811 2686 domain.time = 10. 2812 2687 2813 2688 domain.compute_forcing_terms() 2814 #print domain.quantities['stage'].explicit_update 2815 2816 assert num.allclose(domain.quantities['stage'].explicit_update[1],(3*domain.time+7)/1000)2689 2690 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2691 (3*domain.time+7)/1000) 2817 2692 assert num.allclose(domain.quantities['stage'].explicit_update[0], 0) 2818 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2819 2693 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2820 2694 2821 2695 domain.time = 100. 2822 domain.quantities['stage'].explicit_update[:] = 0.0 # Reset2696 domain.quantities['stage'].explicit_update[:] = 0.0 # Reset 2823 2697 domain.compute_forcing_terms() 2824 #print domain.quantities['stage'].explicit_update 2825 2826 assert num.allclose(domain.quantities['stage'].explicit_update[1],5.0/1000) # Default value2698 2699 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2700 5.0/1000) # Default value 2827 2701 assert num.allclose(domain.quantities['stage'].explicit_update[0], 0) 2828 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2829 2830 2831 2832 2833 2834 2702 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2835 2703 2836 2704 def test_rainfall_forcing_with_evolve(self): … … 2839 2707 Test how forcing terms are called within evolve 2840 2708 """ 2841 2709 2842 2710 # FIXME(Ole): This test is just to experiment 2843 2711 2844 2712 a = [0.0, 0.0] 2845 2713 b = [0.0, 2.0] … … 2850 2718 2851 2719 points = [a, b, c, d, e, f] 2852 #bac, bce, ecf, dbe 2853 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2854 2720 # bac, bce, ecf, dbe 2721 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2855 2722 2856 2723 domain = Domain(points, vertices) 2857 2724 2858 # Flat surface with 1m of water2725 # Flat surface with 1m of water 2859 2726 domain.set_quantity('elevation', 0) 2860 2727 domain.set_quantity('stage', 1.0) … … 2864 2731 domain.set_boundary({'exterior': Br}) 2865 2732 2866 # Setup only one forcing term, time dependent rainfall that expires at t==20 2733 # Setup only one forcing term, time dependent rainfall 2734 # that expires at t==20 2867 2735 from anuga.fit_interpolate.interpolate import Modeltime_too_late 2736 2868 2737 def main_rate(t): 2869 2738 if t > 20: … … 2872 2741 else: 2873 2742 return 3*t + 7 2874 2743 2875 2744 domain.forcing_terms = [] 2876 R = Rainfall(domain, 2877 rate=main_rate, 2878 polygon=[[1,1], [2,1], [2,2], [1,2]], 2879 default_rate=5.0) 2745 R = Rainfall(domain, rate=main_rate, 2746 polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0) 2880 2747 2881 2748 assert num.allclose(R.exchange_area, 1) 2882 2749 2883 2750 domain.forcing_terms.append(R) 2884 2751 2885 2752 for t in domain.evolve(yieldstep=1, finaltime=25): 2886 2753 pass 2887 2888 #print t, domain.quantities['stage'].explicit_update, (3*t+7)/10002889 2890 2754 #FIXME(Ole): A test here is hard because explicit_update also 2891 2755 # receives updates from the flux calculation. 2892 2893 2894 2895 2756 2896 2757 def test_inflow_using_circle(self): … … 2905 2766 2906 2767 points = [a, b, c, d, e, f] 2907 #bac, bce, ecf, dbe 2908 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2909 2768 # bac, bce, ecf, dbe 2769 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2910 2770 2911 2771 domain = Domain(points, vertices) … … 2919 2779 domain.set_boundary({'exterior': Br}) 2920 2780 2921 # Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce) 2781 # Setup only one forcing term, constant inflow of 2 m^3/s 2782 # on a circle affecting triangles #0 and #1 (bac and bce) 2922 2783 domain.forcing_terms = [] 2923 domain.forcing_terms.append( Inflow(domain, rate=2.0, center=(1,1), radius=1) ) 2784 domain.forcing_terms.append(Inflow(domain, rate=2.0, 2785 center=(1,1), radius=1)) 2924 2786 2925 2787 domain.compute_forcing_terms() 2926 #print domain.quantities['stage'].explicit_update 2927 2928 assert num.allclose(domain.quantities['stage'].explicit_update[1],2.0/pi)2929 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)2930 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)2931 2788 2789 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2790 2.0/pi) 2791 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2792 2.0/pi) 2793 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2932 2794 2933 2795 def test_inflow_using_circle_function(self): … … 2942 2804 2943 2805 points = [a, b, c, d, e, f] 2944 #bac, bce, ecf, dbe 2945 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2946 2806 # bac, bce, ecf, dbe 2807 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2947 2808 2948 2809 domain = Domain(points, vertices) … … 2956 2817 domain.set_boundary({'exterior': Br}) 2957 2818 2958 # Setup only one forcing term, time dependent inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce) 2819 # Setup only one forcing term, time dependent inflow of 2 m^3/s 2820 # on a circle affecting triangles #0 and #1 (bac and bce) 2959 2821 domain.forcing_terms = [] 2960 domain.forcing_terms.append( Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) ) 2822 domain.forcing_terms.append(Inflow(domain, rate=lambda t: 2., 2823 center=(1,1), radius=1)) 2961 2824 2962 2825 domain.compute_forcing_terms() 2963 2964 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi) 2965 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi) 2966 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2967 2968 2969 2826 2827 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2828 2.0/pi) 2829 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2830 2.0/pi) 2831 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2970 2832 2971 2833 def test_inflow_catch_too_few_triangles(self): 2972 """test_inflow_catch_too_few_triangles2973 2974 Test that exception is thrown if no triangles are covered by the inflow area2975 2834 """ 2835 Test that exception is thrown if no triangles are covered 2836 by the inflow area 2837 """ 2838 2976 2839 from math import pi, cos, sin 2977 2840 … … 2984 2847 2985 2848 points = [a, b, c, d, e, f] 2986 #bac, bce, ecf, dbe 2987 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2988 2849 # bac, bce, ecf, dbe 2850 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2989 2851 2990 2852 domain = Domain(points, vertices) … … 2998 2860 domain.set_boundary({'exterior': Br}) 2999 2861 3000 # Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)3001 2862 # Setup only one forcing term, constant inflow of 2 m^3/s 2863 # on a circle affecting triangles #0 and #1 (bac and bce) 3002 2864 try: 3003 2865 Inflow(domain, rate=2.0, center=(1,1.1), radius=0.01) … … 3008 2870 raise Exception, msg 3009 2871 3010 3011 3012 3013 2872 def Xtest_inflow_outflow_conservation(self): 3014 """test_inflow_outflow_conservation 3015 3016 Test what happens if water is abstracted from one area and 2873 """ 2874 Test what happens if water is abstracted from one area and 3017 2875 injected into another - especially if there is not enough 3018 water to match the abstraction. 2876 water to match the abstraction. 3019 2877 This tests that the total volume is kept constant under a range of 3020 2878 scenarios. 3021 2879 3022 2880 This test will fail as the problem was only fixed for culverts. 3023 2881 """ 3024 2882 3025 2883 from math import pi, cos, sin 3026 2884 3027 2885 length = 20. 3028 2886 width = 10. 3029 2887 3030 dx = dy = 2 # 1 or 2 OK2888 dx = dy = 2 # 1 or 2 OK 3031 2889 points, vertices, boundary = rectangular_cross(int(length/dx), 3032 2890 int(width/dy), 3033 len1=length, 2891 len1=length, 3034 2892 len2=width) 3035 domain = Domain(points, vertices, boundary) 3036 domain.set_name('test_inflow_conservation') # Output name2893 domain = Domain(points, vertices, boundary) 2894 domain.set_name('test_inflow_conservation') # Output name 3037 2895 domain.set_default_order(2) 3038 3039 2896 3040 2897 # Flat surface with 1m of water … … 3047 2904 domain.set_boundary({'left': Br, 'right': Br, 'bottom': Br, 'top': Br}) 3048 2905 3049 # Setup one forcing term, constant inflow of 2 m^3/s on a circle 2906 # Setup one forcing term, constant inflow of 2 m^3/s on a circle 3050 2907 domain.forcing_terms = [] 3051 domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1)) 2908 domain.forcing_terms.append(Inflow(domain, rate=2.0, 2909 center=(5,5), radius=1)) 3052 2910 3053 2911 domain.compute_forcing_terms() 3054 #print domain.quantities['stage'].explicit_update 3055 2912 3056 2913 # Check that update values are correct 3057 2914 for x in domain.quantities['stage'].explicit_update: 3058 2915 assert num.allclose(x, 2.0/pi) or num.allclose(x, 0.0) 3059 2916 3060 3061 2917 # Check volumes without inflow 3062 domain.forcing_terms = [] 2918 domain.forcing_terms = [] 3063 2919 initial_volume = domain.quantities['stage'].get_integral() 3064 2920 3065 2921 assert num.allclose(initial_volume, width*length*stage) 3066 2922 3067 2923 for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0): 3068 volume = domain.quantities['stage'].get_integral() 3069 assert num.allclose (volume, initial_volume) 3070 3071 2924 volume = domain.quantities['stage'].get_integral() 2925 assert num.allclose(volume, initial_volume) 2926 3072 2927 # Now apply the inflow and check volumes for a range of stage values 3073 2928 for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]: 3074 2929 domain.time = 0.0 3075 domain.set_quantity('stage', stage) 3076 3077 domain.forcing_terms = []3078 domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))2930 domain.set_quantity('stage', stage) 2931 domain.forcing_terms = [] 2932 domain.forcing_terms.append(Inflow(domain, rate=2.0, 2933 center=(5,5), radius=1)) 3079 2934 initial_volume = domain.quantities['stage'].get_integral() 3080 2935 predicted_volume = initial_volume 3081 2936 dt = 0.05 3082 for t in domain.evolve(yieldstep = dt, finaltime =5.0):2937 for t in domain.evolve(yieldstep=dt, finaltime=5.0): 3083 2938 volume = domain.quantities['stage'].get_integral() 3084 3085 assert num.allclose (volume, predicted_volume) 2939 assert num.allclose (volume, predicted_volume) 3086 2940 predicted_volume = predicted_volume + 2.0/pi/100/dt # Why 100? 3087 3088 3089 # Apply equivalent outflow only and check volumesfor a range of stage values2941 2942 # Apply equivalent outflow only and check volumes 2943 # for a range of stage values 3090 2944 for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]: 3091 2945 print stage 3092 2946 3093 2947 domain.time = 0.0 3094 domain.set_quantity('stage', stage) 3095 domain.forcing_terms = [] 3096 domain.forcing_terms.append(Inflow(domain, rate=-2.0, center=(15,5), radius=1)) 2948 domain.set_quantity('stage', stage) 2949 domain.forcing_terms = [] 2950 domain.forcing_terms.append(Inflow(domain, rate=-2.0, 2951 center=(15,5), radius=1)) 3097 2952 initial_volume = domain.quantities['stage'].get_integral() 3098 2953 predicted_volume = initial_volume 3099 2954 dt = 0.05 3100 for t in domain.evolve(yieldstep = dt, finaltime =5.0):2955 for t in domain.evolve(yieldstep=dt, finaltime=5.0): 3101 2956 volume = domain.quantities['stage'].get_integral() 3102 3103 2957 print t, volume, predicted_volume 3104 assert num.allclose (volume, predicted_volume) 3105 predicted_volume = predicted_volume - 2.0/pi/100/dt # Why 100? 3106 3107 3108 # Apply both inflow and outflow and check volumes being constant for a 3109 # range of stage values 3110 for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]: 2958 assert num.allclose (volume, predicted_volume) 2959 predicted_volume = predicted_volume - 2.0/pi/100/dt # Why 100? 2960 2961 # Apply both inflow and outflow and check volumes being constant for a 2962 # range of stage values 2963 for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]: 3111 2964 print stage 3112 2965 3113 2966 domain.time = 0.0 3114 domain.set_quantity('stage', stage) 3115 domain.forcing_terms = [] 3116 domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1)) 3117 domain.forcing_terms.append(Inflow(domain, rate=-2.0, center=(15,5), radius=1)) 2967 domain.set_quantity('stage', stage) 2968 domain.forcing_terms = [] 2969 domain.forcing_terms.append(Inflow(domain, rate=2.0, 2970 center=(5,5), radius=1)) 2971 domain.forcing_terms.append(Inflow(domain, rate=-2.0, 2972 center=(15,5), radius=1)) 3118 2973 initial_volume = domain.quantities['stage'].get_integral() 3119 2974 3120 2975 dt = 0.05 3121 for t in domain.evolve(yieldstep = dt, finaltime =5.0):2976 for t in domain.evolve(yieldstep=dt, finaltime=5.0): 3122 2977 volume = domain.quantities['stage'].get_integral() 3123 2978 3124 2979 print t, volume 3125 assert num.allclose (volume, initial_volume) 3126 3127 3128 2980 assert num.allclose(volume, initial_volume) 3129 2981 3130 2982 ##################################################### 2983 3131 2984 def test_first_order_extrapolator_const_z(self): 3132 3133 2985 a = [0.0, 0.0] 3134 2986 b = [0.0, 2.0] … … 3139 2991 3140 2992 points = [a, b, c, d, e, f] 3141 # bac, bce, ecf,dbe3142 vertices = [ 2993 # bac, bce, ecf, dbe 2994 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 3143 2995 3144 2996 domain = Domain(points, vertices) 3145 val0 = 2. +2.0/33146 val1 = 4. +4.0/33147 val2 = 8. +2.0/33148 val3 = 2. +8.0/33149 3150 zl =zr=-3.75 #Assume constant bed (must be less than stage)3151 domain.set_quantity('elevation', zl*num.ones( (4,3), num.int)) #array default#2997 val0 = 2. + 2.0/3 2998 val1 = 4. + 4.0/3 2999 val2 = 8. + 2.0/3 3000 val3 = 2. + 8.0/3 3001 3002 zl = zr = -3.75 # Assume constant bed (must be less than stage) 3003 domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default# 3152 3004 domain.set_quantity('stage', [[val0, val0-1, val0-2], 3153 3005 [val1, val1+1, val1], … … 3155 3007 [val3-0.5, val3, val3]]) 3156 3008 3157 3158 3159 3009 domain._order_ = 1 3160 3010 domain.distribute_to_vertices_and_edges() … … 3163 3013 C = domain.quantities['stage'].centroid_values 3164 3014 for i in range(3): 3165 assert num.allclose( domain.quantities['stage'].vertex_values[:,i], C)3166 3015 assert num.allclose(domain.quantities['stage'].vertex_values[:,i], 3016 C) 3167 3017 3168 3018 def test_first_order_limiter_variable_z(self): 3169 #Check that first order limiter follows bed_slope 3019 '''Check that first order limiter follows bed_slope''' 3020 3170 3021 from anuga.config import epsilon 3171 3022 3172 3023 a = [0.0, 0.0] 3173 3024 b = [0.0, 2.0] 3174 c = [2.0, 0.0]3025 c = [2.0, 0.0] 3175 3026 d = [0.0, 4.0] 3176 3027 e = [2.0, 2.0] 3177 f = [4.0, 0.0]3028 f = [4.0, 0.0] 3178 3029 3179 3030 points = [a, b, c, d, e, f] 3180 # bac, bce, ecf,dbe3181 vertices = [ 3031 # bac, bce, ecf, dbe 3032 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 3182 3033 3183 3034 domain = Domain(points, vertices) 3184 val0 = 2. +2.0/33185 val1 = 4. +4.0/33186 val2 = 8. +2.0/33187 val3 = 2. +8.0/33035 val0 = 2. + 2.0/3 3036 val1 = 4. + 4.0/3 3037 val2 = 8. + 2.0/3 3038 val3 = 2. + 8.0/3 3188 3039 3189 3040 domain.set_quantity('elevation', [[0,0,0], [6,0,0], … … 3197 3048 L = domain.quantities['stage'].vertex_values 3198 3049 3199 3200 #Check that some stages are not above elevation (within eps) 3201 #- so that the limiter has something to work with 3202 assert not num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon))) 3050 # Check that some stages are not above elevation (within eps) - 3051 # so that the limiter has something to work with 3052 assert not num.alltrue(num.alltrue(num.greater_equal(L, E-epsilon))) 3203 3053 3204 3054 domain._order_ = 1 … … 3206 3056 3207 3057 #Check that all stages are above elevation (within eps) 3208 assert num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon))) 3209 3058 assert num.alltrue(num.alltrue(num.greater_equal(L, E-epsilon))) 3210 3059 3211 3060 ##################################################### 3061 3212 3062 def test_distribute_basic(self): 3213 3063 #Using test data generated by abstract_2d_finite_volumes-2 … … 3222 3072 3223 3073 points = [a, b, c, d, e, f] 3224 # bac, bce, ecf,dbe3225 vertices = [ 3074 # bac, bce, ecf, dbe 3075 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 3226 3076 3227 3077 domain = Domain(points, vertices) … … 3236 3086 L = domain.quantities['stage'].vertex_values 3237 3087 3238 # First order3088 # First order 3239 3089 domain._order_ = 1 3240 3090 domain.distribute_to_vertices_and_edges() 3241 3091 assert num.allclose(L[1], val1) 3242 3092 3243 # Second order3093 # Second order 3244 3094 domain._order_ = 2 3245 domain.beta_w 3246 domain.beta_w_dry 3247 domain.beta_uh 3095 domain.beta_w = 0.9 3096 domain.beta_w_dry = 0.9 3097 domain.beta_uh = 0.9 3248 3098 domain.beta_uh_dry = 0.9 3249 domain.beta_vh 3099 domain.beta_vh = 0.9 3250 3100 domain.beta_vh_dry = 0.9 3251 3101 domain.distribute_to_vertices_and_edges() 3252 3102 assert num.allclose(L[1], [2.2, 4.9, 4.9]) 3253 3254 3255 3103 3256 3104 def test_distribute_away_from_bed(self): … … 3266 3114 3267 3115 points = [a, b, c, d, e, f] 3268 # bac, bce, ecf,dbe3269 vertices = [ 3116 # bac, bce, ecf, dbe 3117 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 3270 3118 3271 3119 domain = Domain(points, vertices) 3272 3120 L = domain.quantities['stage'].vertex_values 3273 3121 3274 def stage(x, y):3122 def stage(x, y): 3275 3123 return x**2 3276 3124 … … 3280 3128 3281 3129 a, b = domain.quantities['stage'].get_gradients() 3282 3130 3283 3131 assert num.allclose(a[1], 3.33333334) 3284 3132 assert num.allclose(b[1], 0.0) … … 3289 3137 3290 3138 domain._order_ = 2 3291 domain.beta_w 3292 domain.beta_w_dry 3293 domain.beta_uh 3139 domain.beta_w = 0.9 3140 domain.beta_w_dry = 0.9 3141 domain.beta_uh = 0.9 3294 3142 domain.beta_uh_dry = 0.9 3295 domain.beta_vh 3143 domain.beta_vh = 0.9 3296 3144 domain.beta_vh_dry = 0.9 3297 3145 domain.distribute_to_vertices_and_edges() 3298 3146 assert num.allclose(L[1], [0.57777777, 2.37777778, 2.37777778]) 3299 3300 3301 3147 3302 3148 def test_distribute_away_from_bed1(self): … … 3312 3158 3313 3159 points = [a, b, c, d, e, f] 3314 # bac, bce, ecf,dbe3315 vertices = [ 3160 # bac, bce, ecf, dbe 3161 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 3316 3162 3317 3163 domain = Domain(points, vertices) 3318 3164 L = domain.quantities['stage'].vertex_values 3319 3165 3320 def stage(x, y):3321 return x**4 +y**23166 def stage(x, y): 3167 return x**4 + y**2 3322 3168 3323 3169 domain.set_quantity('stage', stage, location='centroids') 3324 #print domain.quantities['stage'].centroid_values3325 3170 3326 3171 domain.quantities['stage'].compute_gradients() … … 3334 3179 3335 3180 domain._order_ = 2 3336 domain.beta_w 3337 domain.beta_w_dry 3338 domain.beta_uh 3181 domain.beta_w = 0.9 3182 domain.beta_w_dry = 0.9 3183 domain.beta_uh = 0.9 3339 3184 domain.beta_uh_dry = 0.9 3340 domain.beta_vh 3185 domain.beta_vh = 0.9 3341 3186 domain.beta_vh_dry = 0.9 3342 3187 domain.distribute_to_vertices_and_edges() 3343 3188 assert num.allclose(L[1], [1.07160494, 6.46058131, 7.28262855]) 3344 3189 3345 3346 3347 3190 def test_distribute_near_bed(self): 3348 3349 3191 a = [0.0, 0.0] 3350 3192 b = [0.0, 2.0] … … 3355 3197 3356 3198 points = [a, b, c, d, e, f] 3357 # bac, bce, ecf,dbe3358 vertices = [ 3199 # bac, bce, ecf, dbe 3200 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 3359 3201 3360 3202 domain = Domain(points, vertices) 3361 3203 3362 3363 #Set up for a gradient of (10,0) at mid triangle (bce) 3204 # Set up for a gradient of (10,0) at mid triangle (bce) 3364 3205 def slope(x, y): 3365 3206 return 10*x … … 3371 3212 domain.set_quantity('elevation', slope) 3372 3213 domain.set_quantity('stage', stage, location='centroids') 3373 3374 #print domain.quantities['elevation'].centroid_values3375 #print domain.quantities['stage'].centroid_values3376 3214 3377 3215 E = domain.quantities['elevation'].vertex_values … … 3382 3220 for i in range(len(L)): 3383 3221 volumes.append(num.sum(L[i])/3) 3384 assert num.allclose(volumes[i], domain.quantities['stage'].centroid_values[i])3385 3386 3222 assert num.allclose(volumes[i], 3223 domain.quantities['stage'].centroid_values[i]) 3224 3387 3225 domain._order_ = 1 3388 3226 3389 3227 domain.tight_slope_limiters = 0 3390 3228 domain.distribute_to_vertices_and_edges() 3391 3229 assert num.allclose(L[1], [0.1, 20.1, 20.1]) 3392 3230 for i in range(len(L)): 3393 assert num.allclose(volumes[i], num.sum(L[i])/3) 3394 3395 domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed) 3231 assert num.allclose(volumes[i], num.sum(L[i])/3) 3232 3233 # Allow triangle to be flatter (closer to bed) 3234 domain.tight_slope_limiters = 1 3235 3396 3236 domain.distribute_to_vertices_and_edges() 3397 3237 assert num.allclose(L[1], [0.298, 20.001, 20.001]) 3398 3238 for i in range(len(L)): 3399 assert num.allclose(volumes[i], num.sum(L[i])/3) 3239 assert num.allclose(volumes[i], num.sum(L[i])/3) 3400 3240 3401 3241 domain._order_ = 2 3402 3242 3403 3243 domain.tight_slope_limiters = 0 3404 3244 domain.distribute_to_vertices_and_edges() 3405 assert num.allclose(L[1], [0.1, 20.1, 20.1]) 3245 assert num.allclose(L[1], [0.1, 20.1, 20.1]) 3406 3246 for i in range(len(L)): 3407 assert num.allclose(volumes[i], num.sum(L[i])/3) 3408 3409 domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed) 3247 assert num.allclose(volumes[i], num.sum(L[i])/3) 3248 3249 # Allow triangle to be flatter (closer to bed) 3250 domain.tight_slope_limiters = 1 3251 3410 3252 domain.distribute_to_vertices_and_edges() 3411 3253 assert num.allclose(L[1], [0.298, 20.001, 20.001]) 3412 3254 for i in range(len(L)): 3413 assert num.allclose(volumes[i], num.sum(L[i])/3) 3414 3415 3255 assert num.allclose(volumes[i], num.sum(L[i])/3) 3416 3256 3417 3257 def test_distribute_near_bed1(self): 3418 3419 3258 a = [0.0, 0.0] 3420 3259 b = [0.0, 2.0] … … 3425 3264 3426 3265 points = [a, b, c, d, e, f] 3427 # bac, bce, ecf,dbe3428 vertices = [ 3266 # bac, bce, ecf, dbe 3267 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 3429 3268 3430 3269 domain = Domain(points, vertices) 3431 3270 3432 3433 #Set up for a gradient of (8,2) at mid triangle (bce) 3271 # Set up for a gradient of (8,2) at mid triangle (bce) 3434 3272 def slope(x, y): 3435 return x**4 +y**23273 return x**4 + y**2 3436 3274 3437 3275 h = 0.1 3438 def stage(x, y):3439 return slope(x, y)+h3276 def stage(x, y): 3277 return slope(x, y) + h 3440 3278 3441 3279 domain.set_quantity('elevation', slope) 3442 3280 domain.set_quantity('stage', stage) 3443 3444 #print domain.quantities['elevation'].centroid_values3445 #print domain.quantities['stage'].centroid_values3446 3281 3447 3282 E = domain.quantities['elevation'].vertex_values … … 3452 3287 for i in range(len(L)): 3453 3288 volumes.append(num.sum(L[i])/3) 3454 assert num.allclose(volumes[i], domain.quantities['stage'].centroid_values[i])3455 3456 #print E 3289 assert num.allclose(volumes[i], 3290 domain.quantities['stage'].centroid_values[i]) 3291 3457 3292 domain._order_ = 1 3458 3293 3459 3294 domain.tight_slope_limiters = 0 3460 3295 domain.distribute_to_vertices_and_edges() 3461 assert num.allclose(L[1], [4.1, 16.1, 20.1]) 3296 assert num.allclose(L[1], [4.1, 16.1, 20.1]) 3462 3297 for i in range(len(L)): 3463 3298 assert num.allclose(volumes[i], num.sum(L[i])/3) 3464 3465 3466 domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed) 3299 3300 # Allow triangle to be flatter (closer to bed) 3301 domain.tight_slope_limiters = 1 3302 3467 3303 domain.distribute_to_vertices_and_edges() 3468 3304 assert num.allclose(L[1], [4.2386, 16.0604, 20.001]) 3469 3305 for i in range(len(L)): 3470 assert num.allclose(volumes[i], num.sum(L[i])/3) 3471 3306 assert num.allclose(volumes[i], num.sum(L[i])/3) 3472 3307 3473 3308 domain._order_ = 2 3474 3475 domain.tight_slope_limiters = 0 3309 3310 domain.tight_slope_limiters = 0 3476 3311 domain.distribute_to_vertices_and_edges() 3477 3312 assert num.allclose(L[1], [4.1, 16.1, 20.1]) 3478 3313 for i in range(len(L)): 3479 assert num.allclose(volumes[i], num.sum(L[i])/3) 3480 3481 domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed) 3314 assert num.allclose(volumes[i], num.sum(L[i])/3) 3315 3316 # Allow triangle to be flatter (closer to bed) 3317 domain.tight_slope_limiters = 1 3318 3482 3319 domain.distribute_to_vertices_and_edges() 3483 #print L[1]3484 assert num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\3485 num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\3486 num.allclose(L[1], [4.19351461, 16.10548539, 20.001])# old limiters3487 3320 assert (num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or 3321 num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or 3322 num.allclose(L[1], [4.19351461, 16.10548539, 20.001])) 3323 # old limiters 3324 3488 3325 for i in range(len(L)): 3489 3326 assert num.allclose(volumes[i], num.sum(L[i])/3) 3490 3491 3327 3492 3328 def test_second_order_distribute_real_data(self): … … 3503 3339 3504 3340 points = [a, b, c, d, e, f, g] 3505 # bae, efb, cbf,feg3506 vertices = [ 3341 # bae, efb, cbf, feg 3342 vertices = [[1,0,4], [4,5,1], [2,1,5], [5,4,6]] 3507 3343 3508 3344 domain = Domain(points, vertices) … … 3530 3366 Y = domain.quantities['ymomentum'].vertex_values 3531 3367 3532 #print E3533 3368 domain._order_ = 2 3534 domain.beta_w 3535 domain.beta_w_dry 3536 domain.beta_uh 3369 domain.beta_w = 0.9 3370 domain.beta_w_dry = 0.9 3371 domain.beta_uh = 0.9 3537 3372 domain.beta_uh_dry = 0.9 3538 domain.beta_vh 3373 domain.beta_vh = 0.9 3539 3374 domain.beta_vh_dry = 0.9 3540 3375 3541 3376 # FIXME (Ole): Need tests where this is commented out 3542 domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7) 3543 domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8) 3544 3545 3377 domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7) 3378 domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8) 3379 3546 3380 domain.distribute_to_vertices_and_edges() 3547 3548 #print L[1,:]3549 #print X[1,:]3550 #print Y[1,:]3551 3381 3552 3382 assert num.allclose(L[1,:], [-0.00825735775384, … … 3560 3390 -0.000151505429018]) 3561 3391 3562 3563 3564 3392 def test_balance_deep_and_shallow(self): 3565 3393 """Test that balanced limiters preserve conserved quantites. 3566 3394 This test is using old depth based balanced limiters 3567 3395 """ 3396 3568 3397 import copy 3569 3398 … … 3576 3405 3577 3406 points = [a, b, c, d, e, f] 3578 3579 #bac, bce, ecf, dbe 3580 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 3407 # bac, bce, ecf, dbe 3408 elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 3581 3409 3582 3410 domain = Domain(points, elements) 3583 3411 domain.check_integrity() 3584 3412 3585 # Create a deliberate overshoot3413 # Create a deliberate overshoot 3586 3414 domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 3587 domain.set_quantity('elevation', 0) #Flat bed3415 domain.set_quantity('elevation', 0) # Flat bed 3588 3416 stage = domain.quantities['stage'] 3589 3417 3590 ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy3591 3592 # Limit3593 domain.tight_slope_limiters = 0 3418 ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy 3419 3420 # Limit 3421 domain.tight_slope_limiters = 0 3594 3422 domain.distribute_to_vertices_and_edges() 3595 3423 3596 # Assert that quantities are conserved3424 # Assert that quantities are conserved 3597 3425 for k in range(len(domain)): 3598 assert num.allclose (ref_centroid_values[k], 3599 num.sum(stage.vertex_values[k,:])/3) 3600 3601 3602 #Now try with a non-flat bed - closely hugging initial stage in places 3603 #This will create alphas in the range [0, 0.478260, 1] 3426 assert num.allclose(ref_centroid_values[k], 3427 num.sum(stage.vertex_values[k,:])/3) 3428 3429 # Now try with a non-flat bed - closely hugging initial stage in places 3430 # This will create alphas in the range [0, 0.478260, 1] 3604 3431 domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 3605 3432 domain.set_quantity('elevation', [[0,0,0], 3606 [1.8,1.9,5.9],3607 [4.6,0,0],3608 [0,2,4]])3433 [1.8,1.9,5.9], 3434 [4.6,0,0], 3435 [0,2,4]]) 3609 3436 stage = domain.quantities['stage'] 3610 3437 3611 ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy3612 ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy3613 3614 # Limit3615 domain.tight_slope_limiters = 0 3438 ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy 3439 ref_vertex_values = copy.copy(stage.vertex_values[:]) # Copy 3440 3441 # Limit 3442 domain.tight_slope_limiters = 0 3616 3443 domain.distribute_to_vertices_and_edges() 3617 3444 3618 3619 #Assert that all vertex quantities have changed 3445 # Assert that all vertex quantities have changed 3620 3446 for k in range(len(domain)): 3621 #print ref_vertex_values[k,:], stage.vertex_values[k,:]3622 assert not num.allclose (ref_vertex_values[k,:],stage.vertex_values[k,:])3623 # and assert that quantities are still conserved3447 assert not num.allclose(ref_vertex_values[k,:], 3448 stage.vertex_values[k,:]) 3449 # and assert that quantities are still conserved 3624 3450 for k in range(len(domain)): 3625 assert num.allclose (ref_centroid_values[k], 3626 num.sum(stage.vertex_values[k,:])/3) 3627 3451 assert num.allclose(ref_centroid_values[k], 3452 num.sum(stage.vertex_values[k,:])/3) 3628 3453 3629 3454 # Check actual results 3630 assert num.allclose (stage.vertex_values, 3631 [[2,2,2], 3632 [1.93333333, 2.03333333, 6.03333333], 3633 [6.93333333, 4.53333333, 4.53333333], 3634 [5.33333333, 5.33333333, 5.33333333]]) 3635 3455 assert num.allclose(stage.vertex_values, 3456 [[2,2,2], 3457 [1.93333333, 2.03333333, 6.03333333], 3458 [6.93333333, 4.53333333, 4.53333333], 3459 [5.33333333, 5.33333333, 5.33333333]]) 3636 3460 3637 3461 def test_balance_deep_and_shallow_tight_SL(self): … … 3639 3463 This test is using Tight Slope Limiters 3640 3464 """ 3465 3641 3466 import copy 3642 3467 … … 3649 3474 3650 3475 points = [a, b, c, d, e, f] 3651 3652 #bac, bce, ecf, dbe 3653 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 3476 # bac, bce, ecf, dbe 3477 elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 3654 3478 3655 3479 domain = Domain(points, elements) 3656 3480 domain.check_integrity() 3657 3481 3658 # Create a deliberate overshoot3482 # Create a deliberate overshoot 3659 3483 domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 3660 domain.set_quantity('elevation', 0) #Flat bed3484 domain.set_quantity('elevation', 0) # Flat bed 3661 3485 stage = domain.quantities['stage'] 3662 3486 3663 ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy3664 3665 # Limit3666 domain.tight_slope_limiters = 1 3487 ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy 3488 3489 # Limit 3490 domain.tight_slope_limiters = 1 3667 3491 domain.distribute_to_vertices_and_edges() 3668 3492 3669 # Assert that quantities are conserved3493 # Assert that quantities are conserved 3670 3494 for k in range(len(domain)): 3671 3495 assert num.allclose (ref_centroid_values[k], 3672 3496 num.sum(stage.vertex_values[k,:])/3) 3673 3497 3674 3675 #Now try with a non-flat bed - closely hugging initial stage in places 3676 #This will create alphas in the range [0, 0.478260, 1] 3498 # Now try with a non-flat bed - closely hugging initial stage in places 3499 # This will create alphas in the range [0, 0.478260, 1] 3677 3500 domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 3678 3501 domain.set_quantity('elevation', [[0,0,0], 3679 [1.8,1.9,5.9],3680 [4.6,0,0],3681 [0,2,4]])3502 [1.8,1.9,5.9], 3503 [4.6,0,0], 3504 [0,2,4]]) 3682 3505 stage = domain.quantities['stage'] 3683 3506 3684 ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy3685 ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy3686 3687 # Limit3688 domain.tight_slope_limiters = 1 3507 ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy 3508 ref_vertex_values = copy.copy(stage.vertex_values[:]) # Copy 3509 3510 # Limit 3511 domain.tight_slope_limiters = 1 3689 3512 domain.distribute_to_vertices_and_edges() 3690 3513 3691 3692 #Assert that all vertex quantities have changed 3514 # Assert that all vertex quantities have changed 3693 3515 for k in range(len(domain)): 3694 #print ref_vertex_values[k,:], stage.vertex_values[k,:]3695 assert not num.allclose (ref_vertex_values[k,:],stage.vertex_values[k,:])3696 # and assert that quantities are still conserved3516 assert not num.allclose(ref_vertex_values[k,:], 3517 stage.vertex_values[k,:]) 3518 # and assert that quantities are still conserved 3697 3519 for k in range(len(domain)): 3698 assert num.allclose (ref_centroid_values[k], 3699 num.sum(stage.vertex_values[k,:])/3) 3700 3701 3702 #Also check that Python and C version produce the same 3703 # No longer applicable if tight_slope_limiters == 1 3704 #print stage.vertex_values 3705 #assert allclose (stage.vertex_values, 3706 # [[2,2,2], 3707 # [1.93333333, 2.03333333, 6.03333333], 3708 # [6.93333333, 4.53333333, 4.53333333], 3709 # [5.33333333, 5.33333333, 5.33333333]]) 3710 3711 3520 assert num.allclose(ref_centroid_values[k], 3521 num.sum(stage.vertex_values[k,:])/3) 3712 3522 3713 3523 def test_balance_deep_and_shallow_Froude(self): … … 3716 3526 This test is using tight slope limiters. 3717 3527 """ 3528 3718 3529 import copy 3719 3530 … … 3726 3537 3727 3538 points = [a, b, c, d, e, f] 3728 3729 # bac, bce, ecf, dbe 3730 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 3539 # bac, bce, ecf, dbe 3540 elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 3731 3541 3732 3542 domain = Domain(points, elements) 3733 3543 domain.check_integrity() 3734 3544 domain.tight_slope_limiters = True 3735 domain.use_centroid_velocities = True 3545 domain.use_centroid_velocities = True 3736 3546 3737 3547 # Create non-flat bed - closely hugging initial stage in places … … 3739 3549 domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 3740 3550 domain.set_quantity('elevation', [[0,0,0], 3741 [1.8,1.999,5.999],3742 [4.6,0,0],3743 [0,2,4]])3551 [1.8,1.999,5.999], 3552 [4.6,0,0], 3553 [0,2,4]]) 3744 3554 3745 3555 # Create small momenta, that nonetheless will generate large speeds … … 3748 3558 domain.set_quantity('ymomentum', 0.0890) 3749 3559 3750 3751 3752 3753 3560 stage = domain.quantities['stage'] 3754 3561 elevation = domain.quantities['elevation'] 3755 3562 xmomentum = domain.quantities['xmomentum'] 3756 ymomentum = domain.quantities['ymomentum'] 3563 ymomentum = domain.quantities['ymomentum'] 3757 3564 3758 3565 # Setup triangle #1 to mimick real Froude explosion observed 3759 3566 # in the Onslow example 13 Nov 2007. 3760 3761 3567 stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953] 3762 elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647] 3568 elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647] 3763 3569 xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066] 3764 3570 ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890] 3765 3571 3766 3572 xmomentum.interpolate() 3767 ymomentum.interpolate() 3768 stage.interpolate() 3573 ymomentum.interpolate() 3574 stage.interpolate() 3769 3575 elevation.interpolate() 3770 3576 … … 3780 3586 v = ymomentum/depth 3781 3587 3782 denom = (depth*g)**0.5 3588 denom = (depth*g)**0.5 3783 3589 Fx = u/denom 3784 3590 Fy = v/denom 3785 3786 3591 3787 3592 # Verify against Onslow example (14 Nov 2007) 3788 3593 assert num.allclose(depth.centroid_values[1], 0.278033) … … 3802 3607 assert num.allclose(Fy.centroid_values[1], 0.193924048435) 3803 3608 3804 3805 3609 # But Froude numbers are huge at some vertices and edges 3806 3610 assert num.allclose(Fx.vertex_values[1,:], [-5.85888475e+01, … … 3819 3623 1.06035244e-01, 3820 3624 3.88346947e+02]) 3821 3822 3625 3626 3823 3627 # The task is now to arrange the limiters such that Froude numbers 3824 3628 # remain under control whil at the same time obeying the conservation 3825 3629 # laws. 3826 3827 3828 ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy 3829 ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy 3630 ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy 3631 ref_vertex_values = copy.copy(stage.vertex_values[:]) # Copy 3830 3632 3831 3633 # Limit (and invoke balance_deep_and_shallow) … … 3834 3636 3835 3637 # Redo derived quantities 3836 depth = stage -elevation3638 depth = stage - elevation 3837 3639 u = xmomentum/depth 3838 3640 v = ymomentum/depth … … 3840 3642 # Assert that all vertex velocities stay within one 3841 3643 # order of magnitude of centroid velocities. 3842 #print u.vertex_values[1,:] 3843 #print u.centroid_values[1] 3844 3845 assert num.alltrue(num.absolute(u.vertex_values[1,:]) <= num.absolute(u.centroid_values[1])*10) 3846 assert num.alltrue(num.absolute(v.vertex_values[1,:]) <= num.absolute(v.centroid_values[1])*10) 3847 3848 denom = (depth*g)**0.5 3644 assert num.alltrue(num.absolute(u.vertex_values[1,:]) <= 3645 num.absolute(u.centroid_values[1])*10) 3646 assert num.alltrue(num.absolute(v.vertex_values[1,:]) <= 3647 num.absolute(v.centroid_values[1])*10) 3648 3649 denom = (depth*g)**0.5 3849 3650 Fx = u/denom 3850 3651 Fy = v/denom 3851 3852 3652 3853 3653 # Assert that Froude numbers are less than max value (TBA) 3854 3654 # at vertices, edges and centroids. 3855 3655 from anuga.config import maximum_froude_number 3856 assert num.alltrue(num.absolute(Fx.vertex_values[1,:]) < maximum_froude_number) 3857 assert num.alltrue(num.absolute(Fy.vertex_values[1,:]) < maximum_froude_number) 3858 3656 3657 assert num.alltrue(num.absolute(Fx.vertex_values[1,:]) < 3658 maximum_froude_number) 3659 assert num.alltrue(num.absolute(Fy.vertex_values[1,:]) < 3660 maximum_froude_number) 3859 3661 3860 3662 # Assert that all vertex quantities have changed 3861 3663 for k in range(len(domain)): 3862 #print ref_vertex_values[k,:], stage.vertex_values[k,:] 3863 assert not num.allclose (ref_vertex_values[k,:], 3864 stage.vertex_values[k,:]) 3865 3664 assert not num.allclose(ref_vertex_values[k,:], 3665 stage.vertex_values[k,:]) 3666 3866 3667 # Assert that quantities are still conserved 3867 3668 for k in range(len(domain)): 3868 assert num.allclose (ref_centroid_values[k], 3869 num.sum(stage.vertex_values[k,:])/3) 3870 3871 3872 3669 assert num.allclose(ref_centroid_values[k], 3670 num.sum(stage.vertex_values[k,:])/3) 3671 3873 3672 return 3874 3673 3875 3674 qwidth = 12 3876 for k in [1]: #range(len(domain)): 3877 print 'Triangle %d (C, V, E)' %k 3878 3879 print 'stage'.ljust(qwidth), stage.centroid_values[k],\ 3880 stage.vertex_values[k,:], stage.edge_values[k,:] 3881 print 'elevation'.ljust(qwidth), elevation.centroid_values[k],\ 3882 elevation.vertex_values[k,:], elevation.edge_values[k,:] 3883 print 'depth'.ljust(qwidth), depth.centroid_values[k],\ 3884 depth.vertex_values[k,:], depth.edge_values[k,:] 3885 print 'xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],\ 3886 xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:] 3887 print 'ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],\ 3888 ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:] 3889 print 'u'.ljust(qwidth),u.centroid_values[k],\ 3890 u.vertex_values[k,:], u.edge_values[k,:] 3891 print 'v'.ljust(qwidth), v.centroid_values[k],\ 3892 v.vertex_values[k,:], v.edge_values[k,:] 3893 print 'Fx'.ljust(qwidth), Fx.centroid_values[k],\ 3894 Fx.vertex_values[k,:], Fx.edge_values[k,:] 3895 print 'Fy'.ljust(qwidth), Fy.centroid_values[k],\ 3896 Fy.vertex_values[k,:], Fy.edge_values[k,:] 3897 3898 3899 3900 3901 3675 for k in [1]: # range(len(domain)): 3676 print 'Triangle %d (C, V, E)' % k 3677 3678 print ('stage'.ljust(qwidth), stage.centroid_values[k], 3679 stage.vertex_values[k,:], stage.edge_values[k,:]) 3680 print ('elevation'.ljust(qwidth), elevation.centroid_values[k], 3681 elevation.vertex_values[k,:], elevation.edge_values[k,:]) 3682 print ('depth'.ljust(qwidth), depth.centroid_values[k], 3683 depth.vertex_values[k,:], depth.edge_values[k,:]) 3684 print ('xmomentum'.ljust(qwidth), xmomentum.centroid_values[k], 3685 xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:]) 3686 print ('ymomentum'.ljust(qwidth), ymomentum.centroid_values[k], 3687 ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:]) 3688 print ('u'.ljust(qwidth),u.centroid_values[k], 3689 u.vertex_values[k,:], u.edge_values[k,:]) 3690 print ('v'.ljust(qwidth), v.centroid_values[k], 3691 v.vertex_values[k,:], v.edge_values[k,:]) 3692 print ('Fx'.ljust(qwidth), Fx.centroid_values[k], 3693 Fx.vertex_values[k,:], Fx.edge_values[k,:]) 3694 print ('Fy'.ljust(qwidth), Fy.centroid_values[k], 3695 Fy.vertex_values[k,:], Fy.edge_values[k,:]) 3902 3696 3903 3697 def test_conservation_1(self): … … 3907 3701 initial condition 3908 3702 """ 3703 3909 3704 from mesh_factory import rectangular 3910 3705 3911 # Create basic mesh3706 # Create basic mesh 3912 3707 points, vertices, boundary = rectangular(6, 6) 3913 3708 3914 # Create shallow water domain3709 # Create shallow water domain 3915 3710 domain = Domain(points, vertices, boundary) 3916 3711 domain.smooth = False 3917 domain.default_order =23918 3919 # IC3712 domain.default_order = 2 3713 3714 # IC 3920 3715 def x_slope(x, y): 3921 3716 return x/3 … … 3934 3729 initial_xmom = domain.quantities['xmomentum'].get_integral() 3935 3730 3936 #print initial_xmom 3937 3938 #Evolution 3939 for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0): 3940 volume = domain.quantities['stage'].get_integral() 3941 assert num.allclose (volume, initial_volume) 3731 # Evolution 3732 for t in domain.evolve(yieldstep=0.05, finaltime=5.0): 3733 volume = domain.quantities['stage'].get_integral() 3734 assert num.allclose(volume, initial_volume) 3942 3735 3943 3736 #I don't believe that the total momentum should be the same … … 3949 3742 os.remove(domain.get_name() + '.sww') 3950 3743 3951 3952 3744 def test_conservation_2(self): 3953 3745 """Test that stage is conserved globally … … 3956 3748 initial condition 3957 3749 """ 3750 3958 3751 from mesh_factory import rectangular 3959 3752 3960 # Create basic mesh3753 # Create basic mesh 3961 3754 points, vertices, boundary = rectangular(6, 6) 3962 3755 3963 # Create shallow water domain3756 # Create shallow water domain 3964 3757 domain = Domain(points, vertices, boundary) 3965 3758 domain.smooth = False 3966 domain.default_order =23967 3968 # IC3759 domain.default_order = 2 3760 3761 # IC 3969 3762 def x_slope(x, y): 3970 3763 return x/3 … … 3972 3765 domain.set_quantity('elevation', x_slope) 3973 3766 domain.set_quantity('friction', 0) 3974 domain.set_quantity('stage', 0.4) #Steady3767 domain.set_quantity('stage', 0.4) # Steady 3975 3768 3976 3769 # Boundary conditions (reflective everywhere) … … 3983 3776 initial_xmom = domain.quantities['xmomentum'].get_integral() 3984 3777 3985 #print initial_xmom 3986 3987 #Evolution 3988 for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0): 3989 volume = domain.quantities['stage'].get_integral() 3990 assert num.allclose (volume, initial_volume) 3778 # Evolution 3779 for t in domain.evolve(yieldstep=0.05, finaltime=5.0): 3780 volume = domain.quantities['stage'].get_integral() 3781 assert num.allclose(volume, initial_volume) 3991 3782 3992 3783 #FIXME: What would we expect from momentum … … 4000 3791 """Test that stage is conserved globally 4001 3792 4002 This one uses a larger grid, convoluted bed, reflective b dries and a suitable4003 initial condition3793 This one uses a larger grid, convoluted bed, reflective boundaries 3794 and a suitable initial condition 4004 3795 """ 3796 4005 3797 from mesh_factory import rectangular 4006 3798 4007 # Create basic mesh3799 # Create basic mesh 4008 3800 points, vertices, boundary = rectangular(2, 1) 4009 3801 4010 # Create shallow water domain3802 # Create shallow water domain 4011 3803 domain = Domain(points, vertices, boundary) 4012 3804 domain.smooth = False … … 4014 3806 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 4015 3807 4016 # IC3808 # IC 4017 3809 def x_slope(x, y): 4018 3810 z = 0*x … … 4028 3820 return z 4029 3821 4030 4031 4032 3822 domain.set_quantity('elevation', x_slope) 4033 3823 domain.set_quantity('friction', 0) … … 4044 3834 4045 3835 import copy 4046 ref_centroid_values =\ 4047 copy.copy(domain.quantities['stage'].centroid_values)4048 4049 #print 'ORG', domain.quantities['stage'].centroid_values 3836 3837 ref_centroid_values = copy.copy(domain.quantities['stage'].\ 3838 centroid_values) 3839 4050 3840 domain.distribute_to_vertices_and_edges() 4051 3841 4052 4053 #print domain.quantities['stage'].centroid_values4054 3842 assert num.allclose(domain.quantities['stage'].centroid_values, 4055 3843 ref_centroid_values) 4056 3844 4057 4058 #Check that initial limiter doesn't violate cons quan 3845 # Check that initial limiter doesn't violate cons quan 4059 3846 assert num.allclose(domain.quantities['stage'].get_integral(), 4060 3847 initial_volume) 4061 3848 4062 # Evolution4063 for t in domain.evolve(yieldstep = 0.05, finaltime =10):3849 # Evolution 3850 for t in domain.evolve(yieldstep=0.05, finaltime=10): 4064 3851 volume = domain.quantities['stage'].get_integral() 4065 #print t, volume, initial_volume4066 3852 assert num.allclose (volume, initial_volume) 4067 3853 … … 4071 3857 """Test that stage is conserved globally 4072 3858 4073 This one uses a larger grid, convoluted bed, reflective b dries and a suitable4074 initial condition3859 This one uses a larger grid, convoluted bed, reflective boundaries 3860 and a suitable initial condition 4075 3861 """ 3862 4076 3863 from mesh_factory import rectangular 4077 3864 4078 # Create basic mesh3865 # Create basic mesh 4079 3866 points, vertices, boundary = rectangular(6, 6) 4080 3867 4081 # Create shallow water domain3868 # Create shallow water domain 4082 3869 domain = Domain(points, vertices, boundary) 4083 3870 domain.smooth = False 4084 domain.default_order =24085 4086 # IC3871 domain.default_order = 2 3872 3873 # IC 4087 3874 def x_slope(x, y): 4088 3875 z = 0*x … … 4093 3880 z[i] = -0.5 4094 3881 if 0.5 <= x[i] < 0.7: 4095 #z[i] = 0.3 #OK with beta == 0.24096 z[i] = 0.34 #OK with beta == 0.04097 #z[i] = 0.35 #Fails after 80 timesteps with an error4098 #of the order 1.0e-53882 #z[i] = 0.3 # OK with beta == 0.2 3883 z[i] = 0.34 # OK with beta == 0.0 3884 #z[i] = 0.35 # Fails after 80 timesteps with an error 3885 # of the order 1.0e-5 4099 3886 if 0.7 <= x[i]: 4100 3887 z[i] = x[i]/3 4101 3888 return z 4102 3889 4103 4104 4105 3890 domain.set_quantity('elevation', x_slope) 4106 3891 domain.set_quantity('friction', 0) … … 4117 3902 4118 3903 import copy 4119 ref_centroid_values =\ 4120 copy.copy(domain.quantities['stage'].centroid_values) 4121 4122 #Test limiter by itself