source: branches/numpy_anuga_validation/automated_validation_tests/flow_tests/test_inflow_using_flowline.py @ 6788

"""test_inflow_using_flowline
Test the ability of a flowline to match inflow above the flowline by
creating constant inflow onto a circle at the head of a 20m
wide by 400m long plane dipping at various slopes with a perpendicular flowline and gauge 
downstream of the inflow and a 45 degree flowlines at 200m downstream
"""

verbose = True
import numpy as num

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------
from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
from anuga.shallow_water import Domain
from anuga.shallow_water.shallow_water_domain import Reflective_boundary
from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
from anuga.shallow_water.shallow_water_domain import Transmissive_Momentum_Set_Stage_boundary
from anuga.shallow_water.shallow_water_domain import Inflow
from anuga.shallow_water.data_manager import get_flow_through_cross_section
from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs


#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------

number_of_inflows = 2 # Number of inflows on top of each other
finaltime = 1000.0

length = 400.
width  = 20.
dx = dy = 2          # Resolution: of grid on both axes

points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
                                               len1=length, len2=width)

for mannings_n in [0.150, 0.070, 0.035, 0.012]:
    # Loop over a range of roughnesses  
    
    for slope in [1.0/300, 1.0/150, 1.0/75]:
        # Loop over a range of bedslopes  
        

        domain = Domain(points, vertices, boundary)   
        domain.set_name('inflow_flowline_test')              # Output name


        #----------------------------------------------------------------------
        # Setup initial conditions
        #----------------------------------------------------------------------

        def topography(x, y):
            z=-x * slope
            return z

        domain.set_quantity('elevation', topography)  # Use function for elevation
        domain.set_quantity('friction', mannings_n)   # Constant friction of conc surface   
        domain.set_quantity('stage',
                            expression='elevation')   # Dry initial condition


        #----------------------------------------------------------------------
        # Seup Inflow and compute flow characteristics
        #----------------------------------------------------------------------

        # Fixed Flowrate onto Area 
        fixed_inflow = Inflow(domain,
                              center=(10.0, 10.0),
                              radius=5.00,
                              rate=10.00)   
        
        # Stack this flow
        for i in range(number_of_inflows):
            domain.forcing_terms.append(fixed_inflow)
        
        ref_flow = fixed_inflow.rate*number_of_inflows
        
        # Compute normal depth on plane using Mannings equation
        # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
        normal_depth=(ref_flow*mannings_n/(slope**0.5*width))**0.6
        flow_velocity = ref_flow/(width*normal_depth)
        froude_no = flow_velocity/(9.8*normal_depth)**0.5
        if verbose:
            print
            print 'Slope:', slope, 'Mannings n:', mannings_n, 'Velocity:',flow_velocity,'Froude:',froude_no


        #----------------------------------------------------------------------
        # Setup boundary conditions
        #----------------------------------------------------------------------

        Br = Reflective_boundary(domain)              # Solid reflective wall
        #Bo = Dirichlet_boundary([-(slope*length)+normal_depth, 0, 0])  # Outflow stage at normal depth
        def normal_depth_stage_downstream(t):
            return (-slope*length) + normal_depth 
        Bt = Transmissive_Momentum_Set_Stage_boundary(domain=domain,function=normal_depth_stage_downstream)
        domain.set_boundary({'left': Br, 'right': Bt, 'top': Br, 'bottom': Br})

        #----------------------------------------------------------------------
        # Evolve system through time
        #----------------------------------------------------------------------

        x=200.0      # Reference location for flowlines and gauges
        y=10.00
        
        for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
            if verbose :
                print domain.timestepping_statistics()


        #----------------------------------------------------------------------
        # Compute flow thru flowlines ds of inflow and check against normal depth
        #----------------------------------------------------------------------

        # Square on flowline at 200m
        q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]])
        msg = 'Predicted square on flow was %f, should have been %f' % (q, ref_flow)
        if verbose:
            print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
        assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
        
        # 45 degree flowline at 200m
        q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]])
        msg = 'Predicted 45 flow was %f, should have been %f' % (q, ref_flow)
        if verbose:
            print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
        assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         


       # Stage recorder (gauge) in middle of plane at 200m
        w = domain.get_quantity('stage').get_values(interpolation_points=[[x, y]])[0]
        z = domain.get_quantity('elevation').get_values(interpolation_points=[[x, y]])[0]
        domain_depth = w-z
        
        
        msg = 'Predicted depth of flow at 200m was %f, should have been %f' % (normal_depth, domain_depth)
        if verbose:
            print 'Depth at 200m: ANUGA = %f, Mannings = %f' % (domain_depth, normal_depth)
        assert num.allclose(domain_depth,normal_depth, rtol=1.0e-2), msg