Context Navigation

source: anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py @ 3632

Last change on this file since 3632 was 3632, checked in by ole, 18 years ago
Changed gravity back to 9.8 as unit tests in shallow_water depend on it. Change this later.
File size: 117.4 KB

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

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

Download in other formats: