source: trunk/anuga_validation/validation_tests/periodic_carrier_greenspan.py @ 8427

"""Example of shallow water wave equation analytical solution of the
one-dimensional Carrier and Greenspan wave run-up treated as a two-dimensional solution.

   Copyright 2004
   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
   Geoscience Australia

   Modified by Sudi Mungkasi, Australian National University, 2011
   
Specific methods pertaining to the 2D shallow water equation
are imported from shallow_water
for use with the generic finite volume framework

Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
numerical vector named conserved_quantities.
"""

#-------------------
# Module imports
import sys
import anuga
from math import sqrt, cos, sin, pi
from numpy import asarray
from time import localtime, strftime, gmtime

#-------------------------------------------------------------------------------
# Copy scripts to time stamped output directory and capture screen
# output to file
#-------------------------------------------------------------------------------
time = strftime('%Y%m%d_%H%M%S',localtime())

output_dir = 'periodic_carrier_greenspan_'+time
output_file = 'periodic_carrier_greenspan'

#anuga.copy_code_files(output_dir,__file__)
#start_screen_catcher(output_dir+'_')


#------------------------------------------------------------------------------
# Setup domain
#------------------------------------------------------------------------------
dx = 100.
dy = dx
L = 40000.
W = dx

# structured mesh
points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
domain = anuga.Domain(points, vertices, boundary) 

domain.set_name(output_file)                
domain.set_datadir(output_dir) 

# The following parameters are adapted from Mungkasi & Roberts,
# "Approximations of the Carrier-Greenspan periodic solution"
# International Journal for Numerical Methods in Fluids, 2011
from scipy.special import jn
def j0(x):
    return jn(0.0, x)

def j1(x):
    return jn(1.0, x)
#DIMENSIONAL PARAMETERS
L_0 = 5e3#4               # Length of channel from the origin to shoreline when still
h_0 = 5e2               # Height at origin when the water is still
Tp = 3600.0             # Period of oscillation
a = 250.0                 # Amplitude at origin
g = 9.81                # Gravity
#DIMENSIONLESS PARAMETERS
eps = a/h_0
T = Tp*sqrt(g*h_0)/L_0
A = eps/j0(4.0*pi/T)

#------------------------------------------------------------------------------
# Setup Algorithm
#------------------------------------------------------------------------------
domain.set_default_order(2)
domain.set_timestepping_method(2)
domain.set_beta(0.7)
domain.set_CFL(0.6)


#-----------------------
#Define the boundary condition
#-----------------------
def prescribe(x,t):
    t = t*sqrt(g*h_0)/L_0
    from numpy import zeros
    def fun(q): #Here q=(w, u)
        f = zeros(2)
        f[0] = q[0] + 0.5*q[1]**2.0 - A*j0(4.0*pi/T*(1.0+q[0]-x)**0.5)*cos(2.0*pi/T*(t+q[1]))
        f[1] = q[1] + A*j1(4.0*pi/T*(1.0+q[0]-x)**0.5)*sin(2.0*pi/T*(t+q[1]))/(1+q[0]-x)**0.5
        return f
    from scipy.optimize import fsolve
    w_0, u_0 = fsolve(fun, [0 ,0])   
    return [w_0*h_0,  u_0*sqrt(g*h_0)*h_0,  0.0] #[stage, xmom, ymom]

def boundary_stage(t):
    return prescribe(0.0,t)


#-------------------------------------------------------------------------------
#Initial condition
#-------------------------------------------------------------------------------
#Set bed-elevation and friction(None)
def x_slope(x,y):
    n = x.shape[0]
    z = 0*x
    for i in range(n):
        z[i] = (h_0/L_0)*x[i] - h_0
    return z
domain.set_quantity('elevation', x_slope)

#Set the water depth
def stage(x,y):
    z = x_slope(x,y)
    n = x.shape[0]
    w = 0*x
    for i in range(n):
        w[i] = 0.0
        h = w[i] - z[i]
        if h < 0:
            h = 0
            w[i] = z[i]
    return w
domain.set_quantity('stage', stage)


#-----------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
Br = anuga.Reflective_boundary(domain)
Bt = anuga.Time_boundary(domain, boundary_stage)

# Associate boundary tags with boundary objects
domain.set_boundary({'left': Bt, 'right': Br, 'top': Br, 'bottom': Br})

##visualize = True
##if visualize:
##    from anuga.visualiser import RealtimeVisualiser
##    vis = RealtimeVisualiser(domain)
##    vis.render_quantity_height("elevation", zScale=3.0, offset = 0.01, dynamic=False)
##    vis.render_quantity_height("stage", zScale = 3.0, dynamic=True, opacity = 0.6, wireframe=False)
##    #vis.colour_height_quantity('stage', (lambda q:q['stage'], 1.0, 2.0))
##    vis.colour_height_quantity('stage', (0.4, 0.6, 0.4))
##    vis.start()
##    time.sleep(2.0)

#domain.visualise = True


#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
import time
t0 = time.time()
for t in domain.evolve(yieldstep = 10., finaltime = 5*Tp):
    domain.write_time()
    #print boundary_stage(domain.time)
    #if visualize: vis.update()
    
#if visualize: vis.evolveFinished()    
print 'That took %.2f seconds' %(time.time()-t0)