source: anuga_validation/analytical_solutions/carrier_wave_runup.py @ 7634

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

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

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
from math import sqrt, cos, sin, pi
from numpy import asarray
from time import localtime, strftime, gmtime


#-------------------
# Anuga Imports
#-------------------
from anuga.shallow_water_balanced.swb_domain import Domain, Transmissive_boundary, Reflective_boundary,\
     Dirichlet_boundary, Time_boundary

#from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\
#     Dirichlet_boundary

from anuga.utilities.polygon import inside_polygon, is_inside_triangle
from anuga.abstract_2d_finite_volumes.mesh_factory import strang_mesh


#-------------------
#Convenience functions
#-------------------
def imag(a):
    return a.imag

def real(a):
    return a.real


#--------------------
# Domain
# Strang_domain will search through the file and test to see if there are
# two or three entries. Two entries are for points and three for triangles.

#points, elements = strang_mesh('Run-up.pt')
points, elements = strang_mesh('strang_7389.pt')
domain = Domain(points, elements)

print 'extent ',domain.get_extent()

#----------------
# Order of scheme
# Good compromise between
# limiting and CFL
#---------------
domain.set_default_order(2)
domain.set_timestepping_method(2)
domain.set_beta(0.7)
domain.set_CFL(0.6)


#-------------------
#Set a default tagging
#-------------------
epsilon = 1.0e-12

for id, face in domain.boundary:
    domain.boundary[(id,face)] = 'external'
    if domain.get_vertex_coordinate(id,(face+1)%3)[0] < -200.0+1.0e-10:
        domain.boundary[(id,face)] = 'left'

#---------------------
# Provide file name for storing output
#---------------------
domain.store = True
domain.format = 'sww'

time = strftime('%Y%m%d_%H%M%S',localtime())
output_file= 'carrier_wave_runup_'+time

domain.set_name(output_file)

print "Number of triangles = ", len(domain)

#-----------------------
#Define the boundary condition
#-----------------------
def stage_setup(x,t):
    vh = 0
    alpha = 0.1
    eta = 0.1
    a = 1.5*sqrt(1.+0.9*eta)
    l_0 = 200.
    ii = complex(0,1)
    g = 9.81
    v_0 = sqrt(g*l_0*alpha)
    v1 = 0.

    sigma_max = 100.
    sigma_min = -100.
    for j in range (1,50):
        sigma0 = (sigma_max+sigma_min)/2.
        lambda_prime = 2./a*(t/sqrt(l_0/alpha/g)+v1)
        sigma_prime = sigma0/a
        const = (1.-ii*lambda_prime)**2+sigma_prime**2

        v1 = 8.*eta/a*imag(1./const**(3./2.) \
            -3./4.*(1.-ii*lambda_prime)/const**(5./2.))

        x1 = -v1**2/2.-a**2*sigma_prime**2/16.+eta \
            *real(1.-2.*(5./4.-ii*lambda_prime) \
            /const**(3./2.)+3./2.*(1.-ii*lambda_prime)**2 \
            /const**(5./2.))

        neta1 = x1 + a*a*sigma_prime**2/16.

        v_star1 = v1*v_0
        x_star1 = x1*l_0
        neta_star1 = neta1*alpha*l_0
        stage = neta_star1
        z = stage - x_star1*alpha
        uh = z*v_star1

        if x_star1-x > 0:
            sigma_max = sigma0
        else:
            sigma_min = sigma0

        if abs(abs(sigma0)-100.) < 10:

#     solution does not converge because bed is dry
            stage = 0.
            uh = 0.
            z = 0.

    return [stage, uh, vh]

def boundary_stage(t):
    x = -200
    return stage_setup(x,t)


#---------------------
#Initial condition
#---------------------
print 'Initial condition'
t_star1 = 0.0
slope = -0.1

#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] = -slope*x[i]
    return z

domain.set_quantity('elevation', x_slope)

#Set the water depth
print 'Initial 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], uh, vh = stage_setup(x[i],t_star1)
        h = w[i] - z[i]
        if h < 0:
            h = 0
            w[i] = z[i]
    return w

domain.set_quantity('stage', stage)


#####################
#Set up boundary conditions
Br = Reflective_boundary(domain)
Bw = Time_boundary(domain, boundary_stage)
domain.set_boundary({'left': Bw, 'external': Br})

import time
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

######################
#Evolution

t0 = time.time()


for t in domain.evolve(yieldstep = 1., finaltime = 100):
    domain.write_time()
    print boundary_stage(domain.time)
    if visualize: vis.update()
    
if visualize: vis.evolveFinished()    

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