Context Navigation

source: anuga_core/source/anuga/culvert_flows/test_culvert_class.py @ 7492

Last change on this file since 7492 was 7492, checked in by steve, 15 years ago

Changed the order of update_timestep and compute_forcing_terms.
Can now update the parameter domain.flux_timestep within
the computation of the forcing terms, so as to control stability of
the evolver.

Hopefully this will help in reducing the oscillations seen within the
culvert computation.

File size: 20.8 KB

Line
1	#!/usr/bin/env python
2
3
4	import unittest
5	import os.path
6	import sys
7
8	from anuga.utilities.system_tools import get_pathname_from_package
9	from anuga.utilities.polygon import Polygon_function
10
11	from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
12	from anuga.abstract_2d_finite_volumes.quantity import Quantity
13
14	from anuga.shallow_water import Domain, Reflective_boundary,\
15	Dirichlet_boundary,\
16	Transmissive_boundary, Time_boundary
17
18	from anuga.culvert_flows.culvert_class import Culvert_flow, Culvert_flow_rating, Culvert_flow_energy
19	from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
20
21	from math import pi,pow,sqrt
22
23	import numpy as num
24
25
26	class Test_Culvert(unittest.TestCase):
27	def setUp(self):
28	pass
29
30	def tearDown(self):
31	pass
32
33
34	def test_that_culvert_runs_rating(self):
35	"""test_that_culvert_runs_rating
36
37	This test exercises the culvert and checks values outside rating curve
38	are dealt with
39	"""
40
41	path = get_pathname_from_package('anuga.culvert_flows')
42
43	length = 40.
44	width = 5.
45
46	dx = dy = 1 # Resolution: Length of subdivisions on both axes
47
48	points, vertices, boundary = rectangular_cross(int(length/dx),
49	int(width/dy),
50	len1=length,
51	len2=width)
52	domain = Domain(points, vertices, boundary)
53	domain.set_name('Test_culvert') # Output name
54	domain.set_default_order(2)
55
56
57	#----------------------------------------------------------------------
58	# Setup initial conditions
59	#----------------------------------------------------------------------
60
61	def topography(x, y):
62	"""Set up a weir
63
64	A culvert will connect either side
65	"""
66	# General Slope of Topography
67	z=-x/1000
68
69	N = len(x)
70	for i in range(N):
71
72	# Sloping Embankment Across Channel
73	if 5.0 < x[i] < 10.1:
74	# Cut Out Segment for Culvert face
75	if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0:
76	z[i]=z[i]
77	else:
78	z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face
79	if 10.0 < x[i] < 12.1:
80	z[i] += 2.5 # Flat Crest of Embankment
81	if 12.0 < x[i] < 14.5:
82	# Cut Out Segment for Culvert face
83	if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5:
84	z[i]=z[i]
85	else:
86	z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face
87
88
89	return z
90
91
92	domain.set_quantity('elevation', topography)
93	domain.set_quantity('friction', 0.01) # Constant friction
94	domain.set_quantity('stage',
95	expression='elevation') # Dry initial condition
96
97	filename=os.path.join(path, 'example_rating_curve.csv')
98	culvert = Culvert_flow(domain,
99	culvert_description_filename=filename,
100	end_point0=[9.0, 2.5],
101	end_point1=[13.0, 2.5],
102	use_velocity_head=True,
103	verbose=False)
104
105	domain.forcing_terms.append(culvert)
106
107
108	#-----------------------------------------------------------------------
109	# Setup boundary conditions
110	#-----------------------------------------------------------------------
111
112	# Inflow based on Flow Depth and Approaching Momentum
113	Bi = Dirichlet_boundary([0.0, 0.0, 0.0])
114	Br = Reflective_boundary(domain) # Solid reflective wall
115	Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow
116
117	# Upstream and downstream conditions that will exceed the rating curve
118	# I.e produce delta_h outside the range [0, 10] specified in the the
119	# file example_rating_curve.csv
120	Btus = Time_boundary(domain, lambda t: [100num.sin(2pi*(t-4)/10), 0.0, 0.0])
121	Btds = Time_boundary(domain, lambda t: [-5(num.cos(2pi*(t-4)/20)), 0.0, 0.0])
122	domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
123
124
125	#-----------------------------------------------------------------------
126	# Evolve system through time
127	#-----------------------------------------------------------------------
128
129	min_delta_w = sys.maxint
130	max_delta_w = -min_delta_w
131	for t in domain.evolve(yieldstep = 1, finaltime = 25):
132	delta_w = culvert.inlet.stage - culvert.outlet.stage
133
134	if delta_w > max_delta_w: max_delta_w = delta_w
135	if delta_w < min_delta_w: min_delta_w = delta_w
136
137	pass
138
139	# Check that extreme values in rating curve have been exceeded
140	# so that we know that condition has been exercised
141	assert min_delta_w < 0
142	assert max_delta_w > 10
143
144
145	def test_that_culvert_dry_bed_rating_does_not_produce_flow(self):
146	"""test_that_culvert_in_dry_bed_does_not_produce_flow(self):
147
148	Test that culvert on a sloping dry bed doesn't produce flows
149	although there will be a 'pressure' head due to delta_w > 0
150
151	This one is using the rating curve variant
152	"""
153
154	path = get_pathname_from_package('anuga.culvert_flows')
155
156	length = 40.
157	width = 5.
158
159	dx = dy = 1 # Resolution: Length of subdivisions on both axes
160
161	points, vertices, boundary = rectangular_cross(int(length/dx),
162	int(width/dy),
163	len1=length,
164	len2=width)
165	domain = Domain(points, vertices, boundary)
166	domain.set_name('Test_culvert_dry') # Output name
167	domain.set_default_order(2)
168
169
170	#----------------------------------------------------------------------
171	# Setup initial conditions
172	#----------------------------------------------------------------------
173
174	def topography(x, y):
175	"""Set up a weir
176
177	A culvert will connect either side
178	"""
179	# General Slope of Topography
180	z=-x/1000
181
182	N = len(x)
183	for i in range(N):
184
185	# Sloping Embankment Across Channel
186	if 5.0 < x[i] < 10.1:
187	# Cut Out Segment for Culvert face
188	if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0:
189	z[i]=z[i]
190	else:
191	z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face
192	if 10.0 < x[i] < 12.1:
193	z[i] += 2.5 # Flat Crest of Embankment
194	if 12.0 < x[i] < 14.5:
195	# Cut Out Segment for Culvert face
196	if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5:
197	z[i]=z[i]
198	else:
199	z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face
200
201
202	return z
203
204
205	domain.set_quantity('elevation', topography)
206	domain.set_quantity('friction', 0.01) # Constant friction
207	domain.set_quantity('stage',
208	expression='elevation') # Dry initial condition
209
210
211	filename = os.path.join(path, 'example_rating_curve.csv')
212	culvert = Culvert_flow(domain,
213	culvert_description_filename=filename,
214	end_point0=[9.0, 2.5],
215	end_point1=[13.0, 2.5],
216	verbose=False)
217
218	domain.forcing_terms.append(culvert)
219
220
221	#-----------------------------------------------------------------------
222	# Setup boundary conditions
223	#-----------------------------------------------------------------------
224
225	# Inflow based on Flow Depth and Approaching Momentum
226
227	Br = Reflective_boundary(domain) # Solid reflective wall
228	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
229
230
231	#-----------------------------------------------------------------------
232	# Evolve system through time
233	#-----------------------------------------------------------------------
234
235	ref_volume = domain.get_quantity('stage').get_integral()
236	for t in domain.evolve(yieldstep = 1, finaltime = 25):
237	new_volume = domain.get_quantity('stage').get_integral()
238
239	msg = 'Total volume has changed'
240	assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
241	pass
242
243
244
245	def test_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self):
246	"""test_that_culvert_rating_limits_flow_in_shallow_inlet_condition
247
248	Test that culvert on a sloping dry bed limits flows when very little water
249	is present at inlet
250
251	This one is using the rating curve variant
252	"""
253
254	path = get_pathname_from_package('anuga.culvert_flows')
255
256	length = 40.
257	width = 5.
258
259	dx = dy = 1 # Resolution: Length of subdivisions on both axes
260
261	points, vertices, boundary = rectangular_cross(int(length/dx),
262	int(width/dy),
263	len1=length,
264	len2=width)
265	domain = Domain(points, vertices, boundary)
266	domain.set_name('Test_culvert_shallow') # Output name
267	domain.set_default_order(2)
268
269
270	#----------------------------------------------------------------------
271	# Setup initial conditions
272	#----------------------------------------------------------------------
273
274	def topography(x, y):
275	"""Set up a weir
276
277	A culvert will connect either side
278	"""
279	# General Slope of Topography
280	z=-x/1000
281
282	N = len(x)
283	for i in range(N):
284
285	# Sloping Embankment Across Channel
286	if 5.0 < x[i] < 10.1:
287	# Cut Out Segment for Culvert face
288	if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0:
289	z[i]=z[i]
290	else:
291	z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face
292	if 10.0 < x[i] < 12.1:
293	z[i] += 2.5 # Flat Crest of Embankment
294	if 12.0 < x[i] < 14.5:
295	# Cut Out Segment for Culvert face
296	if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5:
297	z[i]=z[i]
298	else:
299	z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face
300
301
302	return z
303
304
305	domain.set_quantity('elevation', topography)
306	domain.set_quantity('friction', 0.01) # Constant friction
307	domain.set_quantity('stage',
308	expression='elevation + 0.1') # Shallow initial condition
309
310
311	filename = os.path.join(path, 'example_rating_curve.csv')
312	culvert = Culvert_flow(domain,
313	culvert_description_filename=filename,
314	end_point0=[9.0, 2.5],
315	end_point1=[13.0, 2.5],
316	trigger_depth=0.01,
317	verbose=False)
318
319	domain.forcing_terms.append(culvert)
320
321
322	#-----------------------------------------------------------------------
323	# Setup boundary conditions
324	#-----------------------------------------------------------------------
325
326	# Inflow based on Flow Depth and Approaching Momentum
327
328	Br = Reflective_boundary(domain) # Solid reflective wall
329	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
330
331
332
333	#-----------------------------------------------------------------------
334	# Evolve system through time
335	#-----------------------------------------------------------------------
336
337	ref_volume = domain.get_quantity('stage').get_integral()
338	for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
339	new_volume = domain.get_quantity('stage').get_integral()
340
341
342	msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
343	% (new_volume, ref_volume))
344	assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
345
346
347	return
348	# Now try this for a range of other depths
349	for depth in [1.0, 0.5, 0.2, 0.1, 0.05]:
350	domain.set_time(0.0)
351
352	domain.set_quantity('stage',
353	expression='elevation + %f' % depth)
354
355
356	ref_volume = domain.get_quantity('stage').get_integral()
357	for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
358	new_volume = domain.get_quantity('stage').get_integral()
359
360	msg = 'Total volume has changed: Is %.2f m^3 should have been %.2f m^3'\
361	% (new_volume, ref_volume)
362
363	assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
364
365
366
367	def test_that_culvert_dry_bed_boyd_does_not_produce_flow(self):
368	"""test_that_culvert_in_dry_bed_boyd_does_not_produce_flow(self):
369
370	Test that culvert on a sloping dry bed doesn't produce flows
371	although there will be a 'pressure' head due to delta_w > 0
372
373	This one is using the 'Boyd' variant
374	"""
375
376	path = get_pathname_from_package('anuga.culvert_flows')
377
378	length = 40.
379	width = 5.
380
381	dx = dy = 1 # Resolution: Length of subdivisions on both axes
382
383	points, vertices, boundary = rectangular_cross(int(length/dx),
384	int(width/dy),
385	len1=length,
386	len2=width)
387	domain = Domain(points, vertices, boundary)
388	domain.set_name('Test_culvert_dry') # Output name
389	domain.set_default_order(2)
390
391
392	#----------------------------------------------------------------------
393	# Setup initial conditions
394	#----------------------------------------------------------------------
395
396	def topography(x, y):
397	"""Set up a weir
398
399	A culvert will connect either side
400	"""
401	# General Slope of Topography
402	z=-x/1000
403
404	N = len(x)
405	for i in range(N):
406
407	# Sloping Embankment Across Channel
408	if 5.0 < x[i] < 10.1:
409	# Cut Out Segment for Culvert face
410	if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0:
411	z[i]=z[i]
412	else:
413	z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face
414	if 10.0 < x[i] < 12.1:
415	z[i] += 2.5 # Flat Crest of Embankment
416	if 12.0 < x[i] < 14.5:
417	# Cut Out Segment for Culvert face
418	if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5:
419	z[i]=z[i]
420	else:
421	z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face
422
423
424	return z
425
426
427	domain.set_quantity('elevation', topography)
428	domain.set_quantity('friction', 0.01) # Constant friction
429	domain.set_quantity('stage',
430	expression='elevation') # Dry initial condition
431
432
433	filename = os.path.join(path, 'example_rating_curve.csv')
434
435
436	culvert = Culvert_flow(domain,
437	label='Culvert No. 1',
438	description='This culvert is a test unit 1.2m Wide by 0.75m High',
439	end_point0=[9.0, 2.5],
440	end_point1=[13.0, 2.5],
441	width=1.20,height=0.75,
442	culvert_routine=boyd_generalised_culvert_model,
443	number_of_barrels=1,
444	update_interval=2,
445	verbose=True)
446
447	domain.forcing_terms.append(culvert)
448
449
450	#-----------------------------------------------------------------------
451	# Setup boundary conditions
452	#-----------------------------------------------------------------------
453
454	# Inflow based on Flow Depth and Approaching Momentum
455
456	Br = Reflective_boundary(domain) # Solid reflective wall
457	domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
458
459
460	#-----------------------------------------------------------------------
461	# Evolve system through time
462	#-----------------------------------------------------------------------
463
464	ref_volume = domain.get_quantity('stage').get_integral()
465	for t in domain.evolve(yieldstep = 1, finaltime = 25):
466
467	new_volume = domain.get_quantity('stage').get_integral()
468
469	msg = 'Total volume has changed'
470	assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
471	pass
472
473
474
475
476
477	def test_predicted_boyd_flow(self):
478	"""test_predicted_boyd_flow
479
480	Test that flows predicted by the boyd method are consistent with what what
481	is calculated in engineering codes.
482	The data was supplied by Petar Milevski
483	"""
484
485	# FIXME(Ole) this is nowhere near finished
486	path = get_pathname_from_package('anuga.culvert_flows')
487
488	length = 12.
489	width = 5.
490
491	dx = dy = 0.5 # Resolution: Length of subdivisions on both axes
492
493	points, vertices, boundary = rectangular_cross(int(length/dx),
494	int(width/dy),
495	len1=length,
496	len2=width)
497	domain = Domain(points, vertices, boundary)
498
499	domain.set_name('test_culvert') # Output name
500	domain.set_default_order(2)
501
502
503	#----------------------------------------------------------------------
504	# Setup initial conditions
505	#----------------------------------------------------------------------
506
507	def topography(x, y):
508	# General Slope of Topography
509	z=-x/10
510
511	return z
512
513
514	domain.set_quantity('elevation', topography)
515	domain.set_quantity('friction', 0.01) # Constant friction
516	domain.set_quantity('stage', expression='elevation')
517
518
519	Q0 = domain.get_quantity('stage')
520	Q1 = Quantity(domain)
521
522	# Add depths to stage
523	head_water_depth = 0.169
524	tail_water_depth = 0.089
525
526	inlet_poly = [[0,0], [6,0], [6,5], [0,5]]
527	outlet_poly = [[6,0], [12,0], [12,5], [6,5]]
528
529	Q1.set_values(Polygon_function([(inlet_poly, head_water_depth),
530	(outlet_poly, tail_water_depth)]))
531
532	domain.set_quantity('stage', Q0 + Q1)
533
534
535
536	culvert = Culvert_flow(domain,
537	label='Test culvert',
538	description='4 m test culvert',
539	end_point0=[4.0, 2.5],
540	end_point1=[8.0, 2.5],
541	width=1.20,
542	height=0.75,
543	culvert_routine=boyd_generalised_culvert_model,
544	number_of_barrels=1,
545	verbose=True)
546
547
548	domain.forcing_terms.append(culvert)
549
550	# Call
551	culvert(domain)
552
553
554	#-------------------------------------------------------------
555
556	if __name__ == "__main__":
557	suite = unittest.makeSuite(Test_Culvert, 'test')
558	runner = unittest.TextTestRunner() #verbosity=2)
559	runner.run(suite)
560

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

Download in other formats: