Context Navigation

source: inundation/ga/storm_surge/pyvolution/test_shallow_water.py @ 474

Last change on this file since 474 was 443, checked in by ole, 21 years ago
Work towards ensuring that this version (pyvolution mark 3) works as well as previous incarnation (pyvolution mark 2). This is not true yet as some examples e.g. netherlands display 'blow-up'.
File size: 71.6 KB

Rev	Line
[218]	1	#!/usr/bin/env python
	2
	3	import unittest
	4	from math import sqrt, pi
	5
	6	from config import g, epsilon
	7	from Numeric import allclose, array, zeros, ones, Float
	8	from shallow_water import *
	9
	10
	11	class TestCase(unittest.TestCase):
	12	def setUp(self):
	13	pass
	14
	15	def tearDown(self):
	16	pass
	17
[240]	18	def test_rotate(self):
	19	normal = array([0.0,-1.0])
[218]	20
[240]	21	q = array([1.0,2.0,3.0])
	22
	23	r = rotate(q, normal, direction = 1)
	24	assert r[0] == 1
	25	assert r[1] == -3
	26	assert r[2] == 2
	27
	28	w = rotate(r, normal, direction = -1)
	29	assert allclose(w, q)
	30
	31
	32
[218]	33	def test_flux_zero_case(self):
	34	ql = zeros( 3, Float )
	35	qr = zeros( 3, Float )
	36	normal = zeros( 2, Float )
	37	zl = zr = 0.
	38	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	39
	40	assert allclose(flux, [0,0,0])
	41	assert max_speed == 0.
	42
	43	def test_flux_constants(self):
	44	w = 2.0
	45
	46	normal = array([1.,0])
	47	ql = array([w, 0, 0])
	48	qr = array([w, 0, 0])
	49	zl = zr = 0.
	50	h = w - (zl+zr)/2
	51
	52	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	53
	54	assert allclose(flux, [0., 0.5gh**2, 0.])
	55	assert max_speed == sqrt(g*h)
	56
	57
	58	def test_flux1(self):
	59	#Use data from previous version of pyvolution
	60	normal = array([1.,0])
	61	ql = array([-0.2, 2, 3])
	62	qr = array([-0.2, 2, 3])
	63	zl = zr = -0.5
	64	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	65
	66	assert allclose(flux, [2.,13.77433333, 20.])
	67	assert allclose(max_speed, 8.38130948661)
	68
	69
	70	def test_flux2(self):
	71	#Use data from previous version of pyvolution
	72	normal = array([0., -1.])
	73	ql = array([-0.075, 2, 3])
	74	qr = array([-0.075, 2, 3])
	75	zl = zr = -0.375
	76	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	77
	78	assert allclose(flux, [-3.,-20.0, -30.441])
	79	assert allclose(max_speed, 11.7146428199)
	80
	81	def test_flux3(self):
	82	#Use data from previous version of pyvolution
	83	normal = array([-sqrt(2)/2, sqrt(2)/2])
	84	ql = array([-0.075, 2, 3])
	85	qr = array([-0.075, 2, 3])
	86	zl = zr = -0.375
	87	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	88
	89	assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019])
	90	assert allclose(max_speed, 4.0716654239)
	91
	92	def test_flux4(self):
	93	#Use data from previous version of pyvolution
	94	normal = array([-sqrt(2)/2, sqrt(2)/2])
	95	ql = array([-0.34319278, 0.10254161, 0.07273855])
	96	qr = array([-0.30683287, 0.1071986, 0.05930515])
	97	zl = zr = -0.375
	98	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	99
	100	assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364])
	101	assert allclose(max_speed, 1.31414103233)
	102
	103	def test_sw_domain_simple(self):
	104	a = [0.0, 0.0]
	105	b = [0.0, 2.0]
	106	c = [2.0,0.0]
	107	d = [0.0, 4.0]
	108	e = [2.0, 2.0]
	109	f = [4.0,0.0]
	110
	111	points = [a, b, c, d, e, f]
	112	#bac, bce, ecf, dbe, daf, dae
	113	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
	114
	115	domain = Domain(points, vertices)
	116	domain.check_integrity()
	117
	118	for name in ['level', 'xmomentum', 'ymomentum',
	119	'elevation', 'friction']:
	120	assert domain.quantities.has_key(name)
	121
	122
	123	assert domain.get_conserved_quantities(0, edge=1) == 0.
	124
	125
	126	def test_boundary_conditions(self):
	127
	128	a = [0.0, 0.0]
	129	b = [0.0, 2.0]
	130	c = [2.0,0.0]
	131	d = [0.0, 4.0]
	132	e = [2.0, 2.0]
	133	f = [4.0,0.0]
	134
	135	points = [a, b, c, d, e, f]
	136	#bac, bce, ecf, dbe
	137	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
	138	boundary = { (0, 0): 'Third',
	139	(0, 2): 'First',
	140	(2, 0): 'Second',
	141	(2, 1): 'Second',
	142	(3, 1): 'Second',
	143	(3, 2): 'Third'}
	144
	145
	146	domain = Domain(points, vertices, boundary)
	147	domain.check_integrity()
	148
	149
	150	domain.set_quantity('level', [[1,2,3], [5,5,5],
	151	[0,0,9], [-6, 3, 3]])
	152
	153	domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
	154	[3,3,3], [4, 4, 4]])
	155
	156	domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
	157	[30,30,30], [40, 40, 40]])
	158
	159
	160	D = Dirichlet_boundary([5,2,1])
	161	T = Transmissive_boundary(domain)
	162	R = Reflective_boundary(domain)
	163	domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
	164
	165	domain.update_boundary()
	166
	167	#Level
	168	assert domain.quantities['level'].boundary_values[0] == 2.5
	169	assert domain.quantities['level'].boundary_values[0] ==\
	170	domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
	171	assert domain.quantities['level'].boundary_values[1] == 5. #Dirichlet
	172	assert domain.quantities['level'].boundary_values[2] ==\
	173	domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
	174	assert domain.quantities['level'].boundary_values[3] ==\
	175	domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
	176	assert domain.quantities['level'].boundary_values[4] ==\
	177	domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
	178	assert domain.quantities['level'].boundary_values[5] ==\
	179	domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
	180
	181	#Xmomentum
	182	assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
	183	assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
	184	assert domain.quantities['xmomentum'].boundary_values[2] ==\
	185	domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
	186	assert domain.quantities['xmomentum'].boundary_values[3] ==\
	187	domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
	188	assert domain.quantities['xmomentum'].boundary_values[4] ==\
	189	domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
	190	assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective
	191
	192	#Ymomentum
	193	assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
	194	assert domain.quantities['ymomentum'].boundary_values[1] == 1. #Dirichlet
	195	assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
	196	assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
	197	assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
	198	assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
	199
	200
	201
	202
