Context Navigation

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

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

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

Download in other formats: