Context Navigation

test_shallow_water.py @ 7077

Last change on this file since 7077 was 3514, checked in by duncan, 19 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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

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

Context Navigation

source: branches/inundation-numpy-branch/pyvolution/test_shallow_water.py @ 7077

Download in other formats: