source: anuga_work/development/analytical_solutions/Analytical_solution_Yoon_import_mesh.py @ 5029

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

Fixing up directory name

File size: 4.0 KB
Line 
1"""Example of shallow water wave equation analytical solution
2consists of a parabolic profile in a parabolic basin. Analytical
3solutiuon to this problem was derived by Carrier and Greenspan
4and used by Yoon and Chou.
5
6   Copyright 2005
7   Christopher Zoppou, Stephen Roberts, Ole Nielsen
8   ANU, Geoscience Australia
9   
10"""
11
12#---------------
13# Module imports
14import sys
15from os import sep
16sys.path.append('..'+sep+'pyvolution')
17from shallow_water import Domain, Dirichlet_boundary
18from math import sqrt, cos, sin, pi
19from mesh_factory import strang_mesh
20from anuga.pyvolution.util import inside_polygon
21from Numeric import asarray
22from least_squares import Interpolation
23
24#-------------------------------
25# Set up the domain of triangles
26# Strang_domain will search
27# through the file and test to
28# see if there are two or three
29# entries. Two entries are for
30# points and three for triangles.
31points, elements = strang_mesh('yoon_circle.pt')
32domain = Domain(points, elements) 
33
34#----------------
35# Order of scheme
36domain.default_order = 2
37
38domain.smooth = True
39
40#-------------------------------------
41# Provide file name for storing output
42domain.store = False
43domain.format = 'sww'
44domain.filename = 'yoon_mesh_second_order.2'
45print 'Number of triangles = ', len(domain)
46
47#----------------------------------------------------------
48# Decide which quantities are to be stored at each timestep
49domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
50
51#------------------------------------------
52# Reduction operation for get_vertex_values
53from anuga.pyvolution.util import mean
54domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
55
56
57#------------------
58# Initial condition
59print 'Initial condition'
60t = 0.0
61D0 = 1.
62L = 2500.
63R0 = 2000.
64g = 9.81
65
66A = (L**4 - R0**4)/(L**4 + R0**4)
67omega = 2./L*sqrt(2.*g*D0)
68T = pi/omega
69
70#------------------
71# Set bed elevation
72def x_slope(x,y):
73    n = x.shape[0]
74    z = 0*x
75    for i in range(n):
76        r = sqrt(x[i]*x[i] + y[i]*y[i])
77        z[i] = -D0*(1.-r*r/L/L)
78    return z
79domain.set_quantity('elevation', x_slope)
80
81#----------------------------
82# Set the initial water level
83def level(x,y):
84    z = x_slope(x,y)
85    n = x.shape[0]
86    h = 0*x
87    for i in range(n):
88        r = sqrt(x[i]*x[i] + y[i]*y[i])
89        h[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
90                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
91    if h[i] < z[i]:
92        h[i] = z[i]
93    return h   
94domain.set_quantity('stage', level)
95
96
97#---------
98# Boundary
99print 'Boundary conditions'
100domain.set_boundary({'exterior': Dirichlet_boundary([0.0, 0.0, 0.0])})
101
102#---------------------------------------------
103# Find triangle that contains the point points
104# and print to file
105points = [0.,0.]
106for n in range(len(domain.triangles)):
107    t = domain.triangles[n]
108    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
109
110    if inside_polygon(points,tri):
111        print 'Point is within triangle with vertices '+'%s'%tri
112        n_point = n
113t = domain.triangles[n_point]
114tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
115
116filename=domain.filename
117file = open(filename,'w')
118   
119#----------
120# Evolution
121import time
122t0 = time.time()
123for t in domain.evolve(yieldstep = 20.0, finaltime = 3000 ):
124    domain.write_time()
125
126    tri_array = asarray(tri)
127    t_array = asarray([[0,1,2]])
128    interp = Interpolation(tri_array,t_array,[points])
129
130    stage = domain.get_quantity('stage').get_values()[n_point]
131    xmomentum = domain.get_quantity('xmomentum').get_values()[n_point]
132    ymomentum = domain.get_quantity('ymomentum').get_values()[n_point]
133
134    interp_stage = interp.interpolate([[stage[0]],[stage[1]],[stage[2]]])
135    interp_xmomentum = interp.interpolate([[xmomentum[0]],[xmomentum[1]],[xmomentum[2]]])
136    interp_ymomentum = interp.interpolate([[ymomentum[0]],[ymomentum[1]],[ymomentum[2]]])
137   
138    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,interp_stage[0][0],interp_xmomentum[0][0],interp_ymomentum[0][0]) )
139
140file.close()   
141print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.