 Timestamp:
 Apr 1, 2009, 8:41:37 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/anuga_1d/shallow_water_domain_suggestion2.py
r5827 r6694 53 53 54 54 conserved_quantities = ['stage', 'xmomentum'] 55 other_quantities = ['elevation', 'friction', 'height', 'velocity'] 55 evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity'] 56 other_quantities = ['friction'] 56 57 Generic_Domain.__init__(self, coordinates, boundary, 57 conserved_quantities, other_quantities,58 conserved_quantities, evolved_quantities, other_quantities, 58 59 tagged_elements) 59 60 … … 91 92 self.quantities_to_be_stored = ['stage','xmomentum'] 92 93 93 self.__doc__ = 'shallow_water_domain' 94 self.__doc__ = 'shallow_water_domain_suggestion2' 95 self.check_integrity() 94 96 95 97 … … 148 150 msg = 'Second conserved quantity must be "xmomentum"' 149 151 assert self.conserved_quantities[1] == 'xmomentum', msg 152 153 msg = 'First evolved quantity must be "stage"' 154 assert self.evolved_quantities[0] == 'stage', msg 155 msg = 'Second evolved quantity must be "xmomentum"' 156 assert self.evolved_quantities[1] == 'xmomentum', msg 157 msg = 'Third evolved quantity must be "elevation"' 158 assert self.evolved_quantities[2] == 'elevation', msg 159 msg = 'Fourth evolved quantity must be "height"' 160 assert self.evolved_quantities[3] == 'height', msg 161 msg = 'Fifth evolved quantity must be "velocity"' 162 assert self.evolved_quantities[4] == 'velocity', msg 150 163 151 164 def extrapolate_second_order_sw(self): … … 168 181 distribute_to_vertices_and_edges(self) 169 182 170 def evolve(self, yieldstep = None, finaltime = None, duration = None, 171 skip_initial_step = False): 172 """Specialisation of basic evolve method from parent class 173 """ 174 175 #Call check integrity here rather than from user scripts 176 #self.check_integrity() 177 178 #msg = 'Parameter beta_h must be in the interval [0, 1)' 179 #assert 0 <= self.beta_h < 1.0, msg 180 #msg = 'Parameter beta_w must be in the interval [0, 1)' 181 #assert 0 <= self.beta_w < 1.0, msg 182 183 184 #Initial update of vertex and edge values before any storage 185 #and or visualisation 186 187 #self.distribute_to_vertices_and_edges() ????????????????????????????????? 188 189 #Initialise real time viz if requested 190 #if self.visualise is True and self.time == 0.0: 191 # if self.visualiser is None: 192 # self.initialise_visualiser() 193 # 194 # self.visualiser.update_timer() 195 # self.visualiser.setup_all() 196 197 #Store model data, e.g. for visualisation 198 #if self.store is True and self.time == 0.0: 199 # self.initialise_storage() 200 # #print 'Storing results in ' + self.writer.filename 201 #else: 202 # pass 203 # #print 'Results will not be stored.' 204 # #print 'To store results set domain.store = True' 205 # #FIXME: Diagnostic output should be controlled by 206 # # a 'verbose' flag living in domain (or in a parent class) 207 183 def evolve(self, yieldstep = None, finaltime = None, duration = None, skip_initial_step = False): 208 184 #Call basic machinery from parent class 209 for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration, 210 skip_initial_step): 211 #Real time viz 212 # if self.visualise is True: 213 # self.visualiser.update_all() 214 # self.visualiser.update_timer() 215 216 217 #Store model data, e.g. for subsequent visualisation 218 # if self.store is True: 219 # self.store_timestep(self.quantities_to_be_stored) 220 221 #FIXME: Could maybe be taken from specified list 222 #of 'store every step' quantities 223 224 #Pass control on to outer loop for more specific actions 185 for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration, skip_initial_step): 225 186 yield(t) 226 187 … … 352 313 353 314 354 #We have got h and u at vertex, then the following is the calculation of fluxes!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!315 #We have got h and u at vertex, then the following is the calculation of fluxes! 355 316 soundspeed_left = sqrt(g*h_left) 356 317 soundspeed_right = sqrt(g*h_right) … … 656 617 import sys 657 618 from Numeric import zeros, Float 658 from config import epsilon, h0659 619 660 620 N = domain.number_of_elements 661 621 timestep = float(sys.maxint) 662 #epsilon = domain.epsilon622 epsilon = domain.epsilon 663 623 g = domain.g 664 624 neighbours = domain.neighbours … … 667 627 areas = domain.areas 668 628 629 Stage = domain.quantities['stage'] 630 Xmom = domain.quantities['xmomentum'] 631 Bed = domain.quantities['elevation'] 632 Height = domain.quantities['height'] 633 Velocity = domain.quantities['velocity'] 634 635 636 stage_boundary_values = Stage.boundary_values 637 xmom_boundary_values = Xmom.boundary_values 638 bed_boundary_values = Bed.boundary_values 639 height_boundary_values= Height.boundary_values 640 vel_boundary_values = Velocity.boundary_values 641 642 stage_explicit_update = Stage.explicit_update 643 xmom_explicit_update = Xmom.explicit_update 644 bed_explicit_values = Bed.explicit_update 645 height_explicit_values= Height.explicit_update 646 vel_explicit_values = Velocity.explicit_update 647 648 max_speed_array = domain.max_speed_array 649 650 domain.distribute_to_vertices_and_edges() 651 domain.update_boundary() 652 653 stage_V = Stage.vertex_values 654 xmom_V = Xmom.vertex_values 655 bed_V = Bed.vertex_values 656 height_V= Height.vertex_values 657 vel_V = Velocity.vertex_values 658 659 number_of_elements = len(stage_V) 660 661 #from config import epsilon, h0, h_min 662 # ##Check this: 663 #for i in range(N): 664 # height_V[i] = stage_V[i]  bed_V[i] 665 # if height_V[i] <= h_min: #epsilon : 666 # height_V[i] = 0.0 667 # stage_V[i] = bed_V[i] 668 # xmom_V[i] = 0.0 669 # vel_V[i] = 0.0 670 # else: 671 # vel_V[i] = xmom_V[i]/(height_V[i] + h0/height_V[i]) 672 # ##End of part to be checked. 673 674 675 #print 'neighbours=',neighbours 676 #print 'neighbour_vertices',neighbour_vertices 677 #print 'normals=',normals 678 #print 'areas=',areas 679 from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced 680 domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep, 681 epsilon, 682 g, 683 684 neighbours, 685 neighbour_vertices, 686 normals, 687 areas, 688 689 stage_V, 690 xmom_V, 691 bed_V, 692 height_V, 693 vel_V, 694 695 stage_boundary_values, 696 xmom_boundary_values, 697 bed_boundary_values, 698 height_boundary_values, 699 vel_boundary_values, 700 701 stage_explicit_update, 702 xmom_explicit_update, 703 bed_explicit_values, 704 height_explicit_values, 705 vel_explicit_values, 706 707 number_of_elements, 708 max_speed_array) 709 710 711 # ################################### 712 713 714 715 716 717 718 # Module functions for gradient limiting (distribute_to_vertices_and_edges) 719 720 def distribute_to_vertices_and_edges(domain): 721 """Distribution from centroids to vertices specific to the 722 shallow water wave 723 equation. 724 725 It will ensure that h (wz) is always nonnegative even in the 726 presence of steep bedslopes by taking a weighted average between shallow 727 and deep cases. 728 729 In addition, all conserved quantities get distributed as per either a 730 constant (order==1) or a piecewise linear function (order==2). 731 732 FIXME: more explanation about removal of artificial variability etc 733 734 Precondition: 735 All quantities defined at centroids and bed elevation defined at 736 vertices. 737 738 Postcondition 739 Conserved quantities defined at vertices 740 741 """ 742 743 #from config import optimised_gradient_limiter 744 745 #Remove very thin layers of water 746 #protect_against_infinitesimal_and_negative_heights(domain) 747 748 import sys 749 from Numeric import zeros, Float 750 from config import epsilon, h0, h_min 751 752 N = domain.number_of_elements 753 754 #Shortcuts 669 755 Stage = domain.quantities['stage'] 670 Xmom 671 Bed 756 Xmom = domain.quantities['xmomentum'] 757 Bed = domain.quantities['elevation'] 672 758 Height = domain.quantities['height'] 673 759 Velocity = domain.quantities['velocity'] 674 760 761 #Arrays 675 762 w_C = Stage.centroid_values 676 763 uh_C = Xmom.centroid_values 677 z_C = Bed.centroid_values 764 z_C = Bed.centroid_values 678 765 h_C = Height.centroid_values 679 766 u_C = Velocity.centroid_values … … 681 768 for i in range(N): 682 769 h_C[i] = w_C[i]  z_C[i] 683 if h_C[i] < 770 if h_C[i] <= h_min: #epsilon : 684 771 h_C[i] = 0.0 685 772 w_C[i] = z_C[i] 686 773 uh_C[i] = 0.0 687 688 for i in range(len(h_C)): 689 if h_C[i] < epsilon: 690 u_C[i] = 0.0 #Could have been negative 691 h_C[i] = 0.0 774 u_C[i] = 0.0 692 775 else: 693 694 695 for name in [' height', 'velocity']:776 u_C[i] = uh_C[i]/(h_C[i] + h0/h_C[i]) 777 778 for name in ['stage', 'height', 'velocity']: 696 779 Q = domain.quantities[name] 697 780 if domain.order == 1: … … 705 788 raise 'Unknown order' 706 789 707 w_V = domain.quantities['stage'].vertex_values 708 h_V = domain.quantities['height'].vertex_values 709 u_V = domain.quantities['velocity'].vertex_values 710 uh_V = domain.quantities['xmomentum'].vertex_values 711 z_V = w_V  h_V 712 790 stage_V = domain.quantities['stage'].vertex_values 791 bed_V = domain.quantities['elevation'].vertex_values 792 h_V = domain.quantities['height'].vertex_values 793 u_V = domain.quantities['velocity'].vertex_values 794 xmom_V = domain.quantities['xmomentum'].vertex_values 795 796 bed_V[:] = stage_V  h_V 797 xmom_V[:] = u_V * h_V 713 798 714 stage_boundary_values = Stage.boundary_values 715 xmom_boundary_values = Xmom.boundary_values 716 stage_explicit_update = Stage.explicit_update 717 xmom_explicit_update = Xmom.explicit_update 718 max_speed_array = domain.max_speed_array 719 #number_of_elements = len(w_V) 799 return 800 720 801 721 #domain.distribute_to_vertices_and_edges()722 #domain.update_boundary()723 #stage_V = Stage.vertex_values724 #xmom_V = Xmom.vertex_values725 #bed_V = Bed.vertex_values726 727 728 from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced #from comp_flux_ext import compute_fluxes_ext729 730 domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,731 epsilon,732 g,733 neighbours,734 neighbour_vertices,735 normals,736 areas,737 w_V,738 uh_V,739 z_V,740 stage_boundary_values,741 xmom_boundary_values,742 stage_explicit_update,743 xmom_explicit_update,744 N,745 max_speed_array)746 747 # ###################################748 749 750 751 752 753 754 # Module functions for gradient limiting (distribute_to_vertices_and_edges)755 756 def distribute_to_vertices_and_edges(domain):757 """Distribution from centroids to vertices specific to the758 shallow water wave759 equation.760 761 It will ensure that h (wz) is always nonnegative even in the762 presence of steep bedslopes by taking a weighted average between shallow763 and deep cases.764 765 In addition, all conserved quantities get distributed as per either a766 constant (order==1) or a piecewise linear function (order==2).767 768 FIXME: more explanation about removal of artificial variability etc769 770 Precondition:771 All quantities defined at centroids and bed elevation defined at772 vertices.773 774 Postcondition775 Conserved quantities defined at vertices776 777 """778 779 #from config import optimised_gradient_limiter780 781 #Remove very thin layers of water782 protect_against_infinitesimal_and_negative_heights(domain)783 784 import sys785 from Numeric import zeros, Float786 787 N = domain.number_of_elements788 789 #Shortcuts790 Stage = domain.quantities['stage']791 Xmom = domain.quantities['xmomentum']792 Bed = domain.quantities['elevation']793 Height = domain.quantities['height']794 Velocity = domain.quantities['velocity']795 796 #Arrays797 w_C = Stage.centroid_values798 uh_C = Xmom.centroid_values799 z_C = Bed.centroid_values800 h_C = Height.centroid_values801 u_C = Velocity.centroid_values802 803 #print id(h_C)804 for i in range(N):805 h_C[i] = w_C[i]  z_C[i] #This is h at centroids: OK806 807 from config import epsilon, h0808 809 for i in range(len(h_C)):810 if h_C[i] < epsilon:811 u_C[i] = 0.0 #Could have been negative812 h_C[i] = 0.0813 else:814 u_C[i] = uh_C[i]/(h_C[i] + h0/h_C[i])815 816 for name in ['stage', 'velocity', 'height' ]:817 Q = domain.quantities[name]818 if domain.order == 1:819 Q.extrapolate_first_order()820 elif domain.order == 2:821 #print "add extrapolate second order to shallow water"822 #if name != 'height':823 Q.extrapolate_second_order()824 #Q.limit()825 else:826 raise 'Unknown order'827 828 stage = domain.quantities['stage'].vertex_values #This is w at vertices: OK829 h_V = domain.quantities['height'].vertex_values830 bed = stage  h_V831 u_V = domain.quantities['velocity'].vertex_values832 xmom_V = domain.quantities['xmomentum'].vertex_values833 834 #h_V[:] = stage  bed835 xmom_V[:] = u_V * h_V836 return stage, xmom_V, bed, h_V, u_V837 802 # 838 803 def protect_against_infinitesimal_and_negative_heights(domain): … … 1064 1029 #self.xmom = domain.quantities['xmomentum'].edge_values 1065 1030 #self.ymom = domain.quantities['ymomentum'].edge_values 1066 self.normals = domain.normals 1067 self.stage = domain.quantities['stage'].vertex_values 1068 self.xmom = domain.quantities['xmomentum'].vertex_values 1031 self.normals = domain.normals 1032 self.stage = domain.quantities['stage'].vertex_values 1033 self.xmom = domain.quantities['xmomentum'].vertex_values 1034 self.bed = domain.quantities['elevation'].vertex_values 1035 self.height = domain.quantities['height'].vertex_values 1036 self.velocity = domain.quantities['velocity'].vertex_values 1069 1037 1070 1038 from Numeric import zeros, Float 1071 1039 #self.conserved_quantities = zeros(3, Float) 1072 self.conserved_quantities = zeros(2, Float) 1040 #self.conserved_quantities = zeros(2, Float) 1041 self.evolved_quantities = zeros(5, Float) 1073 1042 1074 1043 def __repr__(self): 1075 1044 return 'Reflective_boundary' 1076 1045 1077 1046 """ 1078 1047 def evaluate(self, vol_id, edge_id): 1079 """Reflective boundaries reverses the outward momentum1080 of the volume they serve.1081 """1082 1083 q = self. conserved_quantities1048 #Reflective boundaries reverses the outward momentum 1049 #of the volume they serve. 1050 1051 #q = self.conserved_quantities 1052 q = self.evolved_quantities 1084 1053 q[0] = self.stage[vol_id, edge_id] 1085 1054 q[1] = self.xmom[vol_id, edge_id] … … 1101 1070 1102 1071 return q 1072 """ 1073 def evaluate(self, vol_id, edge_id): 1074 """Reflective boundaries reverses the outward momentum 1075 of the volume they serve. 1076 """ 1077 q = self.evolved_quantities 1078 q[0] = self.stage[vol_id, edge_id] 1079 q[1] = self.xmom[vol_id, edge_id] 1080 q[2] = self.bed[vol_id, edge_id] 1081 q[3] = self.height[vol_id, edge_id] 1082 q[4] = self.velocity[vol_id, edge_id] 1083 return q 1084 1085 1086 1103 1087 1104 1088 class Dirichlet_boundary(Boundary):
Note: See TracChangeset
for help on using the changeset viewer.