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

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

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