Context Navigation

source: anuga_core/source/anuga/shallow_water/test_tsunami_okada.py @ 5612

Last change on this file since 5612 was 5601, checked in by ole, 16 years ago
Cleaned up in tests and made directories independent on where the test is run from.
File size: 9.1 KB

Line
1	import unittest
2	from Numeric import allclose
3	from tsunami_okada import earthquake_tsunami,Okada_func
4
5	class Test_eq(unittest.TestCase):
6	def setUp(self):
7	pass
8
9	def tearDown(self):
10	pass
11
12
13	def test_Okada_func(self):
14	from os import sep, getenv
15	from Numeric import zeros, Float,allclose
16	import sys
17	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
18	from anuga.shallow_water import Domain
19	from anuga.abstract_2d_finite_volumes.quantity import Quantity
20	from anuga.utilities.system_tools import get_pathname_from_package
21	"""
22	Pick the test you want to do; T= 0 test a point source,
23	T= 1 test single rectangular source, T= 2 test multiple
24	rectangular sources
25	"""
26	# Get path where this test is run
27	path = get_pathname_from_package('anuga.shallow_water')
28
29	# Choose what test to proceed
30	T = 1
31
32
33	if T==0:
34	# Fortran output file
35	filename = path+sep+'fullokada_SP.txt'
36
37	# Initial condition of earthquake for multiple source
38	x0 = 7000.0
39	y0 = 10000.0
40	length = 0
41	width =0
42	strike = 0.0
43	depth = 15.0
44	slip = 10.0
45	dip =15.0
46	rake =90.0
47	ns=1
48	NSMAX=1
49	elif T==1:
50	# Fortran output file
51	filename = path+sep+'fullokada_SS.txt'
52
53	# Initial condition of earthquake for multiple source
54	x0 = 7000.0
55	y0 = 10000.0
56	length = 10.0
57	width =6.0
58	strike = 0.0
59	depth = 15.0
60	slip = 10.0
61	dip =15.0
62	rake =90.0
63	ns=1
64	NSMAX=1
65
66	elif T==2:
67
68	# Fortran output file
69	filename = path+sep+'fullokada_MS.txt'
70
71	# Initial condition of earthquake for multiple source
72	x0 = [7000.0,10000.0]
73	y0 = [10000.0,7000.0]
74	length = [10.0,10.0]
75	width =[6.0,6.0]
76	strike = [0.0,0.0]
77	depth = [15.0,15.0]
78	slip = [10.0,10.0]
79	dip = [15.0,15.0]
80	rake = [90.0,90.0]
81	ns=2
82	NSMAX=2
83
84
85
86	# Get output file from original okada fortran script.
87	# Vertical displacement is listed under tmp.
88	polyline_file=open(filename,'r')
89	lines=polyline_file.readlines()
90	polyline_file.close()
91	tmp=[]
92	stage=[]
93	for line in lines [0:]:
94	field = line.split(' ')
95	z=float(field[2])
96	tmp.append(z)
97
98
99
100
101	#create domain
102	dx = dy = 4000
103	l=100000
104	w=100000
105	#create topography
106	def topography(x,y):
107	el=-1000
108	return el
109
110	points, vertices, boundary = rectangular_cross(int(l/dx), int(w/dy),
111	len1=l, len2=w)
112	domain = Domain(points, vertices, boundary)
113	domain.set_name('test')
114	domain.set_quantity('elevation',topography)
115
116	#create variable with elevation data to implement in okada
117	zrec0 = Quantity(domain)
118	zrec0.set_values(0.0)
119	zrec=zrec0.get_vertex_values(xy=True)
120	# call okada
121	Ts= Okada_func(ns=ns, NSMAX=NSMAX,length=length, width=width, dip=dip, \
122	x0=x0, y0=y0, strike=strike, depth=depth, \
123	slip=slip, rake=rake,zrec=zrec)
124
125	#create a variable to store vertical displacement throughout the domain
126	tsunami = Quantity(domain)
127	tsunami.set_values(Ts)
128
129	# get vertical displacement at each point of the domain respecting
130	# original script's order
131	interpolation_points=[]
132	k=0.0
133	for i in range(0,6):
134	for j in range(0,6):
135	p=j*4000
136	Yt=p
137	Xt=k
138	interpolation_points.append([Xt, Yt])
139
140	k=k+4000
141	Z=tsunami.get_values(interpolation_points=interpolation_points,
142	location='edges')
143
144	stage = -Z # FIXME(Ole): Why the sign flip?
145	# Displacement in fortran code is looking downward
146	#print tmp
147	#print 'hello',stage
148	assert allclose(stage,tmp,atol=1.e-3)
149
150	def test_earthquake_tsunami(self):
151	from os import sep, getenv
152	from Numeric import zeros, Float,allclose
153	import sys
154	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
155	from anuga.shallow_water import Domain
156	from anuga.abstract_2d_finite_volumes.quantity import Quantity
157	from anuga.utilities.system_tools import get_pathname_from_package
158	"""
159	Pick the test you want to do; T= 0 test a point source,
160	T= 1 test single rectangular source, T= 2 test multiple
161	rectangular sources
162	"""
163
164	# Get path where this test is run
165	path= get_pathname_from_package('anuga.shallow_water')
166
167	# Choose what test to proceed
168	T=1
169
170	if T==0:
171	# Fortran output file
172	filename = path+sep+'fullokada_SP.txt'
173
174	# Initial condition of earthquake for multiple source
175	x0 = 7000.0
176	y0 = 10000.0
177	length = 0
178	width =0
179	strike = 0.0
180	depth = 15.0
181	slip = 10.0
182	dip =15.0
183	rake =90.0
184	ns=1
185	NSMAX=1
186	elif T==1:
187	# Fortran output file
188	filename = path+sep+'fullokada_SS.txt'
189
190	# Initial condition of earthquake for multiple source
191	x0 = 7000.0
192	y0 = 10000.0
193	length = 10.0
194	width =6.0
195	strike = 0.0
196	depth = 15.0
197	slip = 10.0
198	dip =15.0
199	rake =90.0
200	ns=1
201	NSMAX=1
202
203	elif T==2:
204
205	# Fortran output file
206	filename = path+sep+'fullokada_MS.txt'
207
208	# Initial condition of earthquake for multiple source
209	x0 = [7000.0,10000.0]
210	y0 = [10000.0,7000.0]
211	length = [10.0,10.0]
212	width =[6.0,6.0]
213	strike = [0.0,0.0]
214	depth = [15.0,15.0]
215	slip = [10.0,10.0]
216	dip = [15.0,15.0]
217	rake = [90.0,90.0]
218	ns=2
219	NSMAX=2
220
221
222
223	# Get output file from original okada fortran script.
224	# Vertical displacement is listed under tmp.
225	polyline_file=open(filename,'r')
226	lines=polyline_file.readlines()
227	polyline_file.close()
228	tmp=[]
229	stage=[]
230	for line in lines [0:]:
231	field = line.split(' ')
232	z=float(field[2])
233	tmp.append(z)
234
235
236	# Create domain
237	dx = dy = 4000
238	l=20000
239	w=20000
240
241	# Create topography
242	def topography(x,y):
243	el=-1000
244	return el
245
246	points, vertices, boundary = rectangular_cross(int(l/dx), int(w/dy),
247	len1=l, len2=w)
248	domain = Domain(points, vertices, boundary)
249	domain.set_name('test')
250	domain.set_quantity('elevation',topography)
251	Ts = earthquake_tsunami(ns=ns,NSMAX=NSMAX,length=length, width=width, strike=strike,\
252	depth=depth,dip=dip, xi=x0, yi=y0,z0=0, slip=slip, rake=rake,\
253	domain=domain, verbose=False)
254
255	# Create a variable to store vertical displacement throughout the domain
256	tsunami = Quantity(domain)
257	tsunami.set_values(Ts)
258	interpolation_points=[]
259
260	#k=0.0
261	#for i in range(0,6):
262	# for j in range(0,6):
263	# p=j*4000
264	# Yt=p
265	# Xt=k
266	# Z=tsunami.get_values(interpolation_points=[[Xt,Yt]]
267	# ,location='edges')
268	# stage.append(-Z[0])
269	# k=k+4000
270	#
271	#assert allclose(stage,tmp,atol=1.e-3)
272
273	# Here's a faster way - try that in the first test
274	interpolation_points=[]
275	k=0.0
276	for i in range(0,6):
277	for j in range(0,6):
278	p=j*4000
279	Yt=p
280	Xt=k
281	interpolation_points.append([Xt, Yt])
282
283	k=k+4000
284	Z=tsunami.get_values(interpolation_points=interpolation_points,
285	location='edges')
286
287	stage = -Z # FIXME(Ole): Why the sign flip?
288	# Displacement in fortran code is looking downward
289	#print 'c est fini'
290	#print tmp
291	#print 'hello',stage
292	assert allclose(stage,tmp,atol=1.e-3)
293
294	#-------------------------------------------------------------
295	if __name__ == "__main__":
296	suite = unittest.makeSuite(Test_eq,'test')
297	runner = unittest.TextTestRunner()
298	runner.run(suite)
299

Note: See TracBrowser for help on using the repository browser.

Download in other formats: