source: anuga_work/development/ted_rigby/test_flowhydraulics_in_vee_using_inflow.py @ 7475

"""test_flow_in_vee_using_inflow.py
Test the ability of ANUGA to 
a) maintain the flow rate input into the channel at the upstream end
b) replicate normal depth half way down the channel after half an hour
c) replicate a typical velocity profile across the channel after half an hour

The vee channel is 20m wide by 1000m 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 in the invert of (y=10m)and half way down the channel (x=500m) 
This script explores flow behaviour predicted by ANUGA and compares it with known solutions

This test inputs flow using the 'inflow' function .

"""

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     = 1000.
width      = 20.
side_slope = 1.0
dx = dy    = 1          # Resolution: of points (elev) grid on both axes
ref_flow   = 50.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_flowhydraulics_in_vee_using_inflow')  # 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')   # zero depth initial condition


#------------------------------------------------------------------------------
# Compute normal depth flow characteristics for 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

print 'Mannings solution for normal depth flow in vee is; '
print '    Channel sideslopes = ', side_slope
print '    Channel n          = ', mannings_n
print '    Channel slope      = ', slope
print '    Reference flow     = ', ref_flow
print '    Normal depth       = ', normal_depth
print '    Flow area          = ', area
print '    Wet perimeter      = ', wet_perimeter
print '    Hydraulic radius   = ', hyd_radius
print '    Avg velocity       = ', flow_velocity
print '    Flow width         = ', top_width

#------------------------------------------------------------------------------
# Calculate inflow and append to domain forcing terms
#------------------------------------------------------------------------------

# Fixed Flowrate onto Area 
fixed_inflow = Inflow(domain,
                      center=(10.0, 10.0),
                      radius=5.00,
                      rate=ref_flow)   
domain.forcing_terms.append(fixed_inflow)

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


Br = Reflective_boundary(domain)                       # Solid reflective wall on edge of sides
# set outflow depth to mannings normal depth for flow to avoid drawdown impact upstream
def outflow_depth(t):
    return (-slope*length) + normal_depth 
Bo = Transmissive_Momentum_Set_Stage_boundary(domain=domain,function=outflow_depth)

domain.set_boundary({'left': Br, '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 across flowlines -v- ref_flow 
#------------------------------------------------------------------------------

# Check flow across flowline at 500m (where  it should have normalised)
q=domain.get_flow_through_cross_section([[500.0,0.0],[500.0,width]])
print '90 degree flowline half way down channel: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)

# Check flow across flowline at 998m (wimmediately US of outflow bdry)
q=domain.get_flow_through_cross_section([[998.0,0.0],[998.0,width]])
print '90 degree flowline just US of outflow bdry: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)

#------------------------------------------------------------------------------
# Record temporal flow characteristics just ds inflow bdry
#------------------------------------------------------------------------------

# Store temporal pattern of gauges at head/middle/bottom of channel in csv files 
sww2csv_gauges(sww_file  = 'test_flowhydraulics_in_vee_using_inflow.sww',
               gauge_file= 'gauge_locations.csv',
               out_name  = 'inflow_',
               quantities=['stage','speed','elevation'],
               verbose=True)