source: anuga_work/development/culvert_flow/culvert_boyd_channel.py @ 7244

""" 

Ole Check Culvert Routine from Line 258

Although it is Setup as a Culvert with the Opening presented vertically, 
for now the calculation of flow rate is assuming a horizontal hole in the
ground (Fix this Later)

MOST importantly 2 things...
1. How to use the Create Polygon Routine to enquire Depth ( or later energy)
infront of the Culvert

Done (Ole)

2. How to apply the Culvert velocity and thereby Momentum to the Outlet
Ject presented at the Outlet



Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction

This Version CALCULATES the Culvert Velocity and Uses it to establish 
The Culvert Outlet Momentum

The Aim is to define a Flow Transfer function that Simulates a Culvert
by using the Total Available Energy to Drive the Culvert
as Derived by determining the Difference in Total Energy
between 2 Points, Just Up stream and Just Down Stream of the Culvert
away from the influence of the Flow Abstraction etc..

This example includes a Model Topography that shows a
TYPICAL Headwall Configuration

The aim is to change the Culvert Routine to Model more precisely the
abstraction
from a vertical face.

The inflow must include the impact of Approach velocity.
Similarly the Outflow has MOMENTUM Not just Up welling as in the
Horizontal Style
abstraction

"""

#------------------------------------------------------------------------------
# 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 Inflow, General_forcing

from anuga.culvert_flows.culvert_class import Culvert_flow
from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model


from anuga.utilities.polygon import plot_polygons


from math import pi,sqrt

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

length = 40.
width = 5.

#dx = dy = 1           # Resolution: Length of subdivisions on both axes
#dx = dy = .5           # Resolution: Length of subdivisions on both axes
dx = dy = .25           # Resolution: Length of subdivisions on both axes
#dx = dy = .1           # Resolution: Length of subdivisions on both axes

# OLE.... How do I refine the resolution... in the area where I have the Culvert Opening ???......
#  Can I refine in a X & Y Range ???
points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
                                               len1=length, len2=width)
domain = Domain(points, vertices, boundary)   
domain.set_name('culv_dev_HW_Var_Mom')                 # Output name
domain.set_default_order(2)
domain.H0 = 0.01
domain.tight_slope_limiters = True
domain.set_minimum_storable_height(0.001)





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

# Define the topography (land scape)
def topography(x, y):
    """Set up a weir
    
    A culvert will connect either side of an Embankment with a Headwall type structure
    The aim is for the Model to REALISTICALY model flow through the Culvert
    """
    # General Slope of Topography
    z=-x/100
    floorhgt = 5
    embank_hgt=10
    embank_upslope=embank_hgt/5
    embank_dnslope=embank_hgt/2.5
    # Add bits and Pieces to topography
    N = len(x)
    for i in range(N):

        # Sloping Embankment Across Channel
        
        if 0.0 < x[i] < 7.51:
            z[i]=z[i]+5.0
        if 7.5 < x[i] < 10.1:
            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: # Cut Out Segment for Culvert FACE
               z[i]=z[i]+5.0
            else:
               z[i] +=  embank_upslope*(x[i] -5.0)    # Sloping Segment  U/S Face
        if 10.0 < x[i] < 12.1:
            if  2.0 < y[i]  < 3.0: # Cut Out Segment for Culvert (open Channel)
                #z[i] +=  z[i]+5-(x[i]-10)*2              # Sloping Channel in Embankment
                z[i] +=  embank_hgt                 # Flat Crest of Embankment
            else:
                z[i] +=  embank_hgt                 # Flat Crest of Embankment
        if 12.0 < x[i] < 14.5:
            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5: # Cut Out Segment for Culvert FACE
               z[i]=z[i]
            else:
               z[i] += embank_hgt-embank_dnslope*(x[i] -12.0)       # Sloping D/S Face
      
    return z


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



#------------------------------------------------------------------------------
# Setup culvert inlets and outlets in current topography
#------------------------------------------------------------------------------

# Define culvert inlet and outlets
# NEED TO ADD Mannings N as Fixed Value or Function
# Energy Loss Coefficients as Fixed or Function
# Also Set the Shape & Gap Factors fo rthe Enquiry PolyGons
# ALSO Allow the Invert Level to be provided by the USER


culvert1 = Culvert_flow(domain,
                       label='Culvert No. 1',
                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
                       end_point0=[9.0, 2.5], 
                       end_point1=[13.0, 2.5],
                       width=1.20,height=0.75,
                       culvert_routine=boyd_generalised_culvert_model,        
                       number_of_barrels=2,
                       verbose=True)

culvert2 = Culvert_flow(domain,
                       label='Culvert No. 2',
                       description='This culvert is a circular test with d=1.2m',   
                       end_point0=[9.0, 1.5], 
                       end_point1=[30.0, 3.5],
                       diameter=1.20,
                       invert_level0=7,
                       culvert_routine=boyd_generalised_culvert_model,        
                       verbose=True)
                       
domain.forcing_terms.append(culvert1)
domain.forcing_terms.append(culvert2)


#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
#Bi = Dirichlet_boundary([0.5, 0.0, 0.0])          # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!!
Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
Br = Reflective_boundary(domain)              # Solid reflective wall
Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow

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

#------------------------------------------------------------------------------
# Setup Application of  specialised forcing terms
#------------------------------------------------------------------------------

# This is the new element implemented by Ole to allow direct input of Inflow in m^3/s
fixed_flow = Inflow(domain,
                    rate=20.00,
                    center=(2.1, 2.1),
                    radius=1.261566)   #   Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s

#            flow=file_function('Q/QPMF_Rot_Sub13.tms'))        # Read Time Series in  from File
#             flow=lambda t: min(0.01*t, 0.01942))                             # Time Varying Function   Tap turning up         

domain.forcing_terms.append(fixed_flow)


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


for t in domain.evolve(yieldstep = 0.1, finaltime = 20):
    pass
    #if int(domain.time*100) % 100 == 0:
    #    domain.write_time()
    
    #if domain.get_time() >= 4 and tap.flow != 0.0:
    #    print 'Turning tap off'
    #    tap.flow = 0.0
        
    #if domain.get_time() >= 3 and sink.flow < 0.0:
    #    print 'Turning drain on'
    #    sink.flow = -0.8       
    # Close



#------------------------------------------------------------------------------
# Query output
#------------------------------------------------------------------------------

from anuga.shallow_water.data_manager import get_flow_through_cross_section

swwfilename = domain.get_name()+'.sww'  # Output name from script
print swwfilename

polyline = [[17., 0.], [17., 5.]]

time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True)

from pylab import ion, plot
ion()
plot(time, Q)
raw_input('done')