Changeset 6635


Ignore:
Timestamp:
Mar 26, 2009, 4:40:32 PM (16 years ago)
Author:
ted
Message:

Fixed use of wrong areas in General_forcing and use actual sum of trianle areas instead. Added unit test for sub to super critical regimes.

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

Legend:

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

    r6503 r6635  
    17511751            assert polygon is None, msg
    17521752
    1753             self.exchange_area = radius**2*pi
    17541753
    17551754            # Check that circle center lies within the mesh.
     
    17851784                assert is_inside_polygon(point, bounding_polygon), msg
    17861785               
    1787             # Compute area and check that it is greater than 0   
    1788             self.exchange_area = polygon_area(self.polygon)
    1789            
    1790             msg = 'Polygon %s in forcing term' %(self.polygon)
    1791             msg += ' has area = %f' %self.exchange_area
    1792             assert self.exchange_area > 0.0           
    1793            
    1794                
    1795                
     1786                                 
    17961787
    17971788        # Pointer to update vector
     
    18201811            # Inlet is polygon
    18211812           
    1822             inlet_region = 'polygon=%s, area=%f m^2' %(self.polygon,
    1823                                                        self.exchange_area)
    1824                        
     1813            inlet_region = 'polygon=%s' % (self.polygon)
    18251814            self.exchange_indices = inside_polygon(points, self.polygon)
    18261815           
     
    18341823                msg += 'specified region: %s' %inlet_region
    18351824                raise Exception, msg
     1825
     1826
     1827
     1828        # Compute exchange area as the sum of areas of triangles identified
     1829        # by circle or polygon
     1830        self.exchange_area = 0.0
     1831        for i in self.exchange_indices:
     1832            self.exchange_area += domain.areas[i]
     1833           
     1834
     1835        msg = 'Exchange area in forcing term'
     1836        msg += ' has area = %f' %self.exchange_area
     1837        assert self.exchange_area > 0.0           
     1838           
     1839               
     1840
    18361841           
    18371842        # Check and store default_rate
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r6541 r6635  
    67296729                   
    67306730       
    6731 
     6731    def test_inflow_using_flowline(self):
     6732        """test_inflow_using_flowline
     6733        Test the ability of a flowline to match inflow above the flowline by
     6734        creating constant inflow onto a circle at the head of a 20m
     6735        wide by 300m long plane dipping at 1:300 with a perpendicular flowline and gauge
     6736        downstream of the inflow and a 45 degree flowlines at 60m downstream
     6737        """
     6738
     6739        verbose = False
     6740       
     6741
     6742        #------------------------------------------------------------------------------
     6743        # Import necessary modules
     6744        #------------------------------------------------------------------------------
     6745        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     6746        from anuga.shallow_water import Domain
     6747        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     6748        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     6749        from anuga.shallow_water.shallow_water_domain import Inflow
     6750        from anuga.shallow_water.data_manager import get_flow_through_cross_section
     6751        from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs
     6752
     6753
     6754        #------------------------------------------------------------------------------
     6755        # Setup computational domain
     6756        #------------------------------------------------------------------------------
     6757        number_of_inflows = 2 # Number of inflows on top of each other
     6758
     6759        length = 300.
     6760        width  = 20.
     6761        dx = dy = 2           # Resolution: of grid on both axes
     6762
     6763        points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
     6764                                                       len1=length, len2=width)
     6765
     6766
     6767        for slope in [1.0/300, 1.0/200, 1.0/100, 1.0/50]:
     6768            # Loop over a range of bedslopes representing sub to super critical flows
     6769
     6770            #print slope
     6771            domain = Domain(points, vertices, boundary)   
     6772            domain.set_name('Inflow_flowline_test')              # Output name
     6773           
     6774
     6775            #------------------------------------------------------------------------------
     6776            # Setup initial conditions
     6777            #------------------------------------------------------------------------------
     6778
     6779            def topography(x, y):
     6780                z=-x * slope
     6781                return z
     6782
     6783
     6784            domain.set_quantity('elevation', topography)  # Use function for elevation
     6785            domain.set_quantity('friction', 0.012)        # Constant friction of conc surface   
     6786            domain.set_quantity('stage',
     6787                                expression='elevation')   # Dry initial condition
     6788
     6789
     6790            #------------------------------------------------------------------------------
     6791            # Setup boundary conditions
     6792            #------------------------------------------------------------------------------
     6793
     6794            Br = Reflective_boundary(domain)              # Solid reflective wall
     6795            Bo = Dirichlet_boundary([-100, 0, 0])           # Outflow stsge at -0.9m d=0.1m
     6796
     6797            domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
     6798
     6799            #------------------------------------------------------------------------------
     6800            # Seup Inflow
     6801            #------------------------------------------------------------------------------
     6802
     6803            # Fixed Flowrate onto Area
     6804            fixed_inflow = Inflow(domain,
     6805                                  center=(10.0, 10.0),
     6806                                  radius=5.00,
     6807                                  rate=1.00)   
     6808
     6809            # Stack this flow
     6810            for i in range(number_of_inflows):
     6811                domain.forcing_terms.append(fixed_inflow)
     6812           
     6813            ref_flow = fixed_inflow.rate*number_of_inflows
     6814
     6815            #------------------------------------------------------------------------------
     6816            # Evolve system through time
     6817            #------------------------------------------------------------------------------
     6818
     6819            for t in domain.evolve(yieldstep = 10.0, finaltime = 300):
     6820                pass
     6821                #print t
     6822
     6823            #------------------------------------------------------------------------------
     6824            # Compute flow thru flowlines ds of inflow
     6825            #------------------------------------------------------------------------------
     6826               
     6827            # Square on flowline at 30m
     6828            q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]])
     6829            msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
     6830            #print q, ref_flow
     6831            assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     6832
     6833                       
     6834            # 45 degree flowline at 60m
     6835            q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]])
     6836            msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
     6837            #print q, ref_flow
     6838            assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     6839           
     6840           
     6841       
     6842       
     6843
     6844   
    67326845       
    67336846if __name__ == "__main__":
    6734     suite = unittest.makeSuite(Test_Shallow_Water, 'test')
     6847    suite = unittest.makeSuite(Test_Shallow_Water, 'test_inflow_using_flowline')
    67356848    runner = unittest.TextTestRunner(verbosity=1)   
    67366849    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.