Changeset 215
- Timestamp:
- Aug 24, 2004, 6:11:26 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/domain.py
r212 r215 438 438 timestep = self.timestep 439 439 440 #Reset semi_implicit_updates441 #(explicit updates are already set by compute_fluyzes)442 for i, name in enumerate(self.conserved_quantities):443 Q = self.quantities[name]444 Q.semi_implicit_update[:] = 0.445 446 440 #Compute forcing terms 447 441 self.compute_forcing_terms() 448 442 449 443 #Update conserved_quantities from explicit updates 450 for i, name in enumerate(self.conserved_quantities):444 for name in self.conserved_quantities: 451 445 Q = self.quantities[name] 452 446 Q.update(timestep) 447 448 #Clean up 449 Q.semi_implicit_update[:] = 0.0 450 Q.explicit_update[:] = 0.0 #Not necessary as fluxes will set it 453 451 454 452 -
inundation/ga/storm_surge/pyvolution/quantity.py
r209 r215 328 328 msg = 'Invalid location: %s' %location 329 329 raise msg 330 331 if X is None: 332 msg = 'Given values are None' 333 raise msg 330 334 331 335 import types -
inundation/ga/storm_surge/pyvolution/shallow_water.py
r212 r215 586 586 for k in range(domain.number_of_elements): 587 587 avg_h = sum( h[k,:] )/3 588 588 589 589 #Compute bed slope 590 590 x0, y0, x1, y1, x2, y2 = V[k,:] 591 591 z0, z1, z2 = Elevation.vertex_values[k,:] 592 592 593 zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) 593 594 … … 595 596 Xmom.explicit_update[k] += -g*zx*avg_h 596 597 Ymom.explicit_update[k] += -g*zy*avg_h 597 598 598 599 599 600 def manning_friction(domain): … … 627 628 628 629 ########################### 630 ########################### 631 #Geometries 632 633 634 #FIXME: Rethink this way of creating values. 635 636 637 class Weir: 638 """Set a bathymetry for weir with a hole and a downstream gutter 639 x,y are assumed to be in the unit square 640 """ 641 642 def __init__(self, stage): 643 self.inflow_stage = stage 644 645 def __call__(self, x, y): 646 from Numeric import zeros, Float 647 from math import sqrt 648 649 N = len(x) 650 assert N == len(y) 651 652 z = zeros(N, Float) 653 for i in range(N): 654 z[i] = -x[i]/2 #General slope 655 656 #Flattish bit to the left 657 if x[i] < 0.3: 658 z[i] = -x[i]/10 659 660 #Weir 661 if x[i] >= 0.3 and x[i] < 0.4: 662 z[i] = -x[i]+0.9 663 664 #Dip 665 x0 = 0.6 666 #depth = -1.3 667 depth = -1.0 668 #plateaux = -0.9 669 plateaux = -0.6 670 if y[i] < 0.7: 671 if x[i] > x0 and x[i] < 0.9: 672 z[i] = depth 673 674 #RHS plateaux 675 if x[i] >= 0.9: 676 z[i] = plateaux 677 678 679 elif y[i] >= 0.7 and y[i] < 1.5: 680 #Restrict and deepen 681 if x[i] >= x0 and x[i] < 0.8: 682 z[i] = depth-(y[i]/3-0.3) 683 #z[i] = depth-y[i]/5 684 #z[i] = depth 685 elif x[i] >= 0.8: 686 #RHS plateaux 687 z[i] = plateaux 688 689 elif y[i] >= 1.5: 690 if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2: 691 #Widen up and stay at constant depth 692 z[i] = depth-1.5/5 693 elif x[i] >= 0.8 + (y[i]-1.5)/1.2: 694 #RHS plateaux 695 z[i] = plateaux 696 697 698 #Hole in weir (slightly higher than inflow condition) 699 if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4: 700 z[i] = -x[i]+self.inflow_stage + 0.02 701 702 #Channel behind weir 703 x0 = 0.5 704 if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4: 705 z[i] = -x[i]+self.inflow_stage + 0.02 706 707 if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4: 708 #Flatten it out towards the end 709 z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5 710 711 #Hole to the east 712 x0 = 1.1; y0 = 0.35 713 #if x[i] < -0.2 and y < 0.5: 714 if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2: 715 z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0 716 717 #Tiny channel draining hole 718 if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6: 719 z[i] = -0.9 #North south 720 721 if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65: 722 z[i] = -1.0 + (x[i]-0.9)/3 #East west 723 724 725 726 #Stuff not in use 727 728 #Upward slope at inlet to the north west 729 #if x[i] < 0.0: # and y[i] > 0.5: 730 # #z[i] = -y[i]+0.5 #-x[i]/2 731 # z[i] = x[i]/4 - y[i]**2 + 0.5 732 733 #Hole to the west 734 #x0 = -0.4; y0 = 0.35 # center 735 #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2: 736 # z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2 737 738 739 740 741 742 return z/2 743 744 class Weir_simple: 745 """Set a bathymetry for weir with a hole and a downstream gutter 746 x,y are assumed to be in the unit square 747 """ 748 749 def __init__(self, stage): 750 self.inflow_stage = stage 751 752 def __call__(self, x, y): 753 from Numeric import zeros, Float 754 755 N = len(x) 756 assert N == len(y) 757 758 z = zeros(N, Float) 759 for i in range(N): 760 z[i] = -x[i] #General slope 761 762 #Flat bit to the left 763 if x[i] < 0.3: 764 z[i] = -x[i]/10 #General slope 765 766 #Weir 767 if x[i] > 0.3 and x[i] < 0.4: 768 z[i] = -x[i]+0.9 769 770 #Dip 771 if x[i] > 0.6 and x[i] < 0.9: 772 z[i] = -x[i]-0.5 #-y[i]/5 773 774 #Hole in weir (slightly higher than inflow condition) 775 if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4: 776 z[i] = -x[i]+self.inflow_stage + 0.05 777 778 779 return z/2 780 781 782 629 783 class Constant_height: 630 784 """Set an initial condition with constant water height, e.g 631 785 stage s = z+h 632 786 """ 633 634 #FIXME: Rethink this way of creating values.635 636 787 def __init__(self, W, h): 637 788 self.W = W -
inundation/ga/storm_surge/pyvolution/test_domain.py
r195 r215 165 165 166 166 167 def test_distribute (self):167 def test_distribute_first_order(self): 168 168 """Domain implements a default first order gradient limiter 169 169 """ … … 206 206 207 207 208 #print domain.quantities['level'].centroid_values209 208 domain.distribute_to_vertices_and_edges() 210 209 … … 244 243 245 244 246 domain.set_quantity('level', [[1,2,3], [5,5,5], 247 [0,0,9], [-6, 3, 3]]) 248 249 domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], 250 [3,3,3], [4, 4, 4]]) 251 252 domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], 253 [30,30,30], [40, 40, 40]]) 254 245 domain.set_quantity('level', [1,2,3,4], 'centroids') 246 domain.set_quantity('xmomentum', [1,2,3,4], 'centroids') 247 domain.set_quantity('ymomentum', [1,2,3,4], 'centroids') 248 255 249 256 250 #Assign some values to update vectors 257 #domain.explicit_updates[] 251 #Set explicit_update 252 253 for name in domain.conserved_quantities: 254 domain.quantities[name].explicit_update = array([4.,3.,2.,1.]) 255 domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.]) 256 257 258 #Update with given timestep (assuming no other forcing terms) 259 domain.timestep = 0.1 260 domain.update_conserved_quantities() 261 262 x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) 263 x /= array( [.9,.9,.9,.9] ) 264 265 for name in domain.conserved_quantities: 266 assert allclose(domain.quantities[name].centroid_values, x) 267 258 268 259 269 -
inundation/ga/storm_surge/pyvolution/test_shallow_water.py
r212 r215 534 534 pass 535 535 536 537 536 ##################################################### 537 def test_gravity(self): 538 #Assuming no friction 539 540 from config import g 541 542 a = [0.0, 0.0] 543 b = [0.0, 2.0] 544 c = [2.0, 0.0] 545 d = [0.0, 4.0] 546 e = [2.0, 2.0] 547 f = [4.0, 0.0] 548 549 points = [a, b, c, d, e, f] 550 #bac, bce, ecf, dbe 551 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 552 553 domain = Domain(points, vertices) 554 555 #Set up for a gradient of (3,0) at mid triangle 556 def slope(x, y): 557 return 3*x 558 559 h = 0.1 560 def level(x,y): 561 return slope(x,y)+h 562 563 domain.set_quantity('elevation', slope) 564 domain.set_quantity('level', level) 565 566 for name in domain.conserved_quantities: 567 assert allclose(domain.quantities[name].explicit_update, 0) 568 assert allclose(domain.quantities[name].semi_implicit_update, 0) 569 570 domain.compute_forcing_terms() 571 572 assert allclose(domain.quantities['level'].explicit_update, 0) 573 assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3) 574 assert allclose(domain.quantities['ymomentum'].explicit_update, 0) 575 576 577 578 ##################################################### 538 579 def test_first_order_extrapolator_const_z(self): 580 581 a = [0.0, 0.0] 582 b = [0.0, 2.0] 583 c = [2.0, 0.0] 584 d = [0.0, 4.0] 585 e = [2.0, 2.0] 586 f = [4.0, 0.0] 587 588 points = [a, b, c, d, e, f] 589 #bac, bce, ecf, dbe 590 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 591 592 domain = Domain(points, vertices) 593 val0 = 2.+2.0/3 594 val1 = 4.+4.0/3 595 val2 = 8.+2.0/3 596 val3 = 2.+8.0/3 597 598 zl=zr=-3.75 #Assume constant bed (must be less than level) 599 domain.set_quantity('elevation', zl*ones( (4,3) )) 600 domain.set_quantity('level', [[val0, val0-1, val0-2], 601 [val1, val1+1, val1], 602 [val2, val2-2, val2], 603 [val3-0.5, val3, val3]]) 604 605 606 607 domain.order = 1 608 domain.distribute_to_vertices_and_edges() 609 610 #Check that centroid values were distributed to vertices 611 C = domain.quantities['level'].centroid_values 612 for i in range(3): 613 assert allclose( domain.quantities['level'].vertex_values[:,i], C) 614 615 616 def test_first_order_limiter_variable_z(self): 617 #Check that first order limiter follows bed_slope 618 from Numeric import alltrue, greater_equal 619 from config import epsilon 539 620 540 621 a = [0.0, 0.0] … … 555 636 val3 = 2.+8.0/3 556 637 557 zl=zr=-3.75 #Assume constant bed (must be less than level)558 domain.set_quantity('elevation', zl*ones( (4,3) ))559 domain.set_quantity('level', [[val0, val0-1, val0-2],560 [val1, val1+1, val1],561 [val2, val2-2, val2],562 [val3-0.5, val3, val3]])563 564 565 566 domain.order = 1567 domain.distribute_to_vertices_and_edges()568 569 #Check that centroid values were distributed to vertices570 C = domain.quantities['level'].centroid_values571 for i in range(3):572 assert allclose( domain.quantities['level'].vertex_values[:,i], C)573 574 575 def test_first_order_limiter_variable_z(self):576 #Check that first order limiter follows bed_slope577 from Numeric import alltrue, greater_equal578 from config import epsilon579 580 a = [0.0, 0.0]581 b = [0.0, 2.0]582 c = [2.0,0.0]583 d = [0.0, 4.0]584 e = [2.0, 2.0]585 f = [4.0,0.0]586 587 points = [a, b, c, d, e, f]588 #bac, bce, ecf, dbe589 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]590 591 domain = Domain(points, vertices)592 val0 = 2.+2.0/3593 val1 = 4.+4.0/3594 val2 = 8.+2.0/3595 val3 = 2.+8.0/3596 597 638 domain.set_quantity('elevation', [[0,0,0], [6,0,0], 598 639 [6,6,6], [6,6,6]]) … … 616 657 assert alltrue(alltrue(greater_equal(L,E-epsilon))) 617 658 659 660 ##################################################### 661 def test_distribute_near_bed(self): 662 #Assuming no friction 663 664 from config import g 665 666 a = [0.0, 0.0] 667 b = [0.0, 2.0] 668 c = [2.0, 0.0] 669 d = [0.0, 4.0] 670 e = [2.0, 2.0] 671 f = [4.0, 0.0] 672 673 points = [a, b, c, d, e, f] 674 #bac, bce, ecf, dbe 675 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 676 677 domain = Domain(points, vertices) 678 domain.order = 2 679 680 #Set up for a gradient of (10,0) at mid triangle 681 def slope(x, y): 682 return 10*x 683 684 h = 0.1 685 def level(x,y): 686 return slope(x,y)+h 687 688 domain.set_quantity('elevation', slope) 689 domain.set_quantity('level', level, 'centroids') 690 691 E = domain.quantities['elevation'].vertex_values 692 L = domain.quantities['level'].vertex_values 693 694 #print E 695 domain.distribute_to_vertices_and_edges() 696 #print L 697 698 #FIXME: HERTIL 24/8/4. Write similar test using pyvolution2 and 699 #get test data from that 618 700 619 701 … … 850 932 domain.default_order=2 851 933 852 #Bed-slope and friction 934 #Bed-slope and friction at vertices (and interpolated elsewhere) 853 935 def x_slope(x, y): 854 936 return -x/3
Note: See TracChangeset
for help on using the changeset viewer.