Context Navigation

test_shallow_water.py @ 3383

Last change on this file since 3383 was 2812, checked in by ole, 19 years ago
Forgot to pass caching arguments to Domain
File size: 118.6 KB

Line
1	#!/usr/bin/env python
2
3	import unittest, os
4	from math import sqrt, pi
5
6	from config import g, epsilon
7	from Numeric import allclose, array, zeros, ones, Float, take
8	from utilities.numerical_tools import mean
9
10	from shallow_water import *
11
12
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	def test_flux_zero_case(self):
107	ql = zeros( 3, Float )
108	qr = zeros( 3, Float )
109	normal = zeros( 2, Float )
110	zl = zr = 0.
111	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
112
113	assert allclose(flux, [0,0,0])
114	assert max_speed == 0.
115
116	def test_flux_constants(self):
117	w = 2.0
118
119	normal = array([1.,0])
120	ql = array([w, 0, 0])
121	qr = array([w, 0, 0])
122	zl = zr = 0.
123	h = w - (zl+zr)/2
124
125	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
126
127	assert allclose(flux, [0., 0.5gh**2, 0.])
128	assert max_speed == sqrt(g*h)
129
130	#def test_flux_slope(self):
131	# #FIXME: TODO
132	# w = 2.0
133	#
134	# normal = array([1.,0])
135	# ql = array([w, 0, 0])
136	# qr = array([w, 0, 0])
137	# zl = zr = 0.
138	# h = w - (zl+zr)/2
139	#
140	# flux, max_speed = flux_function(normal, ql, qr, zl, zr)
141	#
142	# assert allclose(flux, [0., 0.5gh**2, 0.])
143	# assert max_speed == sqrt(g*h)
144
145
146	def test_flux1(self):
147	#Use data from previous version of pyvolution
148	normal = array([1.,0])
149	ql = array([-0.2, 2, 3])
150	qr = array([-0.2, 2, 3])
151	zl = zr = -0.5
152	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
153
154	assert allclose(flux, [2.,13.77433333, 20.])
155	assert allclose(max_speed, 8.38130948661)
156
157
158	def test_flux2(self):
159	#Use data from previous version of pyvolution
160	normal = array([0., -1.])
161	ql = array([-0.075, 2, 3])
162	qr = array([-0.075, 2, 3])
163	zl = zr = -0.375
164	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
165
166	assert allclose(flux, [-3.,-20.0, -30.441])
167	assert allclose(max_speed, 11.7146428199)
168
169	def test_flux3(self):
170	#Use data from previous version of pyvolution
171	normal = array([-sqrt(2)/2, sqrt(2)/2])
172	ql = array([-0.075, 2, 3])
173	qr = array([-0.075, 2, 3])
174	zl = zr = -0.375
175	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
176
177	assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
178	assert allclose(max_speed, 4.0716654239)
179
180	def test_flux4(self):
181	#Use data from previous version of pyvolution
182	normal = array([-sqrt(2)/2, sqrt(2)/2])
183	ql = array([-0.34319278, 0.10254161, 0.07273855])
184	qr = array([-0.30683287, 0.1071986, 0.05930515])
185	zl = zr = -0.375
186	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
187
188	assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
189	assert allclose(max_speed, 1.31414103233)
190
191	def test_sw_domain_simple(self):
192	a = [0.0, 0.0]
193	b = [0.0, 2.0]
194	c = [2.0,0.0]
195	d = [0.0, 4.0]
196	e = [2.0, 2.0]
197	f = [4.0,0.0]
198
199	points = [a, b, c, d, e, f]
200	#bac, bce, ecf, dbe, daf, dae
201	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
202
203	domain = Domain(points, vertices)
204	domain.check_integrity()
205
206	for name in ['stage', 'xmomentum', 'ymomentum',
207	'elevation', 'friction']:
208	assert domain.quantities.has_key(name)
209
210
211	assert domain.get_conserved_quantities(0, edge=1) == 0.
212
213
214	def test_boundary_conditions(self):
215
216	a = [0.0, 0.0]
217	b = [0.0, 2.0]
218	c = [2.0,0.0]
219	d = [0.0, 4.0]
220	e = [2.0, 2.0]
221	f = [4.0,0.0]
222
223	points = [a, b, c, d, e, f]
224	#bac, bce, ecf, dbe
225	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
226	boundary = { (0, 0): 'Third',
227	(0, 2): 'First',
228	(2, 0): 'Second',
229	(2, 1): 'Second',
230	(3, 1): 'Second',
231	(3, 2): 'Third'}
232
233
234	domain = Domain(points, vertices, boundary)
235	domain.check_integrity()
236
237
238	domain.set_quantity('stage', [[1,2,3], [5,5,5],
239	[0,0,9], [-6, 3, 3]])
240
241	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
242	[3,3,3], [4, 4, 4]])
243
244	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
245	[30,30,30], [40, 40, 40]])
246
247
248	D = Dirichlet_boundary([5,2,1])
249	T = Transmissive_boundary(domain)
250	R = Reflective_boundary(domain)
251	domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
252
253	domain.update_boundary()
254
255	#Stage
256	assert domain.quantities['stage'].boundary_values[0] == 2.5
257	assert domain.quantities['stage'].boundary_values[0] ==\
258	domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
259	assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
260	assert domain.quantities['stage'].boundary_values[2] ==\
261	domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
262	assert domain.quantities['stage'].boundary_values[3] ==\
263	domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
264	assert domain.quantities['stage'].boundary_values[4] ==\
265	domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
266	assert domain.quantities['stage'].boundary_values[5] ==\
267	domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
268
269	#Xmomentum
270	assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
271	assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
272	assert domain.quantities['xmomentum'].boundary_values[2] ==\
273	domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
274	assert domain.quantities['xmomentum'].boundary_values[3] ==\
275	domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
276	assert domain.quantities['xmomentum'].boundary_values[4] ==\
277	domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
278	assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective
279
280	#Ymomentum
281	assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
282	assert domain.quantities['ymomentum'].boundary_values[1] == 1. #Dirichlet
283	assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
284	assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
285	assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
286	assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
287
288
289	def test_boundary_conditionsII(self):
290
291	a = [0.0, 0.0]
292	b = [0.0, 2.0]
293	c = [2.0,0.0]
294	d = [0.0, 4.0]
295	e = [2.0, 2.0]
296	f = [4.0,0.0]
297
298	points = [a, b, c, d, e, f]
299	#bac, bce, ecf, dbe
300	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
301	boundary = { (0, 0): 'Third',
302	(0, 2): 'First',
303	(2, 0): 'Second',
304	(2, 1): 'Second',
305	(3, 1): 'Second',
306	(3, 2): 'Third',
307	(0, 1): 'Internal'}
308
309
310	domain = Domain(points, vertices, boundary)
311	domain.check_integrity()
312
313
314	domain.set_quantity('stage', [[1,2,3], [5,5,5],
315	[0,0,9], [-6, 3, 3]])
316
317	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
318	[3,3,3], [4, 4, 4]])
319
320	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
321	[30,30,30], [40, 40, 40]])
322
323
324	D = Dirichlet_boundary([5,2,1])
325	T = Transmissive_boundary(domain)
326	R = Reflective_boundary(domain)
327	domain.set_boundary( {'First': D, 'Second': T,
328	'Third': R, 'Internal': None})
329
330	domain.update_boundary()
331	domain.check_integrity()
332
333
334	def test_compute_fluxes0(self):
335	#Do a full triangle and check that fluxes cancel out for
336	#the constant stage case
337
338	a = [0.0, 0.0]
339	b = [0.0, 2.0]
340	c = [2.0,0.0]
341	d = [0.0, 4.0]
342	e = [2.0, 2.0]
343	f = [4.0,0.0]
344
345	points = [a, b, c, d, e, f]
346	#bac, bce, ecf, dbe
347	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
348
349	domain = Domain(points, vertices)
350	domain.set_quantity('stage', [[2,2,2], [2,2,2],
351	[2,2,2], [2,2,2]])
352	domain.check_integrity()
353
354	assert allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
355	assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
356
357	zl=zr=0. #Assume flat bed
358
359	#Flux across right edge of volume 1
360	normal = domain.get_normal(1,0)
361	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
362	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
363	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
364
365	#Check that flux seen from other triangles is inverse
366	tmp = qr; qr=ql; ql=tmp
367	normal = domain.get_normal(2,2)
368	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
369	assert allclose(flux + flux0, 0.)
370
371	#Flux across upper edge of volume 1
372	normal = domain.get_normal(1,1)
373	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
374	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
375	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
376
377	#Check that flux seen from other triangles is inverse
378	tmp = qr; qr=ql; ql=tmp
379	normal = domain.get_normal(3,0)
380	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
381	assert allclose(flux + flux1, 0.)
382
383	#Flux across lower left hypotenuse of volume 1
384	normal = domain.get_normal(1,2)
385	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
386	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
387	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
388
389	#Check that flux seen from other triangles is inverse
390	tmp = qr; qr=ql; ql=tmp
391	normal = domain.get_normal(0,1)
392	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
393	assert allclose(flux + flux2, 0.)
394
395
396	#Scale by edgelengths, add up anc check that total flux is zero
397	e0 = domain.edgelengths[1, 0]
398	e1 = domain.edgelengths[1, 1]
399	e2 = domain.edgelengths[1, 2]
400
401	assert allclose(e0flux0+e1flux1+e2*flux2, 0.)
402
403	#Now check that compute_flux yields zeros as well
404	domain.compute_fluxes()
405
406	for name in ['stage', 'xmomentum', 'ymomentum']:
407	#print name, domain.quantities[name].explicit_update
408	assert allclose(domain.quantities[name].explicit_update[1], 0)
409
410
411
412	def test_compute_fluxes1(self):
413	#Use values from previous version
414
415	a = [0.0, 0.0]
416	b = [0.0, 2.0]
417	c = [2.0,0.0]
418	d = [0.0, 4.0]
419	e = [2.0, 2.0]
420	f = [4.0,0.0]
421
422	points = [a, b, c, d, e, f]
423	#bac, bce, ecf, dbe
424	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
425
426	domain = Domain(points, vertices)
427	val0 = 2.+2.0/3
428	val1 = 4.+4.0/3
429	val2 = 8.+2.0/3
430	val3 = 2.+8.0/3
431
432	domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
433	[val2, val2, val2], [val3, val3, val3]])
434	domain.check_integrity()
435
436	zl=zr=0. #Assume flat bed
437
438	#Flux across right edge of volume 1
439	normal = domain.get_normal(1,0) #Get normal 0 of triangle 1
440	assert allclose(normal, [1, 0])
441
442	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
443	assert allclose(ql, [val1, 0, 0])
444
445	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
446	assert allclose(qr, [val2, 0, 0])
447
448	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
449
450	#Flux across edge in the east direction (as per normal vector)
451	assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
452	assert allclose(max_speed, 9.21592824046)
453
454
455	#Flux across edge in the west direction (opposite sign for xmomentum)
456	normal_opposite = domain.get_normal(2,2) #Get normal 2 of triangle 2
457	assert allclose(normal_opposite, [-1, 0])
458	flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr)
459	assert allclose(flux_opposite, [-15.3598804, -253.71111111, 0.])
460
461
462	#Flux across upper edge of volume 1
463	normal = domain.get_normal(1,1)
464	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
465	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
466	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
467	assert allclose(flux1, [2.4098563, 0., 123.04444444])
468	assert allclose(max_speed, 7.22956891292)
469
470	#Flux across lower left hypotenuse of volume 1
471	normal = domain.get_normal(1,2)
472	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
473	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
474	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
475
476	assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
477	assert allclose(max_speed, 7.22956891292)
478
479	#Scale, add up and check that compute_fluxes is correct for vol 1
480	e0 = domain.edgelengths[1, 0]
481	e1 = domain.edgelengths[1, 1]
482	e2 = domain.edgelengths[1, 2]
483
484	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
485	assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
486
487
488	domain.compute_fluxes()
489
490	#assert allclose(total_flux, domain.explicit_update[1,:])
491	for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
492	assert allclose(total_flux[i],
493	domain.quantities[name].explicit_update[1])
494
495	#assert allclose(domain.explicit_update, [
496	# [0., -69.68888889, -69.68888889],
497	# [-0.68218178, -166.6, -35.93333333],
498	# [-111.77316251, 69.68888889, 0.],
499	# [-35.68522449, 0., 69.68888889]])
500
501	assert allclose(domain.quantities['stage'].explicit_update,
502	[0., -0.68218178, -111.77316251, -35.68522449])
503	assert allclose(domain.quantities['xmomentum'].explicit_update,
504	[-69.68888889, -166.6, 69.68888889, 0])
505	assert allclose(domain.quantities['ymomentum'].explicit_update,
506	[-69.68888889, -35.93333333, 0., 69.68888889])
507
508
509	#assert allclose(domain.quantities[name].explicit_update
510
511
512
513
514
515	def test_compute_fluxes2(self):
516	#Random values, incl momentum
517
518	a = [0.0, 0.0]
519	b = [0.0, 2.0]
520	c = [2.0,0.0]
521	d = [0.0, 4.0]
522	e = [2.0, 2.0]
523	f = [4.0,0.0]
524
525	points = [a, b, c, d, e, f]
526	#bac, bce, ecf, dbe
527	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
528
529	domain = Domain(points, vertices)
530	val0 = 2.+2.0/3
531	val1 = 4.+4.0/3
532	val2 = 8.+2.0/3
533	val3 = 2.+8.0/3
534
535	zl=zr=0 #Assume flat zero bed
536
537	domain.set_quantity('elevation', zl*ones( (4,3) ))
538
539
540	domain.set_quantity('stage', [[val0, val0-1, val0-2],
541	[val1, val1+1, val1],
542	[val2, val2-2, val2],
543	[val3-0.5, val3, val3]])
544
545	domain.set_quantity('xmomentum',
546	[[1, 2, 3], [3, 4, 5],
547	[1, -1, 0], [0, -2, 2]])
548
549	domain.set_quantity('ymomentum',
550	[[1, -1, 0], [0, -3, 2],
551	[0, 1, 0], [-1, 2, 2]])
552
553
554	domain.check_integrity()
555
556
557
558	#Flux across right edge of volume 1
559	normal = domain.get_normal(1,0)
560	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
561	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
562	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
563
564	#Flux across upper edge of volume 1
565	normal = domain.get_normal(1,1)
566	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
567	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
568	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
569
570	#Flux across lower left hypotenuse of volume 1
571	normal = domain.get_normal(1,2)
572	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
573	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
574	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
575
576	#Scale, add up and check that compute_fluxes is correct for vol 1
577	e0 = domain.edgelengths[1, 0]
578	e1 = domain.edgelengths[1, 1]
579	e2 = domain.edgelengths[1, 2]
580
581	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
582
583
584	domain.compute_fluxes()
585	for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
586	assert allclose(total_flux[i],
587	domain.quantities[name].explicit_update[1])
588	#assert allclose(total_flux, domain.explicit_update[1,:])
589
590
591	def test_compute_fluxes3(self):
592	#Random values, incl momentum
593
594	a = [0.0, 0.0]
595	b = [0.0, 2.0]
596	c = [2.0,0.0]
597	d = [0.0, 4.0]
598	e = [2.0, 2.0]
599	f = [4.0,0.0]
600
601	points = [a, b, c, d, e, f]
602	#bac, bce, ecf, dbe
603	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
604
605	domain = Domain(points, vertices)
606	val0 = 2.+2.0/3
607	val1 = 4.+4.0/3
608	val2 = 8.+2.0/3
609	val3 = 2.+8.0/3
610
611	zl=zr=-3.75 #Assume constant bed (must be less than stage)
612	domain.set_quantity('elevation', zl*ones( (4,3) ))
613
614
615	domain.set_quantity('stage', [[val0, val0-1, val0-2],
616	[val1, val1+1, val1],
617	[val2, val2-2, val2],
618	[val3-0.5, val3, val3]])
619
620	domain.set_quantity('xmomentum',
621	[[1, 2, 3], [3, 4, 5],
622	[1, -1, 0], [0, -2, 2]])
623
624	domain.set_quantity('ymomentum',
625	[[1, -1, 0], [0, -3, 2],
626	[0, 1, 0], [-1, 2, 2]])
627
628
629	domain.check_integrity()
630
631
632
633	#Flux across right edge of volume 1
634	normal = domain.get_normal(1,0)
635	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
636	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
637	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
638
639	#Flux across upper edge of volume 1
640	normal = domain.get_normal(1,1)
641	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
642	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
643	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
644
645	#Flux across lower left hypotenuse of volume 1
646	normal = domain.get_normal(1,2)
647	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
648	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
649	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
650
651	#Scale, add up and check that compute_fluxes is correct for vol 1
652	e0 = domain.edgelengths[1, 0]
653	e1 = domain.edgelengths[1, 1]
654	e2 = domain.edgelengths[1, 2]
655
656	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
657
658	domain.compute_fluxes()
659	for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
660	assert allclose(total_flux[i],
661	domain.quantities[name].explicit_update[1])
662
663
664
665	def test_catching_negative_heights(self):
666
667	a = [0.0, 0.0]
668	b = [0.0, 2.0]
669	c = [2.0,0.0]
670	d = [0.0, 4.0]
671	e = [2.0, 2.0]
672	f = [4.0,0.0]
673
674	points = [a, b, c, d, e, f]
675	#bac, bce, ecf, dbe
676	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
677
678	domain = Domain(points, vertices)
679	val0 = 2.+2.0/3
680	val1 = 4.+4.0/3
681	val2 = 8.+2.0/3
682	val3 = 2.+8.0/3
683
684	zl=zr=4 #Too large
685	domain.set_quantity('elevation', zl*ones( (4,3) ))
686	domain.set_quantity('stage', [[val0, val0-1, val0-2],
687	[val1, val1+1, val1],
688	[val2, val2-2, val2],
689	[val3-0.5, val3, val3]])
690
691	#Should fail
692	try:
693	domain.check_integrity()
694	except:
695	pass
696
697
698
699
700	#####################################################
701	def test_initial_condition(self):
702	"""Test that initial condition is output at time == 0
703	"""
704
705	from config import g
706	import copy
707
708	a = [0.0, 0.0]
709	b = [0.0, 2.0]
710	c = [2.0, 0.0]
711	d = [0.0, 4.0]
712	e = [2.0, 2.0]
713	f = [4.0, 0.0]
714
715	points = [a, b, c, d, e, f]
716	#bac, bce, ecf, dbe
717	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
718
719	domain = Domain(points, vertices)
720
721	#Set up for a gradient of (3,0) at mid triangle
722	def slope(x, y):
723	return 3*x
724
725	h = 0.1
726	def stage(x,y):
727	return slope(x,y)+h
728
729	domain.set_quantity('elevation', slope)
730	domain.set_quantity('stage', stage)
731
732	initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
733
734	domain.set_boundary({'exterior': Reflective_boundary(domain)})
735
736	#Evolution
737	for t in domain.evolve(yieldstep = 1.0, finaltime = 2.0):
738	stage = domain.quantities['stage'].vertex_values
739
740	if t == 0.0:
741	assert allclose(stage, initial_stage)
742	else:
743	assert not allclose(stage, initial_stage)
744
745	os.remove(domain.filename + '.sww')
746
747
748
749	#####################################################
750	def test_gravity(self):
751	#Assuming no friction
752
753	from config import g
754
755	a = [0.0, 0.0]
756	b = [0.0, 2.0]
757	c = [2.0, 0.0]
758	d = [0.0, 4.0]
759	e = [2.0, 2.0]
760	f = [4.0, 0.0]
761
762	points = [a, b, c, d, e, f]
763	#bac, bce, ecf, dbe
764	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
765
766	domain = Domain(points, vertices)
767
768	#Set up for a gradient of (3,0) at mid triangle
769	def slope(x, y):
770	return 3*x
771
772	h = 0.1
773	def stage(x,y):
774	return slope(x,y)+h
775
776	domain.set_quantity('elevation', slope)
777	domain.set_quantity('stage', stage)
778
779	for name in domain.conserved_quantities:
780	assert allclose(domain.quantities[name].explicit_update, 0)
781	assert allclose(domain.quantities[name].semi_implicit_update, 0)
782
783	domain.compute_forcing_terms()
784
785	assert allclose(domain.quantities['stage'].explicit_update, 0)
786	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
787	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
788
789
790	def test_manning_friction(self):
791	from config import g
792
793	a = [0.0, 0.0]
794	b = [0.0, 2.0]
795	c = [2.0, 0.0]
796	d = [0.0, 4.0]
797	e = [2.0, 2.0]
798	f = [4.0, 0.0]
799
800	points = [a, b, c, d, e, f]
801	#bac, bce, ecf, dbe
802	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
803
804	domain = Domain(points, vertices)
805
806	#Set up for a gradient of (3,0) at mid triangle
807	def slope(x, y):
808	return 3*x
809
810	h = 0.1
811	def stage(x,y):
812	return slope(x,y)+h
813
814	eta = 0.07
815	domain.set_quantity('elevation', slope)
816	domain.set_quantity('stage', stage)
817	domain.set_quantity('friction', eta)
818
819	for name in domain.conserved_quantities:
820	assert allclose(domain.quantities[name].explicit_update, 0)
821	assert allclose(domain.quantities[name].semi_implicit_update, 0)
822
823	domain.compute_forcing_terms()
824
825	assert allclose(domain.quantities['stage'].explicit_update, 0)
826	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
827	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
828
829	assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
830	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
831	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
832
833	#Create some momentum for friction to work with
834	domain.set_quantity('xmomentum', 1)
835	S = -g * eta2 / h(7.0/3)
836
837	domain.compute_forcing_terms()
838	assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
839	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
840	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
841
842	#A more complex example
843	domain.quantities['stage'].semi_implicit_update[:] = 0.0
844	domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
845	domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
846
847	domain.set_quantity('xmomentum', 3)
848	domain.set_quantity('ymomentum', 4)
849
850	S = -g * eta*2 5 / h**(7.0/3)
851
852
853	domain.compute_forcing_terms()
854
855	assert allclose(domain.quantities['stage'].semi_implicit_update, 0)
856	assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
857	assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
858
859	def test_constant_wind_stress(self):
860	from config import rho_a, rho_w, eta_w
861	from math import pi, cos, sin, sqrt
862
863	a = [0.0, 0.0]
864	b = [0.0, 2.0]
865	c = [2.0, 0.0]
866	d = [0.0, 4.0]
867	e = [2.0, 2.0]
868	f = [4.0, 0.0]
869
870	points = [a, b, c, d, e, f]
871	#bac, bce, ecf, dbe
872	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
873
874
875	domain = Domain(points, vertices)
876
877	#Flat surface with 1m of water
878	domain.set_quantity('elevation', 0)
879	domain.set_quantity('stage', 1.0)
880	domain.set_quantity('friction', 0)
881
882	Br = Reflective_boundary(domain)
883	domain.set_boundary({'exterior': Br})
884
885	#Setup only one forcing term, constant wind stress
886	s = 100
887	phi = 135
888	domain.forcing_terms = []
889	domain.forcing_terms.append( Wind_stress(s, phi) )
890
891	domain.compute_forcing_terms()
892
893
894	const = eta_w*rho_a/rho_w
895
896	#Convert to radians
897	phi = phi*pi/180
898
899	#Compute velocity vector (u, v)
900	u = s*cos(phi)
901	v = s*sin(phi)
902
903	#Compute wind stress
904	S = const * sqrt(u2 + v2)
905
906	assert allclose(domain.quantities['stage'].explicit_update, 0)
907	assert allclose(domain.quantities['xmomentum'].explicit_update, S*u)
908	assert allclose(domain.quantities['ymomentum'].explicit_update, S*v)
909
910
911	def test_variable_wind_stress(self):
912	from config import rho_a, rho_w, eta_w
913	from math import pi, cos, sin, sqrt
914
915	a = [0.0, 0.0]
916	b = [0.0, 2.0]
917	c = [2.0, 0.0]
918	d = [0.0, 4.0]
919	e = [2.0, 2.0]
920	f = [4.0, 0.0]
921
922	points = [a, b, c, d, e, f]
923	#bac, bce, ecf, dbe
924	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
925
926	domain = Domain(points, vertices)
927
928	#Flat surface with 1m of water
929	domain.set_quantity('elevation', 0)
930	domain.set_quantity('stage', 1.0)
931	domain.set_quantity('friction', 0)
932
933	Br = Reflective_boundary(domain)
934	domain.set_boundary({'exterior': Br})
935
936
937	domain.time = 5.54 #Take a random time (not zero)
938
939	#Setup only one forcing term, constant wind stress
940	s = 100
941	phi = 135
942	domain.forcing_terms = []
943	domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
944
945	domain.compute_forcing_terms()
946
947	#Compute reference solution
948	const = eta_w*rho_a/rho_w
949
950	N = domain.number_of_elements
951
952	xc = domain.get_centroid_coordinates()
953	t = domain.time
954
955	x = xc[:,0]
956	y = xc[:,1]
957	s_vec = speed(t,x,y)
958	phi_vec = angle(t,x,y)
959
960
961	for k in range(N):
962	#Convert to radians
963	phi = phi_vec[k]*pi/180
964	s = s_vec[k]
965
966	#Compute velocity vector (u, v)
967	u = s*cos(phi)
968	v = s*sin(phi)
969
970	#Compute wind stress
971	S = const * sqrt(u2 + v2)
972
973	assert allclose(domain.quantities['stage'].explicit_update[k], 0)
974	assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
975	assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
976
977
978
979
980	def test_windfield_from_file(self):
981	from config import rho_a, rho_w, eta_w
982	from math import pi, cos, sin, sqrt
983	from config import time_format
984	from util import file_function
985	import time
986
987
988	a = [0.0, 0.0]
989	b = [0.0, 2.0]
990	c = [2.0, 0.0]
991	d = [0.0, 4.0]
992	e = [2.0, 2.0]
993	f = [4.0, 0.0]
994
995	points = [a, b, c, d, e, f]
996	#bac, bce, ecf, dbe
997	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
998
999	domain = Domain(points, vertices)
1000
1001	#Flat surface with 1m of water
1002	domain.set_quantity('elevation', 0)
1003	domain.set_quantity('stage', 1.0)
1004	domain.set_quantity('friction', 0)
1005
1006	Br = Reflective_boundary(domain)
1007	domain.set_boundary({'exterior': Br})
1008
1009
1010	domain.time = 7 #Take a time that is represented in file (not zero)
1011
1012	#Write wind stress file (ensure that domain.time is covered)
1013	#Take x=1 and y=0
1014	filename = 'test_windstress_from_file'
1015	start = time.mktime(time.strptime('2000', '%Y'))
1016	fid = open(filename + '.txt', 'w')
1017	dt = 1 #One second interval
1018	t = 0.0
1019	while t <= 10.0:
1020	t_string = time.strftime(time_format, time.gmtime(t+start))
1021
1022	fid.write('%s, %f %f\n' %(t_string,
1023	speed(t,[1],[0])[0],
1024	angle(t,[1],[0])[0]))
1025	t += dt
1026
1027	fid.close()
1028
1029
1030	#Convert ASCII file to NetCDF (Which is what we really like!)
1031	from data_manager import timefile2netcdf
1032	timefile2netcdf(filename)
1033	os.remove(filename + '.txt')
1034
1035
1036	#Setup wind stress
1037	F = file_function(filename + '.tms', quantities = ['Attribute0',
1038	'Attribute1'])
1039	os.remove(filename + '.tms')
1040
1041
1042	#print 'F(5)', F(5)
1043
1044	#print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
1045
1046	#print dir(F)
1047	#print F.T
1048	#print F.precomputed_values
1049	#
1050	#F = file_function(filename + '.txt')
1051	#
1052	#print dir(F)
1053	#print F.T
1054	#print F.Q
1055
1056	W = Wind_stress(F)
1057
1058	domain.forcing_terms = []
1059	domain.forcing_terms.append(W)
1060
1061	domain.compute_forcing_terms()
1062
1063	#Compute reference solution
1064	const = eta_w*rho_a/rho_w
1065
1066	N = domain.number_of_elements
1067
1068	t = domain.time
1069
1070	s = speed(t,[1],[0])[0]
1071	phi = angle(t,[1],[0])[0]
1072
1073	#Convert to radians
1074	phi = phi*pi/180
1075
1076
1077	#Compute velocity vector (u, v)
1078	u = s*cos(phi)
1079	v = s*sin(phi)
1080
1081	#Compute wind stress
1082	S = const * sqrt(u2 + v2)
1083
1084	for k in range(N):
1085	assert allclose(domain.quantities['stage'].explicit_update[k], 0)
1086	assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
1087	assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
1088
1089
1090
1091
1092	def test_wind_stress_error_condition(self):
1093	"""Test that windstress reacts properly when forcing functions
1094	are wrong - e.g. returns a scalar
1095	"""
1096
1097	from config import rho_a, rho_w, eta_w
1098	from math import pi, cos, sin, sqrt
1099
1100	a = [0.0, 0.0]
1101	b = [0.0, 2.0]
1102	c = [2.0, 0.0]
1103	d = [0.0, 4.0]
1104	e = [2.0, 2.0]
1105	f = [4.0, 0.0]
1106
1107	points = [a, b, c, d, e, f]
1108	#bac, bce, ecf, dbe
1109	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1110
1111	domain = Domain(points, vertices)
1112
1113	#Flat surface with 1m of water
1114	domain.set_quantity('elevation', 0)
1115	domain.set_quantity('stage', 1.0)
1116	domain.set_quantity('friction', 0)
1117
1118	Br = Reflective_boundary(domain)
1119	domain.set_boundary({'exterior': Br})
1120
1121
1122	domain.time = 5.54 #Take a random time (not zero)
1123
1124	#Setup only one forcing term, bad func
1125	domain.forcing_terms = []
1126
1127	try:
1128	domain.forcing_terms.append(Wind_stress(s = scalar_func,
1129	phi = angle))
1130	except AssertionError:
1131	pass
1132	else:
1133	msg = 'Should have raised exception'
1134	raise msg
1135
1136
1137	try:
1138	domain.forcing_terms.append(Wind_stress(s = speed,
1139	phi = scalar_func))
1140	except AssertionError:
1141	pass
1142	else:
1143	msg = 'Should have raised exception'
1144	raise msg
1145
1146	try:
1147	domain.forcing_terms.append(Wind_stress(s = speed,
1148	phi = 'xx'))
1149	except:
1150	pass
1151	else:
1152	msg = 'Should have raised exception'
1153	raise msg
1154
1155
1156	#####################################################
1157	def test_first_order_extrapolator_const_z(self):
1158
1159	a = [0.0, 0.0]
1160	b = [0.0, 2.0]
1161	c = [2.0, 0.0]
1162	d = [0.0, 4.0]
1163	e = [2.0, 2.0]
1164	f = [4.0, 0.0]
1165
1166	points = [a, b, c, d, e, f]
1167	#bac, bce, ecf, dbe
1168	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1169
1170	domain = Domain(points, vertices)
1171	val0 = 2.+2.0/3
1172	val1 = 4.+4.0/3
1173	val2 = 8.+2.0/3
1174	val3 = 2.+8.0/3
1175
1176	zl=zr=-3.75 #Assume constant bed (must be less than stage)
1177	domain.set_quantity('elevation', zl*ones( (4,3) ))
1178	domain.set_quantity('stage', [[val0, val0-1, val0-2],
1179	[val1, val1+1, val1],
1180	[val2, val2-2, val2],
1181	[val3-0.5, val3, val3]])
1182
1183
1184
1185	domain.order = 1
1186	domain.distribute_to_vertices_and_edges()
1187
1188	#Check that centroid values were distributed to vertices
1189	C = domain.quantities['stage'].centroid_values
1190	for i in range(3):
1191	assert allclose( domain.quantities['stage'].vertex_values[:,i], C)
1192
1193
1194	def test_first_order_limiter_variable_z(self):
1195	#Check that first order limiter follows bed_slope
1196	from Numeric import alltrue, greater_equal
1197	from config import epsilon
1198
1199	a = [0.0, 0.0]
1200	b = [0.0, 2.0]
1201	c = [2.0,0.0]
1202	d = [0.0, 4.0]
1203	e = [2.0, 2.0]
1204	f = [4.0,0.0]
1205
1206	points = [a, b, c, d, e, f]
1207	#bac, bce, ecf, dbe
1208	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1209
1210	domain = Domain(points, vertices)
1211	val0 = 2.+2.0/3
1212	val1 = 4.+4.0/3
1213	val2 = 8.+2.0/3
1214	val3 = 2.+8.0/3
1215
1216	domain.set_quantity('elevation', [[0,0,0], [6,0,0],
1217	[6,6,6], [6,6,6]])
1218	domain.set_quantity('stage', [[val0, val0, val0],
1219	[val1, val1, val1],
1220	[val2, val2, val2],
1221	[val3, val3, val3]])
1222
1223	E = domain.quantities['elevation'].vertex_values
1224	L = domain.quantities['stage'].vertex_values
1225
1226
1227	#Check that some stages are not above elevation (within eps)
1228	#- so that the limiter has something to work with
1229	assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
1230
1231	domain.order = 1
1232	domain.distribute_to_vertices_and_edges()
1233
1234	#Check that all stages are above elevation (within eps)
1235	assert alltrue(alltrue(greater_equal(L,E-epsilon)))
1236
1237
1238	#####################################################
1239	def test_distribute_basic(self):
1240	#Using test data generated by pyvolution-2
1241	#Assuming no friction and flat bed (0.0)
1242
1243	a = [0.0, 0.0]
1244	b = [0.0, 2.0]
1245	c = [2.0, 0.0]
1246	d = [0.0, 4.0]
1247	e = [2.0, 2.0]
1248	f = [4.0, 0.0]
1249
1250	points = [a, b, c, d, e, f]
1251	#bac, bce, ecf, dbe
1252	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1253
1254	domain = Domain(points, vertices)
1255
1256	val0 = 2.
1257	val1 = 4.
1258	val2 = 8.
1259	val3 = 2.
1260
1261	domain.set_quantity('stage', [val0, val1, val2, val3],
1262	location='centroids')
1263	L = domain.quantities['stage'].vertex_values
1264
1265	#First order
1266	domain.order = 1
1267	domain.distribute_to_vertices_and_edges()
1268	assert allclose(L[1], val1)
1269
1270	#Second order
1271	domain.order = 2
1272	domain.distribute_to_vertices_and_edges()
1273	assert allclose(L[1], [2.2, 4.9, 4.9])
1274
1275
1276
1277	def test_distribute_away_from_bed(self):
1278	#Using test data generated by pyvolution-2
1279	#Assuming no friction and flat bed (0.0)
1280
1281	a = [0.0, 0.0]
1282	b = [0.0, 2.0]
1283	c = [2.0, 0.0]
1284	d = [0.0, 4.0]
1285	e = [2.0, 2.0]
1286	f = [4.0, 0.0]
1287
1288	points = [a, b, c, d, e, f]
1289	#bac, bce, ecf, dbe
1290	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1291
1292	domain = Domain(points, vertices)
1293	L = domain.quantities['stage'].vertex_values
1294
1295	def stage(x,y):
1296	return x**2
1297
1298	domain.set_quantity('stage', stage, location='centroids')
1299
1300	a, b = domain.quantities['stage'].compute_gradients()
1301	assert allclose(a[1], 3.33333334)
1302	assert allclose(b[1], 0.0)
1303
1304	domain.order = 1
1305	domain.distribute_to_vertices_and_edges()
1306	assert allclose(L[1], 1.77777778)
1307
1308	domain.order = 2
1309	domain.distribute_to_vertices_and_edges()
1310	assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
1311
1312
1313
1314	def test_distribute_away_from_bed1(self):
1315	#Using test data generated by pyvolution-2
1316	#Assuming no friction and flat bed (0.0)
1317
1318	a = [0.0, 0.0]
1319	b = [0.0, 2.0]
1320	c = [2.0, 0.0]
1321	d = [0.0, 4.0]
1322	e = [2.0, 2.0]
1323	f = [4.0, 0.0]
1324
1325	points = [a, b, c, d, e, f]
1326	#bac, bce, ecf, dbe
1327	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1328
1329	domain = Domain(points, vertices)
1330	L = domain.quantities['stage'].vertex_values
1331
1332	def stage(x,y):
1333	return x4+y2
1334
1335	domain.set_quantity('stage', stage, location='centroids')
1336	#print domain.quantities['stage'].centroid_values
1337
1338	a, b = domain.quantities['stage'].compute_gradients()
1339	assert allclose(a[1], 25.18518519)
1340	assert allclose(b[1], 3.33333333)
1341
1342	domain.order = 1
1343	domain.distribute_to_vertices_and_edges()
1344	assert allclose(L[1], 4.9382716)
1345
1346	domain.order = 2
1347	domain.distribute_to_vertices_and_edges()
1348	assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
1349
1350
1351
1352	def test_distribute_near_bed(self):
1353	#Using test data generated by pyvolution-2
1354	#Assuming no friction and flat bed (0.0)
1355
1356	a = [0.0, 0.0]
1357	b = [0.0, 2.0]
1358	c = [2.0, 0.0]
1359	d = [0.0, 4.0]
1360	e = [2.0, 2.0]
1361	f = [4.0, 0.0]
1362
1363	points = [a, b, c, d, e, f]
1364	#bac, bce, ecf, dbe
1365	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1366
1367	domain = Domain(points, vertices)
1368
1369
1370	#Set up for a gradient of (3,0) at mid triangle
1371	def slope(x, y):
1372	return 10*x
1373
1374	h = 0.1
1375	def stage(x,y):
1376	return slope(x,y)+h
1377
1378	domain.set_quantity('elevation', slope)
1379	domain.set_quantity('stage', stage, location='centroids')
1380
1381	#print domain.quantities['elevation'].centroid_values
1382	#print domain.quantities['stage'].centroid_values
1383
1384	E = domain.quantities['elevation'].vertex_values
1385	L = domain.quantities['stage'].vertex_values
1386
1387	#print E
1388	domain.order = 1
1389	domain.distribute_to_vertices_and_edges()
1390	##assert allclose(L[1], [0.19999999, 20.05, 20.05])
1391	assert allclose(L[1], [0.1, 20.1, 20.1])
1392
1393	domain.order = 2
1394	domain.distribute_to_vertices_and_edges()
1395	assert allclose(L[1], [0.1, 20.1, 20.1])
1396
1397	def test_distribute_near_bed1(self):
1398	#Using test data generated by pyvolution-2
1399	#Assuming no friction and flat bed (0.0)
1400
1401	a = [0.0, 0.0]
1402	b = [0.0, 2.0]
1403	c = [2.0, 0.0]
1404	d = [0.0, 4.0]
1405	e = [2.0, 2.0]
1406	f = [4.0, 0.0]
1407
1408	points = [a, b, c, d, e, f]
1409	#bac, bce, ecf, dbe
1410	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1411
1412	domain = Domain(points, vertices)
1413
1414
1415	#Set up for a gradient of (3,0) at mid triangle
1416	def slope(x, y):
1417	return x4+y2
1418
1419	h = 0.1
1420	def stage(x,y):
1421	return slope(x,y)+h
1422
1423	domain.set_quantity('elevation', slope)
1424	domain.set_quantity('stage', stage)
1425
1426	#print domain.quantities['elevation'].centroid_values
1427	#print domain.quantities['stage'].centroid_values
1428
1429	E = domain.quantities['elevation'].vertex_values
1430	L = domain.quantities['stage'].vertex_values
1431
1432	#print E
1433	domain.order = 1
1434	domain.distribute_to_vertices_and_edges()
1435	##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
1436	assert allclose(L[1], [4.1, 16.1, 20.1])
1437
1438	domain.order = 2
1439	domain.distribute_to_vertices_and_edges()
1440	assert allclose(L[1], [4.1, 16.1, 20.1])
1441
1442
1443
1444	def test_second_order_distribute_real_data(self):
1445	#Using test data generated by pyvolution-2
1446	#Assuming no friction and flat bed (0.0)
1447
1448	a = [0.0, 0.0]
1449	b = [0.0, 1.0/5]
1450	c = [0.0, 2.0/5]
1451	d = [1.0/5, 0.0]
1452	e = [1.0/5, 1.0/5]
1453	f = [1.0/5, 2.0/5]
1454	g = [2.0/5, 2.0/5]
1455
1456	points = [a, b, c, d, e, f, g]
1457	#bae, efb, cbf, feg
1458	vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
1459
1460	domain = Domain(points, vertices)
1461
1462	def slope(x, y):
1463	return -x/3
1464
1465	domain.set_quantity('elevation', slope)
1466	domain.set_quantity('stage',
1467	[0.01298164, 0.00365611,
1468	0.01440365, -0.0381856437096],
1469	location='centroids')
1470	domain.set_quantity('xmomentum',
1471	[0.00670439, 0.01263789,
1472	0.00647805, 0.0178180740668],
1473	location='centroids')
1474	domain.set_quantity('ymomentum',
1475	[-7.23510980e-004, -6.30413883e-005,
1476	6.30413883e-005, 0.000200907255866],
1477	location='centroids')
1478
1479	E = domain.quantities['elevation'].vertex_values
1480	L = domain.quantities['stage'].vertex_values
1481	X = domain.quantities['xmomentum'].vertex_values
1482	Y = domain.quantities['ymomentum'].vertex_values
1483
1484	#print E
1485	domain.order = 2
1486	domain.beta_h = 0.0 #Use first order in h-limiter
1487	domain.distribute_to_vertices_and_edges()
1488
1489	#print L[1,:]
1490	#print X[1,:]
1491	#print Y[1,:]
1492
1493	assert allclose(L[1,:], [-0.00825735775384,
1494	-0.00801881482869,
1495	0.0272445025825])
1496	assert allclose(X[1,:], [0.0143507718962,
1497	0.0142502147066,
1498	0.00931268339717])
1499	assert allclose(Y[1,:], [-0.000117062180693,
1500	7.94434448109e-005,
1501	-0.000151505429018])
1502
1503
1504
1505	def test_balance_deep_and_shallow(self):
1506	"""Test that balanced limiters preserve conserved quantites.
1507	"""
1508	import copy
1509
1510	a = [0.0, 0.0]
1511	b = [0.0, 2.0]
1512	c = [2.0, 0.0]
1513	d = [0.0, 4.0]
1514	e = [2.0, 2.0]
1515	f = [4.0, 0.0]
1516
1517	points = [a, b, c, d, e, f]
1518
1519	#bac, bce, ecf, dbe
1520	elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1521
1522	mesh = Domain(points, elements)
1523	mesh.check_integrity()
1524
1525	#Create a deliberate overshoot
1526	mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1527	mesh.set_quantity('elevation', 0) #Flat bed
1528	stage = mesh.quantities['stage']
1529
1530	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
1531
1532	#Limit
1533	mesh.distribute_to_vertices_and_edges()
1534
1535	#Assert that quantities are conserved
1536	from Numeric import sum
1537	for k in range(mesh.number_of_elements):
1538	assert allclose (ref_centroid_values[k],
1539	sum(stage.vertex_values[k,:])/3)
1540
1541
1542	#Now try with a non-flat bed - closely hugging initial stage in places
1543	#This will create alphas in the range [0, 0.478260, 1]
1544	mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1545	mesh.set_quantity('elevation', [[0,0,0],
1546	[1.8,1.9,5.9],
1547	[4.6,0,0],
1548	[0,2,4]])
1549	stage = mesh.quantities['stage']
1550
1551	ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
1552	ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
1553
1554	#Limit
1555	mesh.distribute_to_vertices_and_edges()
1556
1557
1558	#Assert that all vertex quantities have changed
1559	for k in range(mesh.number_of_elements):
1560	#print ref_vertex_values[k,:], stage.vertex_values[k,:]
1561	assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
1562	#and assert that quantities are still conserved
1563	from Numeric import sum
1564	for k in range(mesh.number_of_elements):
1565	assert allclose (ref_centroid_values[k],
1566	sum(stage.vertex_values[k,:])/3)
1567
1568
1569	#Also check that Python and C version produce the same
1570	assert allclose (stage.vertex_values,
1571	[[2,2,2],
1572	[1.93333333, 2.03333333, 6.03333333],
1573	[6.93333333, 4.53333333, 4.53333333],
1574	[5.33333333, 5.33333333, 5.33333333]])
1575
1576
1577
1578
1579	def test_conservation_1(self):
1580	"""Test that stage is conserved globally
1581
1582	This one uses a flat bed, reflective bdries and a suitable
1583	initial condition
1584	"""
1585	from mesh_factory import rectangular
1586	from shallow_water import Domain, Reflective_boundary,\
1587	Dirichlet_boundary, Constant_height
1588	from Numeric import array
1589
1590	#Create basic mesh
1591	points, vertices, boundary = rectangular(6, 6)
1592
1593	#Create shallow water domain
1594	domain = Domain(points, vertices, boundary)
1595	domain.smooth = False
1596	domain.default_order=2
1597
1598	#IC
1599	def x_slope(x, y):
1600	return x/3
1601
1602	domain.set_quantity('elevation', 0)
1603	domain.set_quantity('friction', 0)
1604	domain.set_quantity('stage', x_slope)
1605
1606	# Boundary conditions (reflective everywhere)
1607	Br = Reflective_boundary(domain)
1608	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1609
1610	domain.check_integrity()
1611
1612	#domain.visualise = True #If you want to take a sticky beak
1613
1614	initial_volume = domain.quantities['stage'].get_integral()
1615	initial_xmom = domain.quantities['xmomentum'].get_integral()
1616
1617	#print initial_xmom
1618
1619	#Evolution
1620	for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
1621	volume = domain.quantities['stage'].get_integral()
1622	assert allclose (volume, initial_volume)
1623
1624	#I don't believe that the total momentum should be the same
1625	#It starts with zero and ends with zero though
1626	#xmom = domain.quantities['xmomentum'].get_integral()
1627	#print xmom
1628	#assert allclose (xmom, initial_xmom)
1629
1630	os.remove(domain.filename + '.sww')
1631
1632
1633	def test_conservation_2(self):
1634	"""Test that stage is conserved globally
1635
1636	This one uses a slopy bed, reflective bdries and a suitable
1637	initial condition
1638	"""
1639	from mesh_factory import rectangular
1640	from shallow_water import Domain, Reflective_boundary,\
1641	Dirichlet_boundary, Constant_height
1642	from Numeric import array
1643
1644	#Create basic mesh
1645	points, vertices, boundary = rectangular(6, 6)
1646
1647	#Create shallow water domain
1648	domain = Domain(points, vertices, boundary)
1649	domain.smooth = False
1650	domain.default_order=2
1651
1652	#IC
1653	def x_slope(x, y):
1654	return x/3
1655
1656	domain.set_quantity('elevation', x_slope)
1657	domain.set_quantity('friction', 0)
1658	domain.set_quantity('stage', 0.4) #Steady
1659
1660	# Boundary conditions (reflective everywhere)
1661	Br = Reflective_boundary(domain)
1662	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1663
1664	domain.check_integrity()
1665
1666	#domain.visualise = True #If you want to take a sticky beak
1667
1668	initial_volume = domain.quantities['stage'].get_integral()
1669	initial_xmom = domain.quantities['xmomentum'].get_integral()
1670
1671	#print initial_xmom
1672
1673	#Evolution
1674	for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
1675	volume = domain.quantities['stage'].get_integral()
1676	assert allclose (volume, initial_volume)
1677
1678	#FIXME: What would we expect from momentum
1679	#xmom = domain.quantities['xmomentum'].get_integral()
1680	#print xmom
1681	#assert allclose (xmom, initial_xmom)
1682
1683	os.remove(domain.filename + '.sww')
1684
1685	def test_conservation_3(self):
1686	"""Test that stage is conserved globally
1687
1688	This one uses a larger grid, convoluted bed, reflective bdries and a suitable
1689	initial condition
1690	"""
1691	from mesh_factory import rectangular
1692	from shallow_water import Domain, Reflective_boundary,\
1693	Dirichlet_boundary, Constant_height
1694	from Numeric import array
1695
1696	#Create basic mesh
1697	points, vertices, boundary = rectangular(2, 1)
1698
1699	#Create shallow water domain
1700	domain = Domain(points, vertices, boundary)
1701	domain.smooth = False
1702	domain.default_order = 2
1703	domain.beta_h = 0.2
1704	domain.set_quantities_to_be_stored(['stage'])
1705
1706	#IC
1707	def x_slope(x, y):
1708	z = 0*x
1709	for i in range(len(x)):
1710	if x[i] < 0.3:
1711	z[i] = x[i]/3
1712	if 0.3 <= x[i] < 0.5:
1713	z[i] = -0.5
1714	if 0.5 <= x[i] < 0.7:
1715	z[i] = 0.39
1716	if 0.7 <= x[i]:
1717	z[i] = x[i]/3
1718	return z
1719
1720
1721
1722	domain.set_quantity('elevation', x_slope)
1723	domain.set_quantity('friction', 0)
1724	domain.set_quantity('stage', 0.4) #Steady
1725
1726	# Boundary conditions (reflective everywhere)
1727	Br = Reflective_boundary(domain)
1728	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1729
1730	domain.check_integrity()
1731
1732	#domain.visualise = True #If you want to take a sticky beak
1733
1734	initial_volume = domain.quantities['stage'].get_integral()
1735	initial_xmom = domain.quantities['xmomentum'].get_integral()
1736
1737	import copy
1738	ref_centroid_values =\
1739	copy.copy(domain.quantities['stage'].centroid_values)
1740
1741	#print 'ORG', domain.quantities['stage'].centroid_values
1742	domain.distribute_to_vertices_and_edges()
1743
1744
1745	#print domain.quantities['stage'].centroid_values
1746	assert allclose(domain.quantities['stage'].centroid_values,
1747	ref_centroid_values)
1748
1749
1750	#Check that initial limiter doesn't violate cons quan
1751	assert allclose (domain.quantities['stage'].get_integral(),
1752	initial_volume)
1753
1754	#Evolution
1755	for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
1756	volume = domain.quantities['stage'].get_integral()
1757	#print t, volume, initial_volume
1758	assert allclose (volume, initial_volume)
1759
1760	os.remove(domain.filename + '.sww')
1761
1762	def test_conservation_4(self):
1763	"""Test that stage is conserved globally
1764
1765	This one uses a larger grid, convoluted bed, reflective bdries and a suitable
1766	initial condition
1767	"""
1768	from mesh_factory import rectangular
1769	from shallow_water import Domain, Reflective_boundary,\
1770	Dirichlet_boundary, Constant_height
1771	from Numeric import array
1772
1773	#Create basic mesh
1774	points, vertices, boundary = rectangular(6, 6)
1775
1776	#Create shallow water domain
1777	domain = Domain(points, vertices, boundary)
1778	domain.smooth = False
1779	domain.default_order=2
1780	domain.beta_h = 0.0
1781
1782	#IC
1783	def x_slope(x, y):
1784	z = 0*x
1785	for i in range(len(x)):
1786	if x[i] < 0.3:
1787	z[i] = x[i]/3
1788	if 0.3 <= x[i] < 0.5:
1789	z[i] = -0.5
1790	if 0.5 <= x[i] < 0.7:
1791	#z[i] = 0.3 #OK with beta == 0.2
1792	z[i] = 0.34 #OK with beta == 0.0
1793	#z[i] = 0.35#Fails after 80 timesteps with an error
1794	#of the order 1.0e-5
1795	if 0.7 <= x[i]:
1796	z[i] = x[i]/3
1797	return z
1798
1799
1800
1801	domain.set_quantity('elevation', x_slope)
1802	domain.set_quantity('friction', 0)
1803	domain.set_quantity('stage', 0.4) #Steady
1804
1805	# Boundary conditions (reflective everywhere)
1806	Br = Reflective_boundary(domain)
1807	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1808
1809	domain.check_integrity()
1810
1811	#domain.visualise = True #If you want to take a sticky beak
1812
1813	initial_volume = domain.quantities['stage'].get_integral()
1814	initial_xmom = domain.quantities['xmomentum'].get_integral()
1815
1816	import copy
1817	ref_centroid_values =\
1818	copy.copy(domain.quantities['stage'].centroid_values)
1819
1820	#Test limiter by itself
1821	domain.distribute_to_vertices_and_edges()
1822
1823	#Check that initial limiter doesn't violate cons quan
1824	assert allclose (domain.quantities['stage'].get_integral(),
1825	initial_volume)
1826	#NOTE: This would fail if any initial stage was less than the
1827	#corresponding bed elevation - but that is reasonable.
1828
1829
1830	#Evolution
1831	for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
1832	volume = domain.quantities['stage'].get_integral()
1833
1834	#print t, volume, initial_volume
1835
1836
1837	#if not allclose (volume, initial_volume):
1838	# print 't==4.05'
1839	# for k in range(domain.number_of_elements):
1840	# pass
1841	# print domain.quantities['stage'].centroid_values[k] -\
1842	# ref_centroid_values[k]
1843
1844	assert allclose (volume, initial_volume)
1845
1846
1847	os.remove(domain.filename + '.sww')
1848
1849
1850	def test_conservation_5(self):
1851	"""Test that momentum is conserved globally in
1852	steady state scenario
1853
1854	This one uses a slopy bed, dirichlet and reflective bdries
1855	"""
1856	from mesh_factory import rectangular
1857	from shallow_water import Domain, Reflective_boundary,\
1858	Dirichlet_boundary, Constant_height
1859	from Numeric import array
1860
1861	#Create basic mesh
1862	points, vertices, boundary = rectangular(6, 6)
1863
1864	#Create shallow water domain
1865	domain = Domain(points, vertices, boundary)
1866	domain.smooth = False
1867	domain.default_order=2
1868
1869	#IC
1870	def x_slope(x, y):
1871	return x/3
1872
1873	domain.set_quantity('elevation', x_slope)
1874	domain.set_quantity('friction', 0)
1875	domain.set_quantity('stage', 0.4) #Steady
1876
1877	# Boundary conditions (reflective everywhere)
1878	Br = Reflective_boundary(domain)
1879	Bleft = Dirichlet_boundary([0.5,0,0])
1880	Bright = Dirichlet_boundary([0.1,0,0])
1881	domain.set_boundary({'left': Bleft, 'right': Bright,
1882	'top': Br, 'bottom': Br})
1883
1884	domain.check_integrity()
1885
1886	#domain.visualise = True #If you want to take a sticky beak
1887
1888	initial_volume = domain.quantities['stage'].get_integral()
1889	initial_xmom = domain.quantities['xmomentum'].get_integral()
1890
1891
1892	#Evolution
1893	for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
1894	stage = domain.quantities['stage'].get_integral()
1895	xmom = domain.quantities['xmomentum'].get_integral()
1896	ymom = domain.quantities['ymomentum'].get_integral()
1897
1898	if allclose(t, 6): #Steady state reached
1899	steady_xmom = domain.quantities['xmomentum'].get_integral()
1900	steady_ymom = domain.quantities['ymomentum'].get_integral()
1901	steady_stage = domain.quantities['stage'].get_integral()
1902
1903	if t > 6:
1904	#print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
1905	assert allclose(xmom, steady_xmom)
1906	assert allclose(ymom, steady_ymom)
1907	assert allclose(stage, steady_stage)
1908
1909
1910	os.remove(domain.filename + '.sww')
1911
1912
1913
1914
1915
1916	def test_conservation_real(self):
1917	"""Test that momentum is conserved globally
1918
1919	Stephen finally made a test that revealed the problem.
1920	This test failed with code prior to 25 July 2005
1921	"""
1922
1923	yieldstep = 0.01
1924	finaltime = 0.05
1925	min_depth = 1.0e-2
1926
1927
1928	import sys
1929	from os import sep; sys.path.append('..'+sep+'pyvolution')
1930	from mesh_factory import rectangular
1931	from shallow_water import Domain, Reflective_boundary
1932
1933
1934	#Create shallow water domain
1935	points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
1936	domain = Domain(points, vertices, boundary)
1937	domain.smooth = False
1938	domain.visualise = False
1939	domain.default_order = 1
1940	domain.minimum_allowed_height = min_depth
1941
1942	# Set initial condition
1943	class Set_IC:
1944	"""Set an initial condition with a constant value, for x0<x<x1
1945	"""
1946
1947	def __init__(self, x0=0.25, x1=0.5, h=1.0):
1948	self.x0 = x0
1949	self.x1 = x1
1950	self.h = h
1951
1952	def __call__(self, x, y):
1953	return self.h*((x>self.x0)&(x<self.x1))
1954
1955
1956	domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
1957
1958
1959	#Boundaries
1960	R = Reflective_boundary(domain)
1961	domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R})
1962
1963	ref = domain.quantities['stage'].get_integral()
1964
1965	# Evolution
1966	for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
1967	pass
1968	#print 'Integral stage = ',\
1969	# domain.quantities['stage'].get_integral(),\
1970	# ' Time = ',domain.time
1971
1972
1973	now = domain.quantities['stage'].get_integral()
1974
1975	msg = 'Stage not conserved: was %f, now %f' %(ref, now)
1976	assert allclose(ref, now), msg
1977
1978	os.remove(domain.filename + '.sww')
1979
1980	def test_second_order_flat_bed_onestep(self):
1981
1982	from mesh_factory import rectangular
1983	from shallow_water import Domain, Reflective_boundary,\
1984	Dirichlet_boundary, Constant_height
1985	from Numeric import array
1986
1987	#Create basic mesh
1988	points, vertices, boundary = rectangular(6, 6)
1989
1990	#Create shallow water domain
1991	domain = Domain(points, vertices, boundary)
1992	domain.smooth = False
1993	domain.default_order=2
1994
1995	# Boundary conditions
1996	Br = Reflective_boundary(domain)
1997	Bd = Dirichlet_boundary([0.1, 0., 0.])
1998	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1999
2000	domain.check_integrity()
2001
2002	#Evolution
2003	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2004	pass# domain.write_time()
2005
2006	#Data from earlier version of pyvolution
2007	assert allclose(domain.min_timestep, 0.0396825396825)
2008	assert allclose(domain.max_timestep, 0.0396825396825)
2009
2010	assert allclose(domain.quantities['stage'].centroid_values[:12],
2011	[0.00171396, 0.02656103, 0.00241523, 0.02656103,
2012	0.00241523, 0.02656103, 0.00241523, 0.02656103,
2013	0.00241523, 0.02656103, 0.00241523, 0.0272623])
2014
2015	domain.distribute_to_vertices_and_edges()
2016	assert allclose(domain.quantities['stage'].vertex_values[:12,0],
2017	[0.0001714, 0.02656103, 0.00024152,
2018	0.02656103, 0.00024152, 0.02656103,
2019	0.00024152, 0.02656103, 0.00024152,
2020	0.02656103, 0.00024152, 0.0272623])
2021
2022	assert allclose(domain.quantities['stage'].vertex_values[:12,1],
2023	[0.00315012, 0.02656103, 0.00024152, 0.02656103,
2024	0.00024152, 0.02656103, 0.00024152, 0.02656103,
2025	0.00024152, 0.02656103, 0.00040506, 0.0272623])
2026
2027	assert allclose(domain.quantities['stage'].vertex_values[:12,2],
2028	[0.00182037, 0.02656103, 0.00676264,
2029	0.02656103, 0.00676264, 0.02656103,
2030	0.00676264, 0.02656103, 0.00676264,
2031	0.02656103, 0.0065991, 0.0272623])
2032
2033	assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
2034	[0.00113961, 0.01302432, 0.00148672,
2035	0.01302432, 0.00148672, 0.01302432,
2036	0.00148672, 0.01302432, 0.00148672 ,
2037	0.01302432, 0.00148672, 0.01337143])
2038
2039	assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
2040	[-2.91240050e-004 , 1.22721531e-004,
2041	-1.22721531e-004, 1.22721531e-004 ,
2042	-1.22721531e-004, 1.22721531e-004,
2043	-1.22721531e-004 , 1.22721531e-004,
2044	-1.22721531e-004, 1.22721531e-004,
2045	-1.22721531e-004, -4.57969873e-005])
2046
2047	os.remove(domain.filename + '.sww')
2048
2049
2050	def test_second_order_flat_bed_moresteps(self):
2051
2052	from mesh_factory import rectangular
2053	from shallow_water import Domain, Reflective_boundary,\
2054	Dirichlet_boundary, Constant_height
2055	from Numeric import array
2056
2057	#Create basic mesh
2058	points, vertices, boundary = rectangular(6, 6)
2059
2060	#Create shallow water domain
2061	domain = Domain(points, vertices, boundary)
2062	domain.smooth = False
2063	domain.default_order=2
2064
2065	# Boundary conditions
2066	Br = Reflective_boundary(domain)
2067	Bd = Dirichlet_boundary([0.1, 0., 0.])
2068	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2069
2070	domain.check_integrity()
2071
2072	#Evolution
2073	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2074	pass
2075
2076	#Data from earlier version of pyvolution
2077	#assert allclose(domain.min_timestep, 0.0396825396825)
2078	#assert allclose(domain.max_timestep, 0.0396825396825)
2079	#print domain.quantities['stage'].centroid_values
2080
2081	os.remove(domain.filename + '.sww')
2082
2083
2084	def test_flatbed_first_order(self):
2085	from mesh_factory import rectangular
2086	from shallow_water import Domain,\
2087	Reflective_boundary, Dirichlet_boundary
2088
2089	from Numeric import array
2090
2091	#Create basic mesh
2092	N = 8
2093	points, vertices, boundary = rectangular(N, N)
2094
2095	#Create shallow water domain
2096	domain = Domain(points, vertices, boundary)
2097	domain.smooth = False
2098	domain.visualise = False
2099	domain.default_order=1
2100
2101	# Boundary conditions
2102	Br = Reflective_boundary(domain)
2103	Bd = Dirichlet_boundary([0.2,0.,0.])
2104
2105	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2106	domain.check_integrity()
2107
2108
2109	#Evolution
2110	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
2111	pass
2112	#domain.write_time()
2113
2114	#FIXME: These numbers were from version before 25/10
2115	#assert allclose(domain.min_timestep, 0.0140413643926)
2116	#assert allclose(domain.max_timestep, 0.0140947355753)
2117
2118	for i in range(3):
2119	#assert allclose(domain.quantities['stage'].edge_values[:4,i],
2120	# [0.10730244,0.12337617,0.11200126,0.12605666])
2121
2122	assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
2123	[0.07610894,0.06901572,0.07284461,0.06819712])
2124
2125	#assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
2126	# [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
2127
2128
2129	os.remove(domain.filename + '.sww')
2130
2131	def test_flatbed_second_order(self):
2132	from mesh_factory import rectangular
2133	from shallow_water import Domain,\
2134	Reflective_boundary, Dirichlet_boundary
2135
2136	from Numeric import array
2137
2138	#Create basic mesh
2139	N = 8
2140	points, vertices, boundary = rectangular(N, N)
2141
2142	#Create shallow water domain
2143	domain = Domain(points, vertices, boundary)
2144	domain.smooth = False
2145	domain.visualise = False
2146	domain.default_order=2
2147	#domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
2148
2149	# Boundary conditions
2150	Br = Reflective_boundary(domain)
2151	Bd = Dirichlet_boundary([0.2,0.,0.])
2152
2153	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2154	domain.check_integrity()
2155
2156	#Evolution
2157	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2158	pass
2159
2160
2161	assert allclose(domain.min_timestep, 0.0210448446782)
2162	assert allclose(domain.max_timestep, 0.0210448446782)
2163
2164	#print domain.quantities['stage'].vertex_values[:4,0]
2165	#print domain.quantities['xmomentum'].vertex_values[:4,0]
2166	#print domain.quantities['ymomentum'].vertex_values[:4,0]
2167
2168	#FIXME: These numbers were from version before 25/10
2169	#assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2170	# [0.00101913,0.05352143,0.00104852,0.05354394])
2171
2172	#FIXME: These numbers were from version before 21/3/6 -
2173	#could be recreated by setting maximum_allowed_speed to 0 maybe
2174	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2175	# [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
2176
2177	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2178	[ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
2179
2180
2181
2182	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2183	# [0.00090581,0.03685719,0.00088303,0.03687313])
2184
2185	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2186	[-0.00139463,0.0006156,-0.00060364,0.00061827])
2187
2188
2189	os.remove(domain.filename + '.sww')
2190
2191
2192	def test_flatbed_second_order_vmax_0(self):
2193	from mesh_factory import rectangular
2194	from shallow_water import Domain,\
2195	Reflective_boundary, Dirichlet_boundary
2196
2197	from Numeric import array
2198
2199	#Create basic mesh
2200	N = 8
2201	points, vertices, boundary = rectangular(N, N)
2202
2203	#Create shallow water domain
2204	domain = Domain(points, vertices, boundary)
2205	domain.smooth = False
2206	domain.visualise = False
2207	domain.default_order=2
2208	domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
2209
2210	# Boundary conditions
2211	Br = Reflective_boundary(domain)
2212	Bd = Dirichlet_boundary([0.2,0.,0.])
2213
2214	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2215	domain.check_integrity()
2216
2217	#Evolution
2218	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2219	pass
2220
2221
2222	assert allclose(domain.min_timestep, 0.0210448446782)
2223	assert allclose(domain.max_timestep, 0.0210448446782)
2224
2225	#FIXME: These numbers were from version before 21/3/6 -
2226	#could be recreated by setting maximum_allowed_speed to 0 maybe
2227	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2228	[ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
2229
2230
2231	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2232	[-0.00139463,0.0006156,-0.00060364,0.00061827])
2233
2234
2235	os.remove(domain.filename + '.sww')
2236
2237
2238
2239	def test_flatbed_second_order_distribute(self):
2240	#Use real data from pyvolution 2
2241	#painfully setup and extracted.
2242	from mesh_factory import rectangular
2243	from shallow_water import Domain,\
2244	Reflective_boundary, Dirichlet_boundary
2245
2246	from Numeric import array
2247
2248	#Create basic mesh
2249	N = 8
2250	points, vertices, boundary = rectangular(N, N)
2251
2252	#Create shallow water domain
2253	domain = Domain(points, vertices, boundary)
2254	domain.smooth = False
2255	domain.visualise = False
2256	domain.default_order=domain.order=2
2257
2258	# Boundary conditions
2259	Br = Reflective_boundary(domain)
2260	Bd = Dirichlet_boundary([0.2,0.,0.])
2261
2262	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2263	domain.check_integrity()
2264
2265
2266
2267	for V in [False, True]:
2268	if V:
2269	#Set centroids as if system had been evolved
2270	L = zeros(2NN, Float)
2271	L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
2272	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2273	1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
2274	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
2275	1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
2276	5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
2277	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2278	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2279	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
2280	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
2281	0.00000000e+000, 5.57305948e-005]
2282
2283	X = zeros(2NN, Float)
2284	X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
2285	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2286	8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
2287	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
2288	8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
2289	3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
2290	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2291	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2292	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
2293	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
2294	0.00000000e+000, 4.57662812e-005]
2295
2296	Y = zeros(2NN, Float)
2297	Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
2298	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2299	-6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
2300	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
2301	-6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
2302	-1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
2303	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2304	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2305	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
2306	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
2307	0.00000000e+000, -2.57635178e-005]
2308
2309
2310	domain.set_quantity('stage', L, location='centroids')
2311	domain.set_quantity('xmomentum', X, location='centroids')
2312	domain.set_quantity('ymomentum', Y, location='centroids')
2313
2314	domain.check_integrity()
2315	else:
2316	#Evolution
2317	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
2318	pass
2319	assert allclose(domain.min_timestep, 0.0210448446782)
2320	assert allclose(domain.max_timestep, 0.0210448446782)
2321
2322
2323	#Centroids were correct but not vertices.
2324	#Hence the check of distribute below.
2325	assert allclose(domain.quantities['stage'].centroid_values[:4],
2326	[0.00721206,0.05352143,0.01009108,0.05354394])
2327
2328	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
2329	[0.00648352,0.03685719,0.00850733,0.03687313])
2330
2331	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2332	[-0.00139463,0.0006156,-0.00060364,0.00061827])
2333
2334	#print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
2335	#print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
2336
2337	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
2338	##print domain.quantities['xmomentum'].centroid_values[17], V
2339	##print
2340	if not V:
2341	#FIXME: These numbers were from version before 21/3/6 -
2342	#could be recreated by setting maximum_allowed_speed to 0 maybe
2343
2344	#assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
2345	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)
2346
2347	else:
2348	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
2349
2350	import copy
2351	XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
2352	assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2353
2354	domain.distribute_to_vertices_and_edges()
2355
2356	#assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
2357
2358	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],
2359	# 0.0)
2360	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.000286060839592)
2361
2362
2363	#FIXME: These numbers were from version before 25/10
2364	#assert allclose(domain.quantities['stage'].vertex_values[:4,0],
2365	# [0.00101913,0.05352143,0.00104852,0.05354394])
2366
2367	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
2368	[-0.00139463,0.0006156,-0.00060364,0.00061827])
2369
2370
2371	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2372	[0.00090581,0.03685719,0.00088303,0.03687313])
2373
2374
2375	#NB NO longer relvant:
2376
2377	#This was the culprit. First triangles vertex 0 had an
2378	#x-momentum of 0.0064835 instead of 0.00090581 and
2379	#third triangle had 0.00850733 instead of 0.00088303
2380	#print domain.quantities['xmomentum'].vertex_values[:4,0]
2381
2382	#print domain.quantities['xmomentum'].vertex_values[:4,0]
2383	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
2384	# [0.00090581,0.03685719,0.00088303,0.03687313])
2385
2386	os.remove(domain.filename + '.sww')
2387
2388
2389
2390	def test_bedslope_problem_first_order(self):
2391
2392	from mesh_factory import rectangular
2393	from shallow_water import Domain, Reflective_boundary, Constant_height
2394	from Numeric import array
2395
2396	#Create basic mesh
2397	points, vertices, boundary = rectangular(6, 6)
2398
2399	#Create shallow water domain
2400	domain = Domain(points, vertices, boundary)
2401	domain.smooth = False
2402	domain.default_order=1
2403
2404	#Bed-slope and friction
2405	def x_slope(x, y):
2406	return -x/3
2407
2408	domain.set_quantity('elevation', x_slope)
2409
2410	# Boundary conditions
2411	Br = Reflective_boundary(domain)
2412	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2413
2414	#Initial condition
2415	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2416	domain.check_integrity()
2417
2418	#Evolution
2419	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2420	pass# domain.write_time()
2421
2422	assert allclose(domain.min_timestep, 0.050010003001)
2423	assert allclose(domain.max_timestep, 0.050010003001)
2424
2425
2426	os.remove(domain.filename + '.sww')
2427
2428	def test_bedslope_problem_first_order_moresteps(self):
2429
2430	from mesh_factory import rectangular
2431	from shallow_water import Domain, Reflective_boundary, Constant_height
2432	from Numeric import array
2433
2434	#Create basic mesh
2435	points, vertices, boundary = rectangular(6, 6)
2436
2437	#Create shallow water domain
2438	domain = Domain(points, vertices, boundary)
2439	domain.smooth = False
2440	domain.default_order=1
2441	domain.beta_h = 0.0 #Use first order in h-limiter
2442
2443	#Bed-slope and friction
2444	def x_slope(x, y):
2445	return -x/3
2446
2447	domain.set_quantity('elevation', x_slope)
2448
2449	# Boundary conditions
2450	Br = Reflective_boundary(domain)
2451	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2452
2453	#Initial condition
2454	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2455	domain.check_integrity()
2456
2457	#Evolution
2458	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2459	pass# domain.write_time()
2460
2461	#Data from earlier version of pyvolution
2462	#print domain.quantities['stage'].centroid_values
2463
2464	assert allclose(domain.quantities['stage'].centroid_values,
2465	[-0.02998628, -0.01520652, -0.03043492,
2466	-0.0149132, -0.03004706, -0.01476251,
2467	-0.0298215, -0.01467976, -0.02988158,
2468	-0.01474662, -0.03036161, -0.01442995,
2469	-0.07624583, -0.06297061, -0.07733792,
2470	-0.06342237, -0.07695439, -0.06289595,
2471	-0.07635559, -0.0626065, -0.07633628,
2472	-0.06280072, -0.07739632, -0.06386738,
2473	-0.12161738, -0.11028239, -0.1223796,
2474	-0.11095953, -0.12189744, -0.11048616,
2475	-0.12074535, -0.10987605, -0.12014311,
2476	-0.10976691, -0.12096859, -0.11087692,
2477	-0.16868259, -0.15868061, -0.16801135,
2478	-0.1588003, -0.16674343, -0.15813323,
2479	-0.16457595, -0.15693826, -0.16281096,
2480	-0.15585154, -0.16283873, -0.15540068,
2481	-0.17450362, -0.19919913, -0.18062882,
2482	-0.19764131, -0.17783111, -0.19407213,
2483	-0.1736915, -0.19053624, -0.17228678,
2484	-0.19105634, -0.17920133, -0.1968828,
2485	-0.14244395, -0.14604641, -0.14473537,
2486	-0.1506107, -0.14510055, -0.14919522,
2487	-0.14175896, -0.14560798, -0.13911658,
2488	-0.14439383, -0.13924047, -0.14829043])
2489
2490	os.remove(domain.filename + '.sww')
2491
2492	def test_bedslope_problem_second_order_one_step(self):
2493
2494	from mesh_factory import rectangular
2495	from shallow_water import Domain, Reflective_boundary, Constant_height
2496	from Numeric import array
2497
2498	#Create basic mesh
2499	points, vertices, boundary = rectangular(6, 6)
2500
2501	#Create shallow water domain
2502	domain = Domain(points, vertices, boundary)
2503	domain.smooth = False
2504	domain.default_order=2
2505
2506	#Bed-slope and friction at vertices (and interpolated elsewhere)
2507	def x_slope(x, y):
2508	return -x/3
2509
2510	domain.set_quantity('elevation', x_slope)
2511
2512	# Boundary conditions
2513	Br = Reflective_boundary(domain)
2514	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2515
2516	#Initial condition
2517	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2518	domain.check_integrity()
2519
2520	assert allclose(domain.quantities['stage'].centroid_values,
2521	[0.01296296, 0.03148148, 0.01296296,
2522	0.03148148, 0.01296296, 0.03148148,
2523	0.01296296, 0.03148148, 0.01296296,
2524	0.03148148, 0.01296296, 0.03148148,
2525	-0.04259259, -0.02407407, -0.04259259,
2526	-0.02407407, -0.04259259, -0.02407407,
2527	-0.04259259, -0.02407407, -0.04259259,
2528	-0.02407407, -0.04259259, -0.02407407,
2529	-0.09814815, -0.07962963, -0.09814815,
2530	-0.07962963, -0.09814815, -0.07962963,
2531	-0.09814815, -0.07962963, -0.09814815,
2532	-0.07962963, -0.09814815, -0.07962963,
2533	-0.1537037 , -0.13518519, -0.1537037,
2534	-0.13518519, -0.1537037, -0.13518519,
2535	-0.1537037 , -0.13518519, -0.1537037,
2536	-0.13518519, -0.1537037, -0.13518519,
2537	-0.20925926, -0.19074074, -0.20925926,
2538	-0.19074074, -0.20925926, -0.19074074,
2539	-0.20925926, -0.19074074, -0.20925926,
2540	-0.19074074, -0.20925926, -0.19074074,
2541	-0.26481481, -0.2462963, -0.26481481,
2542	-0.2462963, -0.26481481, -0.2462963,
2543	-0.26481481, -0.2462963, -0.26481481,
2544	-0.2462963, -0.26481481, -0.2462963])
2545
2546
2547	#print domain.quantities['stage'].extrapolate_second_order()
2548	#domain.distribute_to_vertices_and_edges()
2549	#print domain.quantities['stage'].vertex_values[:,0]
2550
2551	#Evolution
2552	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
2553	#domain.write_time()
2554	pass
2555
2556
2557	#print domain.quantities['stage'].centroid_values
2558	assert allclose(domain.quantities['stage'].centroid_values,
2559	[0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
2560	0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
2561	0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
2562	-0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
2563	-0.04186556, -0.0248011, -0.04186556, -0.02255593,
2564	-0.09966629, -0.08035666, -0.09742112, -0.08035666,
2565	-0.09742112, -0.08035666, -0.09742112, -0.08035666,
2566	-0.09742112, -0.08035666, -0.09742112, -0.07811149,
2567	-0.15522185, -0.13591222, -0.15297667, -0.13591222,
2568	-0.15297667, -0.13591222, -0.15297667, -0.13591222,
2569	-0.15297667, -0.13591222, -0.15297667, -0.13366704,
2570	-0.2107774, -0.19146777, -0.20853223, -0.19146777,
2571	-0.20853223, -0.19146777, -0.20853223, -0.19146777,
2572	-0.20853223, -0.19146777, -0.20853223, -0.1892226,
2573	-0.26120669, -0.24776246, -0.25865535, -0.24776246,
2574	-0.25865535, -0.24776246, -0.25865535, -0.24776246,
2575	-0.25865535, -0.24776246, -0.25865535, -0.24521113])
2576
2577	os.remove(domain.filename + '.sww')
2578
2579	def test_bedslope_problem_second_order_two_steps(self):
2580
2581	from mesh_factory import rectangular
2582	from shallow_water import Domain, Reflective_boundary, Constant_height
2583	from Numeric import array
2584
2585	#Create basic mesh
2586	points, vertices, boundary = rectangular(6, 6)
2587
2588	#Create shallow water domain
2589	domain = Domain(points, vertices, boundary)
2590	domain.smooth = False
2591	domain.default_order=2
2592	domain.beta_h = 0.0 #Use first order in h-limiter
2593
2594	#Bed-slope and friction at vertices (and interpolated elsewhere)
2595	def x_slope(x, y):
2596	return -x/3
2597
2598	domain.set_quantity('elevation', x_slope)
2599
2600	# Boundary conditions
2601	Br = Reflective_boundary(domain)
2602	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2603
2604	#Initial condition
2605	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2606	domain.check_integrity()
2607
2608	assert allclose(domain.quantities['stage'].centroid_values,
2609	[0.01296296, 0.03148148, 0.01296296,
2610	0.03148148, 0.01296296, 0.03148148,
2611	0.01296296, 0.03148148, 0.01296296,
2612	0.03148148, 0.01296296, 0.03148148,
2613	-0.04259259, -0.02407407, -0.04259259,
2614	-0.02407407, -0.04259259, -0.02407407,
2615	-0.04259259, -0.02407407, -0.04259259,
2616	-0.02407407, -0.04259259, -0.02407407,
2617	-0.09814815, -0.07962963, -0.09814815,
2618	-0.07962963, -0.09814815, -0.07962963,
2619	-0.09814815, -0.07962963, -0.09814815,
2620	-0.07962963, -0.09814815, -0.07962963,
2621	-0.1537037 , -0.13518519, -0.1537037,
2622	-0.13518519, -0.1537037, -0.13518519,
2623	-0.1537037 , -0.13518519, -0.1537037,
2624	-0.13518519, -0.1537037, -0.13518519,
2625	-0.20925926, -0.19074074, -0.20925926,
2626	-0.19074074, -0.20925926, -0.19074074,
2627	-0.20925926, -0.19074074, -0.20925926,
2628	-0.19074074, -0.20925926, -0.19074074,
2629	-0.26481481, -0.2462963, -0.26481481,
2630	-0.2462963, -0.26481481, -0.2462963,
2631	-0.26481481, -0.2462963, -0.26481481,
2632	-0.2462963, -0.26481481, -0.2462963])
2633
2634
2635	#print domain.quantities['stage'].extrapolate_second_order()
2636	#domain.distribute_to_vertices_and_edges()
2637	#print domain.quantities['stage'].vertex_values[:,0]
2638
2639	#Evolution
2640	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2641	pass
2642
2643
2644	#Data from earlier version of pyvolution ft=0.1
2645	assert allclose(domain.min_timestep, 0.0376895634803)
2646	assert allclose(domain.max_timestep, 0.0415635655309)
2647
2648
2649	assert allclose(domain.quantities['stage'].centroid_values,
2650	[0.00855788, 0.01575204, 0.00994606, 0.01717072,
2651	0.01005985, 0.01716362, 0.01005985, 0.01716299,
2652	0.01007098, 0.01736248, 0.01216452, 0.02026776,
2653	-0.04462374, -0.02479045, -0.04199789, -0.0229465,
2654	-0.04184033, -0.02295693, -0.04184013, -0.02295675,
2655	-0.04184486, -0.0228168, -0.04028876, -0.02036486,
2656	-0.10029444, -0.08170809, -0.09772846, -0.08021704,
2657	-0.09760006, -0.08022143, -0.09759984, -0.08022124,
2658	-0.09760261, -0.08008893, -0.09603914, -0.07758209,
2659	-0.15584152, -0.13723138, -0.15327266, -0.13572906,
2660	-0.15314427, -0.13573349, -0.15314405, -0.13573331,
2661	-0.15314679, -0.13560104, -0.15158523, -0.13310701,
2662	-0.21208605, -0.19283913, -0.20955631, -0.19134189,
2663	-0.20942821, -0.19134598, -0.20942799, -0.1913458,
2664	-0.20943005, -0.19120952, -0.20781177, -0.18869401,
2665	-0.25384082, -0.2463294, -0.25047649, -0.24464654,
2666	-0.25031159, -0.24464253, -0.25031112, -0.24464253,
2667	-0.25031463, -0.24454764, -0.24885323, -0.24286438])
2668
2669
2670	os.remove(domain.filename + '.sww')
2671
2672	def test_bedslope_problem_second_order_two_yieldsteps(self):
2673
2674	from mesh_factory import rectangular
2675	from shallow_water import Domain, Reflective_boundary, Constant_height
2676	from Numeric import array
2677
2678	#Create basic mesh
2679	points, vertices, boundary = rectangular(6, 6)
2680
2681	#Create shallow water domain
2682	domain = Domain(points, vertices, boundary)
2683	domain.smooth = False
2684	domain.default_order=2
2685	domain.beta_h = 0.0 #Use first order in h-limiter
2686
2687	#Bed-slope and friction at vertices (and interpolated elsewhere)
2688	def x_slope(x, y):
2689	return -x/3
2690
2691	domain.set_quantity('elevation', x_slope)
2692
2693	# Boundary conditions
2694	Br = Reflective_boundary(domain)
2695	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2696
2697	#Initial condition
2698	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2699	domain.check_integrity()
2700
2701	assert allclose(domain.quantities['stage'].centroid_values,
2702	[0.01296296, 0.03148148, 0.01296296,
2703	0.03148148, 0.01296296, 0.03148148,
2704	0.01296296, 0.03148148, 0.01296296,
2705	0.03148148, 0.01296296, 0.03148148,
2706	-0.04259259, -0.02407407, -0.04259259,
2707	-0.02407407, -0.04259259, -0.02407407,
2708	-0.04259259, -0.02407407, -0.04259259,
2709	-0.02407407, -0.04259259, -0.02407407,
2710	-0.09814815, -0.07962963, -0.09814815,
2711	-0.07962963, -0.09814815, -0.07962963,
2712	-0.09814815, -0.07962963, -0.09814815,
2713	-0.07962963, -0.09814815, -0.07962963,
2714	-0.1537037 , -0.13518519, -0.1537037,
2715	-0.13518519, -0.1537037, -0.13518519,
2716	-0.1537037 , -0.13518519, -0.1537037,
2717	-0.13518519, -0.1537037, -0.13518519,
2718	-0.20925926, -0.19074074, -0.20925926,
2719	-0.19074074, -0.20925926, -0.19074074,
2720	-0.20925926, -0.19074074, -0.20925926,
2721	-0.19074074, -0.20925926, -0.19074074,
2722	-0.26481481, -0.2462963, -0.26481481,
2723	-0.2462963, -0.26481481, -0.2462963,
2724	-0.26481481, -0.2462963, -0.26481481,
2725	-0.2462963, -0.26481481, -0.2462963])
2726
2727
2728	#print domain.quantities['stage'].extrapolate_second_order()
2729	#domain.distribute_to_vertices_and_edges()
2730	#print domain.quantities['stage'].vertex_values[:,0]
2731
2732	#Evolution
2733	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): #0.05??
2734	#domain.write_time()
2735	pass
2736
2737
2738
2739	assert allclose(domain.quantities['stage'].centroid_values,
2740	[0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
2741	0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
2742	0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
2743	-0.0229465, -0.04184033, -0.02295693, -0.04184013,
2744	-0.02295675, -0.04184486, -0.0228168, -0.04028876,
2745	-0.02036486, -0.10029444, -0.08170809, -0.09772846,
2746	-0.08021704, -0.09760006, -0.08022143, -0.09759984,
2747	-0.08022124, -0.09760261, -0.08008893, -0.09603914,
2748	-0.07758209, -0.15584152, -0.13723138, -0.15327266,
2749	-0.13572906, -0.15314427, -0.13573349, -0.15314405,
2750	-0.13573331, -0.15314679, -0.13560104, -0.15158523,
2751	-0.13310701, -0.21208605, -0.19283913, -0.20955631,
2752	-0.19134189, -0.20942821, -0.19134598, -0.20942799,
2753	-0.1913458, -0.20943005, -0.19120952, -0.20781177,
2754	-0.18869401, -0.25384082, -0.2463294, -0.25047649,
2755	-0.24464654, -0.25031159, -0.24464253, -0.25031112,
2756	-0.24464253, -0.25031463, -0.24454764, -0.24885323,
2757	-0.24286438])
2758
2759	os.remove(domain.filename + '.sww')
2760
2761	def test_bedslope_problem_second_order_more_steps(self):
2762
2763	from mesh_factory import rectangular
2764	from shallow_water import Domain, Reflective_boundary, Constant_height
2765	from Numeric import array
2766
2767	#Create basic mesh
2768	points, vertices, boundary = rectangular(6, 6)
2769
2770	#Create shallow water domain
2771	domain = Domain(points, vertices, boundary)
2772	domain.smooth = False
2773	domain.default_order=2
2774	domain.beta_h = 0.0 #Use first order in h-limiter
2775
2776	#Bed-slope and friction at vertices (and interpolated elsewhere)
2777	def x_slope(x, y):
2778	return -x/3
2779
2780	domain.set_quantity('elevation', x_slope)
2781
2782	# Boundary conditions
2783	Br = Reflective_boundary(domain)
2784	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2785
2786	#Initial condition
2787	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2788	domain.check_integrity()
2789
2790	assert allclose(domain.quantities['stage'].centroid_values,
2791	[0.01296296, 0.03148148, 0.01296296,
2792	0.03148148, 0.01296296, 0.03148148,
2793	0.01296296, 0.03148148, 0.01296296,
2794	0.03148148, 0.01296296, 0.03148148,
2795	-0.04259259, -0.02407407, -0.04259259,
2796	-0.02407407, -0.04259259, -0.02407407,
2797	-0.04259259, -0.02407407, -0.04259259,
2798	-0.02407407, -0.04259259, -0.02407407,
2799	-0.09814815, -0.07962963, -0.09814815,
2800	-0.07962963, -0.09814815, -0.07962963,
2801	-0.09814815, -0.07962963, -0.09814815,
2802	-0.07962963, -0.09814815, -0.07962963,
2803	-0.1537037 , -0.13518519, -0.1537037,
2804	-0.13518519, -0.1537037, -0.13518519,
2805	-0.1537037 , -0.13518519, -0.1537037,
2806	-0.13518519, -0.1537037, -0.13518519,
2807	-0.20925926, -0.19074074, -0.20925926,
2808	-0.19074074, -0.20925926, -0.19074074,
2809	-0.20925926, -0.19074074, -0.20925926,
2810	-0.19074074, -0.20925926, -0.19074074,
2811	-0.26481481, -0.2462963, -0.26481481,
2812	-0.2462963, -0.26481481, -0.2462963,
2813	-0.26481481, -0.2462963, -0.26481481,
2814	-0.2462963, -0.26481481, -0.2462963])
2815
2816
2817	#print domain.quantities['stage'].extrapolate_second_order()
2818	#domain.distribute_to_vertices_and_edges()
2819	#print domain.quantities['stage'].vertex_values[:,0]
2820
2821	#Evolution
2822	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
2823	pass
2824
2825
2826	assert allclose(domain.quantities['stage'].centroid_values,
2827	[-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
2828	-0.02901975, -0.0141361, -0.02898816, -0.01418135, -0.02961409, -0.01403487,
2829	-0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
2830	-0.07571145, -0.06235231, -0.0756817, -0.06245309, -0.07652292, -0.06289946,
2831	-0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934, -0.1107174,
2832	-0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
2833	-0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493, -0.15697919,
2834	-0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
2835	-0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
2836	-0.1833593, -0.19505871, -0.1818806, -0.19418042, -0.18586159, -0.19576946,
2837	-0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
2838	-0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
2839
2840	assert allclose(domain.quantities['xmomentum'].centroid_values,
2841	[ 0.00831121, 0.00317948, 0.00731797, 0.00334939, 0.00764717, 0.00348053,
2842	0.00788729, 0.00356522, 0.00780649, 0.00341919, 0.00693525, 0.00310375,
2843	0.02166196, 0.01421475, 0.02017737, 0.01316839, 0.02037015, 0.01368659,
2844	0.02106, 0.01399161, 0.02074514, 0.01354935, 0.01887407, 0.0123113,
2845	0.03775083, 0.02855197, 0.03689337, 0.02759782, 0.03732848, 0.02812072,
2846	0.03872545, 0.02913348, 0.03880939, 0.02803804, 0.03546499, 0.0260039,
2847	0.0632131, 0.04730634, 0.0576324, 0.04592336, 0.05790921, 0.04690514,
2848	0.05986467, 0.04871165, 0.06170068, 0.04811572, 0.05657041, 0.04416292,
2849	0.08489642, 0.07188097, 0.07835261, 0.06843406, 0.07986412, 0.0698247,
2850	0.08201071, 0.07216756, 0.08378418, 0.07273624, 0.080399, 0.06645841,
2851	0.01631548, 0.04691608, 0.0206632, 0.044409, 0.02115518, 0.04560305,
2852	0.02160608, 0.04663725, 0.02174734, 0.04795559, 0.02281427, 0.05667111])
2853
2854
2855	assert allclose(domain.quantities['ymomentum'].centroid_values,
2856	[ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
2857	-9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
2858	-1.99445058e-004, -3.28493467e-004, 6.68217349e-005, -8.42042805e-006,
2859	5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
2860	-8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
2861	-2.70160173e-004, -5.35485454e-004, 2.60780997e-004, 3.12145471e-005,
2862	5.16189608e-004, 1.07069062e-004, 9.29989252e-005, -3.71211119e-004,
2863	1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
2864	-4.74025708e-004, -6.94463009e-004, 6.15092843e-005, 2.22106820e-004,
2865	-6.29589294e-004, 2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
2866	-4.17914641e-004, 6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
2867	-1.67424308e-003, -7.39485066e-004, -1.59966988e-003, 5.68262838e-005,
2868	-1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
2869	-1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
2870	-2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
2871	-1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
2872	-6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
2873	-4.50964850e-003, -3.06319963e-003, 6.08950810e-004, -4.79537921e-004])
2874
2875	os.remove(domain.filename + '.sww')
2876
2877
2878	def test_temp_play(self):
2879
2880	from mesh_factory import rectangular
2881	from shallow_water import Domain, Reflective_boundary, Constant_height
2882	from Numeric import array
2883
2884	#Create basic mesh
2885	points, vertices, boundary = rectangular(5, 5)
2886
2887	#Create shallow water domain
2888	domain = Domain(points, vertices, boundary)
2889	domain.smooth = False
2890	domain.default_order=2
2891	domain.beta_h = 0.0 #Use first order in h-limiter
2892
2893	#Bed-slope and friction at vertices (and interpolated elsewhere)
2894	def x_slope(x, y):
2895	return -x/3
2896
2897	domain.set_quantity('elevation', x_slope)
2898
2899	# Boundary conditions
2900	Br = Reflective_boundary(domain)
2901	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
2902
2903	#Initial condition
2904	domain.set_quantity('stage', Constant_height(x_slope, 0.05))
2905	domain.check_integrity()
2906
2907	#Evolution
2908	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
2909	pass
2910
2911	assert allclose(domain.quantities['stage'].centroid_values[:4],
2912	[0.00206836, 0.01296714, 0.00363415, 0.01438924])
2913	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
2914	[0.01360154, 0.00671133, 0.01264578, 0.00648503])
2915	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
2916	[-1.19201077e-003, -7.23647546e-004,
2917	-6.39083123e-005, 6.29815168e-005])
2918
2919	os.remove(domain.filename + '.sww')
2920
2921	def test_complex_bed(self):
2922	#No friction is tested here
2923
2924	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2925	Transmissive_boundary, Time_boundary,\
2926	Weir_simple as Weir, Constant_height
2927
2928	from mesh_factory import rectangular
2929	from Numeric import array
2930
2931	N = 12
2932	points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
2933	origin=(-0.07, 0))
2934
2935
2936	domain = Domain(points, vertices, boundary)
2937	domain.smooth = False
2938	domain.visualise = False
2939	domain.default_order=2
2940
2941
2942	inflow_stage = 0.1
2943	Z = Weir(inflow_stage)
2944	domain.set_quantity('elevation', Z)
2945
2946	Br = Reflective_boundary(domain)
2947	Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
2948	domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
2949
2950	domain.set_quantity('stage', Constant_height(Z, 0.))
2951
2952	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
2953	pass
2954
2955
2956	#print domain.quantities['stage'].centroid_values
2957
2958	#FIXME: These numbers were from version before 25/10
2959	#assert allclose(domain.quantities['stage'].centroid_values,
2960	# [3.95822638e-002, 5.61022588e-002, 4.66437868e-002, 5.73081011e-002,
2961	# 4.72394613e-002, 5.74684939e-002, 4.74309483e-002, 5.77458084e-002,
2962	# 4.80628177e-002, 5.85656225e-002, 4.90498542e-002, 6.02609831e-002,
2963	# 1.18470315e-002, 1.75136443e-002, 1.18035266e-002, 2.15565695e-002,
2964	# 1.31620268e-002, 2.14351640e-002, 1.32351076e-002, 2.15450687e-002,
2965	# 1.36414028e-002, 2.24274619e-002, 1.51689511e-002, 2.21789655e-002,
2966	# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
2967	# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
2968	# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003, 2.73985295e-004,
2969	# 1.86166683e-001, 8.74070369e-002, 1.86166712e-001, 8.74035875e-002,
2970	# 6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
2971	# 6.11666725e-002, 8.73846774e-002, 1.86166697e-001, 8.74171550e-002,
2972	# -4.83333333e-002, 1.18333333e-001, -4.83333333e-002, 1.18333333e-001,
2973	# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
2974	# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002, 1.18333333e-001,
2975	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2976	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2977	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
2978	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2979	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2980	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
2981	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2982	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2983	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
2984	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2985	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2986	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
2987	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2988	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2989	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
2990	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2991	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2992	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
2993	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2994	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
2995	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
2996
2997	os.remove(domain.filename + '.sww')
2998
2999	def test_spatio_temporal_boundary_1(self):
3000	"""Test that boundary values can be read from file and interpolated
3001	in both time and space.
3002
3003	Verify that the same steady state solution is arrived at and that
3004	time interpolation works.
3005
3006	The full solution history is not exactly the same as
3007	file boundary must read and interpolate from smoothed version
3008	as stored in sww.
3009	"""
3010	import time
3011
3012	#Create sww file of simple propagation from left to right
3013	#through rectangular domain
3014
3015	from mesh_factory import rectangular
3016
3017	#Create basic mesh
3018	points, vertices, boundary = rectangular(3, 3)
3019
3020	#Create shallow water domain
3021	domain1 = Domain(points, vertices, boundary)
3022
3023	domain1.reduction = mean
3024	domain1.smooth = False #Exact result
3025
3026	domain1.default_order = 2
3027	domain1.store = True
3028	domain1.set_datadir('.')
3029	domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
3030
3031	#FIXME: This is extremely important!
3032	#How can we test if they weren't stored?
3033	domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
3034
3035
3036	#Bed-slope and friction at vertices (and interpolated elsewhere)
3037	domain1.set_quantity('elevation', 0)
3038	domain1.set_quantity('friction', 0)
3039
3040	# Boundary conditions
3041	Br = Reflective_boundary(domain1)
3042	Bd = Dirichlet_boundary([0.3,0,0])
3043	domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
3044	#Initial condition
3045	domain1.set_quantity('stage', 0)
3046	domain1.check_integrity()
3047
3048	finaltime = 5
3049	#Evolution (full domain - large steps)
3050	for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
3051	pass
3052	#domain1.write_time()
3053
3054	cv1 = domain1.quantities['stage'].centroid_values
3055
3056
3057	#Create a triangle shaped domain (reusing coordinates from domain 1),
3058	#formed from the lower and right hand boundaries and
3059	#the sw-ne diagonal
3060	#from domain 1. Call it domain2
3061
3062	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3063	[2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3064	[1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3065
3066	vertices = [ [1,2,0], [3,4,1], [2,1,4], [4,5,2],
3067	[6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
3068
3069	boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
3070	(4,2):'right', (6,2):'right', (8,2):'right',
3071	(0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
3072
3073	domain2 = Domain(points, vertices, boundary)
3074
3075	domain2.reduction = domain1.reduction
3076	domain2.smooth = False
3077	domain2.default_order = 2
3078
3079	#Bed-slope and friction at vertices (and interpolated elsewhere)
3080	domain2.set_quantity('elevation', 0)
3081	domain2.set_quantity('friction', 0)
3082	domain2.set_quantity('stage', 0)
3083
3084	# Boundary conditions
3085	Br = Reflective_boundary(domain2)
3086	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
3087	domain2)
3088	domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
3089	domain2.check_integrity()
3090
3091
3092
3093	#Evolution (small steps)
3094	for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
3095	pass
3096
3097
3098	#Use output from domain1 as spatio-temporal boundary for domain2
3099	#and verify that results at right hand side are close.
3100
3101	cv2 = domain2.quantities['stage'].centroid_values
3102
3103	#print take(cv1, (12,14,16)) #Right
3104	#print take(cv2, (4,6,8))
3105	#print take(cv1, (0,6,12)) #Bottom
3106	#print take(cv2, (0,1,4))
3107	#print take(cv1, (0,8,16)) #Diag
3108	#print take(cv2, (0,3,8))
3109
3110	assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
3111	assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
3112	assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
3113
3114	#Cleanup
3115	os.remove(domain1.filename + '.' + domain1.format)
3116	os.remove(domain2.filename + '.' + domain2.format)
3117
3118
3119
3120	def test_spatio_temporal_boundary_2(self):
3121	"""Test that boundary values can be read from file and interpolated
3122	in both time and space.
3123	This is a more basic test, verifying that boundary object
3124	produces the expected results
3125
3126
3127	"""
3128	import time
3129
3130	#Create sww file of simple propagation from left to right
3131	#through rectangular domain
3132
3133	from mesh_factory import rectangular
3134
3135	#Create basic mesh
3136	points, vertices, boundary = rectangular(3, 3)
3137
3138	#Create shallow water domain
3139	domain1 = Domain(points, vertices, boundary)
3140
3141	domain1.reduction = mean
3142	domain1.smooth = True #To mimic MOST output
3143
3144	domain1.default_order = 2
3145	domain1.store = True
3146	domain1.set_datadir('.')
3147	domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
3148
3149	#FIXME: This is extremely important!
3150	#How can we test if they weren't stored?
3151	domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
3152
3153
3154	#Bed-slope and friction at vertices (and interpolated elsewhere)
3155	domain1.set_quantity('elevation', 0)
3156	domain1.set_quantity('friction', 0)
3157
3158	# Boundary conditions
3159	Br = Reflective_boundary(domain1)
3160	Bd = Dirichlet_boundary([0.3,0,0])
3161	domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
3162	#Initial condition
3163	domain1.set_quantity('stage', 0)
3164	domain1.check_integrity()
3165
3166	finaltime = 5
3167	#Evolution (full domain - large steps)
3168	for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
3169	pass
3170	#domain1.write_time()
3171
3172
3173	#Create an triangle shaped domain (coinciding with some
3174	#coordinates from domain 1),
3175	#formed from the lower and right hand boundaries and
3176	#the sw-ne diagonal
3177	#from domain 1. Call it domain2
3178
3179	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
3180	[2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
3181	[1,0], [1,1.0/3], [1,2.0/3], [1,1]]
3182
3183	vertices = [ [1,2,0],
3184	[3,4,1], [2,1,4], [4,5,2],
3185	[6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
3186
3187	boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
3188	(4,2):'right', (6,2):'right', (8,2):'right',
3189	(0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
3190
3191	domain2 = Domain(points, vertices, boundary)
3192
3193	domain2.reduction = domain1.reduction
3194	domain2.smooth = False
3195	domain2.default_order = 2
3196
3197	#Bed-slope and friction at vertices (and interpolated elsewhere)
3198	domain2.set_quantity('elevation', 0)
3199	domain2.set_quantity('friction', 0)
3200	domain2.set_quantity('stage', 0)
3201
3202
3203	#Read results for specific timesteps t=1 and t=2
3204	from Scientific.IO.NetCDF import NetCDFFile
3205	fid = NetCDFFile(domain1.filename + '.' + domain1.format)
3206
3207	x = fid.variables['x'][:]
3208	y = fid.variables['y'][:]
3209	s1 = fid.variables['stage'][1,:]
3210	s2 = fid.variables['stage'][2,:]
3211	fid.close()
3212
3213	from Numeric import take, reshape, concatenate
3214	shp = (len(x), 1)
3215	points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
3216	#The diagonal points of domain 1 are 0, 5, 10, 15
3217
3218	#print points[0], points[5], points[10], points[15]
3219	assert allclose( take(points, [0,5,10,15]),
3220	[[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
3221
3222
3223	# Boundary conditions
3224	Br = Reflective_boundary(domain2)
3225	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
3226	domain2)
3227	domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
3228	domain2.check_integrity()
3229
3230	#Test that interpolation points are the mid points of the all boundary
3231	#segments
3232
3233	boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
3234	[1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
3235	[1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
3236
3237	boundary_midpoints.sort()
3238	R = Bf.F.interpolation_points.tolist()
3239	R.sort()
3240	assert allclose(boundary_midpoints, R)
3241
3242	#Check spatially interpolated output at time == 1
3243	domain2.time = 1
3244
3245	#First diagonal midpoint
3246	R0 = Bf.evaluate(0,0)
3247	assert allclose(R0[0], (s1[0] + s1[5])/2)
3248
3249	#Second diagonal midpoint
3250	R0 = Bf.evaluate(3,0)
3251	assert allclose(R0[0], (s1[5] + s1[10])/2)
3252
3253	#First diagonal midpoint
3254	R0 = Bf.evaluate(8,0)
3255	assert allclose(R0[0], (s1[10] + s1[15])/2)
3256
3257	#Check spatially interpolated output at time == 2
3258	domain2.time = 2
3259
3260	#First diagonal midpoint
3261	R0 = Bf.evaluate(0,0)
3262	assert allclose(R0[0], (s2[0] + s2[5])/2)
3263
3264	#Second diagonal midpoint
3265	R0 = Bf.evaluate(3,0)
3266	assert allclose(R0[0], (s2[5] + s2[10])/2)
3267
3268	#First diagonal midpoint
3269	R0 = Bf.evaluate(8,0)
3270	assert allclose(R0[0], (s2[10] + s2[15])/2)
3271
3272
3273	#Now check temporal interpolation
3274
3275	domain2.time = 1 + 2.0/3
3276
3277	#First diagonal midpoint
3278	R0 = Bf.evaluate(0,0)
3279	assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3)
3280
3281	#Second diagonal midpoint
3282	R0 = Bf.evaluate(3,0)
3283	assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3)
3284
3285	#First diagonal midpoint
3286	R0 = Bf.evaluate(8,0)
3287	assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
3288
3289
3290
3291	#Cleanup
3292	os.remove(domain1.filename + '.' + domain1.format)
3293
3294
3295	def test_pmesh2Domain(self):
3296	import os
3297	import tempfile
3298
3299	fileName = tempfile.mktemp(".tsh")
3300	file = open(fileName,"w")
3301	file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
3302	0 0.0 0.0 0.0 0.0 0.01 \n \
3303	1 1.0 0.0 10.0 10.0 0.02 \n \
3304	2 0.0 1.0 0.0 10.0 0.03 \n \
3305	3 0.5 0.25 8.0 12.0 0.04 \n \
3306	# Vert att title \n \
3307	elevation \n \
3308	stage \n \
3309	friction \n \
3310	2 # <triangle #> [<vertex #>] [<neigbouring triangle #>] \n\
3311	0 0 3 2 -1 -1 1 dsg\n\
3312	1 0 1 3 -1 0 -1 ole nielsen\n\
3313	4 # <segment #> <vertex #> <vertex #> [boundary tag] \n\
3314	0 1 0 2 \n\
3315	1 0 2 3 \n\
3316	2 2 3 \n\
3317	3 3 1 1 \n\
3318	3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
3319	0 216.0 -86.0 \n \
3320	1 160.0 -167.0 \n \
3321	2 114.0 -91.0 \n \
3322	3 # <vertex #> <vertex #> [boundary tag] ...Mesh Segments... \n \
3323	0 0 1 0 \n \
3324	1 1 2 0 \n \
3325	2 2 0 0 \n \
3326	0 # <x> <y> ...Mesh Holes... \n \
3327	0 # <x> <y> <attribute>...Mesh Regions... \n \
3328	0 # <x> <y> <attribute>...Mesh Regions, area... \n\
3329	#Geo reference \n \
3330	56 \n \
3331	140 \n \
3332	120 \n")
3333	file.close()
3334
3335	tags = {}
3336	b1 = Dirichlet_boundary(conserved_quantities = array([0.0]))
3337	b2 = Dirichlet_boundary(conserved_quantities = array([1.0]))
3338	b3 = Dirichlet_boundary(conserved_quantities = array([2.0]))
3339	tags["1"] = b1
3340	tags["2"] = b2
3341	tags["3"] = b3
3342
3343	#from pmesh2domain import pmesh_to_domain_instance
3344	#domain = pmesh_to_domain_instance(fileName, Domain)
3345
3346	domain = Domain(mesh_filename=fileName)
3347	#verbose=True, use_cache=True)
3348
3349	#print "domain.tagged_elements", domain.tagged_elements
3350	## check the quantities
3351	#print domain.quantities['elevation'].vertex_values
3352	answer = [[0., 8., 0.],
3353	[0., 10., 8.]]
3354	assert allclose(domain.quantities['elevation'].vertex_values,
3355	answer)
3356
3357	#print domain.quantities['stage'].vertex_values
3358	answer = [[0., 12., 10.],
3359	[0., 10., 12.]]
3360	assert allclose(domain.quantities['stage'].vertex_values,
3361	answer)
3362
3363	#print domain.quantities['friction'].vertex_values
3364	answer = [[0.01, 0.04, 0.03],
3365	[0.01, 0.02, 0.04]]
3366	assert allclose(domain.quantities['friction'].vertex_values,
3367	answer)
3368
3369	#print domain.quantities['friction'].vertex_values
3370	assert allclose(domain.tagged_elements['dsg'][0],0)
3371	assert allclose(domain.tagged_elements['ole nielsen'][0],1)
3372
3373	self.failUnless( domain.boundary[(1, 0)] == '1',
3374	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3375	self.failUnless( domain.boundary[(1, 2)] == '2',
3376	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3377	self.failUnless( domain.boundary[(0, 1)] == '3',
3378	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3379	self.failUnless( domain.boundary[(0, 0)] == 'exterior',
3380	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3381	#print "domain.boundary",domain.boundary
3382	self.failUnless( len(domain.boundary) == 4,
3383	"test_pmesh2Domain Too many boundaries")
3384	#FIXME change to use get_xllcorner
3385	#print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
3386	self.failUnless(domain.geo_reference.xllcorner == 140.0,
3387	"bad geo_referece")
3388	#************
3389
3390
3391	domain = Domain(fileName)
3392
3393	#print "domain.tagged_elements", domain.tagged_elements
3394	## check the quantities
3395	#print domain.quantities['elevation'].vertex_values
3396	answer = [[0., 8., 0.],
3397	[0., 10., 8.]]
3398	assert allclose(domain.quantities['elevation'].vertex_values,
3399	answer)
3400
3401	#print domain.quantities['stage'].vertex_values
3402	answer = [[0., 12., 10.],
3403	[0., 10., 12.]]
3404	assert allclose(domain.quantities['stage'].vertex_values,
3405	answer)
3406
3407	#print domain.quantities['friction'].vertex_values
3408	answer = [[0.01, 0.04, 0.03],
3409	[0.01, 0.02, 0.04]]
3410	assert allclose(domain.quantities['friction'].vertex_values,
3411	answer)
3412
3413	#print domain.quantities['friction'].vertex_values
3414	assert allclose(domain.tagged_elements['dsg'][0],0)
3415	assert allclose(domain.tagged_elements['ole nielsen'][0],1)
3416
3417	self.failUnless( domain.boundary[(1, 0)] == '1',
3418	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3419	self.failUnless( domain.boundary[(1, 2)] == '2',
3420	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3421	self.failUnless( domain.boundary[(0, 1)] == '3',
3422	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3423	self.failUnless( domain.boundary[(0, 0)] == 'exterior',
3424	"test_tags_to_boundaries failed. Single boundary wasn't added.")
3425	#print "domain.boundary",domain.boundary
3426	self.failUnless( len(domain.boundary) == 4,
3427	"test_pmesh2Domain Too many boundaries")
3428	#FIXME change to use get_xllcorner
3429	#print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
3430	self.failUnless(domain.geo_reference.xllcorner == 140.0,
3431	"bad geo_referece")
3432	#************
3433	os.remove(fileName)
3434
3435	#-------------------------------------------------------------
3436	if __name__ == "__main__":
3437	suite = unittest.makeSuite(Test_Shallow_Water,'test')
3438	runner = unittest.TextTestRunner(verbosity=1)
3439	runner.run(suite)

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

Context Navigation

source: inundation/pyvolution/test_shallow_water.py @ 3383

Download in other formats: