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

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

Work on culvert testing - didn't get very far at all.

File size: 34.3 KB
Line 
1from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
2from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
3from anuga.utilities.system_tools import log_to_file
4from anuga.utilities.polygon import inside_polygon
5from anuga.utilities.polygon import is_inside_polygon
6from anuga.utilities.polygon import plot_polygons
7
8from anuga.utilities.numerical_tools import ensure_numeric
9from Numeric import allclose
10
11import sys
12
13class Below_interval(Exception): pass 
14class Above_interval(Exception): pass
15
16
17
18# FIXME(Ole): Write in C and reuse this function by similar code in interpolate.py
19def interpolate_linearly(x, xvec, yvec):
20         
21    # Find appropriate slot           
22    i = 0
23    while x > xvec[i]: i += 1
24
25    if i == 0: 
26        msg = 'Value provided = %.2f, interpolation minimum = %.2f.' %(x, xvec[0])
27        raise Below_interval, msg
28       
29    if i == len(xvec): 
30        msg = 'Value provided = %.2f, interpolation maximum = %.2f.' %(x, xvec[-1])
31        raise Above_interval, msg       
32       
33   
34    x0 = xvec[i-1]
35    x1 = xvec[i]           
36    alpha = (x - x0)/(x1 - x0) 
37           
38    y0 = yvec[i-1]
39    y1 = yvec[i]                       
40    y = alpha*y1 + (1-alpha)*y0
41   
42    return y
43           
44
45           
46def read_culvert_description(culvert_description_filename):           
47   
48    # Read description file
49    fid = open(culvert_description_filename)
50   
51    read_rating_curve_data = False       
52    rating_curve = []
53    for i, line in enumerate(fid.readlines()):
54       
55        if read_rating_curve_data is True:
56            fields = line.split(',')
57            head_difference = float(fields[0].strip())
58            flow_rate = float(fields[1].strip())               
59            barrel_velocity = float(fields[2].strip())
60           
61            rating_curve.append( [head_difference, flow_rate, barrel_velocity] ) 
62       
63        if i == 0:
64            # Header
65            continue
66        if i == 1:
67            # Metadata
68            fields = line.split(',')
69            label=fields[0].strip()
70            type=fields[1].strip().lower()
71            assert type in ['box', 'pipe']
72           
73            width=float(fields[2].strip())
74            height=float(fields[3].strip())
75            length=float(fields[4].strip())
76            number_of_barrels=int(fields[5].strip())
77            #fields[6] refers to losses
78            description=fields[7].strip()               
79               
80        if line.strip() == '': continue # Skip blanks
81
82        if line.startswith('Rating'):
83            read_rating_curve_data = True
84            # Flow data follows
85               
86    fid.close()
87   
88    return label, type, width, height, length, number_of_barrels, description, rating_curve
89   
90
91class Culvert_flow_rating:
92    """Culvert flow - transfer water from one hole to another
93   
94
95    Input: Two points, pipe_size (either diameter or width, height),
96    mannings_rougness,
97    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
98    top-down_blockage_factor and bottom_up_blockage_factor
99   
100    """ 
101
102    def __init__(self,
103                 domain,
104                 culvert_description_filename=None,
105                 end_point0=None, 
106                 end_point1=None,
107                 update_interval=None,
108                 log_file=False,
109                 discharge_hydrograph=False,
110                 verbose=False):
111       
112        from Numeric import sqrt, sum
113
114       
115        label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
116       
117               
118        # Store culvert information
119        self.label = label
120        self.description = description
121        self.culvert_type = type
122        self.rating_curve = ensure_numeric(rating_curve)
123        self.number_of_barrels = number_of_barrels
124
125        if label is None: label = 'culvert_flow'
126        label += '_' + str(id(self)) 
127        self.label = label
128       
129        # File for storing discharge_hydrograph
130        if discharge_hydrograph is True:
131            self.timeseries_filename = label + '_timeseries.csv'
132            fid = open(self.timeseries_filename, 'w')
133            fid.write('time, discharge\n')
134            fid.close()
135
136        # Log file for storing general textual output
137        if log_file is True:       
138            self.log_filename = label + '.log'         
139            log_to_file(self.log_filename, self.label)       
140            log_to_file(self.log_filename, description)
141            log_to_file(self.log_filename, self.culvert_type)       
142
143
144        # Create the fundamental culvert polygons from POLYGON
145        #if self.culvert_type == 'circle':
146        #    # Redefine width and height for use with create_culvert_polygons
147        #    width = height = diameter
148       
149        P = create_culvert_polygons(end_point0,
150                                    end_point1,
151                                    width=width,   
152                                    height=height,
153                                    number_of_barrels=number_of_barrels)
154       
155        if verbose is True:
156            pass
157            #plot_polygons([[end_point0, end_point1],
158            #               P['exchange_polygon0'],
159            #               P['exchange_polygon1'],
160            #               P['enquiry_polygon0'],
161            #               P['enquiry_polygon1']],
162            #              figname='culvert_polygon_output')
163
164                         
165        # Compute the average point for enquiry
166        enquiry_point0 = sum(P['enquiry_polygon0'][:2])/2 
167        enquiry_point1 = sum(P['enquiry_polygon1'][:2])/2         
168       
169        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
170        self.enquiry_indices = []                 
171        for point in self.enquiry_points:
172            # Find nearest centroid
173            N = len(domain)   
174            points = domain.get_centroid_coordinates(absolute=True)
175
176            # Calculate indices in exchange area for this forcing term
177           
178            triangle_id = min_dist = sys.maxint
179            for k in range(N):
180                x, y = points[k,:] # Centroid
181
182                c = point                               
183                distance = (x-c[0])**2+(y-c[1])**2
184                if distance < min_dist:
185                    min_dist = distance
186                    triangle_id = k
187
188                   
189            if triangle_id < sys.maxint:
190                msg = 'found triangle with centroid (%f, %f)'\
191                    %tuple(points[triangle_id, :])
192                msg += ' for point (%f, %f)' %tuple(point)
193               
194                self.enquiry_indices.append(triangle_id)
195            else:
196                msg = 'Triangle not found for point (%f, %f)' %point
197                raise Exception, msg
198       
199                         
200
201        # Check that all polygons lie within the mesh.
202        bounding_polygon = domain.get_boundary_polygon()
203        for key in P.keys():
204            if key in ['exchange_polygon0', 
205                       'exchange_polygon1']:
206                for point in list(P[key]) + self.enquiry_points:
207                    msg = 'Point %s in polygon %s for culvert %s did not'\
208                        %(str(point), key, self.label)
209                    msg += 'fall within the domain boundary.'
210                    assert is_inside_polygon(point, bounding_polygon), msg
211       
212                   
213        # Create inflow object at each end of the culvert.
214        self.openings = []
215        self.openings.append(Inflow(domain,
216                                    polygon=P['exchange_polygon0']))
217
218        self.openings.append(Inflow(domain,
219                                    polygon=P['exchange_polygon1']))                                   
220
221
222
223        dq = domain.quantities                                           
224        for i, opening in enumerate(self.openings):                           
225            #elevation = dq['elevation'].get_values(location='centroids',
226            #                                      interpolation_points=[self.enquiry_points[i]])
227           
228            elevation = dq['elevation'].get_values(location='centroids',
229                                                   indices=[self.enquiry_indices[i]])           
230            opening.elevation = elevation
231
232        # Assume two openings for now: Referred to as 0 and 1
233        assert len(self.openings) == 2
234                           
235        # Determine pipe direction     
236        self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation
237        if delta_z > 0.0:
238            self.inlet = self.openings[0]
239            self.outlet = self.openings[1]
240        else:               
241            self.outlet = self.openings[0]
242            self.inlet = self.openings[1]       
243       
244       
245        # Store basic geometry
246        self.end_points = [end_point0, end_point1]
247        self.vector = P['vector']
248        self.length = P['length']; assert self.length > 0.0
249        if not allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
250            msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length)
251            msg += ' does not match distance between specified'
252            msg += ' end points (%.2f m)' %self.length
253            print msg
254       
255        self.verbose = verbose
256        self.last_update = 0.0 # For use with update_interval       
257        self.update_interval = update_interval
258       
259
260        # Print something
261        if hasattr(self, 'log_filename'):
262            s = 'Culvert Effective Length = %.2f m' %(self.length)
263            log_to_file(self.log_filename, s)
264   
265            s = 'Culvert Direction is %s\n' %str(self.vector)
266            log_to_file(self.log_filename, s)       
267       
268       
269       
270       
271       
272    def __call__(self, domain):
273        from anuga.utilities.numerical_tools import mean
274       
275        from anuga.config import g, epsilon
276        from Numeric import take, sqrt
277        from anuga.config import velocity_protection       
278
279
280        # Time stuff
281        time = domain.get_time()
282       
283       
284        update = False
285        if self.update_interval is None:
286            update = True
287        else:   
288            if time - self.last_update > self.update_interval or time == 0.0:
289                update = True
290
291                               
292        if update is True:
293            self.last_update = time
294            dq = domain.quantities
295                       
296            # Get average water depths at each opening       
297            openings = self.openings   # There are two Opening [0] and [1]
298            for i, opening in enumerate(openings):
299               
300                # Compute mean values of selected quantitites in the
301                # enquiry area in front of the culvert
302                # Stage and velocity comes from enquiry area
303                # and elevation from exchange area
304               
305                stage = dq['stage'].get_values(location='centroids',
306                                               indices=[self.enquiry_indices[i]])
307
308                   
309                # Store current average stage and depth with each opening object
310                opening.depth = stage - opening.elevation
311                opening.stage = stage
312
313               
314
315            #################  End of the FOR loop ################################################
316
317            # We now need to deal with each opening individually
318               
319            # Determine flow direction based on total energy difference
320
321            delta_w = self.inlet.stage - self.outlet.stage
322           
323            if hasattr(self, 'log_filename'):
324                s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time, 
325                                                                           self.inlet.stage,
326                                                                           self.outlet.stage)
327                log_to_file(self.log_filename, s)
328
329
330            if self.inlet.depth <= 0.01:
331                Q = 0.0
332            else:
333                # Calculate discharge for one barrel and set inlet.rate and outlet.rate
334                try:
335                    Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1]) 
336                except Below_interval, e:
337                    Q = self.rating_curve[0,1]             
338                    msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
339                    msg += str(e)
340                    msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
341                        %(Q, self.label)
342                    if hasattr(self, 'log_filename'):                   
343                        log_to_file(self.log_filename, msg)
344                except Above_interval, e:
345                    Q = self.rating_curve[-1,1]             
346                    msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
347                    msg += str(e)
348                    msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
349                        %(Q, self.label)
350                    if hasattr(self, 'log_filename'):                   
351                        log_to_file(self.log_filename, msg)                   
352
353               
354               
355           
356            # Adjust discharge for multiple barrels
357            Q *= self.number_of_barrels
358           
359           
360            self.inlet.rate = -Q
361            self.outlet.rate = Q
362
363            # Log timeseries to file
364            try:
365                fid = open(self.timeseries_filename, 'a')       
366            except:
367                pass
368            else:   
369                fid.write('%.2f, %.2f\n' %(time, Q))
370                fid.close()
371
372
373
374           
375   
376        # Execute flow term for each opening
377        # This is where Inflow objects are evaluated using the last rate that has been calculated
378        #
379        # This will take place at every internal timestep and update the domain
380        for opening in self.openings:
381            opening(domain)
382
383
384
385       
386       
387
388class Culvert_flow_energy:
389    """Culvert flow - transfer water from one hole to another
390   
391    Using Momentum as Calculated by Culvert Flow !!
392    Could be Several Methods Investigated to do This !!!
393
394    2008_May_08
395    To Ole:
396    OK so here we need to get the Polygon Creating code to create
397    polygons for the culvert Based on
398    the two input Points (X0,Y0) and (X1,Y1) - need to be passed
399    to create polygon
400
401    The two centers are now passed on to create_polygon.
402   
403
404    Input: Two points, pipe_size (either diameter or width, height),
405    mannings_rougness,
406    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
407    top-down_blockage_factor and bottom_up_blockage_factor
408   
409   
410    And the Delta H enquiry should be change from Openings in line 412
411    to the enquiry Polygons infront of the culvert
412    At the moment this script uses only Depth, later we can change it to
413    Energy...
414
415    Once we have Delta H can calculate a Flow Rate and from Flow Rate
416    an Outlet Velocity
417    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
418       
419    Invert levels are optional. If left out they will default to the
420    elevation at the opening.
421       
422    """ 
423
424    def __init__(self,
425                 domain,
426                 label=None, 
427                 description=None,
428                 end_point0=None, 
429                 end_point1=None,
430                 width=None,
431                 height=None,
432                 diameter=None,
433                 manning=None,          # Mannings Roughness for Culvert
434                 invert_level0=None,    # Invert level at opening 0
435                 invert_level1=None,    # Invert level at opening 1
436                 loss_exit=None,
437                 loss_entry=None,
438                 loss_bend=None,
439                 loss_special=None,
440                 blockage_topdwn=None,
441                 blockage_bottup=None,
442                 culvert_routine=None,
443                 number_of_barrels=1,
444                 update_interval=None,
445                 verbose=False):
446       
447        from Numeric import sqrt, sum
448
449        # Input check
450        if diameter is not None:
451            self.culvert_type = 'circle'
452            self.diameter = diameter
453            if height is not None or width is not None:
454                msg = 'Either diameter or width&height must be specified, '
455                msg += 'but not both.'
456                raise Exception, msg
457        else:
458            if height is not None:
459                if width is None:
460                    self.culvert_type = 'square'                               
461                    width = height
462                else:
463                    self.culvert_type = 'rectangle'
464            elif width is not None:
465                if height is None:
466                    self.culvert_type = 'square'                               
467                    height = width
468            else:
469                msg = 'Either diameter or width&height must be specified.'
470                raise Exception, msg               
471               
472            if height == width:
473                self.culvert_type = 'square'                                               
474               
475            self.height = height
476            self.width = width
477
478           
479        assert self.culvert_type in ['circle', 'square', 'rectangle']
480       
481        assert number_of_barrels >= 1
482        self.number_of_barrels = number_of_barrels
483       
484       
485        # Set defaults
486        if manning is None: manning = 0.012   # Default roughness for pipe
487        if loss_exit is None: loss_exit = 1.00
488        if loss_entry is None: loss_entry = 0.50
489        if loss_bend is None: loss_bend=0.00
490        if loss_special is None: loss_special=0.00
491        if blockage_topdwn is None: blockage_topdwn=0.00
492        if blockage_bottup is None: blockage_bottup=0.00
493        if culvert_routine is None: 
494            culvert_routine=boyd_generalised_culvert_model
495           
496        if label is None: label = 'culvert_flow'
497        label += '_' + str(id(self)) 
498        self.label = label
499       
500        # File for storing culvert quantities
501        self.timeseries_filename = label + '_timeseries.csv'
502        fid = open(self.timeseries_filename, 'w')
503        fid.write('time, E0, E1, Velocity, Discharge\n')
504        fid.close()
505
506        # Log file for storing general textual output
507        self.log_filename = label + '.log'         
508        log_to_file(self.log_filename, self.label)       
509        log_to_file(self.log_filename, description)
510        log_to_file(self.log_filename, self.culvert_type)       
511
512
513        # Create the fundamental culvert polygons from POLYGON
514        if self.culvert_type == 'circle':
515            # Redefine width and height for use with create_culvert_polygons
516            width = height = diameter
517       
518        P = create_culvert_polygons(end_point0,
519                                    end_point1,
520                                    width=width,   
521                                    height=height,
522                                    number_of_barrels=number_of_barrels)
523       
524        if verbose is True:
525            pass
526            #plot_polygons([[end_point0, end_point1],
527            #               P['exchange_polygon0'],
528            #               P['exchange_polygon1'],
529            #               P['enquiry_polygon0'],
530            #               P['enquiry_polygon1']],
531            #              figname='culvert_polygon_output')
532            #import sys; sys.exit()                           
533
534           
535        # Compute the average point for enquiry
536        enquiry_point0 = sum(P['enquiry_polygon0'][:2])/2 
537        enquiry_point1 = sum(P['enquiry_polygon1'][:2])/2         
538       
539        self.enquiry_points = [enquiry_point0, enquiry_point1]                           
540        self.enquiry_indices = []                 
541        for point in self.enquiry_points:
542            # Find nearest centroid
543            N = len(domain)   
544            points = domain.get_centroid_coordinates(absolute=True)
545
546            # Calculate indices in exchange area for this forcing term
547           
548            triangle_id = min_dist = sys.maxint
549            for k in range(N):
550                x, y = points[k,:] # Centroid
551
552                c = point                               
553                distance = (x-c[0])**2+(y-c[1])**2
554                if distance < min_dist:
555                    min_dist = distance
556                    triangle_id = k
557
558                   
559            if triangle_id < sys.maxint:
560                msg = 'found triangle with centroid (%f, %f)'\
561                    %tuple(points[triangle_id, :])
562                msg += ' for point (%f, %f)' %tuple(point)
563               
564                self.enquiry_indices.append(triangle_id)
565            else:
566                msg = 'Triangle not found for point (%f, %f)' %point
567                raise Exception, msg
568       
569                         
570
571           
572           
573
574        # Check that all polygons lie within the mesh.
575        bounding_polygon = domain.get_boundary_polygon()
576        for key in P.keys():
577            if key in ['exchange_polygon0', 
578                       'exchange_polygon1',
579                       'enquiry_polygon0',
580                       'enquiry_polygon1']:
581                for point in P[key]:
582               
583                    msg = 'Point %s in polygon %s for culvert %s did not'\
584                        %(str(point), key, self.label)
585                    msg += 'fall within the domain boundary.'
586                    assert is_inside_polygon(point, bounding_polygon), msg
587       
588
589        # Create inflow object at each end of the culvert.
590        self.openings = []
591        self.openings.append(Inflow(domain,
592                                    polygon=P['exchange_polygon0']))
593
594        self.openings.append(Inflow(domain,
595                                    polygon=P['exchange_polygon1']))                                   
596
597
598        # Assume two openings for now: Referred to as 0 and 1
599        assert len(self.openings) == 2
600       
601        # Store basic geometry
602        self.end_points = [end_point0, end_point1]
603        self.invert_levels = [invert_level0, invert_level1]               
604        #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
605        self.enquiry_polylines = [P['enquiry_polygon0'][:2], 
606                                  P['enquiry_polygon1'][:2]]
607        self.vector = P['vector']
608        self.length = P['length']; assert self.length > 0.0
609        self.verbose = verbose
610        self.last_time = 0.0       
611        self.last_update = 0.0 # For use with update_interval       
612        self.update_interval = update_interval
613       
614
615        # Store hydraulic parameters
616        self.manning = manning
617        self.loss_exit = loss_exit
618        self.loss_entry = loss_entry
619        self.loss_bend = loss_bend
620        self.loss_special = loss_special
621        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
622        self.blockage_topdwn = blockage_topdwn
623        self.blockage_bottup = blockage_bottup
624
625        # Store culvert routine
626        self.culvert_routine = culvert_routine
627
628       
629        # Create objects to update momentum (a bit crude at this stage)
630
631       
632        xmom0 = General_forcing(domain, 'xmomentum',
633                                polygon=P['exchange_polygon0'])
634
635        xmom1 = General_forcing(domain, 'xmomentum',
636                                polygon=P['exchange_polygon1'])
637
638        ymom0 = General_forcing(domain, 'ymomentum',
639                                polygon=P['exchange_polygon0'])
640
641        ymom1 = General_forcing(domain, 'ymomentum',
642                                polygon=P['exchange_polygon1'])
643
644        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
645       
646
647        # Print something
648        s = 'Culvert Effective Length = %.2f m' %(self.length)
649        log_to_file(self.log_filename, s)
650   
651        s = 'Culvert Direction is %s\n' %str(self.vector)
652        log_to_file(self.log_filename, s)       
653       
654       
655    def __call__(self, domain):
656        from anuga.utilities.numerical_tools import mean
657       
658        from anuga.config import g, epsilon
659        from Numeric import take, sqrt
660        from anuga.config import velocity_protection       
661
662
663        log_filename = self.log_filename
664         
665        # Time stuff
666        time = domain.get_time()
667       
668       
669        update = False
670        if self.update_interval is None:
671            update = True
672        else:   
673            if time - self.last_update > self.update_interval or time == 0.0:
674                update = True
675
676        #print 'call', time, time - self.last_update
677               
678                               
679        if update is True:
680            #print 'Updating', time, time - self.last_update
681            self.last_update = time
682       
683            delta_t = time-self.last_time
684            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
685            log_to_file(log_filename, s)
686       
687            msg = 'Time did not advance'
688            if time > 0.0: assert delta_t > 0.0, msg
689
690
691            # Get average water depths at each opening       
692            openings = self.openings   # There are two Opening [0] and [1]
693            for i, opening in enumerate(openings):
694                dq = domain.quantities
695               
696                # Compute mean values of selected quantitites in the
697                # exchange area in front of the culvert
698                # Stage and velocity comes from enquiry area
699                # and elevation from exchange area
700               
701                stage = dq['stage'].get_values(location='centroids',
702                                               indices=opening.exchange_indices)           
703                w = mean(stage) # Average stage
704
705                # Use invert level instead of elevation if specified
706                invert_level = self.invert_levels[i]
707                if invert_level is not None:
708                    z = invert_level
709                else:
710                    elevation = dq['elevation'].get_values(location='centroids', 
711                                                           indices=opening.exchange_indices)
712                    z = mean(elevation) # Average elevation
713
714                # Estimated depth above the culvert inlet
715                d = w - z  # Used for calculations involving depth
716                if d < 0.0:
717                    # This is possible since w and z are taken at different locations
718                    #msg = 'D < 0.0: %f' %d
719                    #raise msg
720                    d = 0.0
721               
722
723                # Ratio of depth to culvert height.
724                # If ratio > 1 then culvert is running full
725                if self.culvert_type == 'circle':
726                    ratio = d/self.diameter
727                else:   
728                    ratio = d/self.height 
729                opening.ratio = ratio
730                   
731                   
732                # Average measures of energy in front of this opening
733                #polyline = self.enquiry_polylines[i]
734                #opening.total_energy = domain.get_energy_through_cross_section(polyline,
735                #                                                               kind='total')           
736               
737                id = [self.enquiry_indices[i]]
738                stage = dq['stage'].get_values(location='centroids',
739                                               indices=id)
740                elevation = dq['elevation'].get_values(location='centroids',
741                                               indices=id)                                               
742                xmomentum = dq['xmomentum'].get_values(location='centroids',
743                                               indices=id)                                               
744                ymomentum = dq['xmomentum'].get_values(location='centroids',
745                                                       indices=id)                                                                                             
746                depth = stage-elevation
747                if depth > 0.0:
748                    u = xmomentum/(depth + velocity_protection/depth)
749                    v = ymomentum/(depth + velocity_protection/depth)
750                else:
751                    u = v = 0.0
752
753                   
754                opening.total_energy = 0.5*(u*u + v*v)/g + stage
755                # FIXME(Ole): What happens if this is negative?
756                #print 'Et = %.3f m' %opening.total_energy
757
758                # Store current average stage and depth with each opening object
759                opening.depth = d
760                opening.depth_trigger = d # Use this for now
761                opening.stage = w
762                opening.elevation = z
763               
764
765            #################  End of the FOR loop ################################################
766
767            # We now need to deal with each opening individually
768               
769            # Determine flow direction based on total energy difference
770            delta_Et = openings[0].total_energy - openings[1].total_energy
771
772            if delta_Et > 0:
773                #print 'Flow U/S ---> D/S'
774                inlet = openings[0]
775                outlet = openings[1]
776
777                inlet.momentum = self.opening_momentum[0]
778                outlet.momentum = self.opening_momentum[1]
779
780            else:
781                #print 'Flow D/S ---> U/S'
782                inlet = openings[1]
783                outlet = openings[0]
784
785                inlet.momentum = self.opening_momentum[1]
786                outlet.momentum = self.opening_momentum[0]
787               
788                delta_Et = -delta_Et
789
790            self.inlet = inlet
791            self.outlet = outlet
792               
793            msg = 'Total energy difference is negative'
794            assert delta_Et > 0.0, msg
795
796            delta_h = inlet.stage - outlet.stage
797            delta_z = inlet.elevation - outlet.elevation
798            culvert_slope = (delta_z/self.length)
799
800            if culvert_slope < 0.0:
801                # Adverse gradient - flow is running uphill
802                # Flow will be purely controlled by uphill outlet face
803                if self.verbose is True:
804                    print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
805
806
807            s = 'Delta total energy = %.3f' %(delta_Et)
808            log_to_file(log_filename, s)
809
810
811            # Calculate discharge for one barrel and set inlet.rate and outlet.rate
812            Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
813       
814            # Adjust discharge for multiple barrels
815            Q *= self.number_of_barrels
816
817            # Compute barrel momentum
818            barrel_momentum = barrel_velocity*culvert_outlet_depth
819                   
820            s = 'Barrel velocity = %f' %barrel_velocity
821            log_to_file(log_filename, s)
822
823            # Compute momentum vector at outlet
824            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
825               
826            s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
827            log_to_file(log_filename, s)
828
829            # Log timeseries to file
830            fid = open(self.timeseries_filename, 'a')       
831            fid.write('%f, %f, %f, %f, %f\n'\
832                          %(time, 
833                            openings[0].total_energy,
834                            openings[1].total_energy,
835                            barrel_velocity,
836                            Q))
837            fid.close()
838
839            # Update momentum       
840            delta_t = time - self.last_time
841            if delta_t > 0.0:
842                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
843                xmomentum_rate /= delta_t
844                   
845                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
846                ymomentum_rate /= delta_t
847                       
848                s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
849                log_to_file(log_filename, s)                   
850            else:
851                xmomentum_rate = ymomentum_rate = 0.0
852
853
854            # Set momentum rates for outlet jet
855            outlet.momentum[0].rate = xmomentum_rate
856            outlet.momentum[1].rate = ymomentum_rate
857
858            # Remember this value for next step (IMPORTANT)
859            outlet.momentum[0].value = outlet_mom_x
860            outlet.momentum[1].value = outlet_mom_y                   
861
862            if int(domain.time*100) % 100 == 0:
863                s = 'T=%.5f, Culvert Discharge = %.3f f'\
864                    %(time, Q)
865                s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
866                     %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
867                s += ' Momentum rate: (%.4f, %.4f)'\
868                     %(xmomentum_rate, ymomentum_rate)                   
869                s+='Outlet Vel= %.3f'\
870                    %(barrel_velocity)
871                log_to_file(log_filename, s)
872           
873               
874
875
876        # Execute flow term for each opening
877        # This is where Inflow objects are evaluated and update the domain
878        for opening in self.openings:
879            opening(domain)
880
881        # Execute momentum terms
882        # This is where Inflow objects are evaluated and update the domain
883        self.outlet.momentum[0](domain)
884        self.outlet.momentum[1](domain)       
885           
886        # Store value of time #FIXME(Ole): Maybe only every time we update   
887        self.last_time = time
888
889
890Culvert_flow = Culvert_flow_rating       
Note: See TracBrowser for help on using the repository browser.