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

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

Consolidated culverts - tests pass.
Outstanding issues are

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