[240]	203	def test_compute_fluxes0(self):
[218]	204	#Do a full triangle and check that fluxes cancel out for
	205	#the constant level case
	206
	207	a = [0.0, 0.0]
	208	b = [0.0, 2.0]
	209	c = [2.0,0.0]
	210	d = [0.0, 4.0]
	211	e = [2.0, 2.0]
	212	f = [4.0,0.0]
	213
	214	points = [a, b, c, d, e, f]
	215	#bac, bce, ecf, dbe
	216	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	217
	218	domain = Domain(points, vertices)
	219	domain.set_quantity('level', [[2,2,2], [2,2,2],
	220	[2,2,2], [2,2,2]])
	221	domain.check_integrity()
	222
	223	zl=zr=0. #Assume flat bed
	224
	225	#Flux across right edge of volume 1
	226	normal = domain.get_normal(1,0)
	227	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
	228	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
	229	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
	230
	231	#Check that flux seen from other triangles is inverse
	232	tmp = qr; qr=ql; ql=tmp
	233	normal = domain.get_normal(2,2)
	234	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	235	assert allclose(flux + flux0, 0.)
	236
	237	#Flux across upper edge of volume 1
	238	normal = domain.get_normal(1,1)
	239	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
	240	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
	241	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
	242
	243	#Check that flux seen from other triangles is inverse
	244	tmp = qr; qr=ql; ql=tmp
	245	normal = domain.get_normal(3,0)
	246	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	247	assert allclose(flux + flux1, 0.)
	248
	249	#Flux across lower left hypotenuse of volume 1
	250	normal = domain.get_normal(1,2)
	251	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
	252	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
	253	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
	254
	255	#Check that flux seen from other triangles is inverse
	256	tmp = qr; qr=ql; ql=tmp
	257	normal = domain.get_normal(0,1)
	258	flux, max_speed = flux_function(normal, ql, qr, zl, zr)
	259	assert allclose(flux + flux2, 0.)
	260
	261
	262	#Scale by edgelengths, add up anc check that total flux is zero
	263	e0 = domain.edgelengths[1, 0]
	264	e1 = domain.edgelengths[1, 1]
	265	e2 = domain.edgelengths[1, 2]
	266
	267	assert allclose(e0flux0+e1flux1+e2*flux2, 0.)
	268
	269	#Now check that compute_flux yields zeros as well
	270	domain.compute_fluxes()
	271
	272	for name in ['level', 'xmomentum', 'ymomentum']:
	273	#print name, domain.quantities[name].explicit_update
	274	assert allclose(domain.quantities[name].explicit_update[1], 0)
	275
	276
	277
	278	def test_compute_fluxes1(self):
	279	#Use values from previous version
	280
	281	a = [0.0, 0.0]
	282	b = [0.0, 2.0]
	283	c = [2.0,0.0]
	284	d = [0.0, 4.0]
	285	e = [2.0, 2.0]
	286	f = [4.0,0.0]
	287
	288	points = [a, b, c, d, e, f]
	289	#bac, bce, ecf, dbe
	290	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	291
	292	domain = Domain(points, vertices)
	293	val0 = 2.+2.0/3
	294	val1 = 4.+4.0/3
	295	val2 = 8.+2.0/3
	296	val3 = 2.+8.0/3
	297
	298	domain.set_quantity('level', [[val0, val0, val0], [val1, val1, val1],
	299	[val2, val2, val2], [val3, val3, val3]])
	300	domain.check_integrity()
	301
	302	zl=zr=0. #Assume flat bed
	303
	304	#Flux across right edge of volume 1
	305	normal = domain.get_normal(1,0)
	306	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
	307	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
	308	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
	309	assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
	310	assert allclose(max_speed, 9.21592824046)
	311
	312
	313	#Flux across upper edge of volume 1
	314	normal = domain.get_normal(1,1)
	315	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
	316	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
	317	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
	318	assert allclose(flux1, [2.4098563, 0., 123.04444444])
	319	assert allclose(max_speed, 7.22956891292)
	320
	321	#Flux across lower left hypotenuse of volume 1
	322	normal = domain.get_normal(1,2)
	323	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
	324	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
	325	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
	326
	327	assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
	328	assert allclose(max_speed, 7.22956891292)
	329
	330	#Scale, add up and check that compute_fluxes is correct for vol 1
	331	e0 = domain.edgelengths[1, 0]
	332	e1 = domain.edgelengths[1, 1]
	333	e2 = domain.edgelengths[1, 2]
	334
	335	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
	336	assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
	337
	338
	339	domain.compute_fluxes()
	340
	341	#assert allclose(total_flux, domain.explicit_update[1,:])
	342	for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
	343	assert allclose(total_flux[i],
	344	domain.quantities[name].explicit_update[1])
	345
	346	#assert allclose(domain.explicit_update, [
	347	# [0., -69.68888889, -69.68888889],
	348	# [-0.68218178, -166.6, -35.93333333],
	349	# [-111.77316251, 69.68888889, 0.],
	350	# [-35.68522449, 0., 69.68888889]])
	351
	352	assert allclose(domain.quantities['level'].explicit_update,
	353	[0., -0.68218178, -111.77316251, -35.68522449])
	354	assert allclose(domain.quantities['xmomentum'].explicit_update,
	355	[-69.68888889, -166.6, 69.68888889, 0])
	356	assert allclose(domain.quantities['ymomentum'].explicit_update,
	357	[-69.68888889, -35.93333333, 0., 69.68888889])
	358
	359
	360	#assert allclose(domain.quantities[name].explicit_update
	361
	362
	363
	364
	365
	366	def test_compute_fluxes2(self):
	367	#Random values, incl momentum
	368
	369	a = [0.0, 0.0]
	370	b = [0.0, 2.0]
	371	c = [2.0,0.0]
	372	d = [0.0, 4.0]
	373	e = [2.0, 2.0]
	374	f = [4.0,0.0]
	375
	376	points = [a, b, c, d, e, f]
	377	#bac, bce, ecf, dbe
	378	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	379
	380	domain = Domain(points, vertices)
	381	val0 = 2.+2.0/3
	382	val1 = 4.+4.0/3
	383	val2 = 8.+2.0/3
	384	val3 = 2.+8.0/3
	385
	386	zl=zr=0 #Assume flat zero bed
	387
	388	domain.set_quantity('elevation', zl*ones( (4,3) ))
	389
	390
	391	domain.set_quantity('level', [[val0, val0-1, val0-2],
	392	[val1, val1+1, val1],
	393	[val2, val2-2, val2],
	394	[val3-0.5, val3, val3]])
	395
	396	domain.set_quantity('xmomentum',
	397	[[1, 2, 3], [3, 4, 5],
	398	[1, -1, 0], [0, -2, 2]])
	399
	400	domain.set_quantity('ymomentum',
	401	[[1, -1, 0], [0, -3, 2],
	402	[0, 1, 0], [-1, 2, 2]])
	403
	404
	405	domain.check_integrity()
	406
	407
	408
	409	#Flux across right edge of volume 1
	410	normal = domain.get_normal(1,0)
	411	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
	412	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
	413	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
	414
	415	#Flux across upper edge of volume 1
	416	normal = domain.get_normal(1,1)
	417	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
	418	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
	419	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
	420
	421	#Flux across lower left hypotenuse of volume 1
	422	normal = domain.get_normal(1,2)
	423	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
	424	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
	425	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
	426
	427	#Scale, add up and check that compute_fluxes is correct for vol 1
	428	e0 = domain.edgelengths[1, 0]
	429	e1 = domain.edgelengths[1, 1]
	430	e2 = domain.edgelengths[1, 2]
	431
	432	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
	433
	434
	435	domain.compute_fluxes()
	436	for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
	437	assert allclose(total_flux[i],
	438	domain.quantities[name].explicit_update[1])
	439	#assert allclose(total_flux, domain.explicit_update[1,:])
	440
	441
	442	def test_compute_fluxes3(self):
	443	#Random values, incl momentum
	444
	445	a = [0.0, 0.0]
	446	b = [0.0, 2.0]
	447	c = [2.0,0.0]
	448	d = [0.0, 4.0]
	449	e = [2.0, 2.0]
	450	f = [4.0,0.0]
	451
	452	points = [a, b, c, d, e, f]
	453	#bac, bce, ecf, dbe
	454	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	455
	456	domain = Domain(points, vertices)
	457	val0 = 2.+2.0/3
	458	val1 = 4.+4.0/3
	459	val2 = 8.+2.0/3
	460	val3 = 2.+8.0/3
	461
	462	zl=zr=-3.75 #Assume constant bed (must be less than level)
	463	domain.set_quantity('elevation', zl*ones( (4,3) ))
	464
	465
	466	domain.set_quantity('level', [[val0, val0-1, val0-2],
	467	[val1, val1+1, val1],
	468	[val2, val2-2, val2],
	469	[val3-0.5, val3, val3]])
	470
	471	domain.set_quantity('xmomentum',
	472	[[1, 2, 3], [3, 4, 5],
	473	[1, -1, 0], [0, -2, 2]])
	474
	475	domain.set_quantity('ymomentum',
	476	[[1, -1, 0], [0, -3, 2],
	477	[0, 1, 0], [-1, 2, 2]])
	478
	479
	480	domain.check_integrity()
	481
	482
	483
	484	#Flux across right edge of volume 1
	485	normal = domain.get_normal(1,0)
	486	ql = domain.get_conserved_quantities(vol_id=1, edge=0)
	487	qr = domain.get_conserved_quantities(vol_id=2, edge=2)
	488	flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
	489
	490	#Flux across upper edge of volume 1
	491	normal = domain.get_normal(1,1)
	492	ql = domain.get_conserved_quantities(vol_id=1, edge=1)
	493	qr = domain.get_conserved_quantities(vol_id=3, edge=0)
	494	flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
	495
	496	#Flux across lower left hypotenuse of volume 1
	497	normal = domain.get_normal(1,2)
	498	ql = domain.get_conserved_quantities(vol_id=1, edge=2)
	499	qr = domain.get_conserved_quantities(vol_id=0, edge=1)
	500	flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
	501
	502	#Scale, add up and check that compute_fluxes is correct for vol 1
	503	e0 = domain.edgelengths[1, 0]
	504	e1 = domain.edgelengths[1, 1]
	505	e2 = domain.edgelengths[1, 2]
	506
	507	total_flux = -(e0flux0+e1flux1+e2*flux2)/domain.areas[1]
	508
	509	domain.compute_fluxes()
	510	for i, name in enumerate(['level', 'xmomentum', 'ymomentum']):
	511	assert allclose(total_flux[i],
	512	domain.quantities[name].explicit_update[1])
	513
	514
	515
	516	def test_catching_negative_heights(self):
	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=4 #Too large
	536	domain.set_quantity('elevation', zl*ones( (4,3) ))
	537	domain.set_quantity('level', [[val0, val0-1, val0-2],
	538	[val1, val1+1, val1],
	539	[val2, val2-2, val2],
	540	[val3-0.5, val3, val3]])
	541
	542	#Should fail
	543	try:
	544	domain.check_integrity()
	545	except:
	546	pass
	547
	548	#####################################################
	549	def test_gravity(self):
	550	#Assuming no friction
	551
	552	from config import g
	553
	554	a = [0.0, 0.0]
	555	b = [0.0, 2.0]
	556	c = [2.0, 0.0]
	557	d = [0.0, 4.0]
	558	e = [2.0, 2.0]
	559	f = [4.0, 0.0]
	560
	561	points = [a, b, c, d, e, f]
	562	#bac, bce, ecf, dbe
	563	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	564
	565	domain = Domain(points, vertices)
	566
	567	#Set up for a gradient of (3,0) at mid triangle
	568	def slope(x, y):
	569	return 3*x
	570
	571	h = 0.1
	572	def level(x,y):
	573	return slope(x,y)+h
	574
	575	domain.set_quantity('elevation', slope)
	576	domain.set_quantity('level', level)
	577
	578	for name in domain.conserved_quantities:
	579	assert allclose(domain.quantities[name].explicit_update, 0)
	580	assert allclose(domain.quantities[name].semi_implicit_update, 0)
	581
	582	domain.compute_forcing_terms()
	583
	584	assert allclose(domain.quantities['level'].explicit_update, 0)
	585	assert allclose(domain.quantities['xmomentum'].explicit_update, -gh3)
	586	assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
	587
	588
	589
	590	#####################################################
	591	def test_first_order_extrapolator_const_z(self):
	592
	593	a = [0.0, 0.0]
	594	b = [0.0, 2.0]
	595	c = [2.0, 0.0]
	596	d = [0.0, 4.0]
	597	e = [2.0, 2.0]
	598	f = [4.0, 0.0]
	599
	600	points = [a, b, c, d, e, f]
	601	#bac, bce, ecf, dbe
	602	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	603
	604	domain = Domain(points, vertices)
	605	val0 = 2.+2.0/3
	606	val1 = 4.+4.0/3
	607	val2 = 8.+2.0/3
	608	val3 = 2.+8.0/3
	609
	610	zl=zr=-3.75 #Assume constant bed (must be less than level)
	611	domain.set_quantity('elevation', zl*ones( (4,3) ))
	612	domain.set_quantity('level', [[val0, val0-1, val0-2],
	613	[val1, val1+1, val1],
	614	[val2, val2-2, val2],
	615	[val3-0.5, val3, val3]])
	616
	617
	618
	619	domain.order = 1
	620	domain.distribute_to_vertices_and_edges()
	621
	622	#Check that centroid values were distributed to vertices
	623	C = domain.quantities['level'].centroid_values
	624	for i in range(3):
	625	assert allclose( domain.quantities['level'].vertex_values[:,i], C)
	626
	627
	628	def test_first_order_limiter_variable_z(self):
	629	#Check that first order limiter follows bed_slope
	630	from Numeric import alltrue, greater_equal
	631	from config import epsilon
	632
	633	a = [0.0, 0.0]
	634	b = [0.0, 2.0]
	635	c = [2.0,0.0]
	636	d = [0.0, 4.0]
	637	e = [2.0, 2.0]
	638	f = [4.0,0.0]
	639
	640	points = [a, b, c, d, e, f]
	641	#bac, bce, ecf, dbe
	642	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	643
	644	domain = Domain(points, vertices)
	645	val0 = 2.+2.0/3
	646	val1 = 4.+4.0/3
	647	val2 = 8.+2.0/3
	648	val3 = 2.+8.0/3
	649
	650	domain.set_quantity('elevation', [[0,0,0], [6,0,0],
	651	[6,6,6], [6,6,6]])
	652	domain.set_quantity('level', [[val0, val0, val0],
	653	[val1, val1, val1],
	654	[val2, val2, val2],
	655	[val3, val3, val3]])
	656
	657	E = domain.quantities['elevation'].vertex_values
	658	L = domain.quantities['level'].vertex_values
	659
	660
	661	#Check that some levels are not above elevation (within eps)
	662	#- so that the limiter has something to work with
	663	assert not alltrue(alltrue(greater_equal(L,E-epsilon)))
	664
	665	domain.order = 1
	666	domain.distribute_to_vertices_and_edges()
	667
	668	#Check that all levels are above elevation (within eps)
	669	assert alltrue(alltrue(greater_equal(L,E-epsilon)))
	670
	671
	672	#####################################################
	673	def test_distribute_basic(self):
	674	#Using test data generated by pyvolution-2
	675	#Assuming no friction and flat bed (0.0)
	676
	677	a = [0.0, 0.0]
	678	b = [0.0, 2.0]
	679	c = [2.0, 0.0]
	680	d = [0.0, 4.0]
	681	e = [2.0, 2.0]
	682	f = [4.0, 0.0]
	683
	684	points = [a, b, c, d, e, f]
	685	#bac, bce, ecf, dbe
	686	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	687
	688	domain = Domain(points, vertices)
	689
	690	val0 = 2.
	691	val1 = 4.
	692	val2 = 8.
	693	val3 = 2.
	694
	695	domain.set_quantity('level', [val0, val1, val2, val3], 'centroids')
	696	L = domain.quantities['level'].vertex_values
	697
	698	#First order
	699	domain.order = 1
	700	domain.distribute_to_vertices_and_edges()
	701	assert allclose(L[1], val1)
	702
	703	#Second order
	704	domain.order = 2
	705	domain.distribute_to_vertices_and_edges()
	706	assert allclose(L[1], [2.2, 4.9, 4.9])
	707
	708
	709
	710	def test_distribute_away_from_bed(self):
	711	#Using test data generated by pyvolution-2
	712	#Assuming no friction and flat bed (0.0)
	713
	714	a = [0.0, 0.0]
	715	b = [0.0, 2.0]
	716	c = [2.0, 0.0]
	717	d = [0.0, 4.0]
	718	e = [2.0, 2.0]
	719	f = [4.0, 0.0]
	720
	721	points = [a, b, c, d, e, f]
	722	#bac, bce, ecf, dbe
	723	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	724
	725	domain = Domain(points, vertices)
	726	L = domain.quantities['level'].vertex_values
	727
	728	def level(x,y):
	729	return x**2
	730
	731	domain.set_quantity('level', level, 'centroids')
	732
	733	a, b = domain.quantities['level'].compute_gradients()
	734	assert allclose(a[1], 3.33333334)
	735	assert allclose(b[1], 0.0)
	736
	737	domain.order = 1
	738	domain.distribute_to_vertices_and_edges()
	739	assert allclose(L[1], 1.77777778)
	740
	741	domain.order = 2
	742	domain.distribute_to_vertices_and_edges()
	743	assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
	744
	745
	746
	747	def test_distribute_away_from_bed1(self):
	748	#Using test data generated by pyvolution-2
	749	#Assuming no friction and flat bed (0.0)
	750
	751	a = [0.0, 0.0]
	752	b = [0.0, 2.0]
	753	c = [2.0, 0.0]
	754	d = [0.0, 4.0]
	755	e = [2.0, 2.0]
	756	f = [4.0, 0.0]
	757
	758	points = [a, b, c, d, e, f]
	759	#bac, bce, ecf, dbe
	760	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	761
	762	domain = Domain(points, vertices)
	763	L = domain.quantities['level'].vertex_values
	764
	765	def level(x,y):
	766	return x4+y2
	767
	768	domain.set_quantity('level', level, 'centroids')
	769	#print domain.quantities['level'].centroid_values
	770
	771	a, b = domain.quantities['level'].compute_gradients()
	772	assert allclose(a[1], 25.18518519)
	773	assert allclose(b[1], 3.33333333)
	774
	775	domain.order = 1
	776	domain.distribute_to_vertices_and_edges()
	777	assert allclose(L[1], 4.9382716)
	778
	779	domain.order = 2
	780	domain.distribute_to_vertices_and_edges()
	781	assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
	782
	783
	784
	785	def test_distribute_near_bed(self):
	786	#Using test data generated by pyvolution-2
	787	#Assuming no friction and flat bed (0.0)
	788
	789	a = [0.0, 0.0]
	790	b = [0.0, 2.0]
	791	c = [2.0, 0.0]
	792	d = [0.0, 4.0]
	793	e = [2.0, 2.0]
	794	f = [4.0, 0.0]
	795
	796	points = [a, b, c, d, e, f]
	797	#bac, bce, ecf, dbe
	798	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	799
	800	domain = Domain(points, vertices)
	801
	802
	803	#Set up for a gradient of (3,0) at mid triangle
	804	def slope(x, y):
	805	return 10*x
	806
	807	h = 0.1
	808	def level(x,y):
	809	return slope(x,y)+h
	810
	811	domain.set_quantity('elevation', slope)
	812	domain.set_quantity('level', level, 'centroids')
	813
	814	#print domain.quantities['elevation'].centroid_values
	815	#print domain.quantities['level'].centroid_values
	816
	817	E = domain.quantities['elevation'].vertex_values
	818	L = domain.quantities['level'].vertex_values
	819
	820	#print E
	821	domain.order = 1
	822	domain.distribute_to_vertices_and_edges()
