Changeset 3510 for development/pyvolution1d/shallow_water_1d.py
 Timestamp:
 Aug 21, 2006, 10:03:42 AM (17 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

development/pyvolution1d/shallow_water_1d.py
r3425 r3510 64 64 #forcing terms not included in 1d domain ?WHy? 65 65 self.forcing_terms.append(gravity) 66 self.forcing_terms.append(manning_friction)66 #self.forcing_terms.append(manning_friction) 67 67 #print "\nI have Removed forcing terms line 64 1dsw" 68 68 … … 283 283 284 284 285 # Flux computation 286 #def flux_function(normal, ql, qr, zl, zr): 287 def flux_function(normal, ql, qr, zl, zr): 285 def flux_function_unsplit(normal, ql, qr, zl, zr): 288 286 """Compute fluxes between volumes for the shallow water wave equation 289 287 cast in terms of w = h+z using the 'central scheme' as described in … … 376 374 377 375 return edgeflux, max_speed 376 377 def flux_function_split(normal, ql, qr, zl, zr): 378 from config import g, epsilon 379 from math import sqrt 380 from Numeric import array 381 382 #print 'ql',ql 383 384 #Align momentums with xaxis 385 #q_left = rotate(ql, normal, direction = 1) 386 #q_right = rotate(qr, normal, direction = 1) 387 q_left = ql 388 q_left[1] = q_left[1]*normal 389 q_right = qr 390 q_right[1] = q_right[1]*normal 391 392 #z = (zl+zr)/2 #Take average of field values 393 z = 0.5*(zl+zr) #Take average of field values 394 395 w_left = q_left[0] #w=h+z 396 h_left = w_leftz 397 uh_left = q_left[1] 398 399 if h_left < epsilon: 400 u_left = 0.0 #Could have been negative 401 h_left = 0.0 402 else: 403 u_left = uh_left/h_left 404 405 406 w_right = q_right[0] #w=h+z 407 h_right = w_rightz 408 uh_right = q_right[1] 409 410 411 if h_right < epsilon: 412 u_right = 0.0 #Could have been negative 413 h_right = 0.0 414 else: 415 u_right = uh_right/h_right 416 417 #vh_left = q_left[2] 418 #vh_right = q_right[2] 419 420 #soundspeed_left = sqrt(g*h_left) 421 #soundspeed_right = sqrt(g*h_right) 422 423 #Maximal wave speed 424 #s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0) 425 s_max = max(u_left, u_right, 0) 426 427 #Minimal wave speed 428 #s_min = min(u_leftsoundspeed_left, u_rightsoundspeed_right, 0) 429 s_min = min(u_left, u_right, 0) 430 431 #Flux computation 432 433 #flux_left = array([u_left*h_left, 434 # u_left*uh_left + 0.5*g*h_left*h_left]) 435 #flux_right = array([u_right*h_right, 436 # u_right*uh_right + 0.5*g*h_right*h_right]) 437 flux_left = array([u_left*h_left, 438 u_left*uh_left])# + 0.5*g*h_left*h_left]) 439 flux_right = array([u_right*h_right, 440 u_right*uh_right])# + 0.5*g*h_right*h_right]) 441 442 denom = s_maxs_min 443 if denom == 0.0: 444 edgeflux = array([0.0, 0.0]) 445 max_speed = 0.0 446 else: 447 edgeflux = (s_max*flux_left  s_min*flux_right)/denom 448 edgeflux += s_max*s_min*(q_rightq_left)/denom 449 450 edgeflux[1] = edgeflux[1]*normal 451 452 max_speed = max(abs(s_max), abs(s_min)) 453 454 return edgeflux, max_speed 455 378 456 379 457 def compute_fluxes(domain): … … 471 549 #print 'qr',qr 472 550 473 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr) 551 552 if domain.fluxfunc == "unsplit": 553 edgeflux, max_speed = flux_function_unsplit(normal, ql, qr, zl, zr) 554 elif domain.fluxfunc == "split": 555 edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr) 474 556 #print 'edgeflux', edgeflux 475 557 … … 480 562 try: 481 563 #timestep = min(timestep, 0.5*domain.radii[k]/max_speed) 482 #timestep = 0.01483 484 564 timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed) 485 if (timestep < 1e6) & (enter == True):486 #print "domain.order", domain.order487 #domain.write_time()488 print "cell number", k489 print "timestep", timestep490 print 'max_speed', max_speed491 s = domain.quantities['stage']492 s = s.centroid_values493 xm = domain.quantities['xmomentum']494 xm = xm.centroid_values495 print 'h', s[k]496 print 'xm', xm[k]497 print 'u', xm[k]/s[k]498 enter = False499 500 565 except ZeroDivisionError: 501 566 pass … … 504 569 #quantities get updated 505 570 flux /= domain.areas[k] 506 # ADD ABOVE LINE AGAIN 571 507 572 Stage.explicit_update[k] = flux[0] 508 573 Xmom.explicit_update[k] = flux[1] … … 608 673 609 674 #Remove very thin layers of water 610 ##protect_against_infinitesimal_and_negative_heights(domain)675 protect_against_infinitesimal_and_negative_heights(domain) 611 676 612 677 #Extrapolate all conserved quantities … … 1001 1066 #h = Stage.edge_values  Elevation.edge_values 1002 1067 h = Stage.vertex_values  Elevation.vertex_values 1003 v = Elevation.vertex_values 1068 b = Elevation.vertex_values 1069 w = Stage.vertex_values 1004 1070 1005 1071 x = domain.get_vertex_coordinates() … … 1014 1080 x0, x1 = x[k,:] 1015 1081 #z0, z1, z2 = v[k,:] 1016 z0, z1 = v[k,:] 1082 b0, b1 = b[k,:] 1083 1084 w0, w1 = w[k,:] 1085 wx = gradient(x0, x1, w0, w1) 1017 1086 1018 1087 #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) 1019 zx = gradient(x0, x1, z0, z1)1088 bx = gradient(x0, x1, b0, b1) 1020 1089 1021 1090 #Update momentum 1022 xmom[k] += g*zx*avg_h 1091 if domain.fluxfunc == "unsplit": 1092 xmom[k] += g*bx*avg_h 1093 elif domain.fluxfunc == "split": 1094 xmom[k] += g*wx*avg_h 1023 1095 # ymom[k] += g*zy*avg_h 1024 1096
Note: See TracChangeset
for help on using the changeset viewer.