Context Navigation

test_shallow_water.py @ 2950

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

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

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

Context Navigation

source: inundation/pyvolution/test_shallow_water.py @ 2950

Download in other formats: