source: branches/numpy/anuga/culvert_flows/culvert_class.py @ 6533

Last change on this file since 6533 was 6533, checked in by rwilson, 15 years ago

Revert back to 6481, prior to auto-merge of trunk and numpy branch.

File size: 52.8 KB
Line 
1import sys
2
3from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
4from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
5from anuga.utilities.system_tools import log_to_file
6from anuga.utilities.polygon import inside_polygon
7from anuga.utilities.polygon import is_inside_polygon
8from anuga.utilities.polygon import plot_polygons
9
10from anuga.utilities.numerical_tools import mean
11from anuga.utilities.numerical_tools import ensure_numeric, sign
12       
13from anuga.config import g, epsilon
14from anuga.config import minimum_allowed_height, velocity_protection       
15
16import numpy as num
17
18
19class Below_interval(Exception): pass 
20class Above_interval(Exception): pass
21
22# FIXME(Ole): Take a good hard look at logging here
23
24
25# FIXME(Ole): Write in C and reuse this function by similar code
26# in interpolate.py
27def interpolate_linearly(x, xvec, yvec):
28
29    msg = 'Input to function interpolate_linearly could not be converted '
30    msg += 'to numerical scalar: x = %s' % str(x)
31    try:
32        x = float(x)
33    except:
34        raise Exception, msg
35
36
37    # Check bounds
38    if x < xvec[0]: 
39        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
40            % (x, xvec[0])
41        raise Below_interval, msg
42       
43    if x > xvec[-1]: 
44        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
45            %(x, xvec[-1])
46        raise Above_interval, msg       
47       
48       
49    # Find appropriate slot within bounds           
50    i = 0
51    while x > xvec[i]: i += 1
52
53   
54    x0 = xvec[i-1]
55    x1 = xvec[i]           
56    alpha = (x - x0)/(x1 - x0) 
57           
58    y0 = yvec[i-1]
59    y1 = yvec[i]                       
60    y = alpha*y1 + (1-alpha)*y0
61   
62    return y
63           
64
65           
66def read_culvert_description(culvert_description_filename):           
67   
68    # Read description file
69    fid = open(culvert_description_filename)
70   
71    read_rating_curve_data = False       
72    rating_curve = []
73    for i, line in enumerate(fid.readlines()):
74       
75        if read_rating_curve_data is True:
76            fields = line.split(',')
77            head_difference = float(fields[0].strip())
78            flow_rate = float(fields[1].strip())               
79            barrel_velocity = float(fields[2].strip())
80           
81            rating_curve.append([head_difference, flow_rate, barrel_velocity]) 
82       
83        if i == 0:
84            # Header
85            continue
86        if i == 1:
87            # Metadata
88            fields = line.split(',')
89            label=fields[0].strip()
90            type=fields[1].strip().lower()
91            assert type in ['box', 'pipe']
92           
93            width=float(fields[2].strip())
94            height=float(fields[3].strip())
95            length=float(fields[4].strip())
96            number_of_barrels=int(fields[5].strip())
97            #fields[6] refers to losses
98            description=fields[7].strip()               
99               
100        if line.strip() == '': continue # Skip blanks
101
102        if line.startswith('Rating'):
103            read_rating_curve_data = True
104            # Flow data follows
105               
106    fid.close()
107   
108    return label, type, width, height, length, number_of_barrels, description, rating_curve
109   
110
111   
112
113class Culvert_flow_general:
114    """Culvert flow - transfer water from one hole to another
115   
116    This version will work with either rating curve file or with culvert
117    routine.
118   
119    Input: Two points, pipe_size (either diameter or width, height),
120    mannings_rougness,
121    """ 
122
123    def __init__(self,
124                 domain,
125                 culvert_description_filename=None,
126                 culvert_routine=None,
127                 end_point0=None, 
128                 end_point1=None,
129                 enquiry_point0=None, 
130                 enquiry_point1=None,
131                 type='box',
132                 width=None,
133                 height=None,
134                 length=None,
135                 number_of_barrels=1,
136                 trigger_depth=0.01, # Depth below which no flow happens
137                 manning=None,          # Mannings Roughness for Culvert
138                 sum_loss=None,
139                 use_velocity_head=False, 
140                 use_momentum_jet=False, # FIXME(Ole): Not yet implemented
141                 label=None,
142                 description=None,
143                 update_interval=None,
144                 log_file=False,
145                 discharge_hydrograph=False,
146                 verbose=False):
147       
148
149       
150        # Input check
151       
152        if height is None: height = width       
153        self.height = height
154        self.width = width
155
156       
157        assert number_of_barrels >= 1
158        assert use_velocity_head is True or use_velocity_head is False
159       
160        msg = 'Momentum jet not yet moved to general culvert'
161        assert use_momentum_jet is False, msg
162       
163        self.culvert_routine = culvert_routine       
164        self.culvert_description_filename = culvert_description_filename
165        if culvert_description_filename is not None:
166            label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
167            self.rating_curve = ensure_numeric(rating_curve)           
168       
169        self.domain = domain
170        self.trigger_depth = trigger_depth       
171               
172        if manning is None: 
173            self.manning = 0.012   # Default roughness for pipe
174       
175        if sum_loss is None:
176            self.sum_loss = 0.0
177           
178       
179                       
180        # Store culvert information
181        self.label = label
182        self.description = description
183        self.culvert_type = type
184        self.number_of_barrels = number_of_barrels
185       
186        # Store options       
187        self.use_velocity_head = use_velocity_head
188
189        if label is None: label = 'culvert_flow'
190        label += '_' + str(id(self)) 
191        self.label = label
192       
193        # File for storing discharge_hydrograph
194        if discharge_hydrograph is True:
195            self.timeseries_filename = label + '_timeseries.csv'
196            fid = open(self.timeseries_filename, 'w')
197            fid.write('time, discharge\n')
198            fid.close()
199
200        # Log file for storing general textual output
201        if log_file is True:       
202            self.log_filename = label + '.log'         
203            log_to_file(self.log_filename, self.label)       
204            log_to_file(self.log_filename, description)
205            log_to_file(self.log_filename, self.culvert_type)       
206
207
208        # Create the fundamental culvert polygons from polygon
209        P = create_culvert_polygons(end_point0,
210                                    end_point1,
211                                    width=width,   
212                                    height=height,
213                                    number_of_barrels=number_of_barrels)
214        self.culvert_polygons = P
215       
216        # Select enquiry points
217        if enquiry_point0 is None:
218            enquiry_point0 = P['enquiry_point0']
219           
220        if enquiry_point1 is None:
221            enquiry_point1 = P['enquiry_point1']           
222           
223        if verbose is True:
224            pass
225            #plot_polygons([[end_point0, end_point1],
226            #               P['exchange_polygon0'],
227            #               P['exchange_polygon1'],
228            #               [enquiry_point0, 1.005*enquiry_point0],
229            #               [enquiry_point1, 1.005*enquiry_point1]],
230            #              figname='culvert_polygon_output')
231
232           
233           
234        self.enquiry_points = [enquiry_point0, enquiry_point1]
235        self.enquiry_indices = self.get_enquiry_indices()                 
236        self.check_culvert_inside_domain()
237       
238                   
239        # Create inflow object at each end of the culvert.
240        self.openings = []
241        self.openings.append(Inflow(domain,
242                                    polygon=P['exchange_polygon0']))
243        self.openings.append(Inflow(domain,
244                                    polygon=P['exchange_polygon1']))
245                                   
246        # Assume two openings for now: Referred to as 0 and 1
247        assert len(self.openings) == 2
248
249        # Establish initial values at each enquiry point
250        dq = domain.quantities
251        for i, opening in enumerate(self.openings):
252            idx = self.enquiry_indices[i]
253            elevation = dq['elevation'].get_values(location='centroids',
254                                                   indices=[idx])[0]
255            stage = dq['stage'].get_values(location='centroids',
256                                           indices=[idx])[0]
257            opening.elevation = elevation
258            opening.stage = stage
259            opening.depth = stage-elevation
260
261           
262                           
263        # Determine initial pipe direction.
264        # This may change dynamically based on the total energy difference     
265        # Consequently, this may be superfluous
266        delta_z = self.openings[0].elevation - self.openings[1].elevation
267        if delta_z > 0.0:
268            self.inlet = self.openings[0]
269            self.outlet = self.openings[1]
270        else:               
271            self.outlet = self.openings[0]
272            self.inlet = self.openings[1]       
273       
274       
275        # Store basic geometry
276        self.end_points = [end_point0, end_point1]
277        self.vector = P['vector']
278        self.length = P['length']; assert self.length > 0.0
279        if culvert_description_filename is not None:
280            if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
281                msg = 'WARNING: barrel length specified in "%s" (%.2f m)'\
282                    % (culvert_description_filename, 
283                       length)
284                msg += ' does not match distance between specified'
285                msg += ' end points (%.2f m)' %self.length
286                print msg
287       
288        self.verbose = verbose
289
290
291       
292        # For use with update_interval                       
293        self.last_update = 0.0
294        self.update_interval = update_interval
295       
296
297        # Print some diagnostics to log if requested
298        if hasattr(self, 'log_filename'):
299            s = 'Culvert Effective Length = %.2f m' %(self.length)
300            log_to_file(self.log_filename, s)
301   
302            s = 'Culvert Direction is %s\n' %str(self.vector)
303            log_to_file(self.log_filename, s)       
304       
305       
306       
307       
308       
309    def __call__(self, domain):
310
311        # Time stuff
312        time = domain.get_time()
313       
314       
315        update = False
316        if self.update_interval is None:
317            # Use next timestep as has been computed in domain.py       
318            delta_t = domain.timestep           
319            update = True
320        else:   
321            # Use update interval
322            delta_t = self.update_interval           
323            if time - self.last_update > self.update_interval or time == 0.0:
324                update = True
325
326        if hasattr(self, 'log_filename'):           
327            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
328            log_to_file(self.log_filename, s)
329               
330                               
331        if update is True:
332            self.compute_rates(delta_t)
333           
334   
335        # Execute flow term for each opening
336        # This is where Inflow objects are evaluated using the last rate
337        # that has been calculated
338        #
339        # This will take place at every internal timestep and update the domain
340        for opening in self.openings:
341            opening(domain)
342           
343
344
345    def get_enquiry_indices(self):
346        """Get indices for nearest centroids to self.enquiry_points
347        """
348       
349        domain = self.domain
350       
351        enquiry_indices = []                 
352        for point in self.enquiry_points:
353            # Find nearest centroid
354            N = len(domain)   
355            points = domain.get_centroid_coordinates(absolute=True)
356
357            # Calculate indices in exchange area for this forcing term
358           
359            triangle_id = min_dist = sys.maxint
360            for k in range(N):
361                x, y = points[k,:] # Centroid
362
363                c = point                               
364                distance = (x-c[0])**2+(y-c[1])**2
365                if distance < min_dist:
366                    min_dist = distance
367                    triangle_id = k
368
369                   
370            if triangle_id < sys.maxint:
371                msg = 'found triangle with centroid (%f, %f)'\
372                    %tuple(points[triangle_id, :])
373                msg += ' for point (%f, %f)' %tuple(point)
374               
375                enquiry_indices.append(triangle_id)
376            else:
377                msg = 'Triangle not found for point (%f, %f)' %point
378                raise Exception, msg
379       
380        return enquiry_indices
381
382       
383    def check_culvert_inside_domain(self):
384        """Check that all polygons and enquiry points lie within the mesh.
385        """
386        bounding_polygon = self.domain.get_boundary_polygon()
387        P = self.culvert_polygons
388        for key in P.keys():
389            if key in ['exchange_polygon0', 
390                       'exchange_polygon1']:
391                for point in list(P[key]) + self.enquiry_points:
392                    msg = 'Point %s in polygon %s for culvert %s did not'\
393                        %(str(point), key, self.label)
394                    msg += 'fall within the domain boundary.'
395                    assert is_inside_polygon(point, bounding_polygon), msg
396           
397
398    def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
399        """Adjust Q downwards depending on available water at inlet
400        """
401   
402        if delta_t < epsilon:
403            # No need to adjust if time step is very small or zero
404            # In this case the possible flow will be very large
405            # anyway.
406            return Q
407       
408        # Short hands
409        domain = self.domain       
410        dq = domain.quantities               
411        time = domain.get_time()
412        I = self.inlet
413        idx = I.exchange_indices   
414
415        # Find triangle with the smallest depth
416        stage = dq['stage'].get_values(location='centroids', 
417                                               indices=[idx])
418        elevation = dq['elevation'].get_values(location='centroids', 
419                                               indices=[idx])       
420        depth = stage-elevation
421        min_depth = min(depth.flat)
422
423        # Compute possible flow for exchange region based on
424        # triangle with smallest depth
425        max_Q = min_depth*I.exchange_area/delta_t       
426       
427        # Calculate the minimum in absolute terms of
428        # the requsted flow and the possible flow
429        Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
430       
431        if abs(Q_reduced) < abs(Q): 
432            msg = '%.2fs: Requested flow is ' % time
433            msg += 'greater than what is supported by the smallest '
434            msg += 'depth at inlet exchange area:\n        '
435            msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\
436                % (min_depth, I.exchange_area, delta_t)
437            msg += ' = %.2f m^3/s\n        ' % Q_reduced
438            msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
439            if self.verbose is True:
440                print msg
441            if hasattr(self, 'log_filename'):                   
442                log_to_file(self.log_filename, msg)
443       
444        return Q_reduced   
445                       
446           
447    def compute_rates(self, delta_t):
448        """Compute new rates for inlet and outlet
449        """
450
451        # Short hands
452        domain = self.domain       
453        dq = domain.quantities               
454       
455        # Time stuff
456        time = domain.get_time()
457        self.last_update = time
458
459           
460        if hasattr(self, 'log_filename'):
461            log_filename = self.log_filename
462           
463        # Compute stage and energy at the
464        # enquiry points at each end of the culvert
465        openings = self.openings
466        for i, opening in enumerate(openings):
467            idx = self.enquiry_indices[i]               
468           
469            stage = dq['stage'].get_values(location='centroids',
470                                               indices=[idx])[0]
471            depth = h = stage-opening.elevation
472                                                           
473                                               
474            if self.use_velocity_head is True:
475                xmomentum = dq['xmomentum'].get_values(location='centroids',
476                                                       indices=[idx])[0]
477                ymomentum = dq['xmomentum'].get_values(location='centroids',
478                                                       indices=[idx])[0]
479           
480                if h > minimum_allowed_height:
481                    u = xmomentum/(h + velocity_protection/h)
482                    v = ymomentum/(h + velocity_protection/h)
483                else:
484                    u = v = 0.0
485               
486                velocity_head = 0.5*(u*u + v*v)/g   
487            else:
488                velocity_head = 0.0
489           
490            opening.total_energy = velocity_head + stage
491            opening.specific_energy = velocity_head + depth
492            opening.stage = stage
493            opening.depth = depth
494           
495
496        # We now need to deal with each opening individually
497        # Determine flow direction based on total energy difference
498        delta_total_energy = openings[0].total_energy - openings[1].total_energy
499        if delta_total_energy > 0:
500            #print 'Flow U/S ---> D/S'
501            inlet = openings[0]
502            outlet = openings[1]
503        else:
504            #print 'Flow D/S ---> U/S'
505            inlet = openings[1]
506            outlet = openings[0]
507           
508            delta_total_energy = -delta_total_energy
509
510        self.inlet = inlet
511        self.outlet = outlet
512           
513        msg = 'Total energy difference is negative'
514        assert delta_total_energy > 0.0, msg
515
516        # Recompute slope and issue warning if flow is uphill
517        # These values do not enter the computation
518        delta_z = inlet.elevation - outlet.elevation
519        culvert_slope = (delta_z/self.length)
520        if culvert_slope < 0.0:
521            # Adverse gradient - flow is running uphill
522            # Flow will be purely controlled by uphill outlet face
523            if self.verbose is True:
524                print '%.2fs - WARNING: Flow is running uphill.' % time
525           
526        if hasattr(self, 'log_filename'):
527            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
528                %(time, self.inlet.stage, self.outlet.stage)
529            log_to_file(self.log_filename, s)
530            s = 'Delta total energy = %.3f' %(delta_total_energy)
531            log_to_file(log_filename, s)
532
533           
534        # Determine controlling energy (driving head) for culvert
535        if inlet.specific_energy > delta_total_energy:
536            # Outlet control
537            driving_head = delta_total_energy
538        else:
539            # Inlet control
540            driving_head = inlet.specific_energy
541           
542
543
544        if self.inlet.depth <= self.trigger_depth:
545            Q = 0.0
546        else:
547            # Calculate discharge for one barrel and
548            # set inlet.rate and outlet.rate
549           
550            if self.culvert_description_filename is not None:
551                try:
552                    Q = interpolate_linearly(driving_head, 
553                                             self.rating_curve[:,0], 
554                                             self.rating_curve[:,1]) 
555                except Below_interval, e:
556                    Q = self.rating_curve[0,1]             
557                    msg = '%.2fs: ' % time
558                    msg += 'Delta head smaller than rating curve minimum: '
559                    msg += str(e)
560                    msg += '\n        '
561                    msg += 'I will use minimum discharge %.2f m^3/s ' % Q
562                    msg += 'for culvert "%s"' % self.label
563                   
564                    if hasattr(self, 'log_filename'):                   
565                        log_to_file(self.log_filename, msg)
566                except Above_interval, e:
567                    Q = self.rating_curve[-1,1]         
568                    msg = '%.2fs: ' % time                 
569                    msg += 'Delta head greater than rating curve maximum: '
570                    msg += str(e)
571                    msg += '\n        '
572                    msg += 'I will use maximum discharge %.2f m^3/s ' % Q
573                    msg += 'for culvert "%s"' % self.label
574                   
575                    if hasattr(self, 'log_filename'):                   
576                        log_to_file(self.log_filename, msg)
577            else:
578                # User culvert routine
579                Q, barrel_velocity, culvert_outlet_depth =\
580                    self.culvert_routine(self, delta_total_energy, g)
581           
582           
583       
584        # Adjust discharge for multiple barrels
585        Q *= self.number_of_barrels
586       
587
588        Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
589       
590        self.inlet.rate = -Q
591        self.outlet.rate = Q
592
593        # Log timeseries to file
594        try:
595            fid = open(self.timeseries_filename, 'a')       
596        except:
597            pass
598        else:   
599            fid.write('%.2f, %.2f\n' %(time, Q))
600            fid.close()
601
602                           
603   
604class Culvert_flow_rating:
605    """Culvert flow - transfer water from one hole to another
606   
607
608    Input: Two points, pipe_size (either diameter or width, height),
609    mannings_rougness,
610    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
611    top-down_blockage_factor and bottom_up_blockage_factor
612   
613    """ 
614
615    def __init__(self,
616                 domain,
617                 culvert_description_filename=None,
618                 end_point0=None, 
619                 end_point1=None,
620                 enquiry_point0=None, 
621                 enquiry_point1=None,
622                 update_interval=None,
623                 log_file=False,
624                 discharge_hydrograph=False,
625                 verbose=False):
626       
627
628       
629        label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
630       
631               
632        # Store culvert information
633        self.label = label
634        self.description = description
635        self.culvert_type = type
636        self.rating_curve = ensure_numeric(rating_curve)
637        self.number_of_barrels = number_of_barrels
638
639        if label is None: label = 'culvert_flow'
640        label += '_' + str(id(self)) 
641        self.label = label
642       
643        # File for storing discharge_hydrograph
644        if discharge_hydrograph is True:
645            self.timeseries_filename = label + '_timeseries.csv'
646            fid = open(self.timeseries_filename, 'w')
647            fid.write('time, discharge\n')
648            fid.close()
649
650        # Log file for storing general textual output
651        if log_file is True:       
652            self.log_filename = label + '.log'         
653            log_to_file(self.log_filename, self.label)       
654            log_to_file(self.log_filename, description)
655            log_to_file(self.log_filename, self.culvert_type)       
656
657
658        # Create the fundamental culvert polygons from POLYGON
659        #if self.culvert_type == 'circle':
660        #    # Redefine width and height for use with create_culvert_polygons
661        #    width = height = diameter
662       
663        P = create_culvert_polygons(end_point0,
664                                    end_point1,
665                                    width=width,   
666                                    height=height,
667                                    number_of_barrels=number_of_barrels)
668       
669        # Select enquiry points
670        if enquiry_point0 is None:
671            enquiry_point0 = P['enquiry_point0']
672           
673        if enquiry_point1 is None:
674            enquiry_point1 = P['enquiry_point1']           
675           
676        if verbose is True:
677            pass
678            #plot_polygons([[end_point0, end_point1],
679            #               P['exchange_polygon0'],
680            #               P['exchange_polygon1'],
681            #               [enquiry_point0, 1.005*enquiry_point0],
682            #               [enquiry_point1, 1.005*enquiry_point1]],                           
683            #              figname='culvert_polygon_output')
684
685           
686           
687        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
688
689        self.enquiry_indices = []                 
690        for point in self.enquiry_points:
691            # Find nearest centroid
692            N = len(domain)   
693            points = domain.get_centroid_coordinates(absolute=True)
694
695            # Calculate indices in exchange area for this forcing term
696           
697            triangle_id = min_dist = sys.maxint
698            for k in range(N):
699                x, y = points[k,:] # Centroid
700
701                c = point                               
702                distance = (x-c[0])**2+(y-c[1])**2
703                if distance < min_dist:
704                    min_dist = distance
705                    triangle_id = k
706
707                   
708            if triangle_id < sys.maxint:
709                msg = 'found triangle with centroid (%f, %f)'\
710                    %tuple(points[triangle_id, :])
711                msg += ' for point (%f, %f)' %tuple(point)
712               
713                self.enquiry_indices.append(triangle_id)
714            else:
715                msg = 'Triangle not found for point (%f, %f)' %point
716                raise Exception, msg
717       
718                         
719
720        # Check that all polygons lie within the mesh.
721        bounding_polygon = domain.get_boundary_polygon()
722        for key in P.keys():
723            if key in ['exchange_polygon0', 
724                       'exchange_polygon1']:
725                for point in list(P[key]) + self.enquiry_points:
726                    msg = 'Point %s in polygon %s for culvert %s did not'\
727                        %(str(point), key, self.label)
728                    msg += 'fall within the domain boundary.'
729                    assert is_inside_polygon(point, bounding_polygon), msg
730       
731                   
732        # Create inflow object at each end of the culvert.
733        self.openings = []
734        self.openings.append(Inflow(domain,
735                                    polygon=P['exchange_polygon0']))
736
737        self.openings.append(Inflow(domain,
738                                    polygon=P['exchange_polygon1']))                                   
739
740
741
742        dq = domain.quantities                                           
743        for i, opening in enumerate(self.openings):                           
744            elevation = dq['elevation'].get_values(location='centroids',
745                                                   indices=[self.enquiry_indices[i]])           
746            opening.elevation = elevation
747            opening.stage = elevation # Simple assumption that culvert is dry initially
748
749        # Assume two openings for now: Referred to as 0 and 1
750        assert len(self.openings) == 2
751                           
752        # Determine pipe direction     
753        self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation
754        if delta_z > 0.0:
755            self.inlet = self.openings[0]
756            self.outlet = self.openings[1]
757        else:               
758            self.outlet = self.openings[0]
759            self.inlet = self.openings[1]       
760       
761       
762        # Store basic geometry
763        self.end_points = [end_point0, end_point1]
764        self.vector = P['vector']
765        self.length = P['length']; assert self.length > 0.0
766        if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
767            msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length)
768            msg += ' does not match distance between specified'
769            msg += ' end points (%.2f m)' %self.length
770            print msg
771       
772        self.verbose = verbose
773        self.last_update = 0.0 # For use with update_interval       
774        self.last_time = 0.0               
775        self.update_interval = update_interval
776       
777
778        # Print something
779        if hasattr(self, 'log_filename'):
780            s = 'Culvert Effective Length = %.2f m' %(self.length)
781            log_to_file(self.log_filename, s)
782   
783            s = 'Culvert Direction is %s\n' %str(self.vector)
784            log_to_file(self.log_filename, s)       
785       
786       
787       
788       
789       
790    def __call__(self, domain):
791
792        # Time stuff
793        time = domain.get_time()
794       
795       
796        update = False
797        if self.update_interval is None:
798            update = True
799            delta_t = domain.timestep # Next timestep has been computed in domain.py
800        else:   
801            if time - self.last_update > self.update_interval or time == 0.0:
802                update = True
803            delta_t = self.update_interval
804           
805        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
806        if hasattr(self, 'log_filename'):           
807            log_to_file(self.log_filename, s)
808               
809                               
810        if update is True:
811            self.last_update = time
812       
813            dq = domain.quantities
814                       
815            # Get average water depths at each opening       
816            openings = self.openings   # There are two Opening [0] and [1]
817            for i, opening in enumerate(openings):
818               
819                # Compute mean values of selected quantitites in the
820                # enquiry area in front of the culvert
821               
822                stage = dq['stage'].get_values(location='centroids',
823                                               indices=[self.enquiry_indices[i]])
824               
825                # Store current average stage and depth with each opening object
826                opening.depth = stage - opening.elevation
827                opening.stage = stage
828
829               
830
831            #################  End of the FOR loop ################################################
832
833            # We now need to deal with each opening individually
834               
835            # Determine flow direction based on total energy difference
836
837            delta_w = self.inlet.stage - self.outlet.stage
838           
839            if hasattr(self, 'log_filename'):
840                s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time, 
841                                                                           self.inlet.stage,
842                                                                           self.outlet.stage)
843                log_to_file(self.log_filename, s)
844
845
846            if self.inlet.depth <= 0.01:
847                Q = 0.0
848            else:
849                # Calculate discharge for one barrel and set inlet.rate and outlet.rate
850               
851                try:
852                    Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1]) 
853                except Below_interval, e:
854                    Q = self.rating_curve[0,1]             
855                    msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
856                    msg += str(e)
857                    msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
858                        %(Q, self.label)
859                    if hasattr(self, 'log_filename'):                   
860                        log_to_file(self.log_filename, msg)
861                except Above_interval, e:
862                    Q = self.rating_curve[-1,1]         
863                    msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
864                    msg += str(e)
865                    msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
866                        %(Q, self.label)
867                    if hasattr(self, 'log_filename'):                   
868                        log_to_file(self.log_filename, msg)                   
869
870               
871               
872           
873            # Adjust discharge for multiple barrels
874            Q *= self.number_of_barrels
875           
876
877            # Adjust Q downwards depending on available water at inlet
878            stage = self.inlet.get_quantity_values(quantity_name='stage')
879            elevation = self.inlet.get_quantity_values(quantity_name='elevation')
880            depth = stage-elevation
881           
882           
883            V = 0
884            for i, d in enumerate(depth):
885                V += d * domain.areas[i]
886           
887            #Vsimple = mean(depth)*self.inlet.exchange_area # Current volume in exchange area 
888            #print 'Q', Q, 'dt', delta_t, 'Q*dt', Q*delta_t, 'V', V, 'Vsimple', Vsimple
889
890            dt = delta_t           
891            if Q*dt > V:
892           
893                Q_reduced = 0.9*V/dt # Reduce with safety factor
894               
895                msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
896                msg += 'greater than current volume (V) at inlet.\n'
897                msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
898               
899                #print msg
900               
901                if self.verbose is True:
902                    print msg
903                if hasattr(self, 'log_filename'):                   
904                    log_to_file(self.log_filename, msg)                                       
905               
906                Q = Q_reduced
907       
908            self.inlet.rate = -Q
909            self.outlet.rate = Q
910
911            # Log timeseries to file
912            try:
913                fid = open(self.timeseries_filename, 'a')       
914            except:
915                pass
916            else:   
917                fid.write('%.2f, %.2f\n' %(time, Q))
918                fid.close()
919
920            # Store value of time
921            self.last_time = time
922
923           
924   
925        # Execute flow term for each opening
926        # This is where Inflow objects are evaluated using the last rate that has been calculated
927        #
928        # This will take place at every internal timestep and update the domain
929        for opening in self.openings:
930            opening(domain)
931
932
933
934       
935       
936
937class Culvert_flow_energy:
938    """Culvert flow - transfer water from one hole to another
939   
940    Using Momentum as Calculated by Culvert Flow !!
941    Could be Several Methods Investigated to do This !!!
942
943    2008_May_08
944    To Ole:
945    OK so here we need to get the Polygon Creating code to create
946    polygons for the culvert Based on
947    the two input Points (X0,Y0) and (X1,Y1) - need to be passed
948    to create polygon
949
950    The two centers are now passed on to create_polygon.
951   
952
953    Input: Two points, pipe_size (either diameter or width, height),
954    mannings_rougness,
955    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
956    top-down_blockage_factor and bottom_up_blockage_factor
957   
958   
959    And the Delta H enquiry should be change from Openings in line 412
960    to the enquiry Polygons infront of the culvert
961    At the moment this script uses only Depth, later we can change it to
962    Energy...
963
964    Once we have Delta H can calculate a Flow Rate and from Flow Rate
965    an Outlet Velocity
966    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
967       
968    Invert levels are optional. If left out they will default to the
969    elevation at the opening.
970       
971    """ 
972
973    def __init__(self,
974                 domain,
975                 label=None, 
976                 description=None,
977                 end_point0=None, 
978                 end_point1=None,
979                 width=None,
980                 height=None,
981                 diameter=None,
982                 manning=None,          # Mannings Roughness for Culvert
983                 invert_level0=None,    # Invert level at opening 0
984                 invert_level1=None,    # Invert level at opening 1
985                 loss_exit=None,
986                 loss_entry=None,
987                 loss_bend=None,
988                 loss_special=None,
989                 blockage_topdwn=None,
990                 blockage_bottup=None,
991                 culvert_routine=None,
992                 number_of_barrels=1,
993                 enquiry_point0=None, 
994                 enquiry_point1=None,
995                 update_interval=None,
996                 verbose=False):
997       
998        # Input check
999        if diameter is not None:
1000            self.culvert_type = 'circle'
1001            self.diameter = diameter
1002            if height is not None or width is not None:
1003                msg = 'Either diameter or width&height must be specified, '
1004                msg += 'but not both.'
1005                raise Exception, msg
1006        else:
1007            if height is not None:
1008                if width is None:
1009                    self.culvert_type = 'square'                               
1010                    width = height
1011                else:
1012                    self.culvert_type = 'rectangle'
1013            elif width is not None:
1014                if height is None:
1015                    self.culvert_type = 'square'                               
1016                    height = width
1017            else:
1018                msg = 'Either diameter or width&height must be specified.'
1019                raise Exception, msg               
1020               
1021            if height == width:
1022                self.culvert_type = 'square'                                               
1023               
1024            self.height = height
1025            self.width = width
1026
1027           
1028        assert self.culvert_type in ['circle', 'square', 'rectangle']
1029       
1030        assert number_of_barrels >= 1
1031        self.number_of_barrels = number_of_barrels
1032       
1033       
1034        # Set defaults
1035        if manning is None: manning = 0.012   # Default roughness for pipe
1036        if loss_exit is None: loss_exit = 1.00
1037        if loss_entry is None: loss_entry = 0.50
1038        if loss_bend is None: loss_bend=0.00
1039        if loss_special is None: loss_special=0.00
1040        if blockage_topdwn is None: blockage_topdwn=0.00
1041        if blockage_bottup is None: blockage_bottup=0.00
1042        if culvert_routine is None: 
1043            culvert_routine=boyd_generalised_culvert_model
1044           
1045        if label is None: label = 'culvert_flow'
1046        label += '_' + str(id(self)) 
1047        self.label = label
1048       
1049        # File for storing culvert quantities
1050        self.timeseries_filename = label + '_timeseries.csv'
1051        fid = open(self.timeseries_filename, 'w')
1052        fid.write('time, E0, E1, Velocity, Discharge\n')
1053        fid.close()
1054
1055        # Log file for storing general textual output
1056        self.log_filename = label + '.log'         
1057        log_to_file(self.log_filename, self.label)       
1058        log_to_file(self.log_filename, description)
1059        log_to_file(self.log_filename, self.culvert_type)       
1060
1061
1062        # Create the fundamental culvert polygons from POLYGON
1063        if self.culvert_type == 'circle':
1064            # Redefine width and height for use with create_culvert_polygons
1065            width = height = diameter
1066       
1067        P = create_culvert_polygons(end_point0,
1068                                    end_point1,
1069                                    width=width,   
1070                                    height=height,
1071                                    number_of_barrels=number_of_barrels)
1072       
1073        # Select enquiry points
1074        if enquiry_point0 is None:
1075            enquiry_point0 = P['enquiry_point0']
1076           
1077        if enquiry_point1 is None:
1078            enquiry_point1 = P['enquiry_point1']           
1079           
1080        if verbose is True:
1081            pass
1082            #plot_polygons([[end_point0, end_point1],
1083            #               P['exchange_polygon0'],
1084            #               P['exchange_polygon1'],
1085            #               [enquiry_point0, 1.005*enquiry_point0],
1086            #               [enquiry_point1, 1.005*enquiry_point1]],
1087            #              figname='culvert_polygon_output')
1088
1089
1090        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
1091       
1092       
1093        self.enquiry_indices = []                 
1094        for point in self.enquiry_points:
1095            # Find nearest centroid
1096            N = len(domain)   
1097            points = domain.get_centroid_coordinates(absolute=True)
1098
1099            # Calculate indices in exchange area for this forcing term
1100           
1101            triangle_id = min_dist = sys.maxint
1102            for k in range(N):
1103                x, y = points[k,:] # Centroid
1104
1105                c = point                               
1106                distance = (x-c[0])**2+(y-c[1])**2
1107                if distance < min_dist:
1108                    min_dist = distance
1109                    triangle_id = k
1110
1111                   
1112            if triangle_id < sys.maxint:
1113                msg = 'found triangle with centroid (%f, %f)'\
1114                    %tuple(points[triangle_id, :])
1115                msg += ' for point (%f, %f)' %tuple(point)
1116               
1117                self.enquiry_indices.append(triangle_id)
1118            else:
1119                msg = 'Triangle not found for point (%f, %f)' %point
1120                raise Exception, msg
1121       
1122                         
1123
1124           
1125           
1126
1127        # Check that all polygons lie within the mesh.
1128        bounding_polygon = domain.get_boundary_polygon()
1129        for key in P.keys():
1130            if key in ['exchange_polygon0', 
1131                       'exchange_polygon1']:
1132                for point in P[key]:
1133               
1134                    msg = 'Point %s in polygon %s for culvert %s did not'\
1135                        %(str(point), key, self.label)
1136                    msg += 'fall within the domain boundary.'
1137                    assert is_inside_polygon(point, bounding_polygon), msg
1138       
1139
1140        # Create inflow object at each end of the culvert.
1141        self.openings = []
1142        self.openings.append(Inflow(domain,
1143                                    polygon=P['exchange_polygon0']))
1144
1145        self.openings.append(Inflow(domain,
1146                                    polygon=P['exchange_polygon1']))                                   
1147
1148
1149        # Assume two openings for now: Referred to as 0 and 1
1150        assert len(self.openings) == 2
1151       
1152        # Store basic geometry
1153        self.end_points = [end_point0, end_point1]
1154        self.invert_levels = [invert_level0, invert_level1]               
1155        #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
1156        #self.enquiry_polylines = [P['enquiry_polygon0'][:2],
1157        #                          P['enquiry_polygon1'][:2]]
1158        self.vector = P['vector']
1159        self.length = P['length']; assert self.length > 0.0
1160        self.verbose = verbose
1161        self.last_time = 0.0       
1162        self.last_update = 0.0 # For use with update_interval       
1163        self.update_interval = update_interval
1164       
1165
1166        # Store hydraulic parameters
1167        self.manning = manning
1168        self.loss_exit = loss_exit
1169        self.loss_entry = loss_entry
1170        self.loss_bend = loss_bend
1171        self.loss_special = loss_special
1172        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
1173        self.blockage_topdwn = blockage_topdwn
1174        self.blockage_bottup = blockage_bottup
1175
1176        # Store culvert routine
1177        self.culvert_routine = culvert_routine
1178
1179       
1180        # Create objects to update momentum (a bit crude at this stage)
1181
1182       
1183        xmom0 = General_forcing(domain, 'xmomentum',
1184                                polygon=P['exchange_polygon0'])
1185
1186        xmom1 = General_forcing(domain, 'xmomentum',
1187                                polygon=P['exchange_polygon1'])
1188
1189        ymom0 = General_forcing(domain, 'ymomentum',
1190                                polygon=P['exchange_polygon0'])
1191
1192        ymom1 = General_forcing(domain, 'ymomentum',
1193                                polygon=P['exchange_polygon1'])
1194
1195        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
1196       
1197
1198        # Print something
1199        s = 'Culvert Effective Length = %.2f m' %(self.length)
1200        log_to_file(self.log_filename, s)
1201   
1202        s = 'Culvert Direction is %s\n' %str(self.vector)
1203        log_to_file(self.log_filename, s)       
1204       
1205       
1206    def __call__(self, domain):
1207
1208        log_filename = self.log_filename
1209         
1210        # Time stuff
1211        time = domain.get_time()
1212       
1213        # Short hand
1214        dq = domain.quantities
1215               
1216
1217        update = False
1218        if self.update_interval is None:
1219            update = True
1220            delta_t = domain.timestep # Next timestep has been computed in domain.py
1221        else:   
1222            if time - self.last_update > self.update_interval or time == 0.0:
1223                update = True
1224            delta_t = self.update_interval
1225           
1226        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
1227        if hasattr(self, 'log_filename'):           
1228            log_to_file(log_filename, s)
1229               
1230                               
1231        if update is True:
1232            self.last_update = time
1233                       
1234            msg = 'Time did not advance'
1235            if time > 0.0: assert delta_t > 0.0, msg
1236
1237
1238            # Get average water depths at each opening       
1239            openings = self.openings   # There are two Opening [0] and [1]
1240            for i, opening in enumerate(openings):
1241               
1242                # Compute mean values of selected quantitites in the
1243                # exchange area in front of the culvert
1244                     
1245                stage = opening.get_quantity_values(quantity_name='stage')
1246                w = mean(stage) # Average stage
1247
1248                # Use invert level instead of elevation if specified
1249                invert_level = self.invert_levels[i]
1250                if invert_level is not None:
1251                    z = invert_level
1252                else:
1253                    elevation = opening.get_quantity_values(quantity_name='elevation')
1254                    z = mean(elevation) # Average elevation
1255
1256                # Estimated depth above the culvert inlet
1257                d = w - z  # Used for calculations involving depth
1258                if d < 0.0:
1259                    # This is possible since w and z are taken at different locations
1260                    #msg = 'D < 0.0: %f' %d
1261                    #raise msg
1262                    d = 0.0
1263               
1264
1265                # Ratio of depth to culvert height.
1266                # If ratio > 1 then culvert is running full
1267                if self.culvert_type == 'circle':
1268                    ratio = d/self.diameter
1269                else:   
1270                    ratio = d/self.height 
1271                opening.ratio = ratio
1272                   
1273                   
1274                # Average measures of energy in front of this opening
1275               
1276                id = [self.enquiry_indices[i]]
1277                stage = dq['stage'].get_values(location='centroids',
1278                                               indices=id)
1279                elevation = dq['elevation'].get_values(location='centroids',
1280                                                       indices=id)                                               
1281                xmomentum = dq['xmomentum'].get_values(location='centroids',
1282                                                       indices=id)                                               
1283                ymomentum = dq['xmomentum'].get_values(location='centroids',
1284                                                       indices=id)                                                                                             
1285                depth = stage-elevation
1286                if depth > 0.0:
1287                    u = xmomentum/(depth + velocity_protection/depth)
1288                    v = ymomentum/(depth + velocity_protection/depth)
1289                else:
1290                    u = v = 0.0
1291
1292                   
1293                opening.total_energy = 0.5*(u*u + v*v)/g + stage
1294                #print 'Et = %.3f m' %opening.total_energy
1295
1296                # Store current average stage and depth with each opening object
1297                opening.depth = d
1298                opening.depth_trigger = d # Use this for now
1299                opening.stage = w
1300                opening.elevation = z
1301               
1302
1303            #################  End of the FOR loop ################################################
1304
1305            # We now need to deal with each opening individually
1306               
1307            # Determine flow direction based on total energy difference
1308            delta_Et = openings[0].total_energy - openings[1].total_energy
1309
1310            if delta_Et > 0:
1311                #print 'Flow U/S ---> D/S'
1312                inlet = openings[0]
1313                outlet = openings[1]
1314
1315                inlet.momentum = self.opening_momentum[0]
1316                outlet.momentum = self.opening_momentum[1]
1317
1318            else:
1319                #print 'Flow D/S ---> U/S'
1320                inlet = openings[1]
1321                outlet = openings[0]
1322
1323                inlet.momentum = self.opening_momentum[1]
1324                outlet.momentum = self.opening_momentum[0]
1325               
1326                delta_Et = -delta_Et
1327
1328            self.inlet = inlet
1329            self.outlet = outlet
1330               
1331            msg = 'Total energy difference is negative'
1332            assert delta_Et > 0.0, msg
1333
1334            delta_h = inlet.stage - outlet.stage
1335            delta_z = inlet.elevation - outlet.elevation
1336            culvert_slope = (delta_z/self.length)
1337
1338            if culvert_slope < 0.0:
1339                # Adverse gradient - flow is running uphill
1340                # Flow will be purely controlled by uphill outlet face
1341                if self.verbose is True:
1342                    print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
1343
1344
1345            s = 'Delta total energy = %.3f' %(delta_Et)
1346            log_to_file(log_filename, s)
1347
1348
1349            # Calculate discharge for one barrel and set inlet.rate and outlet.rate
1350            Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
1351       
1352            # Adjust discharge for multiple barrels
1353            Q *= self.number_of_barrels
1354
1355            # Compute barrel momentum
1356            barrel_momentum = barrel_velocity*culvert_outlet_depth
1357                   
1358            s = 'Barrel velocity = %f' %barrel_velocity
1359            log_to_file(log_filename, s)
1360
1361            # Compute momentum vector at outlet
1362            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
1363               
1364            s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
1365            log_to_file(log_filename, s)
1366
1367            # Log timeseries to file
1368            fid = open(self.timeseries_filename, 'a')       
1369            fid.write('%f, %f, %f, %f, %f\n'\
1370                          %(time, 
1371                            openings[0].total_energy,
1372                            openings[1].total_energy,
1373                            barrel_velocity,
1374                            Q))
1375            fid.close()
1376
1377            # Update momentum       
1378
1379            if delta_t > 0.0:
1380                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
1381                xmomentum_rate /= delta_t
1382                   
1383                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
1384                ymomentum_rate /= delta_t
1385                       
1386                s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
1387                log_to_file(log_filename, s)                   
1388            else:
1389                xmomentum_rate = ymomentum_rate = 0.0
1390
1391
1392            # Set momentum rates for outlet jet
1393            outlet.momentum[0].rate = xmomentum_rate
1394            outlet.momentum[1].rate = ymomentum_rate
1395
1396            # Remember this value for next step (IMPORTANT)
1397            outlet.momentum[0].value = outlet_mom_x
1398            outlet.momentum[1].value = outlet_mom_y                   
1399
1400            if int(domain.time*100) % 100 == 0:
1401                s = 'T=%.5f, Culvert Discharge = %.3f f'\
1402                    %(time, Q)
1403                s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
1404                     %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
1405                s += ' Momentum rate: (%.4f, %.4f)'\
1406                     %(xmomentum_rate, ymomentum_rate)                   
1407                s+='Outlet Vel= %.3f'\
1408                    %(barrel_velocity)
1409                log_to_file(log_filename, s)
1410           
1411            # Store value of time
1412            self.last_time = time
1413               
1414
1415
1416        # Execute flow term for each opening
1417        # This is where Inflow objects are evaluated and update the domain
1418        for opening in self.openings:
1419            opening(domain)
1420
1421        # Execute momentum terms
1422        # This is where Inflow objects are evaluated and update the domain
1423        self.outlet.momentum[0](domain)
1424        self.outlet.momentum[1](domain)       
1425           
1426
1427
1428Culvert_flow = Culvert_flow_general       
Note: See TracBrowser for help on using the repository browser.