[272]	823	##assert allclose(L[1], [0.19999999, 20.05, 20.05])
	824	assert allclose(L[1], [0.1, 20.1, 20.1])
	825
[218]	826	domain.order = 2
	827	domain.distribute_to_vertices_and_edges()
	828	assert allclose(L[1], [0.1, 20.1, 20.1])
	829
	830	def test_distribute_near_bed1(self):
	831	#Using test data generated by pyvolution-2
	832	#Assuming no friction and flat bed (0.0)
	833
	834	a = [0.0, 0.0]
	835	b = [0.0, 2.0]
	836	c = [2.0, 0.0]
	837	d = [0.0, 4.0]
	838	e = [2.0, 2.0]
	839	f = [4.0, 0.0]
	840
	841	points = [a, b, c, d, e, f]
	842	#bac, bce, ecf, dbe
	843	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	844
	845	domain = Domain(points, vertices)
	846
	847
	848	#Set up for a gradient of (3,0) at mid triangle
	849	def slope(x, y):
	850	return x4+y2
	851
	852	h = 0.1
	853	def level(x,y):
	854	return slope(x,y)+h
	855
	856	domain.set_quantity('elevation', slope)
	857	domain.set_quantity('level', level)
	858
	859	#print domain.quantities['elevation'].centroid_values
	860	#print domain.quantities['level'].centroid_values
	861
	862	E = domain.quantities['elevation'].vertex_values
	863	L = domain.quantities['level'].vertex_values
	864
	865	#print E
