source: anuga_work/development/analytical_solutions/sampson_cross_mesh.py @ 5599

"""Example of shallow water wave equation analytical solution
consists of a flat water surface profile in a parabolic basin
with linear friction. The analytical solution was derived by Sampson in 2002.

   Copyright 2004
   Christopher Zoppou, Stephen Roberts
   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
#

import sys
from os import sep
sys.path.append('..'+sep+'pyvolution')
from anuga.pyvolution.shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\
     Dirichlet_boundary, gravity, linear_friction
from math import sqrt, cos, sin, pi, exp
from anuga.pyvolution.mesh_factory import rectangular_cross
from anuga.pyvolution.quantity import Quantity
from anuga.utilities.polygon import inside_polygon
from Numeric import asarray
from least_squares import Interpolation

#-------------------------------
# Domain
n = 100
m = 100
lenx = 10000.0
leny = 10000.0
origin = (-5000.0, -5000.0)

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

#----------------
# Order of scheme
domain.default_order = 1
domain.smooth = True

domain.quantities['linear_friction'] = Quantity(domain)

#---------------------------------------------------------------------------
# Reconstruct forcing terms with linear friction instead of manning friction
domain.forcing_terms = []
domain.forcing_terms.append(gravity)
domain.forcing_terms.append(linear_friction)

print domain.forcing_terms

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

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

#-----------------------------------------
#Reduction operation for get_vertex_values
from anuga.utilities.numerical_tools import mean
domain.reduction = mean
#domain.reduction = min  #Looks better near steep slopes


#-----------------
#Initial condition
#
print 'Initial condition'
t = 0.0
h0 = 10.
a = 3000.
g = 9.81
tau =0.001
B = 5
A = 0
p = sqrt(8*g*h0/a/a)
s = sqrt(p*p-tau*tau)/2
t = 0.

#------------------------------
#Set bed-elevation and friction
def x_slope(x,y):
    n = x.shape[0]
    z = 0*x
    for i in range(n):
        z[i] = h0 - h0*(1.0 -x[i]*x[i]/a/a - y[i]*y[i]/a/a)
    return z

domain.set_quantity('elevation', x_slope)
domain.set_quantity('linear_friction', tau)

#-------------------
#Set the water stage
def stage(x,y):
    z = x_slope(x,y)
    n = x.shape[0]
    h = 0*x
    for i in range(n):
        h[i] = h0-B*B*exp(-tau*t)/2/g-1/g*(exp(-tau*t/2)*(B*s*cos(s*t) \
               +tau*B/2*sin(s*t)))*x[i] \
          -1/g*(exp(-tau*t/2)*(B*s*sin(s*t)-tau*B/2*cos(s*t)))*y[i]
    if h[i] < z[i]:
        h[i] = z[i]
    return h

domain.set_quantity('stage', stage)


#---------
# Boundary
print 'Boundary conditions'
R = Reflective_boundary(domain)
T = Transmissive_boundary(domain)
D = Dirichlet_boundary([0.0, 0.0, 0.0])
domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D})

#---------------------------------------------
# Find triangle that contains the point points
# and print to file
points = [0.,0.]
for n in range(len(domain.triangles)):
    t = domain.triangles[n]
    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]

    if inside_polygon(points,tri):
        print 'Point is within triangle with vertices '+'%s'%tri
        n_point = n

print 'n_point = ',n_point
t = domain.triangles[n_point]
tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]

filename=domain.filename
file = open(filename,'w')

#----------
# Evolution
import time
t0 = time.time()

Stage     = domain.quantities['stage']
Xmomentum = domain.quantities['xmomentum']
Ymomentum = domain.quantities['ymomentum']

for t in domain.evolve(yieldstep = 20.0, finaltime = 10000.0 ):
    domain.write_time()

    tri_array = asarray(tri)
    t_array = asarray([[0,1,2]])
    interp = Interpolation(tri_array,t_array,[points])


    stage     = Stage.get_values(location='centroids',indices=[n_point])[0]
    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])[0]
    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])[0]
    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )

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