source: anuga_core/source/anuga/culvert_flows/culvert_class.py @ 7340

Last change on this file since 7340 was 7323, checked in by ole, 16 years ago

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