[272]	866	domain.order = 1
[218]	867	domain.distribute_to_vertices_and_edges()
[272]	868	##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
	869	assert allclose(L[1], [4.1, 16.1, 20.1])
[218]	870
	871	domain.order = 2
	872	domain.distribute_to_vertices_and_edges()
	873	assert allclose(L[1], [4.1, 16.1, 20.1])
	874
[229]	875
	876
	877	def test_second_order_distribute_real_data(self):
	878	#Using test data generated by pyvolution-2
	879	#Assuming no friction and flat bed (0.0)
	880
	881	a = [0.0, 0.0]
	882	b = [0.0, 1.0/5]
	883	c = [0.0, 2.0/5]
	884	d = [1.0/5, 0.0]
	885	e = [1.0/5, 1.0/5]
	886	f = [1.0/5, 2.0/5]
	887	g = [2.0/5, 2.0/5]
	888
	889	points = [a, b, c, d, e, f, g]
	890	#bae, efb, cbf, feg
	891	vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
[218]	892
[229]	893	domain = Domain(points, vertices)
	894
	895	def slope(x, y):
	896	return -x/3
	897
	898	domain.set_quantity('elevation', slope)
	899	domain.set_quantity('level',
	900	[0.01298164, 0.00365611, 0.01440365, -0.0381856437096],
	901	'centroids')
	902	domain.set_quantity('xmomentum',
	903	[0.00670439, 0.01263789, 0.00647805, 0.0178180740668],
	904	'centroids')
	905	domain.set_quantity('ymomentum',
	906	[-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866],
	907	'centroids')
	908
	909	E = domain.quantities['elevation'].vertex_values
	910	L = domain.quantities['level'].vertex_values
	911	X = domain.quantities['xmomentum'].vertex_values
	912	Y = domain.quantities['ymomentum'].vertex_values
	913
	914	#print E
	915	domain.order = 2
	916	domain.distribute_to_vertices_and_edges()
	917
	918	#print L[1,:]
	919	#print X[1,:]
	920	#print Y[1,:]
	921
	922	assert allclose(L[1,:], [-0.00825735775384, -0.00801881482869, 0.0272445025825])
	923	assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
	924	assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
	925
	926
[218]	927	def test_second_order_flat_bed_onestep(self):
	928
	929	from mesh_factory import rectangular
	930	from shallow_water import Domain, Reflective_boundary,\
	931	Dirichlet_boundary, Constant_height
	932	from Numeric import array
	933
	934	#Create basic mesh
	935	points, vertices, boundary = rectangular(6, 6)
	936
	937	#Create shallow water domain
	938	domain = Domain(points, vertices, boundary)
	939	domain.smooth = False
	940	domain.default_order=2
	941
	942	# Boundary conditions
	943	Br = Reflective_boundary(domain)
	944	Bd = Dirichlet_boundary([0.1, 0., 0.])
	945	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	946
	947	domain.check_integrity()
	948
	949	#Evolution
	950	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
	951	pass# domain.write_time()
	952
	953	#Data from earlier version of pyvolution
	954	assert allclose(domain.min_timestep, 0.0396825396825)
	955	assert allclose(domain.max_timestep, 0.0396825396825)
	956
	957	assert allclose(domain.quantities['level'].centroid_values[:12],
	958	[0.00171396, 0.02656103, 0.00241523, 0.02656103,
	959	0.00241523, 0.02656103, 0.00241523, 0.02656103,
	960	0.00241523, 0.02656103, 0.00241523, 0.0272623])
	961
	962	domain.distribute_to_vertices_and_edges()
	963	assert allclose(domain.quantities['level'].vertex_values[:12,0],
	964	[0.0001714, 0.02656103, 0.00024152,
	965	0.02656103, 0.00024152, 0.02656103,
	966	0.00024152, 0.02656103, 0.00024152,
	967	0.02656103, 0.00024152, 0.0272623])
	968
	969	assert allclose(domain.quantities['level'].vertex_values[:12,1],
	970	[0.00315012, 0.02656103, 0.00024152, 0.02656103,
	971	0.00024152, 0.02656103, 0.00024152, 0.02656103,
	972	0.00024152, 0.02656103, 0.00040506, 0.0272623])
	973
	974	assert allclose(domain.quantities['level'].vertex_values[:12,2],
	975	[0.00182037, 0.02656103, 0.00676264,
	976	0.02656103, 0.00676264, 0.02656103,
	977	0.00676264, 0.02656103, 0.00676264,
	978	0.02656103, 0.0065991, 0.0272623])
	979
	980	assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
	981	[0.00113961, 0.01302432, 0.00148672,
	982	0.01302432, 0.00148672, 0.01302432,
	983	0.00148672, 0.01302432, 0.00148672 ,
	984	0.01302432, 0.00148672, 0.01337143])
	985
	986	assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
	987	[-2.91240050e-004 , 1.22721531e-004,
	988	-1.22721531e-004, 1.22721531e-004 ,
	989	-1.22721531e-004, 1.22721531e-004,
	990	-1.22721531e-004 , 1.22721531e-004,
	991	-1.22721531e-004, 1.22721531e-004,
	992	-1.22721531e-004, -4.57969873e-005])
	993
	994
	995
	996	def test_second_order_flat_bed_moresteps(self):
	997
	998	from mesh_factory import rectangular
	999	from shallow_water import Domain, Reflective_boundary,\
	1000	Dirichlet_boundary, Constant_height
	1001	from Numeric import array
	1002
	1003	#Create basic mesh
	1004	points, vertices, boundary = rectangular(6, 6)
	1005
	1006	#Create shallow water domain
	1007	domain = Domain(points, vertices, boundary)
	1008	domain.smooth = False
	1009	domain.default_order=2
	1010
	1011	# Boundary conditions
	1012	Br = Reflective_boundary(domain)
	1013	Bd = Dirichlet_boundary([0.1, 0., 0.])
	1014	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	1015
	1016	domain.check_integrity()
	1017
	1018	#Evolution
	1019	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
	1020	pass
	1021
	1022	#Data from earlier version of pyvolution
	1023	#assert allclose(domain.min_timestep, 0.0396825396825)
	1024	#assert allclose(domain.max_timestep, 0.0396825396825)
	1025	#print domain.quantities['level'].centroid_values
	1026
