- Timestamp:
- Sep 20, 2009, 6:35:14 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7452 r7519 147 147 # @param use_cache 148 148 # @param verbose 149 # @param evolved_quantities 149 150 # @param full_send_dict 150 151 # @param ghost_recv_dict … … 163 164 use_cache=False, 164 165 verbose=False, 166 evolved_quantities = None, 167 other_quantities = None, 165 168 full_send_dict=None, 166 169 ghost_recv_dict=None, … … 171 174 172 175 # Define quantities for the shallow_water domain 173 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 174 other_quantities = ['elevation', 'friction'] 176 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 177 178 if other_quantities == None: 179 other_quantities = ['elevation', 'friction'] 175 180 176 181 Generic_Domain.__init__(self, … … 179 184 boundary, 180 185 conserved_quantities, 186 evolved_quantities, 181 187 other_quantities, 182 188 tagged_elements, … … 208 214 self.optimise_dry_cells = optimise_dry_cells 209 215 216 self.use_new_mannings = False 210 217 self.forcing_terms.append(manning_friction_implicit) 211 218 self.forcing_terms.append(gravity) … … 225 232 self.use_centroid_velocities = use_centroid_velocities 226 233 234 235 ## 236 # @brief 237 # @param flag 238 def set_new_mannings_function(self, flag=True): 239 """Cludge to allow unit test to pass, but to 240 also introduce new mannings friction function 241 which takes into account the slope of the bed. 242 The flag is tested in the python wrapper 243 mannings_friction_implicit 244 """ 245 if flag: 246 self.use_new_mannings = True 247 else: 248 self.use_new_mannings = False 249 250 251 ## 252 # @brief 253 # @param flag 254 def set_use_edge_limiter(self, flag=True): 255 """Cludge to allow unit test to pass, but to 256 also introduce new edge limiting. The flag is 257 tested in distribute_to_vertices_and_edges 258 """ 259 if flag: 260 self.use_edge_limiter = True 261 else: 262 self.use_edge_limiter = False 263 264 265 227 266 ## 228 267 # @brief … … 488 527 489 528 490 491 492 ##493 # @brief Get the total flow through an arbitrary poly line.494 # @param polyline Representation of desired cross section.495 # @param verbose True if this method is to be verbose.496 # @note 'polyline' may contain multiple sections allowing complex shapes.497 # @note Assume absolute UTM coordinates.498 def old_get_flow_through_cross_section(self, polyline, verbose=False):499 """Get the total flow through an arbitrary poly line.500 501 This is a run-time equivalent of the function with same name502 in data_manager.py503 504 Input:505 polyline: Representation of desired cross section - it may contain506 multiple sections allowing for complex shapes. Assume507 absolute UTM coordinates.508 Format [[x0, y0], [x1, y1], ...]509 510 Output:511 Q: Total flow [m^3/s] across given segments.512 """513 514 # Find all intersections and associated triangles.515 segments = self.get_intersecting_segments(polyline, use_cache=True,516 verbose=verbose)517 518 # Get midpoints519 midpoints = segment_midpoints(segments)520 521 # Make midpoints Geospatial instances522 midpoints = ensure_geospatial(midpoints, self.geo_reference)523 524 # Compute flow525 if verbose:526 log.critical('Computing flow through specified cross section')527 528 # Get interpolated values529 xmomentum = self.get_quantity('xmomentum')530 ymomentum = self.get_quantity('ymomentum')531 532 uh = xmomentum.get_values(interpolation_points=midpoints,533 use_cache=True)534 vh = ymomentum.get_values(interpolation_points=midpoints,535 use_cache=True)536 537 # Compute and sum flows across each segment538 total_flow = 0539 for i in range(len(uh)):540 # Inner product of momentum vector with segment normal [m^2/s]541 normal = segments[i].normal542 normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]543 544 # Flow across this segment [m^3/s]545 segment_flow = normal_momentum*segments[i].length546 547 # Accumulate548 total_flow += segment_flow549 550 return total_flow551 529 552 530 ## … … 589 567 590 568 return cross_section.get_energy_through_cross_section(kind) 591 592 593 ##594 # @brief595 # @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.600 def old_get_energy_through_cross_section(self, polyline,601 kind='total',602 verbose=False):603 """Obtain average energy head [m] across specified cross section.604 605 Inputs:606 polyline: Representation of desired cross section - it may contain607 multiple sections allowing for complex shapes. Assume608 absolute UTM coordinates.609 Format [[x0, y0], [x1, y1], ...]610 kind: Select which energy to compute.611 Options are 'specific' and 'total' (default)612 613 Output:614 E: Average energy [m] across given segments for all stored times.615 616 The average velocity is computed for each triangle intersected by617 the polyline and averaged weighted by segment lengths.618 619 The typical usage of this function would be to get average energy of620 flow in a channel, and the polyline would then be a cross section621 perpendicular to the flow.622 623 #FIXME (Ole) - need name for this energy reflecting that its dimension624 is [m].625 """626 627 from anuga.config import g, epsilon, velocity_protection as h0628 629 # Find all intersections and associated triangles.630 segments = self.get_intersecting_segments(polyline, use_cache=True,631 verbose=verbose)632 633 # Get midpoints634 midpoints = segment_midpoints(segments)635 636 # Make midpoints Geospatial instances637 midpoints = ensure_geospatial(midpoints, self.geo_reference)638 639 # Compute energy640 if verbose: log.critical('Computing %s energy' % kind)641 642 # Get interpolated values643 stage = self.get_quantity('stage')644 elevation = self.get_quantity('elevation')645 xmomentum = self.get_quantity('xmomentum')646 ymomentum = self.get_quantity('ymomentum')647 648 w = stage.get_values(interpolation_points=midpoints, use_cache=True)649 z = elevation.get_values(interpolation_points=midpoints, use_cache=True)650 uh = xmomentum.get_values(interpolation_points=midpoints,651 use_cache=True)652 vh = ymomentum.get_values(interpolation_points=midpoints,653 use_cache=True)654 h = w-z # Depth655 656 # Compute total length of polyline for use with weighted averages657 total_line_length = 0.0658 for segment in segments:659 total_line_length += segment.length660 661 # Compute and sum flows across each segment662 average_energy = 0.0663 for i in range(len(w)):664 # Average velocity across this segment665 if h[i] > epsilon:666 # Use protection against degenerate velocities667 u = uh[i]/(h[i] + h0/h[i])668 v = vh[i]/(h[i] + h0/h[i])669 else:670 u = v = 0.0671 672 speed_squared = u*u + v*v673 kinetic_energy = 0.5*speed_squared/g674 675 if kind == 'specific':676 segment_energy = h[i] + kinetic_energy677 elif kind == 'total':678 segment_energy = w[i] + kinetic_energy679 else:680 msg = 'Energy kind must be either "specific" or "total".'681 msg += ' I got %s' %kind682 683 # Add to weighted average684 weigth = segments[i].length/total_line_length685 average_energy += segment_energy*weigth686 687 return average_energy688 569 689 570 … … 1866 1747 from shallow_water_ext import gravity as gravity_c 1867 1748 1868 xmom = domain.quantities['xmomentum'].explicit_update1869 ymom = domain.quantities['ymomentum'].explicit_update1749 xmom_update = domain.quantities['xmomentum'].explicit_update 1750 ymom_update = domain.quantities['ymomentum'].explicit_update 1870 1751 1871 1752 stage = domain.quantities['stage'] … … 1878 1759 g = domain.g 1879 1760 1880 gravity_c(g, h, z, x, xmom , ymom) #, 1.0e-6)1761 gravity_c(g, h, z, x, xmom_update, ymom_update) #, 1.0e-6) 1881 1762 1882 1763 ## … … 1889 1770 """ 1890 1771 1891 from shallow_water_ext import manning_friction as manning_friction_c 1772 from shallow_water_ext import manning_friction_old 1773 from shallow_water_ext import manning_friction_new 1892 1774 1893 1775 xmom = domain.quantities['xmomentum'] 1894 1776 ymom = domain.quantities['ymomentum'] 1895 1777 1778 x = domain.get_vertex_coordinates() 1779 1896 1780 w = domain.quantities['stage'].centroid_values 1897 z = domain.quantities['elevation']. centroid_values1781 z = domain.quantities['elevation'].vertex_values 1898 1782 1899 1783 uh = xmom.centroid_values … … 1908 1792 g = domain.g 1909 1793 1910 manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) 1794 if domain.use_new_mannings: 1795 manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update) 1796 else: 1797 manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update) 1798 1911 1799 1912 1800 ## … … 1919 1807 """ 1920 1808 1921 from shallow_water_ext import manning_friction as manning_friction_c 1809 from shallow_water_ext import manning_friction_old 1810 from shallow_water_ext import manning_friction_new 1922 1811 1923 1812 xmom = domain.quantities['xmomentum'] 1924 1813 ymom = domain.quantities['ymomentum'] 1925 1814 1815 x = domain.get_vertex_coordinates() 1816 1926 1817 w = domain.quantities['stage'].centroid_values 1927 z = domain.quantities['elevation']. centroid_values1818 z = domain.quantities['elevation'].vertex_values 1928 1819 1929 1820 uh = xmom.centroid_values … … 1938 1829 g = domain.g 1939 1830 1940 manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) 1831 1832 if domain.use_new_mannings: 1833 manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update) 1834 else: 1835 manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update) 1836 1941 1837 1942 1838
Note: See TracChangeset
for help on using the changeset viewer.