Context Navigation

source: anuga_core/source/anuga/culvert_flows/culvert_class.py @ 6111

Last change on this file since 6111 was 6111, checked in by ole, 15 years ago
Work on culvert testing. New test disabled due to negative total energy issue in Boyd routine.
File size: 34.4 KB

Line
1	from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
2	from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
3	from anuga.utilities.system_tools import log_to_file
4	from anuga.utilities.polygon import inside_polygon
5	from anuga.utilities.polygon import is_inside_polygon
6	from anuga.utilities.polygon import plot_polygons
7
8	from anuga.utilities.numerical_tools import ensure_numeric
9	from Numeric import allclose
10
11	import sys
12
13	class Below_interval(Exception): pass
14	class Above_interval(Exception): pass
15
16
17
18	# FIXME(Ole): Write in C and reuse this function by similar code in interpolate.py
19	def interpolate_linearly(x, xvec, yvec):
20
21	# Check bounds
22	if x < xvec[0]:
23	msg = 'Value provided = %.2f, interpolation minimum = %.2f.' %(x, xvec[0])
24	raise Below_interval, msg
25
26	if x > xvec[-1]:
27	msg = 'Value provided = %.2f, interpolation maximum = %.2f.' %(x, xvec[-1])
28	raise Above_interval, msg
29
30
31	# Find appropriate slot within bounds
32	i = 0
33	while x > xvec[i]: i += 1
34
35
36	x0 = xvec[i-1]
37	x1 = xvec[i]
38	alpha = (x - x0)/(x1 - x0)
39
40	y0 = yvec[i-1]
41	y1 = yvec[i]
42	y = alphay1 + (1-alpha)y0
43
44	return y
45
46
47
48	def read_culvert_description(culvert_description_filename):
49
50	# Read description file
51	fid = open(culvert_description_filename)
52
53	read_rating_curve_data = False
54	rating_curve = []
55	for i, line in enumerate(fid.readlines()):
56
57	if read_rating_curve_data is True:
58	fields = line.split(',')
59	head_difference = float(fields[0].strip())
60	flow_rate = float(fields[1].strip())
61	barrel_velocity = float(fields[2].strip())
62
63	rating_curve.append( [head_difference, flow_rate, barrel_velocity] )
64
65	if i == 0:
66	# Header
67	continue
68	if i == 1:
69	# Metadata
70	fields = line.split(',')
71	label=fields[0].strip()
72	type=fields[1].strip().lower()
73	assert type in ['box', 'pipe']
74
75	width=float(fields[2].strip())
76	height=float(fields[3].strip())
77	length=float(fields[4].strip())
78	number_of_barrels=int(fields[5].strip())
79	#fields[6] refers to losses
80	description=fields[7].strip()
81
82	if line.strip() == '': continue # Skip blanks
83
84	if line.startswith('Rating'):
85	read_rating_curve_data = True
86	# Flow data follows
87
88	fid.close()
89
90	return label, type, width, height, length, number_of_barrels, description, rating_curve
91
92
93	class Culvert_flow_rating:
94	"""Culvert flow - transfer water from one hole to another
95
96
97	Input: Two points, pipe_size (either diameter or width, height),
98	mannings_rougness,
99	inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
100	top-down_blockage_factor and bottom_up_blockage_factor
101
102	"""
103
104	def __init__(self,
105	domain,
106	culvert_description_filename=None,
107	end_point0=None,
108	end_point1=None,
109	update_interval=None,
110	log_file=False,
111	discharge_hydrograph=False,
112	verbose=False):
113
114	from Numeric import sqrt, sum
115
116
117	label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
118
119
120	# Store culvert information
121	self.label = label
122	self.description = description
123	self.culvert_type = type
124	self.rating_curve = ensure_numeric(rating_curve)
125	self.number_of_barrels = number_of_barrels
126
127	if label is None: label = 'culvert_flow'
128	label += '_' + str(id(self))
129	self.label = label
130
131	# File for storing discharge_hydrograph
132	if discharge_hydrograph is True:
133	self.timeseries_filename = label + '_timeseries.csv'
134	fid = open(self.timeseries_filename, 'w')
135	fid.write('time, discharge\n')
136	fid.close()
137
138	# Log file for storing general textual output
139	if log_file is True:
140	self.log_filename = label + '.log'
141	log_to_file(self.log_filename, self.label)
142	log_to_file(self.log_filename, description)
143	log_to_file(self.log_filename, self.culvert_type)
144
145
146	# Create the fundamental culvert polygons from POLYGON
147	#if self.culvert_type == 'circle':
148	# # Redefine width and height for use with create_culvert_polygons
149	# width = height = diameter
150
151	P = create_culvert_polygons(end_point0,
152	end_point1,
153	width=width,
154	height=height,
155	number_of_barrels=number_of_barrels)
156
157	if verbose is True:
158	pass
159	#plot_polygons([[end_point0, end_point1],
160	# P['exchange_polygon0'],
161	# P['exchange_polygon1'],
162	# P['enquiry_polygon0'],
163	# P['enquiry_polygon1']],
164	# figname='culvert_polygon_output')
165
166
167	# Compute the average point for enquiry
168	enquiry_point0 = sum(P['enquiry_polygon0'][:2])/2
169	enquiry_point1 = sum(P['enquiry_polygon1'][:2])/2
170
171	self.enquiry_points = [enquiry_point0, enquiry_point1]
172	self.enquiry_indices = []
173	for point in self.enquiry_points:
174	# Find nearest centroid
175	N = len(domain)
176	points = domain.get_centroid_coordinates(absolute=True)
177
178	# Calculate indices in exchange area for this forcing term
179
180	triangle_id = min_dist = sys.maxint
181	for k in range(N):
182	x, y = points[k,:] # Centroid
183
184	c = point
185	distance = (x-c[0])2+(y-c[1])2
186	if distance < min_dist:
187	min_dist = distance
188	triangle_id = k
189
190
191	if triangle_id < sys.maxint:
192	msg = 'found triangle with centroid (%f, %f)'\
193	%tuple(points[triangle_id, :])
194	msg += ' for point (%f, %f)' %tuple(point)
195
196	self.enquiry_indices.append(triangle_id)
197	else:
198	msg = 'Triangle not found for point (%f, %f)' %point
199	raise Exception, msg
200
201
202
203	# Check that all polygons lie within the mesh.
204	bounding_polygon = domain.get_boundary_polygon()
205	for key in P.keys():
206	if key in ['exchange_polygon0',
207	'exchange_polygon1']:
208	for point in list(P[key]) + self.enquiry_points:
209	msg = 'Point %s in polygon %s for culvert %s did not'\
210	%(str(point), key, self.label)
211	msg += 'fall within the domain boundary.'
212	assert is_inside_polygon(point, bounding_polygon), msg
213
214
215	# Create inflow object at each end of the culvert.
216	self.openings = []
217	self.openings.append(Inflow(domain,
218	polygon=P['exchange_polygon0']))
219
220	self.openings.append(Inflow(domain,
221	polygon=P['exchange_polygon1']))
222
223
224
225	dq = domain.quantities
226	for i, opening in enumerate(self.openings):
227	#elevation = dq['elevation'].get_values(location='centroids',
228	# interpolation_points=[self.enquiry_points[i]])
229
230	elevation = dq['elevation'].get_values(location='centroids',
231	indices=[self.enquiry_indices[i]])
232	opening.elevation = elevation
233	opening.stage = elevation # Simple assumption that culvert is dry initially
234
235	# Assume two openings for now: Referred to as 0 and 1
236	assert len(self.openings) == 2
237
238	# Determine pipe direction
239	self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation
240	if delta_z > 0.0:
241	self.inlet = self.openings[0]
242	self.outlet = self.openings[1]
243	else:
244	self.outlet = self.openings[0]
245	self.inlet = self.openings[1]
246
247
248	# Store basic geometry
249	self.end_points = [end_point0, end_point1]
250	self.vector = P['vector']
251	self.length = P['length']; assert self.length > 0.0
252	if not allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
253	msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length)
254	msg += ' does not match distance between specified'
255	msg += ' end points (%.2f m)' %self.length
256	print msg
257
258	self.verbose = verbose
259	self.last_update = 0.0 # For use with update_interval
260	self.update_interval = update_interval
261
262
263	# Print something
264	if hasattr(self, 'log_filename'):
265	s = 'Culvert Effective Length = %.2f m' %(self.length)
266	log_to_file(self.log_filename, s)
267
268	s = 'Culvert Direction is %s\n' %str(self.vector)
269	log_to_file(self.log_filename, s)
270
271
272
273
274
275	def __call__(self, domain):
276	from anuga.utilities.numerical_tools import mean
277
278	from anuga.config import g, epsilon
279	from Numeric import take, sqrt
280	from anuga.config import velocity_protection
281
282
283	# Time stuff
284	time = domain.get_time()
285
286
287	update = False
288	if self.update_interval is None:
289	update = True
290	else:
291	if time - self.last_update > self.update_interval or time == 0.0:
292	update = True
293
294
295	if update is True:
296	self.last_update = time
297	dq = domain.quantities
298
299	# Get average water depths at each opening
300	openings = self.openings # There are two Opening [0] and [1]
301	for i, opening in enumerate(openings):
302
303	# Compute mean values of selected quantitites in the
304	# enquiry area in front of the culvert
305	# Stage and velocity comes from enquiry area
306	# and elevation from exchange area
307
308	stage = dq['stage'].get_values(location='centroids',
309	indices=[self.enquiry_indices[i]])
310
311	# Store current average stage and depth with each opening object
312	opening.depth = stage - opening.elevation
313	opening.stage = stage
314
315
316
317	################# End of the FOR loop ################################################
318
319	# We now need to deal with each opening individually
320
321	# Determine flow direction based on total energy difference
322
323	delta_w = self.inlet.stage - self.outlet.stage
324
325	if hasattr(self, 'log_filename'):
326	s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time,
327	self.inlet.stage,
328	self.outlet.stage)
329	log_to_file(self.log_filename, s)
330
331
332	if self.inlet.depth <= 0.01:
333	Q = 0.0
334	else:
335	# Calculate discharge for one barrel and set inlet.rate and outlet.rate
336
337	try:
338	Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1])
339	except Below_interval, e:
340	Q = self.rating_curve[0,1]
341	msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
342	msg += str(e)
343	msg += '\n I will use minimum discharge %.2f m^3/s for culvert "%s"'\
344	%(Q, self.label)
345	if hasattr(self, 'log_filename'):
346	log_to_file(self.log_filename, msg)
347	except Above_interval, e:
348	Q = self.rating_curve[-1,1]
349	msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
350	msg += str(e)
351	msg += '\n I will use maximum discharge %.2f m^3/s for culvert "%s"'\
352	%(Q, self.label)
353	if hasattr(self, 'log_filename'):
354	log_to_file(self.log_filename, msg)
355
356
357
358
359	# Adjust discharge for multiple barrels
360	Q *= self.number_of_barrels
361
362
363	self.inlet.rate = -Q
364	self.outlet.rate = Q
365
366	# Log timeseries to file
367	try:
368	fid = open(self.timeseries_filename, 'a')
369	except:
370	pass
371	else:
372	fid.write('%.2f, %.2f\n' %(time, Q))
373	fid.close()
374
375
376
377
378
379	# Execute flow term for each opening
380	# This is where Inflow objects are evaluated using the last rate that has been calculated
381	#
382	# This will take place at every internal timestep and update the domain
383	for opening in self.openings:
384	opening(domain)
385
386
387
388
389
390
391	class Culvert_flow_energy:
392	"""Culvert flow - transfer water from one hole to another
393
394	Using Momentum as Calculated by Culvert Flow !!
395	Could be Several Methods Investigated to do This !!!
396
397	2008_May_08
398	To Ole:
399	OK so here we need to get the Polygon Creating code to create
400	polygons for the culvert Based on
401	the two input Points (X0,Y0) and (X1,Y1) - need to be passed
402	to create polygon
403
404	The two centers are now passed on to create_polygon.
405
406
407	Input: Two points, pipe_size (either diameter or width, height),
408	mannings_rougness,
409	inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
410	top-down_blockage_factor and bottom_up_blockage_factor
411
412
413	And the Delta H enquiry should be change from Openings in line 412
414	to the enquiry Polygons infront of the culvert
415	At the moment this script uses only Depth, later we can change it to
416	Energy...
417
418	Once we have Delta H can calculate a Flow Rate and from Flow Rate
419	an Outlet Velocity
420	The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
421
422	Invert levels are optional. If left out they will default to the
423	elevation at the opening.
424
425	"""
426
427	def __init__(self,
428	domain,
429	label=None,
430	description=None,
431	end_point0=None,
432	end_point1=None,
433	width=None,
434	height=None,
435	diameter=None,
436	manning=None, # Mannings Roughness for Culvert
437	invert_level0=None, # Invert level at opening 0
438	invert_level1=None, # Invert level at opening 1
439	loss_exit=None,
440	loss_entry=None,
441	loss_bend=None,
442	loss_special=None,
443	blockage_topdwn=None,
444	blockage_bottup=None,
445	culvert_routine=None,
446	number_of_barrels=1,
447	update_interval=None,
448	verbose=False):
449
450	from Numeric import sqrt, sum
451
452	# Input check
453	if diameter is not None:
454	self.culvert_type = 'circle'
455	self.diameter = diameter
456	if height is not None or width is not None:
457	msg = 'Either diameter or width&height must be specified, '
458	msg += 'but not both.'
459	raise Exception, msg
460	else:
461	if height is not None:
462	if width is None:
463	self.culvert_type = 'square'
464	width = height
465	else:
466	self.culvert_type = 'rectangle'
467	elif width is not None:
468	if height is None:
469	self.culvert_type = 'square'
470	height = width
471	else:
472	msg = 'Either diameter or width&height must be specified.'
473	raise Exception, msg
474
475	if height == width:
476	self.culvert_type = 'square'
477
478	self.height = height
479	self.width = width
480
481
482	assert self.culvert_type in ['circle', 'square', 'rectangle']
483
484	assert number_of_barrels >= 1
485	self.number_of_barrels = number_of_barrels
486
487
488	# Set defaults
489	if manning is None: manning = 0.012 # Default roughness for pipe
490	if loss_exit is None: loss_exit = 1.00
491	if loss_entry is None: loss_entry = 0.50
492	if loss_bend is None: loss_bend=0.00
493	if loss_special is None: loss_special=0.00
494	if blockage_topdwn is None: blockage_topdwn=0.00
495	if blockage_bottup is None: blockage_bottup=0.00
496	if culvert_routine is None:
497	culvert_routine=boyd_generalised_culvert_model
498
499	if label is None: label = 'culvert_flow'
500	label += '_' + str(id(self))
501	self.label = label
502
503	# File for storing culvert quantities
504	self.timeseries_filename = label + '_timeseries.csv'
505	fid = open(self.timeseries_filename, 'w')
506	fid.write('time, E0, E1, Velocity, Discharge\n')
507	fid.close()
508
509	# Log file for storing general textual output
510	self.log_filename = label + '.log'
511	log_to_file(self.log_filename, self.label)
512	log_to_file(self.log_filename, description)
513	log_to_file(self.log_filename, self.culvert_type)
514
515
516	# Create the fundamental culvert polygons from POLYGON
517	if self.culvert_type == 'circle':
518	# Redefine width and height for use with create_culvert_polygons
519	width = height = diameter
520
521	P = create_culvert_polygons(end_point0,
522	end_point1,
523	width=width,
524	height=height,
525	number_of_barrels=number_of_barrels)
526
527	if verbose is True:
528	pass
529	#plot_polygons([[end_point0, end_point1],
530	# P['exchange_polygon0'],
531	# P['exchange_polygon1'],
532	# P['enquiry_polygon0'],
533	# P['enquiry_polygon1']],
534	# figname='culvert_polygon_output')
535	#import sys; sys.exit()
536
537
538	# Compute the average point for enquiry
539	enquiry_point0 = sum(P['enquiry_polygon0'][:2])/2
540	enquiry_point1 = sum(P['enquiry_polygon1'][:2])/2
541
542	self.enquiry_points = [enquiry_point0, enquiry_point1]
543	self.enquiry_indices = []
544	for point in self.enquiry_points:
545	# Find nearest centroid
546	N = len(domain)
547	points = domain.get_centroid_coordinates(absolute=True)
548
549	# Calculate indices in exchange area for this forcing term
550
551	triangle_id = min_dist = sys.maxint
552	for k in range(N):
553	x, y = points[k,:] # Centroid
554
555	c = point
556	distance = (x-c[0])2+(y-c[1])2
557	if distance < min_dist:
558	min_dist = distance
559	triangle_id = k
560
561
562	if triangle_id < sys.maxint:
563	msg = 'found triangle with centroid (%f, %f)'\
564	%tuple(points[triangle_id, :])
565	msg += ' for point (%f, %f)' %tuple(point)
566
567	self.enquiry_indices.append(triangle_id)
568	else:
569	msg = 'Triangle not found for point (%f, %f)' %point
570	raise Exception, msg
571
572
573
574
575
576
577	# Check that all polygons lie within the mesh.
578	bounding_polygon = domain.get_boundary_polygon()
579	for key in P.keys():
580	if key in ['exchange_polygon0',
581	'exchange_polygon1',
582	'enquiry_polygon0',
583	'enquiry_polygon1']:
584	for point in P[key]:
585
586	msg = 'Point %s in polygon %s for culvert %s did not'\
587	%(str(point), key, self.label)
588	msg += 'fall within the domain boundary.'
589	assert is_inside_polygon(point, bounding_polygon), msg
590
591
592	# Create inflow object at each end of the culvert.
593	self.openings = []
594	self.openings.append(Inflow(domain,
595	polygon=P['exchange_polygon0']))
596
597	self.openings.append(Inflow(domain,
598	polygon=P['exchange_polygon1']))
599
600
601	# Assume two openings for now: Referred to as 0 and 1
602	assert len(self.openings) == 2
603
604	# Store basic geometry
605	self.end_points = [end_point0, end_point1]
606	self.invert_levels = [invert_level0, invert_level1]
607	#self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
608	self.enquiry_polylines = [P['enquiry_polygon0'][:2],
609	P['enquiry_polygon1'][:2]]
610	self.vector = P['vector']
611	self.length = P['length']; assert self.length > 0.0
612	self.verbose = verbose
613	self.last_time = 0.0
614	self.last_update = 0.0 # For use with update_interval
615	self.update_interval = update_interval
616
617
618	# Store hydraulic parameters
619	self.manning = manning
620	self.loss_exit = loss_exit
621	self.loss_entry = loss_entry
622	self.loss_bend = loss_bend
623	self.loss_special = loss_special
624	self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
625	self.blockage_topdwn = blockage_topdwn
626	self.blockage_bottup = blockage_bottup
627
628	# Store culvert routine
629	self.culvert_routine = culvert_routine
630
631
632	# Create objects to update momentum (a bit crude at this stage)
633
634
635	xmom0 = General_forcing(domain, 'xmomentum',
636	polygon=P['exchange_polygon0'])
637
638	xmom1 = General_forcing(domain, 'xmomentum',
639	polygon=P['exchange_polygon1'])
640
641	ymom0 = General_forcing(domain, 'ymomentum',
642	polygon=P['exchange_polygon0'])
643
644	ymom1 = General_forcing(domain, 'ymomentum',
645	polygon=P['exchange_polygon1'])
646
647	self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
648
649
650	# Print something
651	s = 'Culvert Effective Length = %.2f m' %(self.length)
652	log_to_file(self.log_filename, s)
653
654	s = 'Culvert Direction is %s\n' %str(self.vector)
655	log_to_file(self.log_filename, s)
656
657
658	def __call__(self, domain):
659	from anuga.utilities.numerical_tools import mean
660
661	from anuga.config import g, epsilon
662	from Numeric import take, sqrt
663	from anuga.config import velocity_protection
664
665
666	log_filename = self.log_filename
667
668	# Time stuff
669	time = domain.get_time()
670
671
672	update = False
673	if self.update_interval is None:
674	update = True
675	else:
676	if time - self.last_update > self.update_interval or time == 0.0:
677	update = True
678
679	#print 'call', time, time - self.last_update
680
681
682	if update is True:
683	#print 'Updating', time, time - self.last_update
684	self.last_update = time
685
686	delta_t = time-self.last_time
687	s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
688	log_to_file(log_filename, s)
689
690	msg = 'Time did not advance'
691	if time > 0.0: assert delta_t > 0.0, msg
692
693
694	# Get average water depths at each opening
695	openings = self.openings # There are two Opening [0] and [1]
696	for i, opening in enumerate(openings):
697	dq = domain.quantities
698
699	# Compute mean values of selected quantitites in the
700	# exchange area in front of the culvert
701	# Stage and velocity comes from enquiry area
702	# and elevation from exchange area
703
704	stage = dq['stage'].get_values(location='centroids',
705	indices=opening.exchange_indices)
706	w = mean(stage) # Average stage
707
708	# Use invert level instead of elevation if specified
709	invert_level = self.invert_levels[i]
710	if invert_level is not None:
711	z = invert_level
712	else:
713	elevation = dq['elevation'].get_values(location='centroids',
714	indices=opening.exchange_indices)
715	z = mean(elevation) # Average elevation
716
717	# Estimated depth above the culvert inlet
718	d = w - z # Used for calculations involving depth
719	if d < 0.0:
720	# This is possible since w and z are taken at different locations
721	#msg = 'D < 0.0: %f' %d
722	#raise msg
723	d = 0.0
724
725
726	# Ratio of depth to culvert height.
727	# If ratio > 1 then culvert is running full
728	if self.culvert_type == 'circle':
729	ratio = d/self.diameter
730	else:
731	ratio = d/self.height
732	opening.ratio = ratio
733
734
735	# Average measures of energy in front of this opening
736	#polyline = self.enquiry_polylines[i]
737	#opening.total_energy = domain.get_energy_through_cross_section(polyline,
738	# kind='total')
739
740	id = [self.enquiry_indices[i]]
741	stage = dq['stage'].get_values(location='centroids',
742	indices=id)
743	elevation = dq['elevation'].get_values(location='centroids',
744	indices=id)
745	xmomentum = dq['xmomentum'].get_values(location='centroids',
746	indices=id)
747	ymomentum = dq['xmomentum'].get_values(location='centroids',
748	indices=id)
749	depth = stage-elevation
750	if depth > 0.0:
751	u = xmomentum/(depth + velocity_protection/depth)
752	v = ymomentum/(depth + velocity_protection/depth)
753	else:
754	u = v = 0.0
755
756
757	opening.total_energy = 0.5(uu + v*v)/g + stage
758	#print 'Et = %.3f m' %opening.total_energy
759
760	# Store current average stage and depth with each opening object
761	opening.depth = d
762	opening.depth_trigger = d # Use this for now
763	opening.stage = w
764	opening.elevation = z
765
766
767	################# End of the FOR loop ################################################
768
769	# We now need to deal with each opening individually
770
771	# Determine flow direction based on total energy difference
772	delta_Et = openings[0].total_energy - openings[1].total_energy
773
774	if delta_Et > 0:
775	#print 'Flow U/S ---> D/S'
776	inlet = openings[0]
777	outlet = openings[1]
778
779	inlet.momentum = self.opening_momentum[0]
780	outlet.momentum = self.opening_momentum[1]
781
782	else:
783	#print 'Flow D/S ---> U/S'
784	inlet = openings[1]
785	outlet = openings[0]
786
787	inlet.momentum = self.opening_momentum[1]
788	outlet.momentum = self.opening_momentum[0]
789
790	delta_Et = -delta_Et
791
792	self.inlet = inlet
793	self.outlet = outlet
794
795	msg = 'Total energy difference is negative'
796	assert delta_Et > 0.0, msg
797
798	delta_h = inlet.stage - outlet.stage
799	delta_z = inlet.elevation - outlet.elevation
800	culvert_slope = (delta_z/self.length)
801
802	if culvert_slope < 0.0:
803	# Adverse gradient - flow is running uphill
804	# Flow will be purely controlled by uphill outlet face
805	if self.verbose is True:
806	print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
807
808
809	s = 'Delta total energy = %.3f' %(delta_Et)
810	log_to_file(log_filename, s)
811
812
813	# Calculate discharge for one barrel and set inlet.rate and outlet.rate
814	Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
815
816	# Adjust discharge for multiple barrels
817	Q *= self.number_of_barrels
818
819	# Compute barrel momentum
820	barrel_momentum = barrel_velocity*culvert_outlet_depth
821
822	s = 'Barrel velocity = %f' %barrel_velocity
823	log_to_file(log_filename, s)
824
825	# Compute momentum vector at outlet
826	outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
827
828	s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
829	log_to_file(log_filename, s)
830
831	# Log timeseries to file
832	fid = open(self.timeseries_filename, 'a')
833	fid.write('%f, %f, %f, %f, %f\n'\
834	%(time,
835	openings[0].total_energy,
836	openings[1].total_energy,
837	barrel_velocity,
838	Q))
839	fid.close()
840
841	# Update momentum
842	delta_t = time - self.last_time
843	if delta_t > 0.0:
844	xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
845	xmomentum_rate /= delta_t
846
847	ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
848	ymomentum_rate /= delta_t
849
850	s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
851	log_to_file(log_filename, s)
852	else:
853	xmomentum_rate = ymomentum_rate = 0.0
854
855
856	# Set momentum rates for outlet jet
857	outlet.momentum[0].rate = xmomentum_rate
858	outlet.momentum[1].rate = ymomentum_rate
859
860	# Remember this value for next step (IMPORTANT)
861	outlet.momentum[0].value = outlet_mom_x
862	outlet.momentum[1].value = outlet_mom_y
863
864	if int(domain.time*100) % 100 == 0:
865	s = 'T=%.5f, Culvert Discharge = %.3f f'\
866	%(time, Q)
867	s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\
868	%(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
869	s += ' Momentum rate: (%.4f, %.4f)'\
870	%(xmomentum_rate, ymomentum_rate)
871	s+='Outlet Vel= %.3f'\
872	%(barrel_velocity)
873	log_to_file(log_filename, s)
874
875
876
877
878	# Execute flow term for each opening
879	# This is where Inflow objects are evaluated and update the domain
880	for opening in self.openings:
881	opening(domain)
882
883	# Execute momentum terms
884	# This is where Inflow objects are evaluated and update the domain
885	self.outlet.momentum[0](domain)
886	self.outlet.momentum[1](domain)
887
888	# Store value of time #FIXME(Ole): Maybe only every time we update
889	self.last_time = time
890
891
892	Culvert_flow = Culvert_flow_rating

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

Download in other formats: