Context Navigation

source: anuga_core/source/anuga/shallow_water_balanced/test_swb_basic.py @ 7575

Last change on this file since 7575 was 7575, checked in by steve, 14 years ago
Version of compute fluxes with discontinuous bed
File size: 93.3 KB

Line
1	#!/usr/bin/env python
2
3	import unittest, os
4	import os.path
5	from math import pi, sqrt
6	import tempfile
7
8	from anuga.config import g, epsilon
9	from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
10	from anuga.utilities.numerical_tools import mean
11	from anuga.utilities.polygon import is_inside_polygon
12	from anuga.coordinate_transforms.geo_reference import Geo_reference
13	from anuga.abstract_2d_finite_volumes.quantity import Quantity
14	from anuga.geospatial_data.geospatial_data import Geospatial_data
15	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
16
17	from anuga.utilities.system_tools import get_pathname_from_package
18	from swb_domain import *
19
20	import numpy as num
21
22	# Get gateway to C implementation of flux function for direct testing
23	from shallow_water_ext import flux_function_central as flux_function
24
25
26	# For test_fitting_using_shallow_water_domain example
27	def linear_function(point):
28	point = num.array(point)
29	return point[:,0]+point[:,1]
30
31	class Weir:
32	"""Set a bathymetry for weir with a hole and a downstream gutter
33	x,y are assumed to be in the unit square
34	"""
35
36	def __init__(self, stage):
37	self.inflow_stage = stage
38
39	def __call__(self, x, y):
40	N = len(x)
41	assert N == len(y)
42
43	z = num.zeros(N, num.float)
44	for i in range(N):
45	z[i] = -x[i]/2 #General slope
46
47	#Flattish bit to the left
48	if x[i] < 0.3:
49	z[i] = -x[i]/10
50
51	#Weir
52	if x[i] >= 0.3 and x[i] < 0.4:
53	z[i] = -x[i]+0.9
54
55	#Dip
56	x0 = 0.6
57	depth = -1.0
58	plateaux = -0.6
59	if y[i] < 0.7:
60	if x[i] > x0 and x[i] < 0.9:
61	z[i] = depth
62	#RHS plateaux
63	if x[i] >= 0.9:
64	z[i] = plateaux
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	elif x[i] >= 0.8:
70	#RHS plateaux
71	z[i] = plateaux
72	elif y[i] >= 1.5:
73	if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
74	#Widen up and stay at constant depth
75	z[i] = depth-1.5/5
76	elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
77	#RHS plateaux
78	z[i] = plateaux
79
80	#Hole in weir (slightly higher than inflow condition)
81	if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
82	z[i] = -x[i]+self.inflow_stage + 0.02
83
84	#Channel behind weir
85	x0 = 0.5
86	if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
87	z[i] = -x[i]+self.inflow_stage + 0.02
88
89	if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
90	#Flatten it out towards the end
91	z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
92
93	# Hole to the east
94	x0 = 1.1
95	y0 = 0.35
96	if num.sqrt((2(x[i]-x0))2 + (2(y[i]-y0))**2) < 0.2:
97	z[i] = num.sqrt(((x[i]-x0))2 + ((y[i]-y0))2)-1.0
98
99	#Tiny channel draining hole
100	if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
101	z[i] = -0.9 #North south
102
103	if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
104	z[i] = -1.0 + (x[i]-0.9)/3 #East west
105
106	# Stuff not in use
107
108	# Upward slope at inlet to the north west
109	# if x[i] < 0.0: # and y[i] > 0.5:
110	# #z[i] = -y[i]+0.5 #-x[i]/2
111	# z[i] = x[i]/4 - y[i]**2 + 0.5
112
113	# Hole to the west
114	# x0 = -0.4; y0 = 0.35 # center
115	# if sqrt((2(x[i]-x0))2 + (2(y[i]-y0))**2) < 0.2:
116	# z[i] = sqrt(((x[i]-x0))2 + ((y[i]-y0))2)-0.2
117
118	return z/2
119
120
121
122	#########################################################
123
124	class Test_swb_basic(unittest.TestCase):
125	def setUp(self):
126	pass
127
128	def tearDown(self):
129	pass
130
131
132	def test_rotate(self):
133	normal = num.array([0.0, -1.0])
134
135	q = num.array([1.0, 2.0, 3.0])
136
137	r = rotate(q, normal, direction = 1)
138	assert r[0] == 1
139	assert r[1] == -3
140	assert r[2] == 2
141
142	w = rotate(r, normal, direction = -1)
143	assert num.allclose(w, q)
144
145	# Check error check
146	try:
147	rotate(r, num.array([1, 1, 1]))
148	except:
149	pass
150	else:
151	raise Exception, 'Should have raised an exception'
152
153	# Individual flux tests
154	def test_flux_zero_case(self):
155	ql = num.zeros(3, num.float)
156	qr = num.zeros(3, num.float)
157	normal = num.zeros(2, num.float)
158	edgeflux = num.zeros(3, num.float)
159	zl = zr = 0.
160	H0 = 1.0e-3 # As suggested in the manual
161
162	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
163
164	assert num.allclose(edgeflux, [0,0,0])
165	assert max_speed == 0.
166
167	def test_flux_constants(self):
168	w = 2.0
169
170	normal = num.array([1.,0])
171	ql = num.array([w, 0, 0])
172	qr = num.array([w, 0, 0])
173	edgeflux = num.zeros(3, num.float)
174	zl = zr = 0.
175	h = w - (zl+zr)/2
176	H0 = 0.0
177
178	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
179	assert num.allclose(edgeflux, [0., 0.5gh**2, 0.])
180	assert max_speed == num.sqrt(g*h)
181
182	#def test_flux_slope(self):
183	# #FIXME: TODO
184	# w = 2.0
185	#
186	# normal = array([1.,0])
187	# ql = array([w, 0, 0])
188	# qr = array([w, 0, 0])
189	# zl = zr = 0.
190	# h = w - (zl+zr)/2
191	#
192	# flux, max_speed = flux_function(normal, ql, qr, zl, zr)
193	#
194	# assert allclose(flux, [0., 0.5gh**2, 0.])
195	# assert max_speed == sqrt(g*h)
196
197	def test_flux1(self):
198	# Use data from previous version of abstract_2d_finite_volumes
199	normal = num.array([1., 0])
200	ql = num.array([-0.2, 2, 3])
201	qr = num.array([-0.2, 2, 3])
202	zl = zr = -0.5
203	edgeflux = num.zeros(3, num.float)
204
205	H0 = 0.0
206
207	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
208	epsilon, g, H0)
209
210	assert num.allclose(edgeflux, [2., 13.77433333, 20.])
211	assert num.allclose(max_speed, 8.38130948661)
212
213	def test_flux2(self):
214	# Use data from previous version of abstract_2d_finite_volumes
215	normal = num.array([0., -1.])
216	ql = num.array([-0.075, 2, 3])
217	qr = num.array([-0.075, 2, 3])
218	zl = zr = -0.375
219
220	edgeflux = num.zeros(3, num.float)
221	H0 = 0.0
222	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
223	epsilon, g, H0)
224
225	assert num.allclose(edgeflux, [-3., -20.0, -30.441])
226	assert num.allclose(max_speed, 11.7146428199)
227
228	def test_flux3(self):
229	# Use data from previous version of abstract_2d_finite_volumes
230	normal = num.array([-sqrt(2)/2, sqrt(2)/2])
231	ql = num.array([-0.075, 2, 3])
232	qr = num.array([-0.075, 2, 3])
233	zl = zr = -0.375
234
235	edgeflux = num.zeros(3, num.float)
236	H0 = 0.0
237	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
238	epsilon, g, H0)
239
240	assert num.allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019])
241	assert num.allclose(max_speed, 4.0716654239)
242
243	def test_flux4(self):
244	# Use data from previous version of abstract_2d_finite_volumes
245	normal = num.array([-sqrt(2)/2, sqrt(2)/2])
246	ql = num.array([-0.34319278, 0.10254161, 0.07273855])
247	qr = num.array([-0.30683287, 0.1071986, 0.05930515])
248	zl = zr = -0.375
249
250	edgeflux = num.zeros(3, num.float)
251	H0 = 0.0
252	max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
253	epsilon, g, H0)
254
255	assert num.allclose(edgeflux, [-0.04072676, -0.07096636, -0.01604364])
256	assert num.allclose(max_speed, 1.31414103233)
257
258	def test_flux_computation(self):
259	"""test flux calculation (actual C implementation)
260
261	This one tests the constant case where only the pressure term
262	contributes to each edge and cancels out once the total flux has
263	been summed up.
264	"""
265
266	a = [0.0, 0.0]
267	b = [0.0, 2.0]
268	c = [2.0, 0.0]
269	d = [0.0, 4.0]
270	e = [2.0, 2.0]
271	f = [4.0, 0.0]
272
273	points = [a, b, c, d, e, f]
274	# bac, bce, ecf, dbe
275	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
276
277	domain = Domain(points, vertices)
278	domain.check_integrity()
279
280	# The constant case
281	domain.set_quantity('elevation', -1)
282	domain.set_quantity('stage', 1)
283
284	domain.compute_fluxes()
285	# Central triangle
286	assert num.allclose(domain.get_quantity('stage').explicit_update[1], 0)
287
288	# The more general case
289	def surface(x, y):
290	return -x/2
291
292	domain.set_quantity('elevation', -10)
293	domain.set_quantity('stage', surface)
294	domain.set_quantity('xmomentum', 1)
295
296	domain.compute_fluxes()
297
298	#print domain.get_quantity('stage').explicit_update
299	# FIXME (Ole): TODO the general case
300	#assert allclose(domain.get_quantity('stage').explicit_update[1], ...??)
301
302	def test_sw_domain_simple(self):
303	a = [0.0, 0.0]
304	b = [0.0, 2.0]
305	c = [2.0, 0.0]
306	d = [0.0, 4.0]
307	e = [2.0, 2.0]
308	f = [4.0, 0.0]
309
310	points = [a, b, c, d, e, f]
311	# bac, bce, ecf, dbe
312	vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
313
314	#from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain
315	#msg = 'The class %s is not a subclass of the generic domain class %s'\
316	# %(DomainClass, Domain)
317	#assert issubclass(DomainClass, Domain), msg
318
319	domain = Domain(points, vertices)
320	domain.check_integrity()
321
322	for name in ['stage', 'xmomentum', 'ymomentum',
323	'elevation', 'friction']:
324	assert domain.quantities.has_key(name)
325
326	assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
327
328	def xtest_catching_negative_heights(self):
329	#OBSOLETE
330
331	a = [0.0, 0.0]
332	b = [0.0, 2.0]
333	c = [2.0, 0.0]
334	d = [0.0, 4.0]
335	e = [2.0, 2.0]
336	f = [4.0, 0.0]
337
338	points = [a, b, c, d, e, f]
339	# bac, bce, ecf, dbe
340	vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
341
342	domain = Domain(points, vertices)
343	val0 = 2. + 2.0/3
344	val1 = 4. + 4.0/3
345	val2 = 8. + 2.0/3
346	val3 = 2. + 8.0/3
347
348	zl = zr = 4 # Too large
349	domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
350	domain.set_quantity('stage', [[val0, val0-1, val0-2],
351	[val1, val1+1, val1],
352	[val2, val2-2, val2],
353	[val3-0.5, val3, val3]])
354
355	#Should fail
356	try:
357	domain.check_integrity()
358	except:
359	pass
360
361	def test_get_wet_elements(self):
362	a = [0.0, 0.0]
363	b = [0.0, 2.0]
364	c = [2.0, 0.0]
365	d = [0.0, 4.0]
366	e = [2.0, 2.0]
367	f = [4.0, 0.0]
368
369	points = [a, b, c, d, e, f]
370	# bac, bce, ecf, dbe
371	vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
372
373	domain = Domain(points, vertices)
374
375	val0 = 2. + 2.0/3
376	val1 = 4. + 4.0/3
377	val2 = 8. + 2.0/3
378	val3 = 2. + 8.0/3
379
380	zl = zr = 5
381	domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
382	domain.set_quantity('stage', [[val0, val0-1, val0-2],
383	[val1, val1+1, val1],
384	[val2, val2-2, val2],
385	[val3-0.5, val3, val3]])
386
387	domain.check_integrity()
388
389	indices = domain.get_wet_elements()
390	assert num.allclose(indices, [1, 2])
391
392	indices = domain.get_wet_elements(indices=[0, 1, 3])
393	assert num.allclose(indices, [1])
394
395	def test_get_maximum_inundation_1(self):
396	a = [0.0, 0.0]
397	b = [0.0, 2.0]
398	c = [2.0, 0.0]
399	d = [0.0, 4.0]
400	e = [2.0, 2.0]
401	f = [4.0, 0.0]
402
403	points = [a, b, c, d, e, f]
404	# bac, bce, ecf, dbe
405	vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
406
407	domain = Domain(points, vertices)
408
409	domain.set_quantity('elevation', lambda x, y: x+2*y) # 2 4 4 6
410	domain.set_quantity('stage', 3)
411
412	domain.check_integrity()
413
414	indices = domain.get_wet_elements()
415	assert num.allclose(indices, [0])
416
417	q = domain.get_maximum_inundation_elevation()
418	assert num.allclose(q, domain.get_quantity('elevation').\
419	get_values(location='centroids')[0])
420
421	x, y = domain.get_maximum_inundation_location()
422	assert num.allclose([x, y], domain.get_centroid_coordinates()[0])
423
424	def test_get_maximum_inundation_2(self):
425	"""test_get_maximum_inundation_2(self)
426
427	Test multiple wet cells with same elevation
428	"""
429
430	a = [0.0, 0.0]
431	b = [0.0, 2.0]
432	c = [2.0, 0.0]
433	d = [0.0, 4.0]
434	e = [2.0, 2.0]
435	f = [4.0, 0.0]
436
437	points = [a, b, c, d, e, f]
438	# bac, bce, ecf, dbe
439	vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
440
441	domain = Domain(points, vertices)
442
443	domain.set_quantity('elevation', lambda x, y: x+2*y) # 2 4 4 6
444	domain.set_quantity('stage', 4.1)
445
446	domain.check_integrity()
447
448	indices = domain.get_wet_elements()
449	assert num.allclose(indices, [0, 1, 2])
450
451	q = domain.get_maximum_inundation_elevation()
452	assert num.allclose(q, 4)
453
454	x, y = domain.get_maximum_inundation_location()
455	assert num.allclose([x, y], domain.get_centroid_coordinates()[1])
456
457	def test_get_maximum_inundation_3(self):
458	"""test_get_maximum_inundation_3(self)
459
460	Test of real runup example:
461	"""
462
463	from anuga.abstract_2d_finite_volumes.mesh_factory \
464	import rectangular_cross
465
466	initial_runup_height = -0.4
467	final_runup_height = -0.3
468
469	#--------------------------------------------------------------
470	# Setup computational domain
471	#--------------------------------------------------------------
472	N = 5
473	points, vertices, boundary = rectangular_cross(N, N)
474	domain = Domain(points, vertices, boundary)
475	domain.set_maximum_allowed_speed(1.0)
476
477	#--------------------------------------------------------------
478	# Setup initial conditions
479	#--------------------------------------------------------------
480	def topography(x, y):
481	return -x/2 # linear bed slope
482
483	# Use function for elevation
484	domain.set_quantity('elevation', topography)
485	domain.set_quantity('friction', 0.) # Zero friction
486	# Constant negative initial stage
487	domain.set_quantity('stage', initial_runup_height)
488
489	#--------------------------------------------------------------
490	# Setup boundary conditions
491	#--------------------------------------------------------------
492	Br = Reflective_boundary(domain) # Reflective wall
493	Bd = Dirichlet_boundary([final_runup_height, 0, 0]) # Constant inflow
494
495	# All reflective to begin with (still water)
496	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
497
498	#--------------------------------------------------------------
499	# Test initial inundation height
500	#--------------------------------------------------------------
501
502	indices = domain.get_wet_elements()
503	z = domain.get_quantity('elevation').get_values(location='centroids',
504	indices=indices)
505	assert num.alltrue(z < initial_runup_height)
506
507	q = domain.get_maximum_inundation_elevation()
508	# First order accuracy
509	assert num.allclose(q, initial_runup_height, rtol=1.0/N)
510
511	x, y = domain.get_maximum_inundation_location()
512
513	qref = domain.get_quantity('elevation').\
514	get_values(interpolation_points=[[x, y]])
515	assert num.allclose(q, qref)
516
517	wet_elements = domain.get_wet_elements()
518	wet_elevations = domain.get_quantity('elevation').\
519	get_values(location='centroids',
520	indices=wet_elements)
521	assert num.alltrue(wet_elevations < initial_runup_height)
522	assert num.allclose(wet_elevations, z)
523
524	#--------------------------------------------------------------
525	# Let triangles adjust
526	#--------------------------------------------------------------
527	for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
528	pass
529
530	#--------------------------------------------------------------
531	# Test inundation height again
532	#--------------------------------------------------------------
533	indices = domain.get_wet_elements()
534	z = domain.get_quantity('elevation').get_values(location='centroids',
535	indices=indices)
536
537	assert num.alltrue(z < initial_runup_height)
538
539	q = domain.get_maximum_inundation_elevation()
540	# First order accuracy
541	assert num.allclose(q, initial_runup_height, rtol=1.0/N)
542
543	x, y = domain.get_maximum_inundation_location()
544	qref = domain.get_quantity('elevation').\
545	get_values(interpolation_points=[[x, y]])
546	assert num.allclose(q, qref)
547
548	#--------------------------------------------------------------
549	# Update boundary to allow inflow
550	#--------------------------------------------------------------
551	domain.set_boundary({'right': Bd})
552
553	#--------------------------------------------------------------
554	# Evolve system through time
555	#--------------------------------------------------------------
556	for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
557	pass
558
559	#--------------------------------------------------------------
560	# Test inundation height again
561	#--------------------------------------------------------------
562	indices = domain.get_wet_elements()
563	z = domain.get_quantity('elevation').\
564	get_values(location='centroids', indices=indices)
565
566	assert num.alltrue(z < final_runup_height)
567
568	q = domain.get_maximum_inundation_elevation()
569	# First order accuracy
570	assert num.allclose(q, final_runup_height, rtol=1.0/N)
571
572	x, y = domain.get_maximum_inundation_location()
573	qref = domain.get_quantity('elevation').\
574	get_values(interpolation_points=[[x, y]])
575	assert num.allclose(q, qref)
576
577	wet_elements = domain.get_wet_elements()
578	wet_elevations = domain.get_quantity('elevation').\
579	get_values(location='centroids',
580	indices=wet_elements)
581	assert num.alltrue(wet_elevations < final_runup_height)
582	assert num.allclose(wet_elevations, z)
583
584	def test_get_maximum_inundation_from_sww(self):
585	"""test_get_maximum_inundation_from_sww(self)
586
587	Test of get_maximum_inundation_elevation()
588	and get_maximum_inundation_location() from data_manager.py
589
590	This is based on test_get_maximum_inundation_3(self) but works with the
591	stored results instead of with the internal data structure.
592
593	This test uses the underlying get_maximum_inundation_data for tests
594	"""
595
596	from anuga.abstract_2d_finite_volumes.mesh_factory \
597	import rectangular_cross
598	from data_manager import get_maximum_inundation_elevation
599	from data_manager import get_maximum_inundation_location
600	from data_manager import get_maximum_inundation_data
601
602	initial_runup_height = -0.4
603	final_runup_height = -0.3
604
605	#--------------------------------------------------------------
606	# Setup computational domain
607	#--------------------------------------------------------------
608	N = 10
609	points, vertices, boundary = rectangular_cross(N, N)
610	domain = Domain(points, vertices, boundary)
611	domain.set_name('runup_test')
612	#domain.set_maximum_allowed_speed(1.0)
613
614	# FIXME: This works better with old limiters so far
615	#domain.tight_slope_limiters = 0
616
617	#--------------------------------------------------------------
618	# Setup initial conditions
619	#--------------------------------------------------------------
620	def topography(x, y):
621	return -x/2 # linear bed slope
622
623	# Use function for elevation
624	domain.set_quantity('elevation', topography)
625	domain.set_quantity('friction', 0.) # Zero friction
626	# Constant negative initial stage
627	domain.set_quantity('stage', initial_runup_height)
628
629	#--------------------------------------------------------------
630	# Setup boundary conditions
631	#--------------------------------------------------------------
632	Br = Reflective_boundary(domain) # Reflective wall
633	Bd = Dirichlet_boundary([final_runup_height, 0, 0]) # Constant inflow
634
635	# All reflective to begin with (still water)
636	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
637
638	#--------------------------------------------------------------
639	# Test initial inundation height
640	#--------------------------------------------------------------
641	indices = domain.get_wet_elements()
642	z = domain.get_quantity('elevation').\
643	get_values(location='centroids', indices=indices)
644	assert num.alltrue(z < initial_runup_height)
645
646	q_ref = domain.get_maximum_inundation_elevation()
647	# First order accuracy
648	assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
649
650	#--------------------------------------------------------------
651	# Let triangles adjust
652	#--------------------------------------------------------------
653	for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
654	pass
655
656	#--------------------------------------------------------------
657	# Test inundation height again
658	#--------------------------------------------------------------
659	q_ref = domain.get_maximum_inundation_elevation()
660	q = get_maximum_inundation_elevation('runup_test.sww')
661	msg = 'We got %f, should have been %f' % (q, q_ref)
662	assert num.allclose(q, q_ref, rtol=1.0/N), msg
663
664	q = get_maximum_inundation_elevation('runup_test.sww')
665	msg = 'We got %f, should have been %f' % (q, initial_runup_height)
666	assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
667
668	# Test error condition if time interval is out
669	try:
670	q = get_maximum_inundation_elevation('runup_test.sww',
671	time_interval=[2.0, 3.0])
672	except ValueError:
673	pass
674	else:
675	msg = 'should have caught wrong time interval'
676	raise Exception, msg
677
678	# Check correct time interval
679	q, loc = get_maximum_inundation_data('runup_test.sww',
680	time_interval=[0.0, 3.0])
681	msg = 'We got %f, should have been %f' % (q, initial_runup_height)
682	assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
683	assert num.allclose(-loc[0]/2, q) # From topography formula
684
685	#--------------------------------------------------------------
686	# Update boundary to allow inflow
687	#--------------------------------------------------------------
688	domain.set_boundary({'right': Bd})
689
690	#--------------------------------------------------------------
691	# Evolve system through time
692	#--------------------------------------------------------------
693	q_max = None
694	for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
695	skip_initial_step = True):
696	q = domain.get_maximum_inundation_elevation()
697	if q > q_max:
698	q_max = q
699
700	#--------------------------------------------------------------
701	# Test inundation height again
702	#--------------------------------------------------------------
703	indices = domain.get_wet_elements()
704	z = domain.get_quantity('elevation').\
705	get_values(location='centroids', indices=indices)
706
707	assert num.alltrue(z < final_runup_height)
708
709	q = domain.get_maximum_inundation_elevation()
710	# First order accuracy
711	assert num.allclose(q, final_runup_height, rtol=1.0/N)
712
713	q, loc = get_maximum_inundation_data('runup_test.sww',
714	time_interval=[3.0, 3.0])
715	msg = 'We got %f, should have been %f' % (q, final_runup_height)
716	assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
717	assert num.allclose(-loc[0]/2, q) # From topography formula
718
719	q = get_maximum_inundation_elevation('runup_test.sww')
720	loc = get_maximum_inundation_location('runup_test.sww')
721	msg = 'We got %f, should have been %f' % (q, q_max)
722	assert num.allclose(q, q_max, rtol=1.0/N), msg
723	assert num.allclose(-loc[0]/2, q) # From topography formula
724
725	q = get_maximum_inundation_elevation('runup_test.sww',
726	time_interval=[0, 3])
727	msg = 'We got %f, should have been %f' % (q, q_max)
728	assert num.allclose(q, q_max, rtol=1.0/N), msg
729
730	# Check polygon mode
731	# Runup region
732	polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
733	q = get_maximum_inundation_elevation('runup_test.sww',
734	polygon = polygon,
735	time_interval=[0, 3])
736	msg = 'We got %f, should have been %f' % (q, q_max)
737	assert num.allclose(q, q_max, rtol=1.0/N), msg
738
739	# Offshore region
740	polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
741	q, loc = get_maximum_inundation_data('runup_test.sww',
742	polygon = polygon,
743	time_interval=[0, 3])
744	msg = 'We got %f, should have been %f' % (q, -0.475)
745	assert num.allclose(q, -0.475, rtol=1.0/N), msg
746	assert is_inside_polygon(loc, polygon)
747	assert num.allclose(-loc[0]/2, q) # From topography formula
748
749	# Dry region
750	polygon = [[0.0, 0.0], [0.2, 0.0], [0.2, 1.0], [0.0, 1.0]]
751	q, loc = get_maximum_inundation_data('runup_test.sww',
752	polygon = polygon,
753	time_interval=[0, 3])
754	msg = 'We got %s, should have been None' % (q)
755	assert q is None, msg
756	msg = 'We got %s, should have been None' % (loc)
757	assert loc is None, msg
758
759	# Check what happens if no time point is within interval
760	try:
761	q = get_maximum_inundation_elevation('runup_test.sww',
762	time_interval=[2.75, 2.75])
763	except AssertionError:
764	pass
765	else:
766	msg = 'Time interval should have raised an exception'
767	raise Exception, msg
768
769	# Cleanup
770	try:
771	os.remove(domain.get_name() + '.sww')
772	except:
773	pass
774	#FIXME(Ole): Windows won't allow removal of this
775
776
777
778
779
780	def test_another_runup_example(self):
781	"""test_another_runup_example
782
783	Test runup example where actual timeseries at interpolated
784	points are tested.
785	"""
786
787	from anuga.pmesh.mesh_interface import create_mesh_from_regions
788	from anuga.abstract_2d_finite_volumes.mesh_factory \
789	import rectangular_cross
790	from anuga.shallow_water import Domain
791	from anuga.shallow_water import Reflective_boundary
792	from anuga.shallow_water import Dirichlet_boundary
793
794	#-----------------------------------------------------------------
795	# Setup computational domain
796	#-----------------------------------------------------------------
797	points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
798	domain = Domain(points, vertices, boundary) # Create domain
799	domain.set_default_order(2)
800	domain.set_quantities_to_be_stored(None)
801	domain.H0 = 1.0e-3
802
803	#-----------------------------------------------------------------
804	# Setup initial conditions
805	#-----------------------------------------------------------------
806	def topography(x, y):
807	return -x/2 # linear bed slope
808
809	domain.set_quantity('elevation', topography)
810	domain.set_quantity('friction', 0.0)
811	domain.set_quantity('stage', expression='elevation')
812
813	#----------------------------------------------------------------
814	# Setup boundary conditions
815	#----------------------------------------------------------------
816	Br = Reflective_boundary(domain) # Solid reflective wall
817	Bd = Dirichlet_boundary([-0.2, 0., 0.]) # Constant boundary values
818	domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
819
820	#----------------------------------------------------------------
821	# Evolve system through time
822	#----------------------------------------------------------------
823	interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
824	gauge_values = []
825	for _ in interpolation_points:
826	gauge_values.append([])
827
828	time = []
829	for t in domain.evolve(yieldstep=0.1, finaltime=5.0):
830	# Record time series at known points
831	time.append(domain.get_time())
832
833	stage = domain.get_quantity('stage')
834	w = stage.get_values(interpolation_points=interpolation_points)
835
836	for i, _ in enumerate(interpolation_points):
837	gauge_values[i].append(w[i])
838
839	#Reference (nautilus 26/6/2008)
840
841	G0 = [-0.20000000000000001, -0.20000000000000001, -0.20000000000000001, -0.1958465301767274, -0.19059602372752493, -0.18448466250700923, -0.16979321333876071, -0.15976372740651074, -0.1575649333345176, -0.15710373731900584, -0.1530922283220747, -0.18836084336565725, -0.19921529311644628, -0.19923451799698919, -0.19923795414410964, -0.19923178806924047, -0.19925157557666154, -0.19930747801697429, -0.1993266428576112, -0.19932004930281799, -0.19929691326931867, -0.19926285267313795, -0.19916645449780995, -0.1988942593296438, -0.19900620256621993, -0.19914327423060865, -0.19918708440899577, -0.19921557252449132, -0.1992404368022069, -0.19925070370697717, -0.19925975477892274, -0.1992671090445659, -0.19927254203777162, -0.19927631910959256, -0.19927843552031504, -0.19927880339239365, -0.19927763204453783, -0.19927545249577633, -0.19927289590622824, -0.19927076261495152, -0.19926974334736983, -0.19927002562760332, -0.19927138236272213, -0.1992734501064522, -0.19927573251318192, -0.19927778936001547, -0.1992793776883893, -0.19928040577720926, -0.19928092586206753, -0.19928110982948721, -0.19928118887248453]
842
843	G1 = [-0.29999999999999993, -0.29999999999999993, -0.29139135018319512, -0.27257394456094503, -0.24141437432643265, -0.22089173942479151, -0.20796171092975532, -0.19874580192293825, -0.19014580508752857, -0.18421165368665365, -0.18020808282748838, -0.17518824759550247, -0.16436633464497749, -0.18714479115225544, -0.2045242886738807, -0.21011244240826329, -0.21151316017424124, -0.21048112933621732, -0.20772920477355789, -0.20489184334204144, -0.20286043930678221, -0.20094305756540246, -0.19948172752345467, -0.19886725178309209, -0.1986680808256765, -0.19860718133373548, -0.19862076543539733, -0.19888949069732539, -0.19932190310819023, -0.19982482967777454, -0.20036045468470615, -0.20064263130562704, -0.2007255389410077, -0.20068815669152493, -0.20057471332984647, -0.20042203940851802, -0.20026779013499779, -0.20015995671464712, -0.2000684005446525, -0.20001486753189174, -0.20000743467898013, -0.20003739771775905, -0.20008784600912621, -0.20013758305093884, -0.20017277554845025, -0.20018629233766885, -0.20018106288462198, -0.20016648079299326, -0.20015155958426531, -0.20014259747382254, -0.20014096648362509]
844
845
846	G2 = [-0.40000000000000002, -0.38885199453206343, -0.33425057028323962, -0.30154253721772117, -0.27624597383474103, -0.26037811196890087, -0.24415404585285019, -0.22941383121091052, -0.21613496492144549, -0.20418199946908885, -0.19506212965221825, -0.18851924999737435, -0.18271143344718843, -0.16910750701722474, -0.17963775224176018, -0.19442870510406052, -0.20164216917300118, -0.20467219452758528, -0.20608246104917802, -0.20640259931640445, -0.2054749739152594, -0.20380549124050265, -0.20227296931678532, -0.20095834856297176, -0.20000430919304379, -0.19946673053844086, -0.1990733499211611, -0.19882136174363013, -0.19877442300323914, -0.19905182154377868, -0.19943266521643804, -0.19988524183849191, -0.20018905307631765, -0.20031895675727809, -0.20033991149804931, -0.20031574232920274, -0.20027004750680638, -0.20020472427796293, -0.20013382447039607, -0.2000635018536408, -0.20001515436367642, -0.19998427691514989, -0.19997263083178107, -0.19998545383896535, -0.20000134502238734, -0.2000127243362736, -0.20001564474711939, -0.20001267360809977, -0.20002707579781318, -0.20004087961702843, -0.20004212947389177]
847
848	G3 = [-0.45000000000000001, -0.38058172993544187, -0.33756059941741273, -0.31015371357441396, -0.29214769368562965, -0.27545447937118606, -0.25871585649808154, -0.24254276680581988, -0.22758633129006092, -0.21417276895743134, -0.20237184768790789, -0.19369491041576814, -0.18721625333717057, -0.1794243868465818, -0.17052113574042196, -0.18534300640363346, -0.19601184621026671, -0.20185028431829469, -0.20476187496918136, -0.20602933256960082, -0.20598569228739247, -0.20492643756666848, -0.20339695828762758, -0.20196440373022595, -0.20070304052919338, -0.19986227854052355, -0.19933020476408528, -0.19898034831018035, -0.19878317651286193, -0.19886879323961787, -0.19915860801206181, -0.19953675278099042, -0.19992828019602107, -0.20015957043092364, -0.20025268671087426, -0.20028559516444974, -0.20027084877341045, -0.20022991487243985, -0.20016234295579871, -0.20009131445092507, -0.20003149397006523, -0.19998473356473795, -0.19996011913447218, -0.19995647168667186, -0.19996526209120422, -0.19996600297827097, -0.19997268800221216, -0.19998682275066659, -0.20000372259781876, -0.20001628681983963, -0.2000173314086407]
849
850	assert num.allclose(gauge_values[0], G0)
851	assert num.allclose(gauge_values[1], G1)
852	assert num.allclose(gauge_values[2], G2)
853	assert num.allclose(gauge_values[3], G3)
854
855	#####################################################
856
857
858	def test_initial_condition(self):
859	"""test_initial_condition
860
861	Test that initial condition is output at time == 0 and that
862	computed values change as system evolves
863	"""
864
865	from anuga.config import g
866	import copy
867
868	a = [0.0, 0.0]
869	b = [0.0, 2.0]
870	c = [2.0, 0.0]
871	d = [0.0, 4.0]
872	e = [2.0, 2.0]
873	f = [4.0, 0.0]
874
875	points = [a, b, c, d, e, f]
876	# bac, bce, ecf, dbe
877	vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
878
879	domain = Domain(points, vertices)
880
881	#Set up for a gradient of (3,0) at mid triangle (bce)
882	def slope(x, y):
883	return 3*x
884
885	h = 0.1
886	def stage(x, y):
887	return slope(x, y) + h
888
889	domain.set_quantity('elevation', slope)
890	domain.set_quantity('stage', stage)
891
892	# Allow slope limiters to work
893	# (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
894	domain.distribute_to_vertices_and_edges()
895
896	initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
897
898	domain.set_boundary({'exterior': Reflective_boundary(domain)})
899
900	domain.optimise_dry_cells = True
901
902	#Evolution
903	for t in domain.evolve(yieldstep=0.5, finaltime=2.0):
904	stage = domain.quantities['stage'].vertex_values
905
906	if t == 0.0:
907	assert num.allclose(stage, initial_stage)
908	else:
909	assert not num.allclose(stage, initial_stage)
910
911	os.remove(domain.get_name() + '.sww')
912
913	#####################################################
914
915	def test_second_order_flat_bed_onestep(self):
916	from mesh_factory import rectangular
917
918	#Create basic mesh
919	points, vertices, boundary = rectangular(6, 6)
920
921	#Create shallow water domain
922	domain = Domain(points, vertices, boundary)
923	domain.set_default_order(2)
924
925	# Boundary conditions
926	Br = Reflective_boundary(domain)
927	Bd = Dirichlet_boundary([0.1, 0., 0.])
928	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
929
930	domain.check_integrity()
931
932	# Evolution
933	for t in domain.evolve(yieldstep=0.05, finaltime=0.05):
934	pass
935
936
937	# Data from earlier version of abstract_2d_finite_volumes
938	assert num.allclose(domain.recorded_min_timestep, 0.0396825396825) or \
939	num.allclose(domain.recorded_min_timestep, 0.0235282801879)
940
941	assert num.allclose(domain.recorded_max_timestep, 0.0396825396825) or \
942	num.allclose(domain.recorded_max_timestep, 0.0235282801879)
943
944
945
946	assert num.allclose(domain.quantities['stage'].centroid_values[:12],
947	[0.00171396, 0.02656103, 0.00241523, 0.02656103,
948	0.00241523, 0.02656103, 0.00241523, 0.02656103,
949	0.00241523, 0.02656103, 0.00241523, 0.0272623], atol=1.0e-3) or \
950	num.allclose(domain.quantities['stage'].centroid_values[:12],
951	[ 0.00053119, 0.02900893, 0.00077912, 0.02900893,
952	0.00077912, 0.02900893, 0.00077912, 0.02900893,
953	0.00077912, 0.02900893, 0.00077912, 0.02873746], atol=1.0e-3)
954
955	domain.distribute_to_vertices_and_edges()
956
957
958
959	assert num.allclose(domain.quantities['stage'].vertex_values[:12,0],
960	[ -1.96794125e-03, 2.65610347e-02, 0.00000000e+00, 2.65610347e-02,
961	-8.67361738e-19, 2.65610347e-02, 4.33680869e-19, 2.65610347e-02,
962	-2.16840434e-18, 2.65610347e-02, -9.44042339e-05, 2.72623006e-02],
963	atol =1.0e-3) or \
964	num.allclose(domain.quantities['stage'].vertex_values[:12,0],
965	[ -5.51381419e-04, 5.74866732e-02, 1.00006808e-15, 5.72387383e-02,
966	9.99851243e-16, 5.72387383e-02, 1.00050176e-15, 5.72387383e-02,
967	9.99417563e-16, 5.72387383e-02, 1.09882029e-05, 5.66957956e-02],
968	atol=1.0e-3)
969
970
971	assert num.allclose(domain.quantities['stage'].vertex_values[:12,1],
972	[ 5.14188587e-03, 2.65610347e-02, 0.00000000e+00, 2.65610347e-02,
973	8.67361738e-19, 2.65610347e-02, -4.33680869e-19, 2.65610347e-02,
974	1.30104261e-18, 2.65610347e-02, 9.44042339e-05, 2.72623006e-02],
975	atol =1.0e-3) or \
976	num.allclose(domain.quantities['stage'].vertex_values[:12,1],
977	[ 1.59356551e-03, 5.72387383e-02, 1.00006808e-15, 5.72387383e-02,
978	1.00006808e-15, 5.72387383e-02, 9.99634403e-16, 5.72387383e-02,
979	1.00050176e-15, 5.72387383e-02, -1.09882029e-05, 1.47582915e-02],
980	atol =1.0e-3)
981
982	assert num.allclose(domain.quantities['stage'].vertex_values[:12,2],
983	[ 0.00196794, 0.02656103, 0.00724568, 0.02656103,
984	0.00724568, 0.02656103, 0.00724568, 0.02656103,
985	0.00724568, 0.02656103, 0.00724568, 0.0272623 ], atol =1.0e-3) or \
986	num.allclose(domain.quantities['stage'].vertex_values[:12,2],
987	[ 0.00055138, -0.02769862, 0.00233737, -0.02745068,
988	0.00233737, -0.02745068, 0.00233737, -0.02745068,
989	0.00233737, -0.02745068, 0.00233737, 0.01475829], atol =1.0e-3)
990
991
992	assert num.allclose(domain.quantities['xmomentum'].centroid_values[:12],
993	[0.00113961, 0.01302432, 0.00148672,
994	0.01302432, 0.00148672, 0.01302432,
995	0.00148672, 0.01302432, 0.00148672 ,
996	0.01302432, 0.00148672, 0.01337143], atol=1.0e-3) or \
997	num.allclose(domain.quantities['xmomentum'].centroid_values[:12],
998	[ 0.00019529, 0.01425863, 0.00025665,
999	0.01425863, 0.00025665, 0.01425863,
1000	0.00025665, 0.01425863, 0.00025665,
1001	0.01425863, 0.00025665, 0.014423 ], atol=1.0e-3)
1002
1003	assert num.allclose(domain.quantities['ymomentum'].centroid_values[:12],
1004	[-2.91240050e-004 , 1.22721531e-004,
1005	-1.22721531e-004, 1.22721531e-004 ,
1006	-1.22721531e-004, 1.22721531e-004,
1007	-1.22721531e-004 , 1.22721531e-004,
1008	-1.22721531e-004, 1.22721531e-004,
1009	-1.22721531e-004, -4.57969873e-005], atol=1.0e-5) or \
1010	num.allclose(domain.quantities['ymomentum'].centroid_values[:12],
1011	[ -6.38239364e-05, 2.16943067e-05,
1012	-2.16943067e-05, 2.16943067e-05,
1013	-2.16943067e-05, 2.16943067e-05,
1014	-2.16943067e-05, 2.16943067e-05,
1015	-2.16943067e-05, 2.16943067e-05,
1016	-2.16943067e-05, -4.62796434e-04], atol=1.0e-5)
1017
1018	os.remove(domain.get_name() + '.sww')
1019
1020	def test_second_order_flat_bed_moresteps(self):
1021	from mesh_factory import rectangular
1022
1023	# Create basic mesh
1024	points, vertices, boundary = rectangular(6, 6)
1025
1026	# Create shallow water domain
1027	domain = Domain(points, vertices, boundary)
1028	domain.smooth = False
1029	domain.default_order = 2
1030
1031	# Boundary conditions
1032	Br = Reflective_boundary(domain)
1033	Bd = Dirichlet_boundary([0.1, 0., 0.])
1034	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1035
1036	domain.check_integrity()
1037
1038	# Evolution
1039	for t in domain.evolve(yieldstep=0.05, finaltime=0.1):
1040	pass
1041
1042	# Data from earlier version of abstract_2d_finite_volumes
1043	#assert allclose(domain.recorded_min_timestep, 0.0396825396825)
1044	#assert allclose(domain.recorded_max_timestep, 0.0396825396825)
1045	#print domain.quantities['stage'].centroid_values
1046
1047	os.remove(domain.get_name() + '.sww')
1048
1049	def test_flatbed_first_order(self):
1050	from mesh_factory import rectangular
1051
1052	# Create basic mesh
1053	N = 8
1054	points, vertices, boundary = rectangular(N, N)
1055
1056	# Create shallow water domain
1057	domain = Domain(points, vertices, boundary)
1058	domain.smooth = False
1059	domain.default_order = 1
1060	domain.H0 = 1.0e-3 # As suggested in the manual
1061
1062	# Boundary conditions
1063	Br = Reflective_boundary(domain)
1064	Bd = Dirichlet_boundary([0.2, 0., 0.])
1065
1066	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1067	domain.check_integrity()
1068
1069	# Evolution
1070	for t in domain.evolve(yieldstep=0.02, finaltime=0.5):
1071	pass
1072
1073	# FIXME: These numbers were from version before 25/10
1074	#assert allclose(domain.recorded_min_timestep, 0.0140413643926)
1075	#assert allclose(domain.recorded_max_timestep, 0.0140947355753)
1076
1077	for i in range(3):
1078	#assert allclose(domain.quantities['stage'].edge_values[:4,i],
1079	# [0.10730244,0.12337617,0.11200126,0.12605666])
1080	assert num.allclose(domain.quantities['xmomentum'].\
1081	edge_values[:4,i],
1082	[0.07610894,0.06901572,0.07284461,0.06819712])
1083	#assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
1084	# [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
1085
1086	os.remove(domain.get_name() + '.sww')
1087
1088	def test_flatbed_second_order(self):
1089	from mesh_factory import rectangular
1090
1091	# Create basic mesh
1092	N = 8
1093	points, vertices, boundary = rectangular(N, N)
1094
1095	# Create shallow water domain
1096	domain = Domain(points, vertices, boundary)
1097	domain.set_store_vertices_uniquely(True)
1098	domain.set_default_order(2)
1099
1100	# Boundary conditions
1101	Br = Reflective_boundary(domain)
1102	Bd = Dirichlet_boundary([0.2, 0., 0.])
1103
1104	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1105	domain.check_integrity()
1106
1107	# Evolution
1108	for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
1109	pass
1110
1111
1112
1113	msg = 'min step was %f instead of %f' % (domain.recorded_min_timestep,
1114	0.0155604907816)
1115
1116	assert num.allclose(domain.recorded_min_timestep, 0.0155604907816), msg
1117	assert num.allclose(domain.recorded_max_timestep, 0.0155604907816)
1118
1119
1120	assert num.allclose(domain.quantities['stage'].vertex_values[:4,0],
1121	[-0.009, 0.0535, 0.0, 0.0535], atol=1.0e-3) or \
1122	num.allclose(domain.quantities['stage'].vertex_values[:4,0],
1123	[-3.54158995e-03,1.22050959e-01,-2.36227400e-05,1.21501627e-01], atol=1.0e-3)
1124
1125
1126	assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1127	[-0.008, 0.0368, 0.0, 0.0368], atol=1.0e-3) or \
1128	num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1129	[-2.32056226e-03,9.10618822e-02, -1.06135035e-05,9.75175956e-02], atol=1.0e-3)
1130
1131	assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1132	[ 0.002 , 6.0e-04, 0.0, 6.0e-04],
1133	atol=1.0e-3) or \
1134	num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1135	[ 1.43500775e-03, 6.07102924e-05, 1.59329371e-06, 8.44572599e-03],
1136	atol=1.0e-3)
1137
1138	os.remove(domain.get_name() + '.sww')
1139
1140
1141	def test_flatbed_second_order_vmax_0(self):
1142	from mesh_factory import rectangular
1143
1144	# Create basic mesh
1145	N = 8
1146	points, vertices, boundary = rectangular(N, N)
1147
1148	# Create shallow water domain
1149	domain = Domain(points, vertices, boundary)
1150
1151	domain.set_store_vertices_uniquely(True)
1152	domain.set_default_order(2)
1153
1154
1155	# Boundary conditions
1156	Br = Reflective_boundary(domain)
1157	Bd = Dirichlet_boundary([0.2, 0., 0.])
1158
1159	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1160	domain.check_integrity()
1161
1162	# Evolution
1163	for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
1164	pass
1165
1166
1167	assert num.allclose(domain.recorded_min_timestep, 0.0210448446782) or \
1168	num.allclose(domain.recorded_min_timestep, 0.0155604907816)
1169
1170	assert num.allclose(domain.recorded_max_timestep, 0.0210448446782) or \
1171	num.allclose(domain.recorded_min_timestep, 0.0155604907816)
1172
1173
1174	assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1175	[ -2.32056226e-03, 9.10618822e-02, -1.06135035e-05, 9.75175956e-02],
1176	atol=1.0e-3)
1177
1178	assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1179	[ 1.43500775e-03, 6.07102924e-05, 1.59329371e-06, 8.44572599e-03],
1180	atol=1.0e-3)
1181
1182	os.remove(domain.get_name() + '.sww')
1183
1184	def test_flatbed_second_order_distribute(self):
1185	#Use real data from anuga.abstract_2d_finite_volumes 2
1186	#painfully setup and extracted.
1187
1188	from mesh_factory import rectangular
1189
1190	# Create basic mesh
1191	N = 8
1192	points, vertices, boundary = rectangular(N, N)
1193
1194	# Create shallow water domain
1195	domain = Domain(points, vertices, boundary)
1196
1197	domain.set_store_vertices_uniquely(True)
1198	domain.set_default_order(2)
1199
1200	# Boundary conditions
1201	Br = Reflective_boundary(domain)
1202	Bd = Dirichlet_boundary([0.2, 0., 0.])
1203
1204	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
1205	domain.check_integrity()
1206
1207	for V in [False, True]:
1208	if V:
1209	# Set centroids as if system had been evolved
1210	L = num.zeros(2NN, num.float)
1211	L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
1212	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1213	1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
1214	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
1215	1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
1216	5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
1217	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1218	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1219	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
1220	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
1221	0.00000000e+000, 5.57305948e-005]
1222
1223	X = num.zeros(2NN, num.float)
1224	X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
1225	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1226	8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
1227	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
1228	8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
1229	3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
1230	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1231	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1232	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
1233	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
1234	0.00000000e+000, 4.57662812e-005]
1235
1236	Y = num.zeros(2NN, num.float)
1237	Y[:32] = [-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
1238	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1239	-6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
1240	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
1241	-6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
1242	-1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
1243	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1244	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1245	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
1246	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
1247	0.00000000e+000, -2.57635178e-005]
1248
1249	domain.set_quantity('stage', L, location='centroids')
1250	domain.set_quantity('xmomentum', X, location='centroids')
1251	domain.set_quantity('ymomentum', Y, location='centroids')
1252
1253	domain.check_integrity()
1254	else:
1255	# Evolution
1256	for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
1257	pass
1258
1259
1260	assert num.allclose(domain.recorded_min_timestep, 0.0155604907816)
1261	assert num.allclose(domain.recorded_max_timestep, 0.0155604907816)
1262
1263	#print domain.quantities['stage'].centroid_values[:4]
1264	#print domain.quantities['xmomentum'].centroid_values[:4]
1265	#print domain.quantities['ymomentum'].centroid_values[:4]
1266
1267	#Centroids were correct but not vertices.
1268	#Hence the check of distribute below.
1269
1270	if not V:
1271
1272	assert num.allclose(domain.quantities['stage'].centroid_values[:4],
1273	[0.00725574, 0.05350737, 0.01008413, 0.0535293], atol=1.0e-3) or \
1274	num.allclose(domain.quantities['stage'].centroid_values[:4],
1275	[0.00318259, 0.06261678, 0.00420215, 0.06285189], atol=1.0e-3)
1276
1277	assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
1278	[0.00654964, 0.03684904, 0.00852561, 0.03686323],atol=1.0e-3) or \
1279	num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
1280	[0.00218173, 0.04482164, 0.0026334, 0.04491656],atol=1.0e-3)
1281
1282	assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
1283	[-0.00143169, 0.00061027, -0.00062325, 0.00061408],atol=1.0e-3) or \
1284	num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
1285	[-6.46340592e-04,-6.16702557e-05,-2.83424134e-04, 6.48556590e-05],atol=1.0e-3)
1286
1287
1288	assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0,
1289	atol=3.0e-4)
1290	else:
1291	assert num.allclose(domain.quantities['xmomentum'].\
1292	centroid_values[17],
1293	0.00028606084)
1294	return #FIXME - Bailout for V True
1295
1296	import copy
1297
1298	XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
1299	assert num.allclose(domain.quantities['xmomentum'].centroid_values,
1300	XX)
1301
1302	domain.distribute_to_vertices_and_edges()
1303
1304	assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0, atol=3.0e-4)
1305
1306	assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1307	[ 1.84104149e-03, 6.05658846e-04, 1.77092716e-07, 6.10687334e-04],
1308	atol=1.0e-4) or \
1309	num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
1310	[1.43500775e-03, 6.07102924e-05, 1.59329371e-06, 8.44572599e-03],
1311	atol=1.0e-4)
1312
1313	assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1314	[ -8.31184293e-03, 3.68841505e-02, -2.42843889e-06, 3.68900189e-02],
1315	atol=1.0e-4) or \
1316	num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
1317	[-2.32056226e-03, 9.10618822e-02, -1.06135035e-05, 9.75175956e-02],
1318	rtol=1.0e-2)
1319
1320
1321	os.remove(domain.get_name() + '.sww')
1322
1323
1324	def test_bedslope_problem_second_order_more_steps(self):
1325	"""test_bedslope_problem_second_order_more_step
1326
1327	Test shallow water balanced finite volumes
1328	"""
1329
1330	from mesh_factory import rectangular
1331
1332	# Create basic mesh
1333	points, vertices, boundary = rectangular(6, 6)
1334
1335	# Create shallow water domain
1336	domain = Domain(points, vertices, boundary)
1337
1338	domain.set_store_vertices_uniquely(True)
1339	domain.set_default_order(2)
1340
1341	# Bed-slope and friction at vertices (and interpolated elsewhere)
1342	def x_slope(x, y):
1343	return -x/3
1344
1345	domain.set_quantity('elevation', x_slope)
1346
1347	# Boundary conditions
1348	Br = Reflective_boundary(domain)
1349	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1350
1351	# Initial condition
1352	domain.set_quantity('stage', expression='elevation+0.05')
1353	domain.check_integrity()
1354
1355	assert num.allclose(domain.quantities['stage'].centroid_values,
1356	[ 0.01296296, 0.03148148, 0.01296296,
1357	0.03148148, 0.01296296, 0.03148148,
1358	0.01296296, 0.03148148, 0.01296296,
1359	0.03148148, 0.01296296, 0.03148148,
1360	-0.04259259, -0.02407407, -0.04259259,
1361	-0.02407407, -0.04259259, -0.02407407,
1362	-0.04259259, -0.02407407, -0.04259259,
1363	-0.02407407, -0.04259259, -0.02407407,
1364	-0.09814815, -0.07962963, -0.09814815,
1365	-0.07962963, -0.09814815, -0.07962963,
1366	-0.09814815, -0.07962963, -0.09814815,
1367	-0.07962963, -0.09814815, -0.07962963,
1368	-0.1537037 , -0.13518519, -0.1537037,
1369	-0.13518519, -0.1537037, -0.13518519,
1370	-0.1537037 , -0.13518519, -0.1537037,
1371	-0.13518519, -0.1537037, -0.13518519,
1372	-0.20925926, -0.19074074, -0.20925926,
1373	-0.19074074, -0.20925926, -0.19074074,
1374	-0.20925926, -0.19074074, -0.20925926,
1375	-0.19074074, -0.20925926, -0.19074074,
1376	-0.26481481, -0.2462963, -0.26481481,
1377	-0.2462963, -0.26481481, -0.2462963,
1378	-0.26481481, -0.2462963, -0.26481481,
1379	-0.2462963, -0.26481481, -0.2462963])
1380
1381	# Evolution
1382	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
1383	pass
1384
1385
1386	assert num.allclose(domain.quantities['stage'].centroid_values,
1387	[-0.02901283, -0.01619385, -0.03040423, -0.01564474, -0.02936756, -0.01507953,
1388	-0.02858108, -0.01491531, -0.02793549, -0.0147037, -0.02792804, -0.014363,
1389	-0.07794301, -0.05951952, -0.07675098, -0.06182336, -0.07736607, -0.06079504,
1390	-0.07696764, -0.06039043, -0.07708793, -0.0601453, -0.07669911, -0.06020719,
1391	-0.12223185, -0.10857309, -0.12286676, -0.10837591, -0.12386938, -0.10842744,
1392	-0.12363769, -0.10790002, -0.12304837, -0.10871278, -0.12543768, -0.10961026,
1393	-0.15664473, -0.14630207, -0.15838364, -0.14910271, -0.15804002, -0.15029627,
1394	-0.15829717, -0.1503869, -0.15852604, -0.14971109, -0.15856346, -0.15205092,
1395	-0.20900931, -0.19658843, -0.20669607, -0.19558708, -0.20654186, -0.19492423,
1396	-0.20438765, -0.19492931, -0.20644142, -0.19423147, -0.20237449, -0.19198454,
1397	-0.13699658, -0.14209126, -0.13600697, -0.14334968, -0.1347657, -0.14224247,
1398	-0.13442376, -0.14136926, -0.13501004, -0.14339389, -0.13479263, -0.14304073], atol=1.0e-2) or \
1399	num.allclose(domain.quantities['stage'].centroid_values,
1400	[-0.03393968, -0.0166423, -0.03253538, -0.01722023, -0.03270405, -0.01728606,
1401	-0.03277786, -0.0173903, -0.03333736, -0.01743236, -0.03189526, -0.01463918,
1402	-0.07951756, -0.06410763, -0.07847973, -0.06350794, -0.07842429, -0.06240852,
1403	-0.07808697, -0.06255924, -0.07854662, -0.06322442, -0.07867314, -0.06287121,
1404	-0.11533356, -0.10559238, -0.11971301, -0.10742123, -0.1215759 , -0.10830046,
1405	-0.12202867, -0.10831703, -0.122214, -0.10854099, -0.12343779, -0.11035803,
1406	-0.15725714, -0.14300757, -0.15559898, -0.1447275 , -0.15478568, -0.14483551,
1407	-0.15461918, -0.14489704, -0.15462074, -0.14516256, -0.15522298, -0.1452902,
1408	-0.22637615, -0.19192974, -0.20922654, -0.1907441 , -0.20900039, -0.19074809,
1409	-0.20897969, -0.19073365, -0.209195, -0.19071396, -0.20922513, -0.19067714,
1410	-0.11357515, -0.14185801, -0.13224763, -0.14395805, -0.13379438, -0.14497114,
1411	-0.13437773, -0.14536013, -0.13607796, -0.14799629, -0.13148351, -0.15568502], atol=1.0e-1)
1412
1413
1414
1415	assert num.allclose(domain.quantities['xmomentum'].centroid_values,
1416	[ 0.00478273, 0.003297, 0.00471129, 0.00320957, 0.00462171, 0.00320135,
1417	0.00458295, 0.00317193, 0.00451704, 0.00314308, 0.00442684, 0.00320466,
1418	0.01512907, 0.01150756, 0.01604672, 0.01156605, 0.01583911, 0.01135809,
1419	0.01578499, 0.01132479, 0.01543668, 0.01100614, 0.01570445, 0.0120152,
1420	0.04019477, 0.02721469, 0.03509982, 0.02735229, 0.03369315, 0.02727871,
1421	0.03317931, 0.02706421, 0.03332704, 0.02722779, 0.03170258, 0.02556134,
1422	0.07157025, 0.06074271, 0.07249738, 0.05570979, 0.07311261, 0.05428175,
1423	0.07316986, 0.05379702, 0.0719581, 0.05230996, 0.07034837, 0.05468702,
1424	0.08145001, 0.07753479, 0.08148804, 0.08119069, 0.08247295, 0.08134969,
1425	0.0823216, 0.081411, 0.08190964, 0.08151441, 0.08163076, 0.08166174,
1426	0.03680205, 0.0768216, 0.03943625, 0.07791183, 0.03930529, 0.07760588,
1427	0.03949756, 0.07839929, 0.03992892, 0.08001416, 0.04444335, 0.08628738],
1428	atol=1.0e-2) or \
1429	num.allclose(domain.quantities['xmomentum'].centroid_values,
1430	[ 0.00178414, 0.00147791, 0.00373636, 0.00169124, 0.00395649, 0.0014468,
1431	0.00387617, 0.00135572, 0.00338418, 0.00134554, 0.00404961, 0.00252769,
1432	0.01365204, 0.00890416, 0.01381613, 0.00986246, 0.01419385, 0.01145017,
1433	0.01465116, 0.01125933, 0.01407359, 0.01055426, 0.01403563, 0.01095544,
1434	0.04653827, 0.03018236, 0.03709973, 0.0265533 , 0.0337694 , 0.02541724,
1435	0.03304266, 0.02535335, 0.03264548, 0.02484769, 0.03047682, 0.02205757,
1436	0.07400338, 0.06470583, 0.07756503, 0.06098108, 0.07942593, 0.06086531,
1437	0.07977427, 0.06074404, 0.07979513, 0.06019911, 0.07806395, 0.06011152,
1438	0.07305045, 0.07883894, 0.08120393, 0.08166623, 0.08180501, 0.08166251,
1439	0.0818353 , 0.08169641, 0.08173762, 0.08174118, 0.08176467, 0.08181817,
1440	0.01549926, 0.08259719, 0.01835423, 0.07302656, 0.01672924, 0.07198839,
1441	0.01676006, 0.07223233, 0.01775672, 0.07362164, 0.01955846, 0.09361223],
1442	atol=1.0e-2)
1443
1444
1445	assert num.allclose(domain.quantities['ymomentum'].centroid_values,
1446	[ -1.09771684e-05, -2.60328801e-05, -1.03481959e-05, -7.75907380e-05,
1447	-5.00409090e-05, -7.83807512e-05, -3.60509918e-05, -6.19321031e-05,
1448	-1.40041903e-05, -2.95707259e-05, 3.90296618e-06, 1.87556544e-05,
1449	9.27848053e-05, 6.66937557e-07, 1.00653468e-04, 8.24734209e-06,
1450	-1.04548672e-05, -4.40402988e-05, -2.95549946e-05, -1.86360736e-05,
1451	1.12527016e-04, 1.27240681e-04, 2.02147284e-04, 9.18457482e-05,
1452	1.41781748e-03, 7.23407624e-04, 5.09160779e-04, 1.29136939e-04,
1453	-4.70131286e-05, -1.00180290e-04, -1.76806614e-05, -4.19421384e-06,
1454	-6.17759681e-05, -3.02124967e-05, 4.32689360e-04, 5.49717934e-04,
1455	1.15031101e-03, 1.02737170e-03, 5.77937840e-04, 3.36230967e-04,
1456	5.44877516e-04, -7.28594977e-05, 4.60064858e-04, -3.94125434e-05,
1457	7.48242964e-04, 2.88528341e-04, 6.25148041e-05, -1.74477175e-04,
1458	-5.06603166e-05, 7.07720999e-04, -2.04937748e-04, 3.38595573e-05,
1459	-4.64116229e-05, 1.49325340e-04, -2.41342281e-05, 1.83817970e-04,
1460	-1.44417277e-05, 2.47823834e-04, 7.91185571e-05, 1.71615793e-04,
1461	1.56883043e-03, 8.39352974e-04, 3.23353846e-03, 1.70597880e-03,
1462	2.27789107e-03, 1.48928169e-03, 2.09854126e-03, 1.50248643e-03,
1463	2.83029467e-03, 1.09151499e-03, 6.52455118e-03, -2.04468968e-03],
1464	atol=1.0e-3) or \
1465	num.allclose(domain.quantities['ymomentum'].centroid_values,
1466	[ -1.24810991e-04, -3.08228767e-04, -1.56701128e-04, -1.01904208e-04,
1467	-3.36282053e-05, -1.17956840e-04, -3.55986664e-05, -9.38578996e-05,
1468	7.13704069e-05, 2.47022380e-05, 1.71121489e-04, 2.65941677e-04,
1469	6.90055205e-04, 1.99195585e-04, 1.33804448e-04, -1.66563316e-04,
1470	-2.00962830e-04, -3.81664130e-05, -9.50456053e-05, -3.14620186e-06,
1471	1.29388102e-04, 3.16945980e-04, 4.77556581e-04, 2.57217342e-04,
1472	1.42300612e-03, 9.60776359e-04, 5.08941026e-04, 1.06939990e-04,
1473	6.37673950e-05, -2.69783047e-04, -8.55760509e-05, -2.12987309e-04,
1474	-5.86840949e-06, -9.75751293e-05, 8.25447727e-04, 1.14139065e-03,
1475	8.56206468e-04, 3.83113329e-04, 1.75041847e-04, 4.39999200e-04,
1476	3.75156469e-04, 2.48774698e-04, 4.09671654e-04, 2.07125615e-04,
1477	4.59587647e-04, 2.70581830e-04, -1.24082302e-06, -4.29155678e-04,
1478	-9.66841218e-03, 4.93278794e-04, -5.25778806e-06, -4.90396857e-05,
1479	-9.75373988e-06, 7.28023591e-06, -5.20499868e-06, 3.61013683e-05,
1480	-7.54919544e-06, 4.14115771e-05, -1.35778834e-05, -2.23991903e-05,
1481	3.63635844e-02, 5.29865244e-04, 5.13015379e-03, 1.19233296e-03,
1482	4.70681275e-04, 2.62292296e-04, -1.28084045e-04, 7.04826916e-04,
1483	1.50377987e-04, 1.35053814e-03, 1.30710492e-02, 1.93011958e-03],
1484	atol=1.0e-1)
1485
1486
1487
1488	os.remove(domain.get_name() + '.sww')
1489
1490
1491	def test_temp_play(self):
1492	from mesh_factory import rectangular
1493
1494	# Create basic mesh
1495	points, vertices, boundary = rectangular(5, 5)
1496
1497	# Create shallow water domain
1498	domain = Domain(points, vertices, boundary)
1499
1500	domain.set_store_vertices_uniquely(True)
1501	domain.set_default_order(2)
1502
1503
1504
1505	# Bed-slope and friction at vertices (and interpolated elsewhere)
1506	def x_slope(x, y):
1507	return -x/3
1508
1509	domain.set_quantity('elevation', x_slope)
1510
1511	# Boundary conditions
1512	Br = Reflective_boundary(domain)
1513	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
1514
1515	# Initial condition
1516	domain.set_quantity('stage', expression='elevation+0.05')
1517	domain.check_integrity()
1518
1519	# Evolution
1520	for t in domain.evolve(yieldstep=0.05, finaltime=0.1):
1521	pass
1522
1523
1524	assert num.allclose(domain.quantities['stage'].centroid_values[:4],
1525	[ 0.01, 0.015, 0.01, 0.015], atol=1.0e-2)
1526
1527	assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
1528	[ 0.015, 0.01, 0.015, 0.01], atol=1.0e-2)
1529
1530	assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
1531	[ 0.0, 0.0, 0.0, 0.0]
1532	, atol=1.0e-3)
1533
1534	os.remove(domain.get_name() + '.sww')
1535
1536	def test_complex_bed(self):
1537	# No friction is tested here
1538
1539	from mesh_factory import rectangular
1540
1541	N = 12
1542	points, vertices, boundary = rectangular(N, N/2, len1=1.2, len2=0.6,
1543	origin=(-0.07, 0))
1544
1545
1546	domain = Domain(points, vertices, boundary)
1547	domain.smooth = False
1548	domain.set_default_order(2)
1549	domain.set_timestepping_method('rk2')
1550	domain.set_beta(1.0)
1551
1552	inflow_stage = 0.1
1553	Z = Weir(inflow_stage)
1554	domain.set_quantity('elevation', Z)
1555
1556	Br = Reflective_boundary(domain)
1557	Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
1558	domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
1559
1560	domain.set_quantity('stage', expression='elevation')
1561
1562	for t in domain.evolve(yieldstep=0.02, finaltime=0.2):
1563	pass
1564
1565	#FIXME: These numbers were from version before 25/10
1566	#assert allclose(domain.quantities['stage'].centroid_values,
1567	# [3.95822638e-002, 5.61022588e-002, 4.66437868e-002, 5.73081011e-002,
1568	# 4.72394613e-002, 5.74684939e-002, 4.74309483e-002, 5.77458084e-002,
1569	# 4.80628177e-002, 5.85656225e-002, 4.90498542e-002, 6.02609831e-002,
1570	# 1.18470315e-002, 1.75136443e-002, 1.18035266e-002, 2.15565695e-002,
1571	# 1.31620268e-002, 2.14351640e-002, 1.32351076e-002, 2.15450687e-002,
1572	# 1.36414028e-002, 2.24274619e-002, 1.51689511e-002, 2.21789655e-002,
1573	# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
1574	# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
1575	# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003, 2.73985295e-004,
1576	# 1.86166683e-001, 8.74070369e-002, 1.86166712e-001, 8.74035875e-002,
1577	# 6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
1578	# 6.11666725e-002, 8.73846774e-002, 1.86166697e-001, 8.74171550e-002,
1579	# -4.83333333e-002, 1.18333333e-001, -4.83333333e-002, 1.18333333e-001,
1580	# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
1581	# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002, 1.18333333e-001,
1582	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1583	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1584	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
1585	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1586	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1587	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
1588	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1589	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1590	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
1591	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1592	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1593	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
1594	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1595	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1596	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
1597	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1598	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1599	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
1600	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1601	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
1602	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
1603
1604	os.remove(domain.get_name() + '.sww')
1605
1606
1607	def test_tight_slope_limiters(self):
1608	"""Test that new slope limiters (Feb 2007) don't induce extremely
1609	small timesteps. This test actually reveals the problem as it
1610	was in March-April 2007
1611	"""
1612	import time, os
1613	from Scientific.IO.NetCDF import NetCDFFile
1614	from data_manager import extent_sww
1615	from mesh_factory import rectangular_cross
1616
1617	# Create basic mesh
1618	points, vertices, boundary = rectangular_cross(2, 2)
1619
1620	# Create shallow water domain
1621	domain = Domain(points, vertices, boundary)
1622	domain.set_default_order(2)
1623	domain.set_beta(1.0)
1624	domain.set_timestepping_method('euler')
1625	#domain.set_CFL(0.5)
1626
1627
1628	# This will pass
1629	#domain.tight_slope_limiters = 1
1630	#domain.H0 = 0.01
1631
1632	# This will fail
1633	#domain.tight_slope_limiters = 1
1634	#domain.H0 = 0.001
1635
1636	# This will pass provided C extension implements limiting of
1637	# momentum in _compute_speeds
1638	#domain.tight_slope_limiters = 1
1639	#domain.H0 = 0.001
1640	#domain.protect_against_isolated_degenerate_timesteps = True
1641
1642	# Set some field values
1643	domain.set_quantity('elevation', lambda x,y: -x)
1644	domain.set_quantity('friction', 0.03)
1645
1646	# Boundary conditions
1647	B = Transmissive_boundary(domain)
1648	domain.set_boundary({'left': B, 'right': B, 'top': B, 'bottom': B})
1649
1650	# Initial condition - with jumps
1651	bed = domain.quantities['elevation'].vertex_values
1652	stage = num.zeros(bed.shape, num.float)
1653
1654	h = 0.3
1655	for i in range(stage.shape[0]):
1656	if i % 2 == 1:
1657	stage[i,:] = bed[i,:] + h
1658	else:
1659	stage[i,:] = bed[i,:]
1660
1661	domain.set_quantity('stage', stage)
1662
1663	domain.distribute_to_vertices_and_edges()
1664
1665	domain.set_name('tight_limiters')
1666	domain.smooth = True
1667	domain.reduction = mean
1668	domain.set_datadir('.')
1669	domain.smooth = False
1670	domain.store = True
1671
1672	# Evolution
1673	for t in domain.evolve(yieldstep=0.1, finaltime=0.3):
1674	#domain.write_time(track_speeds=True)
1675	stage = domain.quantities['stage'].vertex_values
1676
1677	# Get NetCDF
1678	#fid = NetCDFFile(domain.writer.filename, netcdf_mode_r)
1679	#stage_file = fid.variables['stage']
1680
1681	#fid.close()
1682
1683	os.remove(domain.writer.filename)
1684
1685
1686
1687	def test_pmesh2Domain(self):
1688	import os
1689	import tempfile
1690
1691	fileName = tempfile.mktemp(".tsh")
1692	file = open(fileName, "w")
1693	file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
1694	0 0.0 0.0 0.0 0.0 0.01 \n \
1695	1 1.0 0.0 10.0 10.0 0.02 \n \
1696	2 0.0 1.0 0.0 10.0 0.03 \n \
1697	3 0.5 0.25 8.0 12.0 0.04 \n \
1698	# Vert att title \n \
1699	elevation \n \
1700	stage \n \
1701	friction \n \
1702	2 # <triangle #> [<vertex #>] [<neigbouring triangle #>] \n\
1703	0 0 3 2 -1 -1 1 dsg\n\
1704	1 0 1 3 -1 0 -1 ole nielsen\n\
1705	4 # <segment #> <vertex #> <vertex #> [boundary tag] \n\
1706	0 1 0 2 \n\
1707	1 0 2 3 \n\
1708	2 2 3 \n\
1709	3 3 1 1 \n\
1710	3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
1711	0 216.0 -86.0 \n \
1712	1 160.0 -167.0 \n \
1713	2 114.0 -91.0 \n \
1714	3 # <vertex #> <vertex #> [boundary tag] ...Mesh Segments... \n \
1715	0 0 1 0 \n \
1716	1 1 2 0 \n \
1717	2 2 0 0 \n \
1718	0 # <x> <y> ...Mesh Holes... \n \
1719	0 # <x> <y> <attribute>...Mesh Regions... \n \
1720	0 # <x> <y> <attribute>...Mesh Regions, area... \n\
1721	#Geo reference \n \
1722	56 \n \
1723	140 \n \
1724	120 \n")
1725	file.close()
1726
1727	tags = {}
1728	b1 = Dirichlet_boundary(conserved_quantities = num.array([0.0]))
1729	b2 = Dirichlet_boundary(conserved_quantities = num.array([1.0]))
1730	b3 = Dirichlet_boundary(conserved_quantities = num.array([2.0]))
1731	tags["1"] = b1
1732	tags["2"] = b2
1733	tags["3"] = b3
1734
1735	domain = Domain(mesh_filename=fileName)
1736	# verbose=True, use_cache=True)
1737
1738	## check the quantities
1739	answer = [[0., 8., 0.],
1740	[0., 10., 8.]]
1741	assert num.allclose(domain.quantities['elevation'].vertex_values,
1742	answer)
1743
1744	answer = [[0., 12., 10.],
1745	[0., 10., 12.]]
1746	assert num.allclose(domain.quantities['stage'].vertex_values,
1747	answer)
1748
1749	answer = [[0.01, 0.04, 0.03],
1750	[0.01, 0.02, 0.04]]
1751	assert num.allclose(domain.quantities['friction'].vertex_values,
1752	answer)
1753
1754	tagged_elements = domain.get_tagged_elements()
1755	assert num.allclose(tagged_elements['dsg'][0], 0)
1756	assert num.allclose(tagged_elements['ole nielsen'][0], 1)
1757
1758	msg = "test_tags_to_boundaries failed. Single boundary wasn't added."
1759	self.failUnless( domain.boundary[(1, 0)] == '1', msg)
1760	self.failUnless( domain.boundary[(1, 2)] == '2', msg)
1761	self.failUnless( domain.boundary[(0, 1)] == '3', msg)
1762	self.failUnless( domain.boundary[(0, 0)] == 'exterior', msg)
1763	msg = "test_pmesh2Domain Too many boundaries"
1764	self.failUnless( len(domain.boundary) == 4, msg)
1765
1766	# FIXME change to use get_xllcorner
1767	msg = 'Bad geo-reference'
1768	self.failUnless(domain.geo_reference.xllcorner == 140.0, msg)
1769
1770	domain = Domain(fileName)
1771
1772	answer = [[0., 8., 0.],
1773	[0., 10., 8.]]
1774	assert num.allclose(domain.quantities['elevation'].vertex_values,
1775	answer)
1776
1777	answer = [[0., 12., 10.],
1778	[0., 10., 12.]]
1779	assert num.allclose(domain.quantities['stage'].vertex_values,
1780	answer)
1781
1782	answer = [[0.01, 0.04, 0.03],
1783	[0.01, 0.02, 0.04]]
1784	assert num.allclose(domain.quantities['friction'].vertex_values,
1785	answer)
1786
1787	tagged_elements = domain.get_tagged_elements()
1788	assert num.allclose(tagged_elements['dsg'][0], 0)
1789	assert num.allclose(tagged_elements['ole nielsen'][0], 1)
1790
1791	msg = "test_tags_to_boundaries failed. Single boundary wasn't added."
1792	self.failUnless(domain.boundary[(1, 0)] == '1', msg)
1793	self.failUnless(domain.boundary[(1, 2)] == '2', msg)
1794	self.failUnless(domain.boundary[(0, 1)] == '3', msg)
1795	self.failUnless(domain.boundary[(0, 0)] == 'exterior', msg)
1796	msg = "test_pmesh2Domain Too many boundaries"
1797	self.failUnless(len(domain.boundary) == 4, msg)
1798
1799	# FIXME change to use get_xllcorner
1800	msg = 'Bad geo_reference'
1801	self.failUnless(domain.geo_reference.xllcorner == 140.0, msg)
1802
1803	os.remove(fileName)
1804
1805	def test_get_lone_vertices(self):
1806	a = [0.0, 0.0]
1807	b = [0.0, 2.0]
1808	c = [2.0, 0.0]
1809	d = [0.0, 4.0]
1810	e = [2.0, 2.0]
1811	f = [4.0, 0.0]
1812
1813	points = [a, b, c, d, e, f]
1814	# bac, bce, ecf, dbe
1815	vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1816	boundary = {(0, 0): 'Third',
1817	(0, 2): 'First',
1818	(2, 0): 'Second',
1819	(2, 1): 'Second',
1820	(3, 1): 'Second',
1821	(3, 2): 'Third'}
1822
1823	domain = Domain(points, vertices, boundary)
1824	domain.get_lone_vertices()
1825
1826	def test_fitting_using_shallow_water_domain(self):
1827	#Mesh in zone 56 (absolute coords)
1828
1829	x0 = 314036.58727982
1830	y0 = 6224951.2960092
1831
1832	a = [x0+0.0, y0+0.0]
1833	b = [x0+0.0, y0+2.0]
1834	c = [x0+2.0, y0+0.0]
1835	d = [x0+0.0, y0+4.0]
1836	e = [x0+2.0, y0+2.0]
1837	f = [x0+4.0, y0+0.0]
1838
1839	points = [a, b, c, d, e, f]
1840
1841	# bac, bce, ecf, dbe
1842	elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1843
1844	# absolute going in ..
1845	mesh4 = Domain(points, elements, geo_reference=Geo_reference(56, 0, 0))
1846	mesh4.check_integrity()
1847	quantity = Quantity(mesh4)
1848
1849	# Get (enough) datapoints (relative to georef)
1850	data_points_rel = [[ 0.66666667, 0.66666667],
1851	[ 1.33333333, 1.33333333],
1852	[ 2.66666667, 0.66666667],
1853	[ 0.66666667, 2.66666667],
1854	[ 0.0, 1.0],
1855	[ 0.0, 3.0],
1856	[ 1.0, 0.0],
1857	[ 1.0, 1.0],
1858	[ 1.0, 2.0],
1859	[ 1.0, 3.0],
1860	[ 2.0, 1.0],
1861	[ 3.0, 0.0],
1862	[ 3.0, 1.0]]
1863
1864	data_geo_spatial = Geospatial_data(data_points_rel,
1865	geo_reference=Geo_reference(56,
1866	x0,
1867	y0))
1868	data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
1869	attributes = linear_function(data_points_absolute)
1870	att = 'spam_and_eggs'
1871
1872	# Create .txt file
1873	ptsfile = tempfile.mktemp(".txt")
1874	file = open(ptsfile, "w")
1875	file.write(" x,y," + att + " \n")
1876	for data_point, attribute in map(None, data_points_absolute, attributes):
1877	row = (str(data_point[0]) + ',' +
1878	str(data_point[1]) + ',' +
1879	str(attribute))
1880	file.write(row + "\n")
1881	file.close()
1882
1883	# Check that values can be set from file
1884	quantity.set_values(filename=ptsfile, attribute_name=att, alpha=0)
1885	answer = linear_function(quantity.domain.get_vertex_coordinates())
1886
1887	assert num.allclose(quantity.vertex_values.flat, answer)
1888
1889	# Check that values can be set from file using default attribute
1890	quantity.set_values(filename = ptsfile, alpha = 0)
1891	assert num.allclose(quantity.vertex_values.flat, answer)
1892
1893	# Cleanup
1894	import os
1895	os.remove(ptsfile)
1896
1897	def test_fitting_example_that_crashed(self):
1898	"""This unit test has been derived from a real world example
1899	(the Towradgi '98 rainstorm simulation).
1900
1901	It shows a condition where fitting as called from set_quantity crashes
1902	when ANUGA mesh is reused. The test passes in the case where a new mesh
1903	is created.
1904
1905	See ticket:314
1906	"""
1907
1908	verbose = False
1909
1910	from anuga.shallow_water import Domain
1911	from anuga.pmesh.mesh_interface import create_mesh_from_regions
1912	from anuga.geospatial_data.geospatial_data import Geospatial_data
1913
1914
1915	# Get path where this test is run
1916	path = get_pathname_from_package('anuga.shallow_water')
1917
1918
1919	#----------------------------------------------------------------------
1920	# Create domain
1921	#--------------------------------------------------------------------
1922	W = 303400
1923	N = 6195800
1924	E = 308640
1925	S = 6193120
1926	bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
1927
1928	offending_regions = []
1929
1930	# From culvert 8
1931	offending_regions.append([[307611.43896231, 6193631.6894806],
1932	[307600.11394969, 6193608.2855474],
1933	[307597.41349586, 6193609.59227963],
1934	[307608.73850848, 6193632.99621282]])
1935	offending_regions.append([[307633.69143231, 6193620.9216536],
1936	[307622.36641969, 6193597.5177204],
1937	[307625.06687352, 6193596.21098818],
1938	[307636.39188614, 6193619.61492137]])
1939
1940	# From culvert 9
1941	offending_regions.append([[306326.69660524, 6194818.62900522],
1942	[306324.67939476, 6194804.37099478],
1943	[306323.75856492, 6194804.50127295],
1944	[306325.7757754, 6194818.7592834]])
1945	offending_regions.append([[306365.57160524, 6194813.12900522],
1946	[306363.55439476, 6194798.87099478],
1947	[306364.4752246, 6194798.7407166],
1948	[306366.49243508, 6194812.99872705]])
1949
1950	# From culvert 10
1951	offending_regions.append([[306955.071019428608, 6194465.704096679576],
1952	[306951.616980571358, 6194457.295903320424],
1953	[306950.044491164153, 6194457.941873183474],
1954	[306953.498530021403, 6194466.350066542625]])
1955	offending_regions.append([[307002.540019428649, 6194446.204096679576],
1956	[306999.085980571399, 6194437.795903320424],
1957	[307000.658469978604, 6194437.149933457375],
1958	[307004.112508835853, 6194445.558126816526]])
1959
1960	interior_regions = []
1961	for polygon in offending_regions:
1962	interior_regions.append( [polygon, 100] )
1963
1964	meshname = os.path.join(path, 'offending_mesh.msh')
1965	create_mesh_from_regions(bounding_polygon,
1966	boundary_tags={'south': [0], 'east': [1],
1967	'north': [2], 'west': [3]},
1968	maximum_triangle_area=1000000,
1969	interior_regions=interior_regions,
1970	filename=meshname,
1971	use_cache=False,
1972	verbose=verbose)
1973
1974	domain = Domain(meshname, use_cache=False, verbose=verbose)
1975
1976	#----------------------------------------------------------------------
1977	# Fit data point to mesh
1978	#----------------------------------------------------------------------
1979
1980	points_file = os.path.join(path, 'offending_point.pts')
1981
1982	# Offending point
1983	G = Geospatial_data(data_points=[[306953.344, 6194461.5]],
1984	attributes=[1])
1985	G.export_points_file(points_file)
1986
1987	try:
1988	domain.set_quantity('elevation', filename=points_file,
1989	use_cache=False, verbose=verbose, alpha=0.01)
1990	except RuntimeError, e:
1991	msg = 'Test failed: %s' % str(e)
1992	raise Exception, msg
1993	# clean up in case raise fails
1994	os.remove(meshname)
1995	os.remove(points_file)
1996	else:
1997	os.remove(meshname)
1998	os.remove(points_file)
1999
2000	#finally:
2001	# Cleanup regardless
2002	#FIXME(Ole): Finally does not work like this in python2.4
2003	#FIXME(Ole): Reinstate this when Python2.4 is out of the way
2004	#FIXME(Ole): Python 2.6 apparently introduces something called 'with'
2005	#os.remove(meshname)
2006	#os.remove(points_file)
2007
2008
2009	def test_fitting_example_that_crashed_2(self):
2010	"""test_fitting_example_that_crashed_2
2011
2012	This unit test has been derived from a real world example
2013	(the JJKelly study, by Petar Milevski).
2014
2015	It shows a condition where set_quantity crashes due to AtA
2016	not being built properly
2017
2018	See ticket:314
2019	"""
2020
2021	verbose = False
2022
2023	from anuga.shallow_water import Domain
2024	from anuga.pmesh.mesh_interface import create_mesh_from_regions
2025	from anuga.geospatial_data import Geospatial_data
2026
2027	# Get path where this test is run
2028	path = get_pathname_from_package('anuga.shallow_water')
2029
2030	meshname = os.path.join(path, 'test_mesh.msh')
2031	W = 304180
2032	S = 6185270
2033	E = 307650
2034	N = 6189040
2035	maximum_triangle_area = 1000000
2036
2037	bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
2038
2039	create_mesh_from_regions(bounding_polygon,
2040	boundary_tags={'south': [0],
2041	'east': [1],
2042	'north': [2],
2043	'west': [3]},
2044	maximum_triangle_area=maximum_triangle_area,
2045	filename=meshname,
2046	use_cache=False,
2047	verbose=verbose)
2048
2049	domain = Domain(meshname, use_cache=True, verbose=verbose)
2050
2051	# Large test set revealed one problem
2052	points_file = os.path.join(path, 'test_points_large.csv')
2053	try:
2054	domain.set_quantity('elevation', filename=points_file,
2055	use_cache=False, verbose=verbose)
2056	except AssertionError, e:
2057	msg = 'Test failed: %s' % str(e)
2058	raise Exception, msg
2059	# Cleanup in case this failed
2060	os.remove(meshname)
2061
2062	# Small test set revealed another problem
2063	points_file = os.path.join(path, 'test_points_small.csv')
2064	try:
2065	domain.set_quantity('elevation', filename=points_file,
2066	use_cache=False, verbose=verbose)
2067	except AssertionError, e:
2068	msg = 'Test failed: %s' % str(e)
2069	raise Exception, msg
2070	# Cleanup in case this failed
2071	os.remove(meshname)
2072	else:
2073	os.remove(meshname)
2074
2075
2076
2077
2078	def test_variable_elevation(self):
2079	"""test_variable_elevation
2080
2081	This will test that elevagtion van be stored in sww files
2082	as a time dependent quantity.
2083
2084	It will also chck that storage of other quantities
2085	can be controlled this way.
2086	"""
2087
2088	#---------------------------------------------------------------------
2089	# Import necessary modules
2090	#---------------------------------------------------------------------
2091	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
2092	from anuga.shallow_water import Domain
2093	from anuga.shallow_water import Reflective_boundary
2094	from anuga.shallow_water import Dirichlet_boundary
2095	from anuga.shallow_water import Time_boundary
2096
2097	#---------------------------------------------------------------------
2098	# Setup computational domain
2099	#---------------------------------------------------------------------
2100	length = 8.
2101	width = 6.
2102	dx = dy = 1 # Resolution: Length of subdivisions on both axes
2103
2104	inc = 0.05 # Elevation increment
2105
2106	points, vertices, boundary = rectangular_cross(int(length/dx),
2107	int(width/dy),
2108	len1=length,
2109	len2=width)
2110	domain = Domain(points, vertices, boundary)
2111	domain.set_name('channel_variable_test') # Output name
2112	domain.set_quantities_to_be_stored({'elevation': 2,
2113	'stage': 2})
2114
2115	#---------------------------------------------------------------------
2116	# Setup initial conditions
2117	#---------------------------------------------------------------------
2118
2119	def pole_increment(x,y):
2120	"""This provides a small increment to a pole located mid stream
2121	For use with variable elevation data
2122	"""
2123
2124	z = 0.0*x
2125
2126	N = len(x)
2127	for i in range(N):
2128	# Pole
2129	if (x[i] - 4)2 + (y[i] - 2)2 < 1.0**2:
2130	z[i] += inc
2131	return z
2132
2133	domain.set_quantity('elevation', 0.0) # Flat bed initially
2134	domain.set_quantity('friction', 0.01) # Constant friction
2135	domain.set_quantity('stage', 0.0) # Dry initial condition
2136
2137	#------------------------------------------------------------------
2138	# Setup boundary conditions
2139	#------------------------------------------------------------------
2140	Bi = Dirichlet_boundary([0.4, 0, 0]) # Inflow
2141	Br = Reflective_boundary(domain) # Solid reflective wall
2142	Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow
2143
2144	domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
2145
2146	#-------------------------------------------------------------------
2147	# Evolve system through time
2148	#-------------------------------------------------------------------
2149
2150	for t in domain.evolve(yieldstep=1, finaltime=3.0):
2151	#print domain.timestepping_statistics()
2152
2153	domain.add_quantity('elevation', pole_increment)
2154
2155
2156	# Check that quantities have been stored correctly
2157	from Scientific.IO.NetCDF import NetCDFFile
2158	sww_file = domain.get_name() + '.sww'
2159	fid = NetCDFFile(sww_file)
2160
2161	x = fid.variables['x'][:]
2162	y = fid.variables['y'][:]
2163	stage = fid.variables['stage'][:]
2164	elevation = fid.variables['elevation'][:]
2165	fid.close()
2166
2167	os.remove(sww_file)
2168
2169
2170	assert len(stage.shape) == 2
2171	assert len(elevation.shape) == 2
2172
2173	M, N = stage.shape
2174
2175	for i in range(M):
2176	# For each timestep
2177	assert num.allclose(max(elevation[i,:]), i * inc)
2178
2179
2180
2181	#################################################################################
2182
2183	if __name__ == "__main__":
2184	suite = unittest.makeSuite(Test_swb_basic, 'test')
2185	runner = unittest.TextTestRunner(verbosity=1)
2186	runner.run(suite)

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

Download in other formats: