Context Navigation

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

Last change on this file since 3704 was 3704, checked in by steve, 18 years ago
Changed beta defaults to allow more limiting at wet/dry interface
File size: 128.6 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, alltrue, array, zeros, ones, Float, take
8	from anuga.utilities.numerical_tools import mean
9
10	from shallow_water_domain import *
11	from shallow_water_domain import flux_function_central as flux_function
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 xtest_catching_negative_heights(self):
671
672	#OBSOLETE
673
674	a = [0.0, 0.0]
675	b = [0.0, 2.0]
676	c = [2.0,0.0]
677	d = [0.0, 4.0]
678	e = [2.0, 2.0]
679	f = [4.0,0.0]
680
681	points = [a, b, c, d, e, f]
682	#bac, bce, ecf, dbe
683	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
684
685	domain = Domain(points, vertices)
686	val0 = 2.+2.0/3
687	val1 = 4.+4.0/3
688	val2 = 8.+2.0/3
689	val3 = 2.+8.0/3
690
691	zl=zr=4 #Too large
692	domain.set_quantity('elevation', zl*ones( (4,3) ))
693	domain.set_quantity('stage', [[val0, val0-1, val0-2],
694	[val1, val1+1, val1],
695	[val2, val2-2, val2],
696	[val3-0.5, val3, val3]])
697
698	#Should fail
699	try:
700	domain.check_integrity()
701	except:
702	pass
703
704
705
706	def test_get_wet_elements(self):
707
708	a = [0.0, 0.0]
709	b = [0.0, 2.0]
710	c = [2.0,0.0]
711	d = [0.0, 4.0]
712	e = [2.0, 2.0]
713	f = [4.0,0.0]
714
715	points = [a, b, c, d, e, f]
716	#bac, bce, ecf, dbe
717	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
718
719	domain = Domain(points, vertices)
720	val0 = 2.+2.0/3
721	val1 = 4.+4.0/3
722	val2 = 8.+2.0/3
723	val3 = 2.+8.0/3
724
725	zl=zr=5
726	domain.set_quantity('elevation', zl*ones( (4,3) ))
727	domain.set_quantity('stage', [[val0, val0-1, val0-2],
728	[val1, val1+1, val1],
729	[val2, val2-2, val2],
730	[val3-0.5, val3, val3]])
731
732
733
734	#print domain.get_quantity('elevation').get_values(location='centroids')
735	#print domain.get_quantity('stage').get_values(location='centroids')
736	domain.check_integrity()
737
738	indices = domain.get_wet_elements()
739	assert allclose(indices, [1,2])
740
741	indices = domain.get_wet_elements(indices=[0,1,3])
742	assert allclose(indices, [1])
743
744
745
746	def test_get_maximum_inundation_1(self):
747
748	a = [0.0, 0.0]
749	b = [0.0, 2.0]
750	c = [2.0,0.0]
751	d = [0.0, 4.0]
752	e = [2.0, 2.0]
753	f = [4.0,0.0]
754
755	points = [a, b, c, d, e, f]
756	#bac, bce, ecf, dbe
757	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
758
759	domain = Domain(points, vertices)
760
761	domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
762	domain.set_quantity('stage', 3)
763
764	domain.check_integrity()
765
766	indices = domain.get_wet_elements()
767	assert allclose(indices, [0])
768
769	q = domain.get_maximum_inundation_elevation()
770	assert allclose(q, domain.get_quantity('elevation').get_values(location='centroids')[0])
771
772	x, y = domain.get_maximum_inundation_location()
773	assert allclose([x, y], domain.get_centroid_coordinates()[0])
774
775
776	def test_get_maximum_inundation_2(self):
777	"""test_get_maximum_inundation_2(self)
778
779	Test multiple wet cells with same elevation
780	"""
781
782	a = [0.0, 0.0]
783	b = [0.0, 2.0]
784	c = [2.0,0.0]
785	d = [0.0, 4.0]
786	e = [2.0, 2.0]
787	f = [4.0,0.0]
788
789	points = [a, b, c, d, e, f]
790	#bac, bce, ecf, dbe
791	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
792
793	domain = Domain(points, vertices)
794
795	domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
796	domain.set_quantity('stage', 4.1)
797
798	domain.check_integrity()
799
800	indices = domain.get_wet_elements()
801	assert allclose(indices, [0,1,2])
802
803	q = domain.get_maximum_inundation_elevation()
804	assert allclose(q, 4)
805
806	x, y = domain.get_maximum_inundation_location()
807	assert allclose([x, y], domain.get_centroid_coordinates()[1])
808
809
810	def test_get_maximum_inundation_3(self):
811	"""test_get_maximum_inundation_3(self)
812
813	Test real runup example
814	"""
815
816	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
817
818	initial_runup_height = -0.4
819	final_runup_height = -0.3
820
821
822	#--------------------------------------------------------------
823	# Setup computational domain
824	#--------------------------------------------------------------
825	N = 5
826	points, vertices, boundary = rectangular_cross(N, N)
827	domain = Domain(points, vertices, boundary)
828
829	#--------------------------------------------------------------
830	# Setup initial conditions
831	#--------------------------------------------------------------
832	def topography(x,y):
833	return -x/2 # linear bed slope
834
835
836	domain.set_quantity('elevation', topography) # Use function for elevation
837	domain.set_quantity('friction', 0.) # Zero friction
838	domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
839
840
841	#--------------------------------------------------------------
842	# Setup boundary conditions
843	#--------------------------------------------------------------
844	Br = Reflective_boundary(domain) # Reflective wall
845	Bd = Dirichlet_boundary([final_runup_height, # Constant inflow
846	0,
847	0])
848
849	# All reflective to begin with (still water)
850	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
851
852
853	#--------------------------------------------------------------
854	# Test initial inundation height
855	#--------------------------------------------------------------
856
857	indices = domain.get_wet_elements()
858	z = domain.get_quantity('elevation').\
859	get_values(location='centroids', indices=indices)
860	assert alltrue(z<initial_runup_height)
861
862	q = domain.get_maximum_inundation_elevation()
863	assert allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
864
865	x, y = domain.get_maximum_inundation_location()
866
867	qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
868	assert allclose(q, qref)
869
870
871	wet_elements = domain.get_wet_elements()
872	wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
873	indices=wet_elements)
874	assert alltrue(wet_elevations<initial_runup_height)
875	assert allclose(wet_elevations, z)
876
877
878	#print domain.get_quantity('elevation').get_maximum_value(indices=wet_elements)
879	#print domain.get_quantity('elevation').get_maximum_location(indices=wet_elements)
880	#print domain.get_quantity('elevation').get_maximum_index(indices=wet_elements)
881
882
883	#--------------------------------------------------------------
884	# Let triangles adjust
885	#--------------------------------------------------------------
886	for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
887	pass
888
889
890	#--------------------------------------------------------------
891	# Test inundation height again
892	#--------------------------------------------------------------
893
894	indices = domain.get_wet_elements()
895	z = domain.get_quantity('elevation').\
896	get_values(location='centroids', indices=indices)
897
898	assert alltrue(z<initial_runup_height)
899
900	q = domain.get_maximum_inundation_elevation()
901	assert allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
902
903	x, y = domain.get_maximum_inundation_location()
904	qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
905	assert allclose(q, qref)
906
907
908	#--------------------------------------------------------------
909	# Update boundary to allow inflow
910	#--------------------------------------------------------------
911	domain.modify_boundary({'right': Bd})
912
913
914	#--------------------------------------------------------------
915	# Evolve system through time
916	#--------------------------------------------------------------
917	for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
918	#domain.write_time()
919	pass
920
921	#--------------------------------------------------------------
922	# Test inundation height again
923	#--------------------------------------------------------------
924
925	indices = domain.get_wet_elements()
926	z = domain.get_quantity('elevation').\
927	get_values(location='centroids', indices=indices)
928
929	assert alltrue(z<final_runup_height)
930
931	q = domain.get_maximum_inundation_elevation()
932	assert allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
933
934	x, y = domain.get_maximum_inundation_location()
935	qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
936	assert allclose(q, qref)
937
938
939	wet_elements = domain.get_wet_elements()
940	wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
941	indices=wet_elements)
942	assert alltrue(wet_elevations<final_runup_height)
943	assert allclose(wet_elevations, z)
944
945
946
947	#####################################################
948	def test_initial_condition(self):
949	"""Test that initial condition is output at time == 0
950	"""
951
952	from anuga.config import g
953	import copy
954
955	a = [0.0, 0.0]
956	b = [0.0, 2.0]
957	c = [2.0, 0.0]
958	d = [0.0, 4.0]
959	e = [2.0, 2.0]
960	f = [4.0, 0.0]
961
962	points = [a, b, c, d, e, f]
963	#bac, bce, ecf, dbe
964	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
965
966	domain = Domain(points, vertices)
967
968	#Set up for a gradient of (3,0) at mid triangle
969	def slope(x, y):
970	return 3*x
971
972	h = 0.1
973	def stage(x,y):
974	return slope(x,y)+h
975
976	domain.set_quantity('elevation', slope)
977	domain.set_quantity('stage', stage)
978
979	initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
980
981	domain.set_boundary({'exterior': Reflective_boundary(domain)})
982
983	#Evolution
984	for t in domain.evolve(yieldstep = 1.0, finaltime = 2.0):
985	stage = domain.quantities['stage'].vertex_values
986
987	if t == 0.0:
988	assert allclose(stage, initial_stage)
989	else:
990	assert not allclose(stage, initial_stage)
991
992	os.remove(domain.filename + '.sww')
993
994
995
996	#####################################################
997	def test_gravity(self):
998	#Assuming no friction
999
1000	from anuga.config import g
1001
1002	a = [0.0, 0.0]
1003	b = [0.0, 2.0]
1004	c = [2.0, 0.0]
1005	d = [0.0, 4.0]
1006	e = [2.0, 2.0]
1007	f = [4.0, 0.0]
1008
1009	points = [a, b, c, d, e, f]
1010	#bac, bce, ecf, dbe
1011	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1012
1013	domain = Domain(points, vertices)
1014
1015	#Set up for a gradient of (3,0) at mid triangle
1016	def slope(x, y):
1017	return 3*x
1018
1019	h = 0.1
1020	def stage(x,y):
1021	return slope(x,y)+h
1022
1023	domain.set_quantity('elevation', slope)
1024	domain.set_quantity('stage', stage)
1025
1026	for name in domain.conserved_quantities:
1027	assert allclose(domain.quantities[name].explicit_update, 0)
1028	assert allclose(domain.quantities[name].semi_implicit_update, 0)
1029
1030	domain.compute_forcing_terms()
1031
1032	assert allclose(domain.quantities['stage'].explicit_update, 0)
1033	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
1034	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
1035
1036
1037	def test_manning_friction(self):
1038	from anuga.config import g
1039
1040	a = [0.0, 0.0]
1041	b = [0.0, 2.0]
1042	c = [2.0, 0.0]
1043	d = [0.0, 4.0]
1044	e = [2.0, 2.0]
1045	f = [4.0, 0.0]
1046
1047	points = [a, b, c, d, e, f]
1048	#bac, bce, ecf, dbe
1049	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1050
1051	domain = Domain(points, vertices)
1052
1053	#Set up for a gradient of (3,0) at mid triangle
1054	def slope(x, y):
1055	return 3*x
1056
1057	h = 0.1
1058	def stage(x,y):
1059	return slope(x,y)+h
1060
1061	eta = 0.07
1062	domain.set_quantity('elevation', slope)
1063	domain.set_quantity('stage', stage)
1064	domain.set_quantity('friction', eta)
1065
1066	for name in domain.conserved_quantities:
1067	assert allclose(domain.quantities[name].explicit_update, 0)
1068	assert allclose(domain.quantities[name].semi_implicit_update, 0)
1069
1070	domain.compute_forcing_terms()
1071
1072	assert allclose(domain.quantities['stage'].explicit_update, 0)
1073	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
1074	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
1075
1076	assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1077	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
1078	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
1079
1080	#Create some momentum for friction to work with
1081	domain.set_quantity('xmomentum', 1)
1082	S = -g * eta2 / h(7.0/3)
1083
1084	domain.compute_forcing_terms()
1085	assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1086	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
1087	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
1088
1089	#A more complex example
1090	domain.quantities['stage'].semi_implicit_update[:] = 0.0
1091	domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
1092	domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
1093
1094	domain.set_quantity('xmomentum', 3)
1095	domain.set_quantity('ymomentum', 4)
1096
1097	S = -g * eta*2 5 / h**(7.0/3)
1098
1099
1100	domain.compute_forcing_terms()
1101
1102	assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
1103	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
1104	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
1105
1106	def test_constant_wind_stress(self):
1107	from anuga.config import rho_a, rho_w, eta_w
1108	from math import pi, cos, sin, sqrt
1109
1110	a = [0.0, 0.0]
1111	b = [0.0, 2.0]
1112	c = [2.0, 0.0]
1113	d = [0.0, 4.0]
1114	e = [2.0, 2.0]
1115	f = [4.0, 0.0]
1116
1117	points = [a, b, c, d, e, f]
1118	#bac, bce, ecf, dbe
1119	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1120
1121
1122	domain = Domain(points, vertices)
1123
1124	#Flat surface with 1m of water
1125	domain.set_quantity('elevation', 0)
1126	domain.set_quantity('stage', 1.0)
1127	domain.set_quantity('friction', 0)
1128
1129	Br = Reflective_boundary(domain)
1130	domain.set_boundary({'exterior': Br})
1131
1132	#Setup only one forcing term, constant wind stress
1133	s = 100
1134	phi = 135
1135	domain.forcing_terms = []
1136	domain.forcing_terms.append( Wind_stress(s, phi) )
1137
1138	domain.compute_forcing_terms()
1139
1140
1141	const = eta_w*rho_a/rho_w
1142
1143	#Convert to radians
1144	phi = phi*pi/180
1145
1146	#Compute velocity vector (u, v)
1147	u = s*cos(phi)
1148	v = s*sin(phi)
1149
1150	#Compute wind stress
1151	S = const * sqrt(u2 + v2)
1152
1153	assert allclose(domain.quantities['stage'].explicit_update, 0)
1154	assert allclose(domain.quantities['xmomentum'].explicit_update, S*u)
1155	assert allclose(domain.quantities['ymomentum'].explicit_update, S*v)
1156
1157
1158	def test_variable_wind_stress(self):
1159	from anuga.config import rho_a, rho_w, eta_w
1160	from math import pi, cos, sin, sqrt
1161
1162	a = [0.0, 0.0]
1163	b = [0.0, 2.0]
1164	c = [2.0, 0.0]
1165	d = [0.0, 4.0]
1166	e = [2.0, 2.0]
1167	f = [4.0, 0.0]
1168
1169	points = [a, b, c, d, e, f]
1170	#bac, bce, ecf, dbe
1171	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1172
1173	domain = Domain(points, vertices)
1174
1175	#Flat surface with 1m of water
1176	domain.set_quantity('elevation', 0)
1177	domain.set_quantity('stage', 1.0)
1178	domain.set_quantity('friction', 0)
1179
1180	Br = Reflective_boundary(domain)
1181	domain.set_boundary({'exterior': Br})
1182
1183
1184	domain.time = 5.54 #Take a random time (not zero)
1185
1186	#Setup only one forcing term, constant wind stress
1187	s = 100
1188	phi = 135
1189	domain.forcing_terms = []
1190	domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
1191
1192	domain.compute_forcing_terms()
1193
1194	#Compute reference solution
1195	const = eta_w*rho_a/rho_w
1196
1197	N = domain.number_of_elements
1198
1199	xc = domain.get_centroid_coordinates()
1200	t = domain.time
1201
1202	x = xc[:,0]
1203	y = xc[:,1]
1204	s_vec = speed(t,x,y)
1205	phi_vec = angle(t,x,y)
1206
1207
1208	for k in range(N):
1209	#Convert to radians
1210	phi = phi_vec[k]*pi/180
1211	s = s_vec[k]
1212
1213	#Compute velocity vector (u, v)
1214	u = s*cos(phi)
1215	v = s*sin(phi)
1216
1217	#Compute wind stress
1218	S = const * sqrt(u2 + v2)
1219
1220	assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1221	assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1222	assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1223
1224
1225
1226
1227	def test_windfield_from_file(self):
1228	from anuga.config import rho_a, rho_w, eta_w
1229	from math import pi, cos, sin, sqrt
1230	from anuga.config import time_format
1231	from anuga.abstract_2d_finite_volumes.util import file_function
1232	import time
1233
1234
1235	a = [0.0, 0.0]
1236	b = [0.0, 2.0]
1237	c = [2.0, 0.0]
1238	d = [0.0, 4.0]
1239	e = [2.0, 2.0]
1240	f = [4.0, 0.0]
1241
1242	points = [a, b, c, d, e, f]
1243	#bac, bce, ecf, dbe
1244	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1245
1246	domain = Domain(points, vertices)
1247
1248	#Flat surface with 1m of water
1249	domain.set_quantity('elevation', 0)
1250	domain.set_quantity('stage', 1.0)
1251	domain.set_quantity('friction', 0)
1252
1253	Br = Reflective_boundary(domain)
1254	domain.set_boundary({'exterior': Br})
1255
1256
1257	domain.time = 7 #Take a time that is represented in file (not zero)
1258
1259	#Write wind stress file (ensure that domain.time is covered)
1260	#Take x=1 and y=0
1261	filename = 'test_windstress_from_file'
1262	start = time.mktime(time.strptime('2000', '%Y'))
1263	fid = open(filename + '.txt', 'w')
1264	dt = 1 #One second interval
1265	t = 0.0
1266	while t <= 10.0:
1267	t_string = time.strftime(time_format, time.gmtime(t+start))
1268
1269	fid.write('%s, %f %f\n' %(t_string,
1270	speed(t,[1],[0])[0],
1271	angle(t,[1],[0])[0]))
1272	t += dt
1273
1274	fid.close()
1275
1276
1277	#Convert ASCII file to NetCDF (Which is what we really like!)
1278	from data_manager import timefile2netcdf
1279	timefile2netcdf(filename)
1280	os.remove(filename + '.txt')
1281
1282
1283	#Setup wind stress
1284	F = file_function(filename + '.tms', quantities = ['Attribute0',
1285	'Attribute1'])
1286	os.remove(filename + '.tms')
1287
1288
1289	#print 'F(5)', F(5)
1290
1291	#print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
1292
1293	#print dir(F)
1294	#print F.T
1295	#print F.precomputed_values
1296	#
1297	#F = file_function(filename + '.txt')
1298	#
1299	#print dir(F)
1300	#print F.T
1301	#print F.Q
1302
1303	W = Wind_stress(F)
1304
1305	domain.forcing_terms = []
1306	domain.forcing_terms.append(W)
1307
1308	domain.compute_forcing_terms()
1309
1310	#Compute reference solution
1311	const = eta_w*rho_a/rho_w
1312
1313	N = domain.number_of_elements
1314
1315	t = domain.time
1316
1317	s = speed(t,[1],[0])[0]
1318	phi = angle(t,[1],[0])[0]
1319
1320	#Convert to radians
1321	phi = phi*pi/180
1322
1323
1324	#Compute velocity vector (u, v)
1325	u = s*cos(phi)
1326	v = s*sin(phi)
1327
1328	#Compute wind stress
1329	S = const * sqrt(u2 + v2)
1330
1331	for k in range(N):
1332	assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1333	assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1334	assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1335
1336
1337
1338
1339	def test_wind_stress_error_condition(self):
1340	"""Test that windstress reacts properly when forcing functions
1341	are wrong - e.g. returns a scalar
1342	"""
1343
1344	from anuga.config import rho_a, rho_w, eta_w
1345	from math import pi, cos, sin, sqrt
1346
1347	a = [0.0, 0.0]
1348	b = [0.0, 2.0]
1349	c = [2.0, 0.0]
1350	d = [0.0, 4.0]
1351	e = [2.0, 2.0]
1352	f = [4.0, 0.0]
1353
1354	points = [a, b, c, d, e, f]
1355	#bac, bce, ecf, dbe
1356	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1357
1358	domain = Domain(points, vertices)
1359
1360	#Flat surface with 1m of water
1361	domain.set_quantity('elevation', 0)
1362	domain.set_quantity('stage', 1.0)
1363	domain.set_quantity('friction', 0)
1364
1365	Br = Reflective_boundary(domain)
1366	domain.set_boundary({'exterior': Br})
1367
1368
1369	domain.time = 5.54 #Take a random time (not zero)
1370
1371	#Setup only one forcing term, bad func
1372	domain.forcing_terms = []
1373
1374	try:
1375	domain.forcing_terms.append(Wind_stress(s = scalar_func,
1376	phi = angle))
1377	except AssertionError:
1378	pass
1379	else:
1380	msg = 'Should have raised exception'
1381	raise msg
1382
1383
1384	try:
1385	domain.forcing_terms.append(Wind_stress(s = speed,
1386	phi = scalar_func))
1387	except AssertionError:
1388	pass
1389	else:
1390	msg = 'Should have raised exception'
1391	raise msg
1392
1393	try:
1394	domain.forcing_terms.append(Wind_stress(s = speed,
1395	phi = 'xx'))
1396	except:
1397	pass
1398	else:
1399	msg = 'Should have raised exception'
1400	raise msg
1401
1402
1403	#####################################################
1404	def test_first_order_extrapolator_const_z(self):
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	val0 = 2.+2.0/3
1419	val1 = 4.+4.0/3
1420	val2 = 8.+2.0/3
1421	val3 = 2.+8.0/3
1422
1423	zl=zr=-3.75 #Assume constant bed (must be less than stage)
1424	domain.set_quantity('elevation', zl*ones( (4,3) ))
1425	domain.set_quantity('stage', [[val0, val0-1, val0-2],
1426	[val1, val1+1, val1],
1427	[val2, val2-2, val2],
1428	[val3-0.5, val3, val3]])
1429
1430
1431
1432	domain._order_ = 1
1433	domain.distribute_to_vertices_and_edges()
1434
1435	#Check that centroid values were distributed to vertices
1436	C = domain.quantities['stage'].centroid_values
1437	for i in range(3):
1438	assert allclose( domain.quantities['stage'].vertex_values[:,i], C)
1439
1440
1441	def test_first_order_limiter_variable_z(self):
1442	#Check that first order limiter follows bed_slope
1443	from Numeric import alltrue, greater_equal
1444	from anuga.config import epsilon
1445
1446	a = [0.0, 0.0]
1447	b = [0.0, 2.0]
1448	c = [2.0,0.0]
1449	d = [0.0, 4.0]
1450	e = [2.0, 2.0]
1451	f = [4.0,0.0]
1452
1453	points = [a, b, c, d, e, f]
1454	#bac, bce, ecf, dbe
1455	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1456
1457	domain = Domain(points, vertices)
1458	val0 = 2.+2.0/3
1459	val1 = 4.+4.0/3
1460	val2 = 8.+2.0/3
1461	val3 = 2.+8.0/3
1462
1463	domain.set_quantity('elevation', [[0,0,0], [6,0,0],
1464	[6,6,6], [6,6,6]])
1465	domain.set_quantity('stage', [[val0, val0, val0],
1466	[val1, val1, val1],
1467	[val2, val2, val2],
1468	[val3, val3, val3]])
1469
1470	E = domain.quantities['elevation'].vertex_values
1471	L = domain.quantities['stage'].vertex_values
1472
1473
1474	#Check that some stages are not above elevation (within eps)
1475	#- so that the limiter has something to work with
1476	assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
1477
1478	domain._order_ = 1
1479	domain.distribute_to_vertices_and_edges()
1480
1481	#Check that all stages are above elevation (within eps)
1482	assert alltrue(alltrue(greater_equal(L,E-epsilon)))
1483
1484
1485	#####################################################
1486	def test_distribute_basic(self):
1487	#Using test data generated by abstract_2d_finite_volumes-2
1488	#Assuming no friction and flat bed (0.0)
1489
1490	a = [0.0, 0.0]
1491	b = [0.0, 2.0]
1492	c = [2.0, 0.0]
1493	d = [0.0, 4.0]
1494	e = [2.0, 2.0]
1495	f = [4.0, 0.0]
1496
1497	points = [a, b, c, d, e, f]
1498	#bac, bce, ecf, dbe
1499	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1500
1501	domain = Domain(points, vertices)
1502
1503	val0 = 2.
1504	val1 = 4.
1505	val2 = 8.
1506	val3 = 2.
1507
1508	domain.set_quantity('stage', [val0, val1, val2, val3],
1509	location='centroids')
1510	L = domain.quantities['stage'].vertex_values
1511
1512	#First order
1513	domain._order_ = 1
1514	domain.distribute_to_vertices_and_edges()
1515	assert allclose(L[1], val1)
1516
1517	#Second order
1518	domain._order_ = 2
1519	domain.beta_w = 0.9
1520	domain.beta_w_dry = 0.9
1521	domain.beta_uh = 0.9
1522	domain.beta_uh_dry = 0.9
1523	domain.beta_vh = 0.9
1524	domain.beta_vh_dry = 0.9
1525	domain.distribute_to_vertices_and_edges()
1526	assert allclose(L[1], [2.2, 4.9, 4.9])
1527
1528
1529
1530	def test_distribute_away_from_bed(self):
1531	#Using test data generated by abstract_2d_finite_volumes-2
1532	#Assuming no friction and flat bed (0.0)
1533
1534	a = [0.0, 0.0]
1535	b = [0.0, 2.0]
1536	c = [2.0, 0.0]
1537	d = [0.0, 4.0]
1538	e = [2.0, 2.0]
1539	f = [4.0, 0.0]
1540
1541	points = [a, b, c, d, e, f]
1542	#bac, bce, ecf, dbe
1543	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1544
1545	domain = Domain(points, vertices)
1546	L = domain.quantities['stage'].vertex_values
1547
1548	def stage(x,y):
1549	return x**2
1550
1551	domain.set_quantity('stage', stage, location='centroids')
1552
1553	a, b = domain.quantities['stage'].compute_gradients()
1554	assert allclose(a[1], 3.33333334)
1555	assert allclose(b[1], 0.0)
1556
1557	domain._order_ = 1
1558	domain.distribute_to_vertices_and_edges()
1559	assert allclose(L[1], 1.77777778)
1560
1561	domain._order_ = 2
1562	domain.beta_w = 0.9
1563	domain.beta_w_dry = 0.9
1564	domain.beta_uh = 0.9
1565	domain.beta_uh_dry = 0.9
1566	domain.beta_vh = 0.9
1567	domain.beta_vh_dry = 0.9
1568	domain.distribute_to_vertices_and_edges()
1569	assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
1570
1571
1572
1573	def test_distribute_away_from_bed1(self):
1574	#Using test data generated by abstract_2d_finite_volumes-2
1575	#Assuming no friction and flat bed (0.0)
1576
1577	a = [0.0, 0.0]
1578	b = [0.0, 2.0]
1579	c = [2.0, 0.0]
1580	d = [0.0, 4.0]
1581	e = [2.0, 2.0]
1582	f = [4.0, 0.0]
1583
1584	points = [a, b, c, d, e, f]
1585	#bac, bce, ecf, dbe
1586	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1587
1588	domain = Domain(points, vertices)
1589	L = domain.quantities['stage'].vertex_values
1590
1591	def stage(x,y):
1592	return x4+y2
1593
1594	domain.set_quantity('stage', stage, location='centroids')
1595	#print domain.quantities['stage'].centroid_values
1596
1597	a, b = domain.quantities['stage'].compute_gradients()
1598	assert allclose(a[1], 25.18518519)
1599	assert allclose(b[1], 3.33333333)
1600
1601	domain._order_ = 1
1602	domain.distribute_to_vertices_and_edges()
1603	assert allclose(L[1], 4.9382716)
1604
1605	domain._order_ = 2
1606	domain.beta_w = 0.9
1607	domain.beta_w_dry = 0.9
1608	domain.beta_uh = 0.9
1609	domain.beta_uh_dry = 0.9
1610	domain.beta_vh = 0.9
1611	domain.beta_vh_dry = 0.9
1612	domain.distribute_to_vertices_and_edges()
1613	assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
1614
1615
1616
1617	def test_distribute_near_bed(self):
1618	#Using test data generated by abstract_2d_finite_volumes-2
1619	#Assuming no friction and flat bed (0.0)
1620
1621	a = [0.0, 0.0]
1622	b = [0.0, 2.0]
1623	c = [2.0, 0.0]
1624	d = [0.0, 4.0]
1625	e = [2.0, 2.0]
1626	f = [4.0, 0.0]
1627
1628	points = [a, b, c, d, e, f]
1629	#bac, bce, ecf, dbe
1630	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1631
1632	domain = Domain(points, vertices)
1633
1634
1635	#Set up for a gradient of (3,0) at mid triangle
1636	def slope(x, y):
1637	return 10*x
1638
1639	h = 0.1
1640	def stage(x,y):
1641	return slope(x,y)+h
1642
1643	domain.set_quantity('elevation', slope)
1644	domain.set_quantity('stage', stage, location='centroids')
1645
1646	#print domain.quantities['elevation'].centroid_values
1647	#print domain.quantities['stage'].centroid_values
1648
1649	E = domain.quantities['elevation'].vertex_values
1650	L = domain.quantities['stage'].vertex_values
1651
1652	#print E
1653	domain._order_ = 1
1654	domain.distribute_to_vertices_and_edges()
1655	##assert allclose(L[1], [0.19999999, 20.05, 20.05])
1656	assert allclose(L[1], [0.1, 20.1, 20.1])
1657
1658	domain._order_ = 2
1659	domain.distribute_to_vertices_and_edges()
1660	assert allclose(L[1], [0.1, 20.1, 20.1])
1661
1662	def test_distribute_near_bed1(self):
1663	#Using test data generated by abstract_2d_finite_volumes-2
1664	#Assuming no friction and flat bed (0.0)
1665
1666	a = [0.0, 0.0]
1667	b = [0.0, 2.0]
1668	c = [2.0, 0.0]
1669	d = [0.0, 4.0]
1670	e = [2.0, 2.0]
1671	f = [4.0, 0.0]
1672
1673	points = [a, b, c, d, e, f]
1674	#bac, bce, ecf, dbe
1675	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1676
1677	domain = Domain(points, vertices)
1678
1679
1680	#Set up for a gradient of (3,0) at mid triangle
1681	def slope(x, y):
1682	return x4+y2
1683
1684	h = 0.1
1685	def stage(x,y):
1686	return slope(x,y)+h
1687
1688	domain.set_quantity('elevation', slope)
1689	domain.set_quantity('stage', stage)
1690
1691	#print domain.quantities['elevation'].centroid_values
1692	#print domain.quantities['stage'].centroid_values
1693
1694	E = domain.quantities['elevation'].vertex_values
1695	L = domain.quantities['stage'].vertex_values
1696
1697	#print E
1698	domain._order_ = 1
1699	domain.distribute_to_vertices_and_edges()
1700	##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
1701	assert allclose(L[1], [4.1, 16.1, 20.1])
1702
1703	domain._order_ = 2
1704	domain.distribute_to_vertices_and_edges()
1705	assert allclose(L[1], [4.1, 16.1, 20.1])
1706
1707
1708
1709	def test_second_order_distribute_real_data(self):
1710	#Using test data generated by abstract_2d_finite_volumes-2
1711	#Assuming no friction and flat bed (0.0)
1712
1713	a = [0.0, 0.0]
1714	b = [0.0, 1.0/5]
1715	c = [0.0, 2.0/5]
1716	d = [1.0/5, 0.0]
1717	e = [1.0/5, 1.0/5]
1718	f = [1.0/5, 2.0/5]
1719	g = [2.0/5, 2.0/5]
1720
1721	points = [a, b, c, d, e, f, g]
1722	#bae, efb, cbf, feg
1723	vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
1724
1725	domain = Domain(points, vertices)
1726
1727	def slope(x, y):
1728	return -x/3
1729
1730	domain.set_quantity('elevation', slope)
1731	domain.set_quantity('stage',
1732	[0.01298164, 0.00365611,
1733	0.01440365, -0.0381856437096],
1734	location='centroids')
1735	domain.set_quantity('xmomentum',
1736	[0.00670439, 0.01263789,
1737	0.00647805, 0.0178180740668],
1738	location='centroids')
1739	domain.set_quantity('ymomentum',
1740	[-7.23510980e-004, -6.30413883e-005,
1741	6.30413883e-005, 0.000200907255866],
1742	location='centroids')
1743
1744	E = domain.quantities['elevation'].vertex_values
1745	L = domain.quantities['stage'].vertex_values
1746	X = domain.quantities['xmomentum'].vertex_values
1747	Y = domain.quantities['ymomentum'].vertex_values
1748
1749	#print E
1750	domain._order_ = 2
1751	domain.beta_w = 0.9
1752	domain.beta_w_dry = 0.9
1753	domain.beta_uh = 0.9
1754	domain.beta_uh_dry = 0.9
1755	domain.beta_vh = 0.9
1756	domain.beta_vh_dry = 0.9
1757	domain.beta_h = 0.0 #Use first order in h-limiter
1758	domain.distribute_to_vertices_and_edges()
1759
1760	#print L[1,:]
1761	#print X[1,:]
1762	#print Y[1,:]
1763
1764	assert allclose(L[1,:], [-0.00825735775384,
1765	-0.00801881482869,
1766	0.0272445025825])
1767	assert allclose(X[1,:], [0.0143507718962,
1768	0.0142502147066,
1769	0.00931268339717])
1770	assert allclose(Y[1,:], [-0.000117062180693,
1771	7.94434448109e-005,
1772	-0.000151505429018])
1773
1774
1775
1776	def test_balance_deep_and_shallow(self):
1777	"""Test that balanced limiters preserve conserved quantites.
1778	"""
1779	import copy
1780
1781	a = [0.0, 0.0]
1782	b = [0.0, 2.0]
1783	c = [2.0, 0.0]
1784	d = [0.0, 4.0]
1785	e = [2.0, 2.0]
1786	f = [4.0, 0.0]
1787
1788	points = [a, b, c, d, e, f]
1789
1790	#bac, bce, ecf, dbe
1791	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1792
1793	mesh = Domain(points, elements)
1794	mesh.check_integrity()
1795
1796	#Create a deliberate overshoot
1797	mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1798	mesh.set_quantity('elevation', 0) #Flat bed
1799	stage = mesh.quantities['stage']
1800
1801	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
1802
1803	#Limit
1804	mesh.distribute_to_vertices_and_edges()
1805
1806	#Assert that quantities are conserved
1807	from Numeric import sum
1808	for k in range(mesh.number_of_elements):
1809	assert allclose (ref_centroid_values[k],
1810	sum(stage.vertex_values[k,:])/3)
1811
1812
1813	#Now try with a non-flat bed - closely hugging initial stage in places
1814	#This will create alphas in the range [0, 0.478260, 1]
1815	mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1816	mesh.set_quantity('elevation', [[0,0,0],
1817	[1.8,1.9,5.9],
1818	[4.6,0,0],
1819	[0,2,4]])
1820	stage = mesh.quantities['stage']
1821
1822	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
1823	ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
1824
1825	#Limit
1826	mesh.distribute_to_vertices_and_edges()
1827
1828
1829	#Assert that all vertex quantities have changed
1830	for k in range(mesh.number_of_elements):
1831	#print ref_vertex_values[k,:], stage.vertex_values[k,:]
1832	assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
1833	#and assert that quantities are still conserved
1834	from Numeric import sum
1835	for k in range(mesh.number_of_elements):
1836	assert allclose (ref_centroid_values[k],
1837	sum(stage.vertex_values[k,:])/3)
1838
1839
1840	#Also check that Python and C version produce the same
1841	assert allclose (stage.vertex_values,
1842	[[2,2,2],
1843	[1.93333333, 2.03333333, 6.03333333],
1844	[6.93333333, 4.53333333, 4.53333333],
1845	[5.33333333, 5.33333333, 5.33333333]])
1846
1847
1848
1849
1850	def test_conservation_1(self):
1851	"""Test that stage is conserved globally
1852
1853	This one uses a flat bed, reflective bdries and a suitable
1854	initial condition
1855	"""
1856	from mesh_factory import rectangular
1857	from Numeric import array
1858
1859	#Create basic mesh
1860	points, vertices, boundary = rectangular(6, 6)
1861
1862	#Create shallow water domain
1863	domain = Domain(points, vertices, boundary)
1864	domain.smooth = False
1865	domain.default_order=2
1866
1867	#IC
1868	def x_slope(x, y):
1869	return x/3
1870
1871	domain.set_quantity('elevation', 0)
1872	domain.set_quantity('friction', 0)
1873	domain.set_quantity('stage', x_slope)
1874
1875	# Boundary conditions (reflective everywhere)
1876	Br = Reflective_boundary(domain)
1877	domain.set_boundary({'left': Br, 'right': Br, '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	#print initial_xmom
1887
1888	#Evolution
1889	for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
1890	volume = domain.quantities['stage'].get_integral()
1891	assert allclose (volume, initial_volume)
1892
1893	#I don't believe that the total momentum should be the same
1894	#It starts with zero and ends with zero though
1895	#xmom = domain.quantities['xmomentum'].get_integral()
1896	#print xmom
1897	#assert allclose (xmom, initial_xmom)
1898
1899	os.remove(domain.filename + '.sww')
1900
1901
1902	def test_conservation_2(self):
1903	"""Test that stage is conserved globally
1904
1905	This one uses a slopy bed, reflective bdries and a suitable
1906	initial condition
1907	"""
1908	from mesh_factory import rectangular
1909	from Numeric import array
1910
1911	#Create basic mesh
1912	points, vertices, boundary = rectangular(6, 6)
1913
1914	#Create shallow water domain
1915	domain = Domain(points, vertices, boundary)
1916	domain.smooth = False
1917	domain.default_order=2
1918
1919	#IC
1920	def x_slope(x, y):
1921	return x/3
1922
1923	domain.set_quantity('elevation', x_slope)
1924	domain.set_quantity('friction', 0)
1925	domain.set_quantity('stage', 0.4) #Steady
1926
1927	# Boundary conditions (reflective everywhere)
1928	Br = Reflective_boundary(domain)
1929	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1930
1931	domain.check_integrity()
1932
1933	#domain.visualise = True #If you want to take a sticky beak
1934
1935	initial_volume = domain.quantities['stage'].get_integral()
1936	initial_xmom = domain.quantities['xmomentum'].get_integral()
1937
1938	#print initial_xmom
1939
1940	#Evolution
1941	for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
1942	volume = domain.quantities['stage'].get_integral()
1943	assert allclose (volume, initial_volume)
1944
1945	#FIXME: What would we expect from momentum
1946	#xmom = domain.quantities['xmomentum'].get_integral()
1947	#print xmom
1948	#assert allclose (xmom, initial_xmom)
1949
1950	os.remove(domain.filename + '.sww')
1951
1952	def test_conservation_3(self):
1953	"""Test that stage is conserved globally
1954
1955	This one uses a larger grid, convoluted bed, reflective bdries and a suitable
1956	initial condition
1957	"""
1958	from mesh_factory import rectangular
1959	from Numeric import array
1960
1961	#Create basic mesh
1962	points, vertices, boundary = rectangular(2, 1)
1963
1964	#Create shallow water domain
1965	domain = Domain(points, vertices, boundary)
1966	domain.smooth = False
1967	domain.default_order = 2
1968	domain.beta_h = 0.2
1969	domain.set_quantities_to_be_stored(['stage'])
1970
1971	#IC
1972	def x_slope(x, y):
1973	z = 0*x
1974	for i in range(len(x)):
1975	if x[i] < 0.3:
1976	z[i] = x[i]/3
1977	if 0.3 <= x[i] < 0.5:
1978	z[i] = -0.5
1979	if 0.5 <= x[i] < 0.7:
1980	z[i] = 0.39
1981	if 0.7 <= x[i]:
1982	z[i] = x[i]/3
1983	return z
1984
1985
1986
1987	domain.set_quantity('elevation', x_slope)
1988	domain.set_quantity('friction', 0)
1989	domain.set_quantity('stage', 0.4) #Steady
1990
1991	# Boundary conditions (reflective everywhere)
1992	Br = Reflective_boundary(domain)
1993	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1994
1995	domain.check_integrity()
1996
1997	#domain.visualise = True #If you want to take a sticky beak
1998
1999	initial_volume = domain.quantities['stage'].get_integral()
2000	initial_xmom = domain.quantities['xmomentum'].get_integral()
2001
2002	import copy
2003	ref_centroid_values =\
2004	copy.copy(domain.quantities['stage'].centroid_values)
2005
2006	#print 'ORG', domain.quantities['stage'].centroid_values
2007	domain.distribute_to_vertices_and_edges()
2008
2009
2010	#print domain.quantities['stage'].centroid_values
2011	assert allclose(domain.quantities['stage'].centroid_values,
2012	ref_centroid_values)
2013
2014
2015	#Check that initial limiter doesn't violate cons quan
2016	assert allclose (domain.quantities['stage'].get_integral(),
2017	initial_volume)
2018
2019	#Evolution
2020	for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
2021	volume = domain.quantities['stage'].get_integral()
2022	#print t, volume, initial_volume
2023	assert allclose (volume, initial_volume)
2024
2025	os.remove(domain.filename + '.sww')
2026
2027	def test_conservation_4(self):
2028	"""Test that stage is conserved globally
2029
2030	This one uses a larger grid, convoluted bed, reflective bdries and a suitable
2031	initial condition
2032	"""
2033	from mesh_factory import rectangular
2034	from Numeric import array
2035
2036	#Create basic mesh
2037	points, vertices, boundary = rectangular(6, 6)
2038
2039	#Create shallow water domain
2040	domain = Domain(points, vertices, boundary)
2041	domain.smooth = False
2042	domain.default_order=2
2043	domain.beta_h = 0.0
2044
2045	#IC
2046	def x_slope(x, y):
2047	z = 0*x
2048	for i in range(len(x)):
2049	if x[i] < 0.3:
2050	z[i] = x[i]/3
2051	if 0.3 <= x[i] < 0.5:
2052	z[i] = -0.5
2053	if 0.5 <= x[i] < 0.7:
2054	#z[i] = 0.3 #OK with beta == 0.2
2055	z[i] = 0.34 #OK with beta == 0.0
2056	#z[i] = 0.35#Fails after 80 timesteps with an error
2057	#of the order 1.0e-5
2058	if 0.7 <= x[i]:
2059	z[i] = x[i]/3
2060	return z
2061
2062
2063
2064	domain.set_quantity('elevation', x_slope)
2065	domain.set_quantity('friction', 0)
2066	domain.set_quantity('stage', 0.4) #Steady
2067
2068	# Boundary conditions (reflective everywhere)
2069	Br = Reflective_boundary(domain)
2070	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2071
2072	domain.check_integrity()
2073
2074	#domain.visualise = True #If you want to take a sticky beak
2075
2076	initial_volume = domain.quantities['stage'].get_integral()
2077	initial_xmom = domain.quantities['xmomentum'].get_integral()
2078
2079	import copy
2080	ref_centroid_values =\
2081	copy.copy(domain.quantities['stage'].centroid_values)
2082
2083	#Test limiter by itself
2084	domain.distribute_to_vertices_and_edges()
2085
2086	#Check that initial limiter doesn't violate cons quan
2087	assert allclose (domain.quantities['stage'].get_integral(),
2088	initial_volume)
2089	#NOTE: This would fail if any initial stage was less than the
2090	#corresponding bed elevation - but that is reasonable.
2091
2092
2093	#Evolution
2094	for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
2095	volume = domain.quantities['stage'].get_integral()
2096
2097	#print t, volume, initial_volume
2098
2099
2100	#if not allclose (volume, initial_volume):
2101	# print 't==4.05'
2102	# for k in range(domain.number_of_elements):
2103	# pass
2104	# print domain.quantities['stage'].centroid_values[k] -\
2105	# ref_centroid_values[k]
2106
2107	assert allclose (volume, initial_volume)
2108
2109
2110	os.remove(domain.filename + '.sww')
2111
2112
2113	def test_conservation_5(self):
2114	"""Test that momentum is conserved globally in
2115	steady state scenario
2116
2117	This one uses a slopy bed, dirichlet and reflective bdries
2118	"""
2119	from mesh_factory import rectangular
2120	from Numeric import array
2121
2122	#Create basic mesh
2123	points, vertices, boundary = rectangular(6, 6)
2124
2125	#Create shallow water domain
2126	domain = Domain(points, vertices, boundary)
2127	domain.smooth = False
2128	domain.default_order=2
2129
2130	#IC
2131	def x_slope(x, y):
2132	return x/3
2133
2134	domain.set_quantity('elevation', x_slope)
2135	domain.set_quantity('friction', 0)
2136	domain.set_quantity('stage', 0.4) #Steady
2137
2138	# Boundary conditions (reflective everywhere)
2139	Br = Reflective_boundary(domain)
2140	Bleft = Dirichlet_boundary([0.5,0,0])
2141	Bright = Dirichlet_boundary([0.1,0,0])
2142	domain.set_boundary({'left': Bleft, 'right': Bright,
2143	'top': Br, 'bottom': Br})
2144
2145	domain.check_integrity()
2146
2147	#domain.visualise = True #If you want to take a sticky beak
2148
2149	initial_volume = domain.quantities['stage'].get_integral()
2150	initial_xmom = domain.quantities['xmomentum'].get_integral()
2151
2152
2153	#Evolution
2154	for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
2155	stage = domain.quantities['stage'].get_integral()
2156	xmom = domain.quantities['xmomentum'].get_integral()
2157	ymom = domain.quantities['ymomentum'].get_integral()
2158
2159	if allclose(t, 6): #Steady state reached
2160	steady_xmom = domain.quantities['xmomentum'].get_integral()
2161	steady_ymom = domain.quantities['ymomentum'].get_integral()
2162	steady_stage = domain.quantities['stage'].get_integral()
2163
2164	if t > 6:
2165	#print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
2166	assert allclose(xmom, steady_xmom)
2167	assert allclose(ymom, steady_ymom)
2168	assert allclose(stage, steady_stage)
2169
2170
2171	os.remove(domain.filename + '.sww')
2172
2173
2174
2175
2176
2177	def test_conservation_real(self):
2178	"""Test that momentum is conserved globally
2179
2180	Stephen finally made a test that revealed the problem.
2181	This test failed with code prior to 25 July 2005
2182	"""
2183
2184	yieldstep = 0.01
2185	finaltime = 0.05
2186	min_depth = 1.0e-2
2187
2188
2189	import sys
2190	from os import sep; sys.path.append('..'+sep+'abstract_2d_finite_volumes')
2191	from mesh_factory import rectangular
2192
2193
2194	#Create shallow water domain
2195	points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
2196	domain = Domain(points, vertices, boundary)
2197	domain.smooth = False
2198	domain.visualise = False
2199	domain.default_order = 1
2200	domain.minimum_allowed_height = min_depth
2201
2202	# Set initial condition
2203	class Set_IC:
2204	"""Set an initial condition with a constant value, for x0<x<x1
2205	"""
2206
2207	def __init__(self, x0=0.25, x1=0.5, h=1.0):
2208	self.x0 = x0
2209	self.x1 = x1
2210	self.h = h
2211
2212	def __call__(self, x, y):
2213	return self.h*((x>self.x0)&(x<self.x1))
2214
2215
2216	domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
2217
2218
2219	#Boundaries
2220	R = Reflective_boundary(domain)
2221	domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R})
2222
2223	ref = domain.quantities['stage'].get_integral()
2224
2225	# Evolution
2226	for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
2227	pass
2228	#print 'Integral stage = ',\
2229	# domain.quantities['stage'].get_integral(),\
2230	# ' Time = ',domain.time
2231
2232
2233	now = domain.quantities['stage'].get_integral()
2234
2235	msg = 'Stage not conserved: was %f, now %f' %(ref, now)
2236	assert allclose(ref, now), msg
2237
2238	os.remove(domain.filename + '.sww')
2239
2240	def test_second_order_flat_bed_onestep(self):
2241
2242	from mesh_factory import rectangular
2243	from Numeric import array
2244
2245	#Create basic mesh
2246	points, vertices, boundary = rectangular(6, 6)
2247
2248	#Create shallow water domain
2249	domain = Domain(points, vertices, boundary)
2250	domain.smooth = False
2251	domain.default_order=2
2252	domain.beta_w = 0.9
2253	domain.beta_w_dry = 0.9
2254	domain.beta_uh = 0.9
2255	domain.beta_uh_dry = 0.9
2256	domain.beta_vh = 0.9
2257	domain.beta_vh_dry = 0.9
2258
2259	# Boundary conditions
2260	Br = Reflective_boundary(domain)
2261	Bd = Dirichlet_boundary([0.1, 0., 0.])
2262	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2263
2264	domain.check_integrity()
2265
2266	#Evolution
2267	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2268	pass# domain.write_time()
2269
2270	#Data from earlier version of abstract_2d_finite_volumes
2271	assert allclose(domain.min_timestep, 0.0396825396825)
2272	assert allclose(domain.max_timestep, 0.0396825396825)
2273
2274	assert allclose(domain.quantities['stage'].centroid_values[:12],
2275	[0.00171396, 0.02656103, 0.00241523, 0.02656103,
2276	0.00241523, 0.02656103, 0.00241523, 0.02656103,
2277	0.00241523, 0.02656103, 0.00241523, 0.0272623])
2278
2279	domain.distribute_to_vertices_and_edges()
2280	assert allclose(domain.quantities['stage'].vertex_values[:12,0],
2281	[0.0001714, 0.02656103, 0.00024152,
2282	0.02656103, 0.00024152, 0.02656103,
2283	0.00024152, 0.02656103, 0.00024152,
2284	0.02656103, 0.00024152, 0.0272623])
2285
2286	assert allclose(domain.quantities['stage'].vertex_values[:12,1],
2287	[0.00315012, 0.02656103, 0.00024152, 0.02656103,
2288	0.00024152, 0.02656103, 0.00024152, 0.02656103,
2289	0.00024152, 0.02656103, 0.00040506, 0.0272623])
2290
2291	assert allclose(domain.quantities['stage'].vertex_values[:12,2],
2292	[0.00182037, 0.02656103, 0.00676264,
2293	0.02656103, 0.00676264, 0.02656103,
2294	0.00676264, 0.02656103, 0.00676264,
2295	0.02656103, 0.0065991, 0.0272623])
2296
2297	assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
2298	[0.00113961, 0.01302432, 0.00148672,
2299	0.01302432, 0.00148672, 0.01302432,
2300	0.00148672, 0.01302432, 0.00148672 ,
2301	0.01302432, 0.00148672, 0.01337143])
2302
2303	assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
2304	[-2.91240050e-004 , 1.22721531e-004,
2305	-1.22721531e-004, 1.22721531e-004 ,
2306	-1.22721531e-004, 1.22721531e-004,
2307	-1.22721531e-004 , 1.22721531e-004,
2308	-1.22721531e-004, 1.22721531e-004,
2309	-1.22721531e-004, -4.57969873e-005])
2310
2311	os.remove(domain.filename + '.sww')
2312
2313
2314	def test_second_order_flat_bed_moresteps(self):
2315
2316	from mesh_factory import rectangular
2317	from Numeric import array
2318
2319	#Create basic mesh
2320	points, vertices, boundary = rectangular(6, 6)
2321
2322	#Create shallow water domain
2323	domain = Domain(points, vertices, boundary)
2324	domain.smooth = False
2325	domain.default_order=2
2326
2327	# Boundary conditions
2328	Br = Reflective_boundary(domain)
2329	Bd = Dirichlet_boundary([0.1, 0., 0.])
2330	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2331
2332	domain.check_integrity()
2333
2334	#Evolution
2335	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2336	pass
2337
2338	#Data from earlier version of abstract_2d_finite_volumes
2339	#assert allclose(domain.min_timestep, 0.0396825396825)
2340	#assert allclose(domain.max_timestep, 0.0396825396825)
2341	#print domain.quantities['stage'].centroid_values
2342
2343	os.remove(domain.filename + '.sww')
2344
2345
2346	def test_flatbed_first_order(self):
2347	from mesh_factory import rectangular
2348	from Numeric import array
2349
2350	#Create basic mesh
2351	N = 8
2352	points, vertices, boundary = rectangular(N, N)
2353
2354	#Create shallow water domain
2355	domain = Domain(points, vertices, boundary)
2356	domain.smooth = False
2357	domain.visualise = False
2358	domain.default_order=1
2359
2360	# Boundary conditions
2361	Br = Reflective_boundary(domain)
2362	Bd = Dirichlet_boundary([0.2,0.,0.])
2363
2364	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2365	domain.check_integrity()
2366
2367
2368	#Evolution
2369	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
2370	pass
2371	#domain.write_time()
2372
2373	#FIXME: These numbers were from version before 25/10
2374	#assert allclose(domain.min_timestep, 0.0140413643926)
2375	#assert allclose(domain.max_timestep, 0.0140947355753)
2376
2377	for i in range(3):
2378	#assert allclose(domain.quantities['stage'].edge_values[:4,i],
2379	# [0.10730244,0.12337617,0.11200126,0.12605666])
2380
2381	assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
2382	[0.07610894,0.06901572,0.07284461,0.06819712])
2383
2384	#assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
2385	# [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
2386
2387
2388	os.remove(domain.filename + '.sww')
2389
2390	def test_flatbed_second_order(self):
2391	from mesh_factory import rectangular
2392	from Numeric import array
2393
2394	#Create basic mesh
2395	N = 8
2396	points, vertices, boundary = rectangular(N, N)
2397
2398	#Create shallow water domain
2399	domain = Domain(points, vertices, boundary)
2400	domain.smooth = False
2401	domain.visualise = False
2402	domain.default_order=2
2403	domain.beta_w = 0.9
2404	domain.beta_w_dry = 0.9
2405	domain.beta_uh = 0.9
2406	domain.beta_uh_dry = 0.9
2407	domain.beta_vh = 0.9
2408	domain.beta_vh_dry = 0.9
2409	#domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
2410
2411	# Boundary conditions
2412	Br = Reflective_boundary(domain)
2413	Bd = Dirichlet_boundary([0.2,0.,0.])
2414
2415	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2416	domain.check_integrity()
2417
2418	#Evolution
2419	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2420	pass
2421
2422
2423	assert allclose(domain.min_timestep, 0.0210448446782)
2424	assert allclose(domain.max_timestep, 0.0210448446782)
2425
2426	#print domain.quantities['stage'].vertex_values[:4,0]
2427	#print domain.quantities['xmomentum'].vertex_values[:4,0]
2428	#print domain.quantities['ymomentum'].vertex_values[:4,0]
2429
2430	#FIXME: These numbers were from version before 25/10
2431	#assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2432	# [0.00101913,0.05352143,0.00104852,0.05354394])
2433
2434	#FIXME: These numbers were from version before 21/3/6 -
2435	#could be recreated by setting maximum_allowed_speed to 0 maybe
2436	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2437	# [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
2438
2439	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2440	[ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
2441
2442
2443
2444	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2445	# [0.00090581,0.03685719,0.00088303,0.03687313])
2446
2447	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2448	[-0.00139463,0.0006156,-0.00060364,0.00061827])
2449
2450
2451	os.remove(domain.filename + '.sww')
2452
2453
2454	def test_flatbed_second_order_vmax_0(self):
2455	from mesh_factory import rectangular
2456	from Numeric import array
2457
2458	#Create basic mesh
2459	N = 8
2460	points, vertices, boundary = rectangular(N, N)
2461
2462	#Create shallow water domain
2463	domain = Domain(points, vertices, boundary)
2464	domain.smooth = False
2465	domain.visualise = False
2466	domain.default_order=2
2467	domain.beta_w = 0.9
2468	domain.beta_w_dry = 0.9
2469	domain.beta_uh = 0.9
2470	domain.beta_uh_dry = 0.9
2471	domain.beta_vh = 0.9
2472	domain.beta_vh_dry = 0.9
2473	domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
2474
2475	# Boundary conditions
2476	Br = Reflective_boundary(domain)
2477	Bd = Dirichlet_boundary([0.2,0.,0.])
2478
2479	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2480	domain.check_integrity()
2481
2482	#Evolution
2483	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2484	pass
2485
2486
2487	assert allclose(domain.min_timestep, 0.0210448446782)
2488	assert allclose(domain.max_timestep, 0.0210448446782)
2489
2490	#FIXME: These numbers were from version before 21/3/6 -
2491	#could be recreated by setting maximum_allowed_speed to 0 maybe
2492	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2493	[ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
2494
2495
2496	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2497	[-0.00139463,0.0006156,-0.00060364,0.00061827])
2498
2499
2500	os.remove(domain.filename + '.sww')
2501
2502
2503
2504	def test_flatbed_second_order_distribute(self):
2505	#Use real data from anuga.abstract_2d_finite_volumes 2
2506	#painfully setup and extracted.
2507	from mesh_factory import rectangular
2508	from Numeric import array
2509
2510	#Create basic mesh
2511	N = 8
2512	points, vertices, boundary = rectangular(N, N)
2513
2514	#Create shallow water domain
2515	domain = Domain(points, vertices, boundary)
2516	domain.smooth = False
2517	domain.visualise = False
2518	domain.default_order=domain._order_=2
2519	domain.beta_w = 0.9
2520	domain.beta_w_dry = 0.9
2521	domain.beta_uh = 0.9
2522	domain.beta_uh_dry = 0.9
2523	domain.beta_vh = 0.9
2524	domain.beta_vh_dry = 0.9
2525
2526	# Boundary conditions
2527	Br = Reflective_boundary(domain)
2528	Bd = Dirichlet_boundary([0.2,0.,0.])
2529
2530	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2531	domain.check_integrity()
2532
2533
2534
2535	for V in [False, True]:
2536	if V:
2537	#Set centroids as if system had been evolved
2538	L = zeros(2NN, Float)
2539	L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
2540	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2541	1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
2542	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2543	1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
2544	5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
2545	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2546	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2547	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2548	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2549	0.00000000e+000, 5.57305948e-005]
2550
2551	X = zeros(2NN, Float)
2552	X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
2553	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2554	8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
2555	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2556	8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
2557	3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
2558	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2559	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2560	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2561	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2562	0.00000000e+000, 4.57662812e-005]
2563
2564	Y = zeros(2NN, Float)
2565	Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
2566	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2567	-6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
2568	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2569	-6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
2570	-1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
2571	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2572	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2573	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2574	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2575	0.00000000e+000, -2.57635178e-005]
2576
2577
2578	domain.set_quantity('stage', L, location='centroids')
2579	domain.set_quantity('xmomentum', X, location='centroids')
2580	domain.set_quantity('ymomentum', Y, location='centroids')
2581
2582	domain.check_integrity()
2583	else:
2584	#Evolution
2585	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2586	pass
2587	assert allclose(domain.min_timestep, 0.0210448446782)
2588	assert allclose(domain.max_timestep, 0.0210448446782)
2589
2590
2591	#Centroids were correct but not vertices.
2592	#Hence the check of distribute below.
2593	assert allclose(domain.quantities['stage'].centroid_values[:4],
2594	[0.00721206,0.05352143,0.01009108,0.05354394])
2595
2596	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
2597	[0.00648352,0.03685719,0.00850733,0.03687313])
2598
2599	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2600	[-0.00139463,0.0006156,-0.00060364,0.00061827])
2601
2602	#print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
2603	#print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
2604
2605	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
2606	##print domain.quantities['xmomentum'].centroid_values[17], V
2607	##print
2608	if not V:
2609	#FIXME: These numbers were from version before 21/3/6 -
2610	#could be recreated by setting maximum_allowed_speed to 0 maybe
2611
2612	#assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
2613	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)
2614
2615	else:
2616	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
2617
2618	import copy
2619	XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
2620	assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2621
2622	domain.distribute_to_vertices_and_edges()
2623
2624	#assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2625
2626	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],
2627	# 0.0)
2628	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)
2629
2630
2631	#FIXME: These numbers were from version before 25/10
2632	#assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2633	# [0.00101913,0.05352143,0.00104852,0.05354394])
2634
2635	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2636	[-0.00139463,0.0006156,-0.00060364,0.00061827])
2637
2638
2639	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2640	[0.00090581,0.03685719,0.00088303,0.03687313])
2641
2642
2643	#NB NO longer relvant:
2644
2645	#This was the culprit. First triangles vertex 0 had an
2646	#x-momentum of 0.0064835 instead of 0.00090581 and
2647	#third triangle had 0.00850733 instead of 0.00088303
2648	#print domain.quantities['xmomentum'].vertex_values[:4,0]
2649
2650	#print domain.quantities['xmomentum'].vertex_values[:4,0]
2651	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2652	# [0.00090581,0.03685719,0.00088303,0.03687313])
2653
2654	os.remove(domain.filename + '.sww')
2655
2656
2657
2658	def test_bedslope_problem_first_order(self):
2659
2660	from mesh_factory import rectangular
2661	from Numeric import array
2662
2663	#Create basic mesh
2664	points, vertices, boundary = rectangular(6, 6)
2665
2666	#Create shallow water domain
2667	domain = Domain(points, vertices, boundary)
2668	domain.smooth = False
2669	domain.default_order=1
2670
2671	#Bed-slope and friction
2672	def x_slope(x, y):
2673	return -x/3
2674
2675	domain.set_quantity('elevation', x_slope)
2676
2677	# Boundary conditions
2678	Br = Reflective_boundary(domain)
2679	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2680
2681	#Initial condition
2682	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2683	domain.check_integrity()
2684
2685	#Evolution
2686	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2687	pass# domain.write_time()
2688
2689	assert allclose(domain.min_timestep, 0.050010003001)
2690	assert allclose(domain.max_timestep, 0.050010003001)
2691
2692
2693	os.remove(domain.filename + '.sww')
2694
2695	def test_bedslope_problem_first_order_moresteps(self):
2696
2697	from mesh_factory import rectangular
2698	from Numeric import array
2699
2700	#Create basic mesh
2701	points, vertices, boundary = rectangular(6, 6)
2702
2703	#Create shallow water domain
2704	domain = Domain(points, vertices, boundary)
2705	domain.smooth = False
2706	domain.default_order=1
2707	domain.beta_h = 0.0 #Use first order in h-limiter
2708
2709	#Bed-slope and friction
2710	def x_slope(x, y):
2711	return -x/3
2712
2713	domain.set_quantity('elevation', x_slope)
2714
2715	# Boundary conditions
2716	Br = Reflective_boundary(domain)
2717	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2718
2719	#Initial condition
2720	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2721	domain.check_integrity()
2722
2723	#Evolution
2724	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2725	pass# domain.write_time()
2726
2727	#Data from earlier version of abstract_2d_finite_volumes
2728	#print domain.quantities['stage'].centroid_values
2729
2730	assert allclose(domain.quantities['stage'].centroid_values,
2731	[-0.02998628, -0.01520652, -0.03043492,
2732	-0.0149132, -0.03004706, -0.01476251,
2733	-0.0298215, -0.01467976, -0.02988158,
2734	-0.01474662, -0.03036161, -0.01442995,
2735	-0.07624583, -0.06297061, -0.07733792,
2736	-0.06342237, -0.07695439, -0.06289595,
2737	-0.07635559, -0.0626065, -0.07633628,
2738	-0.06280072, -0.07739632, -0.06386738,
2739	-0.12161738, -0.11028239, -0.1223796,
2740	-0.11095953, -0.12189744, -0.11048616,
2741	-0.12074535, -0.10987605, -0.12014311,
2742	-0.10976691, -0.12096859, -0.11087692,
2743	-0.16868259, -0.15868061, -0.16801135,
2744	-0.1588003, -0.16674343, -0.15813323,
2745	-0.16457595, -0.15693826, -0.16281096,
2746	-0.15585154, -0.16283873, -0.15540068,
2747	-0.17450362, -0.19919913, -0.18062882,
2748	-0.19764131, -0.17783111, -0.19407213,
2749	-0.1736915, -0.19053624, -0.17228678,
2750	-0.19105634, -0.17920133, -0.1968828,
2751	-0.14244395, -0.14604641, -0.14473537,
2752	-0.1506107, -0.14510055, -0.14919522,
2753	-0.14175896, -0.14560798, -0.13911658,
2754	-0.14439383, -0.13924047, -0.14829043])
2755
2756	os.remove(domain.filename + '.sww')
2757
2758	def test_bedslope_problem_second_order_one_step(self):
2759
2760	from mesh_factory import rectangular
2761	from Numeric import array
2762
2763	#Create basic mesh
2764	points, vertices, boundary = rectangular(6, 6)
2765
2766	#Create shallow water domain
2767	domain = Domain(points, vertices, boundary)
2768	domain.smooth = False
2769	domain.default_order=2
2770	domain.beta_w = 0.9
2771	domain.beta_w_dry = 0.9
2772	domain.beta_uh = 0.9
2773	domain.beta_uh_dry = 0.9
2774	domain.beta_vh = 0.9
2775	domain.beta_vh_dry = 0.9
2776
2777	#Bed-slope and friction at vertices (and interpolated elsewhere)
2778	def x_slope(x, y):
2779	return -x/3
2780
2781	domain.set_quantity('elevation', x_slope)
2782
2783	# Boundary conditions
2784	Br = Reflective_boundary(domain)
2785	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2786
2787	#Initial condition
2788	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2789	domain.check_integrity()
2790
2791	assert allclose(domain.quantities['stage'].centroid_values,
2792	[0.01296296, 0.03148148, 0.01296296,
2793	0.03148148, 0.01296296, 0.03148148,
2794	0.01296296, 0.03148148, 0.01296296,
2795	0.03148148, 0.01296296, 0.03148148,
2796	-0.04259259, -0.02407407, -0.04259259,
2797	-0.02407407, -0.04259259, -0.02407407,
2798	-0.04259259, -0.02407407, -0.04259259,
2799	-0.02407407, -0.04259259, -0.02407407,
2800	-0.09814815, -0.07962963, -0.09814815,
2801	-0.07962963, -0.09814815, -0.07962963,
2802	-0.09814815, -0.07962963, -0.09814815,
2803	-0.07962963, -0.09814815, -0.07962963,
2804	-0.1537037 , -0.13518519, -0.1537037,
2805	-0.13518519, -0.1537037, -0.13518519,
2806	-0.1537037 , -0.13518519, -0.1537037,
2807	-0.13518519, -0.1537037, -0.13518519,
2808	-0.20925926, -0.19074074, -0.20925926,
2809	-0.19074074, -0.20925926, -0.19074074,
2810	-0.20925926, -0.19074074, -0.20925926,
2811	-0.19074074, -0.20925926, -0.19074074,
2812	-0.26481481, -0.2462963, -0.26481481,
2813	-0.2462963, -0.26481481, -0.2462963,
2814	-0.26481481, -0.2462963, -0.26481481,
2815	-0.2462963, -0.26481481, -0.2462963])
2816
2817
2818	#print domain.quantities['stage'].extrapolate_second_order()
2819	#domain.distribute_to_vertices_and_edges()
2820	#print domain.quantities['stage'].vertex_values[:,0]
2821
2822	#Evolution
2823	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2824	#domain.write_time()
2825	pass
2826
2827
2828	#print domain.quantities['stage'].centroid_values
2829	assert allclose(domain.quantities['stage'].centroid_values,
2830	[0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
2831	0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
2832	0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
2833	-0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
2834	-0.04186556, -0.0248011, -0.04186556, -0.02255593,
2835	-0.09966629, -0.08035666, -0.09742112, -0.08035666,
2836	-0.09742112, -0.08035666, -0.09742112, -0.08035666,
2837	-0.09742112, -0.08035666, -0.09742112, -0.07811149,
2838	-0.15522185, -0.13591222, -0.15297667, -0.13591222,
2839	-0.15297667, -0.13591222, -0.15297667, -0.13591222,
2840	-0.15297667, -0.13591222, -0.15297667, -0.13366704,
2841	-0.2107774, -0.19146777, -0.20853223, -0.19146777,
2842	-0.20853223, -0.19146777, -0.20853223, -0.19146777,
2843	-0.20853223, -0.19146777, -0.20853223, -0.1892226,
2844	-0.26120669, -0.24776246, -0.25865535, -0.24776246,
2845	-0.25865535, -0.24776246, -0.25865535, -0.24776246,
2846	-0.25865535, -0.24776246, -0.25865535, -0.24521113])
2847
2848	os.remove(domain.filename + '.sww')
2849
2850	def test_bedslope_problem_second_order_two_steps(self):
2851
2852	from mesh_factory import rectangular
2853	from Numeric import array
2854
2855	#Create basic mesh
2856	points, vertices, boundary = rectangular(6, 6)
2857
2858	#Create shallow water domain
2859	domain = Domain(points, vertices, boundary)
2860	domain.smooth = False
2861	domain.default_order=2
2862	domain.beta_w = 0.9
2863	domain.beta_w_dry = 0.9
2864	domain.beta_uh = 0.9
2865	domain.beta_uh_dry = 0.9
2866	domain.beta_vh = 0.9
2867	domain.beta_vh_dry = 0.9
2868	domain.beta_h = 0.0 #Use first order in h-limiter
2869
2870	#Bed-slope and friction at vertices (and interpolated elsewhere)
2871	def x_slope(x, y):
2872	return -x/3
2873
2874	domain.set_quantity('elevation', x_slope)
2875
2876	# Boundary conditions
2877	Br = Reflective_boundary(domain)
2878	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2879
2880	#Initial condition
2881	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2882	domain.check_integrity()
2883
2884	assert allclose(domain.quantities['stage'].centroid_values,
2885	[0.01296296, 0.03148148, 0.01296296,
2886	0.03148148, 0.01296296, 0.03148148,
2887	0.01296296, 0.03148148, 0.01296296,
2888	0.03148148, 0.01296296, 0.03148148,
2889	-0.04259259, -0.02407407, -0.04259259,
2890	-0.02407407, -0.04259259, -0.02407407,
2891	-0.04259259, -0.02407407, -0.04259259,
2892	-0.02407407, -0.04259259, -0.02407407,
2893	-0.09814815, -0.07962963, -0.09814815,
2894	-0.07962963, -0.09814815, -0.07962963,
2895	-0.09814815, -0.07962963, -0.09814815,
2896	-0.07962963, -0.09814815, -0.07962963,
2897	-0.1537037 , -0.13518519, -0.1537037,
2898	-0.13518519, -0.1537037, -0.13518519,
2899	-0.1537037 , -0.13518519, -0.1537037,
2900	-0.13518519, -0.1537037, -0.13518519,
2901	-0.20925926, -0.19074074, -0.20925926,
2902	-0.19074074, -0.20925926, -0.19074074,
2903	-0.20925926, -0.19074074, -0.20925926,
2904	-0.19074074, -0.20925926, -0.19074074,
2905	-0.26481481, -0.2462963, -0.26481481,
2906	-0.2462963, -0.26481481, -0.2462963,
2907	-0.26481481, -0.2462963, -0.26481481,
2908	-0.2462963, -0.26481481, -0.2462963])
2909
2910
2911	#print domain.quantities['stage'].extrapolate_second_order()
2912	#domain.distribute_to_vertices_and_edges()
2913	#print domain.quantities['stage'].vertex_values[:,0]
2914
2915	#Evolution
2916	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2917	pass
2918
2919
2920	#Data from earlier version of abstract_2d_finite_volumes ft=0.1
2921	assert allclose(domain.min_timestep, 0.0376895634803)
2922	assert allclose(domain.max_timestep, 0.0415635655309)
2923
2924
2925	assert allclose(domain.quantities['stage'].centroid_values,
2926	[0.00855788, 0.01575204, 0.00994606, 0.01717072,
2927	0.01005985, 0.01716362, 0.01005985, 0.01716299,
2928	0.01007098, 0.01736248, 0.01216452, 0.02026776,
2929	-0.04462374, -0.02479045, -0.04199789, -0.0229465,
2930	-0.04184033, -0.02295693, -0.04184013, -0.02295675,
2931	-0.04184486, -0.0228168, -0.04028876, -0.02036486,
2932	-0.10029444, -0.08170809, -0.09772846, -0.08021704,
2933	-0.09760006, -0.08022143, -0.09759984, -0.08022124,
2934	-0.09760261, -0.08008893, -0.09603914, -0.07758209,
2935	-0.15584152, -0.13723138, -0.15327266, -0.13572906,
2936	-0.15314427, -0.13573349, -0.15314405, -0.13573331,
2937	-0.15314679, -0.13560104, -0.15158523, -0.13310701,
2938	-0.21208605, -0.19283913, -0.20955631, -0.19134189,
2939	-0.20942821, -0.19134598, -0.20942799, -0.1913458,
2940	-0.20943005, -0.19120952, -0.20781177, -0.18869401,
2941	-0.25384082, -0.2463294, -0.25047649, -0.24464654,
2942	-0.25031159, -0.24464253, -0.25031112, -0.24464253,
2943	-0.25031463, -0.24454764, -0.24885323, -0.24286438])
2944
2945
2946	os.remove(domain.filename + '.sww')
2947
2948	def test_bedslope_problem_second_order_two_yieldsteps(self):
2949
2950	from mesh_factory import rectangular
2951	from Numeric import array
2952
2953	#Create basic mesh
2954	points, vertices, boundary = rectangular(6, 6)
2955
2956	#Create shallow water domain
2957	domain = Domain(points, vertices, boundary)
2958	domain.smooth = False
2959	domain.default_order=2
2960	domain.beta_w = 0.9
2961	domain.beta_w_dry = 0.9
2962	domain.beta_uh = 0.9
2963	domain.beta_uh_dry = 0.9
2964	domain.beta_vh = 0.9
2965	domain.beta_vh_dry = 0.9
2966	domain.beta_h = 0.0 #Use first order in h-limiter
2967
2968	#Bed-slope and friction at vertices (and interpolated elsewhere)
2969	def x_slope(x, y):
2970	return -x/3
2971
2972	domain.set_quantity('elevation', x_slope)
2973
2974	# Boundary conditions
2975	Br = Reflective_boundary(domain)
2976	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2977
2978	#Initial condition
2979	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2980	domain.check_integrity()
2981
2982	assert allclose(domain.quantities['stage'].centroid_values,
2983	[0.01296296, 0.03148148, 0.01296296,
2984	0.03148148, 0.01296296, 0.03148148,
2985	0.01296296, 0.03148148, 0.01296296,
2986	0.03148148, 0.01296296, 0.03148148,
2987	-0.04259259, -0.02407407, -0.04259259,
2988	-0.02407407, -0.04259259, -0.02407407,
2989	-0.04259259, -0.02407407, -0.04259259,
2990	-0.02407407, -0.04259259, -0.02407407,
2991	-0.09814815, -0.07962963, -0.09814815,
2992	-0.07962963, -0.09814815, -0.07962963,
2993	-0.09814815, -0.07962963, -0.09814815,
2994	-0.07962963, -0.09814815, -0.07962963,
2995	-0.1537037 , -0.13518519, -0.1537037,
2996	-0.13518519, -0.1537037, -0.13518519,
2997	-0.1537037 , -0.13518519, -0.1537037,
2998	-0.13518519, -0.1537037, -0.13518519,
2999	-0.20925926, -0.19074074, -0.20925926,
3000	-0.19074074, -0.20925926, -0.19074074,
3001	-0.20925926, -0.19074074, -0.20925926,
3002	-0.19074074, -0.20925926, -0.19074074,
3003	-0.26481481, -0.2462963, -0.26481481,
3004	-0.2462963, -0.26481481, -0.2462963,
3005	-0.26481481, -0.2462963, -0.26481481,
3006	-0.2462963, -0.26481481, -0.2462963])
3007
3008
3009	#print domain.quantities['stage'].extrapolate_second_order()
3010	#domain.distribute_to_vertices_and_edges()
3011	#print domain.quantities['stage'].vertex_values[:,0]
3012
3013	#Evolution
3014	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): #0.05??
3015	#domain.write_time()
3016	pass
3017
3018
3019
3020	assert allclose(domain.quantities['stage'].centroid_values,
3021	[0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
3022	0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
3023	0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
3024	-0.0229465, -0.04184033, -0.02295693, -0.04184013,
3025	-0.02295675, -0.04184486, -0.0228168, -0.04028876,
3026	-0.02036486, -0.10029444, -0.08170809, -0.09772846,
3027	-0.08021704, -0.09760006, -0.08022143, -0.09759984,
3028	-0.08022124, -0.09760261, -0.08008893, -0.09603914,
3029	-0.07758209, -0.15584152, -0.13723138, -0.15327266,
3030	-0.13572906, -0.15314427, -0.13573349, -0.15314405,
3031	-0.13573331, -0.15314679, -0.13560104, -0.15158523,
3032	-0.13310701, -0.21208605, -0.19283913, -0.20955631,
3033	-0.19134189, -0.20942821, -0.19134598, -0.20942799,
3034	-0.1913458, -0.20943005, -0.19120952, -0.20781177,
3035	-0.18869401, -0.25384082, -0.2463294, -0.25047649,
3036	-0.24464654, -0.25031159, -0.24464253, -0.25031112,
3037	-0.24464253, -0.25031463, -0.24454764, -0.24885323,
3038	-0.24286438])
3039
3040	os.remove(domain.filename + '.sww')
3041
3042	def test_bedslope_problem_second_order_more_steps(self):
3043
3044	from mesh_factory import rectangular
3045	from Numeric import array
3046
3047	#Create basic mesh
3048	points, vertices, boundary = rectangular(6, 6)
3049
3050	#Create shallow water domain
3051	domain = Domain(points, vertices, boundary)
3052	domain.smooth = False
3053	domain.default_order=2
3054	domain.beta_w = 0.9
3055	domain.beta_w_dry = 0.9
3056	domain.beta_uh = 0.9
3057	domain.beta_uh_dry = 0.9
3058	domain.beta_vh = 0.9
3059	domain.beta_vh_dry = 0.9
3060	domain.beta_h = 0.0 #Use first order in h-limiter
3061
3062	#Bed-slope and friction at vertices (and interpolated elsewhere)
3063	def x_slope(x, y):
3064	return -x/3
3065
3066	domain.set_quantity('elevation', x_slope)
3067
3068	# Boundary conditions
3069	Br = Reflective_boundary(domain)
3070	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3071
3072	#Initial condition
3073	domain.set_quantity('stage', expression = 'elevation + 0.05')
3074	domain.check_integrity()
3075
3076	assert allclose(domain.quantities['stage'].centroid_values,
3077	[0.01296296, 0.03148148, 0.01296296,
3078	0.03148148, 0.01296296, 0.03148148,
3079	0.01296296, 0.03148148, 0.01296296,
3080	0.03148148, 0.01296296, 0.03148148,
3081	-0.04259259, -0.02407407, -0.04259259,
3082	-0.02407407, -0.04259259, -0.02407407,
3083	-0.04259259, -0.02407407, -0.04259259,
3084	-0.02407407, -0.04259259, -0.02407407,
3085	-0.09814815, -0.07962963, -0.09814815,
3086	-0.07962963, -0.09814815, -0.07962963,
3087	-0.09814815, -0.07962963, -0.09814815,
3088	-0.07962963, -0.09814815, -0.07962963,
3089	-0.1537037 , -0.13518519, -0.1537037,
3090	-0.13518519, -0.1537037, -0.13518519,
3091	-0.1537037 , -0.13518519, -0.1537037,
3092	-0.13518519, -0.1537037, -0.13518519,
3093	-0.20925926, -0.19074074, -0.20925926,
3094	-0.19074074, -0.20925926, -0.19074074,
3095	-0.20925926, -0.19074074, -0.20925926,
3096	-0.19074074, -0.20925926, -0.19074074,
3097	-0.26481481, -0.2462963, -0.26481481,
3098	-0.2462963, -0.26481481, -0.2462963,
3099	-0.26481481, -0.2462963, -0.26481481,
3100	-0.2462963, -0.26481481, -0.2462963])
3101
3102
3103	#print domain.quantities['stage'].extrapolate_second_order()
3104	#domain.distribute_to_vertices_and_edges()
3105	#print domain.quantities['stage'].vertex_values[:,0]
3106
3107	#Evolution
3108	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
3109	pass
3110
3111
3112	assert allclose(domain.quantities['stage'].centroid_values,
3113	[-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
3114	-0.02901975, -0.0141361, -0.02898816, -0.01418135, -0.02961409, -0.01403487,
3115	-0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
3116	-0.07571145, -0.06235231, -0.0756817, -0.06245309, -0.07652292, -0.06289946,
3117	-0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934, -0.1107174,
3118	-0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
3119	-0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493, -0.15697919,
3120	-0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
3121	-0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
3122	-0.1833593, -0.19505871, -0.1818806, -0.19418042, -0.18586159, -0.19576946,
3123	-0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
3124	-0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
3125
3126	assert allclose(domain.quantities['xmomentum'].centroid_values,
3127	[ 0.00831121, 0.00317948, 0.00731797, 0.00334939, 0.00764717, 0.00348053,
3128	0.00788729, 0.00356522, 0.00780649, 0.00341919, 0.00693525, 0.00310375,
3129	0.02166196, 0.01421475, 0.02017737, 0.01316839, 0.02037015, 0.01368659,
3130	0.02106, 0.01399161, 0.02074514, 0.01354935, 0.01887407, 0.0123113,
3131	0.03775083, 0.02855197, 0.03689337, 0.02759782, 0.03732848, 0.02812072,
3132	0.03872545, 0.02913348, 0.03880939, 0.02803804, 0.03546499, 0.0260039,
3133	0.0632131, 0.04730634, 0.0576324, 0.04592336, 0.05790921, 0.04690514,
3134	0.05986467, 0.04871165, 0.06170068, 0.04811572, 0.05657041, 0.04416292,
3135	0.08489642, 0.07188097, 0.07835261, 0.06843406, 0.07986412, 0.0698247,
3136	0.08201071, 0.07216756, 0.08378418, 0.07273624, 0.080399, 0.06645841,
3137	0.01631548, 0.04691608, 0.0206632, 0.044409, 0.02115518, 0.04560305,
3138	0.02160608, 0.04663725, 0.02174734, 0.04795559, 0.02281427, 0.05667111])
3139
3140
3141	assert allclose(domain.quantities['ymomentum'].centroid_values,
3142	[ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
3143	-9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
3144	-1.99445058e-004, -3.28493467e-004, 6.68217349e-005, -8.42042805e-006,
3145	5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
3146	-8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
3147	-2.70160173e-004, -5.35485454e-004, 2.60780997e-004, 3.12145471e-005,
3148	5.16189608e-004, 1.07069062e-004, 9.29989252e-005, -3.71211119e-004,
3149	1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
3150	-4.74025708e-004, -6.94463009e-004, 6.15092843e-005, 2.22106820e-004,
3151	-6.29589294e-004, 2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
3152	-4.17914641e-004, 6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
3153	-1.67424308e-003, -7.39485066e-004, -1.59966988e-003, 5.68262838e-005,
3154	-1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
3155	-1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
3156	-2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
3157	-1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
3158	-6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
3159	-4.50964850e-003, -3.06319963e-003, 6.08950810e-004, -4.79537921e-004])
3160
3161	os.remove(domain.filename + '.sww')
3162
3163
3164	def test_temp_play(self):
3165
3166	from mesh_factory import rectangular
3167	from Numeric import array
3168
3169	#Create basic mesh
3170	points, vertices, boundary = rectangular(5, 5)
3171
3172	#Create shallow water domain
3173	domain = Domain(points, vertices, boundary)
3174	domain.smooth = False
3175	domain.default_order=2
3176	domain.beta_w = 0.9
3177	domain.beta_w_dry = 0.9
3178	domain.beta_uh = 0.9
3179	domain.beta_uh_dry = 0.9
3180	domain.beta_vh = 0.9
3181	domain.beta_vh_dry = 0.9
3182	domain.beta_h = 0.0 #Use first order in h-limiter
3183
3184	#Bed-slope and friction at vertices (and interpolated elsewhere)
3185	def x_slope(x, y):
3186	return -x/3
3187
3188	domain.set_quantity('elevation', x_slope)
3189
3190	# Boundary conditions
3191	Br = Reflective_boundary(domain)
3192	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3193
3194	#Initial condition
3195	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
3196	domain.check_integrity()
3197
3198	#Evolution
3199	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
3200	pass
3201
3202	assert allclose(domain.quantities['stage'].centroid_values[:4],
3203	[0.00206836, 0.01296714, 0.00363415, 0.01438924])
3204	#print domain.quantities['xmomentum'].centroid_values[:4]
3205	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
3206	[0.01360154, 0.00671133, 0.01264578, 0.00648503])
3207	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
3208	[-1.19201077e-003, -7.23647546e-004,
3209	-6.39083123e-005, 6.29815168e-005])
3210
3211	os.remove(domain.filename + '.sww')
3212
3213	def test_complex_bed(self):
3214	#No friction is tested here
3215
3216	from mesh_factory import rectangular
3217	from Numeric import array
3218
3219	N = 12
3220	points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
3221	origin=(-0.07, 0))
3222
3223
3224	domain = Domain(points, vertices, boundary)
3225	domain.smooth = False
3226	domain.visualise = False
3227	domain.default_order=2
3228
3229
3230	inflow_stage = 0.1
3231	Z = Weir(inflow_stage)
3232	domain.set_quantity('elevation', Z)
3233
3234	Br = Reflective_boundary(domain)
3235	Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
3236	domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
3237
3238	domain.set_quantity('stage', Constant_height(Z, 0.))
3239
3240	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
3241	pass
3242
3243
3244	#print domain.quantities['stage'].centroid_values
3245
3246	#FIXME: These numbers were from version before 25/10
3247	#assert allclose(domain.quantities['stage'].centroid_values,
3248	# [3.95822638e-002, 5.61022588e-002, 4.66437868e-002, 5.73081011e-002,
3249	# 4.72394613e-002, 5.74684939e-002, 4.74309483e-002, 5.77458084e-002,
3250	# 4.80628177e-002, 5.85656225e-002, 4.90498542e-002, 6.02609831e-002,
3251	# 1.18470315e-002, 1.75136443e-002, 1.18035266e-002, 2.15565695e-002,
3252	# 1.31620268e-002, 2.14351640e-002, 1.32351076e-002, 2.15450687e-002,
3253	# 1.36414028e-002, 2.24274619e-002, 1.51689511e-002, 2.21789655e-002,
3254	# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
3255	# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
3256	# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003, 2.73985295e-004,
3257	# 1.86166683e-001, 8.74070369e-002, 1.86166712e-001, 8.74035875e-002,
3258	# 6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
3259	# 6.11666725e-002, 8.73846774e-002, 1.86166697e-001, 8.74171550e-002,
3260	# -4.83333333e-002, 1.18333333e-001, -4.83333333e-002, 1.18333333e-001,
3261	# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
3262	# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002, 1.18333333e-001,
3263	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
3264	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
3265	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
3266	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
3267	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
3268	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
3269	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
3270	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
3271	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
3272	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
3273	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
3274	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
3275	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
3276	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
3277	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
3278	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
3279	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
3280	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
3281	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
3282	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
3283	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
3284
3285	os.remove(domain.filename + '.sww')
3286
3287	def test_spatio_temporal_boundary_1(self):
3288	"""Test that boundary values can be read from file and interpolated
3289	in both time and space.
3290
3291	Verify that the same steady state solution is arrived at and that
3292	time interpolation works.
3293
3294	The full solution history is not exactly the same as
3295	file boundary must read and interpolate from smoothed version
3296	as stored in sww.
3297	"""
3298	import time
3299
3300	#Create sww file of simple propagation from left to right
3301	#through rectangular domain
3302
3303	from mesh_factory import rectangular
3304
3305	#Create basic mesh
3306	points, vertices, boundary = rectangular(3, 3)
3307
3308	#Create shallow water domain
3309	domain1 = Domain(points, vertices, boundary)
3310
3311	domain1.reduction = mean
3312	domain1.smooth = False #Exact result
3313
3314	domain1.default_order = 2
3315	domain1.store = True
3316	domain1.set_datadir('.')
3317	domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
3318
3319	#FIXME: This is extremely important!
3320	#How can we test if they weren't stored?
3321	domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
3322
3323
3324	#Bed-slope and friction at vertices (and interpolated elsewhere)
3325	domain1.set_quantity('elevation', 0)
3326	domain1.set_quantity('friction', 0)
3327
3328	# Boundary conditions
3329	Br = Reflective_boundary(domain1)
3330	Bd = Dirichlet_boundary([0.3,0,0])
3331	domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
3332	#Initial condition
3333	domain1.set_quantity('stage', 0)
3334	domain1.check_integrity()
3335
3336	finaltime = 5
3337	#Evolution (full domain - large steps)
3338	for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
3339	pass
3340	#domain1.write_time()
3341
3342	cv1 = domain1.quantities['stage'].centroid_values
3343
3344
3345	#Create a triangle shaped domain (reusing coordinates from domain 1),
3346	#formed from the lower and right hand boundaries and
3347	#the sw-ne diagonal
3348	#from domain 1. Call it domain2
3349
3350	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3351	[2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3352	[1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3353
3354	vertices = [ [1,2,0], [3,4,1], [2,1,4], [4,5,2],
3355	[6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
3356
3357	boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
3358	(4,2):'right', (6,2):'right', (8,2):'right',
3359	(0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
3360
3361	domain2 = Domain(points, vertices, boundary)
3362
3363	domain2.reduction = domain1.reduction
3364	domain2.smooth = False
3365	domain2.default_order = 2
3366
3367	#Bed-slope and friction at vertices (and interpolated elsewhere)
3368	domain2.set_quantity('elevation', 0)
3369	domain2.set_quantity('friction', 0)
3370	domain2.set_quantity('stage', 0)
3371
3372	# Boundary conditions
3373	Br = Reflective_boundary(domain2)
3374	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
3375	domain2)
3376	domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
3377	domain2.check_integrity()
3378
3379
3380
3381	#Evolution (small steps)
3382	for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
3383	pass
3384
3385
3386	#Use output from domain1 as spatio-temporal boundary for domain2
3387	#and verify that results at right hand side are close.
3388
3389	cv2 = domain2.quantities['stage'].centroid_values
3390
3391	#print take(cv1, (12,14,16)) #Right
3392	#print take(cv2, (4,6,8))
3393	#print take(cv1, (0,6,12)) #Bottom
3394	#print take(cv2, (0,1,4))
3395	#print take(cv1, (0,8,16)) #Diag
3396	#print take(cv2, (0,3,8))
3397
3398	assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
3399	assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
3400	assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
3401
3402	#Cleanup
3403	os.remove(domain1.filename + '.' + domain1.format)
3404	os.remove(domain2.filename + '.' + domain2.format)
3405
3406
3407
3408	def test_spatio_temporal_boundary_2(self):
3409	"""Test that boundary values can be read from file and interpolated
3410	in both time and space.
3411	This is a more basic test, verifying that boundary object
3412	produces the expected results
3413
3414
3415	"""
3416	import time
3417
3418	#Create sww file of simple propagation from left to right
3419	#through rectangular domain
3420
3421	from mesh_factory import rectangular
3422
3423	#Create basic mesh
3424	points, vertices, boundary = rectangular(3, 3)
3425
3426	#Create shallow water domain
3427	domain1 = Domain(points, vertices, boundary)
3428
3429	domain1.reduction = mean
3430	domain1.smooth = True #To mimic MOST output
3431
3432	domain1.default_order = 2
3433	domain1.store = True
3434	domain1.set_datadir('.')
3435	domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
3436
3437	#FIXME: This is extremely important!
3438	#How can we test if they weren't stored?
3439	domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
3440
3441
3442	#Bed-slope and friction at vertices (and interpolated elsewhere)
3443	domain1.set_quantity('elevation', 0)
3444	domain1.set_quantity('friction', 0)
3445
3446	# Boundary conditions
3447	Br = Reflective_boundary(domain1)
3448	Bd = Dirichlet_boundary([0.3,0,0])
3449	domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
3450	#Initial condition
3451	domain1.set_quantity('stage', 0)
3452	domain1.check_integrity()
3453
3454	finaltime = 5
3455	#Evolution (full domain - large steps)
3456	for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
3457	pass
3458	#domain1.write_time()
3459
3460
3461	#Create an triangle shaped domain (coinciding with some
3462	#coordinates from domain 1),
3463	#formed from the lower and right hand boundaries and
3464	#the sw-ne diagonal
3465	#from domain 1. Call it domain2
3466
3467	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3468	[2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3469	[1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3470
3471	vertices = [ [1,2,0],
3472	[3,4,1], [2,1,4], [4,5,2],
3473	[6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
3474
3475	boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
3476	(4,2):'right', (6,2):'right', (8,2):'right',
3477	(0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
3478
3479	domain2 = Domain(points, vertices, boundary)
3480
3481	domain2.reduction = domain1.reduction
3482	domain2.smooth = False
3483	domain2.default_order = 2
3484
3485	#Bed-slope and friction at vertices (and interpolated elsewhere)
3486	domain2.set_quantity('elevation', 0)
3487	domain2.set_quantity('friction', 0)
3488	domain2.set_quantity('stage', 0)
3489
3490
3491	#Read results for specific timesteps t=1 and t=2
3492	from Scientific.IO.NetCDF import NetCDFFile
3493	fid = NetCDFFile(domain1.filename + '.' + domain1.format)
3494
3495	x = fid.variables['x'][:]
3496	y = fid.variables['y'][:]
3497	s1 = fid.variables['stage'][1,:]
3498	s2 = fid.variables['stage'][2,:]
3499	fid.close()
3500
3501	from Numeric import take, reshape, concatenate
3502	shp = (len(x), 1)
3503	points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
3504	#The diagonal points of domain 1 are 0, 5, 10, 15
3505
3506	#print points[0], points[5], points[10], points[15]
3507	assert allclose( take(points, [0,5,10,15]),
3508	[[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
3509
3510
3511	# Boundary conditions
3512	Br = Reflective_boundary(domain2)
3513	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
3514	domain2)
3515	domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
3516	domain2.check_integrity()
3517
3518	#Test that interpolation points are the mid points of the all boundary
3519	#segments
3520
3521	boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
3522	[1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
3523	[1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
3524
3525	boundary_midpoints.sort()
3526	R = Bf.F.interpolation_points.tolist()
3527	R.sort()
3528	assert allclose(boundary_midpoints, R)
3529
3530	#Check spatially interpolated output at time == 1
3531	domain2.time = 1
3532
3533	#First diagonal midpoint
3534	R0 = Bf.evaluate(0,0)
3535	assert allclose(R0[0], (s1[0] + s1[5])/2)
3536
3537	#Second diagonal midpoint
3538	R0 = Bf.evaluate(3,0)
3539	assert allclose(R0[0], (s1[5] + s1[10])/2)
3540
3541	#First diagonal midpoint
3542	R0 = Bf.evaluate(8,0)
3543	assert allclose(R0[0], (s1[10] + s1[15])/2)
3544
3545	#Check spatially interpolated output at time == 2
3546	domain2.time = 2
3547
3548	#First diagonal midpoint
3549	R0 = Bf.evaluate(0,0)
3550	assert allclose(R0[0], (s2[0] + s2[5])/2)
3551
3552	#Second diagonal midpoint
3553	R0 = Bf.evaluate(3,0)
3554	assert allclose(R0[0], (s2[5] + s2[10])/2)
3555
3556	#First diagonal midpoint
3557	R0 = Bf.evaluate(8,0)
3558	assert allclose(R0[0], (s2[10] + s2[15])/2)
3559
3560
3561	#Now check temporal interpolation
3562
3563	domain2.time = 1 + 2.0/3
3564
3565	#First diagonal midpoint
3566	R0 = Bf.evaluate(0,0)
3567	assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
3568
3569	#Second diagonal midpoint
3570	R0 = Bf.evaluate(3,0)
3571	assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
3572
3573	#First diagonal midpoint
3574	R0 = Bf.evaluate(8,0)
3575	assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
3576
3577
3578
3579	#Cleanup
3580	os.remove(domain1.filename + '.' + domain1.format)
3581
3582
3583	def test_pmesh2Domain(self):
3584	import os
3585	import tempfile
3586
3587	fileName = tempfile.mktemp(".tsh")
3588	file = open(fileName,"w")
3589	file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
3590	0 0.0 0.0 0.0 0.0 0.01 \n \
3591	1 1.0 0.0 10.0 10.0 0.02 \n \
3592	2 0.0 1.0 0.0 10.0 0.03 \n \
3593	3 0.5 0.25 8.0 12.0 0.04 \n \
3594	# Vert att title \n \
3595	elevation \n \
3596	stage \n \
3597	friction \n \
3598	2 # <triangle #> [<vertex #>] [<neigbouring triangle #>] \n\
3599	0 0 3 2 -1 -1 1 dsg\n\
3600	1 0 1 3 -1 0 -1 ole nielsen\n\
3601	4 # <segment #> <vertex #> <vertex #> [boundary tag] \n\
3602	0 1 0 2 \n\
3603	1 0 2 3 \n\
3604	2 2 3 \n\
3605	3 3 1 1 \n\
3606	3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
3607	0 216.0 -86.0 \n \
3608	1 160.0 -167.0 \n \
3609	2 114.0 -91.0 \n \
3610	3 # <vertex #> <vertex #> [boundary tag] ...Mesh Segments... \n \
3611	0 0 1 0 \n \
3612	1 1 2 0 \n \
3613	2 2 0 0 \n \
3614	0 # <x> <y> ...Mesh Holes... \n \
3615	0 # <x> <y> <attribute>...Mesh Regions... \n \
3616	0 # <x> <y> <attribute>...Mesh Regions, area... \n\
3617	#Geo reference \n \
3618	56 \n \
3619	140 \n \
3620	120 \n")
3621	file.close()
3622
3623	tags = {}
3624	b1 = Dirichlet_boundary(conserved_quantities = array([0.0]))
3625	b2 = Dirichlet_boundary(conserved_quantities = array([1.0]))
3626	b3 = Dirichlet_boundary(conserved_quantities = array([2.0]))
3627	tags["1"] = b1
3628	tags["2"] = b2
3629	tags["3"] = b3
3630
3631	#from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
3632	#domain = pmesh_to_domain_instance(fileName, Domain)
3633
3634	domain = Domain(mesh_filename=fileName)
3635	#verbose=True, use_cache=True)
3636
3637	#print "domain.tagged_elements", domain.tagged_elements
3638	## check the quantities
3639	#print domain.quantities['elevation'].vertex_values
3640	answer = [[0., 8., 0.],
3641	[0., 10., 8.]]
3642	assert allclose(domain.quantities['elevation'].vertex_values,
3643	answer)
3644
3645	#print domain.quantities['stage'].vertex_values
3646	answer = [[0., 12., 10.],
3647	[0., 10., 12.]]
3648	assert allclose(domain.quantities['stage'].vertex_values,
3649	answer)
3650
3651	#print domain.quantities['friction'].vertex_values
3652	answer = [[0.01, 0.04, 0.03],
3653	[0.01, 0.02, 0.04]]
3654	assert allclose(domain.quantities['friction'].vertex_values,
3655	answer)
3656
3657	#print domain.quantities['friction'].vertex_values
3658	assert allclose(domain.tagged_elements['dsg'][0],0)
3659	assert allclose(domain.tagged_elements['ole nielsen'][0],1)
3660
3661	self.failUnless( domain.boundary[(1, 0)] == '1',
3662	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3663	self.failUnless( domain.boundary[(1, 2)] == '2',
3664	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3665	self.failUnless( domain.boundary[(0, 1)] == '3',
3666	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3667	self.failUnless( domain.boundary[(0, 0)] == 'exterior',
3668	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3669	#print "domain.boundary",domain.boundary
3670	self.failUnless( len(domain.boundary) == 4,
3671	"test_pmesh2Domain Too many boundaries")
3672	#FIXME change to use get_xllcorner
3673	#print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
3674	self.failUnless(domain.geo_reference.xllcorner == 140.0,
3675	"bad geo_referece")
3676
3677
3678	#************
3679
3680
3681	domain = Domain(fileName)
3682
3683	#print "domain.tagged_elements", domain.tagged_elements
3684	## check the quantities
3685	#print domain.quantities['elevation'].vertex_values
3686	answer = [[0., 8., 0.],
3687	[0., 10., 8.]]
3688	assert allclose(domain.quantities['elevation'].vertex_values,
3689	answer)
3690
3691	#print domain.quantities['stage'].vertex_values
3692	answer = [[0., 12., 10.],
3693	[0., 10., 12.]]
3694	assert allclose(domain.quantities['stage'].vertex_values,
3695	answer)
3696
3697	#print domain.quantities['friction'].vertex_values
3698	answer = [[0.01, 0.04, 0.03],
3699	[0.01, 0.02, 0.04]]
3700	assert allclose(domain.quantities['friction'].vertex_values,
3701	answer)
3702
3703	#print domain.quantities['friction'].vertex_values
3704	assert allclose(domain.tagged_elements['dsg'][0],0)
3705	assert allclose(domain.tagged_elements['ole nielsen'][0],1)
3706
3707	self.failUnless( domain.boundary[(1, 0)] == '1',
3708	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3709	self.failUnless( domain.boundary[(1, 2)] == '2',
3710	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3711	self.failUnless( domain.boundary[(0, 1)] == '3',
3712	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3713	self.failUnless( domain.boundary[(0, 0)] == 'exterior',
3714	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3715	#print "domain.boundary",domain.boundary
3716	self.failUnless( len(domain.boundary) == 4,
3717	"test_pmesh2Domain Too many boundaries")
3718	#FIXME change to use get_xllcorner
3719	#print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
3720	self.failUnless(domain.geo_reference.xllcorner == 140.0,
3721	"bad geo_referece")
3722	#************
3723	os.remove(fileName)
3724
3725	#-------------------------------------------------------------
3726
3727	if __name__ == "__main__":
3728	suite = unittest.makeSuite(Test_Shallow_Water,'test')
3729	#suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_3')
3730	runner = unittest.TextTestRunner(verbosity=1)
3731	runner.run(suite)

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

Download in other formats: