Context Navigation

test_shallow_water_domain.py @ 6541

Last change on this file since 6541 was 6541, checked in by ole, 14 years ago
Implemented fix to ticket:314.
File size: 246.3 KB

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