source: trunk/anuga_core/source/anuga/structures/culvert_class.py @ 7939

Last change on this file since 7939 was 7939, checked in by steve, 13 years ago

Adding in culvert routines into structures directory
File size: 61.6 KB
RevLine 
[7939]1import sys
2
3from anuga.shallow_water.forcing 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.geometry.polygon import inside_polygon, is_inside_polygon, plot_polygons
7
8from anuga.utilities.numerical_tools import mean
9from anuga.utilities.numerical_tools import ensure_numeric, sign
10       
11from anuga.config import g, epsilon
12from anuga.config import minimum_allowed_height, velocity_protection       
13import anuga.utilities.log as log
14
15import numpy as num
16from math import sqrt
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                 number_of_smoothing_steps=2000,
137                 trigger_depth=0.01, # Depth below which no flow happens
138                 manning=None,          # Mannings Roughness for Culvert
139                 sum_loss=None,
140                 use_velocity_head=False, # FIXME(Ole): Get rid of - always True
141                 use_momentum_jet=False, # FIXME(Ole): Not yet implemented
142                 label=None,
143                 description=None,
144                 update_interval=None,
145                 log_file=False,
146                 discharge_hydrograph=False,
147                 verbose=False):
148       
149
150       
151        # Input check
152       
153        if height is None: height = width
154       
155        assert number_of_barrels >= 1
156        assert use_velocity_head is True or use_velocity_head is False
157       
158        #msg = 'Momentum jet not yet moved to general culvert'
159        #assert use_momentum_jet is False, msg
160        self.use_momentum_jet = use_momentum_jet
161 
162        self.culvert_routine = culvert_routine       
163        self.culvert_description_filename = culvert_description_filename
164        if culvert_description_filename is not None:
165            label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
166            self.rating_curve = ensure_numeric(rating_curve)           
167
168        self.height = height
169        self.width = width
170
171       
172        self.domain = domain
173        self.trigger_depth = trigger_depth       
174               
175        if manning is None: 
176            self.manning = 0.012   # Default roughness for pipe
177       
178        if sum_loss is None:
179            self.sum_loss = 0.0
180           
181       
182                       
183        # Store culvert information
184        self.label = label
185        self.description = description
186        self.culvert_type = type
187        self.number_of_barrels = number_of_barrels
188       
189        # Store options       
190        self.use_velocity_head = use_velocity_head
191
192        if label is None: label = 'culvert_flow'
193        label += '_' + str(id(self)) 
194        self.label = label
195       
196        # File for storing discharge_hydrograph
197        if discharge_hydrograph is True:
198            self.timeseries_filename = label + '_timeseries.csv'
199            fid = open(self.timeseries_filename, 'w')
200            fid.write('time, discharge\n')
201            fid.close()
202
203        # Log file for storing general textual output
204        if log_file is True:       
205            self.log_filename = label + '.log'         
206            log_to_file(self.log_filename, self.label)       
207            log_to_file(self.log_filename, description)
208            log_to_file(self.log_filename, self.culvert_type)       
209        else:
210            self.log_filename = None
211
212
213        # Create the fundamental culvert polygons from polygon
214        P = create_culvert_polygons(end_point0,
215                                    end_point1,
216                                    width=width,   
217                                    height=height,
218                                    number_of_barrels=number_of_barrels)
219        self.culvert_polygons = P
220       
221        # Select enquiry points
222        if enquiry_point0 is None:
223            enquiry_point0 = P['enquiry_point0']
224           
225        if enquiry_point1 is None:
226            enquiry_point1 = P['enquiry_point1']           
227           
228        if verbose is True:
229            pass
230            #plot_polygons([[end_point0, end_point1],
231            #               P['exchange_polygon0'],
232            #               P['exchange_polygon1'],
233            #               [enquiry_point0, 1.005*enquiry_point0],
234            #               [enquiry_point1, 1.005*enquiry_point1]],
235            #              figname='culvert_polygon_output')
236
237           
238           
239        self.enquiry_points = [enquiry_point0, enquiry_point1]
240        self.enquiry_indices = self.get_enquiry_indices()                 
241        self.check_culvert_inside_domain()
242       
243                   
244        # Create inflow object at each end of the culvert.
245        self.openings = []
246        self.openings.append(Inflow(domain,
247                                    polygon=P['exchange_polygon0']))
248        self.openings.append(Inflow(domain,
249                                    polygon=P['exchange_polygon1']))
250                                   
251        # Assume two openings for now: Referred to as 0 and 1
252        assert len(self.openings) == 2
253
254        # Establish initial values at each enquiry point
255        dq = domain.quantities
256        for i, opening in enumerate(self.openings):
257            idx = self.enquiry_indices[i]
258            elevation = dq['elevation'].get_values(location='centroids',
259                                                   indices=[idx])[0]
260            stage = dq['stage'].get_values(location='centroids',
261                                           indices=[idx])[0]
262            opening.elevation = elevation
263            opening.stage = stage
264            opening.depth = stage-elevation
265
266           
267                           
268        # Determine initial pipe direction.
269        # This may change dynamically based on the total energy difference     
270        # Consequently, this may be superfluous
271        delta_z = self.openings[0].elevation - self.openings[1].elevation
272        if delta_z > 0.0:
273            self.inlet = self.openings[0]
274            self.outlet = self.openings[1]
275        else:               
276            self.outlet = self.openings[0]
277            self.inlet = self.openings[1]       
278       
279       
280        # Store basic geometry
281        self.end_points = [end_point0, end_point1]
282        self.vector = P['vector']
283        self.length = P['length']; assert self.length > 0.0
284        if culvert_description_filename is not None:
285            if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
286                msg = 'WARNING: barrel length specified in "%s" (%.2f m)'\
287                    % (culvert_description_filename, 
288                       length)
289                msg += ' does not match distance between specified'
290                msg += ' end points (%.2f m)' %self.length
291                log.critical(msg)
292       
293        self.verbose = verbose
294
295        # Circular index for flow averaging in culvert
296        self.N = N = number_of_smoothing_steps
297        self.Q_list = [0]*N
298        self.i = i
299       
300        # For use with update_interval                       
301        self.last_update = 0.0
302        self.update_interval = update_interval
303       
304        # Create objects to update momentum (a bit crude at this stage). This is used with the momentum jet.
305        xmom0 = General_forcing(domain, 'xmomentum',
306                                polygon=P['exchange_polygon0'])
307
308        xmom1 = General_forcing(domain, 'xmomentum',
309                                polygon=P['exchange_polygon1'])
310
311        ymom0 = General_forcing(domain, 'ymomentum',
312                                polygon=P['exchange_polygon0'])
313
314        ymom1 = General_forcing(domain, 'ymomentum',
315                                polygon=P['exchange_polygon1'])
316
317        self.opening_momentum = [[xmom0, ymom0], [xmom1, ymom1]]
318
319
320
321        # Print some diagnostics to log if requested
322        if self.log_filename is not None:
323            s = 'Culvert Effective Length = %.2f m' %(self.length)
324            log_to_file(self.log_filename, s)
325   
326            s = 'Culvert Direction is %s\n' %str(self.vector)
327            log_to_file(self.log_filename, s)       
328       
329       
330       
331       
332       
333    def __call__(self, domain):
334
335        # Time stuff
336        time = domain.get_time()
337       
338       
339        update = False
340        if self.update_interval is None:
341            # Use next timestep as has been computed in domain.py       
342            delta_t = domain.timestep           
343            update = True
344        else:   
345            # Use update interval
346            delta_t = self.update_interval           
347            if time - self.last_update > self.update_interval or time == 0.0:
348                update = True
349           
350        if self.log_filename is not None:       
351            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
352            log_to_file(self.log_filename, s)
353               
354                               
355        if update is True:
356            self.compute_rates(delta_t)
357           
358   
359        # Execute flow term for each opening
360        # This is where Inflow objects are evaluated using the last rate
361        # that has been calculated
362        #
363        # This will take place at every internal timestep and update the domain
364        for opening in self.openings:
365            opening(domain)
366           
367
368
369    def get_enquiry_indices(self):
370        """Get indices for nearest centroids to self.enquiry_points
371        """
372       
373        domain = self.domain
374       
375        enquiry_indices = []                 
376        for point in self.enquiry_points:
377            # Find nearest centroid
378            N = len(domain)   
379            points = domain.get_centroid_coordinates(absolute=True)
380
381            # Calculate indices in exchange area for this forcing term
382           
383            triangle_id = min_dist = sys.maxint
384            for k in range(N):
385                x, y = points[k,:] # Centroid
386
387                c = point                               
388                distance = (x-c[0])**2+(y-c[1])**2
389                if distance < min_dist:
390                    min_dist = distance
391                    triangle_id = k
392
393                   
394            if triangle_id < sys.maxint:
395                msg = 'found triangle with centroid (%f, %f)'\
396                    %tuple(points[triangle_id, :])
397                msg += ' for point (%f, %f)' %tuple(point)
398               
399                enquiry_indices.append(triangle_id)
400            else:
401                msg = 'Triangle not found for point (%f, %f)' %point
402                raise Exception, msg
403       
404        return enquiry_indices
405
406       
407    def check_culvert_inside_domain(self):
408        """Check that all polygons and enquiry points lie within the mesh.
409        """
410        bounding_polygon = self.domain.get_boundary_polygon()
411        P = self.culvert_polygons
412        for key in P.keys():
413            if key in ['exchange_polygon0', 
414                       'exchange_polygon1']:
415                for point in list(P[key]) + self.enquiry_points:
416                    msg = 'Point %s in polygon %s for culvert %s did not'\
417                        %(str(point), key, self.label)
418                    msg += 'fall within the domain boundary.'
419                    assert is_inside_polygon(point, bounding_polygon), msg
420           
421
422    def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
423        """Adjust Q downwards depending on available water at inlet
424
425           This is a critical step in modelling bridges and Culverts
426           the predicted flow through a structure based on an abstract
427           algorithm can at times request for water that is simply not
428           available due to any number of constrictions that limit the
429           flow approaching the structure In order to ensure that
430           there is adequate flow available certain checks are
431           required There needs to be a check using the Static Water
432           Volume sitting infront of the structure, In addition if the
433           water is moving the available water will be larger than the
434           static volume
435           
436           NOTE To temporarily switch this off for Debugging purposes
437           rem out line in function def compute_rates below
438        """
439   
440        if delta_t < epsilon:
441            # No need to adjust if time step is very small or zero
442            # In this case the possible flow will be very large
443            # anyway.
444            return Q
445       
446        # Short hands
447        domain = self.domain       
448        dq = domain.quantities               
449        time = domain.get_time()
450        I = self.inlet
451        idx = I.exchange_indices   
452
453        # Find triangle with the smallest depth
454        stage = dq['stage'].get_values(location='centroids', 
455                                               indices=[idx])
456        elevation = dq['elevation'].get_values(location='centroids', 
457                                               indices=[idx])       
458        depth = stage-elevation
459        min_depth = min(depth.flat)  # This may lead to errors if edge of area is at a higher level !!!!
460        avg_depth = mean(depth.flat) # Yes, but this one violates the conservation unit tests
461
462
463
464        # FIXME (Ole): If you want these, use log.critical() and
465        # make the statements depend on verbose
466        #print I.depth
467        #print I.velocity
468        #print self.width
469
470        # max_Q Based on Volume Calcs
471
472
473        depth_term = min_depth*I.exchange_area/delta_t
474        if min_depth < 0.2:
475            # Only add velocity term in shallow waters (< 20 cm)
476            # This is a little ad hoc, but maybe it is reasonable
477            velocity_term = self.width*min_depth*I.velocity
478        else:
479            velocity_term = 0.0
480
481        # This one takes approaching water into account   
482        max_Q = max(velocity_term, depth_term)
483
484        # This one preserves Volume
485        #max_Q = depth_term
486
487
488        if self.verbose is True:
489            log.critical('Max_Q = %f' % max_Q)           
490            msg = 'Width = %.2fm, Depth at inlet = %.2f m, Velocity = %.2f m/s.      ' % (self.width, I.depth, I.velocity)
491            msg += 'Max Q = %.2f m^3/s' %(max_Q)
492            log.critical(msg)
493
494        if self.log_filename is not None:               
495            log_to_file(self.log_filename, msg)
496        # New Procedure for assessing the flow available to the Culvert
497        # This routine uses the GET FLOW THROUGH CROSS SECTION
498        #   Need to check Several Polyline however
499        # Firstly 3 sides of the exchange Poly
500        # then only the Line Directly infront of the Polygon
501        # Access polygon Points from   self.inlet.polygon
502     
503        #  The Following computes the flow crossing over 3 sides of the exchange polygon for the structure
504        # Clearly the flow in the culvert can not be more than that flowing toward it through the exhange polygon
505       
506        #q1 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][1:3]) # First Side Segment
507        #q2 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][2:])   # Second Face Segment
508        #q3 =domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'].take([3,0], axis=0)) # Third Side Segment
509        # q4 = domain.get_flow_through_cross_section([self.culvert_polygons['exchange_polygon0'][1:4]][0])
510        #q4=max(q1,0.0)+max(q2,0.0)+max(q3,0.0)
511        # To use only the Flow crossing the 3 sides of the Exchange Polygon use the following Line Only
512        #max_Q=max(q1,q2,q3,q4)
513        # Try Simple Smoothing using Average of 2 approaches
514        #max_Q=(max(q1,q2,q3,q4)+max_Q)/2.0
515        # Calculate the minimum in absolute terms of
516        # the requsted flow and the possible flow
517        Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
518        if self.verbose is True:
519            msg = 'Initial Q Reduced = %.2f m3/s.      ' % (Q_reduced)
520            log.critical(msg)
521
522        if self.log_filename is not None:               
523            log_to_file(self.log_filename, msg)
524        # Now Keep Rolling Average of Computed Discharge to Reduce / Remove Oscillations
525        #  can use delta_t if we want to averageover a time frame for example
526        # N = 5.0/delta_t  Will provide the average over 5 seconds
527
528        self.i=(self.i+1)%self.N
529        self.Q_list[self.i]=Q_reduced
530        Q_reduced = sum(self.Q_list)/len(self.Q_list)
531
532        if self.verbose is True:
533            msg = 'Final Q Reduced = %.2f m3/s.      ' % (Q_reduced)
534            log.critical(msg)
535
536        if self.log_filename is not None:               
537            log_to_file(self.log_filename, msg)
538
539
540        if abs(Q_reduced) < abs(Q): 
541            msg = '%.2fs: Requested flow is ' % time
542            msg += 'greater than what is supported by the smallest '
543            msg += 'depth at inlet exchange area:\n        '
544            msg += 'inlet exchange area: %.2f '% (I.exchange_area) 
545            msg += 'velocity at inlet :%.2f '% (I.velocity)
546            msg += 'Vel* Exch Area = : %.2f '% (I.velocity*avg_depth*self.width)
547            msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\
548                % (avg_depth, I.exchange_area, delta_t)
549            msg += ' = %.2f m^3/s\n        ' % Q_reduced
550            msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
551            msg += 'Note calculate max_Q from V %.2f m^3/s ' % (max_Q)
552            if self.verbose is True:
553                log.critical(msg)
554               
555            if self.log_filename is not None:               
556                log_to_file(self.log_filename, msg)
557       
558        return Q_reduced   
559                       
560           
561    def compute_rates(self, delta_t):
562        """Compute new rates for inlet and outlet
563        """
564
565        # Short hands
566        domain = self.domain       
567        dq = domain.quantities               
568       
569        # Time stuff
570        time = domain.get_time()
571        self.last_update = time
572
573           
574        if hasattr(self, 'log_filename'):
575            log_filename = self.log_filename
576           
577        # Compute stage, energy and velocity at the
578        # enquiry points at each end of the culvert
579        openings = self.openings
580        for i, opening in enumerate(openings):
581            idx = self.enquiry_indices[i]               
582           
583            stage = dq['stage'].get_values(location='centroids',
584                                           indices=[idx])[0]
585            depth = h = stage-opening.elevation
586                                                           
587           
588            # Get velocity                                 
589            xmomentum = dq['xmomentum'].get_values(location='centroids',
590                                                   indices=[idx])[0]
591            ymomentum = dq['xmomentum'].get_values(location='centroids',
592                                                   indices=[idx])[0]
593
594            if h > minimum_allowed_height:
595                u = xmomentum/(h + velocity_protection/h)
596                v = ymomentum/(h + velocity_protection/h)
597            else:
598                u = v = 0.0
599               
600            v_squared = u*u + v*v
601           
602            if self.use_velocity_head is True:
603                velocity_head = 0.5*v_squared/g   
604            else:
605                velocity_head = 0.0
606           
607            opening.total_energy = velocity_head + stage
608            opening.specific_energy = velocity_head + depth
609            opening.stage = stage
610            opening.depth = depth
611            opening.velocity = sqrt(v_squared)
612           
613
614        # We now need to deal with each opening individually
615        # Determine flow direction based on total energy difference
616        delta_total_energy = openings[0].total_energy - openings[1].total_energy
617        if delta_total_energy > 0:
618            inlet = openings[0]
619            outlet = openings[1]
620
621            # FIXME: I think this whole momentum jet thing could be a bit more elegant
622            inlet.momentum = self.opening_momentum[0]
623            outlet.momentum = self.opening_momentum[1]
624        else:
625            inlet = openings[1]
626            outlet = openings[0]
627           
628            inlet.momentum = self.opening_momentum[1]
629            outlet.momentum = self.opening_momentum[0]
630
631            delta_total_energy = -delta_total_energy
632
633        self.inlet = inlet
634        self.outlet = outlet
635           
636        msg = 'Total energy difference is negative'
637        assert delta_total_energy >= 0.0, msg
638
639        # Recompute slope and issue warning if flow is uphill
640        # These values do not enter the computation
641        delta_z = inlet.elevation - outlet.elevation
642        culvert_slope = (delta_z/self.length)
643        if culvert_slope < 0.0:
644            # Adverse gradient - flow is running uphill
645            # Flow will be purely controlled by uphill outlet face
646            if self.verbose is True:
647                log.critical('%.2fs - WARNING: Flow is running uphill.' % time)
648           
649        if self.log_filename is not None:
650            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
651                %(time, self.inlet.stage, self.outlet.stage)
652            log_to_file(self.log_filename, s)
653            s = 'Delta total energy = %.3f' %(delta_total_energy)
654            log_to_file(log_filename, s)
655
656           
657        # Determine controlling energy (driving head) for culvert
658        if inlet.specific_energy > delta_total_energy:
659            # Outlet control
660            driving_head = delta_total_energy
661        else:
662            # Inlet control
663            driving_head = inlet.specific_energy
664           
665
666
667        if self.inlet.depth <= self.trigger_depth:
668            Q = 0.0
669        else:
670            # Calculate discharge for one barrel and
671            # set inlet.rate and outlet.rate
672           
673            if self.culvert_description_filename is not None:
674                try:
675                    Q = interpolate_linearly(driving_head, 
676                                             self.rating_curve[:,0], 
677                                             self.rating_curve[:,1]) 
678                except Below_interval, e:
679                    Q = self.rating_curve[0,1]             
680                    msg = '%.2fs: ' % time
681                    msg += 'Delta head smaller than rating curve minimum: '
682                    msg += str(e)
683                    msg += '\n        '
684                    msg += 'I will use minimum discharge %.2f m^3/s ' % Q
685                    msg += 'for culvert "%s"' % self.label
686                   
687                    if hasattr(self, 'log_filename'):                   
688                        log_to_file(self.log_filename, msg)
689                except Above_interval, e:
690                    Q = self.rating_curve[-1,1]         
691                    msg = '%.2fs: ' % time                 
692                    msg += 'Delta head greater than rating curve maximum: '
693                    msg += str(e)
694                    msg += '\n        '
695                    msg += 'I will use maximum discharge %.2f m^3/s ' % Q
696                    msg += 'for culvert "%s"' % self.label
697                   
698                    if self.log_filename is not None:                   
699                        log_to_file(self.log_filename, msg)
700            else:
701                # User culvert routine
702                Q, barrel_velocity, culvert_outlet_depth =\
703                    self.culvert_routine(inlet.depth,
704                                         outlet.depth,
705                                         inlet.velocity,
706                                         outlet.velocity,
707                                         inlet.specific_energy, 
708                                         delta_total_energy, 
709                                         g,
710                                         culvert_length=self.length,
711                                         culvert_width=self.width,
712                                         culvert_height=self.height,
713                                         culvert_type=self.culvert_type,
714                                         manning=self.manning,
715                                         sum_loss=self.sum_loss,
716                                         log_filename=self.log_filename)
717           
718           
719       
720        # Adjust discharge for multiple barrels
721        Q *= self.number_of_barrels
722
723        # Adjust discharge for available water at the inlet
724        Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
725       
726        self.inlet.rate = -Q
727        self.outlet.rate = Q
728
729
730        # Momentum jet stuff
731        if self.use_momentum_jet is True:
732
733
734            # Compute barrel momentum
735            barrel_momentum = barrel_velocity*culvert_outlet_depth
736
737            if self.log_filename is not None:                                   
738                s = 'Barrel velocity = %f' %barrel_velocity
739                log_to_file(self.log_filename, s)
740
741            # Compute momentum vector at outlet
742            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
743               
744            if self.log_filename is not None:               
745                s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
746                log_to_file(self.log_filename, s)
747
748
749            # Update momentum       
750            if delta_t > 0.0:
751                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
752                xmomentum_rate /= delta_t
753                   
754                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
755                ymomentum_rate /= delta_t
756                       
757                if self.log_filename is not None:               
758                    s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
759                    log_to_file(self.log_filename, s)                   
760            else:
761                xmomentum_rate = ymomentum_rate = 0.0
762
763
764            # Set momentum rates for outlet jet
765            outlet.momentum[0].rate = xmomentum_rate
766            outlet.momentum[1].rate = ymomentum_rate
767
768            # Remember this value for next step (IMPORTANT)
769            outlet.momentum[0].value = outlet_mom_x
770            outlet.momentum[1].value = outlet_mom_y                   
771
772            if int(domain.time*100) % 100 == 0:
773
774                if self.log_filename is not None:               
775                    s = 'T=%.5f, Culvert Discharge = %.3f f'\
776                        %(time, Q)
777                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
778                        %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
779                    s += ' Momentum rate: (%.4f, %.4f)'\
780                        %(xmomentum_rate, ymomentum_rate)                   
781                    s+='Outlet Vel= %.3f'\
782                        %(barrel_velocity)
783                    log_to_file(self.log_filename, s)
784
785
786            # Execute momentum terms
787            # This is where Inflow objects are evaluated and update the domain
788                self.outlet.momentum[0](domain)
789                self.outlet.momentum[1](domain)       
790           
791
792           
793        # Log timeseries to file
794        try:
795            fid = open(self.timeseries_filename, 'a')       
796        except:
797            pass
798        else:   
799            fid.write('%.2f, %.2f\n' %(time, Q))
800            fid.close()
801
802        # Store value of time
803        self.last_time = time
804
805
806           
807
808
809
810                           
811# OBSOLETE (Except for momentum jet in Culvert_flow_energy)   
812class Culvert_flow_rating:
813    """Culvert flow - transfer water from one hole to another
814   
815
816    Input: Two points, pipe_size (either diameter or width, height),
817    mannings_rougness,
818    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
819    top-down_blockage_factor and bottom_up_blockage_factor
820   
821    """ 
822
823    def __init__(self,
824                 domain,
825                 culvert_description_filename=None,
826                 end_point0=None, 
827                 end_point1=None,
828                 enquiry_point0=None, 
829                 enquiry_point1=None,
830                 update_interval=None,
831                 log_file=False,
832                 discharge_hydrograph=False,
833                 verbose=False):
834       
835
836       
837        label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
838       
839               
840        # Store culvert information
841        self.label = label
842        self.description = description
843        self.culvert_type = type
844        self.rating_curve = ensure_numeric(rating_curve)
845        self.number_of_barrels = number_of_barrels
846
847        if label is None: label = 'culvert_flow'
848        label += '_' + str(id(self)) 
849        self.label = label
850       
851        # File for storing discharge_hydrograph
852        if discharge_hydrograph is True:
853            self.timeseries_filename = label + '_timeseries.csv'
854            fid = open(self.timeseries_filename, 'w')
855            fid.write('time, discharge\n')
856            fid.close()
857
858        # Log file for storing general textual output
859        if log_file is True:       
860            self.log_filename = label + '.log'         
861            log_to_file(self.log_filename, self.label)       
862            log_to_file(self.log_filename, description)
863            log_to_file(self.log_filename, self.culvert_type)       
864
865
866        # Create the fundamental culvert polygons from POLYGON
867        #if self.culvert_type == 'circle':
868        #    # Redefine width and height for use with create_culvert_polygons
869        #    width = height = diameter
870       
871        P = create_culvert_polygons(end_point0,
872                                    end_point1,
873                                    width=width,   
874                                    height=height,
875                                    number_of_barrels=number_of_barrels)
876       
877        # Select enquiry points
878        if enquiry_point0 is None:
879            enquiry_point0 = P['enquiry_point0']
880           
881        if enquiry_point1 is None:
882            enquiry_point1 = P['enquiry_point1']           
883           
884        if verbose is True:
885            pass
886            #plot_polygons([[end_point0, end_point1],
887            #               P['exchange_polygon0'],
888            #               P['exchange_polygon1'],
889            #               [enquiry_point0, 1.005*enquiry_point0],
890            #               [enquiry_point1, 1.005*enquiry_point1]],                           
891            #              figname='culvert_polygon_output')
892
893           
894           
895        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
896
897        self.enquiry_indices = []                 
898        for point in self.enquiry_points:
899            # Find nearest centroid
900            N = len(domain)   
901            points = domain.get_centroid_coordinates(absolute=True)
902
903            # Calculate indices in exchange area for this forcing term
904           
905            triangle_id = min_dist = sys.maxint
906            for k in range(N):
907                x, y = points[k,:] # Centroid
908
909                c = point                               
910                distance = (x-c[0])**2+(y-c[1])**2
911                if distance < min_dist:
912                    min_dist = distance
913                    triangle_id = k
914
915                   
916            if triangle_id < sys.maxint:
917                msg = 'found triangle with centroid (%f, %f)'\
918                    %tuple(points[triangle_id, :])
919                msg += ' for point (%f, %f)' %tuple(point)
920               
921                self.enquiry_indices.append(triangle_id)
922            else:
923                msg = 'Triangle not found for point (%f, %f)' %point
924                raise Exception, msg
925       
926                         
927
928        # Check that all polygons lie within the mesh.
929        bounding_polygon = domain.get_boundary_polygon()
930        for key in P.keys():
931            if key in ['exchange_polygon0', 
932                       'exchange_polygon1']:
933                for point in list(P[key]) + self.enquiry_points:
934                    msg = 'Point %s in polygon %s for culvert %s did not'\
935                        %(str(point), key, self.label)
936                    msg += 'fall within the domain boundary.'
937                    assert is_inside_polygon(point, bounding_polygon), msg
938       
939                   
940        # Create inflow object at each end of the culvert.
941        self.openings = []
942        self.openings.append(Inflow(domain,
943                                    polygon=P['exchange_polygon0']))
944
945        self.openings.append(Inflow(domain,
946                                    polygon=P['exchange_polygon1']))                                   
947
948
949
950        dq = domain.quantities                                           
951        for i, opening in enumerate(self.openings):                           
952            elevation = dq['elevation'].get_values(location='centroids',
953                                                   indices=[self.enquiry_indices[i]])           
954            opening.elevation = elevation
955            opening.stage = elevation # Simple assumption that culvert is dry initially
956
957        # Assume two openings for now: Referred to as 0 and 1
958        assert len(self.openings) == 2
959                           
960        # Determine pipe direction     
961        self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation
962        if delta_z > 0.0:
963            self.inlet = self.openings[0]
964            self.outlet = self.openings[1]
965        else:               
966            self.outlet = self.openings[0]
967            self.inlet = self.openings[1]       
968       
969       
970        # Store basic geometry
971        self.end_points = [end_point0, end_point1]
972        self.vector = P['vector']
973        self.length = P['length']; assert self.length > 0.0
974        if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
975            msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length)
976            msg += ' does not match distance between specified'
977            msg += ' end points (%.2f m)' %self.length
978            log.critical(msg)
979       
980        self.verbose = verbose
981        self.last_update = 0.0 # For use with update_interval       
982        self.last_time = 0.0               
983        self.update_interval = update_interval
984       
985
986        # Print something
987        if hasattr(self, 'log_filename'):
988            s = 'Culvert Effective Length = %.2f m' %(self.length)
989            log_to_file(self.log_filename, s)
990   
991            s = 'Culvert Direction is %s\n' %str(self.vector)
992            log_to_file(self.log_filename, s)       
993       
994       
995       
996       
997       
998    def __call__(self, domain):
999
1000        # Time stuff
1001        time = domain.get_time()
1002       
1003       
1004        update = False
1005        if self.update_interval is None:
1006            update = True
1007            delta_t = domain.timestep # Next timestep has been computed in domain.py
1008        else:   
1009            if time - self.last_update > self.update_interval or time == 0.0:
1010                update = True
1011            delta_t = self.update_interval
1012           
1013        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
1014        if hasattr(self, 'log_filename'):           
1015            log_to_file(self.log_filename, s)
1016               
1017                               
1018        if update is True:
1019            self.last_update = time
1020       
1021            dq = domain.quantities
1022                       
1023            # Get average water depths at each opening       
1024            openings = self.openings   # There are two Opening [0] and [1]
1025            for i, opening in enumerate(openings):
1026               
1027                # Compute mean values of selected quantitites in the
1028                # enquiry area in front of the culvert
1029               
1030                stage = dq['stage'].get_values(location='centroids',
1031                                               indices=[self.enquiry_indices[i]])
1032               
1033                # Store current average stage and depth with each opening object
1034                opening.depth = stage - opening.elevation
1035                opening.stage = stage
1036
1037               
1038
1039            #################  End of the FOR loop ################################################
1040
1041            # We now need to deal with each opening individually
1042               
1043            # Determine flow direction based on total energy difference
1044
1045            delta_w = self.inlet.stage - self.outlet.stage
1046           
1047            if hasattr(self, 'log_filename'):
1048                s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time, 
1049                                                                           self.inlet.stage,
1050                                                                           self.outlet.stage)
1051                log_to_file(self.log_filename, s)
1052
1053
1054            if self.inlet.depth <= 0.01:
1055                Q = 0.0
1056            else:
1057                # Calculate discharge for one barrel and set inlet.rate and outlet.rate
1058               
1059                try:
1060                    Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1]) 
1061                except Below_interval, e:
1062                    Q = self.rating_curve[0,1]             
1063                    msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
1064                    msg += str(e)
1065                    msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
1066                        %(Q, self.label)
1067                    if hasattr(self, 'log_filename'):                   
1068                        log_to_file(self.log_filename, msg)
1069                except Above_interval, e:
1070                    Q = self.rating_curve[-1,1]         
1071                    msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
1072                    msg += str(e)
1073                    msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
1074                        %(Q, self.label)
1075                    if hasattr(self, 'log_filename'):                   
1076                        log_to_file(self.log_filename, msg)                   
1077
1078               
1079               
1080           
1081            # Adjust discharge for multiple barrels
1082            Q *= self.number_of_barrels
1083           
1084
1085            # Adjust Q downwards depending on available water at inlet
1086            stage = self.inlet.get_quantity_values(quantity_name='stage')
1087            elevation = self.inlet.get_quantity_values(quantity_name='elevation')
1088            depth = stage-elevation
1089           
1090           
1091            V = 0
1092            for i, d in enumerate(depth):
1093                V += d * domain.areas[i]
1094           
1095            dt = delta_t           
1096            if Q*dt > V:
1097           
1098                Q_reduced = 0.9*V/dt # Reduce with safety factor
1099               
1100                msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
1101                msg += 'greater than current volume (V) at inlet.\n'
1102                msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
1103               
1104                if self.verbose is True:
1105                    log.critical(msg)
1106                if hasattr(self, 'log_filename'):                   
1107                    log_to_file(self.log_filename, msg)                                       
1108               
1109                Q = Q_reduced
1110       
1111            self.inlet.rate = -Q
1112            self.outlet.rate = Q
1113
1114            # Log timeseries to file
1115            try:
1116                fid = open(self.timeseries_filename, 'a')       
1117            except:
1118                pass
1119            else:   
1120                fid.write('%.2f, %.2f\n' %(time, Q))
1121                fid.close()
1122
1123            # Store value of time
1124            self.last_time = time
1125
1126           
1127   
1128        # Execute flow term for each opening
1129        # This is where Inflow objects are evaluated using the last rate that has been calculated
1130        #
1131        # This will take place at every internal timestep and update the domain
1132        for opening in self.openings:
1133            opening(domain)
1134
1135
1136
1137       
1138       
1139
1140class Culvert_flow_energy:
1141    """Culvert flow - transfer water from one hole to another
1142   
1143    Using Momentum as Calculated by Culvert Flow !!
1144    Could be Several Methods Investigated to do This !!!
1145
1146    2008_May_08
1147    To Ole:
1148    OK so here we need to get the Polygon Creating code to create
1149    polygons for the culvert Based on
1150    the two input Points (X0,Y0) and (X1,Y1) - need to be passed
1151    to create polygon
1152
1153    The two centers are now passed on to create_polygon.
1154   
1155
1156    Input: Two points, pipe_size (either diameter or width, height),
1157    mannings_rougness,
1158    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
1159    top-down_blockage_factor and bottom_up_blockage_factor
1160   
1161   
1162    And the Delta H enquiry should be change from Openings in line 412
1163    to the enquiry Polygons infront of the culvert
1164    At the moment this script uses only Depth, later we can change it to
1165    Energy...
1166
1167    Once we have Delta H can calculate a Flow Rate and from Flow Rate
1168    an Outlet Velocity
1169    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
1170       
1171    Invert levels are optional. If left out they will default to the
1172    elevation at the opening.
1173       
1174    """ 
1175
1176    def __init__(self,
1177                 domain,
1178                 label=None, 
1179                 description=None,
1180                 end_point0=None, 
1181                 end_point1=None,
1182                 width=None,
1183                 height=None,
1184                 diameter=None,
1185                 manning=None,          # Mannings Roughness for Culvert
1186                 invert_level0=None,    # Invert level at opening 0
1187                 invert_level1=None,    # Invert level at opening 1
1188                 loss_exit=None,
1189                 loss_entry=None,
1190                 loss_bend=None,
1191                 loss_special=None,
1192                 blockage_topdwn=None,
1193                 blockage_bottup=None,
1194                 culvert_routine=None,
1195                 number_of_barrels=1,
1196                 enquiry_point0=None, 
1197                 enquiry_point1=None,
1198                 update_interval=None,
1199                 verbose=False):
1200       
1201        # Input check
1202        if diameter is not None:
1203            self.culvert_type = 'circle'
1204            self.diameter = diameter
1205            if height is not None or width is not None:
1206                msg = 'Either diameter or width&height must be specified, '
1207                msg += 'but not both.'
1208                raise Exception, msg
1209        else:
1210            if height is not None:
1211                if width is None:
1212                    self.culvert_type = 'square'                               
1213                    width = height
1214                else:
1215                    self.culvert_type = 'rectangle'
1216            elif width is not None:
1217                if height is None:
1218                    self.culvert_type = 'square'                               
1219                    height = width
1220            else:
1221                msg = 'Either diameter or width&height must be specified.'
1222                raise Exception, msg               
1223               
1224            if height == width:
1225                self.culvert_type = 'square'                                               
1226               
1227            self.height = height
1228            self.width = width
1229
1230           
1231        assert self.culvert_type in ['circle', 'square', 'rectangle']
1232       
1233        assert number_of_barrels >= 1
1234        self.number_of_barrels = number_of_barrels
1235       
1236       
1237        # Set defaults
1238        if manning is None: manning = 0.012   # Default roughness for pipe
1239        if loss_exit is None: loss_exit = 1.00
1240        if loss_entry is None: loss_entry = 0.50
1241        if loss_bend is None: loss_bend=0.00
1242        if loss_special is None: loss_special=0.00
1243        if blockage_topdwn is None: blockage_topdwn=0.00
1244        if blockage_bottup is None: blockage_bottup=0.00
1245        if culvert_routine is None: 
1246            culvert_routine=boyd_generalised_culvert_model
1247           
1248        if label is None: label = 'culvert_flow'
1249        label += '_' + str(id(self)) 
1250        self.label = label
1251       
1252        # File for storing culvert quantities
1253        self.timeseries_filename = label + '_timeseries.csv'
1254        fid = open(self.timeseries_filename, 'w')
1255        fid.write('time, E0, E1, Velocity, Discharge\n')
1256        fid.close()
1257
1258        # Log file for storing general textual output
1259        self.log_filename = label + '.log'         
1260        log_to_file(self.log_filename, self.label)       
1261        log_to_file(self.log_filename, description)
1262        log_to_file(self.log_filename, self.culvert_type)       
1263
1264
1265        # Create the fundamental culvert polygons from POLYGON
1266        if self.culvert_type == 'circle':
1267            # Redefine width and height for use with create_culvert_polygons
1268            width = height = diameter
1269       
1270        P = create_culvert_polygons(end_point0,
1271                                    end_point1,
1272                                    width=width,   
1273                                    height=height,
1274                                    number_of_barrels=number_of_barrels)
1275       
1276        # Select enquiry points
1277        if enquiry_point0 is None:
1278            enquiry_point0 = P['enquiry_point0']
1279           
1280        if enquiry_point1 is None:
1281            enquiry_point1 = P['enquiry_point1']           
1282           
1283        if verbose is True:
1284            pass
1285            #plot_polygons([[end_point0, end_point1],
1286            #               P['exchange_polygon0'],
1287            #               P['exchange_polygon1'],
1288            #               [enquiry_point0, 1.005*enquiry_point0],
1289            #               [enquiry_point1, 1.005*enquiry_point1]],
1290            #              figname='culvert_polygon_output')
1291
1292
1293        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
1294       
1295       
1296        self.enquiry_indices = []                 
1297        for point in self.enquiry_points:
1298            # Find nearest centroid
1299            N = len(domain)   
1300            points = domain.get_centroid_coordinates(absolute=True)
1301
1302            # Calculate indices in exchange area for this forcing term
1303           
1304            triangle_id = min_dist = sys.maxint
1305            for k in range(N):
1306                x, y = points[k,:] # Centroid
1307
1308                c = point                               
1309                distance = (x-c[0])**2+(y-c[1])**2
1310                if distance < min_dist:
1311                    min_dist = distance
1312                    triangle_id = k
1313
1314                   
1315            if triangle_id < sys.maxint:
1316                msg = 'found triangle with centroid (%f, %f)'\
1317                    %tuple(points[triangle_id, :])
1318                msg += ' for point (%f, %f)' %tuple(point)
1319               
1320                self.enquiry_indices.append(triangle_id)
1321            else:
1322                msg = 'Triangle not found for point (%f, %f)' %point
1323                raise Exception, msg
1324       
1325                         
1326
1327           
1328           
1329
1330        # Check that all polygons lie within the mesh.
1331        bounding_polygon = domain.get_boundary_polygon()
1332        for key in P.keys():
1333            if key in ['exchange_polygon0', 
1334                       'exchange_polygon1']:
1335                for point in P[key]:
1336               
1337                    msg = 'Point %s in polygon %s for culvert %s did not'\
1338                        %(str(point), key, self.label)
1339                    msg += 'fall within the domain boundary.'
1340                    assert is_inside_polygon(point, bounding_polygon), msg
1341       
1342
1343        # Create inflow object at each end of the culvert.
1344        self.openings = []
1345        self.openings.append(Inflow(domain,
1346                                    polygon=P['exchange_polygon0']))
1347
1348        self.openings.append(Inflow(domain,
1349                                    polygon=P['exchange_polygon1']))                                   
1350
1351
1352        # Assume two openings for now: Referred to as 0 and 1
1353        assert len(self.openings) == 2
1354       
1355        # Store basic geometry
1356        self.end_points = [end_point0, end_point1]
1357        self.invert_levels = [invert_level0, invert_level1]               
1358        #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
1359        #self.enquiry_polylines = [P['enquiry_polygon0'][:2],
1360        #                          P['enquiry_polygon1'][:2]]
1361        self.vector = P['vector']
1362        self.length = P['length']; assert self.length > 0.0
1363        self.verbose = verbose
1364        self.last_time = 0.0       
1365        self.last_update = 0.0 # For use with update_interval       
1366        self.update_interval = update_interval
1367       
1368
1369        # Store hydraulic parameters
1370        self.manning = manning
1371        self.loss_exit = loss_exit
1372        self.loss_entry = loss_entry
1373        self.loss_bend = loss_bend
1374        self.loss_special = loss_special
1375        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
1376        self.blockage_topdwn = blockage_topdwn
1377        self.blockage_bottup = blockage_bottup
1378
1379        # Store culvert routine
1380        self.culvert_routine = culvert_routine
1381
1382       
1383        # Create objects to update momentum (a bit crude at this stage)
1384        xmom0 = General_forcing(domain, 'xmomentum',
1385                                polygon=P['exchange_polygon0'])
1386
1387        xmom1 = General_forcing(domain, 'xmomentum',
1388                                polygon=P['exchange_polygon1'])
1389
1390        ymom0 = General_forcing(domain, 'ymomentum',
1391                                polygon=P['exchange_polygon0'])
1392
1393        ymom1 = General_forcing(domain, 'ymomentum',
1394                                polygon=P['exchange_polygon1'])
1395
1396        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
1397       
1398
1399        # Print something
1400        s = 'Culvert Effective Length = %.2f m' %(self.length)
1401        log_to_file(self.log_filename, s)
1402   
1403        s = 'Culvert Direction is %s\n' %str(self.vector)
1404        log_to_file(self.log_filename, s)       
1405       
1406       
1407    def __call__(self, domain):
1408
1409        log_filename = self.log_filename
1410         
1411        # Time stuff
1412        time = domain.get_time()
1413       
1414        # Short hand
1415        dq = domain.quantities
1416               
1417
1418        update = False
1419        if self.update_interval is None:
1420            update = True
1421            delta_t = domain.timestep # Next timestep has been computed in domain.py
1422        else:   
1423            if time - self.last_update > self.update_interval or time == 0.0:
1424                update = True
1425            delta_t = self.update_interval
1426           
1427        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
1428        if hasattr(self, 'log_filename'):           
1429            log_to_file(log_filename, s)
1430               
1431                               
1432        if update is True:
1433            self.last_update = time
1434                       
1435            msg = 'Time did not advance'
1436            if time > 0.0: assert delta_t > 0.0, msg
1437
1438
1439            # Get average water depths at each opening       
1440            openings = self.openings   # There are two Opening [0] and [1]
1441            for i, opening in enumerate(openings):
1442               
1443                # Compute mean values of selected quantitites in the
1444                # exchange area in front of the culvert
1445                     
1446                stage = opening.get_quantity_values(quantity_name='stage')
1447                w = mean(stage) # Average stage
1448
1449                # Use invert level instead of elevation if specified
1450                invert_level = self.invert_levels[i]
1451                if invert_level is not None:
1452                    z = invert_level
1453                else:
1454                    elevation = opening.get_quantity_values(quantity_name='elevation')
1455                    z = mean(elevation) # Average elevation
1456
1457                # Estimated depth above the culvert inlet
1458                d = w - z  # Used for calculations involving depth
1459                if d < 0.0:
1460                    # This is possible since w and z are taken at different locations
1461                    #msg = 'D < 0.0: %f' %d
1462                    #raise msg
1463                    d = 0.0
1464               
1465
1466                # Ratio of depth to culvert height.
1467                # If ratio > 1 then culvert is running full
1468                if self.culvert_type == 'circle':
1469                    ratio = d/self.diameter
1470                else:   
1471                    ratio = d/self.height 
1472                opening.ratio = ratio
1473                   
1474                   
1475                # Average measures of energy in front of this opening
1476               
1477                id = [self.enquiry_indices[i]]
1478                stage = dq['stage'].get_values(location='centroids',
1479                                               indices=id)
1480                elevation = dq['elevation'].get_values(location='centroids',
1481                                                       indices=id)                                               
1482                xmomentum = dq['xmomentum'].get_values(location='centroids',
1483                                                       indices=id)                                               
1484                ymomentum = dq['xmomentum'].get_values(location='centroids',
1485                                                       indices=id)                                                                                             
1486                depth = stage-elevation
1487                if depth > 0.0:
1488                    u = xmomentum/(depth + velocity_protection/depth)
1489                    v = ymomentum/(depth + velocity_protection/depth)
1490                else:
1491                    u = v = 0.0
1492
1493                   
1494                opening.total_energy = 0.5*(u*u + v*v)/g + stage
1495
1496                # Store current average stage and depth with each opening object
1497                opening.depth = d
1498                opening.depth_trigger = d # Use this for now
1499                opening.stage = w
1500                opening.elevation = z
1501               
1502
1503            #################  End of the FOR loop ################################################
1504
1505            # We now need to deal with each opening individually
1506               
1507            # Determine flow direction based on total energy difference
1508            delta_Et = openings[0].total_energy - openings[1].total_energy
1509
1510            if delta_Et > 0:
1511                inlet = openings[0]
1512                outlet = openings[1]
1513
1514                inlet.momentum = self.opening_momentum[0]
1515                outlet.momentum = self.opening_momentum[1]
1516
1517            else:
1518                inlet = openings[1]
1519                outlet = openings[0]
1520
1521                inlet.momentum = self.opening_momentum[1]
1522                outlet.momentum = self.opening_momentum[0]
1523               
1524                delta_Et = -delta_Et
1525
1526            self.inlet = inlet
1527            self.outlet = outlet
1528               
1529            msg = 'Total energy difference is negative'
1530            assert delta_Et > 0.0, msg
1531
1532            delta_h = inlet.stage - outlet.stage
1533            delta_z = inlet.elevation - outlet.elevation
1534            culvert_slope = (delta_z/self.length)
1535
1536            if culvert_slope < 0.0:
1537                # Adverse gradient - flow is running uphill
1538                # Flow will be purely controlled by uphill outlet face
1539                if self.verbose is True:
1540                    log.critical('WARNING: Flow is running uphill. Watch Out! '
1541                                 'inlet.elevation=%s, outlet.elevation%s'
1542                                 % (str(inlet.elevation), str(outlet.elevation)))
1543
1544
1545            s = 'Delta total energy = %.3f' %(delta_Et)
1546            log_to_file(log_filename, s)
1547
1548
1549            # Calculate discharge for one barrel and set inlet.rate and outlet.rate
1550            Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
1551       
1552            # Adjust discharge for multiple barrels
1553            Q *= self.number_of_barrels
1554
1555            # Compute barrel momentum
1556            barrel_momentum = barrel_velocity*culvert_outlet_depth
1557                   
1558            s = 'Barrel velocity = %f' %barrel_velocity
1559            log_to_file(log_filename, s)
1560
1561            # Compute momentum vector at outlet
1562            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
1563               
1564            s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
1565            log_to_file(log_filename, s)
1566
1567            # Log timeseries to file
1568            fid = open(self.timeseries_filename, 'a')       
1569            fid.write('%f, %f, %f, %f, %f\n'\
1570                          %(time, 
1571                            openings[0].total_energy,
1572                            openings[1].total_energy,
1573                            barrel_velocity,
1574                            Q))
1575            fid.close()
1576
1577            # Update momentum       
1578            if delta_t > 0.0:
1579                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
1580                xmomentum_rate /= delta_t
1581                   
1582                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
1583                ymomentum_rate /= delta_t
1584                       
1585                s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
1586                log_to_file(log_filename, s)                   
1587            else:
1588                xmomentum_rate = ymomentum_rate = 0.0
1589
1590
1591            # Set momentum rates for outlet jet
1592            outlet.momentum[0].rate = xmomentum_rate
1593            outlet.momentum[1].rate = ymomentum_rate
1594
1595            # Remember this value for next step (IMPORTANT)
1596            outlet.momentum[0].value = outlet_mom_x
1597            outlet.momentum[1].value = outlet_mom_y                   
1598
1599            if int(domain.time*100) % 100 == 0:
1600                s = 'T=%.5f, Culvert Discharge = %.3f f'\
1601                    %(time, Q)
1602                s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
1603                     %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
1604                s += ' Momentum rate: (%.4f, %.4f)'\
1605                     %(xmomentum_rate, ymomentum_rate)                   
1606                s+='Outlet Vel= %.3f'\
1607                    %(barrel_velocity)
1608                log_to_file(log_filename, s)
1609           
1610            # Store value of time
1611            self.last_time = time
1612               
1613
1614
1615        # Execute flow term for each opening
1616        # This is where Inflow objects are evaluated and update the domain
1617        for opening in self.openings:
1618            opening(domain)
1619
1620        # Execute momentum terms
1621        # This is where Inflow objects are evaluated and update the domain
1622        self.outlet.momentum[0](domain)
1623        self.outlet.momentum[1](domain)       
1624           
1625
1626
1627Culvert_flow = Culvert_flow_general       
Note: See TracBrowser for help on using the repository browser.