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

"""ANUGA simulation of simple rip current.

Source: Geometry and wave properties loosely based on those presented in
OBSERVATIONS OF LABORATORY RIP CURRENTS by Brian K. Sapp,
School of Civil and Environmental Engineering
Georgia Institute of Technology
May 2006

I will need to make a version which has the exact same geometry as the Georgia Tech wavetank
if we wish to use a comparison to the results of this study as ANUGA validation as i played
with the geometry somewhat as i completed this model.
"""

#------------------------------------------------------------------------------
# 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 anuga.shallow_water.data_manager import get_mesh_and_quantities_from_file
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, sets the shore at the location defined previously
    z=0.05*(y-location_of_shore)

    #Creates a steeper slope close to the seaward boundary giving a region of deepwater
    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)
    #Creates a steeper slope close to the landward boundary, simulating a beach etc
    #(helps to prevent too much reflection of wave energy off the landward boundary)         
    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)
    
    # Sets up the left hand side of the sandbank
    #amount which it deviates from parallel with the beach is controlled by 'bank_slope'
    #width of the channel (the gap between the two segments of the sandbank) is controlled by 'halfchannelwidth'
    #the height of the sandbar is controlled by 'sandbar'
    #'steepness' provides control over the slope of the soundbar (smaller values give a more rounded shape, if too small will produce peaks and troughs)
    
    for i in range(N):
        ymin = -bank_slope*x[i] + 112 
        ymax = -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)
    
    # Sets up the right hand side of the sandbank
    #changing the sign in y min and y max allows the two halves of the sandbank to form a v shape
    
    for i in range(N):
        ymin = -bank_slope*(x[i]-length/2) - bank_slope*length/2 + 112
        ymax = -bank_slope*(x[i]-length/2) - 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)

# 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

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

n_time_intervals = simulation_length/timestep
print 'There were %i time steps' % n_time_intervals

U = u/n_time_intervals
V = v/n_time_intervals
X = gauges[:,0] 
Y = gauges[:,1] 


print 'Computation took %.2f seconds' % (time.time()-t0)

key_auto_length = (max(V))/5     #makes the key vector a sensible length not sure how to label it with the correct value though

from matplotlib import *
from pylab import *

figure()
Q = quiver(X,Y,U,V)
qk = quiverkey(Q,0.8,0.05,key_auto_length,r'$unknown \frac{m}{s}$',labelpos='E',        #need to get the label to show the value of key_auto_length
               coordinates='figure', fontproperties={'weight': 'bold'})
axis([-10,length + 10, -10, width +10])
title('Simulation of a Rip-Current, Average Velocity Vector Field')

axhline(y=25,color='b')
axhline(y=(location_of_shore),color='r')

x1 = arange(0,55,1)
y1 = -(bank_slope)*x1 + 112
y12 = -(bank_slope)*x1 + 124

x2 = arange(65,length,1)
y2 = -(bank_slope)*x2 + 112
y22 = -(bank_slope)*x2 + 124

plot(x1,y1,x1,y12,x2,y2,x2,y22,color='g')
show()