source: anuga_work/development/analytical_solutions/oblique_shock.py @ 5162

"""
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('..'+sep+'pyvolution')


from anuga.shallow_water import Domain
from anuga.shallow_water import Transmissive_boundary, Reflective_boundary,\
     Dirichlet_boundary
from anuga.visualiser import RealtimeVisualiser
from math import sqrt, cos, sin, pi
from mesh_factory import oblique_cross


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

points, elements, boundary = oblique_cross(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
domain.visualise=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('stage', expression = 'elevation + %d'%h )

#---------------------------------
# Setup visualization
#---------------------------------
vis = RealtimeVisualiser(domain)
vis.render_quantity_height("elevation", dynamic=False)
vis.render_quantity_height("stage", zScale=10, dynamic=True)
vis.colour_height_quantity('stage', (0.75, 0.5, 0.5))
vis.start()

######################
#Evolution
import time
t0 = time.time()
for t in domain.evolve(yieldstep = 0.5, finaltime = 50):
    domain.write_time()
    vis.update()

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


vis.evolveFinished()
vis.join()

#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")