Context Navigation

test_shallow_water_domain.py @ 6258

Last change on this file since 6258 was 6258, checked in by ole, 14 years ago
Made second order accuracy the default and tested.
File size: 243.8 KB

Line
1	#!/usr/bin/env python
2
3	import unittest, os
4	from math import pi, sqrt
5	import tempfile
6
7	from anuga.config import g, epsilon
8	from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
9	import Numeric as num
10	from anuga.utilities.numerical_tools import mean
11	from anuga.utilities.polygon import is_inside_polygon
12	from anuga.coordinate_transforms.geo_reference import Geo_reference
13	from anuga.abstract_2d_finite_volumes.quantity import Quantity
14	from anuga.geospatial_data.geospatial_data import Geospatial_data
15	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
16	from shallow_water_domain import *
17
18	# Get gateway to C implementation of flux function for direct testing
19	from shallow_water_ext import flux_function_central as flux_function
20
21	# For test_fitting_using_shallow_water_domain example
22	def linear_function(point):
23	point = num.array(point)
24	return point[:,0]+point[:,1]
25
26	class Weir:
27	"""Set a bathymetry for weir with a hole and a downstream gutter
28	x,y are assumed to be in the unit square
29	"""
30
31	def __init__(self, stage):
32	self.inflow_stage = stage
33
34	def __call__(self, x, y):
35
36	N = len(x)
37	assert N == len(y)
38
39	z = num.zeros(N, num.Float)
40	for i in range(N):
41	z[i] = -x[i]/2 #General slope
42
43	#Flattish bit to the left
44	if x[i] < 0.3:
45	z[i] = -x[i]/10
46
47	#Weir
48	if x[i] >= 0.3 and x[i] < 0.4:
49	z[i] = -x[i]+0.9
50
51	#Dip
52	x0 = 0.6
53	#depth = -1.3
54	depth = -1.0
55	#plateaux = -0.9
56	plateaux = -0.6
57	if y[i] < 0.7:
58	if x[i] > x0 and x[i] < 0.9:
59	z[i] = depth
60
61	#RHS plateaux
62	if x[i] >= 0.9:
63	z[i] = plateaux
64
65
66	elif y[i] >= 0.7 and y[i] < 1.5:
67	#Restrict and deepen
68	if x[i] >= x0 and x[i] < 0.8:
69	z[i] = depth-(y[i]/3-0.3)
70	#z[i] = depth-y[i]/5
71	#z[i] = depth
72	elif x[i] >= 0.8:
73	#RHS plateaux
74	z[i] = plateaux
75
76	elif y[i] >= 1.5:
77	if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
78	#Widen up and stay at constant depth
79	z[i] = depth-1.5/5
80	elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
81	#RHS plateaux
82	z[i] = plateaux
83
84
85	#Hole in weir (slightly higher than inflow condition)
86	if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
87	z[i] = -x[i]+self.inflow_stage + 0.02
88
89	#Channel behind weir
90	x0 = 0.5
91	if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
92	z[i] = -x[i]+self.inflow_stage + 0.02
93
94	if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
95	#Flatten it out towards the end
96	z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
97
98	#Hole to the east
99	x0 = 1.1; y0 = 0.35
100	#if x[i] < -0.2 and y < 0.5:
101	if num.sqrt((2(x[i]-x0))2 + (2(y[i]-y0))**2) < 0.2:
102	z[i] = num.sqrt(((x[i]-x0))2 + ((y[i]-y0))2)-1.0
103
104	#Tiny channel draining hole
105	if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
106	z[i] = -0.9 #North south
107
108	if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
109	z[i] = -1.0 + (x[i]-0.9)/3 #East west
110
111
112
113	#Stuff not in use
114
115	#Upward slope at inlet to the north west
116	#if x[i] < 0.0: # and y[i] > 0.5:
117	# #z[i] = -y[i]+0.5 #-x[i]/2
118	# z[i] = x[i]/4 - y[i]**2 + 0.5
119
120	#Hole to the west
121	#x0 = -0.4; y0 = 0.35 # center
122	#if sqrt((2(x[i]-x0))2 + (2(y[i]-y0))**2) < 0.2:
123	# z[i] = sqrt(((x[i]-x0))2 + ((y[i]-y0))2)-0.2
124
125
126
127
128
129	return z/2
130
131	class Weir_simple:
132	"""Set a bathymetry for weir with a hole and a downstream gutter
133	x,y are assumed to be in the unit square
134	"""
135
136	def __init__(self, stage):
137	self.inflow_stage = stage
138
139	def __call__(self, x, y):
140
141	N = len(x)
142	assert N == len(y)
143
144	z = num.zeros(N, num.Float)
145	for i in range(N):
146	z[i] = -x[i] #General slope
147
148	#Flat bit to the left
149	if x[i] < 0.3:
150	z[i] = -x[i]/10 #General slope
151
152	#Weir
153	if x[i] > 0.3 and x[i] < 0.4:
154	z[i] = -x[i]+0.9
155
156	#Dip
157	if x[i] > 0.6 and x[i] < 0.9:
158	z[i] = -x[i]-0.5 #-y[i]/5
159
160	#Hole in weir (slightly higher than inflow condition)
161	if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
162	z[i] = -x[i]+self.inflow_stage + 0.05
163
164
165	return z/2
166
167
168
169
170	#Variable windfield implemented using functions
171	def speed(t,x,y):
172	"""Large speeds halfway between center and edges
173	Low speeds at center and edges
174	"""
175
176	from math import exp, cos, pi
177
178	x = num.array(x)
179	y = num.array(y)
180
181	N = len(x)
182	s = 0*x #New array
183
184	for k in range(N):
185
186	r = num.sqrt(x[k]2 + y[k]2)
187
188	factor = exp( -(r-0.15)**2 )
189
190	s[k] = 4000 * factor * (cos(t2pi/150) + 2)
191
192	return s
193
194
195	def scalar_func(t,x,y):
196	"""Function that returns a scalar.
197	Used to test error message when Numeric array is expected
198	"""
199
200	return 17.7
201
202
203	def angle(t,x,y):
204	"""Rotating field
205	"""
206	from math import atan, pi
207
208	x = num.array(x)
209	y = num.array(y)
210
211	N = len(x)
212	a = 0*x #New array
213
214	for k in range(N):
215	r = num.sqrt(x[k]2 + y[k]2)
216
217	angle = atan(y[k]/x[k])
218
219	if x[k] < 0:
220	angle+=pi #atan in ]-pi/2; pi/2[
221
222	#Take normal direction
223	angle -= pi/2
224
225	#Ensure positive radians
226	if angle < 0:
227	angle += 2*pi
228
229	a[k] = angle/pi*180
230
231	return a
232
233
234	class Test_Shallow_Water(unittest.TestCase):
235	def setUp(self):
236	pass
237
238	def tearDown(self):
239	pass
240
241	def test_rotate(self):
242	normal = num.array([0.0,-1.0])
243
244	q = num.array([1.0,2.0,3.0])
245
246	r = rotate(q, normal, direction = 1)
247	assert r[0] == 1
248	assert r[1] == -3
249	assert r[2] == 2
250
251	w = rotate(r, normal, direction = -1)
252	assert num.allclose(w, q)
253
254	#Check error check
255	try:
256	rotate(r, num.array([1,1,1], num.Int)) #array default#
257	except:
258	pass
259	else:
260	raise 'Should have raised an exception'
261
262
263	# Individual flux tests
264	def test_flux_zero_case(self):
265	ql = num.zeros( 3, num.Float )
266	qr = num.zeros( 3, num.Float )
267	normal = num.zeros( 2, num.Float )
268	edgeflux = num.zeros( 3, num.Float )
269	zl = zr = 0.
270	H0 = 0.0
271
272	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
273
274	assert num.allclose(edgeflux, [0,0,0])
275	assert max_speed == 0.
276
277	def test_flux_constants(self):
278	w = 2.0
279
280	normal = num.array([1.,0])
281	ql = num.array([w, 0, 0])
282	qr = num.array([w, 0, 0])
283	edgeflux = num.zeros(3, num.Float)
284	zl = zr = 0.
285	h = w - (zl+zr)/2
286	H0 = 0.0
287
288	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
289	assert num.allclose(edgeflux, [0., 0.5gh**2, 0.])
290	assert max_speed == num.sqrt(g*h)
291
292	#def test_flux_slope(self):
293	# #FIXME: TODO
294	# w = 2.0
295	#
296	# normal = array([1.,0])
297	# ql = array([w, 0, 0])
298	# qr = array([w, 0, 0])
299	# zl = zr = 0.
300	# h = w - (zl+zr)/2
301	#
302	# flux, max_speed = flux_function(normal, ql, qr, zl, zr)
303	#
304	# assert allclose(flux, [0., 0.5gh**2, 0.])
305	# assert max_speed == sqrt(g*h)
306
307
308	def test_flux1(self):
309	#Use data from previous version of abstract_2d_finite_volumes
310	normal = num.array([1.,0])
311	ql = num.array([-0.2, 2, 3])
312	qr = num.array([-0.2, 2, 3])
313	zl = zr = -0.5
314	edgeflux = num.zeros(3, num.Float)
315
316	H0 = 0.0
317
318	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
319
320	assert num.allclose(edgeflux, [2.,13.77433333, 20.])
321	assert num.allclose(max_speed, 8.38130948661)
322
323
324	def test_flux2(self):
325	#Use data from previous version of abstract_2d_finite_volumes
326	normal = num.array([0., -1.])
327	ql = num.array([-0.075, 2, 3])
328	qr = num.array([-0.075, 2, 3])
329	zl = zr = -0.375
330
331	edgeflux = num.zeros(3, num.Float)
332	H0 = 0.0
333	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
334
335	assert num.allclose(edgeflux, [-3.,-20.0, -30.441])
336	assert num.allclose(max_speed, 11.7146428199)
337
338	def test_flux3(self):
339	#Use data from previous version of abstract_2d_finite_volumes
340	normal = num.array([-sqrt(2)/2, sqrt(2)/2])
341	ql = num.array([-0.075, 2, 3])
342	qr = num.array([-0.075, 2, 3])
343	zl = zr = -0.375
344
345	edgeflux = num.zeros(3, num.Float)
346	H0 = 0.0
347	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
348
349	assert num.allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019])
350	assert num.allclose(max_speed, 4.0716654239)
351
352	def test_flux4(self):
353	#Use data from previous version of abstract_2d_finite_volumes
354	normal = num.array([-sqrt(2)/2, sqrt(2)/2])
355	ql = num.array([-0.34319278, 0.10254161, 0.07273855])
356	qr = num.array([-0.30683287, 0.1071986, 0.05930515])
357	zl = zr = -0.375
358
359	edgeflux = num.zeros(3, num.Float)
360	H0 = 0.0
361	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
362
363	assert num.allclose(edgeflux, [-0.04072676, -0.07096636, -0.01604364])
364	assert num.allclose(max_speed, 1.31414103233)
365
366	def test_flux_computation(self):
367	"""test_flux_computation - test flux calculation (actual C implementation)
368	This one tests the constant case where only the pressure term contributes to each edge and cancels out
369	once the total flux has been summed up.
370	"""
371
372	a = [0.0, 0.0]
373	b = [0.0, 2.0]
374	c = [2.0,0.0]
375	d = [0.0, 4.0]
376	e = [2.0, 2.0]
377	f = [4.0,0.0]
378
379	points = [a, b, c, d, e, f]
380	#bac, bce, ecf, dbe, daf, dae
381	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
382
383	domain = Domain(points, vertices)
384	domain.check_integrity()
385
386	# The constant case
387	domain.set_quantity('elevation', -1)
388	domain.set_quantity('stage', 1)
389
390	domain.compute_fluxes()
391	assert num.allclose(domain.get_quantity('stage').explicit_update[1], 0) # Central triangle
392
393
394	# The more general case
395	def surface(x,y):
396	return -x/2
397
398	domain.set_quantity('elevation', -10)
399	domain.set_quantity('stage', surface)
400	domain.set_quantity('xmomentum', 1)
401
402	domain.compute_fluxes()
403
404	#print domain.get_quantity('stage').explicit_update
405	# FIXME (Ole): TODO the general case
406	#assert allclose(domain.get_quantity('stage').explicit_update[1], ........??)
407
408
409
410	def test_sw_domain_simple(self):
411	a = [0.0, 0.0]
412	b = [0.0, 2.0]
413	c = [2.0,0.0]
414	d = [0.0, 4.0]
415	e = [2.0, 2.0]
416	f = [4.0,0.0]
417
418	points = [a, b, c, d, e, f]
419	#bac, bce, ecf, dbe, daf, dae
420	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
421
422
423	#from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain
424	#msg = 'The class %s is not a subclass of the generic domain class %s'\
425	# %(DomainClass, Domain)
426	#assert issubclass(DomainClass, Domain), msg
427
428	domain = Domain(points, vertices)
429	domain.check_integrity()
430
431	for name in ['stage', 'xmomentum', 'ymomentum',
432	'elevation', 'friction']:
433	assert domain.quantities.has_key(name)
434
435
436	assert domain.get_conserved_quantities(0, edge=1) == 0.
437
438
439	def test_boundary_conditions(self):
440
441	a = [0.0, 0.0]
442	b = [0.0, 2.0]
443	c = [2.0,0.0]
444	d = [0.0, 4.0]
445	e = [2.0, 2.0]
446	f = [4.0,0.0]
447
448	points = [a, b, c, d, e, f]
449	#bac, bce, ecf, dbe
450	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
451	boundary = { (0, 0): 'Third',
452	(0, 2): 'First',
453	(2, 0): 'Second',
454	(2, 1): 'Second',
455	(3, 1): 'Second',
456	(3, 2): 'Third'}
457
458
459	domain = Domain(points, vertices, boundary)
460	domain.check_integrity()
461
462
463	domain.set_quantity('stage', [[1,2,3], [5,5,5],
464	[0,0,9], [-6, 3, 3]])
465
466	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
467	[3,3,3], [4, 4, 4]])
468
469	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
470	[30,30,30], [40, 40, 40]])
471
472
473	D = Dirichlet_boundary([5,2,1])
474	T = Transmissive_boundary(domain)
475	R = Reflective_boundary(domain)
476	domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
477
478	domain.update_boundary()
479
480	#Stage
481	assert domain.quantities['stage'].boundary_values[0] == 2.5
482	assert domain.quantities['stage'].boundary_values[0] ==\
483	domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
484	assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
485	assert domain.quantities['stage'].boundary_values[2] ==\
486	domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
487	assert domain.quantities['stage'].boundary_values[3] ==\
488	domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
489	assert domain.quantities['stage'].boundary_values[4] ==\
490	domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
491	assert domain.quantities['stage'].boundary_values[5] ==\
492	domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
493
494	#Xmomentum
495	assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
496	assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
497	assert domain.quantities['xmomentum'].boundary_values[2] ==\
498	domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
499	assert domain.quantities['xmomentum'].boundary_values[3] ==\
500	domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
501	assert domain.quantities['xmomentum'].boundary_values[4] ==\
502	domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
503	assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective
504
505	#Ymomentum
506	assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
507	assert domain.quantities['ymomentum'].boundary_values[1] == 1. #Dirichlet
508	assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
509	assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
510	assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
511	assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
512
513
514	def test_boundary_conditionsII(self):
515
516	a = [0.0, 0.0]
517	b = [0.0, 2.0]
518	c = [2.0,0.0]
519	d = [0.0, 4.0]
520	e = [2.0, 2.0]
521	f = [4.0,0.0]
522
523	points = [a, b, c, d, e, f]
524	#bac, bce, ecf, dbe
525	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
526	boundary = { (0, 0): 'Third',
527	(0, 2): 'First',
528	(2, 0): 'Second',
529	(2, 1): 'Second',
530	(3, 1): 'Second',
531	(3, 2): 'Third',
532	(0, 1): 'Internal'}
533
534
535	domain = Domain(points, vertices, boundary)
536	domain.check_integrity()
537
538
539	domain.set_quantity('stage', [[1,2,3], [5,5,5],
540	[0,0,9], [-6, 3, 3]])
541
542	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
543	[3,3,3], [4, 4, 4]])
544
545	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
546	[30,30,30], [40, 40, 40]])
547
548
549	D = Dirichlet_boundary([5,2,1])
550	T = Transmissive_boundary(domain)
551	R = Reflective_boundary(domain)
552	domain.set_boundary( {'First': D, 'Second': T,
553	'Third': R, 'Internal': None})
554
555	domain.update_boundary()
556	domain.check_integrity()
557
558
559
560
561	def test_boundary_conditionsIII(self):
562	"""test_boundary_conditionsIII
563
564	Test Transmissive_stage_zero_momentum_boundary
565	"""
566
567	a = [0.0, 0.0]
568	b = [0.0, 2.0]
569	c = [2.0,0.0]
570	d = [0.0, 4.0]
571	e = [2.0, 2.0]
572	f = [4.0,0.0]
573
574	points = [a, b, c, d, e, f]
575	#bac, bce, ecf, dbe
576	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
577	boundary = { (0, 0): 'Third',
578	(0, 2): 'First',
579	(2, 0): 'Second',
580	(2, 1): 'Second',
581	(3, 1): 'Second',
582	(3, 2): 'Third'}
583
584
585	domain = Domain(points, vertices, boundary)
586	domain.check_integrity()
587
588
589	domain.set_quantity('stage', [[1,2,3], [5,5,5],
590	[0,0,9], [-6, 3, 3]])
591
592	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
593	[3,3,3], [4, 4, 4]])
594
595	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
596	[30,30,30], [40, 40, 40]])
597
598
599	D = Dirichlet_boundary([5,2,1])
600	T = Transmissive_stage_zero_momentum_boundary(domain)
601	R = Reflective_boundary(domain)
602	domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
603
604	domain.update_boundary()
605
606	# Stage
607	assert domain.quantities['stage'].boundary_values[0] == 2.5
608	assert domain.quantities['stage'].boundary_values[0] ==\
609	domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
610	assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
611	assert domain.quantities['stage'].boundary_values[2] ==\
612	domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
613	assert domain.quantities['stage'].boundary_values[3] ==\
614	domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
615	assert domain.quantities['stage'].boundary_values[4] ==\
616	domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
617	assert domain.quantities['stage'].boundary_values[5] ==\
618	domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
619
620	# Xmomentum
621	assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
622	assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
623	assert domain.quantities['xmomentum'].boundary_values[2] == 0.0
624	assert domain.quantities['xmomentum'].boundary_values[3] == 0.0
625	assert domain.quantities['xmomentum'].boundary_values[4] == 0.0
626	assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective
627
628	# Ymomentum
629	assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
630	assert domain.quantities['ymomentum'].boundary_values[1] == 1. #Dirichlet
631	assert domain.quantities['ymomentum'].boundary_values[2] == 0.0
632	assert domain.quantities['ymomentum'].boundary_values[3] == 0.0
633	assert domain.quantities['ymomentum'].boundary_values[4] == 0.0
634	assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
635
636
637
638
639	def test_boundary_condition_time(self):
640	"""test_boundary_condition_time
641
642	This tests that boundary conditions are evaluated
643	using the right time from domain.
644
645	"""
646
647	# Setup computational domain
648	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
649
650
651	#--------------------------------------------------------------
652	# Setup computational domain
653	#--------------------------------------------------------------
654	N = 5
655	points, vertices, boundary = rectangular_cross(N, N)
656	domain = Domain(points, vertices, boundary)
657
658	#--------------------------------------------------------------
659	# Setup initial conditions
660	#--------------------------------------------------------------
661	domain.set_quantity('elevation', 0.0)
662	domain.set_quantity('friction', 0.0)
663	domain.set_quantity('stage', 0.0)
664
665
666	#--------------------------------------------------------------
667	# Setup boundary conditions
668	#--------------------------------------------------------------
669	Bt = Time_boundary(domain=domain, # Time dependent boundary
670	f=lambda t: [t, 0.0, 0.0])
671
672	Br = Reflective_boundary(domain) # Reflective wall
673
674	domain.set_boundary({'left': Bt, 'right': Br, 'top': Br, 'bottom': Br})
675
676	for t in domain.evolve(yieldstep = 10, finaltime = 20.0):
677	q = Bt.evaluate()
678
679	# FIXME (Ole): This test would not have passed in
680	# changeset:5846.
681	msg = 'Time boundary not evaluated correctly'
682	assert num.allclose(t, q[0]), msg
683
684
685
686	def test_compute_fluxes0(self):
687	# Do a full triangle and check that fluxes cancel out for
688	# the constant stage case
689
690	a = [0.0, 0.0]
691	b = [0.0, 2.0]
692	c = [2.0,0.0]
693	d = [0.0, 4.0]
694	e = [2.0, 2.0]
695	f = [4.0,0.0]
696
697	points = [a, b, c, d, e, f]
698	#bac, bce, ecf, dbe
699	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
700
701	domain = Domain(points, vertices)
702	domain.set_quantity('stage', [[2,2,2], [2,2,2],
703	[2,2,2], [2,2,2]])
704	domain.check_integrity()
705
706	assert num.allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
707	assert num.allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
708
709	zl=zr=0. # Assume flat bed
710
711	edgeflux = num.zeros(3, num.Float)
712	edgeflux0 = num.zeros(3, num.Float)
713	edgeflux1 = num.zeros(3, num.Float)
714	edgeflux2 = num.zeros(3, num.Float)
715	H0 = 0.0
716
717	# Flux across right edge of volume 1
718	normal = domain.get_normal(1,0)
719	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
720	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
721	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)
722
723	# Check that flux seen from other triangles is inverse
724	tmp = qr; qr=ql; ql=tmp
725	normal = domain.get_normal(2,2)
726	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
727
728	assert num.allclose(edgeflux0 + edgeflux, 0.)
729
730	# Flux across upper edge of volume 1
731	normal = domain.get_normal(1,1)
732	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
733	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
734	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)
735
736	# Check that flux seen from other triangles is inverse
737	tmp = qr; qr=ql; ql=tmp
738	normal = domain.get_normal(3,0)
739	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
740
741	assert num.allclose(edgeflux1 + edgeflux, 0.)
742
743
744	# Flux across lower left hypotenuse of volume 1
745	normal = domain.get_normal(1,2)
746	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
747	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
748	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)
749
750	# Check that flux seen from other triangles is inverse
751	tmp = qr; qr=ql; ql=tmp
752	normal = domain.get_normal(0,1)
753	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
754	assert num.allclose(edgeflux2 + edgeflux, 0.)
755
756
757	# Scale by edgelengths, add up anc check that total flux is zero
758	e0 = domain.edgelengths[1, 0]
759	e1 = domain.edgelengths[1, 1]
760	e2 = domain.edgelengths[1, 2]
761
762	assert num.allclose(e0edgeflux0+e1edgeflux1+e2*edgeflux2, 0.)
763
764	# Now check that compute_flux yields zeros as well
765	domain.compute_fluxes()
766
767	for name in ['stage', 'xmomentum', 'ymomentum']:
768	#print name, domain.quantities[name].explicit_update
769	assert num.allclose(domain.quantities[name].explicit_update[1], 0)
770
771
772
773	def test_compute_fluxes1(self):
774	#Use values from previous version
775
776	a = [0.0, 0.0]
777	b = [0.0, 2.0]
778	c = [2.0,0.0]
779	d = [0.0, 4.0]
780	e = [2.0, 2.0]
781	f = [4.0,0.0]
782
783	points = [a, b, c, d, e, f]
784	#bac, bce, ecf, dbe
785	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
786
787	domain = Domain(points, vertices)
788	val0 = 2.+2.0/3
789	val1 = 4.+4.0/3
790	val2 = 8.+2.0/3
791	val3 = 2.+8.0/3
792
793	domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
794	[val2, val2, val2], [val3, val3, val3]])
795	domain.check_integrity()
796
797	zl=zr=0. #Assume flat bed
798
799	edgeflux = num.zeros(3, num.Float)
800	edgeflux0 = num.zeros(3, num.Float)
801	edgeflux1 = num.zeros(3, num.Float)
802	edgeflux2 = num.zeros(3, num.Float)
803	H0 = 0.0
804
805
806	# Flux across right edge of volume 1
807	normal = domain.get_normal(1,0) #Get normal 0 of triangle 1
808	assert num.allclose(normal, [1, 0])
809
810	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
811	assert num.allclose(ql, [val1, 0, 0])
812
813	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
814	assert num.allclose(qr, [val2, 0, 0])
815
816	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)
817
818	# Flux across edge in the east direction (as per normal vector)
819	assert num.allclose(edgeflux0, [-15.3598804, 253.71111111, 0.])
820	assert num.allclose(max_speed, 9.21592824046)
821
822
823	#Flux across edge in the west direction (opposite sign for xmomentum)
824	normal_opposite = domain.get_normal(2,2) #Get normal 2 of triangle 2
825	assert num.allclose(normal_opposite, [-1, 0])
826
827	max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
828	#flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr)
829	assert num.allclose(edgeflux, [-15.3598804, -253.71111111, 0.])
830
831
832	#Flux across upper edge of volume 1
833	normal = domain.get_normal(1,1)
834	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
835	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
836	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)
837
838	assert num.allclose(edgeflux1, [2.4098563, 0., 123.04444444])
839	assert num.allclose(max_speed, 7.22956891292)
840
841	#Flux across lower left hypotenuse of volume 1
842	normal = domain.get_normal(1,2)
843	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
844	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
845	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)
846
847	assert num.allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738])
848	assert num.allclose(max_speed, 7.22956891292)
849
850	#Scale, add up and check that compute_fluxes is correct for vol 1
851	e0 = domain.edgelengths[1, 0]
852	e1 = domain.edgelengths[1, 1]
853	e2 = domain.edgelengths[1, 2]
854
855	total_flux = -(e0edgeflux0+e1edgeflux1+e2*edgeflux2)/domain.areas[1]
856	assert num.allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
857
858
859	domain.compute_fluxes()
860
861	#assert num.allclose(total_flux, domain.explicit_update[1,:])
862	for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
863	assert num.allclose(total_flux[i],
864	domain.quantities[name].explicit_update[1])
865
866	#assert allclose(domain.explicit_update, [
867	# [0., -69.68888889, -69.68888889],
868	# [-0.68218178, -166.6, -35.93333333],
869	# [-111.77316251, 69.68888889, 0.],
870	# [-35.68522449, 0., 69.68888889]])
871
872	assert num.allclose(domain.quantities['stage'].explicit_update,
873	[0., -0.68218178, -111.77316251, -35.68522449])
874	assert num.allclose(domain.quantities['xmomentum'].explicit_update,
875	[-69.68888889, -166.6, 69.68888889, 0])
876	assert num.allclose(domain.quantities['ymomentum'].explicit_update,
877	[-69.68888889, -35.93333333, 0., 69.68888889])
878
879
880	#assert allclose(domain.quantities[name].explicit_update
881
882
883
884
885
886	def test_compute_fluxes2(self):
887	#Random values, incl momentum
888
889	a = [0.0, 0.0]
890	b = [0.0, 2.0]
891	c = [2.0,0.0]
892	d = [0.0, 4.0]
893	e = [2.0, 2.0]
894	f = [4.0,0.0]
895
896	points = [a, b, c, d, e, f]
897	#bac, bce, ecf, dbe
898	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
899
900	domain = Domain(points, vertices)
901	val0 = 2.+2.0/3
902	val1 = 4.+4.0/3
903	val2 = 8.+2.0/3
904	val3 = 2.+8.0/3
905
906	zl=zr=0 #Assume flat zero bed
907	edgeflux = num.zeros(3, num.Float)
908	edgeflux0 = num.zeros(3, num.Float)
909	edgeflux1 = num.zeros(3, num.Float)
910	edgeflux2 = num.zeros(3, num.Float)
911	H0 = 0.0
912
913
914	domain.set_quantity('elevation', zl*num.ones( (4,3) ))
915
916
917	domain.set_quantity('stage', [[val0, val0-1, val0-2],
918	[val1, val1+1, val1],
919	[val2, val2-2, val2],
920	[val3-0.5, val3, val3]])
921
922	domain.set_quantity('xmomentum',
923	[[1, 2, 3], [3, 4, 5],
924	[1, -1, 0], [0, -2, 2]])
925
926	domain.set_quantity('ymomentum',
927	[[1, -1, 0], [0, -3, 2],
928	[0, 1, 0], [-1, 2, 2]])
929
930
931	domain.check_integrity()
932
933
934
935	#Flux across right edge of volume 1
936	normal = domain.get_normal(1,0)
937	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
938	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
939	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)
940
941	#Flux across upper edge of volume 1
942	normal = domain.get_normal(1,1)
943	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
944	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
945	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)
946
947	#Flux across lower left hypotenuse of volume 1
948	normal = domain.get_normal(1,2)
949	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
950	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
951	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)
952
953	#Scale, add up and check that compute_fluxes is correct for vol 1
954	e0 = domain.edgelengths[1, 0]
955	e1 = domain.edgelengths[1, 1]
956	e2 = domain.edgelengths[1, 2]
957
958	total_flux = -(e0edgeflux0+e1edgeflux1+e2*edgeflux2)/domain.areas[1]
959
960
961	domain.compute_fluxes()
962	for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
963	assert num.allclose(total_flux[i],
964	domain.quantities[name].explicit_update[1])
965	#assert allclose(total_flux, domain.explicit_update[1,:])
966
967
968	# FIXME (Ole): Need test like this for fluxes in very shallow water.
969	def test_compute_fluxes3(self):
970	#Random values, incl momentum
971
972	a = [0.0, 0.0]
973	b = [0.0, 2.0]
974	c = [2.0,0.0]
975	d = [0.0, 4.0]
976	e = [2.0, 2.0]
977	f = [4.0,0.0]
978
979	points = [a, b, c, d, e, f]
980	#bac, bce, ecf, dbe
981	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
982
983	domain = Domain(points, vertices)
984	val0 = 2.+2.0/3
985	val1 = 4.+4.0/3
986	val2 = 8.+2.0/3
987	val3 = 2.+8.0/3
988
989	zl=zr=-3.75 #Assume constant bed (must be less than stage)
990	domain.set_quantity('elevation', zl*num.ones( (4,3) ))
991
992
993	edgeflux = num.zeros(3, num.Float)
994	edgeflux0 = num.zeros(3, num.Float)
995	edgeflux1 = num.zeros(3, num.Float)
996	edgeflux2 = num.zeros(3, num.Float)
997	H0 = 0.0
998
999
1000
1001	domain.set_quantity('stage', [[val0, val0-1, val0-2],
1002	[val1, val1+1, val1],
1003	[val2, val2-2, val2],
1004	[val3-0.5, val3, val3]])
1005
1006	domain.set_quantity('xmomentum',
1007	[[1, 2, 3], [3, 4, 5],
1008	[1, -1, 0], [0, -2, 2]])
1009
1010	domain.set_quantity('ymomentum',
1011	[[1, -1, 0], [0, -3, 2],
1012	[0, 1, 0], [-1, 2, 2]])
1013
1014
1015	domain.check_integrity()
1016
1017
1018
1019	#Flux across right edge of volume 1
1020	normal = domain.get_normal(1,0)
1021	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
1022	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
1023	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)
1024
1025	#Flux across upper edge of volume 1
1026	normal = domain.get_normal(1,1)
1027	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
1028	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
1029	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)
1030
1031	#Flux across lower left hypotenuse of volume 1
1032	normal = domain.get_normal(1,2)
1033	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
1034	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
1035	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)
1036
1037	#Scale, add up and check that compute_fluxes is correct for vol 1
1038	e0 = domain.edgelengths[1, 0]
1039	e1 = domain.edgelengths[1, 1]
1040	e2 = domain.edgelengths[1, 2]
1041
1042	total_flux = -(e0edgeflux0+e1edgeflux1+e2*edgeflux2)/domain.areas[1]
1043
1044	domain.compute_fluxes()
1045	for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
1046	assert num.allclose(total_flux[i],
1047	domain.quantities[name].explicit_update[1])
1048
1049
1050
1051	def xtest_catching_negative_heights(self):
1052
1053	#OBSOLETE
1054
1055	a = [0.0, 0.0]
1056	b = [0.0, 2.0]
1057	c = [2.0,0.0]
1058	d = [0.0, 4.0]
1059	e = [2.0, 2.0]
1060	f = [4.0,0.0]
1061
1062	points = [a, b, c, d, e, f]
1063	#bac, bce, ecf, dbe
1064	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1065
1066	domain = Domain(points, vertices)
1067	val0 = 2.+2.0/3
1068	val1 = 4.+4.0/3
1069	val2 = 8.+2.0/3
1070	val3 = 2.+8.0/3
1071
1072	zl=zr=4 #Too large
1073	domain.set_quantity('elevation', zl*num.ones( (4,3) ))
1074	domain.set_quantity('stage', [[val0, val0-1, val0-2],
1075	[val1, val1+1, val1],
1076	[val2, val2-2, val2],
1077	[val3-0.5, val3, val3]])
1078
1079	#Should fail
1080	try:
1081	domain.check_integrity()
1082	except:
1083	pass
1084
1085
1086
1087	def test_get_wet_elements(self):
1088
1089	a = [0.0, 0.0]
1090	b = [0.0, 2.0]
1091	c = [2.0,0.0]
1092	d = [0.0, 4.0]
1093	e = [2.0, 2.0]
1094	f = [4.0,0.0]
1095
1096	points = [a, b, c, d, e, f]
1097	#bac, bce, ecf, dbe
1098	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1099
1100	domain = Domain(points, vertices)
1101	val0 = 2.+2.0/3
1102	val1 = 4.+4.0/3
1103	val2 = 8.+2.0/3
1104	val3 = 2.+8.0/3
1105
1106	zl=zr=5
1107	domain.set_quantity('elevation', zl*num.ones( (4,3) ))
1108	domain.set_quantity('stage', [[val0, val0-1, val0-2],
1109	[val1, val1+1, val1],
1110	[val2, val2-2, val2],
1111	[val3-0.5, val3, val3]])
1112
1113
1114
1115	#print domain.get_quantity('elevation').get_values(location='centroids')
1116	#print domain.get_quantity('stage').get_values(location='centroids')
1117	domain.check_integrity()
1118
1119	indices = domain.get_wet_elements()
1120	assert num.allclose(indices, [1,2])
1121
1122	indices = domain.get_wet_elements(indices=[0,1,3])
1123	assert num.allclose(indices, [1])
1124
1125
1126
1127	def test_get_maximum_inundation_1(self):
1128
1129	a = [0.0, 0.0]
1130	b = [0.0, 2.0]
1131	c = [2.0,0.0]
1132	d = [0.0, 4.0]
1133	e = [2.0, 2.0]
1134	f = [4.0,0.0]
1135
1136	points = [a, b, c, d, e, f]
1137	#bac, bce, ecf, dbe
1138	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1139
1140	domain = Domain(points, vertices)
1141
1142	domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
1143	domain.set_quantity('stage', 3)
1144
1145	domain.check_integrity()
1146
1147	indices = domain.get_wet_elements()
1148	assert num.allclose(indices, [0])
1149
1150	q = domain.get_maximum_inundation_elevation()
1151	assert num.allclose(q, domain.get_quantity('elevation').get_values(location='centroids')[0])
1152
1153	x, y = domain.get_maximum_inundation_location()
1154	assert num.allclose([x, y], domain.get_centroid_coordinates()[0])
1155
1156
1157	def test_get_maximum_inundation_2(self):
1158	"""test_get_maximum_inundation_2(self)
1159
1160	Test multiple wet cells with same elevation
1161	"""
1162
1163	a = [0.0, 0.0]
1164	b = [0.0, 2.0]
1165	c = [2.0,0.0]
1166	d = [0.0, 4.0]
1167	e = [2.0, 2.0]
1168	f = [4.0,0.0]
1169
1170	points = [a, b, c, d, e, f]
1171	#bac, bce, ecf, dbe
1172	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1173
1174	domain = Domain(points, vertices)
1175
1176	domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
1177	domain.set_quantity('stage', 4.1)
1178
1179	domain.check_integrity()
1180
1181	indices = domain.get_wet_elements()
1182	assert num.allclose(indices, [0,1,2])
1183
1184	q = domain.get_maximum_inundation_elevation()
1185	assert num.allclose(q, 4)
1186
1187	x, y = domain.get_maximum_inundation_location()
1188	assert num.allclose([x, y], domain.get_centroid_coordinates()[1])
1189
1190
1191	def test_get_maximum_inundation_3(self):
1192	"""test_get_maximum_inundation_3(self)
1193
1194	Test of real runup example:
1195	"""
1196
1197	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
1198
1199	initial_runup_height = -0.4
1200	final_runup_height = -0.3
1201
1202
1203	#--------------------------------------------------------------
1204	# Setup computational domain
1205	#--------------------------------------------------------------
1206	N = 5
1207	points, vertices, boundary = rectangular_cross(N, N)
1208	domain = Domain(points, vertices, boundary)
1209	domain.set_maximum_allowed_speed(1.0)
1210
1211	#--------------------------------------------------------------
1212	# Setup initial conditions
1213	#--------------------------------------------------------------
1214	def topography(x,y):
1215	return -x/2 # linear bed slope
1216
1217
1218	domain.set_quantity('elevation', topography) # Use function for elevation
1219	domain.set_quantity('friction', 0.) # Zero friction
1220	domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
1221
1222
1223	#--------------------------------------------------------------
1224	# Setup boundary conditions
1225	#--------------------------------------------------------------
1226	Br = Reflective_boundary(domain) # Reflective wall
1227	Bd = Dirichlet_boundary([final_runup_height, # Constant inflow
1228	0,
1229	0])
1230
1231	# All reflective to begin with (still water)
1232	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1233
1234
1235	#--------------------------------------------------------------
1236	# Test initial inundation height
1237	#--------------------------------------------------------------
1238
1239	indices = domain.get_wet_elements()
1240	z = domain.get_quantity('elevation').\
1241	get_values(location='centroids', indices=indices)
1242	assert num.alltrue(z<initial_runup_height)
1243
1244	q = domain.get_maximum_inundation_elevation()
1245	assert num.allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
1246
1247	x, y = domain.get_maximum_inundation_location()
1248
1249	qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
1250	assert num.allclose(q, qref)
1251
1252
1253	wet_elements = domain.get_wet_elements()
1254	wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
1255	indices=wet_elements)
1256	assert num.alltrue(wet_elevations<initial_runup_height)
1257	assert num.allclose(wet_elevations, z)
1258
1259
1260	#print domain.get_quantity('elevation').get_maximum_value(indices=wet_elements)
1261	#print domain.get_quantity('elevation').get_maximum_location(indices=wet_elements)
1262	#print domain.get_quantity('elevation').get_maximum_index(indices=wet_elements)
1263
1264
1265	#--------------------------------------------------------------
1266	# Let triangles adjust
1267	#--------------------------------------------------------------
1268	for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
1269	pass
1270
1271
1272	#--------------------------------------------------------------
1273	# Test inundation height again
1274	#--------------------------------------------------------------
1275
1276	indices = domain.get_wet_elements()
1277	z = domain.get_quantity('elevation').\
1278	get_values(location='centroids', indices=indices)
1279
1280	assert num.alltrue(z<initial_runup_height)
1281
1282	q = domain.get_maximum_inundation_elevation()
1283	assert num.allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
1284
1285	x, y = domain.get_maximum_inundation_location()
1286	qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
1287	assert num.allclose(q, qref)
1288
1289
1290	#--------------------------------------------------------------
1291	# Update boundary to allow inflow
1292	#--------------------------------------------------------------
1293	domain.set_boundary({'right': Bd})
1294
1295
1296	#--------------------------------------------------------------
1297	# Evolve system through time
1298	#--------------------------------------------------------------
1299	for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
1300	#print domain.timestepping_statistics(track_speeds=True)
1301	#domain.write_time()
1302	pass
1303
1304	#--------------------------------------------------------------
1305	# Test inundation height again
1306	#--------------------------------------------------------------
1307
1308	indices = domain.get_wet_elements()
1309	z = domain.get_quantity('elevation').\
1310	get_values(location='centroids', indices=indices)
1311
1312	assert num.alltrue(z<final_runup_height)
1313
1314	q = domain.get_maximum_inundation_elevation()
1315	assert num.allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
1316
1317	x, y = domain.get_maximum_inundation_location()
1318	qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
1319	assert num.allclose(q, qref)
1320
1321
1322	wet_elements = domain.get_wet_elements()
1323	wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
1324	indices=wet_elements)
1325	assert num.alltrue(wet_elevations<final_runup_height)
1326	assert num.allclose(wet_elevations, z)
1327
1328
1329
1330	def test_get_maximum_inundation_from_sww(self):
1331	"""test_get_maximum_inundation_from_sww(self)
1332
1333	Test of get_maximum_inundation_elevation()
1334	and get_maximum_inundation_location() from data_manager.py
1335
1336	This is based on test_get_maximum_inundation_3(self) but works with the
1337	stored results instead of with the internal data structure.
1338
1339	This test uses the underlying get_maximum_inundation_data for tests
1340	"""
1341
1342	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
1343	from data_manager import get_maximum_inundation_elevation
1344	from data_manager import get_maximum_inundation_location
1345	from data_manager import get_maximum_inundation_data
1346
1347
1348	initial_runup_height = -0.4
1349	final_runup_height = -0.3
1350
1351
1352	#--------------------------------------------------------------
1353	# Setup computational domain
1354	#--------------------------------------------------------------
1355	N = 10
1356	points, vertices, boundary = rectangular_cross(N, N)
1357	domain = Domain(points, vertices, boundary)
1358	domain.set_name('runup_test')
1359	domain.set_maximum_allowed_speed(1.0)
1360
1361	domain.tight_slope_limiters = 0 # FIXME: This works better with old limiters so far
1362
1363	#--------------------------------------------------------------
1364	# Setup initial conditions
1365	#--------------------------------------------------------------
1366	def topography(x,y):
1367	return -x/2 # linear bed slope
1368
1369
1370	domain.set_quantity('elevation', topography) # Use function for elevation
1371	domain.set_quantity('friction', 0.) # Zero friction
1372	domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
1373
1374
1375	#--------------------------------------------------------------
1376	# Setup boundary conditions
1377	#--------------------------------------------------------------
1378	Br = Reflective_boundary(domain) # Reflective wall
1379	Bd = Dirichlet_boundary([final_runup_height, # Constant inflow
1380	0,
1381	0])
1382
1383	# All reflective to begin with (still water)
1384	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1385
1386
1387	#--------------------------------------------------------------
1388	# Test initial inundation height
1389	#--------------------------------------------------------------
1390
1391	indices = domain.get_wet_elements()
1392	z = domain.get_quantity('elevation').\
1393	get_values(location='centroids', indices=indices)
1394	assert num.alltrue(z<initial_runup_height)
1395
1396	q_ref = domain.get_maximum_inundation_elevation()
1397	assert num.allclose(q_ref, initial_runup_height, rtol = 1.0/N) # First order accuracy
1398
1399
1400	#--------------------------------------------------------------
1401	# Let triangles adjust
1402	#--------------------------------------------------------------
1403	for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
1404	pass
1405
1406
1407	#--------------------------------------------------------------
1408	# Test inundation height again
1409	#--------------------------------------------------------------
1410
1411	q_ref = domain.get_maximum_inundation_elevation()
1412	q = get_maximum_inundation_elevation('runup_test.sww')
1413	msg = 'We got %f, should have been %f' %(q, q_ref)
1414	assert num.allclose(q, q_ref, rtol=1.0/N), msg
1415
1416
1417	q = get_maximum_inundation_elevation('runup_test.sww')
1418	msg = 'We got %f, should have been %f' %(q, initial_runup_height)
1419	assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
1420
1421
1422	# Test error condition if time interval is out
1423	try:
1424	q = get_maximum_inundation_elevation('runup_test.sww',
1425	time_interval=[2.0, 3.0])
1426	except ValueError:
1427	pass
1428	else:
1429	msg = 'should have caught wrong time interval'
1430	raise Exception, msg
1431
1432	# Check correct time interval
1433	q, loc = get_maximum_inundation_data('runup_test.sww',
1434	time_interval=[0.0, 3.0])
1435	msg = 'We got %f, should have been %f' %(q, initial_runup_height)
1436	assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
1437	assert num.allclose(-loc[0]/2, q) # From topography formula
1438
1439
1440	#--------------------------------------------------------------
1441	# Update boundary to allow inflow
1442	#--------------------------------------------------------------
1443	domain.set_boundary({'right': Bd})
1444
1445
1446	#--------------------------------------------------------------
1447	# Evolve system through time
1448	#--------------------------------------------------------------
1449	q_max = None
1450	for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
1451	skip_initial_step = True):
1452	q = domain.get_maximum_inundation_elevation()
1453	if q > q_max: q_max = q
1454
1455
1456	#--------------------------------------------------------------
1457	# Test inundation height again
1458	#--------------------------------------------------------------
1459
1460	indices = domain.get_wet_elements()
1461	z = domain.get_quantity('elevation').\
1462	get_values(location='centroids', indices=indices)
1463
1464	assert num.alltrue(z<final_runup_height)
1465
1466	q = domain.get_maximum_inundation_elevation()
1467	assert num.allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
1468
1469	q, loc = get_maximum_inundation_data('runup_test.sww', time_interval=[3.0, 3.0])
1470	msg = 'We got %f, should have been %f' %(q, final_runup_height)
1471	assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
1472	#print 'loc', loc, q
1473	assert num.allclose(-loc[0]/2, q) # From topography formula
1474
1475	q = get_maximum_inundation_elevation('runup_test.sww')
1476	loc = get_maximum_inundation_location('runup_test.sww')
1477	msg = 'We got %f, should have been %f' %(q, q_max)
1478	assert num.allclose(q, q_max, rtol=1.0/N), msg
1479	#print 'loc', loc, q
1480	assert num.allclose(-loc[0]/2, q) # From topography formula
1481
1482
1483
1484	q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[0, 3])
1485	msg = 'We got %f, should have been %f' %(q, q_max)
1486	assert num.allclose(q, q_max, rtol=1.0/N), msg
1487
1488
1489	# Check polygon mode
1490	polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]] # Runup region
1491	q = get_maximum_inundation_elevation('runup_test.sww',
1492	polygon = polygon,
1493	time_interval=[0, 3])
1494	msg = 'We got %f, should have been %f' %(q, q_max)
1495	assert num.allclose(q, q_max, rtol=1.0/N), msg
1496
1497
1498	polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]] # Offshore region
1499	q, loc = get_maximum_inundation_data('runup_test.sww',
1500	polygon = polygon,
1501	time_interval=[0, 3])
1502	msg = 'We got %f, should have been %f' %(q, -0.475)
1503	assert num.allclose(q, -0.475, rtol=1.0/N), msg
1504	assert is_inside_polygon(loc, polygon)
1505	assert num.allclose(-loc[0]/2, q) # From topography formula
1506
1507
1508	polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]] # Dry region
1509	q, loc = get_maximum_inundation_data('runup_test.sww',
1510	polygon = polygon,
1511	time_interval=[0, 3])
1512	msg = 'We got %s, should have been None' %(q)
1513	assert q is None, msg
1514	msg = 'We got %s, should have been None' %(loc)
1515	assert loc is None, msg
1516
1517	# Check what happens if no time point is within interval
1518	try:
1519	q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[2.75, 2.75])
1520	except AssertionError:
1521	pass
1522	else:
1523	msg = 'Time interval should have raised an exception'
1524	raise msg
1525
1526	# Cleanup
1527	try:
1528	os.remove(domain.get_name() + '.' + domain.format)
1529	except:
1530	pass
1531	#FIXME(Ole): Windows won't allow removal of this
1532
1533
1534
1535	def test_get_flow_through_cross_section_with_geo(self):
1536	"""test_get_flow_through_cross_section(self):
1537
1538	Test that the total flow through a cross section can be
1539	correctly obtained at run-time from the ANUGA domain.
1540
1541	This test creates a flat bed with a known flow through it and tests
1542	that the function correctly returns the expected flow.
1543
1544	The specifics are
1545	e = -1 m
1546	u = 2 m/s
1547	h = 2 m
1548	w = 3 m (width of channel)
1549
1550	q = uhw = 12 m^3/s
1551
1552
1553	This run tries it with georeferencing and with elevation = -1
1554
1555	"""
1556
1557	import time, os
1558	from Scientific.IO.NetCDF import NetCDFFile
1559
1560	# Setup
1561	from mesh_factory import rectangular
1562
1563	# Create basic mesh (20m x 3m)
1564	width = 3
1565	length = 20
1566	t_end = 1
1567	points, vertices, boundary = rectangular(length, width,
1568	length, width)
1569
1570	# Create shallow water domain
1571	domain = Domain(points, vertices, boundary,
1572	geo_reference=Geo_reference(56,308500,6189000))
1573
1574	domain.default_order = 2
1575	domain.set_quantities_to_be_stored(None)
1576
1577
1578	e = -1.0
1579	w = 1.0
1580	h = w-e
1581	u = 2.0
1582	uh = u*h
1583
1584	Br = Reflective_boundary(domain) # Side walls
1585	Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet:
1586
1587
1588	# Initial conditions
1589	domain.set_quantity('elevation', e)
1590	domain.set_quantity('stage', w)
1591	domain.set_quantity('xmomentum', uh)
1592	domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
1593
1594
1595	# Interpolation points down the middle
1596	I = [[0, width/2.],
1597	[length/2., width/2.],
1598	[length, width/2.]]
1599	interpolation_points = domain.geo_reference.get_absolute(I)
1600
1601	# Shortcuts to quantites
1602	stage = domain.get_quantity('stage')
1603	xmomentum = domain.get_quantity('xmomentum')
1604	ymomentum = domain.get_quantity('ymomentum')
1605
1606	for t in domain.evolve(yieldstep=0.1, finaltime = t_end):
1607	# Check that quantities are they should be in the interior
1608
1609	w_t = stage.get_values(interpolation_points)
1610	uh_t = xmomentum.get_values(interpolation_points)
1611	vh_t = ymomentum.get_values(interpolation_points)
1612
1613	assert num.allclose(w_t, w)
1614	assert num.allclose(uh_t, uh)
1615	assert num.allclose(vh_t, 0.0)
1616
1617
1618	# Check flows through the middle
1619	for i in range(5):
1620	x = length/2. + i*0.23674563 # Arbitrary
1621	cross_section = [[x, 0], [x, width]]
1622
1623	cross_section = domain.geo_reference.get_absolute(cross_section)
1624	Q = domain.get_flow_through_cross_section(cross_section,
1625	verbose=False)
1626
1627	assert num.allclose(Q, uh*width)
1628
1629
1630
1631	def test_get_energy_through_cross_section_with_geo(self):
1632	"""test_get_energy_through_cross_section(self):
1633
1634	Test that the total and specific energy through a cross section can be
1635	correctly obtained at run-time from the ANUGA domain.
1636
1637	This test creates a flat bed with a known flow through it and tests
1638	that the function correctly returns the expected energies.
1639
1640	The specifics are
1641	e = -1 m
1642	u = 2 m/s
1643	h = 2 m
1644	w = 3 m (width of channel)
1645
1646	q = uhw = 12 m^3/s
1647
1648
1649	This run tries it with georeferencing and with elevation = -1
1650
1651	"""
1652
1653	import time, os
1654	from Scientific.IO.NetCDF import NetCDFFile
1655
1656	# Setup
1657	from mesh_factory import rectangular
1658
1659	# Create basic mesh (20m x 3m)
1660	width = 3
1661	length = 20
1662	t_end = 1
1663	points, vertices, boundary = rectangular(length, width,
1664	length, width)
1665
1666	# Create shallow water domain
1667	domain = Domain(points, vertices, boundary,
1668	geo_reference=Geo_reference(56,308500,6189000))
1669
1670	domain.default_order = 2
1671	domain.set_quantities_to_be_stored(None)
1672
1673
1674	e = -1.0
1675	w = 1.0
1676	h = w-e
1677	u = 2.0
1678	uh = u*h
1679
1680	Br = Reflective_boundary(domain) # Side walls
1681	Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet:
1682
1683
1684	# Initial conditions
1685	domain.set_quantity('elevation', e)
1686	domain.set_quantity('stage', w)
1687	domain.set_quantity('xmomentum', uh)
1688	domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
1689
1690
1691	# Interpolation points down the middle
1692	I = [[0, width/2.],
1693	[length/2., width/2.],
1694	[length, width/2.]]
1695	interpolation_points = domain.geo_reference.get_absolute(I)
1696
1697	# Shortcuts to quantites
1698	stage = domain.get_quantity('stage')
1699	xmomentum = domain.get_quantity('xmomentum')
1700	ymomentum = domain.get_quantity('ymomentum')
1701
1702	for t in domain.evolve(yieldstep=0.1, finaltime = t_end):
1703	# Check that quantities are they should be in the interior
1704
1705	w_t = stage.get_values(interpolation_points)
1706	uh_t = xmomentum.get_values(interpolation_points)
1707	vh_t = ymomentum.get_values(interpolation_points)
1708
1709	assert num.allclose(w_t, w)
1710	assert num.allclose(uh_t, uh)
1711	assert num.allclose(vh_t, 0.0)
1712
1713
1714	# Check energies through the middle
1715	for i in range(5):
1716	x = length/2. + i*0.23674563 # Arbitrary
1717	cross_section = [[x, 0], [x, width]]
1718
1719	cross_section = domain.geo_reference.get_absolute(cross_section)
1720	Es = domain.get_energy_through_cross_section(cross_section,
1721	kind='specific',
1722	verbose=False)
1723
1724	assert num.allclose(Es, h + 0.5uu/g)
1725
1726	Et = domain.get_energy_through_cross_section(cross_section,
1727	kind='total',
1728	verbose=False)
1729	assert num.allclose(Et, w + 0.5uu/g)
1730
1731
1732
1733
1734
1735	def test_another_runup_example(self):
1736	"""test_another_runup_example
1737
1738	Test runup example where actual timeseries at interpolated
1739	points are tested.
1740	"""
1741
1742	#-----------------------------------------------------------------
1743	# Import necessary modules
1744	#-----------------------------------------------------------------
1745
1746	from anuga.pmesh.mesh_interface import create_mesh_from_regions
1747	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
1748	from anuga.shallow_water import Domain
1749	from anuga.shallow_water import Reflective_boundary
1750	from anuga.shallow_water import Dirichlet_boundary
1751
1752
1753	#-----------------------------------------------------------------
1754	# Setup computational domain
1755	#-----------------------------------------------------------------
1756	points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
1757	domain = Domain(points, vertices, boundary) # Create domain
1758	domain.set_default_order(1)
1759	domain.set_quantities_to_be_stored(None)
1760	domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
1761
1762	# FIXME (Ole): Need tests where this is commented out
1763	domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
1764	domain.H0 = 0 # Backwards compatibility (6/2/7)
1765	domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
1766
1767
1768	#-----------------------------------------------------------------
1769	# Setup initial conditions
1770	#-----------------------------------------------------------------
1771
1772	def topography(x,y):
1773	return -x/2 # linear bed slope
1774
1775	domain.set_quantity('elevation', topography)
1776	domain.set_quantity('friction', 0.0)
1777	domain.set_quantity('stage', expression='elevation')
1778
1779
1780	#----------------------------------------------------------------
1781	# Setup boundary conditions
1782	#----------------------------------------------------------------
1783
1784	Br = Reflective_boundary(domain) # Solid reflective wall
1785	Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
1786	domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
1787
1788
1789	#----------------------------------------------------------------
1790	# Evolve system through time
1791	#----------------------------------------------------------------
1792
1793	interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
1794	gauge_values = []
1795	for _ in interpolation_points:
1796	gauge_values.append([])
1797
1798	time = []
1799	for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
1800	# Record time series at known points
1801	time.append(domain.get_time())
1802
1803	stage = domain.get_quantity('stage')
1804	w = stage.get_values(interpolation_points=interpolation_points)
1805
1806	for i, _ in enumerate(interpolation_points):
1807	gauge_values[i].append(w[i])
1808
1809
1810	#print
1811	#print time
1812	#print
1813	#for i, (x,y) in enumerate(interpolation_points):
1814	# print i, gauge_values[i]
1815	# print
1816
1817	#Reference (nautilus 26/6/2008)
1818
1819	G0 = [-0.20000000000000001, -0.20000000000000001, -0.19920600846161715, -0.19153647344085376, -0.19127622768281194, -0.1770671909675095, -0.16739412133181927, -0.16196038919122191, -0.15621633053131384, -0.15130021599977705, -0.13930978857215484, -0.19349274358263582, -0.19975307598803765, -0.19999897143103357, -0.1999999995532111, -0.19999999999949952, -0.19999999999949952, -0.19999999999949952, -0.19997270012494556, -0.19925805948554556, -0.19934513778450533, -0.19966484196394893, -0.1997352860102084, -0.19968260481750394, -0.19980280797303882, -0.19998804881822749, -0.19999999778075916, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167]
1820
1821	G1 = [-0.29999999999999993, -0.29999588068034899, -0.29250047332330331, -0.28335081844518584, -0.26142206997410805, -0.22656028856329835, -0.21224087216745585, -0.19934324109114465, -0.1889857939783175, -0.18146311603911383, -0.17401078727434263, -0.15419361061257214, -0.16225060576782063, -0.19010941396999181, -0.20901161407004412, -0.21670683975774699, -0.21771386270738891, -0.21481284465869752, -0.21063120869004387, -0.20669243364582401, -0.20320707386714859, -0.19984087691926442, -0.19725417448019505, -0.19633783049069981, -0.19650494599999785, -0.19708316838336942, -0.19779309449413818, -0.19853070294429562, -0.19917342167307153, -0.19964814677795845, -0.19991627610824922, -0.20013162970144974, -0.20029864969405509, -0.20036259676501131, -0.20030682824965193, -0.20016105135750167, -0.19997664501985918, -0.19980185871568762, -0.19966836175417696, -0.19958856744312226, -0.19955954696194517, -0.19956950051110917, -0.19960377086336181, -0.19964885299433241, -0.19969427478531132, -0.19973301547655564, -0.19976121574277764, -0.19977765285688653, -0.19978315117522441, -0.19977994634841739, -0.19977101394878494]
1822
1823	G2 = [-0.40000000000000002, -0.39077401254732241, -0.33350466136630474, -0.29771023004255281, -0.27605439066140897, -0.25986156218997497, -0.24502185018573647, -0.231792624329521, -0.21981564668803993, -0.20870707082936543, -0.19877739883776599, -0.18980922837977957, -0.17308011674005838, -0.16306400164013773, -0.17798470933304333, -0.1929554075869116, -0.20236705191987037, -0.20695767560655007, -0.20841025876092567, -0.20792102174869989, -0.20655350005579293, -0.20492002526259828, -0.20310627026780645, -0.20105983335287836, -0.19937394565794653, -0.19853917506699659, -0.19836389977624452, -0.19850305023602796, -0.19877764028836831, -0.19910928131034669, -0.19943705712418805, -0.19970344172958865, -0.19991076989870474, -0.20010020127747646, -0.20025937787100062, -0.20035087292905965, -0.20035829921463297, -0.20029606557316171, -0.20019606915365515, -0.20009096093399206, -0.20000371608204368, -0.19994495432920584, -0.19991535665176338, -0.19990981826533513, -0.19992106419898723, -0.19994189853516578, -0.19996624091229293, -0.19998946016985167, -0.20000842303470234, -0.20002144460718174, -0.20002815561337187]
1824
1825	G3 = [-0.45000000000000001, -0.37631169657400332, -0.33000044342859486, -0.30586045469008522, -0.28843572253009941, -0.27215308978603808, -0.25712951540331219, -0.2431608296216613, -0.23032023651386374, -0.2184546873456619, -0.20735123704254332, -0.19740397194806389, -0.1859829564064375, -0.16675980728362105, -0.16951575032846536, -0.1832860872609344, -0.19485758939241243, -0.20231368291811427, -0.20625610376074754, -0.20758116241495619, -0.20721445402086161, -0.20603406830353785, -0.20450262808396991, -0.2026769581185151, -0.2007401212066364, -0.19931160535777592, -0.19863606301128725, -0.19848511940572691, -0.19860091042948352, -0.19885490669377764, -0.19916542732701112, -0.19946678238611959, -0.19971209594104697, -0.19991912886512292, -0.2001058430788881, -0.20024959409472989, -0.20032160254609382, -0.20031583165752354, -0.20025051539293123, -0.2001556115816068, -0.20005952955420872, -0.1999814429561611, -0.19992977821558131, -0.19990457708664208, -0.19990104785490476, -0.19991257153954825, -0.19993258231880562, -0.19995548502882532, -0.19997700760919687, -0.19999429663503748, -0.20000588800248761]
1826
1827	#FIXME (DSG):This is a hack so the anuga install, not precompiled
1828	# works on DSG's win2000, python 2.3
1829	#The problem is the gauge_values[X] are 52 long, not 51.
1830	#
1831	# This was probably fixed by Stephen in changeset:3804
1832	#if len(gauge_values[0]) == 52: gauge_values[0].pop()
1833	#if len(gauge_values[1]) == 52: gauge_values[1].pop()
1834	#if len(gauge_values[2]) == 52: gauge_values[2].pop()
1835	#if len(gauge_values[3]) == 52: gauge_values[3].pop()
1836
1837	## print len(G0), len(gauge_values[0])
1838	## print len(G1), len(gauge_values[1])
1839
1840	#print gauge_values[3]
1841	#print G0[:4]
1842	#print array(gauge_values[0])-array(G0)
1843
1844
1845	assert num.allclose(gauge_values[0], G0)
1846	assert num.allclose(gauge_values[1], G1)
1847	assert num.allclose(gauge_values[2], G2)
1848	assert num.allclose(gauge_values[3], G3)
1849
1850
1851
1852
1853
1854
1855
1856	#####################################################
1857
1858	def test_flux_optimisation(self):
1859	"""test_flux_optimisation
1860	Test that fluxes are correctly computed using
1861	dry and still cell exclusions
1862	"""
1863
1864	from anuga.config import g
1865	import copy
1866
1867	a = [0.0, 0.0]
1868	b = [0.0, 2.0]
1869	c = [2.0, 0.0]
1870	d = [0.0, 4.0]
1871	e = [2.0, 2.0]
1872	f = [4.0, 0.0]
1873
1874	points = [a, b, c, d, e, f]
1875	#bac, bce, ecf, dbe
1876	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1877
1878	domain = Domain(points, vertices)
1879
1880	#Set up for a gradient of (3,0) at mid triangle (bce)
1881	def slope(x, y):
1882	return 3*x
1883
1884	h = 0.1
1885	def stage(x,y):
1886	return slope(x,y)+h
1887
1888	domain.set_quantity('elevation', slope)
1889	domain.set_quantity('stage', stage)
1890
1891	# Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
1892	domain.distribute_to_vertices_and_edges()
1893
1894	initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
1895
1896	domain.set_boundary({'exterior': Reflective_boundary(domain)})
1897
1898
1899	# Check that update arrays are initialised to zero
1900	assert num.allclose(domain.get_quantity('stage').explicit_update, 0)
1901	assert num.allclose(domain.get_quantity('xmomentum').explicit_update, 0)
1902	assert num.allclose(domain.get_quantity('ymomentum').explicit_update, 0)
1903
1904
1905	# Get true values
1906	domain.optimise_dry_cells = False
1907	domain.compute_fluxes()
1908	stage_ref = copy.copy(domain.get_quantity('stage').explicit_update)
1909	xmom_ref = copy.copy(domain.get_quantity('xmomentum').explicit_update)
1910	ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update)
1911
1912	# Try with flux optimisation
1913	domain.optimise_dry_cells = True
1914	domain.compute_fluxes()
1915
1916	assert num.allclose(stage_ref, domain.get_quantity('stage').explicit_update)
1917	assert num.allclose(xmom_ref, domain.get_quantity('xmomentum').explicit_update)
1918	assert num.allclose(ymom_ref, domain.get_quantity('ymomentum').explicit_update)
1919
1920
1921
1922	def test_initial_condition(self):
1923	"""test_initial_condition
1924	Test that initial condition is output at time == 0 and that
1925	computed values change as system evolves
1926	"""
1927
1928	from anuga.config import g
1929	import copy
1930
1931	a = [0.0, 0.0]
1932	b = [0.0, 2.0]
1933	c = [2.0, 0.0]
1934	d = [0.0, 4.0]
1935	e = [2.0, 2.0]
1936	f = [4.0, 0.0]
1937
1938	points = [a, b, c, d, e, f]
1939	#bac, bce, ecf, dbe
1940	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1941
1942	domain = Domain(points, vertices)
1943
1944	#Set up for a gradient of (3,0) at mid triangle (bce)
1945	def slope(x, y):
1946	return 3*x
1947
1948	h = 0.1
1949	def stage(x,y):
1950	return slope(x,y)+h
1951
1952	domain.set_quantity('elevation', slope)
1953	domain.set_quantity('stage', stage)
1954
1955	# Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
1956	domain.distribute_to_vertices_and_edges()
1957
1958	initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
1959
1960	domain.set_boundary({'exterior': Reflective_boundary(domain)})
1961
1962	domain.optimise_dry_cells = True
1963	#Evolution
1964	for t in domain.evolve(yieldstep = 0.5, finaltime = 2.0):
1965	stage = domain.quantities['stage'].vertex_values
1966
1967	if t == 0.0:
1968	assert num.allclose(stage, initial_stage)
1969	else:
1970	assert not num.allclose(stage, initial_stage)
1971
1972
1973	os.remove(domain.get_name() + '.sww')
1974
1975
1976
1977	#####################################################
1978	def test_gravity(self):
1979	#Assuming no friction
1980
1981	from anuga.config import g
1982
1983	a = [0.0, 0.0]
1984	b = [0.0, 2.0]
1985	c = [2.0, 0.0]
1986	d = [0.0, 4.0]
1987	e = [2.0, 2.0]
1988	f = [4.0, 0.0]
1989
1990	points = [a, b, c, d, e, f]
1991	#bac, bce, ecf, dbe
1992	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1993
1994	domain = Domain(points, vertices)
1995
1996	#Set up for a gradient of (3,0) at mid triangle (bce)
1997	def slope(x, y):
1998	return 3*x
1999
2000	h = 0.1
2001	def stage(x,y):
2002	return slope(x,y)+h
2003
2004	domain.set_quantity('elevation', slope)
2005	domain.set_quantity('stage', stage)
2006
2007	for name in domain.conserved_quantities:
2008	assert num.allclose(domain.quantities[name].explicit_update, 0)
2009	assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
2010
2011	domain.compute_forcing_terms()
2012
2013	assert num.allclose(domain.quantities['stage'].explicit_update, 0)
2014	assert num.allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
2015	assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
2016
2017
2018	def test_manning_friction(self):
2019	from anuga.config import g
2020
2021	a = [0.0, 0.0]
2022	b = [0.0, 2.0]
2023	c = [2.0, 0.0]
2024	d = [0.0, 4.0]
2025	e = [2.0, 2.0]
2026	f = [4.0, 0.0]
2027
2028	points = [a, b, c, d, e, f]
2029	#bac, bce, ecf, dbe
2030	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2031
2032	domain = Domain(points, vertices)
2033
2034	#Set up for a gradient of (3,0) at mid triangle (bce)
2035	def slope(x, y):
2036	return 3*x
2037
2038	h = 0.1
2039	def stage(x,y):
2040	return slope(x,y)+h
2041
2042	eta = 0.07
2043	domain.set_quantity('elevation', slope)
2044	domain.set_quantity('stage', stage)
2045	domain.set_quantity('friction', eta)
2046
2047	for name in domain.conserved_quantities:
2048	assert num.allclose(domain.quantities[name].explicit_update, 0)
2049	assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
2050
2051	domain.compute_forcing_terms()
2052
2053	assert num.allclose(domain.quantities['stage'].explicit_update, 0)
2054	assert num.allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
2055	assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
2056
2057	assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
2058	assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
2059	assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
2060
2061	#Create some momentum for friction to work with
2062	domain.set_quantity('xmomentum', 1)
2063	S = -g * eta2 / h(7.0/3)
2064
2065	domain.compute_forcing_terms()
2066	assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
2067	assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
2068	assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
2069
2070	#A more complex example
2071	domain.quantities['stage'].semi_implicit_update[:] = 0.0
2072	domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
2073	domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
2074
2075	domain.set_quantity('xmomentum', 3)
2076	domain.set_quantity('ymomentum', 4)
2077
2078	S = -g * eta*2 5 / h**(7.0/3)
2079
2080
2081	domain.compute_forcing_terms()
2082
2083	assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
2084	assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
2085	assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
2086
2087	def test_constant_wind_stress(self):
2088	from anuga.config import rho_a, rho_w, eta_w
2089	from math import pi, cos, sin
2090
2091	a = [0.0, 0.0]
2092	b = [0.0, 2.0]
2093	c = [2.0, 0.0]
2094	d = [0.0, 4.0]
2095	e = [2.0, 2.0]
2096	f = [4.0, 0.0]
2097
2098	points = [a, b, c, d, e, f]
2099	#bac, bce, ecf, dbe
2100	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2101
2102
2103	domain = Domain(points, vertices)
2104
2105	#Flat surface with 1m of water
2106	domain.set_quantity('elevation', 0)
2107	domain.set_quantity('stage', 1.0)
2108	domain.set_quantity('friction', 0)
2109
2110	Br = Reflective_boundary(domain)
2111	domain.set_boundary({'exterior': Br})
2112
2113	#Setup only one forcing term, constant wind stress
2114	s = 100
2115	phi = 135
2116	domain.forcing_terms = []
2117	domain.forcing_terms.append( Wind_stress(s, phi) )
2118
2119	domain.compute_forcing_terms()
2120
2121
2122	const = eta_w*rho_a/rho_w
2123
2124	#Convert to radians
2125	phi = phi*pi/180
2126
2127	#Compute velocity vector (u, v)
2128	u = s*cos(phi)
2129	v = s*sin(phi)
2130
2131	#Compute wind stress
2132	S = const * num.sqrt(u2 + v2)
2133
2134	assert num.allclose(domain.quantities['stage'].explicit_update, 0)
2135	assert num.allclose(domain.quantities['xmomentum'].explicit_update, S*u)
2136	assert num.allclose(domain.quantities['ymomentum'].explicit_update, S*v)
2137
2138
2139	def test_variable_wind_stress(self):
2140	from anuga.config import rho_a, rho_w, eta_w
2141	from math import pi, cos, sin
2142
2143	a = [0.0, 0.0]
2144	b = [0.0, 2.0]
2145	c = [2.0, 0.0]
2146	d = [0.0, 4.0]
2147	e = [2.0, 2.0]
2148	f = [4.0, 0.0]
2149
2150	points = [a, b, c, d, e, f]
2151	#bac, bce, ecf, dbe
2152	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2153
2154	domain = Domain(points, vertices)
2155
2156	#Flat surface with 1m of water
2157	domain.set_quantity('elevation', 0)
2158	domain.set_quantity('stage', 1.0)
2159	domain.set_quantity('friction', 0)
2160
2161	Br = Reflective_boundary(domain)
2162	domain.set_boundary({'exterior': Br})
2163
2164
2165	domain.time = 5.54 #Take a random time (not zero)
2166
2167	#Setup only one forcing term, constant wind stress
2168	s = 100
2169	phi = 135
2170	domain.forcing_terms = []
2171	domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
2172
2173	domain.compute_forcing_terms()
2174
2175	#Compute reference solution
2176	const = eta_w*rho_a/rho_w
2177
2178	N = len(domain) # number_of_triangles
2179
2180	xc = domain.get_centroid_coordinates()
2181	t = domain.time
2182
2183	x = xc[:,0]
2184	y = xc[:,1]
2185	s_vec = speed(t,x,y)
2186	phi_vec = angle(t,x,y)
2187
2188
2189	for k in range(N):
2190	#Convert to radians
2191	phi = phi_vec[k]*pi/180
2192	s = s_vec[k]
2193
2194	#Compute velocity vector (u, v)
2195	u = s*cos(phi)
2196	v = s*sin(phi)
2197
2198	#Compute wind stress
2199	S = const * num.sqrt(u2 + v2)
2200
2201	assert num.allclose(domain.quantities['stage'].explicit_update[k], 0)
2202	assert num.allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
2203	assert num.allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
2204
2205
2206
2207
2208
2209
2210	def test_windfield_from_file(self):
2211	from anuga.config import rho_a, rho_w, eta_w
2212	from math import pi, cos, sin
2213	from anuga.config import time_format
2214	from anuga.abstract_2d_finite_volumes.util import file_function
2215	import time
2216
2217
2218	a = [0.0, 0.0]
2219	b = [0.0, 2.0]
2220	c = [2.0, 0.0]
2221	d = [0.0, 4.0]
2222	e = [2.0, 2.0]
2223	f = [4.0, 0.0]
2224
2225	points = [a, b, c, d, e, f]
2226	#bac, bce, ecf, dbe
2227	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2228
2229	domain = Domain(points, vertices)
2230
2231	#Flat surface with 1m of water
2232	domain.set_quantity('elevation', 0)
2233	domain.set_quantity('stage', 1.0)
2234	domain.set_quantity('friction', 0)
2235
2236	Br = Reflective_boundary(domain)
2237	domain.set_boundary({'exterior': Br})
2238
2239
2240	domain.time = 7 #Take a time that is represented in file (not zero)
2241
2242	#Write wind stress file (ensure that domain.time is covered)
2243	#Take x=1 and y=0
2244	filename = 'test_windstress_from_file'
2245	start = time.mktime(time.strptime('2000', '%Y'))
2246	fid = open(filename + '.txt', 'w')
2247	dt = 1 #One second interval
2248	t = 0.0
2249	while t <= 10.0:
2250	t_string = time.strftime(time_format, time.gmtime(t+start))
2251
2252	fid.write('%s, %f %f\n' %(t_string,
2253	speed(t,[1],[0])[0],
2254	angle(t,[1],[0])[0]))
2255	t += dt
2256
2257	fid.close()
2258
2259
2260	#Convert ASCII file to NetCDF (Which is what we really like!)
2261	from data_manager import timefile2netcdf
2262	timefile2netcdf(filename)
2263	os.remove(filename + '.txt')
2264
2265
2266	#Setup wind stress
2267	F = file_function(filename + '.tms', quantities = ['Attribute0',
2268	'Attribute1'])
2269	os.remove(filename + '.tms')
2270
2271
2272	#print 'F(5)', F(5)
2273
2274	#print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
2275
2276	#print dir(F)
2277	#print F.T
2278	#print F.precomputed_values
2279	#
2280	#F = file_function(filename + '.txt')
2281	#
2282	#print dir(F)
2283	#print F.T
2284	#print F.Q
2285
2286	W = Wind_stress(F)
2287
2288	domain.forcing_terms = []
2289	domain.forcing_terms.append(W)
2290
2291	domain.compute_forcing_terms()
2292
2293	#Compute reference solution
2294	const = eta_w*rho_a/rho_w
2295
2296	N = len(domain) # number_of_triangles
2297
2298	t = domain.time
2299
2300	s = speed(t,[1],[0])[0]
2301	phi = angle(t,[1],[0])[0]
2302
2303	#Convert to radians
2304	phi = phi*pi/180
2305
2306
2307	#Compute velocity vector (u, v)
2308	u = s*cos(phi)
2309	v = s*sin(phi)
2310
2311	#Compute wind stress
2312	S = const * num.sqrt(u2 + v2)
2313
2314	for k in range(N):
2315	assert num.allclose(domain.quantities['stage'].explicit_update[k], 0)
2316	assert num.allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
2317	assert num.allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
2318
2319
2320	def test_windfield_from_file_seconds(self):
2321	from anuga.config import rho_a, rho_w, eta_w
2322	from math import pi, cos, sin
2323	from anuga.config import time_format
2324	from anuga.abstract_2d_finite_volumes.util import file_function
2325	import time
2326
2327
2328	a = [0.0, 0.0]
2329	b = [0.0, 2.0]
2330	c = [2.0, 0.0]
2331	d = [0.0, 4.0]
2332	e = [2.0, 2.0]
2333	f = [4.0, 0.0]
2334
2335	points = [a, b, c, d, e, f]
2336	#bac, bce, ecf, dbe
2337	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2338
2339	domain = Domain(points, vertices)
2340
2341	#Flat surface with 1m of water
2342	domain.set_quantity('elevation', 0)
2343	domain.set_quantity('stage', 1.0)
2344	domain.set_quantity('friction', 0)
2345
2346	Br = Reflective_boundary(domain)
2347	domain.set_boundary({'exterior': Br})
2348
2349
2350	domain.time = 7 #Take a time that is represented in file (not zero)
2351
2352	#Write wind stress file (ensure that domain.time is covered)
2353	#Take x=1 and y=0
2354	filename = 'test_windstress_from_file'
2355	start = time.mktime(time.strptime('2000', '%Y'))
2356	fid = open(filename + '.txt', 'w')
2357	dt = 0.5 #1 #One second interval
2358	t = 0.0
2359	while t <= 10.0:
2360	fid.write('%s, %f %f\n' %(str(t),
2361	speed(t,[1],[0])[0],
2362	angle(t,[1],[0])[0]))
2363	t += dt
2364
2365	fid.close()
2366
2367
2368	#Convert ASCII file to NetCDF (Which is what we really like!)
2369	from data_manager import timefile2netcdf
2370	timefile2netcdf(filename, time_as_seconds=True)
2371	os.remove(filename + '.txt')
2372
2373
2374	#Setup wind stress
2375	F = file_function(filename + '.tms', quantities = ['Attribute0',
2376	'Attribute1'])
2377	os.remove(filename + '.tms')
2378
2379
2380	#print 'F(5)', F(5)
2381
2382	#print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
2383
2384	#print dir(F)
2385	#print F.T
2386	#print F.precomputed_values
2387	#
2388	#F = file_function(filename + '.txt')
2389	#
2390	#print dir(F)
2391	#print F.T
2392	#print F.Q
2393
2394	W = Wind_stress(F)
2395
2396	domain.forcing_terms = []
2397	domain.forcing_terms.append(W)
2398
2399	domain.compute_forcing_terms()
2400
2401	#Compute reference solution
2402	const = eta_w*rho_a/rho_w
2403
2404	N = len(domain) # number_of_triangles
2405
2406	t = domain.time
2407
2408	s = speed(t,[1],[0])[0]
2409	phi = angle(t,[1],[0])[0]
2410
2411	#Convert to radians
2412	phi = phi*pi/180
2413
2414
2415	#Compute velocity vector (u, v)
2416	u = s*cos(phi)
2417	v = s*sin(phi)
2418
2419	#Compute wind stress
2420	S = const * num.sqrt(u2 + v2)
2421
2422	for k in range(N):
2423	assert num.allclose(domain.quantities['stage'].explicit_update[k], 0)
2424	assert num.allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
2425	assert num.allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
2426
2427
2428
2429
2430	def test_wind_stress_error_condition(self):
2431	"""Test that windstress reacts properly when forcing functions
2432	are wrong - e.g. returns a scalar
2433	"""
2434
2435	from anuga.config import rho_a, rho_w, eta_w
2436	from math import pi, cos, sin
2437
2438	a = [0.0, 0.0]
2439	b = [0.0, 2.0]
2440	c = [2.0, 0.0]
2441	d = [0.0, 4.0]
2442	e = [2.0, 2.0]
2443	f = [4.0, 0.0]
2444
2445	points = [a, b, c, d, e, f]
2446	#bac, bce, ecf, dbe
2447	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2448
2449	domain = Domain(points, vertices)
2450
2451	#Flat surface with 1m of water
2452	domain.set_quantity('elevation', 0)
2453	domain.set_quantity('stage', 1.0)
2454	domain.set_quantity('friction', 0)
2455
2456	Br = Reflective_boundary(domain)
2457	domain.set_boundary({'exterior': Br})
2458
2459
2460	domain.time = 5.54 #Take a random time (not zero)
2461
2462	#Setup only one forcing term, bad func
2463	domain.forcing_terms = []
2464
2465	try:
2466	domain.forcing_terms.append(Wind_stress(s = scalar_func,
2467	phi = angle))
2468	except AssertionError:
2469	pass
2470	else:
2471	msg = 'Should have raised exception'
2472	raise msg
2473
2474
2475	try:
2476	domain.forcing_terms.append(Wind_stress(s = speed,
2477	phi = scalar_func))
2478	except AssertionError:
2479	pass
2480	else:
2481	msg = 'Should have raised exception'
2482	raise msg
2483
2484	try:
2485	domain.forcing_terms.append(Wind_stress(s = speed,
2486	phi = 'xx'))
2487	except:
2488	pass
2489	else:
2490	msg = 'Should have raised exception'
2491	raise msg
2492
2493
2494
2495	def test_rainfall(self):
2496	from math import pi, cos, sin
2497
2498	a = [0.0, 0.0]
2499	b = [0.0, 2.0]
2500	c = [2.0, 0.0]
2501	d = [0.0, 4.0]
2502	e = [2.0, 2.0]
2503	f = [4.0, 0.0]
2504
2505	points = [a, b, c, d, e, f]
2506	#bac, bce, ecf, dbe
2507	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2508
2509
2510	domain = Domain(points, vertices)
2511
2512	#Flat surface with 1m of water
2513	domain.set_quantity('elevation', 0)
2514	domain.set_quantity('stage', 1.0)
2515	domain.set_quantity('friction', 0)
2516
2517	Br = Reflective_boundary(domain)
2518	domain.set_boundary({'exterior': Br})
2519
2520	# Setup only one forcing term, constant rainfall
2521	domain.forcing_terms = []
2522	domain.forcing_terms.append( Rainfall(domain, rate=2.0) )
2523
2524	domain.compute_forcing_terms()
2525	assert num.allclose(domain.quantities['stage'].explicit_update, 2.0/1000)
2526
2527
2528
2529	def test_rainfall_restricted_by_polygon(self):
2530	from math import pi, cos, sin
2531
2532	a = [0.0, 0.0]
2533	b = [0.0, 2.0]
2534	c = [2.0, 0.0]
2535	d = [0.0, 4.0]
2536	e = [2.0, 2.0]
2537	f = [4.0, 0.0]
2538
2539	points = [a, b, c, d, e, f]
2540	#bac, bce, ecf, dbe
2541	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2542
2543
2544	domain = Domain(points, vertices)
2545
2546	#Flat surface with 1m of water
2547	domain.set_quantity('elevation', 0)
2548	domain.set_quantity('stage', 1.0)
2549	domain.set_quantity('friction', 0)
2550
2551	Br = Reflective_boundary(domain)
2552	domain.set_boundary({'exterior': Br})
2553
2554	# Setup only one forcing term, constant rainfall restricted to a polygon enclosing triangle #1 (bce)
2555	domain.forcing_terms = []
2556	R = Rainfall(domain,
2557	rate=2.0,
2558	polygon = [[1,1], [2,1], [2,2], [1,2]])
2559
2560	assert num.allclose(R.exchange_area, 1)
2561
2562	domain.forcing_terms.append(R)
2563
2564	domain.compute_forcing_terms()
2565	#print domain.quantities['stage'].explicit_update
2566
2567	assert num.allclose(domain.quantities['stage'].explicit_update[1],
2568	2.0/1000)
2569	assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2570	assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2571
2572
2573
2574	def test_time_dependent_rainfall_restricted_by_polygon(self):
2575
2576	a = [0.0, 0.0]
2577	b = [0.0, 2.0]
2578	c = [2.0, 0.0]
2579	d = [0.0, 4.0]
2580	e = [2.0, 2.0]
2581	f = [4.0, 0.0]
2582
2583	points = [a, b, c, d, e, f]
2584	#bac, bce, ecf, dbe
2585	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2586
2587
2588	domain = Domain(points, vertices)
2589
2590	#Flat surface with 1m of water
2591	domain.set_quantity('elevation', 0)
2592	domain.set_quantity('stage', 1.0)
2593	domain.set_quantity('friction', 0)
2594
2595	Br = Reflective_boundary(domain)
2596	domain.set_boundary({'exterior': Br})
2597
2598	# Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce)
2599	domain.forcing_terms = []
2600	R = Rainfall(domain,
2601	rate=lambda t: 3*t + 7,
2602	polygon = [[1,1], [2,1], [2,2], [1,2]])
2603
2604	assert num.allclose(R.exchange_area, 1)
2605
2606	domain.forcing_terms.append(R)
2607
2608
2609	domain.time = 10.
2610
2611	domain.compute_forcing_terms()
2612	#print domain.quantities['stage'].explicit_update
2613
2614	assert num.allclose(domain.quantities['stage'].explicit_update[1],
2615	(3*domain.time+7)/1000)
2616	assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2617	assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2618
2619
2620
2621
2622	def test_time_dependent_rainfall_using_starttime(self):
2623
2624	rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.Float)
2625
2626	a = [0.0, 0.0]
2627	b = [0.0, 2.0]
2628	c = [2.0, 0.0]
2629	d = [0.0, 4.0]
2630	e = [2.0, 2.0]
2631	f = [4.0, 0.0]
2632
2633	points = [a, b, c, d, e, f]
2634	#bac, bce, ecf, dbe
2635	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2636
2637
2638	domain = Domain(points, vertices)
2639
2640	#Flat surface with 1m of water
2641	domain.set_quantity('elevation', 0)
2642	domain.set_quantity('stage', 1.0)
2643	domain.set_quantity('friction', 0)
2644
2645	Br = Reflective_boundary(domain)
2646	domain.set_boundary({'exterior': Br})
2647
2648	# Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce)
2649	domain.forcing_terms = []
2650	R = Rainfall(domain,
2651	rate=lambda t: 3*t + 7,
2652	polygon=rainfall_poly)
2653
2654	assert num.allclose(R.exchange_area, 1)
2655
2656	domain.forcing_terms.append(R)
2657
2658	# This will test that time used in the forcing function takes
2659	# startime into account.
2660	domain.starttime = 5.0
2661
2662	domain.time = 7.
2663
2664	domain.compute_forcing_terms()
2665	#print domain.quantities['stage'].explicit_update
2666
2667	#print domain.get_time()
2668	assert num.allclose(domain.quantities['stage'].explicit_update[1],
2669	(3*domain.get_time()+7)/1000)
2670	assert num.allclose(domain.quantities['stage'].explicit_update[1],
2671	(3*(domain.time + domain.starttime)+7)/1000)
2672
2673	# Using internal time her should fail
2674	assert not num.allclose(domain.quantities['stage'].explicit_update[1],
2675	(3*domain.time+7)/1000)
2676
2677	assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2678	assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2679
2680
2681
2682
2683	def test_time_dependent_rainfall_using_georef(self):
2684	"""test_time_dependent_rainfall_using_georef
2685
2686	This will also test the General forcing term using georef
2687	"""
2688
2689	#Mesh in zone 56 (absolute coords)
2690
2691	x0 = 314036.58727982
2692	y0 = 6224951.2960092
2693
2694
2695	rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.Float)
2696	rainfall_poly += [x0, y0]
2697
2698	a = [0.0, 0.0]
2699	b = [0.0, 2.0]
2700	c = [2.0, 0.0]
2701	d = [0.0, 4.0]
2702	e = [2.0, 2.0]
2703	f = [4.0, 0.0]
2704
2705	points = [a, b, c, d, e, f]
2706	#bac, bce, ecf, dbe
2707	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2708
2709
2710	domain = Domain(points, vertices,
2711	geo_reference = Geo_reference(56, x0, y0))
2712
2713	#Flat surface with 1m of water
2714	domain.set_quantity('elevation', 0)
2715	domain.set_quantity('stage', 1.0)
2716	domain.set_quantity('friction', 0)
2717
2718	Br = Reflective_boundary(domain)
2719	domain.set_boundary({'exterior': Br})
2720
2721	# Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce)
2722	domain.forcing_terms = []
2723	R = Rainfall(domain,
2724	rate=lambda t: 3*t + 7,
2725	polygon=rainfall_poly)
2726
2727	assert num.allclose(R.exchange_area, 1)
2728
2729	domain.forcing_terms.append(R)
2730
2731	# This will test that time used in the forcing function takes
2732	# startime into account.
2733	domain.starttime = 5.0
2734
2735	domain.time = 7.
2736
2737	domain.compute_forcing_terms()
2738	#print domain.quantities['stage'].explicit_update
2739
2740	#print domain.get_time()
2741	assert num.allclose(domain.quantities['stage'].explicit_update[1],
2742	(3*domain.get_time()+7)/1000)
2743	assert num.allclose(domain.quantities['stage'].explicit_update[1],
2744	(3*(domain.time + domain.starttime)+7)/1000)
2745
2746	# Using internal time her should fail
2747	assert not num.allclose(domain.quantities['stage'].explicit_update[1],
2748	(3*domain.time+7)/1000)
2749
2750	assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2751	assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2752
2753
2754
2755
2756
2757
2758	def test_time_dependent_rainfall_restricted_by_polygon_with_default(self):
2759	"""test_time_dependent_rainfall_restricted_by_polygon_with_default
2760
2761	Test that default rainfall can be used when given rate runs out of data.
2762	"""
2763	a = [0.0, 0.0]
2764	b = [0.0, 2.0]
2765	c = [2.0, 0.0]
2766	d = [0.0, 4.0]
2767	e = [2.0, 2.0]
2768	f = [4.0, 0.0]
2769
2770	points = [a, b, c, d, e, f]
2771	#bac, bce, ecf, dbe
2772	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2773
2774
2775	domain = Domain(points, vertices)
2776
2777	#Flat surface with 1m of water
2778	domain.set_quantity('elevation', 0)
2779	domain.set_quantity('stage', 1.0)
2780	domain.set_quantity('friction', 0)
2781
2782	Br = Reflective_boundary(domain)
2783	domain.set_boundary({'exterior': Br})
2784
2785	# Setup only one forcing term, time dependent rainfall that expires at t==20
2786	from anuga.fit_interpolate.interpolate import Modeltime_too_late
2787	def main_rate(t):
2788	if t > 20:
2789	msg = 'Model time exceeded.'
2790	raise Modeltime_too_late, msg
2791	else:
2792	return 3*t + 7
2793
2794	domain.forcing_terms = []
2795	R = Rainfall(domain,
2796	rate=main_rate,
2797	polygon = [[1,1], [2,1], [2,2], [1,2]],
2798	default_rate=5.0)
2799
2800	assert num.allclose(R.exchange_area, 1)
2801
2802	domain.forcing_terms.append(R)
2803
2804
2805	domain.time = 10.
2806
2807	domain.compute_forcing_terms()
2808	#print domain.quantities['stage'].explicit_update
2809
2810	assert num.allclose(domain.quantities['stage'].explicit_update[1], (3*domain.time+7)/1000)
2811	assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2812	assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2813
2814
2815	domain.time = 100.
2816	domain.quantities['stage'].explicit_update[:] = 0.0 # Reset
2817	domain.compute_forcing_terms()
2818	#print domain.quantities['stage'].explicit_update
2819
2820	assert num.allclose(domain.quantities['stage'].explicit_update[1], 5.0/1000) # Default value
2821	assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
2822	assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2823
2824
2825
2826
2827
2828
2829
2830	def test_rainfall_forcing_with_evolve(self):
2831	"""test_rainfall_forcing_with_evolve
2832
2833	Test how forcing terms are called within evolve
2834	"""
2835
2836	# FIXME(Ole): This test is just to experiment
2837
2838	a = [0.0, 0.0]
2839	b = [0.0, 2.0]
2840	c = [2.0, 0.0]
2841	d = [0.0, 4.0]
2842	e = [2.0, 2.0]
2843	f = [4.0, 0.0]
2844
2845	points = [a, b, c, d, e, f]
2846	#bac, bce, ecf, dbe
2847	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2848
2849
2850	domain = Domain(points, vertices)
2851
2852	#Flat surface with 1m of water
2853	domain.set_quantity('elevation', 0)
2854	domain.set_quantity('stage', 1.0)
2855	domain.set_quantity('friction', 0)
2856
2857	Br = Reflective_boundary(domain)
2858	domain.set_boundary({'exterior': Br})
2859
2860	# Setup only one forcing term, time dependent rainfall that expires at t==20
2861	from anuga.fit_interpolate.interpolate import Modeltime_too_late
2862	def main_rate(t):
2863	if t > 20:
2864	msg = 'Model time exceeded.'
2865	raise Modeltime_too_late, msg
2866	else:
2867	return 3*t + 7
2868
2869	domain.forcing_terms = []
2870	R = Rainfall(domain,
2871	rate=main_rate,
2872	polygon=[[1,1], [2,1], [2,2], [1,2]],
2873	default_rate=5.0)
2874
2875	assert num.allclose(R.exchange_area, 1)
2876
2877	domain.forcing_terms.append(R)
2878
2879	for t in domain.evolve(yieldstep=1, finaltime=25):
2880	pass
2881
2882	#print t, domain.quantities['stage'].explicit_update, (3*t+7)/1000
2883
2884	#FIXME(Ole): A test here is hard because explicit_update also
2885	# receives updates from the flux calculation.
2886
2887
2888
2889
2890	def test_inflow_using_circle(self):
2891	from math import pi, cos, sin
2892
2893	a = [0.0, 0.0]
2894	b = [0.0, 2.0]
2895	c = [2.0, 0.0]
2896	d = [0.0, 4.0]
2897	e = [2.0, 2.0]
2898	f = [4.0, 0.0]
2899
2900	points = [a, b, c, d, e, f]
2901	#bac, bce, ecf, dbe
2902	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2903
2904
2905	domain = Domain(points, vertices)
2906
2907	# Flat surface with 1m of water
2908	domain.set_quantity('elevation', 0)
2909	domain.set_quantity('stage', 1.0)
2910	domain.set_quantity('friction', 0)
2911
2912	Br = Reflective_boundary(domain)
2913	domain.set_boundary({'exterior': Br})
2914
2915	# Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
2916	domain.forcing_terms = []
2917	domain.forcing_terms.append( Inflow(domain, rate=2.0, center=(1,1), radius=1) )
2918
2919	domain.compute_forcing_terms()
2920	#print domain.quantities['stage'].explicit_update
2921
2922	assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
2923	assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
2924	assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2925
2926
2927	def test_inflow_using_circle_function(self):
2928	from math import pi, cos, sin
2929
2930	a = [0.0, 0.0]
2931	b = [0.0, 2.0]
2932	c = [2.0, 0.0]
2933	d = [0.0, 4.0]
2934	e = [2.0, 2.0]
2935	f = [4.0, 0.0]
2936
2937	points = [a, b, c, d, e, f]
2938	#bac, bce, ecf, dbe
2939	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2940
2941
2942	domain = Domain(points, vertices)
2943
2944	# Flat surface with 1m of water
2945	domain.set_quantity('elevation', 0)
2946	domain.set_quantity('stage', 1.0)
2947	domain.set_quantity('friction', 0)
2948
2949	Br = Reflective_boundary(domain)
2950	domain.set_boundary({'exterior': Br})
2951
2952	# Setup only one forcing term, time dependent inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
2953	domain.forcing_terms = []
2954	domain.forcing_terms.append( Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) )
2955
2956	domain.compute_forcing_terms()
2957
2958	assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
2959	assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
2960	assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
2961
2962
2963
2964
2965	def test_inflow_catch_too_few_triangles(self):
2966	"""test_inflow_catch_too_few_triangles
2967
2968	Test that exception is thrown if no triangles are covered by the inflow area
2969	"""
2970	from math import pi, cos, sin
2971
2972	a = [0.0, 0.0]
2973	b = [0.0, 2.0]
2974	c = [2.0, 0.0]
2975	d = [0.0, 4.0]
2976	e = [2.0, 2.0]
2977	f = [4.0, 0.0]
2978
2979	points = [a, b, c, d, e, f]
2980	#bac, bce, ecf, dbe
2981	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2982
2983
2984	domain = Domain(points, vertices)
2985
2986	# Flat surface with 1m of water
2987	domain.set_quantity('elevation', 0)
2988	domain.set_quantity('stage', 1.0)
2989	domain.set_quantity('friction', 0)
2990
2991	Br = Reflective_boundary(domain)
2992	domain.set_boundary({'exterior': Br})
2993
2994	# Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
2995
2996	try:
2997	Inflow(domain, rate=2.0, center=(1,1.1), radius=0.01)
2998	except:
2999	pass
3000	else:
3001	msg = 'Should have raised exception'
3002	raise Exception, msg
3003
3004
3005
3006
3007	def Xtest_inflow_outflow_conservation(self):
3008	"""test_inflow_outflow_conservation
3009
3010	Test what happens if water is abstracted from one area and
3011	injected into another - especially if there is not enough
3012	water to match the abstraction.
3013	This tests that the total volume is kept constant under a range of
3014	scenarios.
3015
3016	This test will fail as the problem was only fixed for culverts.
3017	"""
3018
3019	from math import pi, cos, sin
3020
3021	length = 20.
3022	width = 10.
3023
3024	dx = dy = 2 # 1 or 2 OK
3025	points, vertices, boundary = rectangular_cross(int(length/dx),
3026	int(width/dy),
3027	len1=length,
3028	len2=width)
3029	domain = Domain(points, vertices, boundary)
3030	domain.set_name('test_inflow_conservation') # Output name
3031	domain.set_default_order(2)
3032
3033
3034	# Flat surface with 1m of water
3035	stage = 1.0
3036	domain.set_quantity('elevation', 0)
3037	domain.set_quantity('stage', stage)
3038	domain.set_quantity('friction', 0)
3039
3040	Br = Reflective_boundary(domain)
3041	domain.set_boundary({'left': Br, 'right': Br, 'bottom': Br, 'top': Br})
3042
3043	# Setup one forcing term, constant inflow of 2 m^3/s on a circle
3044	domain.forcing_terms = []
3045	domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))
3046
3047	domain.compute_forcing_terms()
3048	#print domain.quantities['stage'].explicit_update
3049
3050	# Check that update values are correct
3051	for x in domain.quantities['stage'].explicit_update:
3052	assert num.allclose(x, 2.0/pi) or num.allclose(x, 0.0)
3053
3054
3055	# Check volumes without inflow
3056	domain.forcing_terms = []
3057	initial_volume = domain.quantities['stage'].get_integral()
3058
3059	assert num.allclose(initial_volume, widthlengthstage)
3060
3061	for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
3062	volume = domain.quantities['stage'].get_integral()
3063	assert num.allclose (volume, initial_volume)
3064
3065
3066	# Now apply the inflow and check volumes for a range of stage values
3067	for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
3068	domain.time = 0.0
3069	domain.set_quantity('stage', stage)
3070
3071	domain.forcing_terms = []
3072	domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))
3073	initial_volume = domain.quantities['stage'].get_integral()
3074	predicted_volume = initial_volume
3075	dt = 0.05
3076	for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
3077	volume = domain.quantities['stage'].get_integral()
3078
3079	assert num.allclose (volume, predicted_volume)
3080	predicted_volume = predicted_volume + 2.0/pi/100/dt # Why 100?
3081
3082
3083	# Apply equivalent outflow only and check volumes for a range of stage values
3084	for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
3085	print stage
3086
3087	domain.time = 0.0
3088	domain.set_quantity('stage', stage)
3089	domain.forcing_terms = []
3090	domain.forcing_terms.append(Inflow(domain, rate=-2.0, center=(15,5), radius=1))
3091	initial_volume = domain.quantities['stage'].get_integral()
3092	predicted_volume = initial_volume
3093	dt = 0.05
3094	for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
3095	volume = domain.quantities['stage'].get_integral()
3096
3097	print t, volume, predicted_volume
3098	assert num.allclose (volume, predicted_volume)
3099	predicted_volume = predicted_volume - 2.0/pi/100/dt # Why 100?
3100
3101
3102	# Apply both inflow and outflow and check volumes being constant for a
3103	# range of stage values
3104	for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
3105	print stage
3106
3107	domain.time = 0.0
3108	domain.set_quantity('stage', stage)
3109	domain.forcing_terms = []
3110	domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))
3111	domain.forcing_terms.append(Inflow(domain, rate=-2.0, center=(15,5), radius=1))
3112	initial_volume = domain.quantities['stage'].get_integral()
3113
3114	dt = 0.05
3115	for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
3116	volume = domain.quantities['stage'].get_integral()
3117
3118	print t, volume
3119	assert num.allclose (volume, initial_volume)
3120
3121
3122
3123
3124	#####################################################
3125	def test_first_order_extrapolator_const_z(self):
3126
3127	a = [0.0, 0.0]
3128	b = [0.0, 2.0]
3129	c = [2.0, 0.0]
3130	d = [0.0, 4.0]
3131	e = [2.0, 2.0]
3132	f = [4.0, 0.0]
3133
3134	points = [a, b, c, d, e, f]
3135	#bac, bce, ecf, dbe
3136	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3137
3138	domain = Domain(points, vertices)
3139	val0 = 2.+2.0/3
3140	val1 = 4.+4.0/3
3141	val2 = 8.+2.0/3
3142	val3 = 2.+8.0/3
3143
3144	zl=zr=-3.75 #Assume constant bed (must be less than stage)
3145	domain.set_quantity('elevation', zl*num.ones( (4,3) ))
3146	domain.set_quantity('stage', [[val0, val0-1, val0-2],
3147	[val1, val1+1, val1],
3148	[val2, val2-2, val2],
3149	[val3-0.5, val3, val3]])
3150
3151
3152
3153	domain._order_ = 1
3154	domain.distribute_to_vertices_and_edges()
3155
3156	#Check that centroid values were distributed to vertices
3157	C = domain.quantities['stage'].centroid_values
3158	for i in range(3):
3159	assert num.allclose( domain.quantities['stage'].vertex_values[:,i], C)
3160
3161
3162	def test_first_order_limiter_variable_z(self):
3163	#Check that first order limiter follows bed_slope
3164	from anuga.config import epsilon
3165
3166	a = [0.0, 0.0]
3167	b = [0.0, 2.0]
3168	c = [2.0,0.0]
3169	d = [0.0, 4.0]
3170	e = [2.0, 2.0]
3171	f = [4.0,0.0]
3172
3173	points = [a, b, c, d, e, f]
3174	#bac, bce, ecf, dbe
3175	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3176
3177	domain = Domain(points, vertices)
3178	val0 = 2.+2.0/3
3179	val1 = 4.+4.0/3
3180	val2 = 8.+2.0/3
3181	val3 = 2.+8.0/3
3182
3183	domain.set_quantity('elevation', [[0,0,0], [6,0,0],
3184	[6,6,6], [6,6,6]])
3185	domain.set_quantity('stage', [[val0, val0, val0],
3186	[val1, val1, val1],
3187	[val2, val2, val2],
3188	[val3, val3, val3]])
3189
3190	E = domain.quantities['elevation'].vertex_values
3191	L = domain.quantities['stage'].vertex_values
3192
3193
3194	#Check that some stages are not above elevation (within eps)
3195	#- so that the limiter has something to work with
3196	assert not num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon)))
3197
3198	domain._order_ = 1
3199	domain.distribute_to_vertices_and_edges()
3200
3201	#Check that all stages are above elevation (within eps)
3202	assert num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon)))
3203
3204
3205	#####################################################
3206	def test_distribute_basic(self):
3207	#Using test data generated by abstract_2d_finite_volumes-2
3208	#Assuming no friction and flat bed (0.0)
3209
3210	a = [0.0, 0.0]
3211	b = [0.0, 2.0]
3212	c = [2.0, 0.0]
3213	d = [0.0, 4.0]
3214	e = [2.0, 2.0]
3215	f = [4.0, 0.0]
3216
3217	points = [a, b, c, d, e, f]
3218	#bac, bce, ecf, dbe
3219	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3220
3221	domain = Domain(points, vertices)
3222
3223	val0 = 2.
3224	val1 = 4.
3225	val2 = 8.
3226	val3 = 2.
3227
3228	domain.set_quantity('stage', [val0, val1, val2, val3],
3229	location='centroids')
3230	L = domain.quantities['stage'].vertex_values
3231
3232	#First order
3233	domain._order_ = 1
3234	domain.distribute_to_vertices_and_edges()
3235	assert num.allclose(L[1], val1)
3236
3237	#Second order
3238	domain._order_ = 2
3239	domain.beta_w = 0.9
3240	domain.beta_w_dry = 0.9
3241	domain.beta_uh = 0.9
3242	domain.beta_uh_dry = 0.9
3243	domain.beta_vh = 0.9
3244	domain.beta_vh_dry = 0.9
3245	domain.distribute_to_vertices_and_edges()
3246	assert num.allclose(L[1], [2.2, 4.9, 4.9])
3247
3248
3249
3250	def test_distribute_away_from_bed(self):
3251	#Using test data generated by abstract_2d_finite_volumes-2
3252	#Assuming no friction and flat bed (0.0)
3253
3254	a = [0.0, 0.0]
3255	b = [0.0, 2.0]
3256	c = [2.0, 0.0]
3257	d = [0.0, 4.0]
3258	e = [2.0, 2.0]
3259	f = [4.0, 0.0]
3260
3261	points = [a, b, c, d, e, f]
3262	#bac, bce, ecf, dbe
3263	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3264
3265	domain = Domain(points, vertices)
3266	L = domain.quantities['stage'].vertex_values
3267
3268	def stage(x,y):
3269	return x**2
3270
3271	domain.set_quantity('stage', stage, location='centroids')
3272
3273	domain.quantities['stage'].compute_gradients()
3274
3275	a, b = domain.quantities['stage'].get_gradients()
3276
3277	assert num.allclose(a[1], 3.33333334)
3278	assert num.allclose(b[1], 0.0)
3279
3280	domain._order_ = 1
3281	domain.distribute_to_vertices_and_edges()
3282	assert num.allclose(L[1], 1.77777778)
3283
3284	domain._order_ = 2
3285	domain.beta_w = 0.9
3286	domain.beta_w_dry = 0.9
3287	domain.beta_uh = 0.9
3288	domain.beta_uh_dry = 0.9
3289	domain.beta_vh = 0.9
3290	domain.beta_vh_dry = 0.9
3291	domain.distribute_to_vertices_and_edges()
3292	assert num.allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
3293
3294
3295
3296	def test_distribute_away_from_bed1(self):
3297	#Using test data generated by abstract_2d_finite_volumes-2
3298	#Assuming no friction and flat bed (0.0)
3299
3300	a = [0.0, 0.0]
3301	b = [0.0, 2.0]
3302	c = [2.0, 0.0]
3303	d = [0.0, 4.0]
3304	e = [2.0, 2.0]
3305	f = [4.0, 0.0]
3306
3307	points = [a, b, c, d, e, f]
3308	#bac, bce, ecf, dbe
3309	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3310
3311	domain = Domain(points, vertices)
3312	L = domain.quantities['stage'].vertex_values
3313
3314	def stage(x,y):
3315	return x4+y2
3316
3317	domain.set_quantity('stage', stage, location='centroids')
3318	#print domain.quantities['stage'].centroid_values
3319
3320	domain.quantities['stage'].compute_gradients()
3321	a, b = domain.quantities['stage'].get_gradients()
3322	assert num.allclose(a[1], 25.18518519)
3323	assert num.allclose(b[1], 3.33333333)
3324
3325	domain._order_ = 1
3326	domain.distribute_to_vertices_and_edges()
3327	assert num.allclose(L[1], 4.9382716)
3328
3329	domain._order_ = 2
3330	domain.beta_w = 0.9
3331	domain.beta_w_dry = 0.9
3332	domain.beta_uh = 0.9
3333	domain.beta_uh_dry = 0.9
3334	domain.beta_vh = 0.9
3335	domain.beta_vh_dry = 0.9
3336	domain.distribute_to_vertices_and_edges()
3337	assert num.allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
3338
3339
3340
3341	def test_distribute_near_bed(self):
3342
3343	a = [0.0, 0.0]
3344	b = [0.0, 2.0]
3345	c = [2.0, 0.0]
3346	d = [0.0, 4.0]
3347	e = [2.0, 2.0]
3348	f = [4.0, 0.0]
3349
3350	points = [a, b, c, d, e, f]
3351	#bac, bce, ecf, dbe
3352	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3353
3354	domain = Domain(points, vertices)
3355
3356
3357	#Set up for a gradient of (10,0) at mid triangle (bce)
3358	def slope(x, y):
3359	return 10*x
3360
3361	h = 0.1
3362	def stage(x, y):
3363	return slope(x, y) + h
3364
3365	domain.set_quantity('elevation', slope)
3366	domain.set_quantity('stage', stage, location='centroids')
3367
3368	#print domain.quantities['elevation'].centroid_values
3369	#print domain.quantities['stage'].centroid_values
3370
3371	E = domain.quantities['elevation'].vertex_values
3372	L = domain.quantities['stage'].vertex_values
3373
3374	# Get reference values
3375	volumes = []
3376	for i in range(len(L)):
3377	volumes.append(num.sum(L[i])/3)
3378	assert num.allclose(volumes[i], domain.quantities['stage'].centroid_values[i])
3379
3380
3381	domain._order_ = 1
3382
3383	domain.tight_slope_limiters = 0
3384	domain.distribute_to_vertices_and_edges()
3385	assert num.allclose(L[1], [0.1, 20.1, 20.1])
3386	for i in range(len(L)):
3387	assert num.allclose(volumes[i], num.sum(L[i])/3)
3388
3389	domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
3390	domain.distribute_to_vertices_and_edges()
3391	assert num.allclose(L[1], [0.298, 20.001, 20.001])
3392	for i in range(len(L)):
3393	assert num.allclose(volumes[i], num.sum(L[i])/3)
3394
3395	domain._order_ = 2
3396
3397	domain.tight_slope_limiters = 0
3398	domain.distribute_to_vertices_and_edges()
3399	assert num.allclose(L[1], [0.1, 20.1, 20.1])
3400	for i in range(len(L)):
3401	assert num.allclose(volumes[i], num.sum(L[i])/3)
3402
3403	domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
3404	domain.distribute_to_vertices_and_edges()
3405	assert num.allclose(L[1], [0.298, 20.001, 20.001])
3406	for i in range(len(L)):
3407	assert num.allclose(volumes[i], num.sum(L[i])/3)
3408
3409
3410
3411	def test_distribute_near_bed1(self):
3412
3413	a = [0.0, 0.0]
3414	b = [0.0, 2.0]
3415	c = [2.0, 0.0]
3416	d = [0.0, 4.0]
3417	e = [2.0, 2.0]
3418	f = [4.0, 0.0]
3419
3420	points = [a, b, c, d, e, f]
3421	#bac, bce, ecf, dbe
3422	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
3423
3424	domain = Domain(points, vertices)
3425
3426
3427	#Set up for a gradient of (8,2) at mid triangle (bce)
3428	def slope(x, y):
3429	return x4+y2
3430
3431	h = 0.1
3432	def stage(x,y):
3433	return slope(x,y)+h
3434
3435	domain.set_quantity('elevation', slope)
3436	domain.set_quantity('stage', stage)
3437
3438	#print domain.quantities['elevation'].centroid_values
3439	#print domain.quantities['stage'].centroid_values
3440
3441	E = domain.quantities['elevation'].vertex_values
3442	L = domain.quantities['stage'].vertex_values
3443
3444	# Get reference values
3445	volumes = []
3446	for i in range(len(L)):
3447	volumes.append(num.sum(L[i])/3)
3448	assert num.allclose(volumes[i], domain.quantities['stage'].centroid_values[i])
3449
3450	#print E
3451	domain._order_ = 1
3452
3453	domain.tight_slope_limiters = 0
3454	domain.distribute_to_vertices_and_edges()
3455	assert num.allclose(L[1], [4.1, 16.1, 20.1])
3456	for i in range(len(L)):
3457	assert num.allclose(volumes[i], num.sum(L[i])/3)
3458
3459
3460	domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
3461	domain.distribute_to_vertices_and_edges()
3462	assert num.allclose(L[1], [4.2386, 16.0604, 20.001])
3463	for i in range(len(L)):
3464	assert num.allclose(volumes[i], num.sum(L[i])/3)
3465
3466
3467	domain._order_ = 2
3468
3469	domain.tight_slope_limiters = 0
3470	domain.distribute_to_vertices_and_edges()
3471	assert num.allclose(L[1], [4.1, 16.1, 20.1])
3472	for i in range(len(L)):
3473	assert num.allclose(volumes[i], num.sum(L[i])/3)
3474
3475	domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
3476	domain.distribute_to_vertices_and_edges()
3477	#print L[1]
3478	assert num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\
3479	num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\
3480	num.allclose(L[1], [4.19351461, 16.10548539, 20.001]) # old limiters
3481
3482	for i in range(len(L)):
3483	assert num.allclose(volumes[i], num.sum(L[i])/3)
3484
3485
3486	def test_second_order_distribute_real_data(self):
3487	#Using test data generated by abstract_2d_finite_volumes-2
3488	#Assuming no friction and flat bed (0.0)
3489
3490	a = [0.0, 0.0]
3491	b = [0.0, 1.0/5]
3492	c = [0.0, 2.0/5]
3493	d = [1.0/5, 0.0]
3494	e = [1.0/5, 1.0/5]
3495	f = [1.0/5, 2.0/5]
3496	g = [2.0/5, 2.0/5]
3497
3498	points = [a, b, c, d, e, f, g]
3499	#bae, efb, cbf, feg
3500	vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
3501
3502	domain = Domain(points, vertices)
3503
3504	def slope(x, y):
3505	return -x/3
3506
3507	domain.set_quantity('elevation', slope)
3508	domain.set_quantity('stage',
3509	[0.01298164, 0.00365611,
3510	0.01440365, -0.0381856437096],
3511	location='centroids')
3512	domain.set_quantity('xmomentum',
3513	[0.00670439, 0.01263789,
3514	0.00647805, 0.0178180740668],
3515	location='centroids')
3516	domain.set_quantity('ymomentum',
3517	[-7.23510980e-004, -6.30413883e-005,
3518	6.30413883e-005, 0.000200907255866],
3519	location='centroids')
3520
3521	E = domain.quantities['elevation'].vertex_values
3522	L = domain.quantities['stage'].vertex_values
3523	X = domain.quantities['xmomentum'].vertex_values
3524	Y = domain.quantities['ymomentum'].vertex_values
3525
3526	#print E
3527	domain._order_ = 2
3528	domain.beta_w = 0.9
3529	domain.beta_w_dry = 0.9
3530	domain.beta_uh = 0.9
3531	domain.beta_uh_dry = 0.9
3532	domain.beta_vh = 0.9
3533	domain.beta_vh_dry = 0.9
3534
3535	# FIXME (Ole): Need tests where this is commented out
3536	domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
3537	domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
3538
3539
3540	domain.distribute_to_vertices_and_edges()
3541
3542	#print L[1,:]
3543	#print X[1,:]
3544	#print Y[1,:]
3545
3546	assert num.allclose(L[1,:], [-0.00825735775384,
3547	-0.00801881482869,
3548	0.0272445025825])
3549	assert num.allclose(X[1,:], [0.0143507718962,
3550	0.0142502147066,
3551	0.00931268339717])
3552	assert num.allclose(Y[1,:], [-0.000117062180693,
3553	7.94434448109e-005,
3554	-0.000151505429018])
3555
3556
3557
3558	def test_balance_deep_and_shallow(self):
3559	"""Test that balanced limiters preserve conserved quantites.
3560	This test is using old depth based balanced limiters
3561	"""
3562	import copy
3563
3564	a = [0.0, 0.0]
3565	b = [0.0, 2.0]
3566	c = [2.0, 0.0]
3567	d = [0.0, 4.0]
3568	e = [2.0, 2.0]
3569	f = [4.0, 0.0]
3570
3571	points = [a, b, c, d, e, f]
3572
3573	#bac, bce, ecf, dbe
3574	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
3575
3576	domain = Domain(points, elements)
3577	domain.check_integrity()
3578
3579	#Create a deliberate overshoot
3580	domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3581	domain.set_quantity('elevation', 0) #Flat bed
3582	stage = domain.quantities['stage']
3583
3584	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
3585
3586	#Limit
3587	domain.tight_slope_limiters = 0
3588	domain.distribute_to_vertices_and_edges()
3589
3590	#Assert that quantities are conserved
3591	for k in range(len(domain)):
3592	assert num.allclose (ref_centroid_values[k],
3593	num.sum(stage.vertex_values[k,:])/3)
3594
3595
3596	#Now try with a non-flat bed - closely hugging initial stage in places
3597	#This will create alphas in the range [0, 0.478260, 1]
3598	domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3599	domain.set_quantity('elevation', [[0,0,0],
3600	[1.8,1.9,5.9],
3601	[4.6,0,0],
3602	[0,2,4]])
3603	stage = domain.quantities['stage']
3604
3605	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
3606	ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
3607
3608	#Limit
3609	domain.tight_slope_limiters = 0
3610	domain.distribute_to_vertices_and_edges()
3611
3612
3613	#Assert that all vertex quantities have changed
3614	for k in range(len(domain)):
3615	#print ref_vertex_values[k,:], stage.vertex_values[k,:]
3616	assert not num.allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
3617	#and assert that quantities are still conserved
3618	for k in range(len(domain)):
3619	assert num.allclose (ref_centroid_values[k],
3620	num.sum(stage.vertex_values[k,:])/3)
3621
3622
3623	# Check actual results
3624	assert num.allclose (stage.vertex_values,
3625	[[2,2,2],
3626	[1.93333333, 2.03333333, 6.03333333],
3627	[6.93333333, 4.53333333, 4.53333333],
3628	[5.33333333, 5.33333333, 5.33333333]])
3629
3630
3631	def test_balance_deep_and_shallow_tight_SL(self):
3632	"""Test that balanced limiters preserve conserved quantites.
3633	This test is using Tight Slope Limiters
3634	"""
3635	import copy
3636
3637	a = [0.0, 0.0]
3638	b = [0.0, 2.0]
3639	c = [2.0, 0.0]
3640	d = [0.0, 4.0]
3641	e = [2.0, 2.0]
3642	f = [4.0, 0.0]
3643
3644	points = [a, b, c, d, e, f]
3645
3646	#bac, bce, ecf, dbe
3647	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
3648
3649	domain = Domain(points, elements)
3650	domain.check_integrity()
3651
3652	#Create a deliberate overshoot
3653	domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3654	domain.set_quantity('elevation', 0) #Flat bed
3655	stage = domain.quantities['stage']
3656
3657	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
3658
3659	#Limit
3660	domain.tight_slope_limiters = 1
3661	domain.distribute_to_vertices_and_edges()
3662
3663	#Assert that quantities are conserved
3664	for k in range(len(domain)):
3665	assert num.allclose (ref_centroid_values[k],
3666	num.sum(stage.vertex_values[k,:])/3)
3667
3668
3669	#Now try with a non-flat bed - closely hugging initial stage in places
3670	#This will create alphas in the range [0, 0.478260, 1]
3671	domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3672	domain.set_quantity('elevation', [[0,0,0],
3673	[1.8,1.9,5.9],
3674	[4.6,0,0],
3675	[0,2,4]])
3676	stage = domain.quantities['stage']
3677
3678	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
3679	ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
3680
3681	#Limit
3682	domain.tight_slope_limiters = 1
3683	domain.distribute_to_vertices_and_edges()
3684
3685
3686	#Assert that all vertex quantities have changed
3687	for k in range(len(domain)):
3688	#print ref_vertex_values[k,:], stage.vertex_values[k,:]
3689	assert not num.allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
3690	#and assert that quantities are still conserved
3691	for k in range(len(domain)):
3692	assert num.allclose (ref_centroid_values[k],
3693	num.sum(stage.vertex_values[k,:])/3)
3694
3695
3696	#Also check that Python and C version produce the same
3697	# No longer applicable if tight_slope_limiters == 1
3698	#print stage.vertex_values
3699	#assert allclose (stage.vertex_values,
3700	# [[2,2,2],
3701	# [1.93333333, 2.03333333, 6.03333333],
3702	# [6.93333333, 4.53333333, 4.53333333],
3703	# [5.33333333, 5.33333333, 5.33333333]])
3704
3705
3706
3707	def test_balance_deep_and_shallow_Froude(self):
3708	"""Test that balanced limiters preserve conserved quantites -
3709	and also that excessive Froude numbers are dealt with.
3710	This test is using tight slope limiters.
3711	"""
3712	import copy
3713
3714	a = [0.0, 0.0]
3715	b = [0.0, 2.0]
3716	c = [2.0, 0.0]
3717	d = [0.0, 4.0]
3718	e = [2.0, 2.0]
3719	f = [4.0, 0.0]
3720
3721	points = [a, b, c, d, e, f]
3722
3723	# bac, bce, ecf, dbe
3724	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
3725
3726	domain = Domain(points, elements)
3727	domain.check_integrity()
3728	domain.tight_slope_limiters = True
3729	domain.use_centroid_velocities = True
3730
3731	# Create non-flat bed - closely hugging initial stage in places
3732	# This will create alphas in the range [0, 0.478260, 1]
3733	domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
3734	domain.set_quantity('elevation', [[0,0,0],
3735	[1.8,1.999,5.999],
3736	[4.6,0,0],
3737	[0,2,4]])
3738
3739	# Create small momenta, that nonetheless will generate large speeds
3740	# due to shallow depth at isolated vertices
3741	domain.set_quantity('xmomentum', -0.0058)
3742	domain.set_quantity('ymomentum', 0.0890)
3743
3744
3745
3746
3747	stage = domain.quantities['stage']
3748	elevation = domain.quantities['elevation']
3749	xmomentum = domain.quantities['xmomentum']
3750	ymomentum = domain.quantities['ymomentum']
3751
3752	# Setup triangle #1 to mimick real Froude explosion observed
3753	# in the Onslow example 13 Nov 2007.
3754
3755	stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953]
3756	elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647]
3757	xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066]
3758	ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890]
3759
3760	xmomentum.interpolate()
3761	ymomentum.interpolate()
3762	stage.interpolate()
3763	elevation.interpolate()
3764
3765	# Verify interpolation
3766	assert num.allclose(stage.centroid_values[1], 1.5233)
3767	assert num.allclose(elevation.centroid_values[1], 1.2452667)
3768	assert num.allclose(xmomentum.centroid_values[1], -0.0058)
3769	assert num.allclose(ymomentum.centroid_values[1], 0.089)
3770
3771	# Derived quantities
3772	depth = stage-elevation
3773	u = xmomentum/depth
3774	v = ymomentum/depth
3775
3776	denom = (depthg)*0.5
3777	Fx = u/denom
3778	Fy = v/denom
3779
3780
3781	# Verify against Onslow example (14 Nov 2007)
3782	assert num.allclose(depth.centroid_values[1], 0.278033)
3783	assert num.allclose(u.centroid_values[1], -0.0208608)
3784	assert num.allclose(v.centroid_values[1], 0.3201055)
3785
3786	assert num.allclose(denom.centroid_values[1],
3787	num.sqrt(depth.centroid_values[1]*g))
3788
3789	assert num.allclose(u.centroid_values[1]/denom.centroid_values[1],
3790	-0.012637746977)
3791	assert num.allclose(Fx.centroid_values[1],
3792	u.centroid_values[1]/denom.centroid_values[1])
3793
3794	# Check that Froude numbers are small at centroids.
3795	assert num.allclose(Fx.centroid_values[1], -0.012637746977)
3796	assert num.allclose(Fy.centroid_values[1], 0.193924048435)
3797
3798
3799	# But Froude numbers are huge at some vertices and edges
3800	assert num.allclose(Fx.vertex_values[1,:], [-5.85888475e+01,
3801	-1.27775313e+01,
3802	-2.78511420e-03])
3803
3804	assert num.allclose(Fx.edge_values[1,:], [-6.89150773e-03,
3805	-7.38672488e-03,
3806	-2.35626238e+01])
3807
3808	assert num.allclose(Fy.vertex_values[1,:], [8.99035764e+02,
3809	2.27440057e+02,
3810	3.75568430e-02])
3811
3812	assert num.allclose(Fy.edge_values[1,:], [1.05748998e-01,
3813	1.06035244e-01,
3814	3.88346947e+02])
3815
3816
3817	# The task is now to arrange the limiters such that Froude numbers
3818	# remain under control whil at the same time obeying the conservation
3819	# laws.
3820
3821
3822	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
3823	ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
3824
3825	# Limit (and invoke balance_deep_and_shallow)
3826	domain.tight_slope_limiters = 1
3827	domain.distribute_to_vertices_and_edges()
3828
3829	# Redo derived quantities
3830	depth = stage-elevation
3831	u = xmomentum/depth
3832	v = ymomentum/depth
3833
3834	# Assert that all vertex velocities stay within one
3835	# order of magnitude of centroid velocities.
3836	#print u.vertex_values[1,:]
3837	#print u.centroid_values[1]
3838
3839	assert num.alltrue(num.absolute(u.vertex_values[1,:]) <= num.absolute(u.centroid_values[1])*10)
3840	assert num.alltrue(num.absolute(v.vertex_values[1,:]) <= num.absolute(v.centroid_values[1])*10)
3841
3842	denom = (depthg)*0.5
3843	Fx = u/denom
3844	Fy = v/denom
3845
3846
3847	# Assert that Froude numbers are less than max value (TBA)
3848	# at vertices, edges and centroids.
3849	from anuga.config import maximum_froude_number
3850	assert num.alltrue(num.absolute(Fx.vertex_values[1,:]) < maximum_froude_number)
3851	assert num.alltrue(num.absolute(Fy.vertex_values[1,:]) < maximum_froude_number)
3852
3853
3854	# Assert that all vertex quantities have changed
3855	for k in range(len(domain)):
3856	#print ref_vertex_values[k,:], stage.vertex_values[k,:]
3857	assert not num.allclose (ref_vertex_values[k,:],
3858	stage.vertex_values[k,:])
3859
3860	# Assert that quantities are still conserved
3861	for k in range(len(domain)):
3862	assert num.allclose (ref_centroid_values[k],
3863	num.sum(stage.vertex_values[k,:])/3)
3864
3865
3866
3867	return
3868
3869	qwidth = 12
3870	for k in [1]: #range(len(domain)):
3871	print 'Triangle %d (C, V, E)' %k
3872
3873	print 'stage'.ljust(qwidth), stage.centroid_values[k],\
3874	stage.vertex_values[k,:], stage.edge_values[k,:]
3875	print 'elevation'.ljust(qwidth), elevation.centroid_values[k],\
3876	elevation.vertex_values[k,:], elevation.edge_values[k,:]
3877	print 'depth'.ljust(qwidth), depth.centroid_values[k],\
3878	depth.vertex_values[k,:], depth.edge_values[k,:]
3879	print 'xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],\
3880	xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:]
3881	print 'ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],\
3882	ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:]
3883	print 'u'.ljust(qwidth),u.centroid_values[k],\
3884	u.vertex_values[k,:], u.edge_values[k,:]
3885	print 'v'.ljust(qwidth), v.centroid_values[k],\
3886	v.vertex_values[k,:], v.edge_values[k,:]
3887	print 'Fx'.ljust(qwidth), Fx.centroid_values[k],\
3888	Fx.vertex_values[k,:], Fx.edge_values[k,:]
3889	print 'Fy'.ljust(qwidth), Fy.centroid_values[k],\
3890	Fy.vertex_values[k,:], Fy.edge_values[k,:]
3891
3892
3893
3894
3895
3896
3897	def test_conservation_1(self):
3898	"""Test that stage is conserved globally
3899
3900	This one uses a flat bed, reflective bdries and a suitable
3901	initial condition
3902	"""
3903	from mesh_factory import rectangular
3904
3905	#Create basic mesh
3906	points, vertices, boundary = rectangular(6, 6)
3907
3908	#Create shallow water domain
3909	domain = Domain(points, vertices, boundary)
3910	domain.smooth = False
3911	domain.default_order=2
3912
3913	#IC
3914	def x_slope(x, y):
3915	return x/3
3916
3917	domain.set_quantity('elevation', 0)
3918	domain.set_quantity('friction', 0)
3919	domain.set_quantity('stage', x_slope)
3920
3921	# Boundary conditions (reflective everywhere)
3922	Br = Reflective_boundary(domain)
3923	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3924
3925	domain.check_integrity()
3926
3927	initial_volume = domain.quantities['stage'].get_integral()
3928	initial_xmom = domain.quantities['xmomentum'].get_integral()
3929
3930	#print initial_xmom
3931
3932	#Evolution
3933	for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
3934	volume = domain.quantities['stage'].get_integral()
3935	assert num.allclose (volume, initial_volume)
3936
3937	#I don't believe that the total momentum should be the same
3938	#It starts with zero and ends with zero though
3939	#xmom = domain.quantities['xmomentum'].get_integral()
3940	#print xmom
3941	#assert allclose (xmom, initial_xmom)
3942
3943	os.remove(domain.get_name() + '.sww')
3944
3945
3946	def test_conservation_2(self):
3947	"""Test that stage is conserved globally
3948
3949	This one uses a slopy bed, reflective bdries and a suitable
3950	initial condition
3951	"""
3952	from mesh_factory import rectangular
3953
3954	#Create basic mesh
3955	points, vertices, boundary = rectangular(6, 6)
3956
3957	#Create shallow water domain
3958	domain = Domain(points, vertices, boundary)
3959	domain.smooth = False
3960	domain.default_order=2
3961
3962	#IC
3963	def x_slope(x, y):
3964	return x/3
3965
3966	domain.set_quantity('elevation', x_slope)
3967	domain.set_quantity('friction', 0)
3968	domain.set_quantity('stage', 0.4) #Steady
3969
3970	# Boundary conditions (reflective everywhere)
3971	Br = Reflective_boundary(domain)
3972	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
3973
3974	domain.check_integrity()
3975
3976	initial_volume = domain.quantities['stage'].get_integral()
3977	initial_xmom = domain.quantities['xmomentum'].get_integral()
3978
3979	#print initial_xmom
3980
3981	#Evolution
3982	for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
3983	volume = domain.quantities['stage'].get_integral()
3984	assert num.allclose (volume, initial_volume)
3985
3986	#FIXME: What would we expect from momentum
3987	#xmom = domain.quantities['xmomentum'].get_integral()
3988	#print xmom
3989	#assert allclose (xmom, initial_xmom)
3990
3991	os.remove(domain.get_name() + '.sww')
3992
3993	def test_conservation_3(self):
3994	"""Test that stage is conserved globally
3995
3996	This one uses a larger grid, convoluted bed, reflective bdries and a suitable
3997	initial condition
3998	"""
3999	from mesh_factory import rectangular
4000
4001	#Create basic mesh
4002	points, vertices, boundary = rectangular(2, 1)
4003
4004	#Create shallow water domain
4005	domain = Domain(points, vertices, boundary)
4006	domain.smooth = False
4007	domain.default_order = 2
4008	domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
4009
4010	#IC
4011	def x_slope(x, y):
4012	z = 0*x
4013	for i in range(len(x)):
4014	if x[i] < 0.3:
4015	z[i] = x[i]/3
4016	if 0.3 <= x[i] < 0.5:
4017	z[i] = -0.5
4018	if 0.5 <= x[i] < 0.7:
4019	z[i] = 0.39
4020	if 0.7 <= x[i]:
4021	z[i] = x[i]/3
4022	return z
4023
4024
4025
4026	domain.set_quantity('elevation', x_slope)
4027	domain.set_quantity('friction', 0)
4028	domain.set_quantity('stage', 0.4) #Steady
4029
4030	# Boundary conditions (reflective everywhere)
4031	Br = Reflective_boundary(domain)
4032	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4033
4034	domain.check_integrity()
4035
4036	initial_volume = domain.quantities['stage'].get_integral()
4037	initial_xmom = domain.quantities['xmomentum'].get_integral()
4038
4039	import copy
4040	ref_centroid_values =\
4041	copy.copy(domain.quantities['stage'].centroid_values)
4042
4043	#print 'ORG', domain.quantities['stage'].centroid_values
4044	domain.distribute_to_vertices_and_edges()
4045
4046
4047	#print domain.quantities['stage'].centroid_values
4048	assert num.allclose(domain.quantities['stage'].centroid_values,
4049	ref_centroid_values)
4050
4051
4052	#Check that initial limiter doesn't violate cons quan
4053	assert num.allclose(domain.quantities['stage'].get_integral(),
4054	initial_volume)
4055
4056	#Evolution
4057	for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
4058	volume = domain.quantities['stage'].get_integral()
4059	#print t, volume, initial_volume
4060	assert num.allclose (volume, initial_volume)
4061
4062	os.remove(domain.get_name() + '.sww')
4063
4064	def test_conservation_4(self):
4065	"""Test that stage is conserved globally
4066
4067	This one uses a larger grid, convoluted bed, reflective bdries and a suitable
4068	initial condition
4069	"""
4070	from mesh_factory import rectangular
4071
4072	#Create basic mesh
4073	points, vertices, boundary = rectangular(6, 6)
4074
4075	#Create shallow water domain
4076	domain = Domain(points, vertices, boundary)
4077	domain.smooth = False
4078	domain.default_order=2
4079
4080	#IC
4081	def x_slope(x, y):
4082	z = 0*x
4083	for i in range(len(x)):
4084	if x[i] < 0.3:
4085	z[i] = x[i]/3
4086	if 0.3 <= x[i] < 0.5:
4087	z[i] = -0.5
4088	if 0.5 <= x[i] < 0.7:
4089	#z[i] = 0.3 #OK with beta == 0.2
4090	z[i] = 0.34 #OK with beta == 0.0
4091	#z[i] = 0.35#Fails after 80 timesteps with an error
4092	#of the order 1.0e-5
4093	if 0.7 <= x[i]:
4094	z[i] = x[i]/3
4095	return z
4096
4097
4098
4099	domain.set_quantity('elevation', x_slope)
4100	domain.set_quantity('friction', 0)
4101	domain.set_quantity('stage', 0.4) #Steady
4102
4103	# Boundary conditions (reflective everywhere)
4104	Br = Reflective_boundary(domain)
4105	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4106
4107	domain.check_integrity()
4108
4109	initial_volume = domain.quantities['stage'].get_integral()
4110	initial_xmom = domain.quantities['xmomentum'].get_integral()
4111
4112	import copy
4113	ref_centroid_values =\
4114	copy.copy(domain.quantities['stage'].centroid_values)
4115
4116	#Test limiter by itself
4117	domain.distribute_to_vertices_and_edges()
4118
4119	#Check that initial limiter doesn't violate cons quan
4120	assert num.allclose (domain.quantities['stage'].get_integral(),
4121	initial_volume)
4122	#NOTE: This would fail if any initial stage was less than the
4123	#corresponding bed elevation - but that is reasonable.
4124
4125
4126	#Evolution
4127	for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
4128	volume = domain.quantities['stage'].get_integral()
4129
4130	#print t, volume, initial_volume
4131
4132	assert num.allclose (volume, initial_volume)
4133
4134
4135	os.remove(domain.get_name() + '.sww')
4136
4137
4138	def test_conservation_5(self):
4139	"""Test that momentum is conserved globally in
4140	steady state scenario
4141
4142	This one uses a slopy bed, dirichlet and reflective bdries
4143	"""
4144	from mesh_factory import rectangular
4145
4146	# Create basic mesh
4147	points, vertices, boundary = rectangular(6, 6)
4148
4149	# Create shallow water domain
4150	domain = Domain(points, vertices, boundary)
4151	domain.smooth = False
4152	domain.default_order = 2
4153
4154	# IC
4155	def x_slope(x, y):
4156	return x/3
4157
4158	domain.set_quantity('elevation', x_slope)
4159	domain.set_quantity('friction', 0)
4160	domain.set_quantity('stage', 0.4) # Steady
4161
4162	# Boundary conditions (reflective everywhere)
4163	Br = Reflective_boundary(domain)
4164	Bleft = Dirichlet_boundary([0.5,0,0])
4165	Bright = Dirichlet_boundary([0.1,0,0])
4166	domain.set_boundary({'left': Bleft, 'right': Bright,
4167	'top': Br, 'bottom': Br})
4168
4169	domain.check_integrity()
4170
4171	initial_volume = domain.quantities['stage'].get_integral()
4172	initial_xmom = domain.quantities['xmomentum'].get_integral()
4173
4174
4175	# Evolution
4176	for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
4177	stage = domain.quantities['stage'].get_integral()
4178	xmom = domain.quantities['xmomentum'].get_integral()
4179	ymom = domain.quantities['ymomentum'].get_integral()
4180
4181	if num.allclose(t, 6): # Steady state reached
4182	steady_xmom = domain.quantities['xmomentum'].get_integral()
4183	steady_ymom = domain.quantities['ymomentum'].get_integral()
4184	steady_stage = domain.quantities['stage'].get_integral()
4185
4186	if t > 6:
4187	#print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
4188	msg = 'xmom=%.2f, steady_xmom=%.2f' %(xmom, steady_xmom)
4189	assert num.allclose(xmom, steady_xmom), msg
4190	assert num.allclose(ymom, steady_ymom)
4191	assert num.allclose(stage, steady_stage)
4192
4193
4194	os.remove(domain.get_name() + '.sww')
4195
4196
4197
4198
4199
4200	def test_conservation_real(self):
4201	"""Test that momentum is conserved globally
4202
4203	Stephen finally made a test that revealed the problem.
4204	This test failed with code prior to 25 July 2005
4205	"""
4206
4207	yieldstep = 0.01
4208	finaltime = 0.05
4209	min_depth = 1.0e-2
4210
4211
4212	import sys
4213	from os import sep; sys.path.append('..'+sep+'abstract_2d_finite_volumes')
4214	from mesh_factory import rectangular
4215
4216
4217	#Create shallow water domain
4218	points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
4219	domain = Domain(points, vertices, boundary)
4220	domain.smooth = False
4221	domain.default_order = 1
4222	domain.minimum_allowed_height = min_depth
4223
4224	# Set initial condition
4225	class Set_IC:
4226	"""Set an initial condition with a constant value, for x0<x<x1
4227	"""
4228
4229	def __init__(self, x0=0.25, x1=0.5, h=1.0):
4230	self.x0 = x0
4231	self.x1 = x1
4232	self.h = h
4233
4234	def __call__(self, x, y):
4235	return self.h*((x>self.x0)&(x<self.x1))
4236
4237
4238	domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
4239
4240
4241	#Boundaries
4242	R = Reflective_boundary(domain)
4243	domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R})
4244
4245	ref = domain.quantities['stage'].get_integral()
4246
4247	# Evolution
4248	for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
4249	pass
4250	#print 'Integral stage = ',\
4251	# domain.quantities['stage'].get_integral(),\
4252	# ' Time = ',domain.time
4253
4254
4255	now = domain.quantities['stage'].get_integral()
4256
4257	msg = 'Stage not conserved: was %f, now %f' %(ref, now)
4258	assert num.allclose(ref, now), msg
4259
4260	os.remove(domain.get_name() + '.sww')
4261
4262	def test_second_order_flat_bed_onestep(self):
4263
4264	from mesh_factory import rectangular
4265
4266	#Create basic mesh
4267	points, vertices, boundary = rectangular(6, 6)
4268
4269	#Create shallow water domain
4270	domain = Domain(points, vertices, boundary)
4271	domain.smooth = False
4272	domain.default_order = 2
4273	domain.beta_w = 0.9
4274	domain.beta_w_dry = 0.9
4275	domain.beta_uh = 0.9
4276	domain.beta_uh_dry = 0.9
4277	domain.beta_vh = 0.9
4278	domain.beta_vh_dry = 0.9
4279	domain.H0 = 0 # Backwards compatibility (6/2/7)
4280
4281	# Boundary conditions
4282	Br = Reflective_boundary(domain)
4283	Bd = Dirichlet_boundary([0.1, 0., 0.])
4284	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4285
4286	domain.check_integrity()
4287
4288	# Evolution
4289	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
4290	pass# domain.write_time()
4291
4292	# Data from earlier version of abstract_2d_finite_volumes
4293	assert num.allclose(domain.min_timestep, 0.0396825396825)
4294	assert num.allclose(domain.max_timestep, 0.0396825396825)
4295
4296	assert num.allclose(domain.quantities['stage'].centroid_values[:12],
4297	[0.00171396, 0.02656103, 0.00241523, 0.02656103,
4298	0.00241523, 0.02656103, 0.00241523, 0.02656103,
4299	0.00241523, 0.02656103, 0.00241523, 0.0272623])
4300
4301	domain.distribute_to_vertices_and_edges()
4302
4303	assert num.allclose(domain.quantities['stage'].vertex_values[:12,0],
4304	[0.0001714, 0.02656103, 0.00024152,
4305	0.02656103, 0.00024152, 0.02656103,
4306	0.00024152, 0.02656103, 0.00024152,
4307	0.02656103, 0.00024152, 0.0272623])
4308
4309	assert num.allclose(domain.quantities['stage'].vertex_values[:12,1],
4310	[0.00315012, 0.02656103, 0.00024152, 0.02656103,
4311	0.00024152, 0.02656103, 0.00024152, 0.02656103,
4312	0.00024152, 0.02656103, 0.00040506, 0.0272623])
4313
4314	assert num.allclose(domain.quantities['stage'].vertex_values[:12,2],
4315	[0.00182037, 0.02656103, 0.00676264,
4316	0.02656103, 0.00676264, 0.02656103,
4317	0.00676264, 0.02656103, 0.00676264,
4318	0.02656103, 0.0065991, 0.0272623])
4319
4320	assert num.allclose(domain.quantities['xmomentum'].centroid_values[:12],
4321	[0.00113961, 0.01302432, 0.00148672,
4322	0.01302432, 0.00148672, 0.01302432,
4323	0.00148672, 0.01302432, 0.00148672 ,
4324	0.01302432, 0.00148672, 0.01337143])
4325
4326	assert num.allclose(domain.quantities['ymomentum'].centroid_values[:12],
4327	[-2.91240050e-004 , 1.22721531e-004,
4328	-1.22721531e-004, 1.22721531e-004 ,
4329	-1.22721531e-004, 1.22721531e-004,
4330	-1.22721531e-004 , 1.22721531e-004,
4331	-1.22721531e-004, 1.22721531e-004,
4332	-1.22721531e-004, -4.57969873e-005])
4333
4334	os.remove(domain.get_name() + '.sww')
4335
4336
4337	def test_second_order_flat_bed_moresteps(self):
4338
4339	from mesh_factory import rectangular
4340
4341	#Create basic mesh
4342	points, vertices, boundary = rectangular(6, 6)
4343
4344	#Create shallow water domain
4345	domain = Domain(points, vertices, boundary)
4346	domain.smooth = False
4347	domain.default_order=2
4348
4349	# Boundary conditions
4350	Br = Reflective_boundary(domain)
4351	Bd = Dirichlet_boundary([0.1, 0., 0.])
4352	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4353
4354	domain.check_integrity()
4355
4356	#Evolution
4357	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
4358	pass
4359
4360	#Data from earlier version of abstract_2d_finite_volumes
4361	#assert allclose(domain.min_timestep, 0.0396825396825)
4362	#assert allclose(domain.max_timestep, 0.0396825396825)
4363	#print domain.quantities['stage'].centroid_values
4364
4365	os.remove(domain.get_name() + '.sww')
4366
4367
4368	def test_flatbed_first_order(self):
4369	from mesh_factory import rectangular
4370
4371	#Create basic mesh
4372	N = 8
4373	points, vertices, boundary = rectangular(N, N)
4374
4375	#Create shallow water domain
4376	domain = Domain(points, vertices, boundary)
4377	domain.smooth = False
4378	domain.default_order=1
4379	domain.H0 = 0 # Backwards compatibility (6/2/7)
4380
4381	# Boundary conditions
4382	Br = Reflective_boundary(domain)
4383	Bd = Dirichlet_boundary([0.2,0.,0.])
4384
4385	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4386	domain.check_integrity()
4387
4388
4389	#Evolution
4390	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
4391	pass
4392	#domain.write_time()
4393
4394	#FIXME: These numbers were from version before 25/10
4395	#assert allclose(domain.min_timestep, 0.0140413643926)
4396	#assert allclose(domain.max_timestep, 0.0140947355753)
4397
4398	for i in range(3):
4399	#assert allclose(domain.quantities['stage'].edge_values[:4,i],
4400	# [0.10730244,0.12337617,0.11200126,0.12605666])
4401
4402	assert num.allclose(domain.quantities['xmomentum'].edge_values[:4,i],
4403	[0.07610894,0.06901572,0.07284461,0.06819712])
4404
4405	#assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
4406	# [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
4407
4408
4409	os.remove(domain.get_name() + '.sww')
4410
4411	def test_flatbed_second_order(self):
4412	from mesh_factory import rectangular
4413
4414	#Create basic mesh
4415	N = 8
4416	points, vertices, boundary = rectangular(N, N)
4417
4418	#Create shallow water domain
4419	domain = Domain(points, vertices, boundary)
4420	domain.smooth = False
4421	domain.default_order=2
4422	domain.beta_w = 0.9
4423	domain.beta_w_dry = 0.9
4424	domain.beta_uh = 0.9
4425	domain.beta_uh_dry = 0.9
4426	domain.beta_vh = 0.9
4427	domain.beta_vh_dry = 0.9
4428	#domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
4429	domain.H0 = 0 # Backwards compatibility (6/2/7)
4430	domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
4431	domain.set_maximum_allowed_speed(1.0)
4432
4433	# Boundary conditions
4434	Br = Reflective_boundary(domain)
4435	Bd = Dirichlet_boundary([0.2,0.,0.])
4436
4437	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4438	domain.check_integrity()
4439
4440	# Evolution
4441	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
4442	pass
4443
4444	msg = 'min step was %f instead of %f' %(domain.min_timestep,
4445	0.0210448446782)
4446
4447	assert num.allclose(domain.min_timestep, 0.0210448446782), msg
4448	assert num.allclose(domain.max_timestep, 0.0210448446782)
4449
4450	#print domain.quantities['stage'].vertex_values[:4,0]
4451	#print domain.quantities['xmomentum'].vertex_values[:4,0]
4452	#print domain.quantities['ymomentum'].vertex_values[:4,0]
4453
4454	#FIXME: These numbers were from version before 25/10
4455	#assert allclose(domain.quantities['stage'].vertex_values[:4,0],
4456	# [0.00101913,0.05352143,0.00104852,0.05354394])
4457
4458	#FIXME: These numbers were from version before 21/3/6 -
4459	#could be recreated by setting maximum_allowed_speed to 0 maybe
4460	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4461	# [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
4462
4463	assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4464	[ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
4465
4466
4467
4468	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4469	# [0.00090581,0.03685719,0.00088303,0.03687313])
4470
4471	assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
4472	[-0.00139463,0.0006156,-0.00060364,0.00061827])
4473
4474
4475	os.remove(domain.get_name() + '.sww')
4476
4477
4478	def test_flatbed_second_order_vmax_0(self):
4479	from mesh_factory import rectangular
4480
4481	#Create basic mesh
4482	N = 8
4483	points, vertices, boundary = rectangular(N, N)
4484
4485	#Create shallow water domain
4486	domain = Domain(points, vertices, boundary)
4487	domain.smooth = False
4488	domain.default_order=2
4489	domain.beta_w = 0.9
4490	domain.beta_w_dry = 0.9
4491	domain.beta_uh = 0.9
4492	domain.beta_uh_dry = 0.9
4493	domain.beta_vh = 0.9
4494	domain.beta_vh_dry = 0.9
4495	domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
4496	domain.H0 = 0 # Backwards compatibility (6/2/7)
4497	domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
4498
4499	# Boundary conditions
4500	Br = Reflective_boundary(domain)
4501	Bd = Dirichlet_boundary([0.2,0.,0.])
4502
4503	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4504	domain.check_integrity()
4505
4506	#Evolution
4507	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
4508	pass
4509
4510
4511	assert num.allclose(domain.min_timestep, 0.0210448446782)
4512	assert num.allclose(domain.max_timestep, 0.0210448446782)
4513
4514	#FIXME: These numbers were from version before 21/3/6 -
4515	#could be recreated by setting maximum_allowed_speed to 0 maybe
4516	assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4517	[ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
4518
4519
4520	assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
4521	[-0.00139463,0.0006156,-0.00060364,0.00061827])
4522
4523
4524	os.remove(domain.get_name() + '.sww')
4525
4526
4527
4528	def test_flatbed_second_order_distribute(self):
4529	#Use real data from anuga.abstract_2d_finite_volumes 2
4530	#painfully setup and extracted.
4531	from mesh_factory import rectangular
4532
4533	#Create basic mesh
4534	N = 8
4535	points, vertices, boundary = rectangular(N, N)
4536
4537	#Create shallow water domain
4538	domain = Domain(points, vertices, boundary)
4539	domain.smooth = False
4540	domain.default_order=domain._order_=2
4541	domain.beta_w = 0.9
4542	domain.beta_w_dry = 0.9
4543	domain.beta_uh = 0.9
4544	domain.beta_uh_dry = 0.9
4545	domain.beta_vh = 0.9
4546	domain.beta_vh_dry = 0.9
4547	domain.H0 = 0 # Backwards compatibility (6/2/7)
4548	domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
4549	domain.set_maximum_allowed_speed(1.0)
4550
4551	# Boundary conditions
4552	Br = Reflective_boundary(domain)
4553	Bd = Dirichlet_boundary([0.2,0.,0.])
4554
4555	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
4556	domain.check_integrity()
4557
4558
4559
4560	for V in [False, True]:
4561	if V:
4562	#Set centroids as if system had been evolved
4563	L = num.zeros(2NN, num.Float)
4564	L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
4565	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
4566	1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
4567	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
4568	1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
4569	5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
4570	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
4571	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
4572	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
4573	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
4574	0.00000000e+000, 5.57305948e-005]
4575
4576	X = num.zeros(2NN, num.Float)
4577	X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
4578	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
4579	8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
4580	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
4581	8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
4582	3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
4583	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
4584	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
4585	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
4586	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
4587	0.00000000e+000, 4.57662812e-005]
4588
4589	Y = num.zeros(2NN, num.Float)
4590	Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
4591	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
4592	-6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
4593	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
4594	-6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
4595	-1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
4596	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
4597	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
4598	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
4599	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
4600	0.00000000e+000, -2.57635178e-005]
4601
4602
4603	domain.set_quantity('stage', L, location='centroids')
4604	domain.set_quantity('xmomentum', X, location='centroids')
4605	domain.set_quantity('ymomentum', Y, location='centroids')
4606
4607	domain.check_integrity()
4608	else:
4609	#Evolution
4610	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
4611	pass
4612	assert num.allclose(domain.min_timestep, 0.0210448446782)
4613	assert num.allclose(domain.max_timestep, 0.0210448446782)
4614
4615
4616	#Centroids were correct but not vertices.
4617	#Hence the check of distribute below.
4618	assert num.allclose(domain.quantities['stage'].centroid_values[:4],
4619	[0.00721206,0.05352143,0.01009108,0.05354394])
4620
4621	assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
4622	[0.00648352,0.03685719,0.00850733,0.03687313])
4623
4624	assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
4625	[-0.00139463,0.0006156,-0.00060364,0.00061827])
4626
4627	#print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
4628	#print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
4629
4630	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
4631	##print domain.quantities['xmomentum'].centroid_values[17], V
4632	##print
4633	if not V:
4634	#FIXME: These numbers were from version before 21/3/6 -
4635	#could be recreated by setting maximum_allowed_speed to 0 maybe
4636
4637	#assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
4638	assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)
4639
4640	else:
4641	assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
4642
4643	import copy
4644	XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
4645	assert num.allclose(domain.quantities['xmomentum'].centroid_values, XX)
4646
4647	domain.distribute_to_vertices_and_edges()
4648
4649	#assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
4650
4651	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],
4652	# 0.0)
4653	assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)
4654
4655
4656	#FIXME: These numbers were from version before 25/10
4657	#assert allclose(domain.quantities['stage'].vertex_values[:4,0],
4658	# [0.00101913,0.05352143,0.00104852,0.05354394])
4659
4660	assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
4661	[-0.00139463,0.0006156,-0.00060364,0.00061827])
4662
4663
4664	assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4665	[0.00090581,0.03685719,0.00088303,0.03687313])
4666
4667
4668	#NB NO longer relvant:
4669
4670	#This was the culprit. First triangles vertex 0 had an
4671	#x-momentum of 0.0064835 instead of 0.00090581 and
4672	#third triangle had 0.00850733 instead of 0.00088303
4673	#print domain.quantities['xmomentum'].vertex_values[:4,0]
4674
4675	#print domain.quantities['xmomentum'].vertex_values[:4,0]
4676	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
4677	# [0.00090581,0.03685719,0.00088303,0.03687313])
4678
4679	os.remove(domain.get_name() + '.sww')
4680
4681
4682
4683	def test_bedslope_problem_first_order(self):
4684
4685	from mesh_factory import rectangular
4686
4687	#Create basic mesh
4688	points, vertices, boundary = rectangular(6, 6)
4689
4690	#Create shallow water domain
4691	domain = Domain(points, vertices, boundary)
4692	domain.smooth = False
4693	domain.default_order = 1
4694
4695	#Bed-slope and friction
4696	def x_slope(x, y):
4697	return -x/3
4698
4699	domain.set_quantity('elevation', x_slope)
4700
4701	# Boundary conditions
4702	Br = Reflective_boundary(domain)
4703	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4704
4705	#Initial condition
4706	#domain.set_quantity('stage', Constant_height(x_slope, 0.05))
4707	domain.set_quantity('stage', expression='elevation+0.05')
4708	domain.check_integrity()
4709
4710	#Evolution
4711	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
4712	pass# domain.write_time()
4713
4714	# FIXME (Ole): Need some other assertion here!
4715	#print domain.min_timestep, domain.max_timestep
4716	#assert allclose(domain.min_timestep, 0.050010003001)
4717	#assert allclose(domain.max_timestep, 0.050010003001)
4718
4719
4720	os.remove(domain.get_name() + '.sww')
4721
4722	def test_bedslope_problem_first_order_moresteps(self):
4723
4724	from mesh_factory import rectangular
4725
4726	#Create basic mesh
4727	points, vertices, boundary = rectangular(6, 6)
4728
4729	#Create shallow water domain
4730	domain = Domain(points, vertices, boundary)
4731	domain.smooth = False
4732	domain.default_order = 1
4733
4734	# FIXME (Ole): Need tests where these two are commented out
4735	domain.H0 = 0 # Backwards compatibility (6/2/7)
4736	domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
4737	domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
4738
4739	#Bed-slope and friction
4740	def x_slope(x, y):
4741	return -x/3
4742
4743	domain.set_quantity('elevation', x_slope)
4744
4745	# Boundary conditions
4746	Br = Reflective_boundary(domain)
4747	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4748
4749	#Initial condition
4750	domain.set_quantity('stage', expression='elevation+0.05')
4751	domain.check_integrity()
4752
4753	#Evolution
4754	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
4755	pass# domain.write_time()
4756
4757	#Data from earlier version of abstract_2d_finite_volumes
4758	#print domain.quantities['stage'].centroid_values
4759
4760	assert num.allclose(domain.quantities['stage'].centroid_values,
4761	[-0.02998628, -0.01520652, -0.03043492,
4762	-0.0149132, -0.03004706, -0.01476251,
4763	-0.0298215, -0.01467976, -0.02988158,
4764	-0.01474662, -0.03036161, -0.01442995,
4765	-0.07624583, -0.06297061, -0.07733792,
4766	-0.06342237, -0.07695439, -0.06289595,
4767	-0.07635559, -0.0626065, -0.07633628,
4768	-0.06280072, -0.07739632, -0.06386738,
4769	-0.12161738, -0.11028239, -0.1223796,
4770	-0.11095953, -0.12189744, -0.11048616,
4771	-0.12074535, -0.10987605, -0.12014311,
4772	-0.10976691, -0.12096859, -0.11087692,
4773	-0.16868259, -0.15868061, -0.16801135,
4774	-0.1588003, -0.16674343, -0.15813323,
4775	-0.16457595, -0.15693826, -0.16281096,
4776	-0.15585154, -0.16283873, -0.15540068,
4777	-0.17450362, -0.19919913, -0.18062882,
4778	-0.19764131, -0.17783111, -0.19407213,
4779	-0.1736915, -0.19053624, -0.17228678,
4780	-0.19105634, -0.17920133, -0.1968828,
4781	-0.14244395, -0.14604641, -0.14473537,
4782	-0.1506107, -0.14510055, -0.14919522,
4783	-0.14175896, -0.14560798, -0.13911658,
4784	-0.14439383, -0.13924047, -0.14829043])
4785
4786	os.remove(domain.get_name() + '.sww')
4787
4788	def test_bedslope_problem_second_order_one_step(self):
4789
4790	from mesh_factory import rectangular
4791
4792	#Create basic mesh
4793	points, vertices, boundary = rectangular(6, 6)
4794
4795	#Create shallow water domain
4796	domain = Domain(points, vertices, boundary)
4797	domain.smooth = False
4798	domain.default_order=2
4799	domain.beta_w = 0.9
4800	domain.beta_w_dry = 0.9
4801	domain.beta_uh = 0.9
4802	domain.beta_uh_dry = 0.9
4803	domain.beta_vh = 0.9
4804	domain.beta_vh_dry = 0.9
4805
4806
4807	# FIXME (Ole): Need tests where this is commented out
4808	domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
4809	domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
4810
4811	#Bed-slope and friction at vertices (and interpolated elsewhere)
4812	def x_slope(x, y):
4813	return -x/3
4814
4815	domain.set_quantity('elevation', x_slope)
4816
4817	# Boundary conditions
4818	Br = Reflective_boundary(domain)
4819	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4820
4821	#Initial condition
4822	domain.set_quantity('stage', expression='elevation+0.05')
4823	domain.check_integrity()
4824
4825	assert num.allclose(domain.quantities['stage'].centroid_values,
4826	[ 0.01296296, 0.03148148, 0.01296296,
4827	0.03148148, 0.01296296, 0.03148148,
4828	0.01296296, 0.03148148, 0.01296296,
4829	0.03148148, 0.01296296, 0.03148148,
4830	-0.04259259, -0.02407407, -0.04259259,
4831	-0.02407407, -0.04259259, -0.02407407,
4832	-0.04259259, -0.02407407, -0.04259259,
4833	-0.02407407, -0.04259259, -0.02407407,
4834	-0.09814815, -0.07962963, -0.09814815,
4835	-0.07962963, -0.09814815, -0.07962963,
4836	-0.09814815, -0.07962963, -0.09814815,
4837	-0.07962963, -0.09814815, -0.07962963,
4838	-0.1537037, -0.13518519, -0.1537037,
4839	-0.13518519, -0.1537037, -0.13518519,
4840	-0.1537037 , -0.13518519, -0.1537037,
4841	-0.13518519, -0.1537037, -0.13518519,
4842	-0.20925926, -0.19074074, -0.20925926,
4843	-0.19074074, -0.20925926, -0.19074074,
4844	-0.20925926, -0.19074074, -0.20925926,
4845	-0.19074074, -0.20925926, -0.19074074,
4846	-0.26481481, -0.2462963, -0.26481481,
4847	-0.2462963, -0.26481481, -0.2462963,
4848	-0.26481481, -0.2462963, -0.26481481,
4849	-0.2462963, -0.26481481, -0.2462963])
4850
4851
4852	#print domain.quantities['stage'].extrapolate_second_order()
4853	#domain.distribute_to_vertices_and_edges()
4854	#print domain.quantities['stage'].vertex_values[:,0]
4855
4856	#Evolution
4857	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
4858	#domain.write_time()
4859	pass
4860
4861
4862	#print domain.quantities['stage'].centroid_values
4863	assert num.allclose(domain.quantities['stage'].centroid_values,
4864	[ 0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
4865	0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
4866	0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
4867	-0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
4868	-0.04186556, -0.0248011, -0.04186556, -0.02255593,
4869	-0.09966629, -0.08035666, -0.09742112, -0.08035666,
4870	-0.09742112, -0.08035666, -0.09742112, -0.08035666,
4871	-0.09742112, -0.08035666, -0.09742112, -0.07811149,
4872	-0.15522185, -0.13591222, -0.15297667, -0.13591222,
4873	-0.15297667, -0.13591222, -0.15297667, -0.13591222,
4874	-0.15297667, -0.13591222, -0.15297667, -0.13366704,
4875	-0.2107774, -0.19146777, -0.20853223, -0.19146777,
4876	-0.20853223, -0.19146777, -0.20853223, -0.19146777,
4877	-0.20853223, -0.19146777, -0.20853223, -0.1892226,
4878	-0.26120669, -0.24776246, -0.25865535, -0.24776246,
4879	-0.25865535, -0.24776246, -0.25865535, -0.24776246,
4880	-0.25865535, -0.24776246, -0.25865535, -0.24521113])
4881
4882	os.remove(domain.get_name() + '.sww')
4883
4884	def test_bedslope_problem_second_order_two_steps(self):
4885
4886	from mesh_factory import rectangular
4887
4888	#Create basic mesh
4889	points, vertices, boundary = rectangular(6, 6)
4890
4891	#Create shallow water domain
4892	domain = Domain(points, vertices, boundary)
4893	domain.smooth = False
4894	domain.default_order=2
4895	domain.beta_w = 0.9
4896	domain.beta_w_dry = 0.9
4897	domain.beta_uh = 0.9
4898	domain.beta_uh_dry = 0.9
4899	domain.beta_vh = 0.9
4900	domain.beta_vh_dry = 0.9
4901
4902	# FIXME (Ole): Need tests where this is commented out
4903	domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
4904	domain.H0 = 0 # Backwards compatibility (6/2/7)
4905	domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
4906
4907
4908	#Bed-slope and friction at vertices (and interpolated elsewhere)
4909	def x_slope(x, y):
4910	return -x/3
4911
4912	domain.set_quantity('elevation', x_slope)
4913
4914	# Boundary conditions
4915	Br = Reflective_boundary(domain)
4916	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
4917
4918	#Initial condition
4919	domain.set_quantity('stage', expression='elevation+0.05')
4920	domain.check_integrity()
4921
4922	assert num.allclose(domain.quantities['stage'].centroid_values,
4923	[ 0.01296296, 0.03148148, 0.01296296,
4924	0.03148148, 0.01296296, 0.03148148,
4925	0.01296296, 0.03148148, 0.01296296,
4926	0.03148148, 0.01296296, 0.03148148,
4927	-0.04259259, -0.02407407, -0.04259259,
4928	-0.02407407, -0.04259259, -0.02407407,
4929	-0.04259259, -0.02407407, -0.04259259,
4930	-0.02407407, -0.04259259, -0.02407407,
4931	-0.09814815, -0.07962963, -0.09814815,
4932	-0.07962963, -0.09814815, -0.07962963,
4933	-0.09814815, -0.07962963, -0.09814815,
4934	-0.07962963, -0.09814815, -0.07962963,
4935	-0.1537037 , -0.13518519, -0.1537037,
4936	-0.13518519, -0.1537037, -0.13518519,
4937	-0.1537037 , -0.13518519, -0.1537037,
4938	-0.13518519, -0.1537037, -0.13518519,
4939	-0.20925926, -0.19074074, -0.20925926,
4940	-0.19074074, -0.20925926, -0.19074074,
4941	-0.20925926, -0.19074074, -0.20925926,
4942	-0.19074074, -0.20925926, -0.19074074,
4943	-0.26481481, -0.2462963, -0.26481481,
4944	-0.2462963, -0.26481481, -0.2462963,
4945	-0.26481481, -0.2462963, -0.26481481,
4946	-0.2462963, -0.26481481, -0.2462963])
4947
4948
4949	#print domain.quantities['stage'].extrapolate_second_order()
4950	#domain.distribute_to_vertices_and_edges()
4951	#print domain.quantities['stage'].vertex_values[:,0]
4952
4953	#Evolution
4954	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
4955	pass
4956
4957
4958	#Data from earlier version of abstract_2d_finite_volumes ft=0.1
4959	assert num.allclose(domain.min_timestep, 0.0376895634803)
4960	assert num.allclose(domain.max_timestep, 0.0415635655309)
4961
4962
4963	assert num.allclose(domain.quantities['stage'].centroid_values,
4964	[ 0.00855788, 0.01575204, 0.00994606, 0.01717072,
4965	0.01005985, 0.01716362, 0.01005985, 0.01716299,
4966	0.01007098, 0.01736248, 0.01216452, 0.02026776,
4967	-0.04462374, -0.02479045, -0.04199789, -0.0229465,
4968	-0.04184033, -0.02295693, -0.04184013, -0.02295675,
4969	-0.04184486, -0.0228168, -0.04028876, -0.02036486,
4970	-0.10029444, -0.08170809, -0.09772846, -0.08021704,
4971	-0.09760006, -0.08022143, -0.09759984, -0.08022124,
4972	-0.09760261, -0.08008893, -0.09603914, -0.07758209,
4973	-0.15584152, -0.13723138, -0.15327266, -0.13572906,
4974	-0.15314427, -0.13573349, -0.15314405, -0.13573331,
4975	-0.15314679, -0.13560104, -0.15158523, -0.13310701,
4976	-0.21208605, -0.19283913, -0.20955631, -0.19134189,
4977	-0.20942821, -0.19134598, -0.20942799, -0.1913458,
4978	-0.20943005, -0.19120952, -0.20781177, -0.18869401,
4979	-0.25384082, -0.2463294, -0.25047649, -0.24464654,
4980	-0.25031159, -0.24464253, -0.25031112, -0.24464253,
4981	-0.25031463, -0.24454764, -0.24885323, -0.24286438])
4982
4983
4984	os.remove(domain.get_name() + '.sww')
4985
4986	def test_bedslope_problem_second_order_two_yieldsteps(self):
4987
4988	from mesh_factory import rectangular
4989
4990	#Create basic mesh
4991	points, vertices, boundary = rectangular(6, 6)
4992
4993	#Create shallow water domain
4994	domain = Domain(points, vertices, boundary)
4995	domain.smooth = False
4996	domain.default_order=2
4997	domain.beta_w = 0.9
4998	domain.beta_w_dry = 0.9
4999	domain.beta_uh = 0.9
5000	domain.beta_uh_dry = 0.9
5001	domain.beta_vh = 0.9
5002	domain.beta_vh_dry = 0.9
5003
5004	# FIXME (Ole): Need tests where this is commented out
5005	domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
5006	domain.H0 = 0 # Backwards compatibility (6/2/7)
5007	domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
5008
5009
5010	#Bed-slope and friction at vertices (and interpolated elsewhere)
5011	def x_slope(x, y):
5012	return -x/3
5013
5014	domain.set_quantity('elevation', x_slope)
5015
5016	# Boundary conditions
5017	Br = Reflective_boundary(domain)
5018	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
5019
5020	#Initial condition
5021	domain.set_quantity('stage', expression='elevation+0.05')
5022	domain.check_integrity()
5023
5024	assert num.allclose(domain.quantities['stage'].centroid_values,
5025	[ 0.01296296, 0.03148148, 0.01296296,
5026	0.03148148, 0.01296296, 0.03148148,
5027	0.01296296, 0.03148148, 0.01296296,
5028	0.03148148, 0.01296296, 0.03148148,
5029	-0.04259259, -0.02407407, -0.04259259,
5030	-0.02407407, -0.04259259, -0.02407407,
5031	-0.04259259, -0.02407407, -0.04259259,
5032	-0.02407407, -0.04259259, -0.02407407,
5033	-0.09814815, -0.07962963, -0.09814815,
5034	-0.07962963, -0.09814815, -0.07962963,
5035	-0.09814815, -0.07962963, -0.09814815,
5036	-0.07962963, -0.09814815, -0.07962963,
5037	-0.1537037 , -0.13518519, -0.1537037,
5038	-0.13518519, -0.1537037, -0.13518519,
5039	-0.1537037 , -0.13518519, -0.1537037,
5040	-0.13518519, -0.1537037, -0.13518519,
5041	-0.20925926, -0.19074074, -0.20925926,
5042	-0.19074074, -0.20925926, -0.19074074,
5043	-0.20925926, -0.19074074, -0.20925926,
5044	-0.19074074, -0.20925926, -0.19074074,
5045	-0.26481481, -0.2462963, -0.26481481,
5046	-0.2462963, -0.26481481, -0.2462963,
5047	-0.26481481, -0.2462963, -0.26481481,
5048	-0.2462963, -0.26481481, -0.2462963])
5049
5050
5051	#print domain.quantities['stage'].extrapolate_second_order()
5052	#domain.distribute_to_vertices_and_edges()
5053	#print domain.quantities['stage'].vertex_values[:,0]
5054
5055	#Evolution
5056	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): #0.05??
5057	#domain.write_time()
5058	pass
5059
5060
5061
5062	assert num.allclose(domain.quantities['stage'].centroid_values,
5063	[ 0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
5064	0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
5065	0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
5066	-0.0229465, -0.04184033, -0.02295693, -0.04184013,
5067	-0.02295675, -0.04184486, -0.0228168, -0.04028876,
5068	-0.02036486, -0.10029444, -0.08170809, -0.09772846,
5069	-0.08021704, -0.09760006, -0.08022143, -0.09759984,
5070	-0.08022124, -0.09760261, -0.08008893, -0.09603914,
5071	-0.07758209, -0.15584152, -0.13723138, -0.15327266,
5072	-0.13572906, -0.15314427, -0.13573349, -0.15314405,
5073	-0.13573331, -0.15314679, -0.13560104, -0.15158523,
5074	-0.13310701, -0.21208605, -0.19283913, -0.20955631,
5075	-0.19134189, -0.20942821, -0.19134598, -0.20942799,
5076	-0.1913458, -0.20943005, -0.19120952, -0.20781177,
5077	-0.18869401, -0.25384082, -0.2463294, -0.25047649,
5078	-0.24464654, -0.25031159, -0.24464253, -0.25031112,
5079	-0.24464253, -0.25031463, -0.24454764, -0.24885323,
5080	-0.24286438])
5081
5082	os.remove(domain.get_name() + '.sww')
5083
5084	def test_bedslope_problem_second_order_more_steps(self):
5085
5086	from mesh_factory import rectangular
5087
5088	#Create basic mesh
5089	points, vertices, boundary = rectangular(6, 6)
5090
5091	#Create shallow water domain
5092	domain = Domain(points, vertices, boundary)
5093	domain.smooth = False
5094	domain.default_order=2
5095	domain.beta_w = 0.9
5096	domain.beta_w_dry = 0.9
5097	domain.beta_uh = 0.9
5098	domain.beta_uh_dry = 0.9
5099	domain.beta_vh = 0.9
5100	domain.beta_vh_dry = 0.9
5101
5102
5103	# FIXME (Ole): Need tests where these two are commented out
5104	domain.H0 = 0 # Backwards compatibility (6/2/7)
5105	domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
5106	domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
5107
5108
5109
5110	#Bed-slope and friction at vertices (and interpolated elsewhere)
5111	def x_slope(x, y):
5112	return -x/3
5113
5114	domain.set_quantity('elevation', x_slope)
5115
5116	# Boundary conditions
5117	Br = Reflective_boundary(domain)
5118	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
5119
5120	#Initial condition
5121	domain.set_quantity('stage', expression = 'elevation + 0.05')
5122	domain.check_integrity()
5123
5124	assert num.allclose(domain.quantities['stage'].centroid_values,
5125	[ 0.01296296, 0.03148148, 0.01296296,
5126	0.03148148, 0.01296296, 0.03148148,
5127	0.01296296, 0.03148148, 0.01296296,
5128	0.03148148, 0.01296296, 0.03148148,
5129	-0.04259259, -0.02407407, -0.04259259,
5130	-0.02407407, -0.04259259, -0.02407407,
5131	-0.04259259, -0.02407407, -0.04259259,
5132	-0.02407407, -0.04259259, -0.02407407,
5133	-0.09814815, -0.07962963, -0.09814815,
5134	-0.07962963, -0.09814815, -0.07962963,
5135	-0.09814815, -0.07962963, -0.09814815,
5136	-0.07962963, -0.09814815, -0.07962963,
5137	-0.1537037 , -0.13518519, -0.1537037,
5138	-0.13518519, -0.1537037, -0.13518519,
5139	-0.1537037 , -0.13518519, -0.1537037,
5140	-0.13518519, -0.1537037, -0.13518519,
5141	-0.20925926, -0.19074074, -0.20925926,
5142	-0.19074074, -0.20925926, -0.19074074,
5143	-0.20925926, -0.19074074, -0.20925926,
5144	-0.19074074, -0.20925926, -0.19074074,
5145	-0.26481481, -0.2462963, -0.26481481,
5146	-0.2462963, -0.26481481, -0.2462963,
5147	-0.26481481, -0.2462963, -0.26481481,
5148	-0.2462963, -0.26481481, -0.2462963])
5149
5150
5151	#print domain.quantities['stage'].extrapolate_second_order()
5152	#domain.distribute_to_vertices_and_edges()
5153	#print domain.quantities['stage'].vertex_values[:,0]
5154
5155	#Evolution
5156	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
5157
5158	# Check that diagnostics works
5159	msg = domain.timestepping_statistics(track_speeds=True)
5160	#FIXME(Ole): One might check the contents of msg here.
5161
5162
5163
5164	assert num.allclose(domain.quantities['stage'].centroid_values,
5165	[-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
5166	-0.02901975, -0.0141361, -0.02898816, -0.01418135, -0.02961409, -0.01403487,
5167	-0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
5168	-0.07571145, -0.06235231, -0.0756817, -0.06245309, -0.07652292, -0.06289946,
5169	-0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934, -0.1107174,
5170	-0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
5171	-0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493, -0.15697919,
5172	-0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423