source: anuga_work/development/analytical solutions/Analytical_solution_Sampson_cross_mesh.py @ 3970

Last change on this file since 3970 was 3514, checked in by duncan, 19 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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.