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

Last change on this file since 5030 was 5030, checked in by steve, 16 years ago

Changed a few file names (shortened)

File size: 4.7 KB
Line 
1"""Example of shallow water wave equation analytical solution
2consists of a flat water surface profile in a parabolic basin
3with linear friction. The analytical solution was derived by Sampson in 2002.
4
5   Copyright 2004
6   Christopher Zoppou, Stephen Roberts
7   ANU
8
9Specific methods pertaining to the 2D shallow water equation
10are imported from shallow_water
11for use with the generic finite volume framework
12
13Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
14numerical vector named conserved_quantities.
15"""
16
17######################
18# Module imports
19#
20
21import sys
22from os import sep
23sys.path.append('..'+sep+'pyvolution')
24from anuga.pyvolution.shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\
25     Dirichlet_boundary, gravity, linear_friction
26from math import sqrt, cos, sin, pi, exp
27from anuga.pyvolution.mesh_factory import rectangular_cross
28from anuga.pyvolution.quantity import Quantity
29from anuga.utilities.polygon import inside_polygon
30from Numeric import asarray
31from least_squares import Interpolation
32
33#-------------------------------
34# Domain
35n = 100
36m = 100
37lenx = 10000.0
38leny = 10000.0
39origin = (-5000.0, -5000.0)
40
41points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
42domain = Domain(points, elements, boundary)
43
44#----------------
45# Order of scheme
46domain.default_order = 1
47domain.smooth = True
48
49domain.quantities['linear_friction'] = Quantity(domain)
50
51#---------------------------------------------------------------------------
52# Reconstruct forcing terms with linear friction instead of manning friction
53domain.forcing_terms = []
54domain.forcing_terms.append(gravity)
55domain.forcing_terms.append(linear_friction)
56
57print domain.forcing_terms
58
59#-------------------------------------
60# Provide file name for storing output
61domain.store = True
62domain.format = 'sww'
63domain.filename = 'sampson_first_order_cross'
64
65print 'Number of triangles = ', len(domain)
66
67#-----------------------------------------
68#Reduction operation for get_vertex_values
69from anuga.utilities.numerical_tools import mean
70domain.reduction = mean
71#domain.reduction = min  #Looks better near steep slopes
72
73
74#-----------------
75#Initial condition
76#
77print 'Initial condition'
78t = 0.0
79h0 = 10.
80a = 3000.
81g = 9.81
82tau =0.001
83B = 5
84A = 0
85p = sqrt(8*g*h0/a/a)
86s = sqrt(p*p-tau*tau)/2
87t = 0.
88
89#------------------------------
90#Set bed-elevation and friction
91def x_slope(x,y):
92    n = x.shape[0]
93    z = 0*x
94    for i in range(n):
95        z[i] = h0 - h0*(1.0 -x[i]*x[i]/a/a - y[i]*y[i]/a/a)
96    return z
97
98domain.set_quantity('elevation', x_slope)
99domain.set_quantity('linear_friction', tau)
100
101#-------------------
102#Set the water stage
103def stage(x,y):
104    z = x_slope(x,y)
105    n = x.shape[0]
106    h = 0*x
107    for i in range(n):
108        h[i] = h0-B*B*exp(-tau*t)/2/g-1/g*(exp(-tau*t/2)*(B*s*cos(s*t) \
109               +tau*B/2*sin(s*t)))*x[i] \
110          -1/g*(exp(-tau*t/2)*(B*s*sin(s*t)-tau*B/2*cos(s*t)))*y[i]
111    if h[i] < z[i]:
112        h[i] = z[i]
113    return h
114
115domain.set_quantity('stage', stage)
116
117
118#---------
119# Boundary
120print 'Boundary conditions'
121R = Reflective_boundary(domain)
122T = Transmissive_boundary(domain)
123D = Dirichlet_boundary([0.0, 0.0, 0.0])
124domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D})
125
126#---------------------------------------------
127# Find triangle that contains the point points
128# and print to file
129points = [0.,0.]
130for n in range(len(domain.triangles)):
131    t = domain.triangles[n]
132    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
133
134    if inside_polygon(points,tri):
135        print 'Point is within triangle with vertices '+'%s'%tri
136        n_point = n
137
138print 'n_point = ',n_point
139t = domain.triangles[n_point]
140tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
141
142filename=domain.filename
143file = open(filename,'w')
144
145#----------
146# Evolution
147import time
148t0 = time.time()
149
150Stage     = domain.quantities['stage']
151Xmomentum = domain.quantities['xmomentum']
152Ymomentum = domain.quantities['ymomentum']
153
154for t in domain.evolve(yieldstep = 20.0, finaltime = 10000.0 ):
155    domain.write_time()
156
157    tri_array = asarray(tri)
158    t_array = asarray([[0,1,2]])
159    interp = Interpolation(tri_array,t_array,[points])
160
161
162    stage     = Stage.get_values(location='centroids',indices=[n_point])[0]
163    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])[0]
164    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])[0]
165    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )
166
167file.close()
168print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.