source: trunk/anuga_validation/validation_tests/thacker_parabolic/yoon_parabolic_basin_cross_mesh.py @ 8074

Last change on this file since 8074 was 8074, checked in by steve, 13 years ago

Changed name

File size: 4.7 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 Thacker
4and used by Yoon and Chou.
5
6   Copyright 2005
7   Christopher Zoppou, Stephen Roberts, ANU, Geoscience Australia
8
9"""
10
11#---------------
12# Module imports
13#from anuga.shallow_water_balanced.swb_domain import Domain, Transmissive_boundary, Reflective_boundary,\
14#     Dirichlet_boundary
15
16from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\
17     Dirichlet_boundary
18from math import sqrt, cos, sin, pi
19from anuga.interface import rectangular_cross
20from anuga.utilities.polygon import inside_polygon, is_inside_triangle
21from numpy import asarray
22
23
24#-------------------------------
25# Domain
26n = 200
27m = 200
28lenx = 8000.0
29leny = 8000.0
30origin = (-4000.0, -4000.0)
31
32points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
33domain = Domain(points, elements, boundary)
34
35#----------------
36# Order of scheme
37# Good compromise between
38# limiting and CFL
39#---------------
40domain.set_default_order(2)
41domain.set_timestepping_method(2)
42domain.set_beta(0.7)
43domain.set_CFL(0.6)
44
45domain.smooth = True
46
47#-------------------------------------
48# Provide file name for storing output
49domain.store = False
50domain.format = 'sww'
51domain.set_name('yoon_mesh_second_order_cross')
52print 'Number of triangles = ', len(domain)
53
54#----------------------------------------------------------
55# Decide which quantities are to be stored at each timestep
56domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
57
58#------------------------------------------
59# Reduction operation for get_vertex_values
60from anuga.utilities.numerical_tools import mean
61domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
62
63#------------------
64# Initial condition
65print 'Initial condition'
66t = 0.0
67D0 = 1.
68L = 2500.
69R0 = 2000.
70g = 9.81
71
72A = (L**4 - R0**4)/(L**4 + R0**4)
73omega = 2./L*sqrt(2.*g*D0)
74T = pi/omega
75
76#------------------
77# Set bed elevation
78def x_slope(x,y):
79    n = x.shape[0]
80    z = 0*x
81    for i in range(n):
82        r = sqrt(x[i]*x[i] + y[i]*y[i])
83        z[i] = -D0*(1.-r*r/L/L)
84    return z
85domain.set_quantity('elevation', x_slope)
86
87#----------------------------
88# Set the initial water level
89def level(x,y):
90    z = x_slope(x,y)
91    n = x.shape[0]
92    h = 0*x
93    for i in range(n):
94        r = sqrt(x[i]*x[i] + y[i]*y[i])
95        h[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
96                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
97    if h[i] < z[i]:
98        h[i] = z[i]
99    return h
100domain.set_quantity('stage', level)
101
102#---------
103# Boundary
104print 'Boundary conditions'
105R = Reflective_boundary(domain)
106T = Transmissive_boundary(domain)
107D = Dirichlet_boundary([0.0, 0.0, 0.0])
108domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D})
109
110#---------------------------------------------
111# Find triangle that contains the point points
112# and print to file
113
114
115
116points = (0.0, 0.0)
117for n in range(len(domain.triangles)):
118    tri = domain.get_vertex_coordinates(n)
119
120    if is_inside_triangle(points,tri):
121        #print 'Point is within triangle with vertices '+'%s'%tri
122        n_point = n
123
124print 'n_point = ',n_point
125t = domain.triangles[n_point]
126print 't = ', t
127tri = domain.get_vertex_coordinates(n)
128
129filename=domain.get_name()
130file = open(filename,'w')
131
132#----------
133# Evolution
134import time
135t0 = time.time()
136
137time_array = []
138stage_array = []
139Stage     = domain.quantities['stage']
140Xmomentum = domain.quantities['xmomentum']
141Ymomentum = domain.quantities['ymomentum']
142
143for t in domain.evolve(yieldstep = 20.0, finaltime = 17700.0 ):
144    domain.write_time()
145
146    #tri_array = asarray(tri)
147    #t_array = asarray([[0,1,2]])
148    #interp = Interpolation(tri_array,t_array,[points])
149
150
151    stage     = Stage.get_values(location='centroids',indices=[n_point])
152    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])
153    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])
154    #print '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum)
155    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )
156
157    time_array.append(t)
158    stage_array.append(stage)
159
160file.close()
161print 'That took %.2f seconds' %(time.time()-t0)
162
163
164from pylab import *
165ion()
166hold(False)
167plot(time_array, stage_array, 'r.-')
168#title('Gauge %s' %name)
169xlabel('time(s)')
170ylabel('stage (m)')
171#legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
172#savefig(name, dpi = 300)
173
174#raw_input('Next')
175show()
176
177
Note: See TracBrowser for help on using the repository browser.