Changeset 7043


Ignore:
Timestamp:
May 15, 2009, 2:33:56 PM (11 years ago)
Author:
ole
Message:

Added work on Inflow_boundary. This is only a start and the test
is disabled.

Location:
anuga_core/source/anuga/shallow_water
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r6928 r7043  
    121121from warnings import warn
    122122
    123 
     123import math
    124124
    125125#
     
    14011401   
    14021402
    1403        
     1403
     1404   
     1405class Inflow_boundary(Boundary):
     1406    """Apply given flow in m^3/s to boundary segment.
     1407    Depth and momentum is derived using Manning's formula.
     1408
     1409    Underlying domain must be specified when boundary is instantiated
     1410    """
     1411   
     1412    # FIXME (Ole): This is work in progress and definitely not finished.
     1413    # The associated test has been disabled
     1414
     1415    def __init__(self, domain=None, rate=0.0):
     1416        Boundary.__init__(self)
     1417
     1418        if domain is None:
     1419            msg = 'Domain must be specified for '
     1420            msg += 'Inflow boundary'
     1421            raise Exception, msg
     1422
     1423        self.domain = domain
     1424       
     1425        # FIXME(Ole): Allow rate to be time dependent as well
     1426        self.rate = rate
     1427        self.tag = None # Placeholder for tag associated with this object.
     1428
     1429    def __repr__(self):
     1430        return 'Inflow_boundary(%s)' %self.domain
     1431
     1432    def evaluate(self, vol_id, edge_id):
     1433        """Apply inflow rate at each edge of this boundary
     1434        """
     1435       
     1436        # First find all segments having the same tag is vol_id, edge_id
     1437        # This will be done the first time evaluate is called.
     1438        if self.tag is None:
     1439            boundary = self.domain.boundary
     1440            self.tag = boundary[(vol_id, edge_id)]       
     1441           
     1442            # Find total length of boundary with this tag
     1443            length = 0.0
     1444            for v_id, e_id in boundary:
     1445                if self.tag == boundary[(v_id, e_id)]:
     1446                    length += self.domain.mesh.get_edgelength(v_id, e_id)           
     1447
     1448            self.length = length
     1449            self.average_momentum = self.rate/length
     1450           
     1451           
     1452        # Average momentum has now been established across this boundary
     1453        # Compute momentum in the inward normal direction
     1454       
     1455        inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id)       
     1456        xmomentum, ymomentum = self.average_momentum * inward_normal
     1457           
     1458        # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S)
     1459        # Where v is velocity, n is manning's coefficient, h is depth and S is the slope into the domain.
     1460        # Let mu be the momentum (vh), then this equation becomes: mu = 1/n h^{5/3} sqrt(S)
     1461        # from which we can isolate depth to get
     1462        # h = (mu n/sqrt(S) )^{3/5}
     1463       
     1464        slope = 0 # get gradient for this triangle dot normal
     1465       
     1466        # get manning coef from this triangle
     1467        friction = self.domain.get_quantity('friction').get_values(location='edges',
     1468                                                                   indices=[vol_id])[0]
     1469        mannings_n = friction[edge_id]
     1470
     1471        if slope > epsilon and mannings_n > epsilon:
     1472            depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), 3.0/5)
     1473        else:
     1474            depth = 1.0
     1475           
     1476        # Elevation on this edge   
     1477       
     1478        z = self.domain.get_quantity('elevation').get_values(location='edges',
     1479                                                             indices=[vol_id])[0]
     1480        elevation = z[edge_id]
     1481           
     1482        # Assign conserved quantities and return
     1483        q = num.array([elevation + depth, xmomentum, ymomentum], num.Float)
     1484        return q
     1485
     1486
     1487       
     1488   
     1489           
    14041490       
    14051491class Field_boundary(Boundary):
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r6763 r7043  
    30313031       
    30323032        This test will fail as the problem was only fixed for culverts.
     3033        FIXME(Ole): It'd be nice to turn it into a volume conservation test
    30333034        """
    30343035       
     
    68856886           
    68866887           
    6887                                
     6888       
     6889       
    68886890       
    68896891    def test_inflow_using_flowline(self):
     
    69566958
    69576959                #--------------------------------------------------------------
    6958                 # Seup Inflow
     6960                # Setup Inflow
    69596961                #--------------------------------------------------------------
    69606962
     
    70357037
    70367038       
    7037        
     7039
     7040       
     7041                               
     7042       
     7043    def Xtest_inflow_boundary_using_flowline(self):
     7044        """test_inflow_boundary_using_flowline
     7045        Test the ability of a flowline to match inflow above the flowline by
     7046        creating constant inflow into the boundary at the head of a 20m
     7047        wide by 300m long plane dipping at various slopes with a
     7048        perpendicular flowline and gauge downstream of the inflow and
     7049        a 45 degree flowlines at 200m downstream
     7050       
     7051       
     7052        """
     7053
     7054        # FIXME (Ole): Work in progress
     7055       
     7056        verbose = False
     7057       
     7058
     7059        #----------------------------------------------------------------------
     7060        # Import necessary modules
     7061        #----------------------------------------------------------------------
     7062        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     7063        from anuga.shallow_water import Domain
     7064        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     7065        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     7066        from anuga.shallow_water.shallow_water_domain import Inflow_boundary
     7067        from anuga.shallow_water.data_manager import get_flow_through_cross_section
     7068        from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs
     7069
     7070
     7071        #----------------------------------------------------------------------
     7072        # Setup computational domain
     7073        #----------------------------------------------------------------------
     7074
     7075        finaltime = 500 #700.0 # If this is too short, steady state will not be achieved
     7076
     7077        length = 250.
     7078        width  = 20.
     7079        dx = dy = 5          # Resolution: of grid on both axes
     7080
     7081        points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
     7082                                                       len1=length, len2=width)
     7083
     7084        for mannings_n in [0.1, 0.01]:
     7085            # Loop over a range of roughnesses             
     7086           
     7087            for slope in [1.0/300, 1.0/100]:
     7088                # Loop over a range of bedslopes representing sub to super critical flows
     7089
     7090
     7091                domain = Domain(points, vertices, boundary)   
     7092                domain.set_name('inflow_boundary_flowline_test')
     7093               
     7094
     7095                #-------------------------------------------------------------
     7096                # Setup initial conditions
     7097                #-------------------------------------------------------------
     7098
     7099                def topography(x, y):
     7100                    z=-x * slope
     7101                    return z
     7102
     7103                domain.set_quantity('elevation', topography)
     7104                domain.set_quantity('friction', mannings_n)
     7105                domain.set_quantity('stage',
     7106                                    expression='elevation')
     7107
     7108
     7109                   
     7110                #--------------------------------------------------------------
     7111                # Setup boundary conditions
     7112                #--------------------------------------------------------------
     7113               
     7114
     7115
     7116                ref_flow = 10.00
     7117
     7118                # Compute normal depth on plane using Mannings equation
     7119                # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
     7120                normal_depth=(ref_flow*mannings_n/(slope**0.5*width))**0.6
     7121                if verbose:
     7122                    print
     7123                    print 'Slope:', slope, 'Mannings n:', mannings_n
     7124                   
     7125               
     7126               
     7127                Bi = Inflow_boundary(domain, rate=ref_flow)
     7128
     7129                Br = Reflective_boundary(domain)
     7130               
     7131                # Define downstream boundary based on predicted depth
     7132                def normal_depth_stage_downstream(t):
     7133                    return (-slope*length) + normal_depth
     7134               
     7135                Bt = Transmissive_momentum_set_stage_boundary(domain=domain,
     7136                                                              function=normal_depth_stage_downstream)
     7137               
     7138
     7139               
     7140
     7141                domain.set_boundary({'left': Bi,
     7142                                     'right': Bt,
     7143                                     'top': Br,
     7144                                     'bottom': Br})
     7145
     7146
     7147
     7148                #--------------------------------------------------------------
     7149                # Evolve system through time
     7150                #--------------------------------------------------------------
     7151
     7152
     7153                for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
     7154                    pass
     7155                    #if verbose :
     7156                    #    print domain.timestepping_statistics()
     7157                    #    print domain.volumetric_balance_statistics()                                   
     7158
     7159
     7160                #--------------------------------------------------------------
     7161                # Compute flow thru flowlines ds of inflow
     7162                #--------------------------------------------------------------
     7163                   
     7164                # Square on flowline at 200m
     7165                q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]])
     7166                msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
     7167                if verbose:
     7168                    print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
     7169                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     7170
     7171                           
     7172                # 45 degree flowline at 200m
     7173                q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]])
     7174                msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
     7175                if verbose:
     7176                    print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
     7177                   
     7178                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     7179
     7180       
     7181
     7182       
     7183               
     7184       
     7185               
    70387186    def Xtest_friction_dependent_flow_using_flowline(self):
    70397187        """test_friction_dependent_flow_using_flowline
Note: See TracChangeset for help on using the changeset viewer.