Changeset 234
- Timestamp:
- Aug 27, 2004, 4:33:24 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 6 edited
-
HUSK.txt (modified) (1 diff)
-
config.py (modified) (1 diff)
-
domain.py (modified) (1 diff)
-
quantity.py (modified) (3 diffs)
-
shallow_water.py (modified) (9 diffs)
-
show_balanced_limiters.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/HUSK.txt
r232 r234 1 2 Perhaps control old/new style with minimum_allowed height 3 1 4 Test friction term! 2 5 -
inundation/ga/storm_surge/pyvolution/config.py
r221 r234 5 5 6 6 default_boundary_tag = 'exterior' 7 8 newstyle = False #Flase meqans that we are equivalent to previous version pyvolution2 7 9 8 10 -
inundation/ga/storm_surge/pyvolution/domain.py
r232 r234 70 70 71 71 #Defaults 72 from config import max_smallsteps, beta 72 from config import max_smallsteps, beta, newstyle 73 73 self.beta = beta 74 self.newstyle = newstyle 74 75 self.default_order = 1 75 76 self.order = self.default_order -
inundation/ga/storm_surge/pyvolution/quantity.py
r229 r234 183 183 x2, y2 = self.domain.centroids[k2] #V2 centroid 184 184 185 #if k < 4:186 # print 'k, k0, k1, k2 = %d, %d, %d, %d' %(k, k0, k1, k2)187 # #print 'k = %d: (%f, %f), (%f, %f), (%f, %f)'\188 # # %(k, x0, y0, x1, y1, x2, y2)189 # print 'k = %d: (%f, %f, %f)' %(k, q0, q1, q2)190 191 185 #Gradient 192 186 a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) … … 220 214 dq[:,i] = qv[:,i] - qc 221 215 222 #Find min and max at this and neighbour's centroids216 #Find min and max of this and neighbour's centroid values 223 217 qmax = zeros(qc.shape, Float) 224 218 qmin = zeros(qc.shape, Float) … … 275 269 276 270 a, b = self.compute_gradients() 277 ##for k in range(4):278 ## print 'Gradients %d: %16.12f, %16.12f' %(k, a[k], b[k])279 271 280 272 V = self.domain.get_vertex_coordinates() -
inundation/ga/storm_surge/pyvolution/shallow_water.py
r232 r234 272 272 def distribute_to_vertices_and_edges(domain): 273 273 274 ##print 'Calling distrib within SW'274 protect_against_negative_heights_centroid(domain) 275 275 if domain.order == 1: 276 #FIXME: This should be cleaned up, but we try to follow277 #pyvolution 2 strictly for now278 protect_against_negative_heights_centroid(domain)279 protect_against_infinitesimal_heights_centroid(domain)280 276 extrapolate_first_order(domain) 281 277 elif domain.order == 2: 282 #protect_against_negative_heights_centroid(domain) 283 extrapolate_second_order(domain) 278 extrapolate_second_order(domain) 284 279 else: 285 280 raise 'Unknown order' 286 281 287 282 #Compute edge values 288 283 for name in domain.conserved_quantities: … … 290 285 Q.interpolate_from_vertices_to_edges() 291 286 292 293 def protect_against_infinitesimal_heights_centroid(domain): 294 """Adjust height and momentum at centroid if height is less than287 def protect_against_infinitesimal_heights(domain): 288 """Protect against infinitesimal heights and high speeds 289 Adjust height and momentum at centroid if height is less than 295 290 minimum allowed height 296 291 """ 297 #FIXME: Used only in first order 298 299 #Water levels at centroids 292 293 #FIXME: May be unnecessary 294 #Shortcuts 295 wv = domain.quantities['level'].vertex_values 296 zv = domain.quantities['elevation'].vertex_values 300 297 wc = domain.quantities['level'].centroid_values 301 302 #Bed elevations at centroids303 298 zc = domain.quantities['elevation'].centroid_values 304 305 #Water depths at centroids 299 hv = wv-zv 306 300 hc = wc - zc 307 301 308 #Momentums at centroids 302 for k in range(domain.number_of_elements): 303 hmax = max(hv[k,:]) 304 305 if hmax < domain.minimum_allowed_height: 306 307 #Reset negative heights to bed elevation 308 if hc[k] < 0.0: 309 print 'C' 310 wc[k] = zc[k] 311 for i in range(3): 312 if hv[k,i] < 0.0: 313 print 'V' 314 wv[k,i] = zv[k,i] 315 316 317 def protect_against_negative_heights_centroid(domain): 318 """Protect against infinitesimal heights and associated high velocities 319 """ 320 321 #FIXME: Change name appropriately 322 323 #Shortcuts 324 wc = domain.quantities['level'].centroid_values 325 zc = domain.quantities['elevation'].centroid_values 309 326 xmomc = domain.quantities['xmomentum'].centroid_values 310 ymomc = domain.quantities['ymomentum'].centroid_values 311 327 ymomc = domain.quantities['ymomentum'].centroid_values 328 hc = wc - zc #Water depths at centroids 329 330 #Update 312 331 for k in range(domain.number_of_elements): 313 #Protect against infinitesimal heights and high velocities 314 if hc[k] < domain.minimum_allowed_height: 315 #Control level and height 332 333 334 if domain.newstyle: 335 if hc[k] < domain.minimum_allowed_height: 336 if hc[k] < 0.0: 337 #Control level and height 338 wc[k] = zc[k]; hc[k] = 0.0 339 340 #Control momentum 341 xmomc[k] = ymomc[k] = 0.0 342 else: 316 343 if hc[k] < 0.0: 344 #Control level and height 317 345 wc[k] = zc[k]; hc[k] = 0.0 318 346 319 #Control momentum 320 xmomc[k] = ymomc[k] = 0.0 321 322 323 def protect_against_negative_heights_centroid(domain): 324 """Adjust height and momentum at centroid if height is less than zero 325 """ 326 327 328 #Water levels at centroids 329 wc = domain.quantities['level'].centroid_values 330 331 #Bed elevations at centroids 332 zc = domain.quantities['elevation'].centroid_values 333 334 #Water depths at centroids 335 hc = wc - zc 336 337 #Momentums at centroids 338 xmomc = domain.quantities['xmomentum'].centroid_values 339 ymomc = domain.quantities['ymomentum'].centroid_values 340 341 for k in range(domain.number_of_elements): 342 #Protect against infinitesimal heights and high velocities 343 if hc[k] < 0.0: 344 #Control level and height 345 wc[k] = zc[k]; hc[k] = 0.0 346 347 #Control momentum 348 xmomc[k] = ymomc[k] = 0.0 349 350 347 #Control momentum 348 xmomc[k] = ymomc[k] = 0.0 351 349 352 350 … … 367 365 """ 368 366 367 #Shortcuts 368 wc = domain.quantities['level'].centroid_values 369 zc = domain.quantities['elevation'].centroid_values 370 xmomc = domain.quantities['xmomentum'].centroid_values 371 ymomc = domain.quantities['ymomentum'].centroid_values 372 hc = wc - zc #Water depths at centroids 373 374 375 if not domain.newstyle: 376 #Protect against infinitesimal heights and high speeds at_centroid 377 #FIXME: Try without when using the second balance function 378 379 for k in range(domain.number_of_elements): 380 #Protect against infinitesimal heights and high velocities 381 if hc[k] < domain.minimum_allowed_height: 382 #Control level and height 383 if hc[k] < 0.0: 384 wc[k] = zc[k]; hc[k] = 0.0# 385 386 #Control momentum 387 xmomc[k] = ymomc[k] = 0.0 388 369 389 370 390 #Update conserved quantities using straight first order … … 373 393 Q.extrapolate_first_order() 374 394 375 376 377 #Water levels at centroids 378 wc = domain.quantities['level'].centroid_values 379 380 #Bed elevations at centroids 381 zc = domain.quantities['elevation'].centroid_values 382 383 #Water depths at centroids 384 hc = wc - zc 385 386 #Water levels at vertices 387 wv = domain.quantities['level'].vertex_values 388 389 #Bed elevations at vertices 390 zv = domain.quantities['elevation'].vertex_values 391 392 #Water depths at vertices 393 hv = wv-zv 394 395 396 #Computed weighted balance between constant levels and and 397 #levels parallel to the bed elevation. 398 for k in range(domain.number_of_elements): 399 400 #Compute maximal variation in bed elevation 401 z_range = max(abs(zv[k,0]-zc[k]), 402 abs(zv[k,1]-zc[k]), 403 abs(zv[k,2]-zc[k])) 404 405 406 #Weighted balance between stage parallel to bed elevation 407 #(wvi = zvi + hc) and constant stage (wvi = wc = zc+hc) 408 #where i=0,1,2 denotes the vertex ids 409 # 410 #It follows that 411 # wvi = (1-alpha)*(zvi+hc) + alpha*(zc+hc) = 412 # (1-alpha)*zvi + alpha*zc + hc = 413 # zvi + hc + alpha*(zc-zvi) 414 # 415 #where alpha in [0,1] and defined as the ratio between hc and 416 #the maximal difference from zc to zv0, zv1 and zv2 417 # 418 #Mathematically the following can be continued on using hc as 419 # wvi = 420 # zvi + hc + alpha*(zc+hc-zvi-hc) = 421 # zvi + hc + alpha*(hvi-hc) 422 #since hvi = zc+hc-zvi in the constant case 423 424 425 if z_range > 0.0: 426 alpha = min(hc[k]/z_range, 1.0) 427 else: 428 alpha = 1.0 429 430 #Update water levels 431 for i in range(3): 432 #FIXME: Use the original first-order one first, then switch 433 wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i]) 434 #wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k]) 435 436 #FIXME: What about alpha weighting of momentum?? 395 396 397 if domain.newstyle: 398 balance_deep_and_shallow(domain) #This will be better 399 else: 400 old_first_order_balancing_of_deep_and_shallow(domain) #Tests will pass 401 402 437 403 438 404 … … 460 426 """ 461 427 462 #FIXME: first and second order might merge 463 464 #Update conserved quantities using straight second order 428 429 #Update conserved quantities using straight second order with limiter 465 430 for name in domain.conserved_quantities: 466 431 Q = domain.quantities[name] 467 432 Q.extrapolate_second_order() 468 469 #print 'y1', Q.vertex_values[1,:] #OK 470 471 #FIXME - like pyvolution 2 ....................... 472 protect_against_negative_heights_centroid(domain) 473 474 #print 'y1', Q.vertex_values[1,:] #OK 475 476 for name in domain.conserved_quantities: 477 Q = domain.quantities[name] 478 Q.limit() 479 480 481 #print 'y1', Q.vertex_values[1,:] #OK 482 483 #Water levels at centroids 433 Q.limit() 434 435 436 balance_deep_and_shallow(domain) 437 ##protect_against_infinitesimal_heights(domain) #FIXME: Probably not necessary 438 439 440 441 def balance_deep_and_shallow(domain): 442 #FIXME: Grap math comments from old_first_order_balance....(below) 443 444 #Shortcuts 484 445 wc = domain.quantities['level'].centroid_values 485 486 #Bed elevations at centroids487 446 zc = domain.quantities['elevation'].centroid_values 488 489 #Water depths at centroids490 447 hc = wc - zc 491 492 #Water levels at vertices 448 493 449 wv = domain.quantities['level'].vertex_values 494 495 #Bed elevations at vertices496 450 zv = domain.quantities['elevation'].vertex_values 497 498 #Water depths at vertices499 451 hv = wv-zv 500 452 501 502 503 453 #Computed linear combination between constant levels and and 504 454 #levels parallel to the bed elevation. … … 530 480 alpha = 1.0 531 481 532 ##if k==1: print 'alpha', alpha #OK533 534 482 #Update water levels 535 483 # (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = … … 551 499 # xmomc and ymomc (shallow) and momentum 552 500 # from extrapolator xmomv and ymomv (deep). 553 554 501 xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]; 555 502 ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]; 556 503 557 #print 'y1', Q.vertex_values[1,:] #OK558 559 #Finally, protect against infinitesimal heights and high speeds560 #Water depths at vertices561 hv = wv-zv562 hc = wc-zc563 for k in range(domain.number_of_elements):564 hmax = max(hv[k,:])565 566 if hmax < domain.minimum_allowed_height:567 #Reset negative heights to bed elevation568 if hc[k] < 0.0:569 wc[k] = zc[k]570 for i in range(3):571 if hv[k,i] < 0.0:572 wv[k,i] = zv[k,i]573 574 504 575 505 … … 674 604 675 605 #Update momentum 676 Xmom. explicit_update[k] += S*uh677 Ymom. explicit_update[k] += S*vh606 Xmom.semi_implicit_update[k] += S*uh 607 Ymom.semi_implicit_update[k] += S*vh 678 608 679 609 … … 847 777 return self.W(x,y) + self.h 848 778 779 #STUFF 780 781 782 def old_first_order_balancing_of_deep_and_shallow(domain): 783 #FIXME: Will be obsolete, but keep comments from this one. 784 785 #Computed weighted balance between constant levels and and 786 #levels parallel to the bed elevation. 787 wc = domain.quantities['level'].centroid_values 788 zc = domain.quantities['elevation'].centroid_values 789 xmomc = domain.quantities['xmomentum'].centroid_values 790 ymomc = domain.quantities['ymomentum'].centroid_values 791 hc = wc - zc #Water depths at centroids 792 793 wv = domain.quantities['level'].vertex_values 794 zv = domain.quantities['elevation'].vertex_values 795 796 for k in range(domain.number_of_elements): 797 798 #Compute maximal variation in bed elevation 799 z_range = max(abs(zv[k,0]-zc[k]), 800 abs(zv[k,1]-zc[k]), 801 abs(zv[k,2]-zc[k])) 802 803 804 #Weighted balance between stage parallel to bed elevation 805 #(wvi = zvi + hc) and constant stage (wvi = wc = zc+hc) 806 #where i=0,1,2 denotes the vertex ids 807 # 808 #It follows that 809 # wvi = (1-alpha)*(zvi+hc) + alpha*(zc+hc) = 810 # (1-alpha)*zvi + alpha*zc + hc = 811 # zvi + hc + alpha*(zc-zvi) 812 # 813 #where alpha in [0,1] and defined as the ratio between hc and 814 #the maximal difference from zc to zv0, zv1 and zv2 815 # 816 #Mathematically the following could be continued on using hc as 817 # wvi = 818 # zvi + hc + alpha*(zc+hc-zvi-hc) = 819 # zvi + hc + alpha*(hvi-hc) 820 #since hvi = zc+hc-zvi in the constant case 821 822 823 if z_range > 0.0: 824 alpha = min(hc[k]/z_range, 1.0) 825 else: 826 alpha = 1.0 827 828 #Update water levels 829 for i in range(3): 830 #FIXME: Use the original first-order one first, then switch 831 wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i]) 832 #wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k]) #Requires updated hv!! 833 834 #FIXME: What about alpha weighting of momentum?? 849 835 850 836 -
inundation/ga/storm_surge/pyvolution/show_balanced_limiters.py
r233 r234 36 36 domain.visualise = False 37 37 domain.visualise = True 38 domain.default_order=2 38 domain.default_order = 2 39 domain.newstyle = True #True is better 39 40 40 41 … … 46 47 print 'Field values' 47 48 domain.set_quantity('elevation', Z) 48 #domain.set_quantity('friction', manning)49 domain.set_quantity('friction', manning) 49 50 50 51 … … 77 78 78 79 from Numeric import allclose 80 79 81 #Evolve 80 for t in domain.evolve(yieldstep = 0.02, finaltime = 0.7): 81 #for t in domain.evolve(yieldstep = 0.1, finaltime = 20): 82 for t in domain.evolve(yieldstep = 0.1, finaltime = 20): 82 83 domain.write_time() 83 print domain.quantities['level'].centroid_values[:6] 84 85 if domain.time == 0.2: 86 print 'Testing' 87 assert allclose(domain.quantities['level'].centroid_values, 88 [ 3.96042007e-002, 5.61081828e-002, 4.66578380e-002, 5.73165329e-002, 89 4.72542001e-002, 5.74770060e-002, 4.74459150e-002, 5.77550805e-002, 90 4.80791695e-002, 5.85754074e-002, 4.90681598e-002, 6.02682664e-002, 91 1.16686827e-002, 1.75422685e-002, 1.17014731e-002, 2.15810992e-002, 92 1.30549421e-002, 2.14416681e-002, 1.31212934e-002, 2.15486277e-002, 93 1.34996488e-002, 2.24053139e-002, 1.50195194e-002, 2.22306851e-002, 94 -7.26762170e-003, -1.35071582e-003, -7.88143638e-003, -2.18165245e-003, 95 -7.81749271e-003, -1.06732807e-003, -7.76399231e-003, -1.00580353e-003, 96 -7.81877765e-003, -9.81086203e-004, -7.42500800e-003, -2.41412070e-004, 97 1.86244453e-001, 8.79324341e-002, 1.86232625e-001, 8.78313615e-002, 98 6.12537452e-002, -3.73125664e-002, -6.37550753e-002, -3.73269705e-002, 99 6.12145063e-002, 8.77700677e-002, 1.86257693e-001, 8.79121535e-002, 100 -4.83328632e-002, 1.18336527e-001, -4.83328400e-002, 1.18337005e-001, 101 -4.83328850e-002, -6.65776472e-003, -1.73331646e-001, -1.31654218e-001, 102 -1.73332232e-001, -6.66097985e-003, -4.83323869e-002, 1.18339536e-001, 103 -2.48333331e-001, -2.31666623e-001, -2.48333332e-001, -2.31666628e-001, 104 -2.48333332e-001, -2.31666627e-001, -2.48333330e-001, -2.31666575e-001, 105 -2.48333330e-001, -2.31666597e-001, -2.48333329e-001, -2.31666584e-001, 106 -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001, 107 -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001, 108 -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001, 109 -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001, 110 -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001, 111 -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001, 112 -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001, 113 -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001, 114 -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001, 115 -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001, 116 -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001, 117 -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001, 118 -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001, 119 -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001, 120 -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001, 121 -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001, 122 -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001, 123 -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001]) 124 84 #print domain.quantities['level'].centroid_values[:6] 125 85 print 'Done' 126 86
Note: See TracChangeset
for help on using the changeset viewer.
