Context Navigation

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

Last change on this file since 947 was 907, checked in by ole, 20 years ago
Addressed problem with artificial momentum generated by discontinuous water depths in the presence of steep slopes. Now very shallow water is limited with a separate h-limiter controlled by beta_h (see config.py) for details.
File size: 95.9 KB

Line
1	#!/usr/bin/env python
2
3	import unittest, os
4	from math import sqrt, pi
5
6	from config import g, epsilon
7	from Numeric import allclose, array, zeros, ones, Float, take
8	from shallow_water import *
9
10
11
12	#Variable windfield implemented using functions
13	def speed(t,x,y):
14	"""Large speeds halfway between center and edges
15	Low speeds at center and edges
16	"""
17
18	from math import sqrt, exp, cos, pi
19
20	x = array(x)
21	y = array(y)
22
23	N = len(x)
24	s = 0*x #New array
25
26	for k in range(N):
27
28	r = sqrt(x[k]2 + y[k]2)
29
30	factor = exp( -(r-0.15)**2 )
31
32	s[k] = 4000 * factor * (cos(t2pi/150) + 2)
33
34	return s
35
36
37	def scalar_func(t,x,y):
38	"""Function that returns a scalar.
39	Used to test error message when Numeric array is expected
40	"""
41
42	return 17.7
43
44
45	def angle(t,x,y):
46	"""Rotating field
47	"""
48	from math import sqrt, atan, pi
49
50	x = array(x)
51	y = array(y)
52
53	N = len(x)
54	a = 0*x #New array
55
56	for k in range(N):
57	r = sqrt(x[k]2 + y[k]2)
58
59	angle = atan(y[k]/x[k])
60
61	if x[k] < 0:
62	angle+=pi #atan in ]-pi/2; pi/2[
63
64	#Take normal direction
65	angle -= pi/2
66
67	#Ensure positive radians
68	if angle < 0:
69	angle += 2*pi
70
71	a[k] = angle/pi*180
72
73	return a
74
75
76	class TestCase(unittest.TestCase):
77	def setUp(self):
78	pass
79
80	def tearDown(self):
81	pass
82
83	def test_rotate(self):
84	normal = array([0.0,-1.0])
85
86	q = array([1.0,2.0,3.0])
87
88	r = rotate(q, normal, direction = 1)
89	assert r[0] == 1
90	assert r[1] == -3
91	assert r[2] == 2
92
93	w = rotate(r, normal, direction = -1)
94	assert allclose(w, q)
95
96	#Check error check
97	try:
98	rotate(r, array([1,1,1]) )
99	except:
100	pass
101	else:
102	raise 'Should have raised an exception'
103
104	def test_flux_zero_case(self):
105	ql = zeros( 3, Float )
106	qr = zeros( 3, Float )
107	normal = zeros( 2, Float )
108	zl = zr = 0.
109	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
110
111	assert allclose(flux, [0,0,0])
112	assert max_speed == 0.
113
114	def test_flux_constants(self):
115	w = 2.0
116
117	normal = array([1.,0])
118	ql = array([w, 0, 0])
119	qr = array([w, 0, 0])
120	zl = zr = 0.
121	h = w - (zl+zr)/2
122
123	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
124
125	assert allclose(flux, [0., 0.5gh**2, 0.])
126	assert max_speed == sqrt(g*h)
127
128	#def test_flux_slope(self):
129	# #FIXME: TODO
130	# w = 2.0
131	#
132	# normal = array([1.,0])
133	# ql = array([w, 0, 0])
134	# qr = array([w, 0, 0])
135	# zl = zr = 0.
136	# h = w - (zl+zr)/2
137	#
138	# flux, max_speed = flux_function(normal, ql, qr, zl, zr)
139	#
140	# assert allclose(flux, [0., 0.5gh**2, 0.])
141	# assert max_speed == sqrt(g*h)
142
143
144	def test_flux1(self):
145	#Use data from previous version of pyvolution
146	normal = array([1.,0])
147	ql = array([-0.2, 2, 3])
148	qr = array([-0.2, 2, 3])
149	zl = zr = -0.5
150	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
151
152	assert allclose(flux, [2.,13.77433333, 20.])
153	assert allclose(max_speed, 8.38130948661)
154
155
156	def test_flux2(self):
157	#Use data from previous version of pyvolution
158	normal = array([0., -1.])
159	ql = array([-0.075, 2, 3])
160	qr = array([-0.075, 2, 3])
161	zl = zr = -0.375
162	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
163
164	assert allclose(flux, [-3.,-20.0, -30.441])
165	assert allclose(max_speed, 11.7146428199)
166
167	def test_flux3(self):
168	#Use data from previous version of pyvolution
169	normal = array([-sqrt(2)/2, sqrt(2)/2])
170	ql = array([-0.075, 2, 3])
171	qr = array([-0.075, 2, 3])
172	zl = zr = -0.375
173	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
174
175	assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
176	assert allclose(max_speed, 4.0716654239)
177
178	def test_flux4(self):
179	#Use data from previous version of pyvolution
180	normal = array([-sqrt(2)/2, sqrt(2)/2])
181	ql = array([-0.34319278, 0.10254161, 0.07273855])
182	qr = array([-0.30683287, 0.1071986, 0.05930515])
183	zl = zr = -0.375
184	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
185
186	assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
187	assert allclose(max_speed, 1.31414103233)
188
189	def test_sw_domain_simple(self):
190	a = [0.0, 0.0]
191	b = [0.0, 2.0]
192	c = [2.0,0.0]
193	d = [0.0, 4.0]
194	e = [2.0, 2.0]
195	f = [4.0,0.0]
196
197	points = [a, b, c, d, e, f]
198	#bac, bce, ecf, dbe, daf, dae
199	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
200
201	domain = Domain(points, vertices)
202	domain.check_integrity()
203
204	for name in ['stage', 'xmomentum', 'ymomentum',
205	'elevation', 'friction']:
206	assert domain.quantities.has_key(name)
207
208
209	assert domain.get_conserved_quantities(0, edge=1) == 0.
210
211
212	def test_boundary_conditions(self):
213
214	a = [0.0, 0.0]
215	b = [0.0, 2.0]
216	c = [2.0,0.0]
217	d = [0.0, 4.0]
218	e = [2.0, 2.0]
219	f = [4.0,0.0]
220
221	points = [a, b, c, d, e, f]
222	#bac, bce, ecf, dbe
223	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
224	boundary = { (0, 0): 'Third',
225	(0, 2): 'First',
226	(2, 0): 'Second',
227	(2, 1): 'Second',
228	(3, 1): 'Second',
229	(3, 2): 'Third'}
230
231
232	domain = Domain(points, vertices, boundary)
233	domain.check_integrity()
234
235
236	domain.set_quantity('stage', [[1,2,3], [5,5,5],
237	[0,0,9], [-6, 3, 3]])
238
239	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
240	[3,3,3], [4, 4, 4]])
241
242	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
243	[30,30,30], [40, 40, 40]])
244
245
246	D = Dirichlet_boundary([5,2,1])
247	T = Transmissive_boundary(domain)
248	R = Reflective_boundary(domain)
249	domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
250
251	domain.update_boundary()
252
253	#Stage
254	assert domain.quantities['stage'].boundary_values[0] == 2.5
255	assert domain.quantities['stage'].boundary_values[0] ==\
256	domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
257	assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
258	assert domain.quantities['stage'].boundary_values[2] ==\
259	domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
260	assert domain.quantities['stage'].boundary_values[3] ==\
261	domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
262	assert domain.quantities['stage'].boundary_values[4] ==\
263	domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
264	assert domain.quantities['stage'].boundary_values[5] ==\
265	domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
266
267	#Xmomentum
268	assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
269	assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
270	assert domain.quantities['xmomentum'].boundary_values[2] ==\
271	domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
272	assert domain.quantities['xmomentum'].boundary_values[3] ==\
273	domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
274	assert domain.quantities['xmomentum'].boundary_values[4] ==\
275	domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
276	assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective
277
278	#Ymomentum
279	assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
280	assert domain.quantities['ymomentum'].boundary_values[1] == 1. #Dirichlet
281	assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
282	assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
283	assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
284	assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
285
286
287	def test_boundary_conditionsII(self):
288
289	a = [0.0, 0.0]
290	b = [0.0, 2.0]
291	c = [2.0,0.0]
292	d = [0.0, 4.0]
293	e = [2.0, 2.0]
294	f = [4.0,0.0]
295
296	points = [a, b, c, d, e, f]
297	#bac, bce, ecf, dbe
298	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
299	boundary = { (0, 0): 'Third',
300	(0, 2): 'First',
301	(2, 0): 'Second',
302	(2, 1): 'Second',
303	(3, 1): 'Second',
304	(3, 2): 'Third',
305	(0, 1): 'Internal'}
306
307
308	domain = Domain(points, vertices, boundary)
309	domain.check_integrity()
310
311
312	domain.set_quantity('stage', [[1,2,3], [5,5,5],
313	[0,0,9], [-6, 3, 3]])
314
315	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
316	[3,3,3], [4, 4, 4]])
317
318	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
319	[30,30,30], [40, 40, 40]])
320
321
322	D = Dirichlet_boundary([5,2,1])
323	T = Transmissive_boundary(domain)
324	R = Reflective_boundary(domain)
325	domain.set_boundary( {'First': D, 'Second': T,
326	'Third': R, 'Internal': None})
327
328	domain.update_boundary()
329	domain.check_integrity()
330
331
332	def test_compute_fluxes0(self):
333	#Do a full triangle and check that fluxes cancel out for
334	#the constant stage case
335
336	a = [0.0, 0.0]
337	b = [0.0, 2.0]
338	c = [2.0,0.0]
339	d = [0.0, 4.0]
340	e = [2.0, 2.0]
341	f = [4.0,0.0]
342
343	points = [a, b, c, d, e, f]
344	#bac, bce, ecf, dbe
345	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
346
347	domain = Domain(points, vertices)
348	domain.set_quantity('stage', [[2,2,2], [2,2,2],
349	[2,2,2], [2,2,2]])
350	domain.check_integrity()
351
352	assert allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
353	assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
354
355	zl=zr=0. #Assume flat bed
356
357	#Flux across right edge of volume 1
358	normal = domain.get_normal(1,0)
359	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
360	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
361	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
362
363	#Check that flux seen from other triangles is inverse
364	tmp = qr; qr=ql; ql=tmp
365	normal = domain.get_normal(2,2)
366	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
367	assert allclose(flux + flux0, 0.)
368
369	#Flux across upper edge of volume 1
370	normal = domain.get_normal(1,1)
371	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
372	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
373	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
374
375	#Check that flux seen from other triangles is inverse
376	tmp = qr; qr=ql; ql=tmp
377	normal = domain.get_normal(3,0)
378	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
379	assert allclose(flux + flux1, 0.)
380
381	#Flux across lower left hypotenuse of volume 1
382	normal = domain.get_normal(1,2)
383	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
384	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
385	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
386
387	#Check that flux seen from other triangles is inverse
388	tmp = qr; qr=ql; ql=tmp
389	normal = domain.get_normal(0,1)
390	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
391	assert allclose(flux + flux2, 0.)
392
393
394	#Scale by edgelengths, add up anc check that total flux is zero
395	e0 = domain.edgelengths[1, 0]
396	e1 = domain.edgelengths[1, 1]
397	e2 = domain.edgelengths[1, 2]
398
399	assert allclose(e0flux0+e1flux1+e2*flux2, 0.)
400
401	#Now check that compute_flux yields zeros as well
402	domain.compute_fluxes()
403
404	for name in ['stage', 'xmomentum', 'ymomentum']:
405	#print name, domain.quantities[name].explicit_update
406	assert allclose(domain.quantities[name].explicit_update[1], 0)
407
408
409
410	def test_compute_fluxes1(self):
411	#Use values from previous version
412
413	a = [0.0, 0.0]
414	b = [0.0, 2.0]
415	c = [2.0,0.0]
416	d = [0.0, 4.0]
417	e = [2.0, 2.0]
418	f = [4.0,0.0]
419
420	points = [a, b, c, d, e, f]
421	#bac, bce, ecf, dbe
422	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
423
424	domain = Domain(points, vertices)
425	val0 = 2.+2.0/3
426	val1 = 4.+4.0/3
427	val2 = 8.+2.0/3
428	val3 = 2.+8.0/3
429
430	domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
431	[val2, val2, val2], [val3, val3, val3]])
432	domain.check_integrity()
433
434	zl=zr=0. #Assume flat bed
435
436	#Flux across right edge of volume 1
437	normal = domain.get_normal(1,0)
438	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
439	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
440	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
441	assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
442	assert allclose(max_speed, 9.21592824046)
443
444
445	#Flux across upper edge of volume 1
446	normal = domain.get_normal(1,1)
447	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
448	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
449	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
450	assert allclose(flux1, [2.4098563, 0., 123.04444444])
451	assert allclose(max_speed, 7.22956891292)
452
453	#Flux across lower left hypotenuse of volume 1
454	normal = domain.get_normal(1,2)
455	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
456	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
457	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
458
459	assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
460	assert allclose(max_speed, 7.22956891292)
461
462	#Scale, add up and check that compute_fluxes is correct for vol 1
463	e0 = domain.edgelengths[1, 0]
464	e1 = domain.edgelengths[1, 1]
465	e2 = domain.edgelengths[1, 2]
466
467	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
468	assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
469
470
471	domain.compute_fluxes()
472
473	#assert allclose(total_flux, domain.explicit_update[1,:])
474	for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
475	assert allclose(total_flux[i],
476	domain.quantities[name].explicit_update[1])
477
478	#assert allclose(domain.explicit_update, [
479	# [0., -69.68888889, -69.68888889],
480	# [-0.68218178, -166.6, -35.93333333],
481	# [-111.77316251, 69.68888889, 0.],
482	# [-35.68522449, 0., 69.68888889]])
483
484	assert allclose(domain.quantities['stage'].explicit_update,
485	[0., -0.68218178, -111.77316251, -35.68522449])
486	assert allclose(domain.quantities['xmomentum'].explicit_update,
487	[-69.68888889, -166.6, 69.68888889, 0])
488	assert allclose(domain.quantities['ymomentum'].explicit_update,
489	[-69.68888889, -35.93333333, 0., 69.68888889])
490
491
492	#assert allclose(domain.quantities[name].explicit_update
493
494
495
496
497
498	def test_compute_fluxes2(self):
499	#Random values, incl momentum
500
501	a = [0.0, 0.0]
502	b = [0.0, 2.0]
503	c = [2.0,0.0]
504	d = [0.0, 4.0]
505	e = [2.0, 2.0]
506	f = [4.0,0.0]
507
508	points = [a, b, c, d, e, f]
509	#bac, bce, ecf, dbe
510	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
511
512	domain = Domain(points, vertices)
513	val0 = 2.+2.0/3
514	val1 = 4.+4.0/3
515	val2 = 8.+2.0/3
516	val3 = 2.+8.0/3
517
518	zl=zr=0 #Assume flat zero bed
519
520	domain.set_quantity('elevation', zl*ones( (4,3) ))
521
522
523	domain.set_quantity('stage', [[val0, val0-1, val0-2],
524	[val1, val1+1, val1],
525	[val2, val2-2, val2],
526	[val3-0.5, val3, val3]])
527
528	domain.set_quantity('xmomentum',
529	[[1, 2, 3], [3, 4, 5],
530	[1, -1, 0], [0, -2, 2]])
531
532	domain.set_quantity('ymomentum',
533	[[1, -1, 0], [0, -3, 2],
534	[0, 1, 0], [-1, 2, 2]])
535
536
537	domain.check_integrity()
538
539
540
541	#Flux across right edge of volume 1
542	normal = domain.get_normal(1,0)
543	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
544	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
545	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
546
547	#Flux across upper edge of volume 1
548	normal = domain.get_normal(1,1)
549	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
550	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
551	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
552
553	#Flux across lower left hypotenuse of volume 1
554	normal = domain.get_normal(1,2)
555	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
556	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
557	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
558
559	#Scale, add up and check that compute_fluxes is correct for vol 1
560	e0 = domain.edgelengths[1, 0]
561	e1 = domain.edgelengths[1, 1]
562	e2 = domain.edgelengths[1, 2]
563
564	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
565
566
567	domain.compute_fluxes()
568	for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
569	assert allclose(total_flux[i],
570	domain.quantities[name].explicit_update[1])
571	#assert allclose(total_flux, domain.explicit_update[1,:])
572
573
574	def test_compute_fluxes3(self):
575	#Random values, incl momentum
576
577	a = [0.0, 0.0]
578	b = [0.0, 2.0]
579	c = [2.0,0.0]
580	d = [0.0, 4.0]
581	e = [2.0, 2.0]
582	f = [4.0,0.0]
583
584	points = [a, b, c, d, e, f]
585	#bac, bce, ecf, dbe
586	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
587
588	domain = Domain(points, vertices)
589	val0 = 2.+2.0/3
590	val1 = 4.+4.0/3
591	val2 = 8.+2.0/3
592	val3 = 2.+8.0/3
593
594	zl=zr=-3.75 #Assume constant bed (must be less than stage)
595	domain.set_quantity('elevation', zl*ones( (4,3) ))
596
597
598	domain.set_quantity('stage', [[val0, val0-1, val0-2],
599	[val1, val1+1, val1],
600	[val2, val2-2, val2],
601	[val3-0.5, val3, val3]])
602
603	domain.set_quantity('xmomentum',
604	[[1, 2, 3], [3, 4, 5],
605	[1, -1, 0], [0, -2, 2]])
606
607	domain.set_quantity('ymomentum',
608	[[1, -1, 0], [0, -3, 2],
609	[0, 1, 0], [-1, 2, 2]])
610
611
612	domain.check_integrity()
613
614
615
616	#Flux across right edge of volume 1
617	normal = domain.get_normal(1,0)
618	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
619	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
620	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
621
622	#Flux across upper edge of volume 1
623	normal = domain.get_normal(1,1)
624	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
625	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
626	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
627
628	#Flux across lower left hypotenuse of volume 1
629	normal = domain.get_normal(1,2)
630	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
631	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
632	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
633
634	#Scale, add up and check that compute_fluxes is correct for vol 1
635	e0 = domain.edgelengths[1, 0]
636	e1 = domain.edgelengths[1, 1]
637	e2 = domain.edgelengths[1, 2]
638
639	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
640
641	domain.compute_fluxes()
642	for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
643	assert allclose(total_flux[i],
644	domain.quantities[name].explicit_update[1])
645
646
647
648	def test_catching_negative_heights(self):
649
650	a = [0.0, 0.0]
651	b = [0.0, 2.0]
652	c = [2.0,0.0]
653	d = [0.0, 4.0]
654	e = [2.0, 2.0]
655	f = [4.0,0.0]
656
657	points = [a, b, c, d, e, f]
658	#bac, bce, ecf, dbe
659	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
660
661	domain = Domain(points, vertices)
662	val0 = 2.+2.0/3
663	val1 = 4.+4.0/3
664	val2 = 8.+2.0/3
665	val3 = 2.+8.0/3
666
667	zl=zr=4 #Too large
668	domain.set_quantity('elevation', zl*ones( (4,3) ))
669	domain.set_quantity('stage', [[val0, val0-1, val0-2],
670	[val1, val1+1, val1],
671	[val2, val2-2, val2],
672	[val3-0.5, val3, val3]])
673
674	#Should fail
675	try:
676	domain.check_integrity()
677	except:
678	pass
679
680
681
682
683	#####################################################
684	def test_initial_condition(self):
685	"""Test that initial condition is output at time == 0
686	"""
687
688	from config import g
689	import copy
690
691	a = [0.0, 0.0]
692	b = [0.0, 2.0]
693	c = [2.0, 0.0]
694	d = [0.0, 4.0]
695	e = [2.0, 2.0]
696	f = [4.0, 0.0]
697
698	points = [a, b, c, d, e, f]
699	#bac, bce, ecf, dbe
700	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
701
702	domain = Domain(points, vertices)
703
704	#Set up for a gradient of (3,0) at mid triangle
705	def slope(x, y):
706	return 3*x
707
708	h = 0.1
709	def stage(x,y):
710	return slope(x,y)+h
711
712	domain.set_quantity('elevation', slope)
713	domain.set_quantity('stage', stage)
714
715	initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
716
717	domain.set_boundary({'exterior': Reflective_boundary(domain)})
718
719	#Evolution
720	for t in domain.evolve(yieldstep = 1.0, finaltime = 2.0):
721	stage = domain.quantities['stage'].vertex_values
722
723	if t == 0.0:
724	assert allclose(stage, initial_stage)
725	else:
726	assert not allclose(stage, initial_stage)
727
728
729
730
731	#####################################################
732	def test_gravity(self):
733	#Assuming no friction
734
735	from config import g
736
737	a = [0.0, 0.0]
738	b = [0.0, 2.0]
739	c = [2.0, 0.0]
740	d = [0.0, 4.0]
741	e = [2.0, 2.0]
742	f = [4.0, 0.0]
743
744	points = [a, b, c, d, e, f]
745	#bac, bce, ecf, dbe
746	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
747
748	domain = Domain(points, vertices)
749
750	#Set up for a gradient of (3,0) at mid triangle
751	def slope(x, y):
752	return 3*x
753
754	h = 0.1
755	def stage(x,y):
756	return slope(x,y)+h
757
758	domain.set_quantity('elevation', slope)
759	domain.set_quantity('stage', stage)
760
761	for name in domain.conserved_quantities:
762	assert allclose(domain.quantities[name].explicit_update, 0)
763	assert allclose(domain.quantities[name].semi_implicit_update, 0)
764
765	domain.compute_forcing_terms()
766
767	assert allclose(domain.quantities['stage'].explicit_update, 0)
768	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
769	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
770
771
772	def test_manning_friction(self):
773	from config import g
774
775	a = [0.0, 0.0]
776	b = [0.0, 2.0]
777	c = [2.0, 0.0]
778	d = [0.0, 4.0]
779	e = [2.0, 2.0]
780	f = [4.0, 0.0]
781
782	points = [a, b, c, d, e, f]
783	#bac, bce, ecf, dbe
784	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
785
786	domain = Domain(points, vertices)
787
788	#Set up for a gradient of (3,0) at mid triangle
789	def slope(x, y):
790	return 3*x
791
792	h = 0.1
793	def stage(x,y):
794	return slope(x,y)+h
795
796	eta = 0.07
797	domain.set_quantity('elevation', slope)
798	domain.set_quantity('stage', stage)
799	domain.set_quantity('friction', eta)
800
801	for name in domain.conserved_quantities:
802	assert allclose(domain.quantities[name].explicit_update, 0)
803	assert allclose(domain.quantities[name].semi_implicit_update, 0)
804
805	domain.compute_forcing_terms()
806
807	assert allclose(domain.quantities['stage'].explicit_update, 0)
808	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
809	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
810
811	assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
812	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
813	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
814
815	#Create some momentum for friction to work with
816	domain.set_quantity('xmomentum', 1)
817	S = -g * eta2 / h(7.0/3)
818
819	domain.compute_forcing_terms()
820	assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
821	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
822	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
823
824	#A more complex example
825	domain.quantities['stage'].semi_implicit_update[:] = 0.0
826	domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
827	domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
828
829	domain.set_quantity('xmomentum', 3)
830	domain.set_quantity('ymomentum', 4)
831
832	S = -g * eta*2 5 / h**(7.0/3)
833
834
835	domain.compute_forcing_terms()
836
837	assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
838	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
839	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
840
841	def test_constant_wind_stress(self):
842	from config import rho_a, rho_w, eta_w
843	from math import pi, cos, sin, sqrt
844
845	a = [0.0, 0.0]
846	b = [0.0, 2.0]
847	c = [2.0, 0.0]
848	d = [0.0, 4.0]
849	e = [2.0, 2.0]
850	f = [4.0, 0.0]
851
852	points = [a, b, c, d, e, f]
853	#bac, bce, ecf, dbe
854	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
855
856
857	domain = Domain(points, vertices)
858
859	#Flat surface with 1m of water
860	domain.set_quantity('elevation', 0)
861	domain.set_quantity('stage', 1.0)
862	domain.set_quantity('friction', 0)
863
864	Br = Reflective_boundary(domain)
865	domain.set_boundary({'exterior': Br})
866
867	#Setup only one forcing term, constant wind stress
868	s = 100
869	phi = 135
870	domain.forcing_terms = []
871	domain.forcing_terms.append( Wind_stress(s, phi) )
872
873	domain.compute_forcing_terms()
874
875
876	const = eta_w*rho_a/rho_w
877
878	#Convert to radians
879	phi = phi*pi/180
880
881	#Compute velocity vector (u, v)
882	u = s*cos(phi)
883	v = s*sin(phi)
884
885	#Compute wind stress
886	S = const * sqrt(u2 + v2)
887
888	assert allclose(domain.quantities['stage'].explicit_update, 0)
889	assert allclose(domain.quantities['xmomentum'].explicit_update, S*u)
890	assert allclose(domain.quantities['ymomentum'].explicit_update, S*v)
891
892
893	def test_variable_wind_stress(self):
894	from config import rho_a, rho_w, eta_w
895	from math import pi, cos, sin, sqrt
896
897	a = [0.0, 0.0]
898	b = [0.0, 2.0]
899	c = [2.0, 0.0]
900	d = [0.0, 4.0]
901	e = [2.0, 2.0]
902	f = [4.0, 0.0]
903
904	points = [a, b, c, d, e, f]
905	#bac, bce, ecf, dbe
906	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
907
908	domain = Domain(points, vertices)
909
910	#Flat surface with 1m of water
911	domain.set_quantity('elevation', 0)
912	domain.set_quantity('stage', 1.0)
913	domain.set_quantity('friction', 0)
914
915	Br = Reflective_boundary(domain)
916	domain.set_boundary({'exterior': Br})
917
918
919	domain.time = 5.54 #Take a random time (not zero)
920
921	#Setup only one forcing term, constant wind stress
922	s = 100
923	phi = 135
924	domain.forcing_terms = []
925	domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
926
927	domain.compute_forcing_terms()
928
929	#Compute reference solution
930	const = eta_w*rho_a/rho_w
931
932	N = domain.number_of_elements
933
934	xc = domain.get_centroid_coordinates()
935	t = domain.time
936
937	x = xc[:,0]
938	y = xc[:,1]
939	s_vec = speed(t,x,y)
940	phi_vec = angle(t,x,y)
941
942
943	for k in range(N):
944	#Convert to radians
945	phi = phi_vec[k]*pi/180
946	s = s_vec[k]
947
948	#Compute velocity vector (u, v)
949	u = s*cos(phi)
950	v = s*sin(phi)
951
952	#Compute wind stress
953	S = const * sqrt(u2 + v2)
954
955	assert allclose(domain.quantities['stage'].explicit_update[k], 0)
956	assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
957	assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
958
959
960
961
962	def test_windfield_from_file(self):
963	from config import rho_a, rho_w, eta_w
964	from math import pi, cos, sin, sqrt
965	from config import time_format
966	from util import file_function
967	import time
968
969
970	a = [0.0, 0.0]
971	b = [0.0, 2.0]
972	c = [2.0, 0.0]
973	d = [0.0, 4.0]
974	e = [2.0, 2.0]
975	f = [4.0, 0.0]
976
977	points = [a, b, c, d, e, f]
978	#bac, bce, ecf, dbe
979	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
980
981	domain = Domain(points, vertices)
982
983	#Flat surface with 1m of water
984	domain.set_quantity('elevation', 0)
985	domain.set_quantity('stage', 1.0)
986	domain.set_quantity('friction', 0)
987
988	Br = Reflective_boundary(domain)
989	domain.set_boundary({'exterior': Br})
990
991
992	domain.time = 7 #Take a time that is represented in file (not zero)
993
994	#Write wind stress file (ensure that domain.time is covered)
995	#Take x=1 and y=0
996	filename = 'test_windstress_from_file.txt'
997	start = time.mktime(time.strptime('2000', '%Y'))
998	fid = open(filename, 'w')
999	dt = 1 #One second interval
1000	t = 0.0
1001	while t <= 10.0:
1002	t_string = time.strftime(time_format, time.gmtime(t+start))
1003
1004	fid.write('%s, %f %f\n' %(t_string,
1005	speed(t,[1],[0])[0],
1006	angle(t,[1],[0])[0]))
1007	t += dt
1008
1009	fid.close()
1010
1011	#Setup wind stress
1012	F = file_function(filename)
1013	W = Wind_stress(F)
1014	domain.forcing_terms = []
1015	domain.forcing_terms.append(W)
1016
1017	domain.compute_forcing_terms()
1018
1019	#Compute reference solution
1020	const = eta_w*rho_a/rho_w
1021
1022	N = domain.number_of_elements
1023
1024	t = domain.time
1025
1026	s = speed(t,[1],[0])[0]
1027	phi = angle(t,[1],[0])[0]
1028
1029	#Convert to radians
1030	phi = phi*pi/180
1031
1032
1033	#Compute velocity vector (u, v)
1034	u = s*cos(phi)
1035	v = s*sin(phi)
1036
1037	#Compute wind stress
1038	S = const * sqrt(u2 + v2)
1039
1040	for k in range(N):
1041	assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1042	assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1043	assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1044
1045	os.remove(filename)
1046
1047	def test_wind_stress_error_condition(self):
1048	"""Test that windstress reacts properly when forcing functions
1049	are wrong - e.g. returns a scalar
1050	"""
1051
1052	from config import rho_a, rho_w, eta_w
1053	from math import pi, cos, sin, sqrt
1054
1055	a = [0.0, 0.0]
1056	b = [0.0, 2.0]
1057	c = [2.0, 0.0]
1058	d = [0.0, 4.0]
1059	e = [2.0, 2.0]
1060	f = [4.0, 0.0]
1061
1062	points = [a, b, c, d, e, f]
1063	#bac, bce, ecf, dbe
1064	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1065
1066	domain = Domain(points, vertices)
1067
1068	#Flat surface with 1m of water
1069	domain.set_quantity('elevation', 0)
1070	domain.set_quantity('stage', 1.0)
1071	domain.set_quantity('friction', 0)
1072
1073	Br = Reflective_boundary(domain)
1074	domain.set_boundary({'exterior': Br})
1075
1076
1077	domain.time = 5.54 #Take a random time (not zero)
1078
1079	#Setup only one forcing term, bad func
1080	domain.forcing_terms = []
1081
1082	try:
1083	domain.forcing_terms.append(Wind_stress(s = scalar_func,
1084	phi = angle))
1085	except AssertionError:
1086	pass
1087	else:
1088	msg = 'Should have raised exception'
1089	raise msg
1090
1091
1092	try:
1093	domain.forcing_terms.append(Wind_stress(s = speed,
1094	phi = scalar_func))
1095	except AssertionError:
1096	pass
1097	else:
1098	msg = 'Should have raised exception'
1099	raise msg
1100
1101	try:
1102	domain.forcing_terms.append(Wind_stress(s = speed,
1103	phi = 'xx'))
1104	except:
1105	pass
1106	else:
1107	msg = 'Should have raised exception'
1108	raise msg
1109
1110
1111	#####################################################
1112	def test_first_order_extrapolator_const_z(self):
1113
1114	a = [0.0, 0.0]
1115	b = [0.0, 2.0]
1116	c = [2.0, 0.0]
1117	d = [0.0, 4.0]
1118	e = [2.0, 2.0]
1119	f = [4.0, 0.0]
1120
1121	points = [a, b, c, d, e, f]
1122	#bac, bce, ecf, dbe
1123	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1124
1125	domain = Domain(points, vertices)
1126	val0 = 2.+2.0/3
1127	val1 = 4.+4.0/3
1128	val2 = 8.+2.0/3
1129	val3 = 2.+8.0/3
1130
1131	zl=zr=-3.75 #Assume constant bed (must be less than stage)
1132	domain.set_quantity('elevation', zl*ones( (4,3) ))
1133	domain.set_quantity('stage', [[val0, val0-1, val0-2],
1134	[val1, val1+1, val1],
1135	[val2, val2-2, val2],
1136	[val3-0.5, val3, val3]])
1137
1138
1139
1140	domain.order = 1
1141	domain.distribute_to_vertices_and_edges()
1142
1143	#Check that centroid values were distributed to vertices
1144	C = domain.quantities['stage'].centroid_values
1145	for i in range(3):
1146	assert allclose( domain.quantities['stage'].vertex_values[:,i], C)
1147
1148
1149	def test_first_order_limiter_variable_z(self):
1150	#Check that first order limiter follows bed_slope
1151	from Numeric import alltrue, greater_equal
1152	from config import epsilon
1153
1154	a = [0.0, 0.0]
1155	b = [0.0, 2.0]
1156	c = [2.0,0.0]
1157	d = [0.0, 4.0]
1158	e = [2.0, 2.0]
1159	f = [4.0,0.0]
1160
1161	points = [a, b, c, d, e, f]
1162	#bac, bce, ecf, dbe
1163	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1164
1165	domain = Domain(points, vertices)
1166	val0 = 2.+2.0/3
1167	val1 = 4.+4.0/3
1168	val2 = 8.+2.0/3
1169	val3 = 2.+8.0/3
1170
1171	domain.set_quantity('elevation', [[0,0,0], [6,0,0],
1172	[6,6,6], [6,6,6]])
1173	domain.set_quantity('stage', [[val0, val0, val0],
1174	[val1, val1, val1],
1175	[val2, val2, val2],
1176	[val3, val3, val3]])
1177
1178	E = domain.quantities['elevation'].vertex_values
1179	L = domain.quantities['stage'].vertex_values
1180
1181
1182	#Check that some stages are not above elevation (within eps)
1183	#- so that the limiter has something to work with
1184	assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
1185
1186	domain.order = 1
1187	domain.distribute_to_vertices_and_edges()
1188
1189	#Check that all stages are above elevation (within eps)
1190	assert alltrue(alltrue(greater_equal(L,E-epsilon)))
1191
1192
1193	#####################################################
1194	def test_distribute_basic(self):
1195	#Using test data generated by pyvolution-2
1196	#Assuming no friction and flat bed (0.0)
1197
1198	a = [0.0, 0.0]
1199	b = [0.0, 2.0]
1200	c = [2.0, 0.0]
1201	d = [0.0, 4.0]
1202	e = [2.0, 2.0]
1203	f = [4.0, 0.0]
1204
1205	points = [a, b, c, d, e, f]
1206	#bac, bce, ecf, dbe
1207	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1208
1209	domain = Domain(points, vertices)
1210
1211	val0 = 2.
1212	val1 = 4.
1213	val2 = 8.
1214	val3 = 2.
1215
1216	domain.set_quantity('stage', [val0, val1, val2, val3], 'centroids')
1217	L = domain.quantities['stage'].vertex_values
1218
1219	#First order
1220	domain.order = 1
1221	domain.distribute_to_vertices_and_edges()
1222	assert allclose(L[1], val1)
1223
1224	#Second order
1225	domain.order = 2
1226	domain.distribute_to_vertices_and_edges()
1227	assert allclose(L[1], [2.2, 4.9, 4.9])
1228
1229
1230
1231	def test_distribute_away_from_bed(self):
1232	#Using test data generated by pyvolution-2
1233	#Assuming no friction and flat bed (0.0)
1234
1235	a = [0.0, 0.0]
1236	b = [0.0, 2.0]
1237	c = [2.0, 0.0]
1238	d = [0.0, 4.0]
1239	e = [2.0, 2.0]
1240	f = [4.0, 0.0]
1241
1242	points = [a, b, c, d, e, f]
1243	#bac, bce, ecf, dbe
1244	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1245
1246	domain = Domain(points, vertices)
1247	L = domain.quantities['stage'].vertex_values
1248
1249	def stage(x,y):
1250	return x**2
1251
1252	domain.set_quantity('stage', stage, 'centroids')
1253
1254	a, b = domain.quantities['stage'].compute_gradients()
1255	assert allclose(a[1], 3.33333334)
1256	assert allclose(b[1], 0.0)
1257
1258	domain.order = 1
1259	domain.distribute_to_vertices_and_edges()
1260	assert allclose(L[1], 1.77777778)
1261
1262	domain.order = 2
1263	domain.distribute_to_vertices_and_edges()
1264	assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
1265
1266
1267
1268	def test_distribute_away_from_bed1(self):
1269	#Using test data generated by pyvolution-2
1270	#Assuming no friction and flat bed (0.0)
1271
1272	a = [0.0, 0.0]
1273	b = [0.0, 2.0]
1274	c = [2.0, 0.0]
1275	d = [0.0, 4.0]
1276	e = [2.0, 2.0]
1277	f = [4.0, 0.0]
1278
1279	points = [a, b, c, d, e, f]
1280	#bac, bce, ecf, dbe
1281	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1282
1283	domain = Domain(points, vertices)
1284	L = domain.quantities['stage'].vertex_values
1285
1286	def stage(x,y):
1287	return x4+y2
1288
1289	domain.set_quantity('stage', stage, 'centroids')
1290	#print domain.quantities['stage'].centroid_values
1291
1292	a, b = domain.quantities['stage'].compute_gradients()
1293	assert allclose(a[1], 25.18518519)
1294	assert allclose(b[1], 3.33333333)
1295
1296	domain.order = 1
1297	domain.distribute_to_vertices_and_edges()
1298	assert allclose(L[1], 4.9382716)
1299
1300	domain.order = 2
1301	domain.distribute_to_vertices_and_edges()
1302	assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
1303
1304
1305
1306	def test_distribute_near_bed(self):
1307	#Using test data generated by pyvolution-2
1308	#Assuming no friction and flat bed (0.0)
1309
1310	a = [0.0, 0.0]
1311	b = [0.0, 2.0]
1312	c = [2.0, 0.0]
1313	d = [0.0, 4.0]
1314	e = [2.0, 2.0]
1315	f = [4.0, 0.0]
1316
1317	points = [a, b, c, d, e, f]
1318	#bac, bce, ecf, dbe
1319	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1320
1321	domain = Domain(points, vertices)
1322
1323
1324	#Set up for a gradient of (3,0) at mid triangle
1325	def slope(x, y):
1326	return 10*x
1327
1328	h = 0.1
1329	def stage(x,y):
1330	return slope(x,y)+h
1331
1332	domain.set_quantity('elevation', slope)
1333	domain.set_quantity('stage', stage, 'centroids')
1334
1335	#print domain.quantities['elevation'].centroid_values
1336	#print domain.quantities['stage'].centroid_values
1337
1338	E = domain.quantities['elevation'].vertex_values
1339	L = domain.quantities['stage'].vertex_values
1340
1341	#print E
1342	domain.order = 1
1343	domain.distribute_to_vertices_and_edges()
1344	##assert allclose(L[1], [0.19999999, 20.05, 20.05])
1345	assert allclose(L[1], [0.1, 20.1, 20.1])
1346
1347	domain.order = 2
1348	domain.distribute_to_vertices_and_edges()
1349	assert allclose(L[1], [0.1, 20.1, 20.1])
1350
1351	def test_distribute_near_bed1(self):
1352	#Using test data generated by pyvolution-2
1353	#Assuming no friction and flat bed (0.0)
1354
1355	a = [0.0, 0.0]
1356	b = [0.0, 2.0]
1357	c = [2.0, 0.0]
1358	d = [0.0, 4.0]
1359	e = [2.0, 2.0]
1360	f = [4.0, 0.0]
1361
1362	points = [a, b, c, d, e, f]
1363	#bac, bce, ecf, dbe
1364	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1365
1366	domain = Domain(points, vertices)
1367
1368
1369	#Set up for a gradient of (3,0) at mid triangle
1370	def slope(x, y):
1371	return x4+y2
1372
1373	h = 0.1
1374	def stage(x,y):
1375	return slope(x,y)+h
1376
1377	domain.set_quantity('elevation', slope)
1378	domain.set_quantity('stage', stage)
1379
1380	#print domain.quantities['elevation'].centroid_values
1381	#print domain.quantities['stage'].centroid_values
1382
1383	E = domain.quantities['elevation'].vertex_values
1384	L = domain.quantities['stage'].vertex_values
1385
1386	#print E
1387	domain.order = 1
1388	domain.distribute_to_vertices_and_edges()
1389	##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
1390	assert allclose(L[1], [4.1, 16.1, 20.1])
1391
1392	domain.order = 2
1393	domain.distribute_to_vertices_and_edges()
1394	assert allclose(L[1], [4.1, 16.1, 20.1])
1395
1396
1397
1398	def test_second_order_distribute_real_data(self):
1399	#Using test data generated by pyvolution-2
1400	#Assuming no friction and flat bed (0.0)
1401
1402	a = [0.0, 0.0]
1403	b = [0.0, 1.0/5]
1404	c = [0.0, 2.0/5]
1405	d = [1.0/5, 0.0]
1406	e = [1.0/5, 1.0/5]
1407	f = [1.0/5, 2.0/5]
1408	g = [2.0/5, 2.0/5]
1409
1410	points = [a, b, c, d, e, f, g]
1411	#bae, efb, cbf, feg
1412	vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
1413
1414	domain = Domain(points, vertices)
1415
1416	def slope(x, y):
1417	return -x/3
1418
1419	domain.set_quantity('elevation', slope)
1420	domain.set_quantity('stage',
1421	[0.01298164, 0.00365611, 0.01440365, -0.0381856437096],
1422	'centroids')
1423	domain.set_quantity('xmomentum',
1424	[0.00670439, 0.01263789, 0.00647805, 0.0178180740668],
1425	'centroids')
1426	domain.set_quantity('ymomentum',
1427	[-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866],
1428	'centroids')
1429
1430	E = domain.quantities['elevation'].vertex_values
1431	L = domain.quantities['stage'].vertex_values
1432	X = domain.quantities['xmomentum'].vertex_values
1433	Y = domain.quantities['ymomentum'].vertex_values
1434
1435	#print E
1436	domain.order = 2
1437	domain.beta_h = 0.0 #Use first order in h-limiter
1438	domain.distribute_to_vertices_and_edges()
1439
1440	#print L[1,:]
1441	#print X[1,:]
1442	#print Y[1,:]
1443
1444	assert allclose(L[1,:], [-0.00825735775384, -0.00801881482869, 0.0272445025825])
1445	assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
1446	assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
1447
1448
1449	def test_second_order_flat_bed_onestep(self):
1450
1451	from mesh_factory import rectangular
1452	from shallow_water import Domain, Reflective_boundary,\
1453	Dirichlet_boundary, Constant_height
1454	from Numeric import array
1455
1456	#Create basic mesh
1457	points, vertices, boundary = rectangular(6, 6)
1458
1459	#Create shallow water domain
1460	domain = Domain(points, vertices, boundary)
1461	domain.smooth = False
1462	domain.default_order=2
1463
1464	# Boundary conditions
1465	Br = Reflective_boundary(domain)
1466	Bd = Dirichlet_boundary([0.1, 0., 0.])
1467	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1468
1469	domain.check_integrity()
1470
1471	#Evolution
1472	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1473	pass# domain.write_time()
1474
1475	#Data from earlier version of pyvolution
1476	assert allclose(domain.min_timestep, 0.0396825396825)
1477	assert allclose(domain.max_timestep, 0.0396825396825)
1478
1479	assert allclose(domain.quantities['stage'].centroid_values[:12],
1480	[0.00171396, 0.02656103, 0.00241523, 0.02656103,
1481	0.00241523, 0.02656103, 0.00241523, 0.02656103,
1482	0.00241523, 0.02656103, 0.00241523, 0.0272623])
1483
1484	domain.distribute_to_vertices_and_edges()
1485	assert allclose(domain.quantities['stage'].vertex_values[:12,0],
1486	[0.0001714, 0.02656103, 0.00024152,
1487	0.02656103, 0.00024152, 0.02656103,
1488	0.00024152, 0.02656103, 0.00024152,
1489	0.02656103, 0.00024152, 0.0272623])
1490
1491	assert allclose(domain.quantities['stage'].vertex_values[:12,1],
1492	[0.00315012, 0.02656103, 0.00024152, 0.02656103,
1493	0.00024152, 0.02656103, 0.00024152, 0.02656103,
1494	0.00024152, 0.02656103, 0.00040506, 0.0272623])
1495
1496	assert allclose(domain.quantities['stage'].vertex_values[:12,2],
1497	[0.00182037, 0.02656103, 0.00676264,
1498	0.02656103, 0.00676264, 0.02656103,
1499	0.00676264, 0.02656103, 0.00676264,
1500	0.02656103, 0.0065991, 0.0272623])
1501
1502	assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
1503	[0.00113961, 0.01302432, 0.00148672,
1504	0.01302432, 0.00148672, 0.01302432,
1505	0.00148672, 0.01302432, 0.00148672 ,
1506	0.01302432, 0.00148672, 0.01337143])
1507
1508	assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
1509	[-2.91240050e-004 , 1.22721531e-004,
1510	-1.22721531e-004, 1.22721531e-004 ,
1511	-1.22721531e-004, 1.22721531e-004,
1512	-1.22721531e-004 , 1.22721531e-004,
1513	-1.22721531e-004, 1.22721531e-004,
1514	-1.22721531e-004, -4.57969873e-005])
1515
1516
1517
1518	def test_second_order_flat_bed_moresteps(self):
1519
1520	from mesh_factory import rectangular
1521	from shallow_water import Domain, Reflective_boundary,\
1522	Dirichlet_boundary, Constant_height
1523	from Numeric import array
1524
1525	#Create basic mesh
1526	points, vertices, boundary = rectangular(6, 6)
1527
1528	#Create shallow water domain
1529	domain = Domain(points, vertices, boundary)
1530	domain.smooth = False
1531	domain.default_order=2
1532
1533	# Boundary conditions
1534	Br = Reflective_boundary(domain)
1535	Bd = Dirichlet_boundary([0.1, 0., 0.])
1536	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1537
1538	domain.check_integrity()
1539
1540	#Evolution
1541	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
1542	pass
1543
1544	#Data from earlier version of pyvolution
1545	#assert allclose(domain.min_timestep, 0.0396825396825)
1546	#assert allclose(domain.max_timestep, 0.0396825396825)
1547	#print domain.quantities['stage'].centroid_values
1548
1549
1550	def test_flatbed_first_order(self):
1551	from mesh_factory import rectangular
1552	from shallow_water import Domain,\
1553	Reflective_boundary, Dirichlet_boundary
1554
1555	from Numeric import array
1556
1557	#Create basic mesh
1558	N = 8
1559	points, vertices, boundary = rectangular(N, N)
1560
1561	#Create shallow water domain
1562	domain = Domain(points, vertices, boundary)
1563	domain.smooth = False
1564	domain.visualise = False
1565	domain.default_order=1
1566
1567	# Boundary conditions
1568	Br = Reflective_boundary(domain)
1569	Bd = Dirichlet_boundary([0.2,0.,0.])
1570
1571	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1572	domain.check_integrity()
1573
1574
1575	#Evolution
1576	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
1577	pass
1578	#domain.write_time()
1579
1580	#FIXME: These numbers were from version before 25/10
1581	#assert allclose(domain.min_timestep, 0.0140413643926)
1582	#assert allclose(domain.max_timestep, 0.0140947355753)
1583
1584	for i in range(3):
1585	#assert allclose(domain.quantities['stage'].edge_values[:4,i],
1586	# [0.10730244,0.12337617,0.11200126,0.12605666])
1587
1588	assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
1589	[0.07610894,0.06901572,0.07284461,0.06819712])
1590
1591	#assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
1592	# [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
1593
1594	def test_flatbed_second_order(self):
1595	from mesh_factory import rectangular
1596	from shallow_water import Domain,\
1597	Reflective_boundary, Dirichlet_boundary
1598
1599	from Numeric import array
1600
1601	#Create basic mesh
1602	N = 8
1603	points, vertices, boundary = rectangular(N, N)
1604
1605	#Create shallow water domain
1606	domain = Domain(points, vertices, boundary)
1607	domain.smooth = False
1608	domain.visualise = False
1609	domain.default_order=2
1610	#domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
1611
1612	# Boundary conditions
1613	Br = Reflective_boundary(domain)
1614	Bd = Dirichlet_boundary([0.2,0.,0.])
1615
1616	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1617	domain.check_integrity()
1618
1619	#Evolution
1620	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1621	pass
1622
1623
1624	assert allclose(domain.min_timestep, 0.0210448446782)
1625	assert allclose(domain.max_timestep, 0.0210448446782)
1626
1627	#print domain.quantities['stage'].vertex_values[:4,0]
1628	#print domain.quantities['xmomentum'].vertex_values[:4,0]
1629	#print domain.quantities['ymomentum'].vertex_values[:4,0]
1630
1631	#FIXME: These numbers were from version before 25/10
1632	#assert allclose(domain.quantities['stage'].vertex_values[:4,0],
1633	# [0.00101913,0.05352143,0.00104852,0.05354394])
1634
1635	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1636	[ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
1637
1638	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1639	# [0.00090581,0.03685719,0.00088303,0.03687313])
1640
1641	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1642	[-0.00139463,0.0006156,-0.00060364,0.00061827])
1643
1644
1645
1646	def test_flatbed_second_order_distribute(self):
1647	#Use real data from pyvolution 2
1648	#painfully setup and extracted.
1649	from mesh_factory import rectangular
1650	from shallow_water import Domain,\
1651	Reflective_boundary, Dirichlet_boundary
1652
1653	from Numeric import array
1654
1655	#Create basic mesh
1656	N = 8
1657	points, vertices, boundary = rectangular(N, N)
1658
1659	#Create shallow water domain
1660	domain = Domain(points, vertices, boundary)
1661	domain.smooth = False
1662	domain.visualise = False
1663	domain.default_order=domain.order=2
1664
1665	# Boundary conditions
1666	Br = Reflective_boundary(domain)
1667	Bd = Dirichlet_boundary([0.2,0.,0.])
1668
1669	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1670	domain.check_integrity()
1671
1672
1673
1674	for V in [False, True]:
1675	if V:
1676	#Set centroids as if system had been evolved
1677	L = zeros(2NN, Float)
1678	L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
1679	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1680	1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
1681	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1682	1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
1683	5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
1684	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1685	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1686	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1687	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1688	0.00000000e+000, 5.57305948e-005]
1689
1690	X = zeros(2NN, Float)
1691	X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
1692	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1693	8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
1694	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1695	8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
1696	3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
1697	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1698	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1699	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1700	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1701	0.00000000e+000, 4.57662812e-005]
1702
1703	Y = zeros(2NN, Float)
1704	Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
1705	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1706	-6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
1707	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1708	-6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
1709	-1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
1710	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1711	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1712	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1713	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1714	0.00000000e+000, -2.57635178e-005]
1715
1716
1717	domain.set_quantity('stage', L, 'centroids')
1718	domain.set_quantity('xmomentum', X, 'centroids')
1719	domain.set_quantity('ymomentum', Y, 'centroids')
1720
1721	domain.check_integrity()
1722	else:
1723	#Evolution
1724	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
1725	pass
1726	assert allclose(domain.min_timestep, 0.0210448446782)
1727	assert allclose(domain.max_timestep, 0.0210448446782)
1728
1729
1730	#Centroids were correct but not vertices.
1731	#Hence the check of distribute below.
1732	assert allclose(domain.quantities['stage'].centroid_values[:4],
1733	[0.00721206,0.05352143,0.01009108,0.05354394])
1734
1735	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
1736	[0.00648352,0.03685719,0.00850733,0.03687313])
1737
1738	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
1739	[-0.00139463,0.0006156,-0.00060364,0.00061827])
1740
1741	#print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
1742	#print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
1743
1744	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
1745	##print domain.quantities['xmomentum'].centroid_values[17], V
1746	##print
1747	if not V:
1748	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
1749	else:
1750	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
1751
1752	import copy
1753	XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
1754	assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
1755
1756	domain.distribute_to_vertices_and_edges()
1757
1758	#assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
1759
1760	assert allclose(domain.quantities['xmomentum'].centroid_values[17],
1761	0.0)
1762
1763
1764	#FIXME: These numbers were from version before 25/10
1765	#assert allclose(domain.quantities['stage'].vertex_values[:4,0],
1766	# [0.00101913,0.05352143,0.00104852,0.05354394])
1767
1768	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1769	[-0.00139463,0.0006156,-0.00060364,0.00061827])
1770
1771
1772	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1773	[0.00064835,0.03685719,0.00085073,0.03687313])
1774
1775
1776	#NB NO longer relvant:
1777
1778	#This was the culprit. First triangles vertex 0 had an
1779	#x-momentum of 0.0064835 instead of 0.00090581 and
1780	#third triangle had 0.00850733 instead of 0.00088303
1781	#print domain.quantities['xmomentum'].vertex_values[:4,0]
1782
1783	#print domain.quantities['xmomentum'].vertex_values[:4,0]
1784	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1785	# [0.00090581,0.03685719,0.00088303,0.03687313])
1786
1787
1788
1789
1790
1791	def test_bedslope_problem_first_order(self):
1792
1793	from mesh_factory import rectangular
1794	from shallow_water import Domain, Reflective_boundary, Constant_height
1795	from Numeric import array
1796
1797	#Create basic mesh
1798	points, vertices, boundary = rectangular(6, 6)
1799
1800	#Create shallow water domain
1801	domain = Domain(points, vertices, boundary)
1802	domain.smooth = False
1803	domain.default_order=1
1804
1805	#Bed-slope and friction
1806	def x_slope(x, y):
1807	return -x/3
1808
1809	domain.set_quantity('elevation', x_slope)
1810
1811	# Boundary conditions
1812	Br = Reflective_boundary(domain)
1813	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1814
1815	#Initial condition
1816	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
1817	domain.check_integrity()
1818
1819	#Evolution
1820	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1821	pass# domain.write_time()
1822
1823	assert allclose(domain.min_timestep, 0.050010003001)
1824	assert allclose(domain.max_timestep, 0.050010003001)
1825
1826
1827	def test_bedslope_problem_first_order_moresteps(self):
1828
1829	from mesh_factory import rectangular
1830	from shallow_water import Domain, Reflective_boundary, Constant_height
1831	from Numeric import array
1832
1833	#Create basic mesh
1834	points, vertices, boundary = rectangular(6, 6)
1835
1836	#Create shallow water domain
1837	domain = Domain(points, vertices, boundary)
1838	domain.smooth = False
1839	domain.default_order=1
1840	domain.beta_h = 0.0 #Use first order in h-limiter
1841
1842	#Bed-slope and friction
1843	def x_slope(x, y):
1844	return -x/3
1845
1846	domain.set_quantity('elevation', x_slope)
1847
1848	# Boundary conditions
1849	Br = Reflective_boundary(domain)
1850	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1851
1852	#Initial condition
1853	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
1854	domain.check_integrity()
1855
1856	#Evolution
1857	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
1858	pass# domain.write_time()
1859
1860	#Data from earlier version of pyvolution
1861	#print domain.quantities['stage'].centroid_values
1862
1863	assert allclose(domain.quantities['stage'].centroid_values,
1864	[-0.02998628, -0.01520652, -0.03043492,
1865	-0.0149132, -0.03004706, -0.01476251,
1866	-0.0298215, -0.01467976, -0.02988158,
1867	-0.01474662, -0.03036161, -0.01442995,
1868	-0.07624583, -0.06297061, -0.07733792,
1869	-0.06342237, -0.07695439, -0.06289595,
1870	-0.07635559, -0.0626065, -0.07633628,
1871	-0.06280072, -0.07739632, -0.06386738,
1872	-0.12161738, -0.11028239, -0.1223796,
1873	-0.11095953, -0.12189744, -0.11048616,
1874	-0.12074535, -0.10987605, -0.12014311,
1875	-0.10976691, -0.12096859, -0.11087692,
1876	-0.16868259, -0.15868061, -0.16801135,
1877	-0.1588003, -0.16674343, -0.15813323,
1878	-0.16457595, -0.15693826, -0.16281096,
1879	-0.15585154, -0.16283873, -0.15540068,
1880	-0.17450362, -0.19919913, -0.18062882,
1881	-0.19764131, -0.17783111, -0.19407213,
1882	-0.1736915, -0.19053624, -0.17228678,
1883	-0.19105634, -0.17920133, -0.1968828,
1884	-0.14244395, -0.14604641, -0.14473537,
1885	-0.1506107, -0.14510055, -0.14919522,
1886	-0.14175896, -0.14560798, -0.13911658,
1887	-0.14439383, -0.13924047, -0.14829043])
1888
1889
1890	def test_bedslope_problem_second_order_one_step(self):
1891
1892	from mesh_factory import rectangular
1893	from shallow_water import Domain, Reflective_boundary, Constant_height
1894	from Numeric import array
1895
1896	#Create basic mesh
1897	points, vertices, boundary = rectangular(6, 6)
1898
1899	#Create shallow water domain
1900	domain = Domain(points, vertices, boundary)
1901	domain.smooth = False
1902	domain.default_order=2
1903
1904	#Bed-slope and friction at vertices (and interpolated elsewhere)
1905	def x_slope(x, y):
1906	return -x/3
1907
1908	domain.set_quantity('elevation', x_slope)
1909
1910	# Boundary conditions
1911	Br = Reflective_boundary(domain)
1912	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1913
1914	#Initial condition
1915	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
1916	domain.check_integrity()
1917
1918	assert allclose(domain.quantities['stage'].centroid_values,
1919	[0.01296296, 0.03148148, 0.01296296,
1920	0.03148148, 0.01296296, 0.03148148,
1921	0.01296296, 0.03148148, 0.01296296,
1922	0.03148148, 0.01296296, 0.03148148,
1923	-0.04259259, -0.02407407, -0.04259259,
1924	-0.02407407, -0.04259259, -0.02407407,
1925	-0.04259259, -0.02407407, -0.04259259,
1926	-0.02407407, -0.04259259, -0.02407407,
1927	-0.09814815, -0.07962963, -0.09814815,
1928	-0.07962963, -0.09814815, -0.07962963,
1929	-0.09814815, -0.07962963, -0.09814815,
1930	-0.07962963, -0.09814815, -0.07962963,
1931	-0.1537037 , -0.13518519, -0.1537037,
1932	-0.13518519, -0.1537037, -0.13518519,
1933	-0.1537037 , -0.13518519, -0.1537037,
1934	-0.13518519, -0.1537037, -0.13518519,
1935	-0.20925926, -0.19074074, -0.20925926,
1936	-0.19074074, -0.20925926, -0.19074074,
1937	-0.20925926, -0.19074074, -0.20925926,
1938	-0.19074074, -0.20925926, -0.19074074,
1939	-0.26481481, -0.2462963, -0.26481481,
1940	-0.2462963, -0.26481481, -0.2462963,
1941	-0.26481481, -0.2462963, -0.26481481,
1942	-0.2462963, -0.26481481, -0.2462963])
1943
1944
1945	#print domain.quantities['stage'].extrapolate_second_order()
1946	#domain.distribute_to_vertices_and_edges()
1947	#print domain.quantities['stage'].vertex_values[:,0]
1948
1949	#Evolution
1950	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
1951	#domain.write_time()
1952	pass
1953
1954
1955	#print domain.quantities['stage'].centroid_values
1956	assert allclose(domain.quantities['stage'].centroid_values,
1957	[0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
1958	0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
1959	0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
1960	-0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
1961	-0.04186556, -0.0248011, -0.04186556, -0.02255593,
1962	-0.09966629, -0.08035666, -0.09742112, -0.08035666,
1963	-0.09742112, -0.08035666, -0.09742112, -0.08035666,
1964	-0.09742112, -0.08035666, -0.09742112, -0.07811149,
1965	-0.15522185, -0.13591222, -0.15297667, -0.13591222,
1966	-0.15297667, -0.13591222, -0.15297667, -0.13591222,
1967	-0.15297667, -0.13591222, -0.15297667, -0.13366704,
1968	-0.2107774, -0.19146777, -0.20853223, -0.19146777,
1969	-0.20853223, -0.19146777, -0.20853223, -0.19146777,
1970	-0.20853223, -0.19146777, -0.20853223, -0.1892226,
1971	-0.26120669, -0.24776246, -0.25865535, -0.24776246,
1972	-0.25865535, -0.24776246, -0.25865535, -0.24776246,
1973	-0.25865535, -0.24776246, -0.25865535, -0.24521113])
1974
1975
1976
1977	def test_bedslope_problem_second_order_two_steps(self):
1978
1979	from mesh_factory import rectangular
1980	from shallow_water import Domain, Reflective_boundary, Constant_height
1981	from Numeric import array
1982
1983	#Create basic mesh
1984	points, vertices, boundary = rectangular(6, 6)
1985
1986	#Create shallow water domain
1987	domain = Domain(points, vertices, boundary)
1988	domain.smooth = False
1989	domain.default_order=2
1990	domain.beta_h = 0.0 #Use first order in h-limiter
1991
1992	#Bed-slope and friction at vertices (and interpolated elsewhere)
1993	def x_slope(x, y):
1994	return -x/3
1995
1996	domain.set_quantity('elevation', x_slope)
1997
1998	# Boundary conditions
1999	Br = Reflective_boundary(domain)
2000	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2001
2002	#Initial condition
2003	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2004	domain.check_integrity()
2005
2006	assert allclose(domain.quantities['stage'].centroid_values,
2007	[0.01296296, 0.03148148, 0.01296296,
2008	0.03148148, 0.01296296, 0.03148148,
2009	0.01296296, 0.03148148, 0.01296296,
2010	0.03148148, 0.01296296, 0.03148148,
2011	-0.04259259, -0.02407407, -0.04259259,
2012	-0.02407407, -0.04259259, -0.02407407,
2013	-0.04259259, -0.02407407, -0.04259259,
2014	-0.02407407, -0.04259259, -0.02407407,
2015	-0.09814815, -0.07962963, -0.09814815,
2016	-0.07962963, -0.09814815, -0.07962963,
2017	-0.09814815, -0.07962963, -0.09814815,
2018	-0.07962963, -0.09814815, -0.07962963,
2019	-0.1537037 , -0.13518519, -0.1537037,
2020	-0.13518519, -0.1537037, -0.13518519,
2021	-0.1537037 , -0.13518519, -0.1537037,
2022	-0.13518519, -0.1537037, -0.13518519,
2023	-0.20925926, -0.19074074, -0.20925926,
2024	-0.19074074, -0.20925926, -0.19074074,
2025	-0.20925926, -0.19074074, -0.20925926,
2026	-0.19074074, -0.20925926, -0.19074074,
2027	-0.26481481, -0.2462963, -0.26481481,
2028	-0.2462963, -0.26481481, -0.2462963,
2029	-0.26481481, -0.2462963, -0.26481481,
2030	-0.2462963, -0.26481481, -0.2462963])
2031
2032
2033	#print domain.quantities['stage'].extrapolate_second_order()
2034	#domain.distribute_to_vertices_and_edges()
2035	#print domain.quantities['stage'].vertex_values[:,0]
2036
2037	#Evolution
2038	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2039	pass
2040
2041
2042	#Data from earlier version of pyvolution ft=0.1
2043	assert allclose(domain.min_timestep, 0.0376895634803)
2044	assert allclose(domain.max_timestep, 0.0415635655309)
2045
2046
2047	assert allclose(domain.quantities['stage'].centroid_values,
2048	[0.00855788, 0.01575204, 0.00994606, 0.01717072,
2049	0.01005985, 0.01716362, 0.01005985, 0.01716299,
2050	0.01007098, 0.01736248, 0.01216452, 0.02026776,
2051	-0.04462374, -0.02479045, -0.04199789, -0.0229465,
2052	-0.04184033, -0.02295693, -0.04184013, -0.02295675,
2053	-0.04184486, -0.0228168, -0.04028876, -0.02036486,
2054	-0.10029444, -0.08170809, -0.09772846, -0.08021704,
2055	-0.09760006, -0.08022143, -0.09759984, -0.08022124,
2056	-0.09760261, -0.08008893, -0.09603914, -0.07758209,
2057	-0.15584152, -0.13723138, -0.15327266, -0.13572906,
2058	-0.15314427, -0.13573349, -0.15314405, -0.13573331,
2059	-0.15314679, -0.13560104, -0.15158523, -0.13310701,
2060	-0.21208605, -0.19283913, -0.20955631, -0.19134189,
2061	-0.20942821, -0.19134598, -0.20942799, -0.1913458,
2062	-0.20943005, -0.19120952, -0.20781177, -0.18869401,
2063	-0.25384082, -0.2463294, -0.25047649, -0.24464654,
2064	-0.25031159, -0.24464253, -0.25031112, -0.24464253,
2065	-0.25031463, -0.24454764, -0.24885323, -0.24286438])
2066
2067
2068
2069
2070	def test_bedslope_problem_second_order_two_yieldsteps(self):
2071
2072	from mesh_factory import rectangular
2073	from shallow_water import Domain, Reflective_boundary, Constant_height
2074	from Numeric import array
2075
2076	#Create basic mesh
2077	points, vertices, boundary = rectangular(6, 6)
2078
2079	#Create shallow water domain
2080	domain = Domain(points, vertices, boundary)
2081	domain.smooth = False
2082	domain.default_order=2
2083	domain.beta_h = 0.0 #Use first order in h-limiter
2084
2085	#Bed-slope and friction at vertices (and interpolated elsewhere)
2086	def x_slope(x, y):
2087	return -x/3
2088
2089	domain.set_quantity('elevation', x_slope)
2090
2091	# Boundary conditions
2092	Br = Reflective_boundary(domain)
2093	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2094
2095	#Initial condition
2096	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2097	domain.check_integrity()
2098
2099	assert allclose(domain.quantities['stage'].centroid_values,
2100	[0.01296296, 0.03148148, 0.01296296,
2101	0.03148148, 0.01296296, 0.03148148,
2102	0.01296296, 0.03148148, 0.01296296,
2103	0.03148148, 0.01296296, 0.03148148,
2104	-0.04259259, -0.02407407, -0.04259259,
2105	-0.02407407, -0.04259259, -0.02407407,
2106	-0.04259259, -0.02407407, -0.04259259,
2107	-0.02407407, -0.04259259, -0.02407407,
2108	-0.09814815, -0.07962963, -0.09814815,
2109	-0.07962963, -0.09814815, -0.07962963,
2110	-0.09814815, -0.07962963, -0.09814815,
2111	-0.07962963, -0.09814815, -0.07962963,
2112	-0.1537037 , -0.13518519, -0.1537037,
2113	-0.13518519, -0.1537037, -0.13518519,
2114	-0.1537037 , -0.13518519, -0.1537037,
2115	-0.13518519, -0.1537037, -0.13518519,
2116	-0.20925926, -0.19074074, -0.20925926,
2117	-0.19074074, -0.20925926, -0.19074074,
2118	-0.20925926, -0.19074074, -0.20925926,
2119	-0.19074074, -0.20925926, -0.19074074,
2120	-0.26481481, -0.2462963, -0.26481481,
2121	-0.2462963, -0.26481481, -0.2462963,
2122	-0.26481481, -0.2462963, -0.26481481,
2123	-0.2462963, -0.26481481, -0.2462963])
2124
2125
2126	#print domain.quantities['stage'].extrapolate_second_order()
2127	#domain.distribute_to_vertices_and_edges()
2128	#print domain.quantities['stage'].vertex_values[:,0]
2129
2130	#Evolution
2131	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): #0.05??
2132	#domain.write_time()
2133	pass
2134
2135
2136
2137	assert allclose(domain.quantities['stage'].centroid_values,
2138	[0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
2139	0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
2140	0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
2141	-0.0229465, -0.04184033, -0.02295693, -0.04184013,
2142	-0.02295675, -0.04184486, -0.0228168, -0.04028876,
2143	-0.02036486, -0.10029444, -0.08170809, -0.09772846,
2144	-0.08021704, -0.09760006, -0.08022143, -0.09759984,
2145	-0.08022124, -0.09760261, -0.08008893, -0.09603914,
2146	-0.07758209, -0.15584152, -0.13723138, -0.15327266,
2147	-0.13572906, -0.15314427, -0.13573349, -0.15314405,
2148	-0.13573331, -0.15314679, -0.13560104, -0.15158523,
2149	-0.13310701, -0.21208605, -0.19283913, -0.20955631,
2150	-0.19134189, -0.20942821, -0.19134598, -0.20942799,
2151	-0.1913458, -0.20943005, -0.19120952, -0.20781177,
2152	-0.18869401, -0.25384082, -0.2463294, -0.25047649,
2153	-0.24464654, -0.25031159, -0.24464253, -0.25031112,
2154	-0.24464253, -0.25031463, -0.24454764, -0.24885323,
2155	-0.24286438])
2156
2157
2158
2159	def test_bedslope_problem_second_order_more_steps(self):
2160
2161	from mesh_factory import rectangular
2162	from shallow_water import Domain, Reflective_boundary, Constant_height
2163	from Numeric import array
2164
2165	#Create basic mesh
2166	points, vertices, boundary = rectangular(6, 6)
2167
2168	#Create shallow water domain
2169	domain = Domain(points, vertices, boundary)
2170	domain.smooth = False
2171	domain.default_order=2
2172	domain.beta_h = 0.0 #Use first order in h-limiter
2173
2174	#Bed-slope and friction at vertices (and interpolated elsewhere)
2175	def x_slope(x, y):
2176	return -x/3
2177
2178	domain.set_quantity('elevation', x_slope)
2179
2180	# Boundary conditions
2181	Br = Reflective_boundary(domain)
2182	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2183
2184	#Initial condition
2185	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2186	domain.check_integrity()
2187
2188	assert allclose(domain.quantities['stage'].centroid_values,
2189	[0.01296296, 0.03148148, 0.01296296,
2190	0.03148148, 0.01296296, 0.03148148,
2191	0.01296296, 0.03148148, 0.01296296,
2192	0.03148148, 0.01296296, 0.03148148,
2193	-0.04259259, -0.02407407, -0.04259259,
2194	-0.02407407, -0.04259259, -0.02407407,
2195	-0.04259259, -0.02407407, -0.04259259,
2196	-0.02407407, -0.04259259, -0.02407407,
2197	-0.09814815, -0.07962963, -0.09814815,
2198	-0.07962963, -0.09814815, -0.07962963,
2199	-0.09814815, -0.07962963, -0.09814815,
2200	-0.07962963, -0.09814815, -0.07962963,
2201	-0.1537037 , -0.13518519, -0.1537037,
2202	-0.13518519, -0.1537037, -0.13518519,
2203	-0.1537037 , -0.13518519, -0.1537037,
2204	-0.13518519, -0.1537037, -0.13518519,
2205	-0.20925926, -0.19074074, -0.20925926,
2206	-0.19074074, -0.20925926, -0.19074074,
2207	-0.20925926, -0.19074074, -0.20925926,
2208	-0.19074074, -0.20925926, -0.19074074,
2209	-0.26481481, -0.2462963, -0.26481481,
2210	-0.2462963, -0.26481481, -0.2462963,
2211	-0.26481481, -0.2462963, -0.26481481,
2212	-0.2462963, -0.26481481, -0.2462963])
2213
2214
2215	#print domain.quantities['stage'].extrapolate_second_order()
2216	#domain.distribute_to_vertices_and_edges()
2217	#print domain.quantities['stage'].vertex_values[:,0]
2218
2219	#Evolution
2220	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2221	pass
2222
2223
2224	assert allclose(domain.quantities['stage'].centroid_values,
2225	[-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
2226	-0.02901975, -0.0141361, -0.02898816, -0.01418135, -0.02961409, -0.01403487,
2227	-0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
2228	-0.07571145, -0.06235231, -0.0756817, -0.06245309, -0.07652292, -0.06289946,
2229	-0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934, -0.1107174,
2230	-0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
2231	-0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493, -0.15697919,
2232	-0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
2233	-0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
2234	-0.1833593, -0.19505871, -0.1818806, -0.19418042, -0.18586159, -0.19576946,
2235	-0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
2236	-0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
2237
2238	assert allclose(domain.quantities['xmomentum'].centroid_values,
2239	[ 0.00831121, 0.00317948, 0.00731797, 0.00334939, 0.00764717, 0.00348053,
2240	0.00788729, 0.00356522, 0.00780649, 0.00341919, 0.00693525, 0.00310375,
2241	0.02166196, 0.01421475, 0.02017737, 0.01316839, 0.02037015, 0.01368659,
2242	0.02106, 0.01399161, 0.02074514, 0.01354935, 0.01887407, 0.0123113,
2243	0.03775083, 0.02855197, 0.03689337, 0.02759782, 0.03732848, 0.02812072,
2244	0.03872545, 0.02913348, 0.03880939, 0.02803804, 0.03546499, 0.0260039,
2245	0.0632131, 0.04730634, 0.0576324, 0.04592336, 0.05790921, 0.04690514,
2246	0.05986467, 0.04871165, 0.06170068, 0.04811572, 0.05657041, 0.04416292,
2247	0.08489642, 0.07188097, 0.07835261, 0.06843406, 0.07986412, 0.0698247,
2248	0.08201071, 0.07216756, 0.08378418, 0.07273624, 0.080399, 0.06645841,
2249	0.01631548, 0.04691608, 0.0206632, 0.044409, 0.02115518, 0.04560305,
2250	0.02160608, 0.04663725, 0.02174734, 0.04795559, 0.02281427, 0.05667111])
2251
2252
2253	assert allclose(domain.quantities['ymomentum'].centroid_values,
2254	[ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
2255	-9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
2256	-1.99445058e-004, -3.28493467e-004, 6.68217349e-005, -8.42042805e-006,
2257	5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
2258	-8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
2259	-2.70160173e-004, -5.35485454e-004, 2.60780997e-004, 3.12145471e-005,
2260	5.16189608e-004, 1.07069062e-004, 9.29989252e-005, -3.71211119e-004,
2261	1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
2262	-4.74025708e-004, -6.94463009e-004, 6.15092843e-005, 2.22106820e-004,
2263	-6.29589294e-004, 2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
2264	-4.17914641e-004, 6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
2265	-1.67424308e-003, -7.39485066e-004, -1.59966988e-003, 5.68262838e-005,
2266	-1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
2267	-1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
2268	-2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
2269	-1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
2270	-6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
2271	-4.50964850e-003, -3.06319963e-003, 6.08950810e-004, -4.79537921e-004])
2272
2273
2274
2275
2276	def test_temp_play(self):
2277
2278	from mesh_factory import rectangular
2279	from shallow_water import Domain, Reflective_boundary, Constant_height
2280	from Numeric import array
2281
2282	#Create basic mesh
2283	points, vertices, boundary = rectangular(5, 5)
2284
2285	#Create shallow water domain
2286	domain = Domain(points, vertices, boundary)
2287	domain.smooth = False
2288	domain.default_order=2
2289	domain.beta_h = 0.0 #Use first order in h-limiter
2290
2291	#Bed-slope and friction at vertices (and interpolated elsewhere)
2292	def x_slope(x, y):
2293	return -x/3
2294
2295	domain.set_quantity('elevation', x_slope)
2296
2297	# Boundary conditions
2298	Br = Reflective_boundary(domain)
2299	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2300
2301	#Initial condition
2302	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2303	domain.check_integrity()
2304
2305	#Evolution
2306	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2307	pass
2308
2309	assert allclose(domain.quantities['stage'].centroid_values[:4],
2310	[0.00206836, 0.01296714, 0.00363415, 0.01438924])
2311	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
2312	[0.01360154, 0.00671133, 0.01264578, 0.00648503])
2313	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2314	[-1.19201077e-003, -7.23647546e-004,
2315	-6.39083123e-005, 6.29815168e-005])
2316
2317
2318	def test_complex_bed(self):
2319	#No friction is tested here
2320
2321	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2322	Transmissive_boundary, Time_boundary,\
2323	Weir_simple as Weir, Constant_height
2324
2325	from mesh_factory import rectangular
2326	from Numeric import array
2327
2328	N = 12
2329	points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
2330	origin=(-0.07, 0))
2331
2332
2333	domain = Domain(points, vertices, boundary)
2334	domain.smooth = False
2335	domain.visualise = False
2336	domain.default_order=2
2337
2338
2339	inflow_stage = 0.1
2340	Z = Weir(inflow_stage)
2341	domain.set_quantity('elevation', Z)
2342
2343	Br = Reflective_boundary(domain)
2344	Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
2345	domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
2346
2347	domain.set_quantity('stage', Constant_height(Z, 0.))
2348
2349	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
2350	pass
2351
2352
2353	#print domain.quantities['stage'].centroid_values
2354
2355	#FIXME: These numbers were from version before 25/10
2356	#assert allclose(domain.quantities['stage'].centroid_values,
2357	# [3.95822638e-002, 5.61022588e-002, 4.66437868e-002, 5.73081011e-002,
2358	# 4.72394613e-002, 5.74684939e-002, 4.74309483e-002, 5.77458084e-002,
2359	# 4.80628177e-002, 5.85656225e-002, 4.90498542e-002, 6.02609831e-002,
2360	# 1.18470315e-002, 1.75136443e-002, 1.18035266e-002, 2.15565695e-002,
2361	# 1.31620268e-002, 2.14351640e-002, 1.32351076e-002, 2.15450687e-002,
2362	# 1.36414028e-002, 2.24274619e-002, 1.51689511e-002, 2.21789655e-002,
2363	# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
2364	# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
2365	# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003, 2.73985295e-004,
2366	# 1.86166683e-001, 8.74070369e-002, 1.86166712e-001, 8.74035875e-002,
2367	# 6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
2368	# 6.11666725e-002, 8.73846774e-002, 1.86166697e-001, 8.74171550e-002,
2369	# -4.83333333e-002, 1.18333333e-001, -4.83333333e-002, 1.18333333e-001,
2370	# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
2371	# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002, 1.18333333e-001,
2372	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2373	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2374	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2375	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2376	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2377	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2378	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2379	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2380	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2381	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2382	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2383	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2384	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2385	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2386	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2387	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2388	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2389	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2390	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2391	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2392	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
2393
2394
2395
2396	def test_spatio_temporal_boundary(self):
2397	"""Test that boundary values can be read from file and interpolated
2398	in both time and space.
2399
2400	Verify that the same steady state solution is arrived at and that
2401	time interpolation works.
2402
2403	The full solution history is not exactly the same as
2404	file boundary must read and interpolate from smoothed version
2405	as stored in sww.
2406	"""
2407	import time
2408
2409	#Create sww file of simple propagation from left to right
2410	#through rectangular domain
2411
2412	from mesh_factory import rectangular
2413
2414	#Create basic mesh
2415	points, vertices, boundary = rectangular(3, 3)
2416
2417	#Create shallow water domain
2418	domain1 = Domain(points, vertices, boundary)
2419
2420	from util import mean
2421	domain1.reduction = mean
2422	domain1.smooth = False #Exact result
2423
2424	domain1.default_order = 2
2425	domain1.store = True
2426	domain1.set_datadir('.')
2427	domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
2428
2429	#FIXME: This is extremely important!
2430	#How can we test if they weren't stored?
2431	domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
2432
2433
2434	#Bed-slope and friction at vertices (and interpolated elsewhere)
2435	domain1.set_quantity('elevation', 0)
2436	domain1.set_quantity('friction', 0)
2437
2438	# Boundary conditions
2439	Br = Reflective_boundary(domain1)
2440	Bd = Dirichlet_boundary([0.3,0,0])
2441	domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
2442	#Initial condition
2443	domain1.set_quantity('stage', 0)
2444	domain1.check_integrity()
2445
2446	finaltime = 5
2447	#Evolution (full domain - large steps)
2448	for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
2449	pass
2450	#domain1.write_time()
2451
2452	cv1 = domain1.quantities['stage'].centroid_values
2453
2454
2455	#Create an triangle shaped domain (reusing coordinates from domain 1),
2456	#formed from the lower and right hand boundaries and
2457	#the sw-ne diagonal
2458	#from domain 1. Call it domain2
2459
2460	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
2461	[2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
2462	[1,0], [1,1.0/3], [1,2.0/3], [1,1]]
2463
2464	vertices = [ [1,2,0],
2465	[3,4,1], [2,1,4], [4,5,2],
2466	[6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
2467
2468	boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
2469	(4,2):'right', (6,2):'right', (8,2):'right',
2470	(0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
2471
2472	domain2 = Domain(points, vertices, boundary)
2473
2474	domain2.reduction = domain1.reduction
2475	domain2.smooth = False
2476	domain2.default_order = 2
2477
2478	#Bed-slope and friction at vertices (and interpolated elsewhere)
2479	domain2.set_quantity('elevation', 0)
2480	domain2.set_quantity('friction', 0)
2481	domain2.set_quantity('stage', 0)
2482
2483	# Boundary conditions
2484	Br = Reflective_boundary(domain2)
2485	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
2486	domain2)
2487	domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
2488	domain2.check_integrity()
2489
2490
2491
2492	#Evolution (small steps)
2493	for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
2494	pass
2495
2496
2497	#Use output from domain1 as spatio-temporal boundary for domain2
2498	#and verify that results at right hand side are close.
2499
2500	cv2 = domain2.quantities['stage'].centroid_values
2501
2502	#print take(cv1, (12,14,16)) #Right
2503	#print take(cv2, (4,6,8))
2504	#print take(cv1, (0,6,12)) #Bottom
2505	#print take(cv2, (0,1,4))
2506	#print take(cv1, (0,8,16)) #Diag
2507	#print take(cv2, (0,3,8))
2508
2509	assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
2510	assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
2511	assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
2512
2513	#Cleanup
2514	os.remove(domain1.filename + '.' + domain1.format)
2515
2516
2517
2518	def test_spatio_temporal_boundary_2(self):
2519	"""Test that boundary values can be read from file and interpolated
2520	in both time and space.
2521	This is a more basic test, verifying that boundary object
2522	produces the expected results
2523
2524
2525	"""
2526	import time
2527
2528	#Create sww file of simple propagation from left to right
2529	#through rectangular domain
2530
2531	from mesh_factory import rectangular
2532
2533	#Create basic mesh
2534	points, vertices, boundary = rectangular(3, 3)
2535
2536	#Create shallow water domain
2537	domain1 = Domain(points, vertices, boundary)
2538
2539	from util import mean
2540	domain1.reduction = mean
2541	domain1.smooth = True #To mimic MOST output
2542
2543	domain1.default_order = 2
2544	domain1.store = True
2545	domain1.set_datadir('.')
2546	domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
2547
2548	#FIXME: This is extremely important!
2549	#How can we test if they weren't stored?
2550	domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
2551
2552
2553	#Bed-slope and friction at vertices (and interpolated elsewhere)
2554	domain1.set_quantity('elevation', 0)
2555	domain1.set_quantity('friction', 0)
2556
2557	# Boundary conditions
2558	Br = Reflective_boundary(domain1)
2559	Bd = Dirichlet_boundary([0.3,0,0])
2560	domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
2561	#Initial condition
2562	domain1.set_quantity('stage', 0)
2563	domain1.check_integrity()
2564
2565	finaltime = 5
2566	#Evolution (full domain - large steps)
2567	for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
2568	pass
2569	#domain1.write_time()
2570
2571
2572	#Create an triangle shaped domain (coinciding with some
2573	#coordinates from domain 1),
2574	#formed from the lower and right hand boundaries and
2575	#the sw-ne diagonal
2576	#from domain 1. Call it domain2
2577
2578	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
2579	[2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
2580	[1,0], [1,1.0/3], [1,2.0/3], [1,1]]
2581
2582	vertices = [ [1,2,0],
2583	[3,4,1], [2,1,4], [4,5,2],
2584	[6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
2585
2586	boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
2587	(4,2):'right', (6,2):'right', (8,2):'right',
2588	(0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
2589
2590	domain2 = Domain(points, vertices, boundary)
2591
2592	domain2.reduction = domain1.reduction
2593	domain2.smooth = False
2594	domain2.default_order = 2
2595
2596	#Bed-slope and friction at vertices (and interpolated elsewhere)
2597	domain2.set_quantity('elevation', 0)
2598	domain2.set_quantity('friction', 0)
2599	domain2.set_quantity('stage', 0)
2600
2601
2602	#Read results for specific timesteps t=1 and t=2
2603	from Scientific.IO.NetCDF import NetCDFFile
2604	fid = NetCDFFile(domain1.filename + '.' + domain1.format)
2605
2606	x = fid.variables['x'][:]
2607	y = fid.variables['y'][:]
2608	s1 = fid.variables['stage'][1,:]
2609	s2 = fid.variables['stage'][2,:]
2610	fid.close()
2611
2612	from Numeric import take, reshape, concatenate
2613	shp = (len(x), 1)
2614	points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
2615	#The diagonal points of domain 1 are 0, 5, 10, 15
2616
2617	#print points[0], points[5], points[10], points[15]
2618	assert allclose( take(points, [0,5,10,15]),
2619	[[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
2620
2621
2622	# Boundary conditions
2623	Br = Reflective_boundary(domain2)
2624	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
2625	domain2)
2626	domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
2627	domain2.check_integrity()
2628
2629	#Test that interpolation points are the mid points of the all boundary
2630	#segments
2631
2632	boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
2633	[1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
2634	[1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
2635
2636	boundary_midpoints.sort()
2637	R = Bf.F.interpolation_points.tolist()
2638	R.sort()
2639	assert allclose(boundary_midpoints, R)
2640
2641	#Check spatially interpolated output at time == 1
2642	domain2.time = 1
2643
2644	#First diagonal midpoint
2645	R0 = Bf.evaluate(0,0)
2646	assert allclose(R0[0], (s1[0] + s1[5])/2)
2647
2648	#Second diagonal midpoint
2649	R0 = Bf.evaluate(3,0)
2650	assert allclose(R0[0], (s1[5] + s1[10])/2)
2651
2652	#First diagonal midpoint
2653	R0 = Bf.evaluate(8,0)
2654	assert allclose(R0[0], (s1[10] + s1[15])/2)
2655
2656	#Check spatially interpolated output at time == 2
2657	domain2.time = 2
2658
2659	#First diagonal midpoint
2660	R0 = Bf.evaluate(0,0)
2661	assert allclose(R0[0], (s2[0] + s2[5])/2)
2662
2663	#Second diagonal midpoint
2664	R0 = Bf.evaluate(3,0)
2665	assert allclose(R0[0], (s2[5] + s2[10])/2)
2666
2667	#First diagonal midpoint
2668	R0 = Bf.evaluate(8,0)
2669	assert allclose(R0[0], (s2[10] + s2[15])/2)
2670
2671
2672	#Now check temporal interpolation
2673
2674	domain2.time = 1 + 2.0/3
2675
2676	#First diagonal midpoint
2677	R0 = Bf.evaluate(0,0)
2678	assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
2679
2680	#Second diagonal midpoint
2681	R0 = Bf.evaluate(3,0)
2682	assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
2683
2684	#First diagonal midpoint
2685	R0 = Bf.evaluate(8,0)
2686	assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
2687
2688
2689
2690	#Cleanup
2691	os.remove(domain1.filename + '.' + domain1.format)
2692	#-------------------------------------------------------------
2693	if __name__ == "__main__":
2694	suite = unittest.makeSuite(TestCase,'test')
2695	#suite = unittest.makeSuite(TestCase,'test_spatio_temporal_boundary_2')
2696	runner = unittest.TextTestRunner()
2697	runner.run(suite)
2698

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

Download in other formats: