source: trunk/anuga_work/development/gareth/tests/urban_flow/ideal_urban.py @ 8547

"""Idealised urban flow example using ANUGA
Water flowing over a floodplain with buildings 
"""

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------
# Import standard shallow water domain and standard boundaries.
import anuga
import numpy
#from anuga import create_domain_from_regions as create_domain_from_regions
#from anuga.structures.inlet_operator import Inlet_operator

#from balanced_dev import *
#from balanced_dev import create_domain_from_regions as create_domain_from_regions
from anuga_tsunami import *
from anuga_tsunami import create_domain_from_regions as create_domain_from_regions
#------------------------------------------------------------------------------
# Useful parameters for controlling this case
#------------------------------------------------------------------------------

floodplain_length = 35.80 # Model domain length
floodplain_width = 3.60 # Model domain width
floodplain_slope = 0.0 #1./3000.
#chan_initial_depth = 0.8 # Initial depth of water in the channel
#chan_bankfull_depth = 1.0 # Bankfull depth of the channel
#chan_width = 10.0 # Bankfull width of the channel
#bankwidth = 2. # Width of the bank regions -- note that these protrude into the channel
man_n=0.01 # Manning's n
l0 = 0.05 #0.2501 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)

#assert chan_width < floodplain_width, \
#        ' ERROR: Channel width is greater than floodplain width'

#assert bankwidth < chan_width/2., \
#        'ERROR: The bank width must be less than half the channel width'

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

# Define boundary polygon -- in this case clockwise around the proposed boundary
boundary_polygon = [ [0.,0.], 
                     [0., floodplain_width], 
                     [floodplain_length, floodplain_width], 
                     [floodplain_length, 0.], 
                     ]


# Define domain with appropriate boundary conditions
domain = create_domain_from_regions( boundary_polygon, 
                                   boundary_tags={'left': [0],
                                                  'top': [1],
                                                  'right': [2],
                                                  'bottom': [3]},
                                   maximum_triangle_area = 1.0*l0*l0,
                                   minimum_triangle_angle = 28,
                                   mesh_filename = 'urban_flow1.msh',
                                   use_cache=False,
                                   verbose=True)


domain.set_name('urban_flow0p05') # Output name
#domain.extrapolate_velocity_second_order=False
#------------------------------------------------------------------------------
#
# Setup initial conditions
#
#------------------------------------------------------------------------------

## Snippet of code to test the topography function
## We define x,y, then look at z(x,y) to see if it is correct (I use
## pyplot.scatter for this)
## Predefine x,y to be long enough
#increm=0.05
#x = numpy.ones(int(((floodplain_length)/increm+1.0)*((floodplain_width)/increm+1.0) ))
#y = numpy.ones(int(((floodplain_length)/increm+1.0)*((floodplain_width)/increm+1.0) ))
## Fill out x,y
#counter=0
#for i in numpy.arange(0.0,floodplain_length, increm):
#   for j in numpy.arange(0.0, floodplain_width, increm): 
#       x[counter]=i
#       y[counter]=j
#       counter=counter+1
#
#z = topography(x,y)

# Function for topography, based on S. Soares Frazao, Noel and Zech,
# Experiments of dam break flow in the presence of obstables.
# J. Hydraulic Research
def topography(x, y):
    #import pdb
    #pdb.set_trace()
    # Define 'bare-earth' topography
    structure_height=0.8
    elev1= 0.0*y #-y*floodplain_slope 
    # Define constriction
    elev2 = (x>6.9)*(x<7.7)*(y<1.3)*structure_height + \
            (x>6.9)*(x<7.7)*(3.6-y<1.3)*structure_height

    # Add sloping edges
    elev2b = (0.34 - y)*(0.155/0.34)*(y<0.34) + (0.155 - (floodplain_width-y)*(0.155/0.34))*(y>=floodplain_width-0.34)

    # Define building.
    # Do this by using rotated coordinates.
    alpha = 64.0/180.0*numpy.pi
    short_side = 0.4
    long_side=0.8
    # The lower left building corner is at (11.10, 1.75)
    ll = (11.10, 1.75)
    # The lower right building corner is at (11.10, 1.75) + 0.4*(numpy.sin(alpha), -numpy.cos(alpha))
    lr = (ll[0] + short_side*numpy.sin(alpha), ll[1] -short_side*numpy.cos(alpha))
    # The upper left building corner is at (11.10, 1.75) + 0.8*(numpy.cos(64./180.), numpy.sin(64./180.)) 
    ul = (ll[0] + long_side*numpy.cos(alpha), ll[1] + long_side*numpy.sin(alpha))
    c1 = y-numpy.tan(alpha)*x # Coordinate parallel to the long side of the building
    c2 = y +1.0/(numpy.tan(alpha))*x # Coordinate parallel to the short side of the building
    elev3 = (c1 < (ll[1] - numpy.tan(alpha)*ll[0]) )*\
            (c1 > (lr[1] - numpy.tan(alpha)*lr[0]) )*\
            (c2 < (ul[1] + 1.0/numpy.tan(alpha)*ul[0]) )*\
            (c2 > (ll[1] + 1.0/numpy.tan(alpha)*ll[0]) )*\
            structure_height
    return elev1 + elev2 + elev2b + elev3

#Function for stage
def stagetopo(x,y):
    return  (x<7.7)*0.39 + 0.01

domain.set_quantity('elevation', topography) # Use function for elevation
domain.set_quantity('friction', man_n) # Constant friction
domain.set_quantity('stage', stagetopo) # Use function for stage


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

# Trial a set_stage_transmissive_momentum inflow boundary
def inflow_stage_boundary(t):
    return chan_initial_depth - chan_bankfull_depth
# transmissive_momentum_set_stage
Bin_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain=domain, function = inflow_stage_boundary) 

# Standard boundary conditions
Br = anuga.Reflective_boundary(domain) # Solid reflective wall
Bt = anuga.Transmissive_boundary(domain) # Transmissive boundary
# Outflow boundary conditions
def outflow_stage_boundary(t):
    return -floodplain_length*floodplain_slope \
            + chan_initial_depth - chan_bankfull_depth
Bout_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain, function = outflow_stage_boundary) 

# Define boundary conditions
domain.set_boundary({'left': Br, 
                     'right': Br, 
                     'top': Br, 
                     'bottom': Br})


# To monitor the boundary conditions, let's find centroid points representing sub and super critical flows.
#inpt = numpy.argmin( ( abs(domain.centroid_coordinates[:,0] - floodplain_width/2.) + \
#        abs(domain.centroid_coordinates[:,1] - 0.0)  ) )
#outpt = numpy.argmin( ( abs(domain.centroid_coordinates[:, 0] - floodplain_width/2.) +\
#        abs(domain.centroid_coordinates[:, 1] - floodplain_length)  ) )

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

for t in domain.evolve(yieldstep=0.3, finaltime=40.0):
    print domain.timestepping_statistics()