source: branches/numpy_anuga_validation/analytical solutions/Analytical solution_oblique_shock_01.py @ 7080

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

from anuga.shallow_water import Domain
from anuga.shallow_water import Transmissive_boundary, Reflective_boundary,\
     Dirichlet_boundary , Transmissive_Momentum_Set_Stage_boundary

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


######################
# Domain
#
leny = 30.
lenx = 40.
f = 2
n = 50*f
m = 80*f
theta = 25
h_bc = 1.5
p_bc = 5.0

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


print 'Number of Elements = ',domain.number_of_elements
# Order of solver
domain.default_order=2
domain.beta_w = 0.8

# Store output
domain.store=True

# Provide file name for storing output
domain.set_name('oblique_cross_%g_%g_%g_%g_%g' % (n,m,theta,h_bc,p_bc))
print domain.get_name()

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


# Visualization smoothing
domain.smooth=True
domain.visualise=False

#######################
#Bed-slope and friction
class ConstantFunctionT:

    def __init__(self,value):
        self.value = value

    def __call__(self,t):
        return self.value


shock_hb = ConstantFunctionT(h_bc)
shock_hh = ConstantFunctionT(h_bc)
shock_pb = ConstantFunctionT(p_bc)
shock_pp = ConstantFunctionT(p_bc)


def x_slope(x, y):
    return 0*x



def shock_h(x,y):
    n = x.shape[0]
    w = 0*x
    for i in range(n):
        if x[i]<0:
            w[i] = shock_hh(0.0)
        else:
            w[i] = 1.0
    return w


def shock_p(x,y):
    n = x.shape[0]
    w = 0*x
    for i in range(n):
        if x[i]<0:
            w[i] = shock_pp(0.0)
        else:
            w[i] = 0.0
    return w





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([shock_hb(0.0) , shock_pb(0.0), 0.0])
Bts = Transmissive_Momentum_Set_Stage_boundary(domain, shock_hb)

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

######################
#Initial condition

domain.set_quantity('stage', shock_h )
domain.set_quantity('xmomentum',shock_p )

class SetValueWhere:

    def __init__(self,quantity,value=0,xrange=0.0):
        self.quantity =  quantity
        self.value = value
        self.xrange = xrange

    def __call__(self,x,y):
        n = x.shape[0]
        w = self.quantity.centroid_values
        for i in range(n):
            if x[i]<self.xrange:
                w[i] = self.value
        return w


Stage = domain.quantities['stage']
Xmom =  domain.quantities['xmomentum']
Ymom =  domain.quantities['ymomentum']

#for id, face in domain.boundary:
#    print '(',id,',',face,') ',domain.boundary[(id,face)]

######################
#Evolution
import time
t0 = time.time()
for t in domain.evolve(yieldstep = 0.1, finaltime = 20.0):
    domain.write_time()
    id = 3399
    print Stage.get_values(location='centroids',indices=[id])
    print Xmom.get_values(location='centroids',indices=[id])
    print Ymom.get_values(location='centroids',indices=[id])
    vstage =  Stage.get_values(location='centroids',indices=[id])
    vxmom  =  Xmom.get_values(location='centroids',indices=[id])
    id = 12719
    print Stage.get_values(location='centroids',indices=[id])
    print Xmom.get_values(location='centroids',indices=[id])
    print Ymom.get_values(location='centroids',indices=[id])
    tclean = 2.0
#    if t > tclean-0.09 and t < tclean+0.01:
#        print 'Cleaning up Initial profile'
#        setstage = SetValueWhere(quantity=Stage,value=vstage[0],xrange=10)
#        setxmom  = SetValueWhere(quantity=Xmom,value=vxmom[0],xrange=10)
#        Stage.set_values(setstage,location='centroids')
#        Xmom.set_values(setxmom,location='centroids')
#        Dnew = Dirichlet_boundary([vstage[0] , vxmom[0], 0.0])
#        domain.set_boundary({'left': Dnew, 'right': T, 'top': R, 'bottom': R})

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