source: trunk/anuga_core/source/anuga/structures/culvert_operator.py @ 7955

Last change on this file since 7955 was 7955, checked in by steve, 14 years ago

Changing culvert forcing term class to a culvert fractional step operator

File size: 29.2 KB
Line 
1import sys
2
3from anuga.shallow_water.forcing import Inflow, General_forcing
4from anuga.structures.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
7from anuga.geometry.polygon import plot_polygons, polygon_area
8
9
10from anuga.utilities.numerical_tools import mean
11from anuga.utilities.numerical_tools import ensure_numeric, sign
12       
13from anuga.config import g, epsilon
14from anuga.config import minimum_allowed_height, velocity_protection       
15import anuga.utilities.log as log
16
17import numpy as num
18from math import sqrt
19from math import sqrt
20
21class Below_interval(Exception): pass 
22class Above_interval(Exception): pass
23
24# FIXME(Ole): Take a good hard look at logging here
25
26
27# FIXME(Ole): Write in C and reuse this function by similar code
28# in interpolate.py
29def interpolate_linearly(x, xvec, yvec):
30
31    msg = 'Input to function interpolate_linearly could not be converted '
32    msg += 'to numerical scalar: x = %s' % str(x)
33    try:
34        x = float(x)
35    except:
36        raise Exception, msg
37
38
39    # Check bounds
40    if x < xvec[0]: 
41        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
42            % (x, xvec[0])
43        raise Below_interval, msg
44       
45    if x > xvec[-1]: 
46        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
47            %(x, xvec[-1])
48        raise Above_interval, msg       
49       
50       
51    # Find appropriate slot within bounds           
52    i = 0
53    while x > xvec[i]: i += 1
54
55   
56    x0 = xvec[i-1]
57    x1 = xvec[i]           
58    alpha = (x - x0)/(x1 - x0) 
59           
60    y0 = yvec[i-1]
61    y1 = yvec[i]                       
62    y = alpha*y1 + (1-alpha)*y0
63   
64    return y
65           
66
67           
68def read_culvert_description(culvert_description_filename):           
69   
70    # Read description file
71    fid = open(culvert_description_filename)
72   
73    read_rating_curve_data = False       
74    rating_curve = []
75    for i, line in enumerate(fid.readlines()):
76       
77        if read_rating_curve_data is True:
78            fields = line.split(',')
79            head_difference = float(fields[0].strip())
80            flow_rate = float(fields[1].strip())               
81            barrel_velocity = float(fields[2].strip())
82           
83            rating_curve.append([head_difference, flow_rate, barrel_velocity]) 
84       
85        if i == 0:
86            # Header
87            continue
88        if i == 1:
89            # Metadata
90            fields = line.split(',')
91            label=fields[0].strip()
92            type=fields[1].strip().lower()
93            assert type in ['box', 'pipe']
94           
95            width=float(fields[2].strip())
96            height=float(fields[3].strip())
97            length=float(fields[4].strip())
98            number_of_barrels=int(fields[5].strip())
99            #fields[6] refers to losses
100            description=fields[7].strip()               
101               
102        if line.strip() == '': continue # Skip blanks
103
104        if line.startswith('Rating'):
105            read_rating_curve_data = True
106            # Flow data follows
107               
108    fid.close()
109   
110    return label, type, width, height, length, number_of_barrels, description, rating_curve
111   
112
113   
114
115class Generic_box_culvert:
116    """Culvert flow - transfer water from one rectngular box to another.
117    Sets up the geometry of problem
118   
119    This is the base class for culverts. Inherit from this class (and overwrite
120    compute_discharge method for specific subclasses)
121   
122    Input: Two points, pipe_size (either diameter or width, height),
123    mannings_rougness,
124    """ 
125
126    def __init__(self,
127                 domain,
128                 end_point0=None, 
129                 end_point1=None,
130                 enquiry_gap_factor=0.2,
131                 width=None,
132                 height=None,
133                 verbose=False):
134       
135        # Input check
136       
137        if height is None:
138            height = width
139
140        self.height = height
141        self.width = width
142       
143        self.domain = domain
144        self.end_points= [end_point0, end_point1]
145        self.enquiry_gap_factor=enquiry_gap_factor
146        self.verbose=verbose
147       
148        self.filename = None
149       
150
151        # Create the fundamental culvert polygons from polygon
152        self.create_culvert_polygons()
153        self.compute_enquiry_indices()
154        self.check_culvert_inside_domain()
155
156
157        # Establish initial values at each enquiry point
158        dq = domain.quantities
159        for i, opening in enumerate(self.openings):
160            idx = self.enquiry_indices[i]
161            elevation = dq['elevation'].get_values(location='centroids',
162                                                   indices=[idx])[0]
163            stage = dq['stage'].get_values(location='centroids',
164                                           indices=[idx])[0]
165            opening.elevation = elevation
166            opening.stage = stage
167            opening.depth = stage-elevation
168
169           
170                           
171        # Determine initial pipe direction.
172        # This may change dynamically based on the total energy difference     
173        # Consequently, this may be superfluous
174        delta_z = self.openings[0].elevation - self.openings[1].elevation
175        if delta_z > 0.0:
176            self.inlet = self.openings[0]
177            self.outlet = self.openings[1]
178        else:               
179            self.outlet = self.openings[0]
180            self.inlet = self.openings[1]       
181       
182       
183        # Store basic geometry
184        self.end_points = [end_point0, end_point1]
185        self.vector = P['vector']
186        self.length = P['length']; assert self.length > 0.0
187        if culvert_description_filename is not None:
188            if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
189                msg = 'WARNING: barrel length specified in "%s" (%.2f m)'\
190                    % (culvert_description_filename, 
191                       length)
192                msg += ' does not match distance between specified'
193                msg += ' end points (%.2f m)' %self.length
194                log.critical(msg)
195       
196        self.verbose = verbose
197
198        # Circular index for flow averaging in culvert
199        self.N = N = number_of_smoothing_steps
200        self.Q_list = [0]*N
201        self.i = i
202       
203        # For use with update_interval                       
204        self.last_update = 0.0
205        self.update_interval = update_interval
206       
207        # Create objects to update momentum (a bit crude at this stage). This is used with the momentum jet.
208        xmom0 = General_forcing(domain, 'xmomentum',
209                                polygon=P['exchange_polygon0'])
210
211        xmom1 = General_forcing(domain, 'xmomentum',
212                                polygon=P['exchange_polygon1'])
213
214        ymom0 = General_forcing(domain, 'ymomentum',
215                                polygon=P['exchange_polygon0'])
216
217        ymom1 = General_forcing(domain, 'ymomentum',
218                                polygon=P['exchange_polygon1'])
219
220        self.opening_momentum = [[xmom0, ymom0], [xmom1, ymom1]]
221
222
223
224        # Print some diagnostics to log if requested
225        if self.log_filename is not None:
226            s = 'Culvert Effective Length = %.2f m' %(self.length)
227            log_to_file(self.log_filename, s)
228   
229            s = 'Culvert Direction is %s\n' %str(self.vector)
230            log_to_file(self.log_filename, s)       
231
232
233    def set_store_hydrograph_discharge(self,filename=None):
234
235        if filename is None:
236            self.filename = 'culvert_discharge_hydrograph'
237        else:
238            self.filename = filename
239
240        self.discharge_hydrograph = True
241       
242        self.timeseries_filename = self.filename + '_timeseries.csv'
243        fid = open(self.timeseries_filename, 'w')
244        fid.write('time, discharge\n')
245        fid.close()
246
247
248
249
250    def create_culvert_polygons(self):
251        """Create polygons at the end of a culvert inlet and outlet.
252        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
253        for enquiring the total energy at both ends of the culvert and transferring flow.
254        """
255
256        # Calculate geometry
257        x0, y0 = self.end_point0
258        x1, y1 = self.end_point1
259
260        dx = x1-x0
261        dy = y1-y0
262
263        self.culvert_vector = num.array([dx, dy])
264        self.culvert_length = sqrt(num.sum(self.culvert_vector**2))
265
266
267        # Unit direction vector and normal
268        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
269        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
270
271        # Short hands
272        w = 0.5*width*normal # Perpendicular vector of 1/2 width
273        h = height*vector    # Vector of length=height in the
274                             # direction of the culvert
275        gap = (1 + enquiry_gap_factor)*h
276
277
278        self.exchange_polygons = []
279
280        # Build exchange polygon and enquiry point for opening 0
281        p0 = self.end_point0 + w
282        p1 = self.end_point0 - w
283        p2 = p1 - h
284        p3 = p0 - h
285        self.exchange_polygons.append(num.array([p0,p1,p2,p3]))
286        self.enquiry_points= end_point0 - gap
287
288
289        # Build exchange polygon and enquiry point for opening 1
290        p0 = end_point1 + w
291        p1 = end_point1 - w
292        p2 = p1 + h
293        p3 = p0 + h
294        self.exchange_polygon1 = num.array([p0,p1,p2,p3])
295        self.enquiry_point1 = end_point1 + gap
296
297        # Check that enquiry point0 is outside exchange polygon0
298        polygon = self.exchange_polygon0
299        area = polygon_area(polygon)
300
301            msg = 'Polygon %s ' %(polygon)
302            msg += ' has area = %f' % area
303            assert area > 0.0, msg
304
305            for key2 in ['enquiry_point0', 'enquiry_point1']:
306                point = culvert_polygons[key2]
307                msg = 'Enquiry point falls inside an enquiry point.'
308
309                assert not inside_polygon(point, polygon), msg
310
311
312
313        for key1 in ['exchange_polygon0',
314                     'exchange_polygon1']:
315           
316
317        # Return results
318        self.culvert_polygons = culvert_polygons
319
320
321
322       
323    def __call__(self, domain):
324
325        # Time stuff
326        time = domain.get_time()
327       
328       
329        update = False
330        if self.update_interval is None:
331            # Use next timestep as has been computed in domain.py       
332            delta_t = domain.timestep           
333            update = True
334        else:   
335            # Use update interval
336            delta_t = self.update_interval           
337            if time - self.last_update > self.update_interval or time == 0.0:
338                update = True
339           
340        if self.log_filename is not None:       
341            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
342            log_to_file(self.log_filename, s)
343               
344                               
345        if update is True:
346            self.compute_rates(delta_t)
347           
348   
349        # Execute flow term for each opening
350        # This is where Inflow objects are evaluated using the last rate
351        # that has been calculated
352        #
353        # This will take place at every internal timestep and update the domain
354        for opening in self.openings:
355            opening(domain)
356           
357
358
359    def compute_enquiry_indices(self):
360        """Get indices for nearest centroids to self.enquiry_points
361        """
362       
363        domain = self.domain
364       
365        enquiry_indices = []                 
366        for point in self.enquiry_points:
367            # Find nearest centroid
368            N = len(domain)   
369            points = domain.get_centroid_coordinates(absolute=True)
370
371            # Calculate indices in exchange area for this forcing term
372           
373            triangle_id = min_dist = sys.maxint
374            for k in range(N):
375                x, y = points[k,:] # Centroid
376
377                c = point                               
378                distance = (x-c[0])**2+(y-c[1])**2
379                if distance < min_dist:
380                    min_dist = distance
381                    triangle_id = k
382
383                   
384            if triangle_id < sys.maxint:
385                msg = 'found triangle with centroid (%f, %f)'\
386                    %tuple(points[triangle_id, :])
387                msg += ' for point (%f, %f)' %tuple(point)
388               
389                enquiry_indices.append(triangle_id)
390            else:
391                msg = 'Triangle not found for point (%f, %f)' %point
392                raise Exception, msg
393       
394        self.enquiry_indices = enquiry_indices
395
396       
397    def check_culvert_inside_domain(self):
398        """Check that all polygons and enquiry points lie within the mesh.
399        """
400        bounding_polygon = self.domain.get_boundary_polygon()
401        P = self.culvert_polygons
402        for key in P.keys():
403            if key in ['exchange_polygon0', 
404                       'exchange_polygon1']:
405                for point in list(P[key]) + self.enquiry_points:
406                    msg = 'Point %s in polygon %s for culvert %s did not'\
407                        %(str(point), key, self.label)
408                    msg += 'fall within the domain boundary.'
409                    assert is_inside_polygon(point, bounding_polygon), msg
410           
411
412    def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
413        """Adjust Q downwards depending on available water at inlet
414
415           This is a critical step in modelling bridges and Culverts
416           the predicted flow through a structure based on an abstract
417           algorithm can at times request for water that is simply not
418           available due to any number of constrictions that limit the
419           flow approaching the structure In order to ensure that
420           there is adequate flow available certain checks are
421           required There needs to be a check using the Static Water
422           Volume sitting infront of the structure, In addition if the
423           water is moving the available water will be larger than the
424           static volume
425           
426           NOTE To temporarily switch this off for Debugging purposes
427           rem out line in function def compute_rates below
428        """
429   
430        if delta_t < epsilon:
431            # No need to adjust if time step is very small or zero
432            # In this case the possible flow will be very large
433            # anyway.
434            return Q
435       
436        # Short hands
437        domain = self.domain       
438        dq = domain.quantities               
439        time = domain.get_time()
440        I = self.inlet
441        idx = I.exchange_indices   
442
443        # Find triangle with the smallest depth
444        stage = dq['stage'].get_values(location='centroids', 
445                                               indices=[idx])
446        elevation = dq['elevation'].get_values(location='centroids', 
447                                               indices=[idx])       
448        depth = stage-elevation
449        min_depth = min(depth.flat)  # This may lead to errors if edge of area is at a higher level !!!!
450        avg_depth = mean(depth.flat) # Yes, but this one violates the conservation unit tests
451
452
453
454        # FIXME (Ole): If you want these, use log.critical() and
455        # make the statements depend on verbose
456        #print I.depth
457        #print I.velocity
458        #print self.width
459
460        # max_Q Based on Volume Calcs
461
462
463        depth_term = min_depth*I.exchange_area/delta_t
464        if min_depth < 0.2:
465            # Only add velocity term in shallow waters (< 20 cm)
466            # This is a little ad hoc, but maybe it is reasonable
467            velocity_term = self.width*min_depth*I.velocity
468        else:
469            velocity_term = 0.0
470
471        # This one takes approaching water into account   
472        max_Q = max(velocity_term, depth_term)
473
474        # This one preserves Volume
475        #max_Q = depth_term
476
477
478        if self.verbose is True:
479            log.critical('Max_Q = %f' % max_Q)           
480            msg = 'Width = %.2fm, Depth at inlet = %.2f m, Velocity = %.2f m/s.      ' % (self.width, I.depth, I.velocity)
481            msg += 'Max Q = %.2f m^3/s' %(max_Q)
482            log.critical(msg)
483
484        if self.log_filename is not None:               
485            log_to_file(self.log_filename, msg)
486        # New Procedure for assessing the flow available to the Culvert
487        # This routine uses the GET FLOW THROUGH CROSS SECTION
488        #   Need to check Several Polyline however
489        # Firstly 3 sides of the exchange Poly
490        # then only the Line Directly infront of the Polygon
491        # Access polygon Points from   self.inlet.polygon
492     
493        #  The Following computes the flow crossing over 3 sides of the exchange polygon for the structure
494        # Clearly the flow in the culvert can not be more than that flowing toward it through the exhange polygon
495       
496        #q1 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][1:3]) # First Side Segment
497        #q2 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][2:])   # Second Face Segment
498        #q3 =domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'].take([3,0], axis=0)) # Third Side Segment
499        # q4 = domain.get_flow_through_cross_section([self.culvert_polygons['exchange_polygon0'][1:4]][0])
500        #q4=max(q1,0.0)+max(q2,0.0)+max(q3,0.0)
501        # To use only the Flow crossing the 3 sides of the Exchange Polygon use the following Line Only
502        #max_Q=max(q1,q2,q3,q4)
503        # Try Simple Smoothing using Average of 2 approaches
504        #max_Q=(max(q1,q2,q3,q4)+max_Q)/2.0
505        # Calculate the minimum in absolute terms of
506        # the requsted flow and the possible flow
507        Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
508        if self.verbose is True:
509            msg = 'Initial Q Reduced = %.2f m3/s.      ' % (Q_reduced)
510            log.critical(msg)
511
512        if self.log_filename is not None:               
513            log_to_file(self.log_filename, msg)
514        # Now Keep Rolling Average of Computed Discharge to Reduce / Remove Oscillations
515        #  can use delta_t if we want to averageover a time frame for example
516        # N = 5.0/delta_t  Will provide the average over 5 seconds
517
518        self.i=(self.i+1)%self.N
519        self.Q_list[self.i]=Q_reduced
520        Q_reduced = sum(self.Q_list)/len(self.Q_list)
521
522        if self.verbose is True:
523            msg = 'Final Q Reduced = %.2f m3/s.      ' % (Q_reduced)
524            log.critical(msg)
525
526        if self.log_filename is not None:               
527            log_to_file(self.log_filename, msg)
528
529
530        if abs(Q_reduced) < abs(Q): 
531            msg = '%.2fs: Requested flow is ' % time
532            msg += 'greater than what is supported by the smallest '
533            msg += 'depth at inlet exchange area:\n        '
534            msg += 'inlet exchange area: %.2f '% (I.exchange_area) 
535            msg += 'velocity at inlet :%.2f '% (I.velocity)
536            msg += 'Vel* Exch Area = : %.2f '% (I.velocity*avg_depth*self.width)
537            msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\
538                % (avg_depth, I.exchange_area, delta_t)
539            msg += ' = %.2f m^3/s\n        ' % Q_reduced
540            msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
541            msg += 'Note calculate max_Q from V %.2f m^3/s ' % (max_Q)
542            if self.verbose is True:
543                log.critical(msg)
544               
545            if self.log_filename is not None:               
546                log_to_file(self.log_filename, msg)
547       
548        return Q_reduced   
549                       
550           
551    def compute_rates(self, delta_t):
552        """Compute new rates for inlet and outlet
553        """
554
555        # Short hands
556        domain = self.domain       
557        dq = domain.quantities               
558       
559        # Time stuff
560        time = domain.get_time()
561        self.last_update = time
562
563           
564        if hasattr(self, 'log_filename'):
565            log_filename = self.log_filename
566           
567        # Compute stage, energy and velocity at the
568        # enquiry points at each end of the culvert
569        openings = self.openings
570        for i, opening in enumerate(openings):
571            idx = self.enquiry_indices[i]               
572           
573            stage = dq['stage'].get_values(location='centroids',
574                                           indices=[idx])[0]
575            depth = h = stage-opening.elevation
576                                                           
577           
578            # Get velocity                                 
579            xmomentum = dq['xmomentum'].get_values(location='centroids',
580                                                   indices=[idx])[0]
581            ymomentum = dq['xmomentum'].get_values(location='centroids',
582                                                   indices=[idx])[0]
583
584            if h > minimum_allowed_height:
585                u = xmomentum/(h + velocity_protection/h)
586                v = ymomentum/(h + velocity_protection/h)
587            else:
588                u = v = 0.0
589               
590            v_squared = u*u + v*v
591           
592            if self.use_velocity_head is True:
593                velocity_head = 0.5*v_squared/g   
594            else:
595                velocity_head = 0.0
596           
597            opening.total_energy = velocity_head + stage
598            opening.specific_energy = velocity_head + depth
599            opening.stage = stage
600            opening.depth = depth
601            opening.velocity = sqrt(v_squared)
602           
603
604        # We now need to deal with each opening individually
605        # Determine flow direction based on total energy difference
606        delta_total_energy = openings[0].total_energy - openings[1].total_energy
607        if delta_total_energy > 0:
608            inlet = openings[0]
609            outlet = openings[1]
610
611            # FIXME: I think this whole momentum jet thing could be a bit more elegant
612            inlet.momentum = self.opening_momentum[0]
613            outlet.momentum = self.opening_momentum[1]
614        else:
615            inlet = openings[1]
616            outlet = openings[0]
617           
618            inlet.momentum = self.opening_momentum[1]
619            outlet.momentum = self.opening_momentum[0]
620
621            delta_total_energy = -delta_total_energy
622
623        self.inlet = inlet
624        self.outlet = outlet
625           
626        msg = 'Total energy difference is negative'
627        assert delta_total_energy >= 0.0, msg
628
629        # Recompute slope and issue warning if flow is uphill
630        # These values do not enter the computation
631        delta_z = inlet.elevation - outlet.elevation
632        culvert_slope = (delta_z/self.length)
633        if culvert_slope < 0.0:
634            # Adverse gradient - flow is running uphill
635            # Flow will be purely controlled by uphill outlet face
636            if self.verbose is True:
637                log.critical('%.2fs - WARNING: Flow is running uphill.' % time)
638           
639        if self.log_filename is not None:
640            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
641                %(time, self.inlet.stage, self.outlet.stage)
642            log_to_file(self.log_filename, s)
643            s = 'Delta total energy = %.3f' %(delta_total_energy)
644            log_to_file(log_filename, s)
645
646           
647        # Determine controlling energy (driving head) for culvert
648        if inlet.specific_energy > delta_total_energy:
649            # Outlet control
650            driving_head = delta_total_energy
651        else:
652            # Inlet control
653            driving_head = inlet.specific_energy
654           
655
656
657        if self.inlet.depth <= self.trigger_depth:
658            Q = 0.0
659        else:
660            # Calculate discharge for one barrel and
661            # set inlet.rate and outlet.rate
662           
663            if self.culvert_description_filename is not None:
664                try:
665                    Q = interpolate_linearly(driving_head, 
666                                             self.rating_curve[:,0], 
667                                             self.rating_curve[:,1]) 
668                except Below_interval, e:
669                    Q = self.rating_curve[0,1]             
670                    msg = '%.2fs: ' % time
671                    msg += 'Delta head smaller than rating curve minimum: '
672                    msg += str(e)
673                    msg += '\n        '
674                    msg += 'I will use minimum discharge %.2f m^3/s ' % Q
675                    msg += 'for culvert "%s"' % self.label
676                   
677                    if hasattr(self, 'log_filename'):                   
678                        log_to_file(self.log_filename, msg)
679                except Above_interval, e:
680                    Q = self.rating_curve[-1,1]         
681                    msg = '%.2fs: ' % time                 
682                    msg += 'Delta head greater than rating curve maximum: '
683                    msg += str(e)
684                    msg += '\n        '
685                    msg += 'I will use maximum discharge %.2f m^3/s ' % Q
686                    msg += 'for culvert "%s"' % self.label
687                   
688                    if self.log_filename is not None:                   
689                        log_to_file(self.log_filename, msg)
690            else:
691                # User culvert routine
692                Q, barrel_velocity, culvert_outlet_depth =\
693                    self.culvert_routine(inlet.depth,
694                                         outlet.depth,
695                                         inlet.velocity,
696                                         outlet.velocity,
697                                         inlet.specific_energy, 
698                                         delta_total_energy, 
699                                         g,
700                                         culvert_length=self.length,
701                                         culvert_width=self.width,
702                                         culvert_height=self.height,
703                                         culvert_type=self.culvert_type,
704                                         manning=self.manning,
705                                         sum_loss=self.sum_loss,
706                                         log_filename=self.log_filename)
707           
708           
709       
710        # Adjust discharge for multiple barrels
711        Q *= self.number_of_barrels
712
713        # Adjust discharge for available water at the inlet
714        Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
715       
716        self.inlet.rate = -Q
717        self.outlet.rate = Q
718
719
720        # Momentum jet stuff
721        if self.use_momentum_jet is True:
722
723
724            # Compute barrel momentum
725            barrel_momentum = barrel_velocity*culvert_outlet_depth
726
727            if self.log_filename is not None:                                   
728                s = 'Barrel velocity = %f' %barrel_velocity
729                log_to_file(self.log_filename, s)
730
731            # Compute momentum vector at outlet
732            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
733               
734            if self.log_filename is not None:               
735                s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
736                log_to_file(self.log_filename, s)
737
738
739            # Update momentum       
740            if delta_t > 0.0:
741                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
742                xmomentum_rate /= delta_t
743                   
744                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
745                ymomentum_rate /= delta_t
746                       
747                if self.log_filename is not None:               
748                    s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
749                    log_to_file(self.log_filename, s)                   
750            else:
751                xmomentum_rate = ymomentum_rate = 0.0
752
753
754            # Set momentum rates for outlet jet
755            outlet.momentum[0].rate = xmomentum_rate
756            outlet.momentum[1].rate = ymomentum_rate
757
758            # Remember this value for next step (IMPORTANT)
759            outlet.momentum[0].value = outlet_mom_x
760            outlet.momentum[1].value = outlet_mom_y                   
761
762            if int(domain.time*100) % 100 == 0:
763
764                if self.log_filename is not None:               
765                    s = 'T=%.5f, Culvert Discharge = %.3f f'\
766                        %(time, Q)
767                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
768                        %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
769                    s += ' Momentum rate: (%.4f, %.4f)'\
770                        %(xmomentum_rate, ymomentum_rate)                   
771                    s+='Outlet Vel= %.3f'\
772                        %(barrel_velocity)
773                    log_to_file(self.log_filename, s)
774
775
776            # Execute momentum terms
777            # This is where Inflow objects are evaluated and update the domain
778                self.outlet.momentum[0](domain)
779                self.outlet.momentum[1](domain)       
780           
781
782           
783        # Log timeseries to file
784        try:
785            fid = open(self.timeseries_filename, 'a')       
786        except:
787            pass
788        else:   
789            fid.write('%.2f, %.2f\n' %(time, Q))
790            fid.close()
791
792        # Store value of time
793        self.last_time = time
794
795
796           
797
798Culvert_flow = Culvert_flow_general       
Note: See TracBrowser for help on using the repository browser.