Context Navigation

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

Last change on this file since 3632 was 3632, checked in by ole, 18 years ago
Changed gravity back to 9.8 as unit tests in shallow_water depend on it. Change this later.
File size: 117.4 KB

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

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

Download in other formats: