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

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

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