[221]	1027
	1028	def test_flatbed_first_order(self):
	1029	from mesh_factory import rectangular
	1030	from shallow_water import Domain,\
	1031	Reflective_boundary, Dirichlet_boundary
	1032
	1033	from Numeric import array
	1034
	1035	#Create basic mesh
	1036	N = 8
	1037	points, vertices, boundary = rectangular(N, N)
	1038
	1039	#Create shallow water domain
	1040	domain = Domain(points, vertices, boundary)
	1041	domain.smooth = False
	1042	domain.visualise = False
	1043	domain.default_order=1
	1044
	1045	# Boundary conditions
	1046	Br = Reflective_boundary(domain)
	1047	Bd = Dirichlet_boundary([0.2,0.,0.])
[218]	1048
[221]	1049	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	1050	domain.check_integrity()
	1051
	1052
	1053	#Evolution
	1054	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
	1055	pass
	1056	#domain.write_time()
	1057
[443]	1058	#FIXME: These numbers were from version before 25/10
	1059	#assert allclose(domain.min_timestep, 0.0140413643926)
	1060	#assert allclose(domain.max_timestep, 0.0140947355753)
[221]	1061
	1062	for i in range(3):
[443]	1063	#assert allclose(domain.quantities['level'].edge_values[:4,i],
	1064	# [0.10730244,0.12337617,0.11200126,0.12605666])
[221]	1065
	1066	assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
	1067	[0.07610894,0.06901572,0.07284461,0.06819712])
	1068
[443]	1069	#assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
	1070	# [-0.0060238, -0.00157404, -0.00309633, -0.0001637])
[218]	1071
[221]	1072	def test_flatbed_second_order(self):
	1073	from mesh_factory import rectangular
	1074	from shallow_water import Domain,\
	1075	Reflective_boundary, Dirichlet_boundary
	1076
	1077	from Numeric import array
	1078
	1079	#Create basic mesh
	1080	N = 8
	1081	points, vertices, boundary = rectangular(N, N)
	1082
	1083	#Create shallow water domain
	1084	domain = Domain(points, vertices, boundary)
	1085	domain.smooth = False
	1086	domain.visualise = False
	1087	domain.default_order=2
[272]	1088	#domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
[221]	1089
	1090	# Boundary conditions
	1091	Br = Reflective_boundary(domain)
	1092	Bd = Dirichlet_boundary([0.2,0.,0.])
	1093
	1094	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	1095	domain.check_integrity()
	1096
	1097	#Evolution
	1098	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
	1099	pass
	1100
	1101
	1102	assert allclose(domain.min_timestep, 0.0210448446782)
	1103	assert allclose(domain.max_timestep, 0.0210448446782)
	1104
	1105	#print domain.quantities['level'].vertex_values[:4,0]
	1106	#print domain.quantities['xmomentum'].vertex_values[:4,0]
	1107	#print domain.quantities['ymomentum'].vertex_values[:4,0]
[443]	1108
	1109	#FIXME: These numbers were from version before 25/10
	1110	#assert allclose(domain.quantities['level'].vertex_values[:4,0],
	1111	# [0.00101913,0.05352143,0.00104852,0.05354394])
[221]	1112
	1113	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
[272]	1114	[ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
	1115
	1116	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
	1117	# [0.00090581,0.03685719,0.00088303,0.03687313])
	1118
[221]	1119	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
	1120	[-0.00139463,0.0006156,-0.00060364,0.00061827])
	1121
	1122
	1123
	1124	def test_flatbed_second_order_distribute(self):
	1125	#Use real data from pyvolution 2
	1126	#painfully setup and extracted.
	1127	from mesh_factory import rectangular
	1128	from shallow_water import Domain,\
	1129	Reflective_boundary, Dirichlet_boundary
	1130
	1131	from Numeric import array
	1132
	1133	#Create basic mesh
	1134	N = 8
	1135	points, vertices, boundary = rectangular(N, N)
	1136
	1137	#Create shallow water domain
	1138	domain = Domain(points, vertices, boundary)
	1139	domain.smooth = False
	1140	domain.visualise = False
	1141	domain.default_order=domain.order=2
	1142
	1143	# Boundary conditions
	1144	Br = Reflective_boundary(domain)
	1145	Bd = Dirichlet_boundary([0.2,0.,0.])
	1146
	1147	domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
	1148	domain.check_integrity()
	1149
	1150
	1151
	1152	for V in [False, True]:
	1153	if V:
	1154	#Set centroids as if system had been evolved
	1155	L = zeros(2NN, Float)
	1156	L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
	1157	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
	1158	1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
	1159	5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
	1160	1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
	1161	5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
	1162	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
	1163	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
	1164	0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
	1165	4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
	1166	0.00000000e+000, 5.57305948e-005]
	1167
	1168	X = zeros(2NN, Float)
	1169	X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
	1170	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
	1171	8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
	1172	3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
	1173	8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
	1174	3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
	1175	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
	1176	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
	1177	0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
	1178	3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
	1179	0.00000000e+000, 4.57662812e-005]
	1180
	1181	Y = zeros(2NN, Float)
	1182	Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
	1183	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
	1184	-6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
	1185	6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
	1186	-6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
	1187	-1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
	1188	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
	1189	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
	1190	0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
	1191	-2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
	1192	0.00000000e+000, -2.57635178e-005]
	1193
	1194
	1195	domain.set_quantity('level', L, 'centroids')
	1196	domain.set_quantity('xmomentum', X, 'centroids')
	1197	domain.set_quantity('ymomentum', Y, 'centroids')
	1198
	1199	domain.check_integrity()
	1200	else:
	1201	#Evolution
	1202	for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
	1203	pass
	1204	assert allclose(domain.min_timestep, 0.0210448446782)
	1205	assert allclose(domain.max_timestep, 0.0210448446782)
	1206
	1207
	1208	#Centroids were correct but not vertices.
	1209	#Hence the check of distribute below.
	1210	assert allclose(domain.quantities['level'].centroid_values[:4],
	1211	[0.00721206,0.05352143,0.01009108,0.05354394])
	1212
	1213	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
	1214	[0.00648352,0.03685719,0.00850733,0.03687313])
	1215
	1216	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
	1217	[-0.00139463,0.0006156,-0.00060364,0.00061827])
	1218
	1219	#print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
	1220	#print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
	1221
[272]	1222	#assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
	1223	##print domain.quantities['xmomentum'].centroid_values[17], V
	1224	##print
	1225	if not V:
	1226	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
	1227	else:
	1228	assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)
[221]	1229
	1230	import copy
	1231	XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
	1232	assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
	1233
	1234	domain.distribute_to_vertices_and_edges()
	1235
[272]	1236	#assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
[221]	1237
	1238	assert allclose(domain.quantities['xmomentum'].centroid_values[17],
[272]	1239	0.0)
[221]	1240
	1241
[443]	1242	#FIXME: These numbers were from version before 25/10
	1243	#assert allclose(domain.quantities['level'].vertex_values[:4,0],
	1244	# [0.00101913,0.05352143,0.00104852,0.05354394])
[221]	1245
	1246	assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
	1247	[-0.00139463,0.0006156,-0.00060364,0.00061827])
	1248
[272]	1249
	1250	assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
	1251	[0.00064835,0.03685719,0.00085073,0.03687313])
	1252
	1253
	1254	#NB NO longer relvant:
	1255
[221]	1256	#This was the culprit. First triangles vertex 0 had an
	1257	#x-momentum of 0.0064835 instead of 0.00090581 and
	1258	#third triangle had 0.00850733 instead of 0.00088303
	1259	#print domain.quantities['xmomentum'].vertex_values[:4,0]
	1260
[272]	1261	#print domain.quantities['xmomentum'].vertex_values[:4,0]
	1262	#assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
	1263	# [0.00090581,0.03685719,0.00088303,0.03687313])
[221]	1264
	1265
[272]	1266
[221]	1267
	1268
[218]	1269	def test_bedslope_problem_first_order(self):
	1270
	1271	from mesh_factory import rectangular
	1272	from shallow_water import Domain, Reflective_boundary, Constant_height
	1273	from Numeric import array
	1274
	1275	#Create basic mesh
	1276	points, vertices, boundary = rectangular(6, 6)
	1277
	1278	#Create shallow water domain
	1279	domain = Domain(points, vertices, boundary)
	1280	domain.smooth = False
	1281	domain.default_order=1
	1282
	1283	#Bed-slope and friction
	1284	def x_slope(x, y):
	1285	return -x/3
	1286
	1287	domain.set_quantity('elevation', x_slope)
	1288
	1289	# Boundary conditions
	1290	Br = Reflective_boundary(domain)
	1291	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1292
	1293	#Initial condition
	1294	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1295	domain.check_integrity()
	1296
	1297	#Evolution
	1298	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
	1299	pass# domain.write_time()
	1300
[272]	1301	assert allclose(domain.min_timestep, 0.050010003001)
	1302	assert allclose(domain.max_timestep, 0.050010003001)
