Context Navigation

← Previous Change
Next Change →

culvert_operator.py

Timestamp:

Aug 19, 2010, 4:41:10 PM (14 years ago)

Author:

steve

Message:

Changing culvert forcing term class to a culvert fractional step operator

File:

: 1 moved

trunk/anuga_core/source/anuga/structures/culvert_operator.py (moved) (moved from trunk/anuga_core/source/anuga/structures/culvert_class.py) (7 diffs)

Legend:

: Unmodified
: Added
: Removed

trunk/anuga_core/source/anuga/structures/culvert_operator.py

-                      r7952
+                      r7955
 from anuga.structures.culvert_polygons import create_culvert_polygons
 from anuga.utilities.system_tools import log_to_file
+from anuga.geometry.polygon import inside_polygon, is_inside_polygon, plot_polygons
+from anuga.geometry.polygon import inside_polygon, is_inside_polygon
+from anuga.geometry.polygon import plot_polygons, polygon_area
 from anuga.utilities.numerical_tools import mean
 …
+class Culvert_flow_general:
+    """Culvert flow - transfer water from one hole to another
+    This version will work with either rating curve file or with culvert
+    routine.
+class Generic_box_culvert:
+    """Culvert flow - transfer water from one rectngular box to another.
+    Sets up the geometry of problem
+    This is the base class for culverts. Inherit from this class (and overwrite
+    compute_discharge method for specific subclasses)
     Input: Two points, pipe_size (either diameter or width, height),
 …
     def __init__(self,
                  domain,
-                 culvert_description_filename=None,
-                 culvert_routine=None,
                  end_point0=None,
                  end_point1=None,
+                 enquiry_point0=None,
+                 enquiry_point1=None,
+                 type='box',
+                 enquiry_gap_factor=0.2,
                  width=None,
                  height=None,
-                 length=None,
-                 number_of_barrels=1,
-                 number_of_smoothing_steps=2000,
-                 trigger_depth=0.01, # Depth below which no flow happens
-                 manning=None,          # Mannings Roughness for Culvert
-                 sum_loss=None,
-                 use_velocity_head=False, # FIXME(Ole): Get rid of - always True
-                 use_momentum_jet=False, # FIXME(Ole): Not yet implemented
-                 label=None,
-                 description=None,
-                 update_interval=None,
-                 log_file=False,
-                 discharge_hydrograph=False,
                  verbose=False):
         # Input check
+        if height is None: height = width
+        assert number_of_barrels >= 1
+        assert use_velocity_head is True or use_velocity_head is False
+        #msg = 'Momentum jet not yet moved to general culvert'
+        #assert use_momentum_jet is False, msg
+        self.use_momentum_jet = use_momentum_jet
+        self.culvert_routine = culvert_routine
+        self.culvert_description_filename = culvert_description_filename
+        if culvert_description_filename is not None:
+            label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
+            self.rating_curve = ensure_numeric(rating_curve)
+        if height is None:
+            height = width
         self.height = height
         self.width = width
+        self.domain = domain
+        self.trigger_depth = trigger_depth
+        if manning is None:
+            self.manning = 0.012   # Default roughness for pipe
+        if sum_loss is None:
+            self.sum_loss = 0.0
+        # Store culvert information
+        self.label = label
+        self.description = description
+        self.culvert_type = type
+        self.number_of_barrels = number_of_barrels
+        # Store options
+        self.use_velocity_head = use_velocity_head
+        if label is None: label = 'culvert_flow'
+        label += '_' + str(id(self))
+        self.label = label
+        # File for storing discharge_hydrograph
+        if discharge_hydrograph is True:
+            self.timeseries_filename = label + '_timeseries.csv'
+            fid = open(self.timeseries_filename, 'w')
+            fid.write('time, discharge\n')
+            fid.close()
+        # Log file for storing general textual output
+        if log_file is True:
+            self.log_filename = label + '.log'
+            log_to_file(self.log_filename, self.label)
+            log_to_file(self.log_filename, description)
+            log_to_file(self.log_filename, self.culvert_type)
+        else:
+            self.log_filename = None
+        self.domain = domain
+        self.end_points= [end_point0, end_point1]
+        self.enquiry_gap_factor=enquiry_gap_factor
+        self.verbose=verbose
+        self.filename = None
         # Create the fundamental culvert polygons from polygon
