Changeset 6059


Ignore:
Timestamp:
Dec 11, 2008, 1:45:58 PM (15 years ago)
Author:
ole
Message:

Computed energy at point instead of polyline in Culver_flow_energy

Location:
anuga_core/source/anuga/culvert_flows
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

    r6055 r6059  
    2121     Transmissive_boundary, Time_boundary
    2222
    23 from anuga.culvert_flows.culvert_class import Culvert_flow
     23from anuga.culvert_flows.culvert_class import Culvert_flow_rating, Culvert_flow_energy
    2424from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
    2525     
     
    104104
    105105
    106 culvert = Culvert_flow(domain,
     106culvert_rating = Culvert_flow_rating(domain,
    107107                       culvert_description_filename='example_rating_curve.csv',
    108108                       end_point0=[9.0, 2.5],
     
    111111
    112112
    113 #culvert = Culvert_flow(domain,
    114 #                       label='Culvert No. 1',
    115 #                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
    116 #                       end_point0=[9.0, 2.5],
    117 #                       end_point1=[13.0, 2.5],
    118 #                       width=1.20,height=0.75,
    119 #                       culvert_routine=boyd_generalised_culvert_model,       
    120 #                       number_of_barrels=1,
    121 #                       update_interval=2,
    122 #                       verbose=True)
     113culvert_energy = Culvert_flow_energy(domain,
     114                       label='Culvert No. 1',
     115                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
     116                       end_point0=[9.0, 2.5],
     117                       end_point1=[13.0, 2.5],
     118                       width=1.20,height=0.75,
     119                       culvert_routine=boyd_generalised_culvert_model,       
     120                       number_of_barrels=1,
     121                       update_interval=2,
     122                       verbose=True)
    123123
    124 domain.forcing_terms.append(culvert)
     124domain.forcing_terms.append(culvert_energy)
    125125
    126126#------------------------------------------------------------------------------
  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r6055 r6059  
    8989   
    9090
    91 class Culvert_flow:
     91class Culvert_flow_rating:
    9292    """Culvert flow - transfer water from one hole to another
    9393   
     
    528528            #import sys; sys.exit()                           
    529529
     530           
     531        # Compute the average point for enquiry
     532        enquiry_point0 = sum(P['enquiry_polygon0'][:2])/2
     533        enquiry_point1 = sum(P['enquiry_polygon1'][:2])/2         
     534       
     535        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
     536        self.enquiry_indices = []                 
     537        for point in self.enquiry_points:
     538            # Find nearest centroid
     539            N = len(domain)   
     540            points = domain.get_centroid_coordinates(absolute=True)
     541
     542            # Calculate indices in exchange area for this forcing term
     543           
     544            triangle_id = min_dist = sys.maxint
     545            for k in range(N):
     546                x, y = points[k,:] # Centroid
     547
     548                c = point                               
     549                distance = (x-c[0])**2+(y-c[1])**2
     550                if distance < min_dist:
     551                    min_dist = distance
     552                    triangle_id = k
     553
     554                   
     555            if triangle_id < sys.maxint:
     556                msg = 'found triangle with centroid (%f, %f)'\
     557                    %tuple(points[triangle_id, :])
     558                msg += ' for point (%f, %f)' %tuple(point)
     559               
     560                self.enquiry_indices.append(triangle_id)
     561            else:
     562                msg = 'Triangle not found for point (%f, %f)' %point
     563                raise Exception, msg
     564       
     565                         
     566
     567           
     568           
    530569
    531570        # Check that all polygons lie within the mesh.
     
    688727                   
    689728                # Average measures of energy in front of this opening
    690                 polyline = self.enquiry_polylines[i]
    691                 #print 't = %.4f, opening=%d,' %(domain.time, i),
    692                 opening.total_energy = domain.get_energy_through_cross_section(polyline,
    693                                                                                kind='total')           
     729                #polyline = self.enquiry_polylines[i]
     730                #opening.total_energy = domain.get_energy_through_cross_section(polyline,
     731                #                                                               kind='total')           
     732               
     733                id = [self.enquiry_indices[i]]
     734                stage = dq['stage'].get_values(location='centroids',
     735                                               indices=id)
     736                elevation = dq['elevation'].get_values(location='centroids',
     737                                               indices=id)                                               
     738                xmomentum = dq['xmomentum'].get_values(location='centroids',
     739                                               indices=id)                                               
     740                ymomentum = dq['xmomentum'].get_values(location='centroids',
     741                                                       indices=id)                                                                                             
     742                depth = stage-elevation
     743                if depth > 0.0:
     744                    u = xmomentum/(depth + velocity_protection/depth)
     745                    v = ymomentum/(depth + velocity_protection/depth)
     746                else:
     747                    u = v = 0.0
     748
     749                   
     750                opening.total_energy = 0.5*(u*u + v*v)/g + stage
    694751                #print 'Et = %.3f m' %opening.total_energy
    695752
     
    826883
    827884
    828        
     885Culvert_flow = Culvert_flow_rating       
Note: See TracChangeset for help on using the changeset viewer.