[218]	1303
	1304
	1305	def test_bedslope_problem_first_order_moresteps(self):
	1306
	1307	from mesh_factory import rectangular
	1308	from shallow_water import Domain, Reflective_boundary, Constant_height
	1309	from Numeric import array
	1310
	1311	#Create basic mesh
	1312	points, vertices, boundary = rectangular(6, 6)
	1313
	1314	#Create shallow water domain
	1315	domain = Domain(points, vertices, boundary)
	1316	domain.smooth = False
	1317	domain.default_order=1
	1318
	1319	#Bed-slope and friction
	1320	def x_slope(x, y):
	1321	return -x/3
	1322
	1323	domain.set_quantity('elevation', x_slope)
	1324
	1325	# Boundary conditions
	1326	Br = Reflective_boundary(domain)
	1327	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1328
	1329	#Initial condition
	1330	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1331	domain.check_integrity()
	1332
	1333	#Evolution
	1334	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
	1335	pass# domain.write_time()
	1336
	1337	#Data from earlier version of pyvolution
[272]	1338	#print domain.quantities['level'].centroid_values
	1339
[218]	1340	assert allclose(domain.quantities['level'].centroid_values,
[272]	1341	[-0.02998628, -0.01520652, -0.03043492,
	1342	-0.0149132, -0.03004706, -0.01476251,
	1343	-0.0298215, -0.01467976, -0.02988158,
	1344	-0.01474662, -0.03036161, -0.01442995,
	1345	-0.07624583, -0.06297061, -0.07733792,
	1346	-0.06342237, -0.07695439, -0.06289595,
	1347	-0.07635559, -0.0626065, -0.07633628,
	1348	-0.06280072, -0.07739632, -0.06386738,
	1349	-0.12161738, -0.11028239, -0.1223796,
	1350	-0.11095953, -0.12189744, -0.11048616,
	1351	-0.12074535, -0.10987605, -0.12014311,
	1352	-0.10976691, -0.12096859, -0.11087692,
	1353	-0.16868259, -0.15868061, -0.16801135,
	1354	-0.1588003, -0.16674343, -0.15813323,
	1355	-0.16457595, -0.15693826, -0.16281096,
	1356	-0.15585154, -0.16283873, -0.15540068,
	1357	-0.17450362, -0.19919913, -0.18062882,
	1358	-0.19764131, -0.17783111, -0.19407213,
	1359	-0.1736915, -0.19053624, -0.17228678,
	1360	-0.19105634, -0.17920133, -0.1968828,
	1361	-0.14244395, -0.14604641, -0.14473537,
	1362	-0.1506107, -0.14510055, -0.14919522,
	1363	-0.14175896, -0.14560798, -0.13911658,
	1364	-0.14439383, -0.13924047, -0.14829043])
	1365
[218]	1366
[229]	1367	def test_bedslope_problem_second_order_one_step(self):
[218]	1368
	1369	from mesh_factory import rectangular
	1370	from shallow_water import Domain, Reflective_boundary, Constant_height
	1371	from Numeric import array
	1372
	1373	#Create basic mesh
	1374	points, vertices, boundary = rectangular(6, 6)
	1375
	1376	#Create shallow water domain
	1377	domain = Domain(points, vertices, boundary)
	1378	domain.smooth = False
	1379	domain.default_order=2
	1380
	1381	#Bed-slope and friction at vertices (and interpolated elsewhere)
	1382	def x_slope(x, y):
	1383	return -x/3
	1384
	1385	domain.set_quantity('elevation', x_slope)
	1386
	1387	# Boundary conditions
	1388	Br = Reflective_boundary(domain)
	1389	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1390
	1391	#Initial condition
	1392	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1393	domain.check_integrity()
	1394
	1395	assert allclose(domain.quantities['level'].centroid_values,
	1396	[0.01296296, 0.03148148, 0.01296296,
	1397	0.03148148, 0.01296296, 0.03148148,
	1398	0.01296296, 0.03148148, 0.01296296,
	1399	0.03148148, 0.01296296, 0.03148148,
	1400	-0.04259259, -0.02407407, -0.04259259,
	1401	-0.02407407, -0.04259259, -0.02407407,
	1402	-0.04259259, -0.02407407, -0.04259259,
	1403	-0.02407407, -0.04259259, -0.02407407,
	1404	-0.09814815, -0.07962963, -0.09814815,
	1405	-0.07962963, -0.09814815, -0.07962963,
	1406	-0.09814815, -0.07962963, -0.09814815,
	1407	-0.07962963, -0.09814815, -0.07962963,
	1408	-0.1537037 , -0.13518519, -0.1537037,
	1409	-0.13518519, -0.1537037, -0.13518519,
	1410	-0.1537037 , -0.13518519, -0.1537037,
	1411	-0.13518519, -0.1537037, -0.13518519,
	1412	-0.20925926, -0.19074074, -0.20925926,
	1413	-0.19074074, -0.20925926, -0.19074074,
	1414	-0.20925926, -0.19074074, -0.20925926,
	1415	-0.19074074, -0.20925926, -0.19074074,
	1416	-0.26481481, -0.2462963, -0.26481481,
	1417	-0.2462963, -0.26481481, -0.2462963,
	1418	-0.26481481, -0.2462963, -0.26481481,
	1419	-0.2462963, -0.26481481, -0.2462963])
	1420
	1421
	1422	#print domain.quantities['level'].extrapolate_second_order()
	1423	#domain.distribute_to_vertices_and_edges()
	1424	#print domain.quantities['level'].vertex_values[:,0]
	1425
	1426	#Evolution
[229]	1427	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
[218]	1428	#domain.write_time()
	1429	pass
	1430
	1431
[229]	1432	#print domain.quantities['level'].centroid_values
	1433	assert allclose(domain.quantities['level'].centroid_values,
	1434	[0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
	1435	0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
	1436	0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
	1437	-0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
	1438	-0.04186556, -0.0248011, -0.04186556, -0.02255593,
	1439	-0.09966629, -0.08035666, -0.09742112, -0.08035666,
	1440	-0.09742112, -0.08035666, -0.09742112, -0.08035666,
	1441	-0.09742112, -0.08035666, -0.09742112, -0.07811149,
	1442	-0.15522185, -0.13591222, -0.15297667, -0.13591222,
	1443	-0.15297667, -0.13591222, -0.15297667, -0.13591222,
	1444	-0.15297667, -0.13591222, -0.15297667, -0.13366704,
	1445	-0.2107774, -0.19146777, -0.20853223, -0.19146777,
	1446	-0.20853223, -0.19146777, -0.20853223, -0.19146777,
	1447	-0.20853223, -0.19146777, -0.20853223, -0.1892226,
	1448	-0.26120669, -0.24776246, -0.25865535, -0.24776246,
	1449	-0.25865535, -0.24776246, -0.25865535, -0.24776246,
	1450	-0.25865535, -0.24776246, -0.25865535, -0.24521113])
	1451
	1452
	1453
	1454	def test_bedslope_problem_second_order_two_steps(self):
	1455
	1456	from mesh_factory import rectangular
	1457	from shallow_water import Domain, Reflective_boundary, Constant_height
	1458	from Numeric import array
	1459
	1460	#Create basic mesh
	1461	points, vertices, boundary = rectangular(6, 6)
	1462
	1463	#Create shallow water domain
	1464	domain = Domain(points, vertices, boundary)
	1465	domain.smooth = False
	1466	domain.default_order=2
	1467
	1468	#Bed-slope and friction at vertices (and interpolated elsewhere)
	1469	def x_slope(x, y):
	1470	return -x/3
	1471
	1472	domain.set_quantity('elevation', x_slope)
	1473
	1474	# Boundary conditions
	1475	Br = Reflective_boundary(domain)
	1476	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1477
	1478	#Initial condition
	1479	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1480	domain.check_integrity()
	1481
	1482	assert allclose(domain.quantities['level'].centroid_values,
	1483	[0.01296296, 0.03148148, 0.01296296,
	1484	0.03148148, 0.01296296, 0.03148148,
	1485	0.01296296, 0.03148148, 0.01296296,
	1486	0.03148148, 0.01296296, 0.03148148,
	1487	-0.04259259, -0.02407407, -0.04259259,
	1488	-0.02407407, -0.04259259, -0.02407407,
	1489	-0.04259259, -0.02407407, -0.04259259,
	1490	-0.02407407, -0.04259259, -0.02407407,
	1491	-0.09814815, -0.07962963, -0.09814815,
	1492	-0.07962963, -0.09814815, -0.07962963,
	1493	-0.09814815, -0.07962963, -0.09814815,
	1494	-0.07962963, -0.09814815, -0.07962963,
	1495	-0.1537037 , -0.13518519, -0.1537037,
	1496	-0.13518519, -0.1537037, -0.13518519,
	1497	-0.1537037 , -0.13518519, -0.1537037,
	1498	-0.13518519, -0.1537037, -0.13518519,
	1499	-0.20925926, -0.19074074, -0.20925926,
	1500	-0.19074074, -0.20925926, -0.19074074,
	1501	-0.20925926, -0.19074074, -0.20925926,
	1502	-0.19074074, -0.20925926, -0.19074074,
	1503	-0.26481481, -0.2462963, -0.26481481,
	1504	-0.2462963, -0.26481481, -0.2462963,
	1505	-0.26481481, -0.2462963, -0.26481481,
	1506	-0.2462963, -0.26481481, -0.2462963])
	1507
	1508
	1509	#print domain.quantities['level'].extrapolate_second_order()
	1510	#domain.distribute_to_vertices_and_edges()
	1511	#print domain.quantities['level'].vertex_values[:,0]
