Context Navigation

source: branches/source_numpy_conversion/anuga/shallow_water/test_shallow_water_domain.py @ 7177

Last change on this file since 7177 was 5947, checked in by rwilson, 16 years ago
More NumPy? changes.
File size: 227.5 KB

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

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

Download in other formats: