# source:branches/numpy/anuga/culvert_flows/culvert_class.py

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

Back-merge from Numeric trunk.

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