[218]	1512
[229]	1513	#Evolution
	1514	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
	1515	pass
	1516
	1517
	1518	#Data from earlier version of pyvolution ft=0.1
[218]	1519	assert allclose(domain.min_timestep, 0.0376895634803)
	1520	assert allclose(domain.max_timestep, 0.0415635655309)
	1521
	1522
[229]	1523	assert allclose(domain.quantities['level'].centroid_values,
	1524	[0.00855788, 0.01575204, 0.00994606, 0.01717072,
	1525	0.01005985, 0.01716362, 0.01005985, 0.01716299,
	1526	0.01007098, 0.01736248, 0.01216452, 0.02026776,
	1527	-0.04462374, -0.02479045, -0.04199789, -0.0229465,
	1528	-0.04184033, -0.02295693, -0.04184013, -0.02295675,
	1529	-0.04184486, -0.0228168, -0.04028876, -0.02036486,
	1530	-0.10029444, -0.08170809, -0.09772846, -0.08021704,
	1531	-0.09760006, -0.08022143, -0.09759984, -0.08022124,
	1532	-0.09760261, -0.08008893, -0.09603914, -0.07758209,
	1533	-0.15584152, -0.13723138, -0.15327266, -0.13572906,
	1534	-0.15314427, -0.13573349, -0.15314405, -0.13573331,
	1535	-0.15314679, -0.13560104, -0.15158523, -0.13310701,
	1536	-0.21208605, -0.19283913, -0.20955631, -0.19134189,
	1537	-0.20942821, -0.19134598, -0.20942799, -0.1913458,
	1538	-0.20943005, -0.19120952, -0.20781177, -0.18869401,
	1539	-0.25384082, -0.2463294, -0.25047649, -0.24464654,
	1540	-0.25031159, -0.24464253, -0.25031112, -0.24464253,
	1541	-0.25031463, -0.24454764, -0.24885323, -0.24286438])
	1542
	1543
	1544
	1545
	1546	def test_bedslope_problem_second_order_two_yieldsteps(self):
	1547
	1548	from mesh_factory import rectangular
	1549	from shallow_water import Domain, Reflective_boundary, Constant_height
	1550	from Numeric import array
	1551
	1552	#Create basic mesh
	1553	points, vertices, boundary = rectangular(6, 6)
	1554
	1555	#Create shallow water domain
	1556	domain = Domain(points, vertices, boundary)
	1557	domain.smooth = False
	1558	domain.default_order=2
	1559
	1560	#Bed-slope and friction at vertices (and interpolated elsewhere)
	1561	def x_slope(x, y):
	1562	return -x/3
	1563
	1564	domain.set_quantity('elevation', x_slope)
	1565
	1566	# Boundary conditions
	1567	Br = Reflective_boundary(domain)
	1568	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1569
	1570	#Initial condition
	1571	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1572	domain.check_integrity()
	1573
	1574	assert allclose(domain.quantities['level'].centroid_values,
	1575	[0.01296296, 0.03148148, 0.01296296,
	1576	0.03148148, 0.01296296, 0.03148148,
	1577	0.01296296, 0.03148148, 0.01296296,
	1578	0.03148148, 0.01296296, 0.03148148,
	1579	-0.04259259, -0.02407407, -0.04259259,
	1580	-0.02407407, -0.04259259, -0.02407407,
	1581	-0.04259259, -0.02407407, -0.04259259,
	1582	-0.02407407, -0.04259259, -0.02407407,
	1583	-0.09814815, -0.07962963, -0.09814815,
	1584	-0.07962963, -0.09814815, -0.07962963,
	1585	-0.09814815, -0.07962963, -0.09814815,
	1586	-0.07962963, -0.09814815, -0.07962963,
	1587	-0.1537037 , -0.13518519, -0.1537037,
	1588	-0.13518519, -0.1537037, -0.13518519,
	1589	-0.1537037 , -0.13518519, -0.1537037,
	1590	-0.13518519, -0.1537037, -0.13518519,
	1591	-0.20925926, -0.19074074, -0.20925926,
	1592	-0.19074074, -0.20925926, -0.19074074,
	1593	-0.20925926, -0.19074074, -0.20925926,
	1594	-0.19074074, -0.20925926, -0.19074074,
	1595	-0.26481481, -0.2462963, -0.26481481,
	1596	-0.2462963, -0.26481481, -0.2462963,
	1597	-0.26481481, -0.2462963, -0.26481481,
	1598	-0.2462963, -0.26481481, -0.2462963])
	1599
	1600
	1601	#print domain.quantities['level'].extrapolate_second_order()
	1602	#domain.distribute_to_vertices_and_edges()
	1603	#print domain.quantities['level'].vertex_values[:,0]
