Changeset 212
- Timestamp:
- Aug 24, 2004, 4:11:25 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/domain.py
r209 r212 38 38 39 39 40 #FIXME: Move these explanations elsewhere 41 40 42 #Create an empty list for explicit forcing terms 41 43 # … … 51 53 # FIXME: How to call and how function should look 52 54 53 self. explicit_forcing_terms = []55 self.forcing_terms = [] 54 56 55 57 … … 64 66 # q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1}) 65 67 66 self.semi_implicit_forcing_terms = []68 ###self.semi_implicit_forcing_terms = [] 67 69 68 70 … … 216 218 %(self.time, self.min_timestep, self.number_of_steps, 217 219 self.number_of_first_order_steps) 218 else: 220 elif self.min_timestep > self.max_timestep: 221 print 'Time = %.4f, steps=%d (%d)'\ 222 %(self.time, self.number_of_steps, 223 self.number_of_first_order_steps) 224 else: 219 225 print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\ 220 226 %(self.time, self.min_timestep, … … 416 422 """ 417 423 418 for f in self.explicit_forcing_terms: 419 420 pass 421 424 for f in self.forcing_terms: 425 f(self) 422 426 423 427 -
inundation/ga/storm_surge/pyvolution/shallow_water.py
r209 r212 33 33 self.minimum_allowed_height = minimum_allowed_height 34 34 35 35 self.forcing_terms.append(gravity) 36 self.forcing_terms.append(manning_friction) 36 37 37 38 … … 280 281 raise 'Unknown order' 281 282 283 #Compute edge values 284 for name in domain.conserved_quantities: 285 Q = domain.quantities[name] 286 Q.interpolate_from_vertices_to_edges() 287 288 282 289 def protect_against_infinitesimal_heights_centroid(domain): 283 290 """Adjust height and momentum at centroid if height is less than … … 333 340 Q = domain.quantities[name] 334 341 Q.extrapolate_first_order() 335 Q.interpolate_from_vertices_to_edges() 342 336 343 337 344 … … 419 426 """ 420 427 421 #FIXME: first and second order migh merge 422 423 from Numeric import minimum, maximum 424 428 #FIXME: first and second order might merge 429 425 430 #Update conserved quantities using straight second order 426 431 for name in domain.conserved_quantities: 427 432 Q = domain.quantities[name] 428 429 433 Q.extrapolate_second_order() 430 Q.limiter() 431 Q.interpolate_from_vertices_to_edges() 434 Q.limit() 432 435 433 436 … … 469 472 470 473 471 hmin = min imum( hv[k, :] )474 hmin = min( hv[k, :] ) 472 475 473 476 #Create alpha in [0,1], where alpha==0 means using shallow … … 500 503 # from extrapolator xmomv and ymomv (deep). 501 504 502 503 ##f##or i in range(3):504 ##### xmomv[k,i] = (1-alpha)*xmomc[k] + alpha*xmomv[k,i];505 506 505 xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]; 507 506 … … 512 511 hc = wc-zc 513 512 for k in range(domain.number_of_elements): 514 hmax = max imum(hv[k,:])515 516 if hmax < minimum_allowed_height:513 hmax = max(hv[k,:]) 514 515 if hmax < domain.minimum_allowed_height: 517 516 #Reset negative heights to bed elevation 518 517 if hc[k] < 0.0: 519 518 wc[k] = zc[k] 520 ###hc[k] = 0.0521 519 for i in range(3): 522 520 if hv[k,i] < 0.0: 523 521 wv[k,i] = zv[k,i] 524 ##hv0 = 0.0;}525 522 526 523 … … 569 566 ######################### 570 567 #Standard forcing terms: 571 568 # 572 569 def gravity(domain): 573 570 """Implement forcing function for bed slope working with … … 629 626 630 627 628 ########################### 629 class Constant_height: 630 """Set an initial condition with constant water height, e.g 631 stage s = z+h 632 """ 633 634 #FIXME: Rethink this way of creating values. 635 636 def __init__(self, W, h): 637 self.W = W 638 self.h = h 639 640 def __call__(self, x, y): 641 if self.W is None: 642 from Numeric import ones, Float 643 return self.h*ones(len(x), Float) 644 else: 645 return self.W(x,y) + self.h 646 647 648 631 649 ############################################## 632 650 #Initialise module -
inundation/ga/storm_surge/pyvolution/test_quantity.py
r209 r212 217 217 218 218 219 # def test_second_order_extrapolation2(self): 220 221 # initialise_consecutive_datastructure(points=6+4, elements=4) 222 223 # a = Point (0.0, 0.0) 224 # b = Point (0.0, 2.0) 225 # c = Point (2.0, 0.0) 226 # d = Point (0.0, 4.0) 227 # e = Point (2.0, 2.0) 228 # f = Point (4.0, 0.0) 229 230 # #Set up for a gradient of (3,1), f(x) = 3x+y 231 # v1 = Volume(b,a,c,array([2.0+2.0/3,0,0])) 232 # v2 = Volume(b,c,e,array([4.0+4.0/3,0,0])) 233 # v3 = Volume(e,c,f,array([8.0+2.0/3,0,0])) 234 # v4 = Volume(d,b,e,array([2.0+8.0/3,0,0])) 235 236 # #Setup neighbour structure 237 # domain = Domain([v1,v2,v3,v4]) 238 # domain.precompute() 239 # domain.check_integrity() 240 241 # #Lets's check first order first, hey 242 # domain.order = 1 243 # domain.limiter = dummy_limiter 244 # distribute_to_vertices_and_edges(domain) 245 246 # assert allclose(v2.conserved_quantities_vertex0, 247 # v2.conserved_quantities_centroid) 248 # assert allclose(v2.conserved_quantities_vertex1, 249 # v2.conserved_quantities_centroid) 250 # assert allclose(v2.conserved_quantities_vertex2, 251 # v2.conserved_quantities_centroid) 252 253 254 # #Flux across right edge of volume 1 255 # #Outward pointing normal vector 256 # from shallow_water import flux_using_stage as flux_function 257 # normals = Volume.normals 258 259 # normal = Vector.coordinates[normals[1][2]] 260 # ql = Volume.conserved_quantities_face2[1] 261 # qr = Volume.conserved_quantities_face1[0] 262 # fl = array([0.,0.]) 263 # fr = array([0.,0.]) 264 # flux0, max_speed = flux_function(normal, ql, qr, fl, fr) 265 266 # #print flux0, max_speed 267 268 # #print 269 # #print v1.conserved_quantities_face0,\ 270 # # v2.conserved_quantities_face0,\ 271 # # v3.conserved_quantities_face0,\ 272 # # v4.conserved_quantities_face0 273 # #print 274 # #edgelengths = Volume.geometric[:,2:] 275 # #print 276 # #print 277 278 # from python_versions import compute_flux 279 # compute_flux(domain, 100) 280 # F1 = Volume.explicit_update 281 282 # from domain import compute_flux 283 # compute_flux(domain, 100) 284 # F2 = Volume.explicit_update 285 286 # assert allclose(F1, F2) 287 288 # #print F1 289 # #print F2 290 291 292 293 # #Gradient of fitted pwl surface 294 # a, b = compute_gradient(v2.id) 295 296 # assert abs(a[0] - 3.0) < epsilon 297 # assert abs(b[0] - 1.0) < epsilon 298 # #assert qmin[0] == 2.0 + 2.0/3 299 # #assert qmax[0] == 8.0 + 2.0/3 300 301 # #And now for the second order stuff 302 # # - the full second order extrapolation 303 # domain.order = 2 304 # domain.limiter = dummy_limiter 305 # distribute_to_vertices_and_edges(domain) 306 307 # assert allclose(v2.conserved_quantities_vertex0[0], 2.0) 308 # assert allclose(v2.conserved_quantities_vertex1[0], 6.0) 309 # assert allclose(v2.conserved_quantities_vertex2[0], 8.0) 310 219 def test_second_order_extrapolation2(self): 220 quantity = Quantity(self.mesh4) 221 222 #Set up for a gradient of (3,1), f(x) = 3x+y 223 quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3], 224 location = 'centroids') 225 226 #Gradients 227 a, b = quantity.compute_gradients() 228 229 #gradient bewteen t0 and t1 is undefined as det==0 230 assert a[0] == 0.0 231 assert b[0] == 0.0 232 #The others are OK 233 for i in range(1,4): 234 assert allclose(a[i], 3.0) 235 assert allclose(b[i], 1.0) 236 237 238 quantity.extrapolate_second_order() 239 240 assert allclose(quantity.vertex_values[1,0], 2.0) 241 assert allclose(quantity.vertex_values[1,1], 6.0) 242 assert allclose(quantity.vertex_values[1,2], 8.0) 311 243 312 244 -
inundation/ga/storm_surge/pyvolution/test_shallow_water.py
r205 r212 614 614 615 615 #Check that all levels are above elevation (within eps) 616 617 618 assert alltrue(alltrue(greater_equal(L,E-epsilon))) 619 620 616 assert alltrue(alltrue(greater_equal(L,E-epsilon))) 617 618 619 620 def test_second_order_flat_bed_onestep(self): 621 622 from mesh_factory import rectangular 623 from shallow_water import Domain, Reflective_boundary,\ 624 Dirichlet_boundary, Constant_height 625 from Numeric import array 626 627 #Create basic mesh 628 points, vertices, boundary = rectangular(6, 6) 629 630 #Create shallow water domain 631 domain = Domain(points, vertices, boundary) 632 domain.smooth = False 633 domain.default_order=2 634 635 # Boundary conditions 636 Br = Reflective_boundary(domain) 637 Bd = Dirichlet_boundary([0.1, 0., 0.]) 638 domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) 639 640 domain.check_integrity() 641 642 #Evolution 643 for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05): 644 pass# domain.write_time() 645 646 #Data from earlier version of pyvolution 647 assert allclose(domain.min_timestep, 0.0396825396825) 648 assert allclose(domain.max_timestep, 0.0396825396825) 649 650 assert allclose(domain.quantities['level'].centroid_values[:12], 651 [0.00171396, 0.02656103, 0.00241523, 0.02656103, 652 0.00241523, 0.02656103, 0.00241523, 0.02656103, 653 0.00241523, 0.02656103, 0.00241523, 0.0272623]) 654 655 domain.distribute_to_vertices_and_edges() 656 assert allclose(domain.quantities['level'].vertex_values[:12,0], 657 [0.0001714, 0.02656103, 0.00024152, 658 0.02656103, 0.00024152, 0.02656103, 659 0.00024152, 0.02656103, 0.00024152, 660 0.02656103, 0.00024152, 0.0272623]) 661 662 assert allclose(domain.quantities['level'].vertex_values[:12,1], 663 [0.00315012, 0.02656103, 0.00024152, 0.02656103, 664 0.00024152, 0.02656103, 0.00024152, 0.02656103, 665 0.00024152, 0.02656103, 0.00040506, 0.0272623]) 666 667 assert allclose(domain.quantities['level'].vertex_values[:12,2], 668 [0.00182037, 0.02656103, 0.00676264, 669 0.02656103, 0.00676264, 0.02656103, 670 0.00676264, 0.02656103, 0.00676264, 671 0.02656103, 0.0065991, 0.0272623]) 672 673 assert allclose(domain.quantities['xmomentum'].centroid_values[:12], 674 [0.00113961, 0.01302432, 0.00148672, 675 0.01302432, 0.00148672, 0.01302432, 676 0.00148672, 0.01302432, 0.00148672 , 677 0.01302432, 0.00148672, 0.01337143]) 678 679 assert allclose(domain.quantities['ymomentum'].centroid_values[:12], 680 [-2.91240050e-004 , 1.22721531e-004, 681 -1.22721531e-004, 1.22721531e-004 , 682 -1.22721531e-004, 1.22721531e-004, 683 -1.22721531e-004 , 1.22721531e-004, 684 -1.22721531e-004, 1.22721531e-004, 685 -1.22721531e-004, -4.57969873e-005]) 686 687 688 689 def test_second_order_flat_bed_moresteps(self): 690 691 from mesh_factory import rectangular 692 from shallow_water import Domain, Reflective_boundary,\ 693 Dirichlet_boundary, Constant_height 694 from Numeric import array 695 696 #Create basic mesh 697 points, vertices, boundary = rectangular(6, 6) 698 699 #Create shallow water domain 700 domain = Domain(points, vertices, boundary) 701 domain.smooth = False 702 domain.default_order=2 703 704 # Boundary conditions 705 Br = Reflective_boundary(domain) 706 Bd = Dirichlet_boundary([0.1, 0., 0.]) 707 domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) 708 709 domain.check_integrity() 710 711 #Evolution 712 for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): 713 pass 714 715 #Data from earlier version of pyvolution 716 #assert allclose(domain.min_timestep, 0.0396825396825) 717 #assert allclose(domain.max_timestep, 0.0396825396825) 718 #print domain.quantities['level'].centroid_values 719 720 721 722 def test_bedslope_problem_first_order(self): 723 724 from mesh_factory import rectangular 725 from shallow_water import Domain, Reflective_boundary, Constant_height 726 from Numeric import array 727 728 #Create basic mesh 729 points, vertices, boundary = rectangular(6, 6) 730 731 #Create shallow water domain 732 domain = Domain(points, vertices, boundary) 733 domain.smooth = False 734 domain.default_order=1 735 736 #Bed-slope and friction 737 def x_slope(x, y): 738 return -x/3 739 740 domain.set_quantity('elevation', x_slope) 741 742 # Boundary conditions 743 Br = Reflective_boundary(domain) 744 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 745 746 #Initial condition 747 domain.set_quantity('level', Constant_height(x_slope, 0.05)) 748 domain.check_integrity() 749 750 #Evolution 751 for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05): 752 pass# domain.write_time() 753 754 #Data from earlier version of pyvolution 755 assert allclose(domain.min_timestep, 0.0423354136479) 756 assert allclose(domain.max_timestep, 0.0479397182296) 757 assert allclose(domain.quantities['level'].centroid_values, 758 [0.00974056, 0.02131643, 0.01373269, 0.02139684, 759 0.01373269, 0.02139684, 0.01373269, 0.02139684, 760 0.01373269, 0.02139684, 0.01386079, 0.0254346, 761 -0.04562579, -0.025253, -0.04154658, -0.02512009, 762 -0.04154658, -0.02512009, -0.04154658, -0.02512009, 763 -0.04154658, -0.02512009, -0.04140825, -0.0210463, 764 -0.10118135, -0.08080855, -0.09710213, -0.08067564, 765 -0.09710213, -0.08067564, -0.09710213, -0.08067564, 766 -0.09710213, -0.08067564, -0.0969638 , -0.07660185, 767 -0.1567369, -0.13636411, -0.15265769, -0.1362312, 768 -0.15265769, -0.1362312, -0.15265769, -0.1362312, 769 -0.15265769, -0.1362312, -0.15251936, -0.13215741, 770 -0.21229246, -0.19191967, -0.20821325, -0.19178675, 771 -0.20821325, -0.19178675, -0.20821325, -0.19178675, 772 -0.20821325, -0.19178675, -0.20807491, -0.18771296, 773 -0.25875433, -0.24722502, -0.25470346, -0.24709274, 774 -0.25470346, -0.24709274, -0.25470346, -0.24709274, 775 -0.25470346, -0.24709274, -0.25460675, -0.24309962]) 776 777 778 779 def test_bedslope_problem_first_order_moresteps(self): 780 781 from mesh_factory import rectangular 782 from shallow_water import Domain, Reflective_boundary, Constant_height 783 from Numeric import array 784 785 #Create basic mesh 786 points, vertices, boundary = rectangular(6, 6) 787 788 #Create shallow water domain 789 domain = Domain(points, vertices, boundary) 790 domain.smooth = False 791 domain.default_order=1 792 793 #Bed-slope and friction 794 def x_slope(x, y): 795 return -x/3 796 797 domain.set_quantity('elevation', x_slope) 798 799 # Boundary conditions 800 Br = Reflective_boundary(domain) 801 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 802 803 #Initial condition 804 domain.set_quantity('level', Constant_height(x_slope, 0.05)) 805 domain.check_integrity() 806 807 #Evolution 808 for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5): 809 pass# domain.write_time() 810 811 #Data from earlier version of pyvolution 812 assert allclose(domain.quantities['level'].centroid_values, 813 [-0.03558044, -0.01758186, -0.03457868, 814 -0.01757463, -0.03400681, -0.01765817, 815 -0.0334924, -0.01763917, -0.03304351, 816 -0.01745604, -0.03259134, -0.01690249, 817 -0.08343238, -0.06844054, -0.08114626, 818 -0.06776272, -0.07985917, -0.0670156, 819 -0.07857347, -0.06636824, -0.07751765, 820 -0.06571005, -0.07694892, -0.06485356, 821 -0.12816762, -0.11399292, -0.12388354, 822 -0.11298543, -0.12146159, -0.11160487, 823 -0.11890573, -0.11016723, -0.11720749, 824 -0.10922083, -0.11748089, -0.10796284, 825 -0.17324244, -0.15860786, -0.16784961, 826 -0.15785693, -0.16501143, -0.1558837, 827 -0.1618709, -0.15384493, -0.15979744, 828 -0.15251944, -0.15976661, -0.15094149, 829 -0.17055387, -0.1934305, -0.17869999, 830 -0.19402574, -0.17570318, -0.19030006, 831 -0.17114482, -0.18720482, -0.16970896, 832 -0.18836511, -0.17655148, -0.19499397, 833 -0.14139494, -0.14492939, -0.14365011, 834 -0.14969524, -0.14343543, -0.14754936, 835 -0.13944213, -0.14337135, -0.13651514, 836 -0.14183568, -0.13623283, -0.14529803]) 837 838 def test_bedslope_problem_second_order(self): 839 840 from mesh_factory import rectangular 841 from shallow_water import Domain, Reflective_boundary, Constant_height 842 from Numeric import array 843 844 #Create basic mesh 845 points, vertices, boundary = rectangular(6, 6) 846 847 #Create shallow water domain 848 domain = Domain(points, vertices, boundary) 849 domain.smooth = False 850 domain.default_order=2 851 852 #Bed-slope and friction 853 def x_slope(x, y): 854 return -x/3 855 856 domain.set_quantity('elevation', x_slope) 857 858 # Boundary conditions 859 Br = Reflective_boundary(domain) 860 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 861 862 #Initial condition 863 domain.set_quantity('level', Constant_height(x_slope, 0.05)) 864 domain.check_integrity() 865 866 assert allclose(domain.quantities['level'].centroid_values, 867 [0.01296296, 0.03148148, 0.01296296, 868 0.03148148, 0.01296296, 0.03148148, 869 0.01296296, 0.03148148, 0.01296296, 870 0.03148148, 0.01296296, 0.03148148, 871 -0.04259259, -0.02407407, -0.04259259, 872 -0.02407407, -0.04259259, -0.02407407, 873 -0.04259259, -0.02407407, -0.04259259, 874 -0.02407407, -0.04259259, -0.02407407, 875 -0.09814815, -0.07962963, -0.09814815, 876 -0.07962963, -0.09814815, -0.07962963, 877 -0.09814815, -0.07962963, -0.09814815, 878 -0.07962963, -0.09814815, -0.07962963, 879 -0.1537037 , -0.13518519, -0.1537037, 880 -0.13518519, -0.1537037, -0.13518519, 881 -0.1537037 , -0.13518519, -0.1537037, 882 -0.13518519, -0.1537037, -0.13518519, 883 -0.20925926, -0.19074074, -0.20925926, 884 -0.19074074, -0.20925926, -0.19074074, 885 -0.20925926, -0.19074074, -0.20925926, 886 -0.19074074, -0.20925926, -0.19074074, 887 -0.26481481, -0.2462963, -0.26481481, 888 -0.2462963, -0.26481481, -0.2462963, 889 -0.26481481, -0.2462963, -0.26481481, 890 -0.2462963, -0.26481481, -0.2462963]) 891 892 893 #print domain.quantities['level'].extrapolate_second_order() 894 #domain.distribute_to_vertices_and_edges() 895 #print domain.quantities['level'].vertex_values[:,0] 896 897 #Evolution 898 for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): 899 #domain.write_time() 900 pass 901 902 903 #print domain.quantities['level'].centroid_values 904 905 #Data from earlier version of pyvolution 906 assert allclose(domain.min_timestep, 0.0376895634803) 907 assert allclose(domain.max_timestep, 0.0415635655309) 908 909 910 #FIXME: Fails 911 #assert allclose(domain.quantities['level'].centroid_values, 912 # [0.00855788, 0.01575204, 0.00994606, 0.01717072, 913 # 0.01005985, 0.01716362, 0.01005985, 0.01716299, 914 # 0.01007098, 0.01736248, 0.01216452, 0.02026776, 915 # -0.04462374, -0.02479045, -0.04199789, -0.0229465, 916 # -0.04184033, -0.02295693, -0.04184013, -0.02295675, 917 # -0.04184486, -0.0228168, -0.04028876, -0.02036486, 918 # -0.10029444, -0.08170809, -0.09772846, -0.08021704, 919 # -0.09760006, -0.08022143, -0.09759984, -0.08022124, 920 # -0.09760261, -0.08008893, -0.09603914, -0.07758209, 921 # -0.15584152, -0.13723138, -0.15327266, -0.13572906, 922 # -0.15314427, -0.13573349, -0.15314405, -0.13573331, 923 # -0.15314679, -0.13560104, -0.15158523, -0.13310701, 924 # -0.21208605, -0.19283913, -0.20955631, -0.19134189, 925 # -0.20942821, -0.19134598, -0.20942799, -0.1913458, 926 # -0.20943005, -0.19120952, -0.20781177, -0.18869401, 927 # -0.25384082, -0.2463294, -0.25047649, -0.24464654, 928 # -0.25031159, -0.24464253, -0.25031112, -0.24464253, 929 # -0.25031463, -0.24454764, -0.24885323, -0.24286438]) 930 621 931 622 932 #-------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.