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

Last change on this file since 6123 was 6123, checked in by ole, 15 years ago

Implented general inlet/outlet control.
Works with rating curve but not Boyd.

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