+        P = create_culvert_polygons(end_point0,
+                                    end_point1,
+                                    width=width,
+                                    height=height,
+                                    number_of_barrels=number_of_barrels)
+        self.culvert_polygons = P
+        # Select enquiry points
+        if enquiry_point0 is None:
+            enquiry_point0 = P['enquiry_point0']
+        if enquiry_point1 is None:
+            enquiry_point1 = P['enquiry_point1']
+        if verbose is True:
+            pass
+            #plot_polygons([[end_point0, end_point1],
+            #               P['exchange_polygon0'],
+            #               P['exchange_polygon1'],
+            #               [enquiry_point0, 1.005*enquiry_point0],
+            #               [enquiry_point1, 1.005*enquiry_point1]],
+            #              figname='culvert_polygon_output')
+        self.enquiry_points = [enquiry_point0, enquiry_point1]
+        self.enquiry_indices = self.get_enquiry_indices()
+        self.create_culvert_polygons()
+        self.compute_enquiry_indices()
         self.check_culvert_inside_domain()
+        # Create inflow object at each end of the culvert.
+        self.openings = []
+        self.openings.append(Inflow(domain,
+                                    polygon=P['exchange_polygon0']))
+        self.openings.append(Inflow(domain,
+                                    polygon=P['exchange_polygon1']))
+        # Assume two openings for now: Referred to as 0 and 1
+        assert len(self.openings) == 2
         # Establish initial values at each enquiry point
 …
             s = 'Culvert Direction is %s\n' %str(self.vector)
             log_to_file(self.log_filename, s)
+    def set_store_hydrograph_discharge(self,filename=None):
+        if filename is None:
+            self.filename = 'culvert_discharge_hydrograph'
+        else:
+            self.filename = filename
+        self.discharge_hydrograph = True
+        self.timeseries_filename = self.filename + '_timeseries.csv'
+        fid = open(self.timeseries_filename, 'w')
+        fid.write('time, discharge\n')
+        fid.close()
+    def create_culvert_polygons(self):
+        """Create polygons at the end of a culvert inlet and outlet.
+        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
+        for enquiring the total energy at both ends of the culvert and transferring flow.
+        """
+        # Calculate geometry
+        x0, y0 = self.end_point0
+        x1, y1 = self.end_point1
+        dx = x1-x0
+        dy = y1-y0
+        self.culvert_vector = num.array([dx, dy])
+        self.culvert_length = sqrt(num.sum(self.culvert_vector**2))
+        # Unit direction vector and normal
+        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
+        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
+        # Short hands
+        w = 0.5*width*normal # Perpendicular vector of 1/2 width
+        h = height*vector    # Vector of length=height in the
+                             # direction of the culvert
+        gap = (1 + enquiry_gap_factor)*h
+        self.exchange_polygons = []
+        # Build exchange polygon and enquiry point for opening 0
+        p0 = self.end_point0 + w
+        p1 = self.end_point0 - w
+        p2 = p1 - h
+        p3 = p0 - h
+        self.exchange_polygons.append(num.array([p0,p1,p2,p3]))
+        self.enquiry_points= end_point0 - gap
+        # Build exchange polygon and enquiry point for opening 1
+        p0 = end_point1 + w
+        p1 = end_point1 - w
+        p2 = p1 + h
+        p3 = p0 + h
+        self.exchange_polygon1 = num.array([p0,p1,p2,p3])
+        self.enquiry_point1 = end_point1 + gap
+        # Check that enquiry point0 is outside exchange polygon0
+        polygon = self.exchange_polygon0
+        area = polygon_area(polygon)
+            msg = 'Polygon %s ' %(polygon)
+            msg += ' has area = %f' % area
+            assert area > 0.0, msg
+            for key2 in ['enquiry_point0', 'enquiry_point1']:
+                point = culvert_polygons[key2]
+                msg = 'Enquiry point falls inside an enquiry point.'
+                assert not inside_polygon(point, polygon), msg
+        for key1 in ['exchange_polygon0',
+                     'exchange_polygon1']:
+        # Return results
+        self.culvert_polygons = culvert_polygons
     def __call__(self, domain):
 …
     def get_enquiry_indices(self):
