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

"""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    = 1          # 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


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

def topography(x, y):
    z=-x * slope + abs((y-10)* side_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
#------------------------------------------------------------------------------

# 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

#------------------------------------------------------------------------------
# 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 = Transmissive_boundary(domain)                     # momentum as above bdry
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()


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