Changeset 7352
- Timestamp:
- Aug 11, 2009, 4:38:46 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7350 r7352 557 557 # @note 'polyline' may contain multiple sections allowing complex shapes. 558 558 # @note Assume absolute UTM coordinates. 559 def new_get_energy_through_cross_section(self, polyline, 560 kind='total', 561 verbose=False): 562 """Obtain average energy head [m] across specified cross section. 563 564 Inputs: 565 polyline: Representation of desired cross section - it may contain 566 multiple sections allowing for complex shapes. Assume 567 absolute UTM coordinates. 568 Format [[x0, y0], [x1, y1], ...] 569 kind: Select which energy to compute. 570 Options are 'specific' and 'total' (default) 571 572 Output: 573 E: Average energy [m] across given segments for all stored times. 574 575 The average velocity is computed for each triangle intersected by 576 the polyline and averaged weighted by segment lengths. 577 578 The typical usage of this function would be to get average energy of 579 flow in a channel, and the polyline would then be a cross section 580 perpendicular to the flow. 581 582 #FIXME (Ole) - need name for this energy reflecting that its dimension 583 is [m]. 584 """ 585 586 587 588 cross_section = Cross_section(self, polyline, verbose) 589 590 return cross_section.get_energy_through_cross_section() 591 592 593 ## 594 # @brief 595 # @param polyline Representation of desired cross section. 596 # @param kind Select energy type to compute ('specific' or 'total'). 597 # @param verbose True if this method is to be verbose. 598 # @note 'polyline' may contain multiple sections allowing complex shapes. 599 # @note Assume absolute UTM coordinates. 559 600 def get_energy_through_cross_section(self, polyline, 560 601 kind='total', … … 646 687 return average_energy 647 688 689 648 690 ## 649 691 # @brief Run integrity checks on shallow water domain. … … 679 721 else: 680 722 distribute_using_vertex_limiter(self) 723 724 681 725 682 726 ## … … 2774 2818 2775 2819 2820 ## 2821 # @brief calculate current energy flow through cross section 2822 def get_energy_through_cross_section(self): 2823 """Obtain average energy head [m] across specified cross section. 2824 2825 Output: 2826 E: Average energy [m] across given segments for all stored times. 2827 2828 The average velocity is computed for each triangle intersected by 2829 the polyline and averaged weighted by segment lengths. 2830 2831 The typical usage of this function would be to get average energy of 2832 flow in a channel, and the polyline would then be a cross section 2833 perpendicular to the flow. 2834 2835 #FIXME (Ole) - need name for this energy reflecting that its dimension 2836 is [m]. 2837 """ 2838 2839 # Get interpolated values 2840 stage = self.domain.get_quantity('stage') 2841 elevation = self.domain.get_quantity('elevation') 2842 xmomentum = self.domain.get_quantity('xmomentum') 2843 ymomentum = self.domain.get_quantity('ymomentum') 2844 2845 w = stage.get_values(interpolation_points=midpoints, use_cache=True) 2846 z = elevation.get_values(interpolation_points=midpoints, use_cache=True) 2847 uh = xmomentum.get_values(interpolation_points=midpoints, 2848 use_cache=True) 2849 vh = ymomentum.get_values(interpolation_points=midpoints, 2850 use_cache=True) 2851 h = w-z # Depth 2852 2853 # Compute total length of polyline for use with weighted averages 2854 total_line_length = 0.0 2855 for segment in segments: 2856 total_line_length += segment.length 2857 2858 # Compute and sum flows across each segment 2859 average_energy = 0.0 2860 for i in range(len(w)): 2861 # Average velocity across this segment 2862 if h[i] > epsilon: 2863 # Use protection against degenerate velocities 2864 u = uh[i]/(h[i] + h0/h[i]) 2865 v = vh[i]/(h[i] + h0/h[i]) 2866 else: 2867 u = v = 0.0 2868 2869 speed_squared = u*u + v*v 2870 kinetic_energy = 0.5*speed_squared/g 2871 2872 if kind == 'specific': 2873 segment_energy = h[i] + kinetic_energy 2874 elif kind == 'total': 2875 segment_energy = w[i] + kinetic_energy 2876 else: 2877 msg = 'Energy kind must be either "specific" or "total".' 2878 msg += ' I got %s' %kind 2879 2880 # Add to weighted average 2881 weigth = segments[i].length/total_line_length 2882 average_energy += segment_energy*weigth 2883 2884 return average_energy 2885 2776 2886 2777 2887
Note: See TracChangeset
for help on using the changeset viewer.