source: anuga_work/development/classroom/ripcurrent.py @ 7580

"""Simple water flow example using ANUGA
"""

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------
from anuga.interface import create_domain_from_regions
from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
from anuga.shallow_water.shallow_water_domain import Reflective_boundary
from anuga.shallow_water.shallow_water_domain import Time_boundary
from pylab import figure, quiver, show, cos, sin, pi
import numpy
import csv
import time


#------------------------------------------------------------------------------
# Parameters
#------------------------------------------------------------------------------
filename = "WORKING-RIP-LAB_Expt-Geometry_Triangular_Mesh"
location_of_shore = 140 # The position along the y axis of the shorefront
sandbar = 1.2           # Height of sandbar
sealevel = 0            # Height of coast above sea level
steepness = 8000        # Period of sandbar - 
                        # larger number gives smoother slope - longer period 
halfchannelwidth = 5
bank_slope = 0.1
simulation_length = 1
timestep = 1


#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
length = 120
width = 170
seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh
feature_resolution = 1.0
beach_resolution = 10.0

sea_boundary_polygon = [[0,0],[length,0],[length,width],[0,width]]
feature_boundary_polygon = [[0,100],[length,100],[length,150],[0,150]]
beach_interior_polygon = [[0,150],[length,150],[length,width],[0,width]]

meshname = str(filename)+'.msh'

# Interior regions
feature_regions = [[feature_boundary_polygon, feature_resolution],
                   [beach_interior_polygon, beach_resolution]]

domain = create_domain_from_regions(sea_boundary_polygon,
                                    boundary_tags={'bottom': [0],
                                                   'right' : [1],
                                                   'top' :   [2],
                                                   'left':   [3]},
                                    maximum_triangle_area = seafloor_resolution,
                                    mesh_filename = meshname,
                                    interior_regions = feature_regions,
                                    use_cache = True,
                                    verbose = True)
                                    
domain.set_name(filename)                  # Output name
print domain.statistics()


#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
def topography(x,y):
    """Complex topography defined by a function of vectors x and y."""

    # General slope and buildings
    z=0.05*(y-(location_of_shore))
    
    N = len(x)
    for i in range(N):
        if y[i] < 25:
            z[i] = (0.2*(y[i]-25)) + 0.05*(y[i]-(location_of_shore))
    for i in range(N):
        if y[i]>150:
            z[i] = (0.1*(y[i]-150)) + 0.05*(y[i]-(location_of_shore))

    return z


def topography3(x,y):
    z=0*x
    
    N = len(x)
    
    # It would be great with a comment about what this does
    for i in range(N):
        ymin = -1*(bank_slope)*x[i] + 112 
        ymax = -1*(bank_slope)*x[i] + 124
        xmin = 0
        xmax = (length/2)-halfchannelwidth
        if ymin < y[i] < ymax and xmin < x[i]< xmax:
            z[i] += sandbar*((cos((y[i]-118)/steepness)))
    
    # It would be great with a comment about what this does        
    for i in range(N):
        ymin = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+112)
        ymax = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+124)
        xmin = (length/2)+halfchannelwidth
        xmax = 183
        if ymin < y[i] < ymax and xmin < x[i] < xmax:
            z[i] += sandbar*(cos((y[i]-118)/steepness))
            
    return z

domain.set_quantity('elevation', topography)  # Apply base elevation function
domain.add_quantity('elevation', topography3) # Add elevation modification
domain.set_quantity('friction', 0.01)         # Constant friction
domain.set_quantity('stage', 0)               # Constant initial condition at mean sea level


#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
Br = Reflective_boundary(domain)              # Solid reflective wall
Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow

def wave(t):
    """Define wave driving the system
    """
    
    A = 0.4  # Amplitude of wave [m] (wave height)
    T = 5    # Wave period [s]

    if t < 30000000000:
        return [A*sin(2*pi*t/T) + 1, 0, 0] 
    else:
        return [0.0, 0, 0]

Bt = Time_boundary(domain, f=wave)


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


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

# Read in gauge locations for interpolation and convert them to floats
gauge_file = open('New_gauges.csv', 'r')
G = [(float(x[0]), float(x[1])) for x in csv.reader(gauge_file,
                                                    dialect='excel', 
                                                    delimiter=',')]
gauges = numpy.array(G)
number_of_gauges = len(gauges)
print gauges                               

# Allocate space for velocity values
u = numpy.zeros(number_of_gauges)
v = numpy.zeros(number_of_gauges)

t0 = time.time()
for t in domain.evolve(yieldstep = timestep, finaltime = simulation_length):
    print domain.timestepping_statistics()
    S = domain.get_quantity('stage').get_values(interpolation_points=gauges)
    E = domain.get_quantity('elevation').get_values(interpolation_points=gauges)
    depth = S-E

    uh = domain.get_quantity('xmomentum').get_values(interpolation_points=gauges)
    vh = domain.get_quantity('ymomentum').get_values(interpolation_points=gauges)
    u += uh/depth
    v += vh/depth
print 'Evolution took %.2f seconds' % (time.time()-t0)

    
#------------------------------------------------------------------------------
# Post processing
#------------------------------------------------------------------------------    

# Use this

#     from anuga.fit_interpolate.interpolate import Interpolation_function

#     # Get mesh and quantities from sww file
#     X = get_mesh_and_quantities_from_file(filename,
#                                           quantities=quantity_names,
#                                           verbose=verbose)
#     mesh, quantities, time = X

#     # Find all intersections and associated triangles.
#     segments = mesh.get_intersecting_segments(polyline, verbose=verbose)

#     # Get midpoints
#     interpolation_points = segment_midpoints(segments)

#     # Interpolate
#     if verbose:
#         log.critical('Interpolating - total number of interpolation points = %d'
#                      % len(interpolation_points))

#     I = Interpolation_function(time,
#                                quantities,
#                                quantity_names=quantity_names,
#                                vertex_coordinates=mesh.nodes,
#                                triangles=mesh.triangles,
#                                interpolation_points=interpolation_points,
#                                verbose=verbose)

#     return segments, I






n_time_intervals = simulation_length/timestep
u_average = u/n_time_intervals
v_average = v/n_time_intervals

#print "there were", n_time_intervals, "time steps"

#print "sum y velocity", v
#print "average y velocity", v_average
#print "sum x velocity", u
#print "average x velocity", u_average

x_output = file('x_velocity.txt', 'w')
y_output = file('y_velocity.txt', 'w')

#print >> x_output, " "
#print >> y_output, " "

#print >> x_output, u_average
#print >> y_output, v_average


X = gauges[:,0] 
Y = gauges[:,1] 
                              
U = u_average.tolist()
V = v_average.tolist()

#print "U = ", U
#print "U has type", type(U)


figure()
quiver(X,Y,U,V)
show()