Context Navigation

source: inundation/ga/storm_surge/pyvolution/test_shallow_water.py @ 390

Last change on this file since 390 was 272, checked in by ole, 21 years ago
Coded update as C implementation and verified improvement. Removed oldstyle obsolete balanced limiter
File size: 71.3 KB

Line
1	#!/usr/bin/env python
2
3	import unittest
4	from math import sqrt, pi
5
6	from config import g, epsilon
7	from Numeric import allclose, array, zeros, ones, Float
8	from shallow_water import *
9
10
11	class TestCase(unittest.TestCase):
12	def setUp(self):
13	pass
14
15	def tearDown(self):
16	pass
17
18	def test_rotate(self):
19	normal = array([0.0,-1.0])
20
21	q = array([1.0,2.0,3.0])
22
23	r = rotate(q, normal, direction = 1)
24	assert r[0] == 1
25	assert r[1] == -3
26	assert r[2] == 2
27
28	w = rotate(r, normal, direction = -1)
29	assert allclose(w, q)
30
31
32
33	def test_flux_zero_case(self):
34	ql = zeros( 3, Float )
35	qr = zeros( 3, Float )
36	normal = zeros( 2, Float )
37	zl = zr = 0.
38	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
39
40	assert allclose(flux, [0,0,0])
41	assert max_speed == 0.
42
43	def test_flux_constants(self):
44	w = 2.0
45
46	normal = array([1.,0])
47	ql = array([w, 0, 0])
48	qr = array([w, 0, 0])
49	zl = zr = 0.
50	h = w - (zl+zr)/2
51
52	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
53
54	assert allclose(flux, [0., 0.5gh**2, 0.])
55	assert max_speed == sqrt(g*h)
56
57
58	def test_flux1(self):
59	#Use data from previous version of pyvolution
60	normal = array([1.,0])
61	ql = array([-0.2, 2, 3])
62	qr = array([-0.2, 2, 3])
63	zl = zr = -0.5
64	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
65
66	assert allclose(flux, [2.,13.77433333, 20.])
67	assert allclose(max_speed, 8.38130948661)
68
69
70	def test_flux2(self):
71	#Use data from previous version of pyvolution
72	normal = array([0., -1.])
73	ql = array([-0.075, 2, 3])
74	qr = array([-0.075, 2, 3])
75	zl = zr = -0.375
76	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
77
78	assert allclose(flux, [-3.,-20.0, -30.441])
79	assert allclose(max_speed, 11.7146428199)
80
81	def test_flux3(self):
82	#Use data from previous version of pyvolution
83	normal = array([-sqrt(2)/2, sqrt(2)/2])
84	ql = array([-0.075, 2, 3])
85	qr = array([-0.075, 2, 3])
86	zl = zr = -0.375
87	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
88
89	assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
90	assert allclose(max_speed, 4.0716654239)
91
92	def test_flux4(self):
93	#Use data from previous version of pyvolution
94	normal = array([-sqrt(2)/2, sqrt(2)/2])
95	ql = array([-0.34319278, 0.10254161, 0.07273855])
96	qr = array([-0.30683287, 0.1071986, 0.05930515])
97	zl = zr = -0.375
98	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
99
100	assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
101	assert allclose(max_speed, 1.31414103233)
102
103	def test_sw_domain_simple(self):
104	a = [0.0, 0.0]
105	b = [0.0, 2.0]
106	c = [2.0,0.0]
107	d = [0.0, 4.0]
108	e = [2.0, 2.0]
109	f = [4.0,0.0]
110
111	points = [a, b, c, d, e, f]
112	#bac, bce, ecf, dbe, daf, dae
113	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
114
115	domain = Domain(points, vertices)
116	domain.check_integrity()
117
118	for name in ['level', 'xmomentum', 'ymomentum',
119	'elevation', 'friction']:
120	assert domain.quantities.has_key(name)
121
122
123	assert domain.get_conserved_quantities(0, edge=1) == 0.
124
125
126	def test_boundary_conditions(self):
127
128	a = [0.0, 0.0]
129	b = [0.0, 2.0]
130	c = [2.0,0.0]
131	d = [0.0, 4.0]
132	e = [2.0, 2.0]
133	f = [4.0,0.0]
134
135	points = [a, b, c, d, e, f]
136	#bac, bce, ecf, dbe
137	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
138	boundary = { (0, 0): 'Third',
139	(0, 2): 'First',
140	(2, 0): 'Second',
141	(2, 1): 'Second',
142	(3, 1): 'Second',
143	(3, 2): 'Third'}
144
145
146	domain = Domain(points, vertices, boundary)
147	domain.check_integrity()
148
149
150	domain.set_quantity('level', [[1,2,3], [5,5,5],
151	[0,0,9], [-6, 3, 3]])
152
153	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
154	[3,3,3], [4, 4, 4]])
155
156	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
157	[30,30,30], [40, 40, 40]])
158
159
160	D = Dirichlet_boundary([5,2,1])
161	T = Transmissive_boundary(domain)
162	R = Reflective_boundary(domain)
163	domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
164
165	domain.update_boundary()
166
167	#Level
168	assert domain.quantities['level'].boundary_values[0] == 2.5
169	assert domain.quantities['level'].boundary_values[0] ==\
170	domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
171	assert domain.quantities['level'].boundary_values[1] == 5. #Dirichlet
172	assert domain.quantities['level'].boundary_values[2] ==\
173	domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
174	assert domain.quantities['level'].boundary_values[3] ==\
175	domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
176	assert domain.quantities['level'].boundary_values[4] ==\
177	domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
178	assert domain.quantities['level'].boundary_values[5] ==\
179	domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
180
181	#Xmomentum
182	assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
183	assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
184	assert domain.quantities['xmomentum'].boundary_values[2] ==\
185	domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
186	assert domain.quantities['xmomentum'].boundary_values[3] ==\
187	domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
188	assert domain.quantities['xmomentum'].boundary_values[4] ==\
189	domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
190	assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective
191
192	#Ymomentum
193	assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
194	assert domain.quantities['ymomentum'].boundary_values[1] == 1. #Dirichlet
195	assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
196	assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
197	assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
198	assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
199
200
201
202
203	def test_compute_fluxes0(self):
204	#Do a full triangle and check that fluxes cancel out for
205	#the constant level case
206
207	a = [0.0, 0.0]
208	b = [0.0, 2.0]
209	c = [2.0,0.0]
210	d = [0.0, 4.0]
211	e = [2.0, 2.0]
212	f = [4.0,0.0]
213
214	points = [a, b, c, d, e, f]
215	#bac, bce, ecf, dbe
216	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
217
218	domain = Domain(points, vertices)
219	domain.set_quantity('level', [[2,2,2], [2,2,2],
220	[2,2,2], [2,2,2]])
221	domain.check_integrity()
222
223	zl=zr=0. #Assume flat bed
224
225	#Flux across right edge of volume 1
226	normal = domain.get_normal(1,0)
227	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
228	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
229	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
230
231	#Check that flux seen from other triangles is inverse
232	tmp = qr; qr=ql; ql=tmp
233	normal = domain.get_normal(2,2)
234	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
235	assert allclose(flux + flux0, 0.)
236
237	#Flux across upper edge of volume 1
238	normal = domain.get_normal(1,1)
239	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
240	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
241	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
242
243	#Check that flux seen from other triangles is inverse
244	tmp = qr; qr=ql; ql=tmp
245	normal = domain.get_normal(3,0)
246	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
247	assert allclose(flux + flux1, 0.)
248
249	#Flux across lower left hypotenuse of volume 1
250	normal = domain.get_normal(1,2)
251	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
252	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
253	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
254
255	#Check that flux seen from other triangles is inverse
256	tmp = qr; qr=ql; ql=tmp
257	normal = domain.get_normal(0,1)
258	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
259	assert allclose(flux + flux2, 0.)
260
261
262	#Scale by edgelengths, add up anc check that total flux is zero
263	e0 = domain.edgelengths[1, 0]
264	e1 = domain.edgelengths[1, 1]
265	e2 = domain.edgelengths[1, 2]
266
267	assert allclose(e0flux0+e1flux1+e2*flux2, 0.)
268
269	#Now check that compute_flux yields zeros as well
270	domain.compute_fluxes()
271
272	for name in ['level', 'xmomentum', 'ymomentum']:
273	#print name, domain.quantities[name].explicit_update
274	assert allclose(domain.quantities[name].explicit_update[1], 0)
275
276
277
278	def test_compute_fluxes1(self):
279	#Use values from previous version
280
281	a = [0.0, 0.0]
282	b = [0.0, 2.0]
283	c = [2.0,0.0]
284	d = [0.0, 4.0]
285	e = [2.0, 2.0]
286	f = [4.0,0.0]
287
288	points = [a, b, c, d, e, f]
289	#bac, bce, ecf, dbe
290	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
291
292	domain = Domain(points, vertices)
293	val0 = 2.+2.0/3
294	val1 = 4.+4.0/3
295	val2 = 8.+2.0/3
296	val3 = 2.+8.0/3
297
298	domain.set_quantity('level', [[val0, val0, val0], [val1, val1, val1],
299	[val2, val2, val2], [val3, val3, val3]])
300	domain.check_integrity()
301
302	zl=zr=0. #Assume flat bed
303
304	#Flux across right edge of volume 1
305	normal = domain.get_normal(1,0)
306	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
307	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
308	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
309	assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
310	assert allclose(max_speed, 9.21592824046)
311
312
313	#Flux across upper edge of volume 1
314	normal = domain.get_normal(1,1)
315	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
316	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
317	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
318	assert allclose(flux1, [2.4098563, 0., 123.04444444])
319	assert allclose(max_speed, 7.22956891292)
320
321	#Flux across lower left hypotenuse of volume 1
322	normal = domain.get_normal(1,2)
323	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
324	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
325	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
326
327	assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
328	assert allclose(max_speed, 7.22956891292)
329
330	#Scale, add up and check that compute_fluxes is correct for vol 1
331	e0 = domain.edgelengths[1, 0]
332	e1 = domain.edgelengths[1, 1]
333	e2 = domain.edgelengths[1, 2]
334
335	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
336	assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
337
338
339	domain.compute_fluxes()
340
341	#assert allclose(total_flux, domain.explicit_update[1,:])
342	for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
343	assert allclose(total_flux[i],
344	domain.quantities[name].explicit_update[1])
345
346	#assert allclose(domain.explicit_update, [
347	# [0., -69.68888889, -69.68888889],
348	# [-0.68218178, -166.6, -35.93333333],
349	# [-111.77316251, 69.68888889, 0.],
350	# [-35.68522449, 0., 69.68888889]])
351
352	assert allclose(domain.quantities['level'].explicit_update,
353	[0., -0.68218178, -111.77316251, -35.68522449])
354	assert allclose(domain.quantities['xmomentum'].explicit_update,
355	[-69.68888889, -166.6, 69.68888889, 0])
356	assert allclose(domain.quantities['ymomentum'].explicit_update,
357	[-69.68888889, -35.93333333, 0., 69.68888889])
358
359
360	#assert allclose(domain.quantities[name].explicit_update
361
362
363
364
365
366	def test_compute_fluxes2(self):
367	#Random values, incl momentum
368
369	a = [0.0, 0.0]
370	b = [0.0, 2.0]
371	c = [2.0,0.0]
372	d = [0.0, 4.0]
373	e = [2.0, 2.0]
374	f = [4.0,0.0]
375
376	points = [a, b, c, d, e, f]
377	#bac, bce, ecf, dbe
378	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
379
380	domain = Domain(points, vertices)
381	val0 = 2.+2.0/3
382	val1 = 4.+4.0/3
383	val2 = 8.+2.0/3
384	val3 = 2.+8.0/3
385
386	zl=zr=0 #Assume flat zero bed
387
388	domain.set_quantity('elevation', zl*ones( (4,3) ))
389
390
391	domain.set_quantity('level', [[val0, val0-1, val0-2],
392	[val1, val1+1, val1],
393	[val2, val2-2, val2],
394	[val3-0.5, val3, val3]])
395
396	domain.set_quantity('xmomentum',
397	[[1, 2, 3], [3, 4, 5],
398	[1, -1, 0], [0, -2, 2]])
399
400	domain.set_quantity('ymomentum',
401	[[1, -1, 0], [0, -3, 2],
402	[0, 1, 0], [-1, 2, 2]])
403
404
405	domain.check_integrity()
406
407
408
409	#Flux across right edge of volume 1
410	normal = domain.get_normal(1,0)
411	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
412	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
413	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
414
415	#Flux across upper edge of volume 1
416	normal = domain.get_normal(1,1)
417	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
418	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
419	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
420
421	#Flux across lower left hypotenuse of volume 1
422	normal = domain.get_normal(1,2)
423	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
424	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
425	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
426
427	#Scale, add up and check that compute_fluxes is correct for vol 1
428	e0 = domain.edgelengths[1, 0]
429	e1 = domain.edgelengths[1, 1]
430	e2 = domain.edgelengths[1, 2]
431
432	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
433
434
435	domain.compute_fluxes()
436	for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
437	assert allclose(total_flux[i],
438	domain.quantities[name].explicit_update[1])
439	#assert allclose(total_flux, domain.explicit_update[1,:])
440
441
442	def test_compute_fluxes3(self):
443	#Random values, incl momentum
444
445	a = [0.0, 0.0]
446	b = [0.0, 2.0]
447	c = [2.0,0.0]
448	d = [0.0, 4.0]
449	e = [2.0, 2.0]
450	f = [4.0,0.0]
451
452	points = [a, b, c, d, e, f]
453	#bac, bce, ecf, dbe
454	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
455
456	domain = Domain(points, vertices)
457	val0 = 2.+2.0/3
458	val1 = 4.+4.0/3
459	val2 = 8.+2.0/3
460	val3 = 2.+8.0/3
461
462	zl=zr=-3.75 #Assume constant bed (must be less than level)
463	domain.set_quantity('elevation', zl*ones( (4,3) ))
464
465
466	domain.set_quantity('level', [[val0, val0-1, val0-2],
467	[val1, val1+1, val1],
468	[val2, val2-2, val2],
469	[val3-0.5, val3, val3]])
470
471	domain.set_quantity('xmomentum',
472	[[1, 2, 3], [3, 4, 5],
473	[1, -1, 0], [0, -2, 2]])
474
475	domain.set_quantity('ymomentum',
476	[[1, -1, 0], [0, -3, 2],
477	[0, 1, 0], [-1, 2, 2]])
478
479
480	domain.check_integrity()
481
482
483
484	#Flux across right edge of volume 1
485	normal = domain.get_normal(1,0)
486	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
487	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
488	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
489
490	#Flux across upper edge of volume 1
491	normal = domain.get_normal(1,1)
492	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
493	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
494	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
495
496	#Flux across lower left hypotenuse of volume 1
497	normal = domain.get_normal(1,2)
498	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
499	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
500	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
501
502	#Scale, add up and check that compute_fluxes is correct for vol 1
503	e0 = domain.edgelengths[1, 0]
504	e1 = domain.edgelengths[1, 1]
505	e2 = domain.edgelengths[1, 2]
506
507	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
508
509	domain.compute_fluxes()
510	for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
511	assert allclose(total_flux[i],
512	domain.quantities[name].explicit_update[1])
513
514
515
516	def test_catching_negative_heights(self):
517
518	a = [0.0, 0.0]
519	b = [0.0, 2.0]
520	c = [2.0,0.0]
521	d = [0.0, 4.0]
522	e = [2.0, 2.0]
523	f = [4.0,0.0]
524
525	points = [a, b, c, d, e, f]
526	#bac, bce, ecf, dbe
527	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
528
529	domain = Domain(points, vertices)
530	val0 = 2.+2.0/3
531	val1 = 4.+4.0/3
532	val2 = 8.+2.0/3
533	val3 = 2.+8.0/3
534
535	zl=zr=4 #Too large
536	domain.set_quantity('elevation', zl*ones( (4,3) ))
537	domain.set_quantity('level', [[val0, val0-1, val0-2],
538	[val1, val1+1, val1],
539	[val2, val2-2, val2],
540	[val3-0.5, val3, val3]])
541
542	#Should fail
543	try:
544	domain.check_integrity()
545	except:
546	pass
547
548	#####################################################
549	def test_gravity(self):
550	#Assuming no friction
551
552	from config import g
553
554	a = [0.0, 0.0]
555	b = [0.0, 2.0]
556	c = [2.0, 0.0]
557	d = [0.0, 4.0]
558	e = [2.0, 2.0]
559	f = [4.0, 0.0]
560
561	points = [a, b, c, d, e, f]
562	#bac, bce, ecf, dbe
563	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
564
565	domain = Domain(points, vertices)
566
567	#Set up for a gradient of (3,0) at mid triangle
568	def slope(x, y):
569	return 3*x
570
571	h = 0.1
572	def level(x,y):
573	return slope(x,y)+h
574
575	domain.set_quantity('elevation', slope)
576	domain.set_quantity('level', level)
577
578	for name in domain.conserved_quantities:
579	assert allclose(domain.quantities[name].explicit_update, 0)
580	assert allclose(domain.quantities[name].semi_implicit_update, 0)
581
582	domain.compute_forcing_terms()
583
584	assert allclose(domain.quantities['level'].explicit_update, 0)
585	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
586	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
587
588
589
590	#####################################################
591	def test_first_order_extrapolator_const_z(self):
592
593	a = [0.0, 0.0]
594	b = [0.0, 2.0]
595	c = [2.0, 0.0]
596	d = [0.0, 4.0]
597	e = [2.0, 2.0]
598	f = [4.0, 0.0]
599
600	points = [a, b, c, d, e, f]
601	#bac, bce, ecf, dbe
602	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
603
604	domain = Domain(points, vertices)
605	val0 = 2.+2.0/3
606	val1 = 4.+4.0/3
607	val2 = 8.+2.0/3
608	val3 = 2.+8.0/3
609
610	zl=zr=-3.75 #Assume constant bed (must be less than level)
611	domain.set_quantity('elevation', zl*ones( (4,3) ))
612	domain.set_quantity('level', [[val0, val0-1, val0-2],
613	[val1, val1+1, val1],
614	[val2, val2-2, val2],
615	[val3-0.5, val3, val3]])
616
617
618
619	domain.order = 1
620	domain.distribute_to_vertices_and_edges()
621
622	#Check that centroid values were distributed to vertices
623	C = domain.quantities['level'].centroid_values
624	for i in range(3):
625	assert allclose( domain.quantities['level'].vertex_values[:,i], C)
626
627
628	def test_first_order_limiter_variable_z(self):
629	#Check that first order limiter follows bed_slope
630	from Numeric import alltrue, greater_equal
631	from config import epsilon
632
633	a = [0.0, 0.0]
634	b = [0.0, 2.0]
635	c = [2.0,0.0]
636	d = [0.0, 4.0]
637	e = [2.0, 2.0]
638	f = [4.0,0.0]
639
640	points = [a, b, c, d, e, f]
641	#bac, bce, ecf, dbe
642	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
643
644	domain = Domain(points, vertices)
645	val0 = 2.+2.0/3
646	val1 = 4.+4.0/3
647	val2 = 8.+2.0/3
648	val3 = 2.+8.0/3
649
650	domain.set_quantity('elevation', [[0,0,0], [6,0,0],
651	[6,6,6], [6,6,6]])
652	domain.set_quantity('level', [[val0, val0, val0],
653	[val1, val1, val1],
654	[val2, val2, val2],
655	[val3, val3, val3]])
656
657	E = domain.quantities['elevation'].vertex_values
658	L = domain.quantities['level'].vertex_values
659
660
661	#Check that some levels are not above elevation (within eps)
662	#- so that the limiter has something to work with
663	assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
664
665	domain.order = 1
666	domain.distribute_to_vertices_and_edges()
667
668	#Check that all levels are above elevation (within eps)
669	assert alltrue(alltrue(greater_equal(L,E-epsilon)))
670
671
672	#####################################################
673	def test_distribute_basic(self):
674	#Using test data generated by pyvolution-2
675	#Assuming no friction and flat bed (0.0)
676
677	a = [0.0, 0.0]
678	b = [0.0, 2.0]
679	c = [2.0, 0.0]
680	d = [0.0, 4.0]
681	e = [2.0, 2.0]
682	f = [4.0, 0.0]
683
684	points = [a, b, c, d, e, f]
685	#bac, bce, ecf, dbe
686	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
687
688	domain = Domain(points, vertices)
689
690	val0 = 2.
691	val1 = 4.
692	val2 = 8.
693	val3 = 2.
694
695	domain.set_quantity('level', [val0, val1, val2, val3], 'centroids')
696	L = domain.quantities['level'].vertex_values
697
698	#First order
699	domain.order = 1
700	domain.distribute_to_vertices_and_edges()
701	assert allclose(L[1], val1)
702
703	#Second order
704	domain.order = 2
705	domain.distribute_to_vertices_and_edges()
706	assert allclose(L[1], [2.2, 4.9, 4.9])
707
708
709
710	def test_distribute_away_from_bed(self):
711	#Using test data generated by pyvolution-2
712	#Assuming no friction and flat bed (0.0)
713
714	a = [0.0, 0.0]
715	b = [0.0, 2.0]
716	c = [2.0, 0.0]
717	d = [0.0, 4.0]
718	e = [2.0, 2.0]
719	f = [4.0, 0.0]
720
721	points = [a, b, c, d, e, f]
722	#bac, bce, ecf, dbe
723	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
724
725	domain = Domain(points, vertices)
726	L = domain.quantities['level'].vertex_values
727
728	def level(x,y):
729	return x**2
730
731	domain.set_quantity('level', level, 'centroids')
732
733	a, b = domain.quantities['level'].compute_gradients()
734	assert allclose(a[1], 3.33333334)
735	assert allclose(b[1], 0.0)
736
737	domain.order = 1
738	domain.distribute_to_vertices_and_edges()
739	assert allclose(L[1], 1.77777778)
740
741	domain.order = 2
742	domain.distribute_to_vertices_and_edges()
743	assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
744
745
746
747	def test_distribute_away_from_bed1(self):
748	#Using test data generated by pyvolution-2
749	#Assuming no friction and flat bed (0.0)
750
751	a = [0.0, 0.0]
752	b = [0.0, 2.0]
753	c = [2.0, 0.0]
754	d = [0.0, 4.0]
755	e = [2.0, 2.0]
756	f = [4.0, 0.0]
757
758	points = [a, b, c, d, e, f]
759	#bac, bce, ecf, dbe
760	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
761
762	domain = Domain(points, vertices)
763	L = domain.quantities['level'].vertex_values
764
765	def level(x,y):
766	return x4+y2
767
768	domain.set_quantity('level', level, 'centroids')
769	#print domain.quantities['level'].centroid_values
770
771	a, b = domain.quantities['level'].compute_gradients()
772	assert allclose(a[1], 25.18518519)
773	assert allclose(b[1], 3.33333333)
774
775	domain.order = 1
776	domain.distribute_to_vertices_and_edges()
777	assert allclose(L[1], 4.9382716)
778
779	domain.order = 2
780	domain.distribute_to_vertices_and_edges()
781	assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
782
783
784
785	def test_distribute_near_bed(self):
786	#Using test data generated by pyvolution-2
787	#Assuming no friction and flat bed (0.0)
788
789	a = [0.0, 0.0]
790	b = [0.0, 2.0]
791	c = [2.0, 0.0]
792	d = [0.0, 4.0]
793	e = [2.0, 2.0]
794	f = [4.0, 0.0]
795
796	points = [a, b, c, d, e, f]
797	#bac, bce, ecf, dbe
798	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
799
800	domain = Domain(points, vertices)
801
802
803	#Set up for a gradient of (3,0) at mid triangle
804	def slope(x, y):
805	return 10*x
806
807	h = 0.1
808	def level(x,y):
809	return slope(x,y)+h
810
811	domain.set_quantity('elevation', slope)
812	domain.set_quantity('level', level, 'centroids')
813
814	#print domain.quantities['elevation'].centroid_values
815	#print domain.quantities['level'].centroid_values
816
817	E = domain.quantities['elevation'].vertex_values
818	L = domain.quantities['level'].vertex_values
819
820	#print E
821	domain.order = 1
822	domain.distribute_to_vertices_and_edges()
823	##assert allclose(L[1], [0.19999999, 20.05, 20.05])
824	assert allclose(L[1], [0.1, 20.1, 20.1])
825
826	domain.order = 2
827	domain.distribute_to_vertices_and_edges()
828	assert allclose(L[1], [0.1, 20.1, 20.1])
829
830	def test_distribute_near_bed1(self):
831	#Using test data generated by pyvolution-2
832	#Assuming no friction and flat bed (0.0)
833
834	a = [0.0, 0.0]
835	b = [0.0, 2.0]
836	c = [2.0, 0.0]
837	d = [0.0, 4.0]
838	e = [2.0, 2.0]
839	f = [4.0, 0.0]
840
841	points = [a, b, c, d, e, f]
842	#bac, bce, ecf, dbe
843	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
844
845	domain = Domain(points, vertices)
846
847
848	#Set up for a gradient of (3,0) at mid triangle
849	def slope(x, y):
850	return x4+y2
851
852	h = 0.1
853	def level(x,y):
854	return slope(x,y)+h
855
856	domain.set_quantity('elevation', slope)
857	domain.set_quantity('level', level)
858
859	#print domain.quantities['elevation'].centroid_values
860	#print domain.quantities['level'].centroid_values
861
862	E = domain.quantities['elevation'].vertex_values
863	L = domain.quantities['level'].vertex_values
864
865	#print E
866	domain.order = 1
867	domain.distribute_to_vertices_and_edges()
868	##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
869	assert allclose(L[1], [4.1, 16.1, 20.1])
870
871	domain.order = 2
872	domain.distribute_to_vertices_and_edges()
873	assert allclose(L[1], [4.1, 16.1, 20.1])
874
875
876
877	def test_second_order_distribute_real_data(self):
878	#Using test data generated by pyvolution-2
879	#Assuming no friction and flat bed (0.0)
880
881	a = [0.0, 0.0]
882	b = [0.0, 1.0/5]
883	c = [0.0, 2.0/5]
884	d = [1.0/5, 0.0]
885	e = [1.0/5, 1.0/5]
886	f = [1.0/5, 2.0/5]
887	g = [2.0/5, 2.0/5]
888
889	points = [a, b, c, d, e, f, g]
890	#bae, efb, cbf, feg
891	vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
892
893	domain = Domain(points, vertices)
894
895	def slope(x, y):
896	return -x/3
897
898	domain.set_quantity('elevation', slope)
899	domain.set_quantity('level',
900	[0.01298164, 0.00365611, 0.01440365, -0.0381856437096],
901	'centroids')
902	domain.set_quantity('xmomentum',
903	[0.00670439, 0.01263789, 0.00647805, 0.0178180740668],
904	'centroids')
905	domain.set_quantity('ymomentum',
906	[-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866],
907	'centroids')
908
909	E = domain.quantities['elevation'].vertex_values
910	L = domain.quantities['level'].vertex_values
911	X = domain.quantities['xmomentum'].vertex_values
912	Y = domain.quantities['ymomentum'].vertex_values
913
914	#print E
915	domain.order = 2
916	domain.distribute_to_vertices_and_edges()
917
918	#print L[1,:]
919	#print X[1,:]
920	#print Y[1,:]
921
922	assert allclose(L[1,:], [-0.00825735775384, -0.00801881482869, 0.0272445025825])
923	assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
924	assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
925
926
927	def test_second_order_flat_bed_onestep(self):
928
929	from mesh_factory import rectangular
930	from shallow_water import Domain, Reflective_boundary,\
931	Dirichlet_boundary, Constant_height
932	from Numeric import array
933
934	#Create basic mesh
935	points, vertices, boundary = rectangular(6, 6)
936
937	#Create shallow water domain
938	domain = Domain(points, vertices, boundary)
939	domain.smooth = False
940	domain.default_order=2
941
942	# Boundary conditions
943	Br = Reflective_boundary(domain)
944	Bd = Dirichlet_boundary([0.1, 0., 0.])
945	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
946
947	domain.check_integrity()
948
949	#Evolution
950	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
951	pass# domain.write_time()
952
953	#Data from earlier version of pyvolution
954	assert allclose(domain.min_timestep, 0.0396825396825)
955	assert allclose(domain.max_timestep, 0.0396825396825)
956
957	assert allclose(domain.quantities['level'].centroid_values[:12],
958	[0.00171396, 0.02656103, 0.00241523, 0.02656103,
959	0.00241523, 0.02656103, 0.00241523, 0.02656103,
960	0.00241523, 0.02656103, 0.00241523, 0.0272623])
961
962	domain.distribute_to_vertices_and_edges()
963	assert allclose(domain.quantities['level'].vertex_values[:12,0],
964	[0.0001714, 0.02656103, 0.00024152,
965	0.02656103, 0.00024152, 0.02656103,
966	0.00024152, 0.02656103, 0.00024152,
967	0.02656103, 0.00024152, 0.0272623])
968
969	assert allclose(domain.quantities['level'].vertex_values[:12,1],
970	[0.00315012, 0.02656103, 0.00024152, 0.02656103,
971	0.00024152, 0.02656103, 0.00024152, 0.02656103,
972	0.00024152, 0.02656103, 0.00040506, 0.0272623])
973
974	assert allclose(domain.quantities['level'].vertex_values[:12,2],
975	[0.00182037, 0.02656103, 0.00676264,
976	0.02656103, 0.00676264, 0.02656103,
977	0.00676264, 0.02656103, 0.00676264,
978	0.02656103, 0.0065991, 0.0272623])
979
980	assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
981	[0.00113961, 0.01302432, 0.00148672,
982	0.01302432, 0.00148672, 0.01302432,
983	0.00148672, 0.01302432, 0.00148672 ,
984	0.01302432, 0.00148672, 0.01337143])
985
986	assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
987	[-2.91240050e-004 , 1.22721531e-004,
988	-1.22721531e-004, 1.22721531e-004 ,
989	-1.22721531e-004, 1.22721531e-004,
990	-1.22721531e-004 , 1.22721531e-004,
991	-1.22721531e-004, 1.22721531e-004,
992	-1.22721531e-004, -4.57969873e-005])
993
994
995
996	def test_second_order_flat_bed_moresteps(self):
997
998	from mesh_factory import rectangular
999	from shallow_water import Domain, Reflective_boundary,\
1000	Dirichlet_boundary, Constant_height
1001	from Numeric import array
1002
1003	#Create basic mesh
1004	points, vertices, boundary = rectangular(6, 6)
1005
1006	#Create shallow water domain
1007	domain = Domain(points, vertices, boundary)
1008	domain.smooth = False
1009	domain.default_order=2
1010
1011	# Boundary conditions
1012	Br = Reflective_boundary(domain)
1013	Bd = Dirichlet_boundary([0.1, 0., 0.])
1014	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1015
1016	domain.check_integrity()
1017
1018	#Evolution
1019	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
1020	pass
1021
1022	#Data from earlier version of pyvolution
1023	#assert allclose(domain.min_timestep, 0.0396825396825)
1024	#assert allclose(domain.max_timestep, 0.0396825396825)
1025	#print domain.quantities['level'].centroid_values
1026
1027
1028	def test_flatbed_first_order(self):
1029	from mesh_factory import rectangular
1030	from shallow_water import Domain,\
1031	Reflective_boundary, Dirichlet_boundary
1032
1033	from Numeric import array
1034
1035	#Create basic mesh
1036	N = 8
1037	points, vertices, boundary = rectangular(N, N)
1038
1039	#Create shallow water domain
1040	domain = Domain(points, vertices, boundary)
1041	domain.smooth = False
1042	domain.visualise = False
1043	domain.default_order=1
1044
1045	# Boundary conditions
1046	Br = Reflective_boundary(domain)
1047	Bd = Dirichlet_boundary([0.2,0.,0.])
1048
1049	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1050	domain.check_integrity()
1051
1052
1053	#Evolution
1054	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
1055	pass
1056	#domain.write_time()
1057
1058	assert allclose(domain.min_timestep, 0.0140413643926)
1059	assert allclose(domain.max_timestep, 0.0140947355753)
1060
1061	for i in range(3):
1062	assert allclose(domain.quantities['level'].edge_values[:4,i],
1063	[0.10730244,0.12337617,0.11200126,0.12605666])
1064
1065	assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
1066	[0.07610894,0.06901572,0.07284461,0.06819712])
1067
1068	assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
1069	[-0.0060238, -0.00157404, -0.00309633, -0.0001637])
1070
1071	def test_flatbed_second_order(self):
1072	from mesh_factory import rectangular
1073	from shallow_water import Domain,\
1074	Reflective_boundary, Dirichlet_boundary
1075
1076	from Numeric import array
1077
1078	#Create basic mesh
1079	N = 8
1080	points, vertices, boundary = rectangular(N, N)
1081
1082	#Create shallow water domain
1083	domain = Domain(points, vertices, boundary)
1084	domain.smooth = False
1085	domain.visualise = False
1086	domain.default_order=2
1087	#domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
1088
1089	# Boundary conditions
1090	Br = Reflective_boundary(domain)
1091	Bd = Dirichlet_boundary([0.2,0.,0.])
1092
1093	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1094	domain.check_integrity()
1095
1096	#Evolution
1097	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1098	pass
1099
1100
1101	assert allclose(domain.min_timestep, 0.0210448446782)
1102	assert allclose(domain.max_timestep, 0.0210448446782)
1103
1104	#print domain.quantities['level'].vertex_values[:4,0]
1105	#print domain.quantities['xmomentum'].vertex_values[:4,0]
1106	#print domain.quantities['ymomentum'].vertex_values[:4,0]
1107
1108	assert allclose(domain.quantities['level'].vertex_values[:4,0],
1109	[0.00101913,0.05352143,0.00104852,0.05354394])
1110
1111	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1112	[ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
1113
1114	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1115	# [0.00090581,0.03685719,0.00088303,0.03687313])
1116
1117	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1118	[-0.00139463,0.0006156,-0.00060364,0.00061827])
1119
1120
1121
1122	def test_flatbed_second_order_distribute(self):
1123	#Use real data from pyvolution 2
1124	#painfully setup and extracted.
1125	from mesh_factory import rectangular
1126	from shallow_water import Domain,\
1127	Reflective_boundary, Dirichlet_boundary
1128
1129	from Numeric import array
1130
1131	#Create basic mesh
1132	N = 8
1133	points, vertices, boundary = rectangular(N, N)
1134
1135	#Create shallow water domain
1136	domain = Domain(points, vertices, boundary)
1137	domain.smooth = False
1138	domain.visualise = False
1139	domain.default_order=domain.order=2
1140
1141	# Boundary conditions
1142	Br = Reflective_boundary(domain)
1143	Bd = Dirichlet_boundary([0.2,0.,0.])
1144
1145	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1146	domain.check_integrity()
1147
1148
1149
1150	for V in [False, True]:
1151	if V:
1152	#Set centroids as if system had been evolved
1153	L = zeros(2NN, Float)
1154	L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
1155	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1156	1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
1157	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1158	1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
1159	5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
1160	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1161	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1162	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1163	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1164	0.00000000e+000, 5.57305948e-005]
1165
1166	X = zeros(2NN, Float)
1167	X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
1168	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1169	8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
1170	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1171	8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
1172	3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
1173	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1174	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1175	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1176	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1177	0.00000000e+000, 4.57662812e-005]
1178
1179	Y = zeros(2NN, Float)
1180	Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
1181	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1182	-6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
1183	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1184	-6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
1185	-1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
1186	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1187	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1188	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1189	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1190	0.00000000e+000, -2.57635178e-005]
1191
1192
1193	domain.set_quantity('level', L, 'centroids')
1194	domain.set_quantity('xmomentum', X, 'centroids')
1195	domain.set_quantity('ymomentum', Y, 'centroids')
1196
1197	domain.check_integrity()
1198	else:
1199	#Evolution
1200	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1201	pass
1202	assert allclose(domain.min_timestep, 0.0210448446782)
1203	assert allclose(domain.max_timestep, 0.0210448446782)
1204
1205
1206	#Centroids were correct but not vertices.
1207	#Hence the check of distribute below.
1208	assert allclose(domain.quantities['level'].centroid_values[:4],
1209	[0.00721206,0.05352143,0.01009108,0.05354394])
1210
1211	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
1212	[0.00648352,0.03685719,0.00850733,0.03687313])
1213
1214	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
1215	[-0.00139463,0.0006156,-0.00060364,0.00061827])
1216
1217	#print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
1218	#print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
1219
1220	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
1221	##print domain.quantities['xmomentum'].centroid_values[17], V
1222	##print
1223	if not V:
1224	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
1225	else:
1226	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
1227
1228	import copy
1229	XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
1230	assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
1231
1232	domain.distribute_to_vertices_and_edges()
1233
1234	#assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
1235
1236	assert allclose(domain.quantities['xmomentum'].centroid_values[17],
1237	0.0)
1238
1239
1240	assert allclose(domain.quantities['level'].vertex_values[:4,0],
1241	[0.00101913,0.05352143,0.00104852,0.05354394])
1242
1243	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1244	[-0.00139463,0.0006156,-0.00060364,0.00061827])
1245
1246
1247	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1248	[0.00064835,0.03685719,0.00085073,0.03687313])
1249
1250
1251	#NB NO longer relvant:
1252
1253	#This was the culprit. First triangles vertex 0 had an
1254	#x-momentum of 0.0064835 instead of 0.00090581 and
1255	#third triangle had 0.00850733 instead of 0.00088303
1256	#print domain.quantities['xmomentum'].vertex_values[:4,0]
1257
1258	#print domain.quantities['xmomentum'].vertex_values[:4,0]
1259	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1260	# [0.00090581,0.03685719,0.00088303,0.03687313])
1261
1262
1263
1264
1265
1266	def test_bedslope_problem_first_order(self):
1267
1268	from mesh_factory import rectangular
1269	from shallow_water import Domain, Reflective_boundary, Constant_height
1270	from Numeric import array
1271
1272	#Create basic mesh
1273	points, vertices, boundary = rectangular(6, 6)
1274
1275	#Create shallow water domain
1276	domain = Domain(points, vertices, boundary)
1277	domain.smooth = False
1278	domain.default_order=1
1279
1280	#Bed-slope and friction
1281	def x_slope(x, y):
1282	return -x/3
1283
1284	domain.set_quantity('elevation', x_slope)
1285
1286	# Boundary conditions
1287	Br = Reflective_boundary(domain)
1288	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1289
1290	#Initial condition
1291	domain.set_quantity('level', Constant_height(x_slope, 0.05))
1292	domain.check_integrity()
1293
1294	#Evolution
1295	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1296	pass# domain.write_time()
1297
1298	assert allclose(domain.min_timestep, 0.050010003001)
1299	assert allclose(domain.max_timestep, 0.050010003001)
1300
1301
1302	def test_bedslope_problem_first_order_moresteps(self):
1303
1304	from mesh_factory import rectangular
1305	from shallow_water import Domain, Reflective_boundary, Constant_height
1306	from Numeric import array
1307
1308	#Create basic mesh
1309	points, vertices, boundary = rectangular(6, 6)
1310
1311	#Create shallow water domain
1312	domain = Domain(points, vertices, boundary)
1313	domain.smooth = False
1314	domain.default_order=1
1315
1316	#Bed-slope and friction
1317	def x_slope(x, y):
1318	return -x/3
1319
1320	domain.set_quantity('elevation', x_slope)
1321
1322	# Boundary conditions
1323	Br = Reflective_boundary(domain)
1324	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1325
1326	#Initial condition
1327	domain.set_quantity('level', Constant_height(x_slope, 0.05))
1328	domain.check_integrity()
1329
1330	#Evolution
1331	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
1332	pass# domain.write_time()
1333
1334	#Data from earlier version of pyvolution
1335	#print domain.quantities['level'].centroid_values
1336
1337	assert allclose(domain.quantities['level'].centroid_values,
1338	[-0.02998628, -0.01520652, -0.03043492,
1339	-0.0149132, -0.03004706, -0.01476251,
1340	-0.0298215, -0.01467976, -0.02988158,
1341	-0.01474662, -0.03036161, -0.01442995,
1342	-0.07624583, -0.06297061, -0.07733792,
1343	-0.06342237, -0.07695439, -0.06289595,
1344	-0.07635559, -0.0626065, -0.07633628,
1345	-0.06280072, -0.07739632, -0.06386738,
1346	-0.12161738, -0.11028239, -0.1223796,
1347	-0.11095953, -0.12189744, -0.11048616,
1348	-0.12074535, -0.10987605, -0.12014311,
1349	-0.10976691, -0.12096859, -0.11087692,
1350	-0.16868259, -0.15868061, -0.16801135,
1351	-0.1588003, -0.16674343, -0.15813323,
1352	-0.16457595, -0.15693826, -0.16281096,
1353	-0.15585154, -0.16283873, -0.15540068,
1354	-0.17450362, -0.19919913, -0.18062882,
1355	-0.19764131, -0.17783111, -0.19407213,
1356	-0.1736915, -0.19053624, -0.17228678,
1357	-0.19105634, -0.17920133, -0.1968828,
1358	-0.14244395, -0.14604641, -0.14473537,
1359	-0.1506107, -0.14510055, -0.14919522,
1360	-0.14175896, -0.14560798, -0.13911658,
1361	-0.14439383, -0.13924047, -0.14829043])
1362
1363
1364	def test_bedslope_problem_second_order_one_step(self):
1365
1366	from mesh_factory import rectangular
1367	from shallow_water import Domain, Reflective_boundary, Constant_height
1368	from Numeric import array
1369
1370	#Create basic mesh
1371	points, vertices, boundary = rectangular(6, 6)
1372
1373	#Create shallow water domain
1374	domain = Domain(points, vertices, boundary)
1375	domain.smooth = False
1376	domain.default_order=2
1377
1378	#Bed-slope and friction at vertices (and interpolated elsewhere)
1379	def x_slope(x, y):
1380	return -x/3
1381
1382	domain.set_quantity('elevation', x_slope)
1383
1384	# Boundary conditions
1385	Br = Reflective_boundary(domain)
1386	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1387
1388	#Initial condition
1389	domain.set_quantity('level', Constant_height(x_slope, 0.05))
1390	domain.check_integrity()
1391
1392	assert allclose(domain.quantities['level'].centroid_values,
1393	[0.01296296, 0.03148148, 0.01296296,
1394	0.03148148, 0.01296296, 0.03148148,
1395	0.01296296, 0.03148148, 0.01296296,
1396	0.03148148, 0.01296296, 0.03148148,
1397	-0.04259259, -0.02407407, -0.04259259,
1398	-0.02407407, -0.04259259, -0.02407407,
1399	-0.04259259, -0.02407407, -0.04259259,
1400	-0.02407407, -0.04259259, -0.02407407,
1401	-0.09814815, -0.07962963, -0.09814815,
1402	-0.07962963, -0.09814815, -0.07962963,
1403	-0.09814815, -0.07962963, -0.09814815,
1404	-0.07962963, -0.09814815, -0.07962963,
1405	-0.1537037 , -0.13518519, -0.1537037,
1406	-0.13518519, -0.1537037, -0.13518519,
1407	-0.1537037 , -0.13518519, -0.1537037,
1408	-0.13518519, -0.1537037, -0.13518519,
1409	-0.20925926, -0.19074074, -0.20925926,
1410	-0.19074074, -0.20925926, -0.19074074,
1411	-0.20925926, -0.19074074, -0.20925926,
1412	-0.19074074, -0.20925926, -0.19074074,
1413	-0.26481481, -0.2462963, -0.26481481,
1414	-0.2462963, -0.26481481, -0.2462963,
1415	-0.26481481, -0.2462963, -0.26481481,
1416	-0.2462963, -0.26481481, -0.2462963])
1417
1418
1419	#print domain.quantities['level'].extrapolate_second_order()
1420	#domain.distribute_to_vertices_and_edges()
1421	#print domain.quantities['level'].vertex_values[:,0]
1422
1423	#Evolution
1424	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1425	#domain.write_time()
1426	pass
1427
1428
1429	#print domain.quantities['level'].centroid_values
1430	assert allclose(domain.quantities['level'].centroid_values,
1431	[0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
1432	0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
1433	0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
1434	-0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
1435	-0.04186556, -0.0248011, -0.04186556, -0.02255593,
1436	-0.09966629, -0.08035666, -0.09742112, -0.08035666,
1437	-0.09742112, -0.08035666, -0.09742112, -0.08035666,
1438	-0.09742112, -0.08035666, -0.09742112, -0.07811149,
1439	-0.15522185, -0.13591222, -0.15297667, -0.13591222,
1440	-0.15297667, -0.13591222, -0.15297667, -0.13591222,
1441	-0.15297667, -0.13591222, -0.15297667, -0.13366704,
1442	-0.2107774, -0.19146777, -0.20853223, -0.19146777,
1443	-0.20853223, -0.19146777, -0.20853223, -0.19146777,
1444	-0.20853223, -0.19146777, -0.20853223, -0.1892226,
1445	-0.26120669, -0.24776246, -0.25865535, -0.24776246,
1446	-0.25865535, -0.24776246, -0.25865535, -0.24776246,
1447	-0.25865535, -0.24776246, -0.25865535, -0.24521113])
1448
1449
1450
1451	def test_bedslope_problem_second_order_two_steps(self):
1452
1453	from mesh_factory import rectangular
1454	from shallow_water import Domain, Reflective_boundary, Constant_height
1455	from Numeric import array
1456
1457	#Create basic mesh
1458	points, vertices, boundary = rectangular(6, 6)
1459
1460	#Create shallow water domain
1461	domain = Domain(points, vertices, boundary)
1462	domain.smooth = False
1463	domain.default_order=2
1464
1465	#Bed-slope and friction at vertices (and interpolated elsewhere)
1466	def x_slope(x, y):
1467	return -x/3
1468
1469	domain.set_quantity('elevation', x_slope)
1470
1471	# Boundary conditions
1472	Br = Reflective_boundary(domain)
1473	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1474
1475	#Initial condition
1476	domain.set_quantity('level', Constant_height(x_slope, 0.05))
1477	domain.check_integrity()
1478
1479	assert allclose(domain.quantities['level'].centroid_values,
1480	[0.01296296, 0.03148148, 0.01296296,
1481	0.03148148, 0.01296296, 0.03148148,
1482	0.01296296, 0.03148148, 0.01296296,
1483	0.03148148, 0.01296296, 0.03148148,
1484	-0.04259259, -0.02407407, -0.04259259,
1485	-0.02407407, -0.04259259, -0.02407407,
1486	-0.04259259, -0.02407407, -0.04259259,
1487	-0.02407407, -0.04259259, -0.02407407,
1488	-0.09814815, -0.07962963, -0.09814815,
1489	-0.07962963, -0.09814815, -0.07962963,
1490	-0.09814815, -0.07962963, -0.09814815,
1491	-0.07962963, -0.09814815, -0.07962963,
1492	-0.1537037 , -0.13518519, -0.1537037,
1493	-0.13518519, -0.1537037, -0.13518519,
1494	-0.1537037 , -0.13518519, -0.1537037,
1495	-0.13518519, -0.1537037, -0.13518519,
1496	-0.20925926, -0.19074074, -0.20925926,
1497	-0.19074074, -0.20925926, -0.19074074,
1498	-0.20925926, -0.19074074, -0.20925926,
1499	-0.19074074, -0.20925926, -0.19074074,
1500	-0.26481481, -0.2462963, -0.26481481,
1501	-0.2462963, -0.26481481, -0.2462963,
1502	-0.26481481, -0.2462963, -0.26481481,
1503	-0.2462963, -0.26481481, -0.2462963])
1504
1505
1506	#print domain.quantities['level'].extrapolate_second_order()
1507	#domain.distribute_to_vertices_and_edges()
1508	#print domain.quantities['level'].vertex_values[:,0]
1509
1510	#Evolution
1511	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
1512	pass
1513
1514
1515	#Data from earlier version of pyvolution ft=0.1
1516	assert allclose(domain.min_timestep, 0.0376895634803)
1517	assert allclose(domain.max_timestep, 0.0415635655309)
1518
1519
1520	assert allclose(domain.quantities['level'].centroid_values,
1521	[0.00855788, 0.01575204, 0.00994606, 0.01717072,
1522	0.01005985, 0.01716362, 0.01005985, 0.01716299,
1523	0.01007098, 0.01736248, 0.01216452, 0.02026776,
1524	-0.04462374, -0.02479045, -0.04199789, -0.0229465,
1525	-0.04184033, -0.02295693, -0.04184013, -0.02295675,
1526	-0.04184486, -0.0228168, -0.04028876, -0.02036486,
1527	-0.10029444, -0.08170809, -0.09772846, -0.08021704,
1528	-0.09760006, -0.08022143, -0.09759984, -0.08022124,
1529	-0.09760261, -0.08008893, -0.09603914, -0.07758209,
1530	-0.15584152, -0.13723138, -0.15327266, -0.13572906,
1531	-0.15314427, -0.13573349, -0.15314405, -0.13573331,
1532	-0.15314679, -0.13560104, -0.15158523, -0.13310701,
1533	-0.21208605, -0.19283913, -0.20955631, -0.19134189,
1534	-0.20942821, -0.19134598, -0.20942799, -0.1913458,
1535	-0.20943005, -0.19120952, -0.20781177, -0.18869401,
1536	-0.25384082, -0.2463294, -0.25047649, -0.24464654,
1537	-0.25031159, -0.24464253, -0.25031112, -0.24464253,
1538	-0.25031463, -0.24454764, -0.24885323, -0.24286438])
1539
1540
1541
1542
1543	def test_bedslope_problem_second_order_two_yieldsteps(self):
1544
1545	from mesh_factory import rectangular
1546	from shallow_water import Domain, Reflective_boundary, Constant_height
1547	from Numeric import array
1548
1549	#Create basic mesh
1550	points, vertices, boundary = rectangular(6, 6)
1551
1552	#Create shallow water domain
1553	domain = Domain(points, vertices, boundary)
1554	domain.smooth = False
1555	domain.default_order=2
1556
1557	#Bed-slope and friction at vertices (and interpolated elsewhere)
1558	def x_slope(x, y):
1559	return -x/3
1560
1561	domain.set_quantity('elevation', x_slope)
1562
1563	# Boundary conditions
1564	Br = Reflective_boundary(domain)
1565	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1566
1567	#Initial condition
1568	domain.set_quantity('level', Constant_height(x_slope, 0.05))
1569	domain.check_integrity()
1570
1571	assert allclose(domain.quantities['level'].centroid_values,
1572	[0.01296296, 0.03148148, 0.01296296,
1573	0.03148148, 0.01296296, 0.03148148,
1574	0.01296296, 0.03148148, 0.01296296,
1575	0.03148148, 0.01296296, 0.03148148,
1576	-0.04259259, -0.02407407, -0.04259259,
1577	-0.02407407, -0.04259259, -0.02407407,
1578	-0.04259259, -0.02407407, -0.04259259,
1579	-0.02407407, -0.04259259, -0.02407407,
1580	-0.09814815, -0.07962963, -0.09814815,
1581	-0.07962963, -0.09814815, -0.07962963,
1582	-0.09814815, -0.07962963, -0.09814815,
1583	-0.07962963, -0.09814815, -0.07962963,
1584	-0.1537037 , -0.13518519, -0.1537037,
1585	-0.13518519, -0.1537037, -0.13518519,
1586	-0.1537037 , -0.13518519, -0.1537037,
1587	-0.13518519, -0.1537037, -0.13518519,
1588	-0.20925926, -0.19074074, -0.20925926,
1589	-0.19074074, -0.20925926, -0.19074074,
1590	-0.20925926, -0.19074074, -0.20925926,
1591	-0.19074074, -0.20925926, -0.19074074,
1592	-0.26481481, -0.2462963, -0.26481481,
1593	-0.2462963, -0.26481481, -0.2462963,
1594	-0.26481481, -0.2462963, -0.26481481,
1595	-0.2462963, -0.26481481, -0.2462963])
1596
1597
1598	#print domain.quantities['level'].extrapolate_second_order()
1599	#domain.distribute_to_vertices_and_edges()
1600	#print domain.quantities['level'].vertex_values[:,0]
1601
1602	#Evolution
1603	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): #0.05??
1604	#domain.write_time()
1605	pass
1606
1607
1608
1609	assert allclose(domain.quantities['level'].centroid_values,
1610	[0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
1611	0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
1612	0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
1613	-0.0229465, -0.04184033, -0.02295693, -0.04184013,
1614	-0.02295675, -0.04184486, -0.0228168, -0.04028876,
1615	-0.02036486, -0.10029444, -0.08170809, -0.09772846,
1616	-0.08021704, -0.09760006, -0.08022143, -0.09759984,
1617	-0.08022124, -0.09760261, -0.08008893, -0.09603914,
1618	-0.07758209, -0.15584152, -0.13723138, -0.15327266,
1619	-0.13572906, -0.15314427, -0.13573349, -0.15314405,
1620	-0.13573331, -0.15314679, -0.13560104, -0.15158523,
1621	-0.13310701, -0.21208605, -0.19283913, -0.20955631,
1622	-0.19134189, -0.20942821, -0.19134598, -0.20942799,
1623	-0.1913458, -0.20943005, -0.19120952, -0.20781177,
1624	-0.18869401, -0.25384082, -0.2463294, -0.25047649,
1625	-0.24464654, -0.25031159, -0.24464253, -0.25031112,
1626	-0.24464253, -0.25031463, -0.24454764, -0.24885323,
1627	-0.24286438])
1628
1629
1630
1631	def test_bedslope_problem_second_order_more_steps(self):
1632
1633	from mesh_factory import rectangular
1634	from shallow_water import Domain, Reflective_boundary, Constant_height
1635	from Numeric import array
1636
1637	#Create basic mesh
1638	points, vertices, boundary = rectangular(6, 6)
1639
1640	#Create shallow water domain
1641	domain = Domain(points, vertices, boundary)
1642	domain.smooth = False
1643	domain.default_order=2
1644
1645	#Bed-slope and friction at vertices (and interpolated elsewhere)
1646	def x_slope(x, y):
1647	return -x/3
1648
1649	domain.set_quantity('elevation', x_slope)
1650
1651	# Boundary conditions
1652	Br = Reflective_boundary(domain)
1653	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1654
1655	#Initial condition
1656	domain.set_quantity('level', Constant_height(x_slope, 0.05))
1657	domain.check_integrity()
1658
1659	assert allclose(domain.quantities['level'].centroid_values,
1660	[0.01296296, 0.03148148, 0.01296296,
1661	0.03148148, 0.01296296, 0.03148148,
1662	0.01296296, 0.03148148, 0.01296296,
1663	0.03148148, 0.01296296, 0.03148148,
1664	-0.04259259, -0.02407407, -0.04259259,
1665	-0.02407407, -0.04259259, -0.02407407,
1666	-0.04259259, -0.02407407, -0.04259259,
1667	-0.02407407, -0.04259259, -0.02407407,
1668	-0.09814815, -0.07962963, -0.09814815,
1669	-0.07962963, -0.09814815, -0.07962963,
1670	-0.09814815, -0.07962963, -0.09814815,
1671	-0.07962963, -0.09814815, -0.07962963,
1672	-0.1537037 , -0.13518519, -0.1537037,
1673	-0.13518519, -0.1537037, -0.13518519,
1674	-0.1537037 , -0.13518519, -0.1537037,
1675	-0.13518519, -0.1537037, -0.13518519,
1676	-0.20925926, -0.19074074, -0.20925926,
1677	-0.19074074, -0.20925926, -0.19074074,
1678	-0.20925926, -0.19074074, -0.20925926,
1679	-0.19074074, -0.20925926, -0.19074074,
1680	-0.26481481, -0.2462963, -0.26481481,
1681	-0.2462963, -0.26481481, -0.2462963,
1682	-0.26481481, -0.2462963, -0.26481481,
1683	-0.2462963, -0.26481481, -0.2462963])
1684
1685
1686	#print domain.quantities['level'].extrapolate_second_order()
1687	#domain.distribute_to_vertices_and_edges()
1688	#print domain.quantities['level'].vertex_values[:,0]
1689
1690	#Evolution
1691	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
1692	pass
1693
1694
1695	assert allclose(domain.quantities['level'].centroid_values,
1696	[-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
1697	-0.02901975, -0.0141361, -0.02898816, -0.01418135, -0.02961409, -0.01403487,
1698	-0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
1699	-0.07571145, -0.06235231, -0.0756817, -0.06245309, -0.07652292, -0.06289946,
1700	-0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934, -0.1107174,
1701	-0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
1702	-0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493, -0.15697919,
1703	-0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
1704	-0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
1705	-0.1833593, -0.19505871, -0.1818806, -0.19418042, -0.18586159, -0.19576946,
1706	-0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
1707	-0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
1708
1709	assert allclose(domain.quantities['xmomentum'].centroid_values,
1710	[ 0.00831121, 0.00317948, 0.00731797, 0.00334939, 0.00764717, 0.00348053,
1711	0.00788729, 0.00356522, 0.00780649, 0.00341919, 0.00693525, 0.00310375,
1712	0.02166196, 0.01421475, 0.02017737, 0.01316839, 0.02037015, 0.01368659,
1713	0.02106, 0.01399161, 0.02074514, 0.01354935, 0.01887407, 0.0123113,
1714	0.03775083, 0.02855197, 0.03689337, 0.02759782, 0.03732848, 0.02812072,
1715	0.03872545, 0.02913348, 0.03880939, 0.02803804, 0.03546499, 0.0260039,
1716	0.0632131, 0.04730634, 0.0576324, 0.04592336, 0.05790921, 0.04690514,
1717	0.05986467, 0.04871165, 0.06170068, 0.04811572, 0.05657041, 0.04416292,
1718	0.08489642, 0.07188097, 0.07835261, 0.06843406, 0.07986412, 0.0698247,
1719	0.08201071, 0.07216756, 0.08378418, 0.07273624, 0.080399, 0.06645841,
1720	0.01631548, 0.04691608, 0.0206632, 0.044409, 0.02115518, 0.04560305,
1721	0.02160608, 0.04663725, 0.02174734, 0.04795559, 0.02281427, 0.05667111])
1722
1723
1724	assert allclose(domain.quantities['ymomentum'].centroid_values,
1725	[ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
1726	-9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
1727	-1.99445058e-004, -3.28493467e-004, 6.68217349e-005, -8.42042805e-006,
1728	5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
1729	-8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
1730	-2.70160173e-004, -5.35485454e-004, 2.60780997e-004, 3.12145471e-005,
1731	5.16189608e-004, 1.07069062e-004, 9.29989252e-005, -3.71211119e-004,
1732	1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
1733	-4.74025708e-004, -6.94463009e-004, 6.15092843e-005, 2.22106820e-004,
1734	-6.29589294e-004, 2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
1735	-4.17914641e-004, 6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
1736	-1.67424308e-003, -7.39485066e-004, -1.59966988e-003, 5.68262838e-005,
1737	-1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
1738	-1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
1739	-2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
1740	-1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
1741	-6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
1742	-4.50964850e-003, -3.06319963e-003, 6.08950810e-004, -4.79537921e-004])
1743
1744
1745
1746
1747	def test_temp_play(self):
1748
1749	from mesh_factory import rectangular
1750	from shallow_water import Domain, Reflective_boundary, Constant_height
1751	from Numeric import array
1752
1753	#Create basic mesh
1754	points, vertices, boundary = rectangular(5, 5)
1755
1756	#Create shallow water domain
1757	domain = Domain(points, vertices, boundary)
1758	domain.smooth = False
1759	domain.default_order=2
1760
1761	#Bed-slope and friction at vertices (and interpolated elsewhere)
1762	def x_slope(x, y):
1763	return -x/3
1764
1765	domain.set_quantity('elevation', x_slope)
1766
1767	# Boundary conditions
1768	Br = Reflective_boundary(domain)
1769	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1770
1771	#Initial condition
1772	domain.set_quantity('level', Constant_height(x_slope, 0.05))
1773	domain.check_integrity()
1774
1775	#Evolution
1776	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
1777	pass
1778
1779	assert allclose(domain.quantities['level'].centroid_values[:4],
1780	[0.00206836, 0.01296714, 0.00363415, 0.01438924])
1781	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
1782	[0.01360154, 0.00671133, 0.01264578, 0.00648503])
1783	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
1784	[-1.19201077e-003, -7.23647546e-004,
1785	-6.39083123e-005, 6.29815168e-005])
1786
1787
1788	def test_complex_bed(self):
1789	#No friction is tested here
1790
1791	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
1792	Transmissive_boundary, Time_boundary,\
1793	Weir_simple as Weir, Constant_height
1794
1795	from mesh_factory import rectangular
1796	from Numeric import array
1797
1798	N = 12
1799	points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
1800	origin=(-0.07, 0))
1801
1802
1803	domain = Domain(points, vertices, boundary)
1804	domain.smooth = False
1805	domain.visualise = False
1806	domain.default_order=2
1807
1808
1809	inflow_stage = 0.1
1810	Z = Weir(inflow_stage)
1811	domain.set_quantity('elevation', Z)
1812
1813	Br = Reflective_boundary(domain)
1814	Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
1815	domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
1816
1817	domain.set_quantity('level', Constant_height(Z, 0.))
1818
1819	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
1820	pass
1821
1822
1823	#print domain.quantities['level'].centroid_values
1824
1825
1826	assert allclose(domain.quantities['level'].centroid_values,
1827	[3.95822638e-002, 5.61022588e-002, 4.66437868e-002, 5.73081011e-002,
1828	4.72394613e-002, 5.74684939e-002, 4.74309483e-002, 5.77458084e-002,
1829	4.80628177e-002, 5.85656225e-002, 4.90498542e-002, 6.02609831e-002,
1830	1.18470315e-002, 1.75136443e-002, 1.18035266e-002, 2.15565695e-002,
1831	1.31620268e-002, 2.14351640e-002, 1.32351076e-002, 2.15450687e-002,
1832	1.36414028e-002, 2.24274619e-002, 1.51689511e-002, 2.21789655e-002,
1833	-7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
1834	-8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
1835	-8.10292716e-003, -3.88584984e-004, -7.35226124e-003, 2.73985295e-004,
1836	1.86166683e-001, 8.74070369e-002, 1.86166712e-001, 8.74035875e-002,
1837	6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
1838	6.11666725e-002, 8.73846774e-002, 1.86166697e-001, 8.74171550e-002,
1839	-4.83333333e-002, 1.18333333e-001, -4.83333333e-002, 1.18333333e-001,
1840	-4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
1841	-1.73333333e-001, -6.66666667e-003, -4.83333333e-002, 1.18333333e-001,
1842	-2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1843	-2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1844	-2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1845	-4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1846	-4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1847	-4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1848	-5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1849	-5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1850	-5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1851	-6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1852	-6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1853	-6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1854	-5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1855	-5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1856	-5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1857	-4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1858	-4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1859	-4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1860	-5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1861	-5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1862	-5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
1863
1864
1865
1866
1867	#-------------------------------------------------------------
1868	if __name__ == "__main__":
1869	suite = unittest.makeSuite(TestCase,'test')
1870	runner = unittest.TextTestRunner()
1871	runner.run(suite)
1872

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

Download in other formats: