Context Navigation

test_shallow_water_domain.py @ 6553

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