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

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

Developed protection against situation where requested culvert flows exceed
what is possible at the inlet. This works by computing the largest possible flow for each triangle and use that to reduce the overall flow if necessary.

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