[218]	1604
[229]	1605	#Evolution
	1606	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): #0.05??
	1607	#domain.write_time()
	1608	pass
	1609
	1610
	1611
	1612	assert allclose(domain.quantities['level'].centroid_values,
	1613	[0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
	1614	0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
	1615	0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
	1616	-0.0229465, -0.04184033, -0.02295693, -0.04184013,
	1617	-0.02295675, -0.04184486, -0.0228168, -0.04028876,
	1618	-0.02036486, -0.10029444, -0.08170809, -0.09772846,
	1619	-0.08021704, -0.09760006, -0.08022143, -0.09759984,
	1620	-0.08022124, -0.09760261, -0.08008893, -0.09603914,
	1621	-0.07758209, -0.15584152, -0.13723138, -0.15327266,
	1622	-0.13572906, -0.15314427, -0.13573349, -0.15314405,
	1623	-0.13573331, -0.15314679, -0.13560104, -0.15158523,
	1624	-0.13310701, -0.21208605, -0.19283913, -0.20955631,
	1625	-0.19134189, -0.20942821, -0.19134598, -0.20942799,
	1626	-0.1913458, -0.20943005, -0.19120952, -0.20781177,
	1627	-0.18869401, -0.25384082, -0.2463294, -0.25047649,
	1628	-0.24464654, -0.25031159, -0.24464253, -0.25031112,
	1629	-0.24464253, -0.25031463, -0.24454764, -0.24885323,
	1630	-0.24286438])
	1631
	1632
	1633
	1634	def test_bedslope_problem_second_order_more_steps(self):
	1635
	1636	from mesh_factory import rectangular
	1637	from shallow_water import Domain, Reflective_boundary, Constant_height
	1638	from Numeric import array
	1639
	1640	#Create basic mesh
	1641	points, vertices, boundary = rectangular(6, 6)
	1642
	1643	#Create shallow water domain
	1644	domain = Domain(points, vertices, boundary)
	1645	domain.smooth = False
	1646	domain.default_order=2
	1647
	1648	#Bed-slope and friction at vertices (and interpolated elsewhere)
	1649	def x_slope(x, y):
	1650	return -x/3
	1651
	1652	domain.set_quantity('elevation', x_slope)
	1653
	1654	# Boundary conditions
	1655	Br = Reflective_boundary(domain)
	1656	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1657
	1658	#Initial condition
	1659	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1660	domain.check_integrity()
	1661
	1662	assert allclose(domain.quantities['level'].centroid_values,
	1663	[0.01296296, 0.03148148, 0.01296296,
	1664	0.03148148, 0.01296296, 0.03148148,
	1665	0.01296296, 0.03148148, 0.01296296,
	1666	0.03148148, 0.01296296, 0.03148148,
	1667	-0.04259259, -0.02407407, -0.04259259,
	1668	-0.02407407, -0.04259259, -0.02407407,
	1669	-0.04259259, -0.02407407, -0.04259259,
	1670	-0.02407407, -0.04259259, -0.02407407,
	1671	-0.09814815, -0.07962963, -0.09814815,
	1672	-0.07962963, -0.09814815, -0.07962963,
	1673	-0.09814815, -0.07962963, -0.09814815,
	1674	-0.07962963, -0.09814815, -0.07962963,
	1675	-0.1537037 , -0.13518519, -0.1537037,
	1676	-0.13518519, -0.1537037, -0.13518519,
	1677	-0.1537037 , -0.13518519, -0.1537037,
	1678	-0.13518519, -0.1537037, -0.13518519,
	1679	-0.20925926, -0.19074074, -0.20925926,
	1680	-0.19074074, -0.20925926, -0.19074074,
	1681	-0.20925926, -0.19074074, -0.20925926,
	1682	-0.19074074, -0.20925926, -0.19074074,
	1683	-0.26481481, -0.2462963, -0.26481481,
	1684	-0.2462963, -0.26481481, -0.2462963,
	1685	-0.26481481, -0.2462963, -0.26481481,
	1686	-0.2462963, -0.26481481, -0.2462963])
	1687
	1688
	1689	#print domain.quantities['level'].extrapolate_second_order()
	1690	#domain.distribute_to_vertices_and_edges()
	1691	#print domain.quantities['level'].vertex_values[:,0]
	1692
	1693	#Evolution
	1694	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
	1695	pass
	1696
	1697
	1698	assert allclose(domain.quantities['level'].centroid_values,
	1699	[-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
	1700	-0.02901975, -0.0141361, -0.02898816, -0.01418135, -0.02961409, -0.01403487,
	1701	-0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
	1702	-0.07571145, -0.06235231, -0.0756817, -0.06245309, -0.07652292, -0.06289946,
	1703	-0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934, -0.1107174,
	1704	-0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
	1705	-0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493, -0.15697919,
	1706	-0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
	1707	-0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
	1708	-0.1833593, -0.19505871, -0.1818806, -0.19418042, -0.18586159, -0.19576946,
	1709	-0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
	1710	-0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
	1711
	1712	assert allclose(domain.quantities['xmomentum'].centroid_values,
	1713	[ 0.00831121, 0.00317948, 0.00731797, 0.00334939, 0.00764717, 0.00348053,
	1714	0.00788729, 0.00356522, 0.00780649, 0.00341919, 0.00693525, 0.00310375,
	1715	0.02166196, 0.01421475, 0.02017737, 0.01316839, 0.02037015, 0.01368659,
	1716	0.02106, 0.01399161, 0.02074514, 0.01354935, 0.01887407, 0.0123113,
	1717	0.03775083, 0.02855197, 0.03689337, 0.02759782, 0.03732848, 0.02812072,
	1718	0.03872545, 0.02913348, 0.03880939, 0.02803804, 0.03546499, 0.0260039,
	1719	0.0632131, 0.04730634, 0.0576324, 0.04592336, 0.05790921, 0.04690514,
	1720	0.05986467, 0.04871165, 0.06170068, 0.04811572, 0.05657041, 0.04416292,
	1721	0.08489642, 0.07188097, 0.07835261, 0.06843406, 0.07986412, 0.0698247,
	1722	0.08201071, 0.07216756, 0.08378418, 0.07273624, 0.080399, 0.06645841,
	1723	0.01631548, 0.04691608, 0.0206632, 0.044409, 0.02115518, 0.04560305,
	1724	0.02160608, 0.04663725, 0.02174734, 0.04795559, 0.02281427, 0.05667111])
	1725
	1726
	1727	assert allclose(domain.quantities['ymomentum'].centroid_values,
	1728	[ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
	1729	-9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
	1730	-1.99445058e-004, -3.28493467e-004, 6.68217349e-005, -8.42042805e-006,
	1731	5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
	1732	-8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
	1733	-2.70160173e-004, -5.35485454e-004, 2.60780997e-004, 3.12145471e-005,
	1734	5.16189608e-004, 1.07069062e-004, 9.29989252e-005, -3.71211119e-004,
	1735	1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
	1736	-4.74025708e-004, -6.94463009e-004, 6.15092843e-005, 2.22106820e-004,
	1737	-6.29589294e-004, 2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
	1738	-4.17914641e-004, 6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
	1739	-1.67424308e-003, -7.39485066e-004, -1.59966988e-003, 5.68262838e-005,
	1740	-1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
	1741	-1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
	1742	-2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
	1743	-1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
	1744	-6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
	1745	-4.50964850e-003, -3.06319963e-003, 6.08950810e-004, -4.79537921e-004])
	1746
	1747
	1748
	1749
	1750	def test_temp_play(self):
	1751
	1752	from mesh_factory import rectangular
	1753	from shallow_water import Domain, Reflective_boundary, Constant_height
	1754	from Numeric import array
	1755
	1756	#Create basic mesh
	1757	points, vertices, boundary = rectangular(5, 5)
	1758
	1759	#Create shallow water domain
	1760	domain = Domain(points, vertices, boundary)
	1761	domain.smooth = False
	1762	domain.default_order=2
	1763
	1764	#Bed-slope and friction at vertices (and interpolated elsewhere)
	1765	def x_slope(x, y):
	1766	return -x/3
	1767
	1768	domain.set_quantity('elevation', x_slope)
	1769
	1770	# Boundary conditions
	1771	Br = Reflective_boundary(domain)
	1772	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
	1773
	1774	#Initial condition
	1775	domain.set_quantity('level', Constant_height(x_slope, 0.05))
	1776	domain.check_integrity()
	1777
	1778	#Evolution
	1779	for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
	1780	pass
	1781
	1782	assert allclose(domain.quantities['level'].centroid_values[:4],
	1783	[0.00206836, 0.01296714, 0.00363415, 0.01438924])
	1784	assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
	1785	[0.01360154, 0.00671133, 0.01264578, 0.00648503])
	1786	assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
	1787	[-1.19201077e-003, -7.23647546e-004,
	1788	-6.39083123e-005, 6.29815168e-005])
	1789
	1790
[232]	1791	def test_complex_bed(self):
	1792	#No friction is tested here
[229]	1793
[232]	1794	from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
	1795	Transmissive_boundary, Time_boundary,\
	1796	Weir_simple as Weir, Constant_height
	1797
	1798	from mesh_factory import rectangular
	1799	from Numeric import array
[218]	1800
[232]	1801	N = 12
	1802	points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6,
	1803	origin=(-0.07, 0))
	1804
	1805
	1806	domain = Domain(points, vertices, boundary)
	1807	domain.smooth = False
	1808	domain.visualise = False
	1809	domain.default_order=2
	1810
	1811
	1812	inflow_stage = 0.1
	1813	Z = Weir(inflow_stage)
	1814	domain.set_quantity('elevation', Z)
	1815
	1816	Br = Reflective_boundary(domain)
	1817	Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
	1818	domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
	1819
	1820	domain.set_quantity('level', Constant_height(Z, 0.))
	1821
	1822	for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
	1823	pass
	1824
[272]	1825
	1826	#print domain.quantities['level'].centroid_values
	1827
[443]	1828	#FIXME: These numbers were from version before 25/10
	1829	#assert allclose(domain.quantities['level'].centroid_values,
	1830	# [3.95822638e-002, 5.61022588e-002, 4.66437868e-002, 5.73081011e-002,
	1831	# 4.72394613e-002, 5.74684939e-002, 4.74309483e-002, 5.77458084e-002,
	1832	# 4.80628177e-002, 5.85656225e-002, 4.90498542e-002, 6.02609831e-002,
	1833	# 1.18470315e-002, 1.75136443e-002, 1.18035266e-002, 2.15565695e-002,
	1834	# 1.31620268e-002, 2.14351640e-002, 1.32351076e-002, 2.15450687e-002,
	1835	# 1.36414028e-002, 2.24274619e-002, 1.51689511e-002, 2.21789655e-002,
	1836	# -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
	1837	# -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
	1838	# -8.10292716e-003, -3.88584984e-004, -7.35226124e-003, 2.73985295e-004,
	1839	# 1.86166683e-001, 8.74070369e-002, 1.86166712e-001, 8.74035875e-002,
	1840	# 6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
	1841	# 6.11666725e-002, 8.73846774e-002, 1.86166697e-001, 8.74171550e-002,
	1842	# -4.83333333e-002, 1.18333333e-001, -4.83333333e-002, 1.18333333e-001,
	1843	# -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
	1844	# -1.73333333e-001, -6.66666667e-003, -4.83333333e-002, 1.18333333e-001,
	1845	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
	1846	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
	1847	# -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
	1848	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
	1849	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
	1850	# -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
	1851	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
	1852	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
	1853	# -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
	1854	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
	1855	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
	1856	# -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
	1857	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
	1858	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
	1859	# -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
	1860	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
	1861	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
	1862	# -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
	1863	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
	1864	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
	1865	# -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
[232]	1866
[272]	1867
[232]	1868
	1869
[218]	1870	#-------------------------------------------------------------
	1871	if __name__ == "__main__":
	1872	suite = unittest.makeSuite(TestCase,'test')
	1873	runner = unittest.TextTestRunner()
	1874	runner.run(suite)
	1875

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

Download in other formats: