source: inundation/ga/storm_surge/analytical solutions/Analytical solution_oblique_shock.py @ 669

"""Example of shallow water wave equation
consisting of an asymetrical converging channel.

   Copyright 2004
   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
   Geoscience Australia, ANU
   
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 
#

#Were these used?
#import visualise2_chris as visualise
#import Image, ImageGrab

import sys
from os import sep
sys.path.append('..'+s

ep+'pyvolution')


from shallow_water import Domain, Constant_height
from shallow_water import Transmissive_boundary, Reflective_boundary,\
     Dirichlet_boundary

from math import sqrt, cos, sin, pi
from mesh_factory import oblique


######################
# Domain
#
n = 60
m = 80
leny = 30.
lenx = 40.
n = 50
m = 60

points, elements, boundary = oblique(m, n, lenx, leny)
domain = Domain(points, elements, boundary) 

# Order of solver
domain.default_order=2

# Store output
domain.store=True

# Output format
domain.format="sww" #NET.CDF binary format
                    # "dat" for ASCII

# Provide file name for storing output
domain.filename="oblique"

# Visualization smoothing
domain.smooth=True

#######################
#Bed-slope and friction
def x_slope(x, y):
    return 0*x

domain.set_quantity('elevation', x_slope)
domain.set_quantity('friction', 0.0)

######################
# Boundary conditions
#
R = Reflective_boundary(domain)
T = Transmissive_boundary(domain)
D = Dirichlet_boundary([1.0, 8.57, 0.0])

domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R})

######################
#Initial condition
h = 0.5
domain.set_quantity('level', Constant_height(x_slope, h) )


######################
#Evolution
import time
t0 = time.time()
for t in domain.evolve(yieldstep = 0.5, finaltime = 50):
    domain.write_time()
    
print 'That took %.2f seconds' %(time.time()-t0)

#FIXME: Compute average water depth on either side of shock and compare
#to expected values. And also Froude numbers.


#print "saving file?"
#im = ImageGrab.grab()
#im.save("ccube.eps")