source: anuga_work/development/ted_rigby/test_flow_in_vee_basic.py @ 7741

"""test_flow_in_vee.py
Test the ability of ANUGA to 
a) maintain the flow rate input into the channel at the upstream end

The vee channel is 20m wide by 400m long with side slopes of 1:1 
The equations for dn are generalised with side slope as 1 V in Side_Slope H
The test location is half way down the channel X=200m 
This script explores flow  as given by ANUGA and compares it with 
the ref_flow 
"""

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 Transmissive_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
#------------------------------------------------------------------------------

finaltime = 1800.0    # 1/2 hour

length     = 400.
width      = 20.
side_slope = 1.0
dx = dy    = 2          # Resolution: of points (elev) grid on both axes
ref_flow   = 5.0        # Inflow at head of vee
slope      = 1.0/1000
mannings_n = 0.150

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



domain = Domain(points, vertices, boundary)   
domain.set_name('test_flow_in_vee')              # Output name
print domain.statistics()


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

# Compute normal depth, and velocity for the vee using Mannings equation

normal_depth  = ( ref_flow*mannings_n* ( (2.0*( (1.0+(1.0/side_slope)**2.0)**0.5) )**0.6667 ) / (side_slope*(slope**0.5))  )**0.375
area          = (normal_depth**2.0) *side_slope
wet_perimeter = 2*normal_depth*((1+(side_slope)**2.0)**0.5)
hyd_radius    = area/wet_perimeter
flow_velocity = ref_flow/area
top_width     = 2*side_slope*normal_depth

# Compute xmomentum to input into the ANUGA inflow boundary
# assuming momentum is applied uniformly to inflow boundary by ANUGA

#xmomentum = ref_flow/top_width  # this is all very unrealistic but should replicate ref_flow in ANUGA
xmomentum = ref_flow/width       # This only applies where channel is flat. 

print 'ref_flow', ref_flow
print 'xmomentum', xmomentum
print 'normal_depth', normal_depth
print 'width', width



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

def topography(x, y):
    z = -x * slope + abs((y-10)*side_slope)
    
    # Let both ends of domain be square so that the applied xmomentum 
    # represents flow correctly 
    for i in range(len(x)):
        if x[i] < 20 or x[i] > 380:
            z[i] = -x[i] * slope
            
    return z
    
    
def initial_stage(x, y):
    """Set constant depth based on calculated normal flow
    """

    w = -x*slope + normal_depth
    return w
    
    
domain.set_quantity('elevation', topography)  # Use function for elevation
domain.set_quantity('friction', mannings_n)   # Constant friction of conc surface   
domain.set_quantity('stage', initial_stage) 
domain.set_quantity('xmomentum', xmomentum) # Initial momentum commensurate with normal flow



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


Br = Reflective_boundary(domain)                       # Solid reflective wall on edge of sides
Bi = Dirichlet_boundary([normal_depth,xmomentum,0.0])  # inflow at top
Bo = Dirichlet_boundary([-length*slope+normal_depth,xmomentum,0.0])  # Outflow at bottom

domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})

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

for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
    if verbose :
        print domain.timestepping_statistics()
     
     
    # Check flow across flowline at 200m at every yield step
    q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,width]])
    print '90 degree flowline: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)     


#------------------------------------------------------------------------------
# Compute flow in vee cf inflow
#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Compute flow across flowline -v- ref_flow 
#------------------------------------------------------------------------------

# Check flow across flowline at 200m
q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,width]])
print '90 degree flowline: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)