+    def compute_enquiry_indices(self):
         """Get indices for nearest centroids to self.enquiry_points
         """
 …
                 raise Exception, msg
         return enquiry_indices
+        self.enquiry_indices = enquiry_indices
 …
-# OBSOLETE (Except for momentum jet in Culvert_flow_energy)
-class Culvert_flow_rating:
-    """Culvert flow - transfer water from one hole to another
-    Input: Two points, pipe_size (either diameter or width, height),
-    mannings_rougness,
-    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
-    top-down_blockage_factor and bottom_up_blockage_factor
-    """
-    def __init__(self,
-                 domain,
-                 culvert_description_filename=None,
-                 end_point0=None,
-                 end_point1=None,
-                 enquiry_point0=None,
-                 enquiry_point1=None,
-                 update_interval=None,
-                 log_file=False,
-                 discharge_hydrograph=False,
-                 verbose=False):
-        label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
-        # Store culvert information
-        self.label = label
-        self.description = description
-        self.culvert_type = type
-        self.rating_curve = ensure_numeric(rating_curve)
-        self.number_of_barrels = number_of_barrels
-        if label is None: label = 'culvert_flow'
-        label += '_' + str(id(self))
-        self.label = label
-        # File for storing discharge_hydrograph
-        if discharge_hydrograph is True:
-            self.timeseries_filename = label + '_timeseries.csv'
-            fid = open(self.timeseries_filename, 'w')
-            fid.write('time, discharge\n')
-            fid.close()
-        # Log file for storing general textual output
-        if log_file is True:
-            self.log_filename = label + '.log'
-            log_to_file(self.log_filename, self.label)
-            log_to_file(self.log_filename, description)
-            log_to_file(self.log_filename, self.culvert_type)
-        # Create the fundamental culvert polygons from POLYGON
-        #if self.culvert_type == 'circle':
-        #    # Redefine width and height for use with create_culvert_polygons
-        #    width = height = diameter
-        P = create_culvert_polygons(end_point0,
-                                    end_point1,
-                                    width=width,
-                                    height=height,
-                                    number_of_barrels=number_of_barrels)
-        # Select enquiry points
-        if enquiry_point0 is None:
-            enquiry_point0 = P['enquiry_point0']
-        if enquiry_point1 is None:
-            enquiry_point1 = P['enquiry_point1']
-        if verbose is True:
-            pass
-            #plot_polygons([[end_point0, end_point1],
-            #               P['exchange_polygon0'],
-            #               P['exchange_polygon1'],
-            #               [enquiry_point0, 1.005*enquiry_point0],
-            #               [enquiry_point1, 1.005*enquiry_point1]],
-            #              figname='culvert_polygon_output')
-        self.enquiry_points = [enquiry_point0, enquiry_point1]
-        self.enquiry_indices = []
-        for point in self.enquiry_points:
-            # Find nearest centroid
-            N = len(domain)
-            points = domain.get_centroid_coordinates(absolute=True)
-            # Calculate indices in exchange area for this forcing term
-            triangle_id = min_dist = sys.maxint
-            for k in range(N):
-                x, y = points[k,:] # Centroid
-                c = point
-                distance = (x-c[0])**2+(y-c[1])**2
-                if distance < min_dist:
-                    min_dist = distance
-                    triangle_id = k
-            if triangle_id < sys.maxint:
-                msg = 'found triangle with centroid (%f, %f)'\
-                    %tuple(points[triangle_id, :])
-                msg += ' for point (%f, %f)' %tuple(point)
-                self.enquiry_indices.append(triangle_id)
-            else:
-                msg = 'Triangle not found for point (%f, %f)' %point
-                raise Exception, msg
-        # Check that all polygons lie within the mesh.
-        bounding_polygon = domain.get_boundary_polygon()
-        for key in P.keys():
-            if key in ['exchange_polygon0',
-                       'exchange_polygon1']:
-                for point in list(P[key]) + self.enquiry_points:
-                    msg = 'Point %s in polygon %s for culvert %s did not'\
-                        %(str(point), key, self.label)
-                    msg += 'fall within the domain boundary.'
-                    assert is_inside_polygon(point, bounding_polygon), msg
-        # Create inflow object at each end of the culvert.
-        self.openings = []
-        self.openings.append(Inflow(domain,
-                                    polygon=P['exchange_polygon0']))
-        self.openings.append(Inflow(domain,
-                                    polygon=P['exchange_polygon1']))
-        dq = domain.quantities
-        for i, opening in enumerate(self.openings):
-            elevation = dq['elevation'].get_values(location='centroids',
-                                                   indices=[self.enquiry_indices[i]])
-            opening.elevation = elevation
-            opening.stage = elevation # Simple assumption that culvert is dry initially
-        # Assume two openings for now: Referred to as 0 and 1
-        assert len(self.openings) == 2
-        # Determine pipe direction
-        self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation
-        if delta_z > 0.0:
-            self.inlet = self.openings[0]
-            self.outlet = self.openings[1]
-        else:
-            self.outlet = self.openings[0]
-            self.inlet = self.openings[1]
-        # Store basic geometry
-        self.end_points = [end_point0, end_point1]
-        self.vector = P['vector']
-        self.length = P['length']; assert self.length > 0.0
-        if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
-            msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length)
-            msg += ' does not match distance between specified'
-            msg += ' end points (%.2f m)' %self.length
-            log.critical(msg)
-        self.verbose = verbose
-        self.last_update = 0.0 # For use with update_interval
-        self.last_time = 0.0
-        self.update_interval = update_interval
-        # Print something
-        if hasattr(self, 'log_filename'):
-            s = 'Culvert Effective Length = %.2f m' %(self.length)
-            log_to_file(self.log_filename, s)
-            s = 'Culvert Direction is %s\n' %str(self.vector)
-            log_to_file(self.log_filename, s)
-    def __call__(self, domain):
-        # Time stuff
-        time = domain.get_time()
-        update = False
-        if self.update_interval is None:
-            update = True
-            delta_t = domain.timestep # Next timestep has been computed in domain.py
-        else:
-            if time - self.last_update > self.update_interval or time == 0.0:
-                update = True
-            delta_t = self.update_interval
-        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
-        if hasattr(self, 'log_filename'):
-            log_to_file(self.log_filename, s)
-        if update is True:
-            self.last_update = time
-            dq = domain.quantities
-            # Get average water depths at each opening
-            openings = self.openings   # There are two Opening [0] and [1]
-            for i, opening in enumerate(openings):
-                # Compute mean values of selected quantitites in the
-                # enquiry area in front of the culvert
-                stage = dq['stage'].get_values(location='centroids',
-                                               indices=[self.enquiry_indices[i]])
-                # Store current average stage and depth with each opening object
-                opening.depth = stage - opening.elevation
-                opening.stage = stage
-            #################  End of the FOR loop ################################################
-            # We now need to deal with each opening individually
-            # Determine flow direction based on total energy difference
-            delta_w = self.inlet.stage - self.outlet.stage
-            if hasattr(self, 'log_filename'):
-                s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time,
-                                                                           self.inlet.stage,
-                                                                           self.outlet.stage)
-                log_to_file(self.log_filename, s)
-            if self.inlet.depth <= 0.01:
-                Q = 0.0
-            else:
-                # Calculate discharge for one barrel and set inlet.rate and outlet.rate
-                try:
-                    Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1])
-                except Below_interval, e:
-                    Q = self.rating_curve[0,1]
-                    msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
-                    msg += str(e)
-                    msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
-                        %(Q, self.label)
-                    if hasattr(self, 'log_filename'):
-                        log_to_file(self.log_filename, msg)
-                except Above_interval, e:
-                    Q = self.rating_curve[-1,1]
-                    msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
-                    msg += str(e)
-                    msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
-                        %(Q, self.label)
-                    if hasattr(self, 'log_filename'):
-                        log_to_file(self.log_filename, msg)
-            # Adjust discharge for multiple barrels
-            Q *= self.number_of_barrels
-            # Adjust Q downwards depending on available water at inlet
-            stage = self.inlet.get_quantity_values(quantity_name='stage')
-            elevation = self.inlet.get_quantity_values(quantity_name='elevation')
-            depth = stage-elevation
-            V = 0
-            for i, d in enumerate(depth):
-                V += d * domain.areas[i]
-            dt = delta_t
-            if Q*dt > V:
-                Q_reduced = 0.9*V/dt # Reduce with safety factor
-                msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
-                msg += 'greater than current volume (V) at inlet.\n'
-                msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
-                if self.verbose is True:
-                    log.critical(msg)
-                if hasattr(self, 'log_filename'):
-                    log_to_file(self.log_filename, msg)
-                Q = Q_reduced
-            self.inlet.rate = -Q
-            self.outlet.rate = Q
-            # Log timeseries to file
-            try:
-                fid = open(self.timeseries_filename, 'a')
-            except:
-                pass
-            else:
-                fid.write('%.2f, %.2f\n' %(time, Q))
-                fid.close()
-            # Store value of time
-            self.last_time = time
-        # Execute flow term for each opening
-        # This is where Inflow objects are evaluated using the last rate that has been calculated
+        #
-        # This will take place at every internal timestep and update the domain
-        for opening in self.openings:
-            opening(domain)
-class Culvert_flow_energy:
-    """Culvert flow - transfer water from one hole to another
-    Using Momentum as Calculated by Culvert Flow !!
-    Could be Several Methods Investigated to do This !!!
-_May_08
-    To Ole:
-    OK so here we need to get the Polygon Creating code to create
-    polygons for the culvert Based on
-    the two input Points (X0,Y0) and (X1,Y1) - need to be passed
-    to create polygon
-    The two centers are now passed on to create_polygon.
-    Input: Two points, pipe_size (either diameter or width, height),
-    mannings_rougness,
-    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
-    top-down_blockage_factor and bottom_up_blockage_factor
-    And the Delta H enquiry should be change from Openings in line 412
-    to the enquiry Polygons infront of the culvert
-    At the moment this script uses only Depth, later we can change it to
-    Energy...
-    Once we have Delta H can calculate a Flow Rate and from Flow Rate
-    an Outlet Velocity
-    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
-    Invert levels are optional. If left out they will default to the
-    elevation at the opening.
-    """
-    def __init__(self,
-                 domain,
-                 label=None,
-                 description=None,
-                 end_point0=None,
-                 end_point1=None,
-                 width=None,
-                 height=None,
-                 diameter=None,
-                 manning=None,          # Mannings Roughness for Culvert
-                 invert_level0=None,    # Invert level at opening 0
-                 invert_level1=None,    # Invert level at opening 1
-                 loss_exit=None,
-                 loss_entry=None,
-                 loss_bend=None,
-                 loss_special=None,
-                 blockage_topdwn=None,
-                 blockage_bottup=None,
-                 culvert_routine=None,
-                 number_of_barrels=1,
-                 enquiry_point0=None,
-                 enquiry_point1=None,
-                 update_interval=None,
-                 verbose=False):
-        # Input check
-        if diameter is not None:
-            self.culvert_type = 'circle'
-            self.diameter = diameter
-            if height is not None or width is not None:
-                msg = 'Either diameter or width&height must be specified, '
-                msg += 'but not both.'
-                raise Exception, msg
-        else:
-            if height is not None:
-                if width is None:
-                    self.culvert_type = 'square'
-                    width = height
-                else:
-                    self.culvert_type = 'rectangle'
-            elif width is not None:
-                if height is None:
-                    self.culvert_type = 'square'
-                    height = width
-            else:
-                msg = 'Either diameter or width&height must be specified.'
-                raise Exception, msg
-            if height == width:
-                self.culvert_type = 'square'
-            self.height = height
-            self.width = width
-        assert self.culvert_type in ['circle', 'square', 'rectangle']
-        assert number_of_barrels >= 1
-        self.number_of_barrels = number_of_barrels
-        # Set defaults
-        if manning is None: manning = 0.012   # Default roughness for pipe
-        if loss_exit is None: loss_exit = 1.00
-        if loss_entry is None: loss_entry = 0.50
-        if loss_bend is None: loss_bend=0.00
-        if loss_special is None: loss_special=0.00
-        if blockage_topdwn is None: blockage_topdwn=0.00
-        if blockage_bottup is None: blockage_bottup=0.00
-        if culvert_routine is None:
-            culvert_routine=boyd_generalised_culvert_model
-        if label is None: label = 'culvert_flow'
-        label += '_' + str(id(self))
-        self.label = label
-        # File for storing culvert quantities
-        self.timeseries_filename = label + '_timeseries.csv'
-        fid = open(self.timeseries_filename, 'w')
-        fid.write('time, E0, E1, Velocity, Discharge\n')
-        fid.close()
-        # Log file for storing general textual output
-        self.log_filename = label + '.log'
-        log_to_file(self.log_filename, self.label)
-        log_to_file(self.log_filename, description)
-        log_to_file(self.log_filename, self.culvert_type)
-        # Create the fundamental culvert polygons from POLYGON
-        if self.culvert_type == 'circle':
-            # Redefine width and height for use with create_culvert_polygons
-            width = height = diameter
-        P = create_culvert_polygons(end_point0,
-                                    end_point1,
-                                    width=width,
-                                    height=height,
-                                    number_of_barrels=number_of_barrels)
-        # Select enquiry points
-        if enquiry_point0 is None:
-            enquiry_point0 = P['enquiry_point0']
-        if enquiry_point1 is None:
-            enquiry_point1 = P['enquiry_point1']
-        if verbose is True:
-            pass
-            #plot_polygons([[end_point0, end_point1],
-            #               P['exchange_polygon0'],
-            #               P['exchange_polygon1'],
-            #               [enquiry_point0, 1.005*enquiry_point0],
-            #               [enquiry_point1, 1.005*enquiry_point1]],
-            #              figname='culvert_polygon_output')
-        self.enquiry_points = [enquiry_point0, enquiry_point1]
-        self.enquiry_indices = []
-        for point in self.enquiry_points:
-            # Find nearest centroid
-            N = len(domain)
-            points = domain.get_centroid_coordinates(absolute=True)
-            # Calculate indices in exchange area for this forcing term
-            triangle_id = min_dist = sys.maxint
-            for k in range(N):
-                x, y = points[k,:] # Centroid
-                c = point
-                distance = (x-c[0])**2+(y-c[1])**2
-                if distance < min_dist:
-                    min_dist = distance
-                    triangle_id = k
-            if triangle_id < sys.maxint:
-                msg = 'found triangle with centroid (%f, %f)'\
-                    %tuple(points[triangle_id, :])
-                msg += ' for point (%f, %f)' %tuple(point)
-                self.enquiry_indices.append(triangle_id)
-            else:
-                msg = 'Triangle not found for point (%f, %f)' %point
-                raise Exception, msg
-        # Check that all polygons lie within the mesh.
-        bounding_polygon = domain.get_boundary_polygon()
-        for key in P.keys():
-            if key in ['exchange_polygon0',
-                       'exchange_polygon1']:
-                for point in P[key]:
-                    msg = 'Point %s in polygon %s for culvert %s did not'\
-                        %(str(point), key, self.label)
-                    msg += 'fall within the domain boundary.'
-                    assert is_inside_polygon(point, bounding_polygon), msg
-        # Create inflow object at each end of the culvert.
-        self.openings = []
-        self.openings.append(Inflow(domain,
-                                    polygon=P['exchange_polygon0']))
-        self.openings.append(Inflow(domain,
-                                    polygon=P['exchange_polygon1']))
-        # Assume two openings for now: Referred to as 0 and 1
-        assert len(self.openings) == 2
-        # Store basic geometry
-        self.end_points = [end_point0, end_point1]
-        self.invert_levels = [invert_level0, invert_level1]
-        #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
-        #self.enquiry_polylines = [P['enquiry_polygon0'][:2],
-        #                          P['enquiry_polygon1'][:2]]
-        self.vector = P['vector']
-        self.length = P['length']; assert self.length > 0.0
-        self.verbose = verbose
-        self.last_time = 0.0
-        self.last_update = 0.0 # For use with update_interval
-        self.update_interval = update_interval
-        # Store hydraulic parameters
-        self.manning = manning
-        self.loss_exit = loss_exit
-        self.loss_entry = loss_entry
-        self.loss_bend = loss_bend
-        self.loss_special = loss_special
-        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
-        self.blockage_topdwn = blockage_topdwn
-        self.blockage_bottup = blockage_bottup
-        # Store culvert routine
-        self.culvert_routine = culvert_routine
-        # Create objects to update momentum (a bit crude at this stage)
-        xmom0 = General_forcing(domain, 'xmomentum',
-                                polygon=P['exchange_polygon0'])
-        xmom1 = General_forcing(domain, 'xmomentum',
-                                polygon=P['exchange_polygon1'])
-        ymom0 = General_forcing(domain, 'ymomentum',
-                                polygon=P['exchange_polygon0'])
-        ymom1 = General_forcing(domain, 'ymomentum',
-                                polygon=P['exchange_polygon1'])
-        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
-        # Print something
-        s = 'Culvert Effective Length = %.2f m' %(self.length)
-        log_to_file(self.log_filename, s)
-        s = 'Culvert Direction is %s\n' %str(self.vector)
-        log_to_file(self.log_filename, s)
-    def __call__(self, domain):
-        log_filename = self.log_filename
-        # Time stuff
-        time = domain.get_time()
-        # Short hand
-        dq = domain.quantities
-        update = False
-        if self.update_interval is None:
-            update = True
-            delta_t = domain.timestep # Next timestep has been computed in domain.py
-        else:
-            if time - self.last_update > self.update_interval or time == 0.0:
-                update = True
-            delta_t = self.update_interval
-        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
-        if hasattr(self, 'log_filename'):
-            log_to_file(log_filename, s)
-        if update is True:
-            self.last_update = time
-            msg = 'Time did not advance'
-            if time > 0.0: assert delta_t > 0.0, msg
-            # Get average water depths at each opening
-            openings = self.openings   # There are two Opening [0] and [1]
-            for i, opening in enumerate(openings):
-                # Compute mean values of selected quantitites in the
-                # exchange area in front of the culvert
-                stage = opening.get_quantity_values(quantity_name='stage')
-                w = mean(stage) # Average stage
-                # Use invert level instead of elevation if specified
-                invert_level = self.invert_levels[i]
-                if invert_level is not None:
-                    z = invert_level
-                else:
-                    elevation = opening.get_quantity_values(quantity_name='elevation')
-                    z = mean(elevation) # Average elevation
-                # Estimated depth above the culvert inlet
-                d = w - z  # Used for calculations involving depth
-                if d < 0.0:
-                    # This is possible since w and z are taken at different locations
-                    #msg = 'D < 0.0: %f' %d
-                    #raise msg
-                    d = 0.0
-                # Ratio of depth to culvert height.
-                # If ratio > 1 then culvert is running full
-                if self.culvert_type == 'circle':
-                    ratio = d/self.diameter
-                else:
-                    ratio = d/self.height
-                opening.ratio = ratio
-                # Average measures of energy in front of this opening
-                id = [self.enquiry_indices[i]]
-                stage = dq['stage'].get_values(location='centroids',
-                                               indices=id)
-                elevation = dq['elevation'].get_values(location='centroids',
-                                                       indices=id)
-                xmomentum = dq['xmomentum'].get_values(location='centroids',
-                                                       indices=id)
-                ymomentum = dq['xmomentum'].get_values(location='centroids',
-                                                       indices=id)
-                depth = stage-elevation
-                if depth > 0.0:
-                    u = xmomentum/(depth + velocity_protection/depth)
-                    v = ymomentum/(depth + velocity_protection/depth)
-                else:
-                    u = v = 0.0
-                opening.total_energy = 0.5*(u*u + v*v)/g + stage
-                # Store current average stage and depth with each opening object
-                opening.depth = d
-                opening.depth_trigger = d # Use this for now
-                opening.stage = w
-                opening.elevation = z
-            #################  End of the FOR loop ################################################
-            # We now need to deal with each opening individually
-            # Determine flow direction based on total energy difference
-            delta_Et = openings[0].total_energy - openings[1].total_energy
-            if delta_Et > 0:
-                inlet = openings[0]
-                outlet = openings[1]
-                inlet.momentum = self.opening_momentum[0]
-                outlet.momentum = self.opening_momentum[1]
-            else:
-                inlet = openings[1]
-                outlet = openings[0]
-                inlet.momentum = self.opening_momentum[1]
-                outlet.momentum = self.opening_momentum[0]
-                delta_Et = -delta_Et
-            self.inlet = inlet
-            self.outlet = outlet
-            msg = 'Total energy difference is negative'
-            assert delta_Et > 0.0, msg
-            delta_h = inlet.stage - outlet.stage
-            delta_z = inlet.elevation - outlet.elevation
-            culvert_slope = (delta_z/self.length)
-            if culvert_slope < 0.0:
-                # Adverse gradient - flow is running uphill
-                # Flow will be purely controlled by uphill outlet face
-                if self.verbose is True:
-                    log.critical('WARNING: Flow is running uphill. Watch Out! '
-                                 'inlet.elevation=%s, outlet.elevation%s'
-                                 % (str(inlet.elevation), str(outlet.elevation)))
-            s = 'Delta total energy = %.3f' %(delta_Et)
-            log_to_file(log_filename, s)
-            # Calculate discharge for one barrel and set inlet.rate and outlet.rate
-            Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
-            # Adjust discharge for multiple barrels
-            Q *= self.number_of_barrels
-            # Compute barrel momentum
-            barrel_momentum = barrel_velocity*culvert_outlet_depth
-            s = 'Barrel velocity = %f' %barrel_velocity
-            log_to_file(log_filename, s)
-            # Compute momentum vector at outlet
-            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
-            s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
-            log_to_file(log_filename, s)
-            # Log timeseries to file
-            fid = open(self.timeseries_filename, 'a')
-            fid.write('%f, %f, %f, %f, %f\n'\
-                          %(time,
-                            openings[0].total_energy,
-                            openings[1].total_energy,
-                            barrel_velocity,
-                            Q))
-            fid.close()
-            # Update momentum
-            if delta_t > 0.0:
-                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
-                xmomentum_rate /= delta_t
-                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
-                ymomentum_rate /= delta_t
-                s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
-                log_to_file(log_filename, s)
-            else:
-                xmomentum_rate = ymomentum_rate = 0.0
-            # Set momentum rates for outlet jet
-            outlet.momentum[0].rate = xmomentum_rate
-            outlet.momentum[1].rate = ymomentum_rate
-            # Remember this value for next step (IMPORTANT)
-            outlet.momentum[0].value = outlet_mom_x
-            outlet.momentum[1].value = outlet_mom_y
-            if int(domain.time*100) % 100 == 0:
-                s = 'T=%.5f, Culvert Discharge = %.3f f'\
-                    %(time, Q)
-                s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
-                     %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
-                s += ' Momentum rate: (%.4f, %.4f)'\
-                     %(xmomentum_rate, ymomentum_rate)
-                s+='Outlet Vel= %.3f'\
-                    %(barrel_velocity)
-                log_to_file(log_filename, s)
-            # Store value of time
-            self.last_time = time
-        # Execute flow term for each opening
-        # This is where Inflow objects are evaluated and update the domain
-        for opening in self.openings:
-            opening(domain)
-        # Execute momentum terms
-        # This is where Inflow objects are evaluated and update the domain
-        self.outlet.momentum[0](domain)
-        self.outlet.momentum[1](domain)
 Culvert_flow = Culvert_flow_general

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 7955 for trunk/anuga_core/source/anuga/structures/culvert_operator.py

Legend:

trunk/anuga_core/source/anuga/structures/culvert_operator.py

Download in other formats: