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

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

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