source: trunk/anuga_core/source/anuga/structures/culvert_operator.py @ 7962

Last change on this file since 7962 was 7962, checked in by steve, 13 years ago

Got a flow from one region to another.
File size: 28.6 KB
Line 
1import sys
2
3from anuga.shallow_water.forcing import Inflow, General_forcing
4from anuga.utilities.system_tools import log_to_file
5from anuga.geometry.polygon import inside_polygon, is_inside_polygon
6from anuga.geometry.polygon import plot_polygons, polygon_area
7
8
9from anuga.utilities.numerical_tools import mean
10from anuga.utilities.numerical_tools import ensure_numeric, sign
11       
12from anuga.config import g, epsilon
13from anuga.config import minimum_allowed_height, velocity_protection       
14import anuga.utilities.log as log
15
16import numpy as num
17from math import sqrt
18
19class Below_interval(Exception): pass 
20class Above_interval(Exception): pass
21
22   
23
24class Generic_box_culvert:
25    """Culvert flow - transfer water from one rectangular box to another.
26    Sets up the geometry of problem
27   
28    This is the base class for culverts. Inherit from this class (and overwrite
29    compute_discharge method for specific subclasses)
30   
31    Input: Two points, pipe_size (either diameter or width, height),
32    mannings_rougness,
33    """ 
34
35    def __init__(self,
36                 domain,
37                 end_point0=None, 
38                 end_point1=None,
39                 enquiry_gap_factor=0.2,
40                 width=None,
41                 height=None,
42                 verbose=False):
43       
44        # Input check
45       
46        self.domain = domain
47
48        self.domain.set_fractional_step_operator(self)
49
50        self.end_points= [end_point0, end_point1]
51        self.enquiry_gap_factor = enquiry_gap_factor
52       
53        if height is None:
54            height = width
55
56        self.width = width
57        self.height = height
58       
59        self.verbose=verbose
60        self.filename = None
61       
62        # Create the fundamental culvert polygons from polygon
63        self.create_culvert_polygons()
64        self.compute_enquiry_indices()
65        self.check_culvert_inside_domain()
66        self.compute_exchange_triangle_indices()
67
68
69        #self.print_stats()
70
71
72
73    def __call__(self):
74
75
76        # Time stuff
77        time     = self.domain.get_time()
78        timestep = self.domain.get_timestep()
79
80
81        inlet_indices  = self.exchange_triangle_indices[0]
82        outlet_indices = self.exchange_triangle_indices[1]
83
84        areas      = self.domain.areas
85        stage     = self.domain.quantities['stage'].centroid_values
86        elevation = self.domain.quantities['elevation'].centroid_values
87       
88        xmom = self.domain.quantities['xmomentum'].centroid_values
89        ymom = self.domain.quantities['ymomentum'].centroid_values
90
91        # Inlet averages
92        inlet_heights  = stage[inlet_indices]-elevation[inlet_indices]
93        inlet_areas = areas[inlet_indices]
94
95        inlet_water = num.sum(inlet_heights*inlet_areas)
96
97        average_inlet_water = inlet_water/self.exchange_areas[0]
98
99        # Outlet averages
100        outlet_heights  = stage[outlet_indices]-elevation[outlet_indices]
101        outlet_areas = areas[outlet_indices]
102
103        outlet_water = num.sum(outlet_heights*outlet_areas)
104
105        average_outlet_water = outlet_water/self.exchange_areas[1]
106
107
108        # Transfer
109        transfer_water = timestep*inlet_water
110
111        stage[inlet_indices] = elevation[inlet_indices] + average_inlet_water - transfer_water
112        xmom[inlet_indices] = 0.0
113        ymom[inlet_indices] = 0.0
114
115
116        stage[outlet_indices] = elevation[outlet_indices] + average_outlet_water + transfer_water
117        xmom[outlet_indices] = 0.0
118        ymom[outlet_indices] = 0.0
119
120
121    def print_stats(self):
122
123        print '====================================='
124        print 'Generic Culvert Operator'
125        print '====================================='
126        print "enquiry_gap_factor"
127        print self.enquiry_gap_factor
128       
129        for i in [0,1]:
130            print '-------------------------------------'
131            print 'exchange_region %i' % i
132            print '-------------------------------------'
133
134            print 'exchange triangle indices and centres'
135            print self.exchange_triangle_indices[i]
136            print self.domain.get_centroid_coordinates()[self.exchange_triangle_indices[i]]
137       
138            print 'end_point'
139            print self.end_points[i]
140
141
142            print 'exchange_polygon'
143            print self.exchange_polygons[i]
144
145            print 'enquiry_point'
146            print self.enquiry_points[i]
147
148        print '====================================='
149
150
151 
152
153
154    def compute_exchange_triangle_indices(self):
155
156        # Get boundary (in absolute coordinates)
157        domain = self.domain
158        bounding_polygon = domain.get_boundary_polygon()
159        centroids = domain.get_centroid_coordinates(absolute=True)
160        self.exchange_triangle_indices = []
161        self.exchange_areas = []
162
163        for i in [0,1]:
164            exchange_polygon = self.exchange_polygons[i]
165
166            # Check that polygon lies within the mesh.
167            for point in exchange_polygon:
168                msg = 'Point %s in polygon for forcing term' % str(point)
169                msg += ' did not fall within the domain boundary.'
170                assert is_inside_polygon(point, bounding_polygon), msg
171
172           
173            exchange_triangle_indices = inside_polygon(centroids, exchange_polygon)
174
175            if len(exchange_triangle_indices) == 0:
176                region = 'polygon=%s' % (exchange_polygon)
177                msg = 'No triangles have been identified in '
178                msg += 'specified region: %s' % region
179                raise Exception, msg
180
181            # Compute exchange area as the sum of areas of triangles identified
182            # by polygon
183            exchange_area = 0.0
184            for j in exchange_triangle_indices:
185                exchange_area += domain.areas[j]
186
187
188            msg = 'Exchange area %f in culvert' % i
189            msg += ' has area = %f' % exchange_area
190            assert exchange_area > 0.0
191
192            self.exchange_triangle_indices.append(exchange_triangle_indices)
193            self.exchange_areas.append(exchange_area)
194
195
196
197    def set_store_hydrograph_discharge(self, filename=None):
198
199        if filename is None:
200            self.filename = 'culvert_discharge_hydrograph'
201        else:
202            self.filename = filename
203
204        self.discharge_hydrograph = True
205       
206        self.timeseries_filename = self.filename + '_timeseries.csv'
207        fid = open(self.timeseries_filename, 'w')
208        fid.write('time, discharge\n')
209        fid.close()
210
211    def create_culvert_polygons(self):
212        """Create polygons at the end of a culvert inlet and outlet.
213        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
214        for enquiring the total energy at both ends of the culvert and transferring flow.
215        """
216
217        # Calculate geometry
218        x0, y0 = self.end_points[0]
219        x1, y1 = self.end_points[1]
220
221        dx = x1 - x0
222        dy = y1 - y0
223
224        self.culvert_vector = num.array([dx, dy])
225        self.culvert_length = sqrt(num.sum(self.culvert_vector**2))
226        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
227
228        # Unit direction vector and normal
229        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
230        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
231
232        # Short hands
233        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
234        h = self.height*self.culvert_vector    # Vector of length=height in the
235                             # direction of the culvert
236        gap = (1 + self.enquiry_gap_factor)*h
237
238        self.exchange_polygons = []
239        self.enquiry_points = []
240
241        # Build exchange polygon and enquiry points 0 and 1
242        for i in [0, 1]:
243            i0 = (2*i-1)
244            p0 = self.end_points[i] + w
245            p1 = self.end_points[i] - w
246            p2 = p1 + i0*h
247            p3 = p0 + i0*h
248            self.exchange_polygons.append(num.array([p0, p1, p2, p3]))
249            self.enquiry_points.append(self.end_points[i] + i0*gap)
250
251
252
253
254        # Check that enquiry points are outside exchange polygons
255        for i in [0,1]:
256            polygon = self.exchange_polygons[i]
257            # FIXME (SR) Probably should calculate the area of all the triangles
258            # associated with this polygon, as there is likely to be some
259            # inconsistency between triangles and ploygon
260            area = polygon_area(polygon)
261           
262
263            msg = 'Polygon %s ' %(polygon)
264            msg += ' has area = %f' % area
265            assert area > 0.0, msg
266
267            for j in [0,1]:
268                point = self.enquiry_points[j]
269                msg = 'Enquiry point falls inside a culvert polygon.'
270
271                assert not inside_polygon(point, polygon), msg
272
273       
274
275
276    def compute_enquiry_indices(self):
277        """Get indices for nearest centroids to self.enquiry_points
278        """
279       
280        domain = self.domain
281       
282        enquiry_indices = []                 
283        for point in self.enquiry_points:
284            # Find nearest centroid
285            N = len(domain)   
286            points = domain.get_centroid_coordinates(absolute=True)
287
288            # Calculate indices in exchange area for this forcing term
289           
290            triangle_id = min_dist = sys.maxint
291            for k in range(N):
292                x, y = points[k,:] # Centroid
293
294                c = point                               
295                distance = (x-c[0])**2+(y-c[1])**2
296                if distance < min_dist:
297                    min_dist = distance
298                    triangle_id = k
299
300                   
301            if triangle_id < sys.maxint:
302                msg = 'found triangle with centroid (%f, %f)'\
303                    %tuple(points[triangle_id, :])
304                msg += ' for point (%f, %f)' %tuple(point)
305               
306                enquiry_indices.append(triangle_id)
307            else:
308                msg = 'Triangle not found for point (%f, %f)' %point
309                raise Exception, msg
310       
311        self.enquiry_indices = enquiry_indices
312
313       
314    def check_culvert_inside_domain(self):
315        """Check that all polygons and enquiry points lie within the mesh.
316        """
317        bounding_polygon = self.domain.get_boundary_polygon()
318        for i in [0, 1]:
319            for point in list(self.exchange_polygons[i]) + self.enquiry_points:
320                msg = 'Point %s did not '\
321                    %(str(point))
322                msg += 'fall within the domain boundary.'
323                assert is_inside_polygon(point, bounding_polygon), msg
324           
325
326    def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
327        """Adjust Q downwards depending on available water at inlet
328
329           This is a critical step in modelling bridges and Culverts
330           the predicted flow through a structure based on an abstract
331           algorithm can at times request for water that is simply not
332           available due to any number of constrictions that limit the
333           flow approaching the structure In order to ensure that
334           there is adequate flow available certain checks are
335           required There needs to be a check using the Static Water
336           Volume sitting infront of the structure, In addition if the
337           water is moving the available water will be larger than the
338           static volume
339           
340           NOTE To temporarily switch this off for Debugging purposes
341           rem out line in function def compute_rates below
342        """
343   
344        if delta_t < epsilon:
345            # No need to adjust if time step is very small or zero
346            # In this case the possible flow will be very large
347            # anyway.
348            return Q
349       
350        # Short hands
351        domain = self.domain       
352        dq = domain.quantities               
353        time = domain.get_time()
354        I = self.inlet
355        idx = I.exchange_indices   
356
357        # Find triangle with the smallest depth
358        stage = dq['stage'].get_values(location='centroids', 
359                                               indices=[idx])
360        elevation = dq['elevation'].get_values(location='centroids', 
361                                               indices=[idx])       
362        depth = stage-elevation
363        min_depth = min(depth.flat)  # This may lead to errors if edge of area is at a higher level !!!!
364        avg_depth = mean(depth.flat) # Yes, but this one violates the conservation unit tests
365
366
367
368        # FIXME (Ole): If you want these, use log.critical() and
369        # make the statements depend on verbose
370        #print I.depth
371        #print I.velocity
372        #print self.width
373
374        # max_Q Based on Volume Calcs
375
376
377        depth_term = min_depth*I.exchange_area/delta_t
378        if min_depth < 0.2:
379            # Only add velocity term in shallow waters (< 20 cm)
380            # This is a little ad hoc, but maybe it is reasonable
381            velocity_term = self.width*min_depth*I.velocity
382        else:
383            velocity_term = 0.0
384
385        # This one takes approaching water into account   
386        max_Q = max(velocity_term, depth_term)
387
388        # This one preserves Volume
389        #max_Q = depth_term
390
391
392        if self.verbose is True:
393            log.critical('Max_Q = %f' % max_Q)           
394            msg = 'Width = %.2fm, Depth at inlet = %.2f m, Velocity = %.2f m/s.      ' % (self.width, I.depth, I.velocity)
395            msg += 'Max Q = %.2f m^3/s' %(max_Q)
396            log.critical(msg)
397
398        if self.log_filename is not None:               
399            log_to_file(self.log_filename, msg)
400        # New Procedure for assessing the flow available to the Culvert
401        # This routine uses the GET FLOW THROUGH CROSS SECTION
402        #   Need to check Several Polyline however
403        # Firstly 3 sides of the exchange Poly
404        # then only the Line Directly infront of the Polygon
405        # Access polygon Points from   self.inlet.polygon
406     
407        #  The Following computes the flow crossing over 3 sides of the exchange polygon for the structure
408        # Clearly the flow in the culvert can not be more than that flowing toward it through the exhange polygon
409       
410        #q1 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][1:3]) # First Side Segment
411        #q2 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][2:])   # Second Face Segment
412        #q3 =domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'].take([3,0], axis=0)) # Third Side Segment
413        # q4 = domain.get_flow_through_cross_section([self.culvert_polygons['exchange_polygon0'][1:4]][0])
414        #q4=max(q1,0.0)+max(q2,0.0)+max(q3,0.0)
415        # To use only the Flow crossing the 3 sides of the Exchange Polygon use the following Line Only
416        #max_Q=max(q1,q2,q3,q4)
417        # Try Simple Smoothing using Average of 2 approaches
418        #max_Q=(max(q1,q2,q3,q4)+max_Q)/2.0
419        # Calculate the minimum in absolute terms of
420        # the requsted flow and the possible flow
421        Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
422        if self.verbose is True:
423            msg = 'Initial Q Reduced = %.2f m3/s.      ' % (Q_reduced)
424            log.critical(msg)
425
426        if self.log_filename is not None:               
427            log_to_file(self.log_filename, msg)
428        # Now Keep Rolling Average of Computed Discharge to Reduce / Remove Oscillations
429        #  can use delta_t if we want to averageover a time frame for example
430        # N = 5.0/delta_t  Will provide the average over 5 seconds
431
432        self.i=(self.i+1)%self.N
433        self.Q_list[self.i]=Q_reduced
434        Q_reduced = sum(self.Q_list)/len(self.Q_list)
435
436        if self.verbose is True:
437            msg = 'Final Q Reduced = %.2f m3/s.      ' % (Q_reduced)
438            log.critical(msg)
439
440        if self.log_filename is not None:               
441            log_to_file(self.log_filename, msg)
442
443
444        if abs(Q_reduced) < abs(Q): 
445            msg = '%.2fs: Requested flow is ' % time
446            msg += 'greater than what is supported by the smallest '
447            msg += 'depth at inlet exchange area:\n        '
448            msg += 'inlet exchange area: %.2f '% (I.exchange_area) 
449            msg += 'velocity at inlet :%.2f '% (I.velocity)
450            msg += 'Vel* Exch Area = : %.2f '% (I.velocity*avg_depth*self.width)
451            msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\
452                % (avg_depth, I.exchange_area, delta_t)
453            msg += ' = %.2f m^3/s\n        ' % Q_reduced
454            msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
455            msg += 'Note calculate max_Q from V %.2f m^3/s ' % (max_Q)
456            if self.verbose is True:
457                log.critical(msg)
458               
459            if self.log_filename is not None:               
460                log_to_file(self.log_filename, msg)
461       
462        return Q_reduced   
463                       
464           
465    def compute_rates(self, delta_t):
466        """Compute new rates for inlet and outlet
467        """
468
469        # Short hands
470        domain = self.domain       
471        dq = domain.quantities               
472       
473        # Time stuff
474        time = domain.get_time()
475        self.last_update = time
476
477           
478        if hasattr(self, 'log_filename'):
479            log_filename = self.log_filename
480           
481        # Compute stage, energy and velocity at the
482        # enquiry points at each end of the culvert
483        openings = self.openings
484        for i, opening in enumerate(openings):
485            idx = self.enquiry_indices[i]               
486           
487            stage = dq['stage'].get_values(location='centroids',
488                                           indices=[idx])[0]
489            depth = h = stage-opening.elevation
490                                                           
491           
492            # Get velocity                                 
493            xmomentum = dq['xmomentum'].get_values(location='centroids',
494                                                   indices=[idx])[0]
495            ymomentum = dq['xmomentum'].get_values(location='centroids',
496                                                   indices=[idx])[0]
497
498            if h > minimum_allowed_height:
499                u = xmomentum/(h + velocity_protection/h)
500                v = ymomentum/(h + velocity_protection/h)
501            else:
502                u = v = 0.0
503               
504            v_squared = u*u + v*v
505           
506            if self.use_velocity_head is True:
507                velocity_head = 0.5*v_squared/g   
508            else:
509                velocity_head = 0.0
510           
511            opening.total_energy = velocity_head + stage
512            opening.specific_energy = velocity_head + depth
513            opening.stage = stage
514            opening.depth = depth
515            opening.velocity = sqrt(v_squared)
516           
517
518        # We now need to deal with each opening individually
519        # Determine flow direction based on total energy difference
520        delta_total_energy = openings[0].total_energy - openings[1].total_energy
521        if delta_total_energy > 0:
522            inlet = openings[0]
523            outlet = openings[1]
524
525            # FIXME: I think this whole momentum jet thing could be a bit more elegant
526            inlet.momentum = self.opening_momentum[0]
527            outlet.momentum = self.opening_momentum[1]
528        else:
529            inlet = openings[1]
530            outlet = openings[0]
531           
532            inlet.momentum = self.opening_momentum[1]
533            outlet.momentum = self.opening_momentum[0]
534
535            delta_total_energy = -delta_total_energy
536
537        self.inlet = inlet
538        self.outlet = outlet
539           
540        msg = 'Total energy difference is negative'
541        assert delta_total_energy >= 0.0, msg
542
543        # Recompute slope and issue warning if flow is uphill
544        # These values do not enter the computation
545        delta_z = inlet.elevation - outlet.elevation
546        culvert_slope = (delta_z/self.length)
547        if culvert_slope < 0.0:
548            # Adverse gradient - flow is running uphill
549            # Flow will be purely controlled by uphill outlet face
550            if self.verbose is True:
551                log.critical('%.2fs - WARNING: Flow is running uphill.' % time)
552           
553        if self.log_filename is not None:
554            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
555                %(time, self.inlet.stage, self.outlet.stage)
556            log_to_file(self.log_filename, s)
557            s = 'Delta total energy = %.3f' %(delta_total_energy)
558            log_to_file(log_filename, s)
559
560           
561        # Determine controlling energy (driving head) for culvert
562        if inlet.specific_energy > delta_total_energy:
563            # Outlet control
564            driving_head = delta_total_energy
565        else:
566            # Inlet control
567            driving_head = inlet.specific_energy
568           
569
570
571        if self.inlet.depth <= self.trigger_depth:
572            Q = 0.0
573        else:
574            # Calculate discharge for one barrel and
575            # set inlet.rate and outlet.rate
576           
577            if self.culvert_description_filename is not None:
578                try:
579                    Q = interpolate_linearly(driving_head, 
580                                             self.rating_curve[:,0], 
581                                             self.rating_curve[:,1]) 
582                except Below_interval, e:
583                    Q = self.rating_curve[0,1]             
584                    msg = '%.2fs: ' % time
585                    msg += 'Delta head smaller than rating curve minimum: '
586                    msg += str(e)
587                    msg += '\n        '
588                    msg += 'I will use minimum discharge %.2f m^3/s ' % Q
589                    msg += 'for culvert "%s"' % self.label
590                   
591                    if hasattr(self, 'log_filename'):                   
592                        log_to_file(self.log_filename, msg)
593                except Above_interval, e:
594                    Q = self.rating_curve[-1,1]         
595                    msg = '%.2fs: ' % time                 
596                    msg += 'Delta head greater than rating curve maximum: '
597                    msg += str(e)
598                    msg += '\n        '
599                    msg += 'I will use maximum discharge %.2f m^3/s ' % Q
600                    msg += 'for culvert "%s"' % self.label
601                   
602                    if self.log_filename is not None:                   
603                        log_to_file(self.log_filename, msg)
604            else:
605                # User culvert routine
606                Q, barrel_velocity, culvert_outlet_depth =\
607                    self.culvert_routine(inlet.depth,
608                                         outlet.depth,
609                                         inlet.velocity,
610                                         outlet.velocity,
611                                         inlet.specific_energy, 
612                                         delta_total_energy, 
613                                         g,
614                                         culvert_length=self.length,
615                                         culvert_width=self.width,
616                                         culvert_height=self.height,
617                                         culvert_type=self.culvert_type,
618                                         manning=self.manning,
619                                         sum_loss=self.sum_loss,
620                                         log_filename=self.log_filename)
621           
622           
623       
624        # Adjust discharge for multiple barrels
625        Q *= self.number_of_barrels
626
627        # Adjust discharge for available water at the inlet
628        Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
629       
630        self.inlet.rate = -Q
631        self.outlet.rate = Q
632
633
634        # Momentum jet stuff
635        if self.use_momentum_jet is True:
636
637
638            # Compute barrel momentum
639            barrel_momentum = barrel_velocity*culvert_outlet_depth
640
641            if self.log_filename is not None:                                   
642                s = 'Barrel velocity = %f' %barrel_velocity
643                log_to_file(self.log_filename, s)
644
645            # Compute momentum vector at outlet
646            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
647               
648            if self.log_filename is not None:               
649                s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
650                log_to_file(self.log_filename, s)
651
652
653            # Update momentum       
654            if delta_t > 0.0:
655                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
656                xmomentum_rate /= delta_t
657                   
658                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
659                ymomentum_rate /= delta_t
660                       
661                if self.log_filename is not None:               
662                    s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
663                    log_to_file(self.log_filename, s)                   
664            else:
665                xmomentum_rate = ymomentum_rate = 0.0
666
667
668            # Set momentum rates for outlet jet
669            outlet.momentum[0].rate = xmomentum_rate
670            outlet.momentum[1].rate = ymomentum_rate
671
672            # Remember this value for next step (IMPORTANT)
673            outlet.momentum[0].value = outlet_mom_x
674            outlet.momentum[1].value = outlet_mom_y                   
675
676            if int(domain.time*100) % 100 == 0:
677
678                if self.log_filename is not None:               
679                    s = 'T=%.5f, Culvert Discharge = %.3f f'\
680                        %(time, Q)
681                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
682                        %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
683                    s += ' Momentum rate: (%.4f, %.4f)'\
684                        %(xmomentum_rate, ymomentum_rate)                   
685                    s+='Outlet Vel= %.3f'\
686                        %(barrel_velocity)
687                    log_to_file(self.log_filename, s)
688
689
690            # Execute momentum terms
691            # This is where Inflow objects are evaluated and update the domain
692                self.outlet.momentum[0](domain)
693                self.outlet.momentum[1](domain)       
694           
695
696           
697        # Log timeseries to file
698        try:
699            fid = open(self.timeseries_filename, 'a')       
700        except:
701            pass
702        else:   
703            fid.write('%.2f, %.2f\n' %(time, Q))
704            fid.close()
705
706        # Store value of time
707        self.last_time = time
708
709
710# FIXME(Ole): Write in C and reuse this function by similar code
711# in interpolate.py
712def interpolate_linearly(x, xvec, yvec):
713
714    msg = 'Input to function interpolate_linearly could not be converted '
715    msg += 'to numerical scalar: x = %s' % str(x)
716    try:
717        x = float(x)
718    except:
719        raise Exception, msg
720
721
722    # Check bounds
723    if x < xvec[0]:
724        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
725            % (x, xvec[0])
726        raise Below_interval, msg
727
728    if x > xvec[-1]:
729        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
730            %(x, xvec[-1])
731        raise Above_interval, msg
732
733
734    # Find appropriate slot within bounds
735    i = 0
736    while x > xvec[i]: i += 1
737
738
739    x0 = xvec[i-1]
740    x1 = xvec[i]
741    alpha = (x - x0)/(x1 - x0)
742
743    y0 = yvec[i-1]
744    y1 = yvec[i]
745    y = alpha*y1 + (1-alpha)*y0
746
747    return y
748
749
750
751def read_culvert_description(culvert_description_filename):
752
753    # Read description file
754    fid = open(culvert_description_filename)
755
756    read_rating_curve_data = False
757    rating_curve = []
758    for i, line in enumerate(fid.readlines()):
759
760        if read_rating_curve_data is True:
761            fields = line.split(',')
762            head_difference = float(fields[0].strip())
763            flow_rate = float(fields[1].strip())
764            barrel_velocity = float(fields[2].strip())
765
766            rating_curve.append([head_difference, flow_rate, barrel_velocity])
767
768        if i == 0:
769            # Header
770            continue
771        if i == 1:
772            # Metadata
773            fields = line.split(',')
774            label=fields[0].strip()
775            type=fields[1].strip().lower()
776            assert type in ['box', 'pipe']
777
778            width=float(fields[2].strip())
779            height=float(fields[3].strip())
780            length=float(fields[4].strip())
781            number_of_barrels=int(fields[5].strip())
782            #fields[6] refers to losses
783            description=fields[7].strip()
784
785        if line.strip() == '': continue # Skip blanks
786
787        if line.startswith('Rating'):
788            read_rating_curve_data = True
789            # Flow data follows
790
791    fid.close()
792
793    return label, type, width, height, length, number_of_barrels, description, rating_curve
794
Note: See TracBrowser for help on using the repository browser.