source: branches/numpy/anuga/culvert_flows/culvert_class.py @ 6517

Last change on this file since 6517 was 6517, checked in by rwilson, 14 years ago

Hand-merged recent changes in main trunk. Still work to be done in shallow_water.

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