Context Navigation

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

Last change on this file since 4684 was 4684, checked in by ole, 16 years ago
Flux test
File size: 175.2 KB

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

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

Download in other formats: