Context Navigation

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

Last change on this file since 750 was 750, checked in by ole, 20 years ago
Test that initial values are indeed output at t == 0 (just to make sure). Also a check that bounsry is created before starting evolve.
File size: 85.5 KB

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

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

Download in other formats: