Context Navigation

test_shallow_water_domain.py @ 5730

Last change on this file since 5730 was 5730, checked in by ole, 16 years ago
Run time - flow through cross section and some refactoring + clean up
File size: 213.3 KB

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