Context Navigation

source: branches/numpy/anuga/culvert_flows/culvert_class.py @ 6982

Last change on this file since 6982 was 6902, checked in by rwilson, 16 years ago
Back-merge from Numeric trunk to numpy branch.
File size: 53.6 KB

Line
1	import sys
2
3	from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
4	from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
5	from anuga.utilities.system_tools import log_to_file
6	from anuga.utilities.polygon import inside_polygon
7	from anuga.utilities.polygon import is_inside_polygon
8	from anuga.utilities.polygon import plot_polygons
9
10	from anuga.utilities.numerical_tools import mean
11	from anuga.utilities.numerical_tools import ensure_numeric, sign
12
13	from anuga.config import g, epsilon
14	from anuga.config import minimum_allowed_height, velocity_protection
15
16	import numpy as num
17
18
19	class Below_interval(Exception): pass
20	class Above_interval(Exception): pass
21
22	# FIXME(Ole): Take a good hard look at logging here
23
24
25	# FIXME(Ole): Write in C and reuse this function by similar code
26	# in interpolate.py
27	def interpolate_linearly(x, xvec, yvec):
28
29	msg = 'Input to function interpolate_linearly could not be converted '
30	msg += 'to numerical scalar: x = %s' % str(x)
31	try:
32	x = float(x)
33	except:
34	raise Exception, msg
35
36
37	# Check bounds
38	if x < xvec[0]:
39	msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
40	% (x, xvec[0])
41	raise Below_interval, msg
42
43	if x > xvec[-1]:
44	msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
45	%(x, xvec[-1])
46	raise Above_interval, msg
47
48
49	# Find appropriate slot within bounds
50	i = 0
51	while x > xvec[i]: i += 1
52
53
54	x0 = xvec[i-1]
55	x1 = xvec[i]
56	alpha = (x - x0)/(x1 - x0)
57
58	y0 = yvec[i-1]
59	y1 = yvec[i]
60	y = alphay1 + (1-alpha)y0
61
62	return y
63
64
65
66	def read_culvert_description(culvert_description_filename):
67
68	# Read description file
69	fid = open(culvert_description_filename)
70
71	read_rating_curve_data = False
72	rating_curve = []
73	for i, line in enumerate(fid.readlines()):
74
75	if read_rating_curve_data is True:
76	fields = line.split(',')
77	head_difference = float(fields[0].strip())
78	flow_rate = float(fields[1].strip())
79	barrel_velocity = float(fields[2].strip())
80
81	rating_curve.append([head_difference, flow_rate, barrel_velocity])
82
83	if i == 0:
84	# Header
85	continue
86	if i == 1:
87	# Metadata
88	fields = line.split(',')
89	label=fields[0].strip()
90	type=fields[1].strip().lower()
91	assert type in ['box', 'pipe']
92
93	width=float(fields[2].strip())
94	height=float(fields[3].strip())
95	length=float(fields[4].strip())
96	number_of_barrels=int(fields[5].strip())
97	#fields[6] refers to losses
98	description=fields[7].strip()
99
100	if line.strip() == '': continue # Skip blanks
101
102	if line.startswith('Rating'):
103	read_rating_curve_data = True
104	# Flow data follows
105
106	fid.close()
107
108	return label, type, width, height, length, number_of_barrels, description, rating_curve
109
110
111
112
113	class Culvert_flow_general:
114	"""Culvert flow - transfer water from one hole to another
115
116	This version will work with either rating curve file or with culvert
117	routine.
118
119	Input: Two points, pipe_size (either diameter or width, height),
120	mannings_rougness,
121	"""
122
123	def __init__(self,
124	domain,
125	culvert_description_filename=None,
126	culvert_routine=None,
127	end_point0=None,
128	end_point1=None,
129	enquiry_point0=None,
130	enquiry_point1=None,
131	type='box',
132	width=None,
133	height=None,
134	length=None,
135	number_of_barrels=1,
136	trigger_depth=0.01, # Depth below which no flow happens
137	manning=None, # Mannings Roughness for Culvert
138	sum_loss=None,
139	use_velocity_head=False,
140	use_momentum_jet=False, # FIXME(Ole): Not yet implemented
141	label=None,
142	description=None,
143	update_interval=None,
144	log_file=False,
145	discharge_hydrograph=False,
146	verbose=False):
147
148
149
150	# Input check
151
152	if height is None: height = width
153	self.height = height
154	self.width = width
155
156
157	assert number_of_barrels >= 1
158	assert use_velocity_head is True or use_velocity_head is False
159
160	msg = 'Momentum jet not yet moved to general culvert'
161	assert use_momentum_jet is False, msg
162
163	self.culvert_routine = culvert_routine
164	self.culvert_description_filename = culvert_description_filename
165	if culvert_description_filename is not None:
166	label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
167	self.rating_curve = ensure_numeric(rating_curve)
168
169	self.domain = domain
170	self.trigger_depth = trigger_depth
171
172	if manning is None:
173	self.manning = 0.012 # Default roughness for pipe
174
175	if sum_loss is None:
176	self.sum_loss = 0.0
177
178
179
180	# Store culvert information
181	self.label = label
182	self.description = description
183	self.culvert_type = type
184	self.number_of_barrels = number_of_barrels
185
186	# Store options
187	self.use_velocity_head = use_velocity_head
188
189	if label is None: label = 'culvert_flow'
190	label += '_' + str(id(self))
191	self.label = label
192
193	# File for storing discharge_hydrograph
194	if discharge_hydrograph is True:
195	self.timeseries_filename = label + '_timeseries.csv'
196	fid = open(self.timeseries_filename, 'w')
197	fid.write('time, discharge\n')
198	fid.close()
199
200	# Log file for storing general textual output
201	if log_file is True:
202	self.log_filename = label + '.log'
203	log_to_file(self.log_filename, self.label)
204	log_to_file(self.log_filename, description)
205	log_to_file(self.log_filename, self.culvert_type)
206	else:
207	self.log_filename = None
208
209
210	# Create the fundamental culvert polygons from polygon
211	P = create_culvert_polygons(end_point0,
212	end_point1,
213	width=width,
214	height=height,
215	number_of_barrels=number_of_barrels)
216	self.culvert_polygons = P
217
218	# Select enquiry points
219	if enquiry_point0 is None:
220	enquiry_point0 = P['enquiry_point0']
221
222	if enquiry_point1 is None:
223	enquiry_point1 = P['enquiry_point1']
224
225	if verbose is True:
226	pass
227	#plot_polygons([[end_point0, end_point1],
228	# P['exchange_polygon0'],
229	# P['exchange_polygon1'],
230	# [enquiry_point0, 1.005*enquiry_point0],
231	# [enquiry_point1, 1.005*enquiry_point1]],
232	# figname='culvert_polygon_output')
233
234
235
236	self.enquiry_points = [enquiry_point0, enquiry_point1]
237	self.enquiry_indices = self.get_enquiry_indices()
238	self.check_culvert_inside_domain()
239
240
241	# Create inflow object at each end of the culvert.
242	self.openings = []
243	self.openings.append(Inflow(domain,
244	polygon=P['exchange_polygon0']))
245	self.openings.append(Inflow(domain,
246	polygon=P['exchange_polygon1']))
247
248	# Assume two openings for now: Referred to as 0 and 1
249	assert len(self.openings) == 2
250
251	# Establish initial values at each enquiry point
252	dq = domain.quantities
253	for i, opening in enumerate(self.openings):
254	idx = self.enquiry_indices[i]
255	elevation = dq['elevation'].get_values(location='centroids',
256	indices=[idx])[0]
257	stage = dq['stage'].get_values(location='centroids',
258	indices=[idx])[0]
259	opening.elevation = elevation
260	opening.stage = stage
261	opening.depth = stage-elevation
262
263
264
265	# Determine initial pipe direction.
266	# This may change dynamically based on the total energy difference
267	# Consequently, this may be superfluous
268	delta_z = self.openings[0].elevation - self.openings[1].elevation
269	if delta_z > 0.0:
270	self.inlet = self.openings[0]
271	self.outlet = self.openings[1]
272	else:
273	self.outlet = self.openings[0]
274	self.inlet = self.openings[1]
275
276
277	# Store basic geometry
278	self.end_points = [end_point0, end_point1]
279	self.vector = P['vector']
280	self.length = P['length']; assert self.length > 0.0
281	if culvert_description_filename is not None:
282	if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
283	msg = 'WARNING: barrel length specified in "%s" (%.2f m)'\
284	% (culvert_description_filename,
285	length)
286	msg += ' does not match distance between specified'
287	msg += ' end points (%.2f m)' %self.length
288	print msg
289
290	self.verbose = verbose
291
292
293
294	# For use with update_interval
295	self.last_update = 0.0
296	self.update_interval = update_interval
297
298
299	# Print some diagnostics to log if requested
300	if self.log_filename is not None:
301	s = 'Culvert Effective Length = %.2f m' %(self.length)
302	log_to_file(self.log_filename, s)
303
304	s = 'Culvert Direction is %s\n' %str(self.vector)
305	log_to_file(self.log_filename, s)
306
307
308
309
310
311	def __call__(self, domain):
312
313	# Time stuff
314	time = domain.get_time()
315
316
317	update = False
318	if self.update_interval is None:
319	# Use next timestep as has been computed in domain.py
320	delta_t = domain.timestep
321	update = True
322	else:
323	# Use update interval
324	delta_t = self.update_interval
325	if time - self.last_update > self.update_interval or time == 0.0:
326	update = True
327
328	if self.log_filename is not None:
329	s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
330	log_to_file(self.log_filename, s)
331
332
333	if update is True:
334	self.compute_rates(delta_t)
335
336
337	# Execute flow term for each opening
338	# This is where Inflow objects are evaluated using the last rate
339	# that has been calculated
340	#
341	# This will take place at every internal timestep and update the domain
342	for opening in self.openings:
343	opening(domain)
344
345
346
347	def get_enquiry_indices(self):
348	"""Get indices for nearest centroids to self.enquiry_points
349	"""
350
351	domain = self.domain
352
353	enquiry_indices = []
354	for point in self.enquiry_points:
355	# Find nearest centroid
356	N = len(domain)
357	points = domain.get_centroid_coordinates(absolute=True)
358
359	# Calculate indices in exchange area for this forcing term
360
361	triangle_id = min_dist = sys.maxint
362	for k in range(N):
363	x, y = points[k,:] # Centroid
364
365	c = point
366	distance = (x-c[0])2+(y-c[1])2
367	if distance < min_dist:
368	min_dist = distance
369	triangle_id = k
370
371
372	if triangle_id < sys.maxint:
373	msg = 'found triangle with centroid (%f, %f)'\
374	%tuple(points[triangle_id, :])
375	msg += ' for point (%f, %f)' %tuple(point)
376
377	enquiry_indices.append(triangle_id)
378	else:
379	msg = 'Triangle not found for point (%f, %f)' %point
380	raise Exception, msg
381
382	return enquiry_indices
383
384
385	def check_culvert_inside_domain(self):
386	"""Check that all polygons and enquiry points lie within the mesh.
387	"""
388	bounding_polygon = self.domain.get_boundary_polygon()
389	P = self.culvert_polygons
390	for key in P.keys():
391	if key in ['exchange_polygon0',
392	'exchange_polygon1']:
393	for point in list(P[key]) + self.enquiry_points:
394	msg = 'Point %s in polygon %s for culvert %s did not'\
395	%(str(point), key, self.label)
396	msg += 'fall within the domain boundary.'
397	assert is_inside_polygon(point, bounding_polygon), msg
398
399
400	def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
401	"""Adjust Q downwards depending on available water at inlet
402	"""
403
404	if delta_t < epsilon:
405	# No need to adjust if time step is very small or zero
406	# In this case the possible flow will be very large
407	# anyway.
408	return Q
409
410	# Short hands
411	domain = self.domain
412	dq = domain.quantities
413	time = domain.get_time()
414	I = self.inlet
415	idx = I.exchange_indices
416
417	# Find triangle with the smallest depth
418	stage = dq['stage'].get_values(location='centroids',
419	indices=[idx])
420	elevation = dq['elevation'].get_values(location='centroids',
421	indices=[idx])
422	depth = stage-elevation
423	min_depth = min(depth.flat)
424
425	# Compute possible flow for exchange region based on
426	# triangle with smallest depth
427	max_Q = min_depth*I.exchange_area/delta_t
428
429	# Calculate the minimum in absolute terms of
430	# the requsted flow and the possible flow
431	Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
432
433	if abs(Q_reduced) < abs(Q):
434	msg = '%.2fs: Requested flow is ' % time
435	msg += 'greater than what is supported by the smallest '
436	msg += 'depth at inlet exchange area:\n '
437	msg += 'h_mininlet_area/delta_t = %.2f%.2f/%.2f '\
438	% (min_depth, I.exchange_area, delta_t)
439	msg += ' = %.2f m^3/s\n ' % Q_reduced
440	msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
441	if self.verbose is True:
442	print msg
443
444	if self.log_filename is not None:
445	log_to_file(self.log_filename, msg)
446
447	return Q_reduced
448
449
450	def compute_rates(self, delta_t):
451	"""Compute new rates for inlet and outlet
452	"""
453
454	# Short hands
455	domain = self.domain
456	dq = domain.quantities
457
458	# Time stuff
459	time = domain.get_time()
460	self.last_update = time
461
462
463	if hasattr(self, 'log_filename'):
464	log_filename = self.log_filename
465
466	# Compute stage and energy at the
467	# enquiry points at each end of the culvert
468	openings = self.openings
469	for i, opening in enumerate(openings):
470	idx = self.enquiry_indices[i]
471
472	stage = dq['stage'].get_values(location='centroids',
473	indices=[idx])[0]
474	depth = h = stage-opening.elevation
475
476
477	if self.use_velocity_head is True:
478	xmomentum = dq['xmomentum'].get_values(location='centroids',
479	indices=[idx])[0]
480	ymomentum = dq['xmomentum'].get_values(location='centroids',
481	indices=[idx])[0]
482
483	if h > minimum_allowed_height:
484	u = xmomentum/(h + velocity_protection/h)
485	v = ymomentum/(h + velocity_protection/h)
486	else:
487	u = v = 0.0
488
489	velocity_head = 0.5(uu + v*v)/g
490	else:
491	velocity_head = 0.0
492
493	opening.total_energy = velocity_head + stage
494	opening.specific_energy = velocity_head + depth
495	opening.stage = stage
496	opening.depth = depth
497
498
499	# We now need to deal with each opening individually
500	# Determine flow direction based on total energy difference
501	delta_total_energy = openings[0].total_energy - openings[1].total_energy
502	if delta_total_energy > 0:
503	#print 'Flow U/S ---> D/S'
504	inlet = openings[0]
505	outlet = openings[1]
506	else:
507	#print 'Flow D/S ---> U/S'
508	inlet = openings[1]
509	outlet = openings[0]
510
511	delta_total_energy = -delta_total_energy
512
513	self.inlet = inlet
514	self.outlet = outlet
515
516	msg = 'Total energy difference is negative'
517	assert delta_total_energy > 0.0, msg
518
519	# Recompute slope and issue warning if flow is uphill
520	# These values do not enter the computation
521	delta_z = inlet.elevation - outlet.elevation
522	culvert_slope = (delta_z/self.length)
523	if culvert_slope < 0.0:
524	# Adverse gradient - flow is running uphill
525	# Flow will be purely controlled by uphill outlet face
526	if self.verbose is True:
527	print '%.2fs - WARNING: Flow is running uphill.' % time
528
529	if self.log_filename is not None:
530	s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
531	%(time, self.inlet.stage, self.outlet.stage)
532	log_to_file(self.log_filename, s)
533	s = 'Delta total energy = %.3f' %(delta_total_energy)
534	log_to_file(log_filename, s)
535
536
537	# Determine controlling energy (driving head) for culvert
538	if inlet.specific_energy > delta_total_energy:
539	# Outlet control
540	driving_head = delta_total_energy
541	else:
542	# Inlet control
543	driving_head = inlet.specific_energy
544
545
546
547	if self.inlet.depth <= self.trigger_depth:
548	Q = 0.0
549	else:
550	# Calculate discharge for one barrel and
551	# set inlet.rate and outlet.rate
552
553	if self.culvert_description_filename is not None:
554	try:
555	Q = interpolate_linearly(driving_head,
556	self.rating_curve[:,0],
557	self.rating_curve[:,1])
558	except Below_interval, e:
559	Q = self.rating_curve[0,1]
560	msg = '%.2fs: ' % time
561	msg += 'Delta head smaller than rating curve minimum: '
562	msg += str(e)
563	msg += '\n '
564	msg += 'I will use minimum discharge %.2f m^3/s ' % Q
565	msg += 'for culvert "%s"' % self.label
566
567	if hasattr(self, 'log_filename'):
568	log_to_file(self.log_filename, msg)
569	except Above_interval, e:
570	Q = self.rating_curve[-1,1]
571	msg = '%.2fs: ' % time
572	msg += 'Delta head greater than rating curve maximum: '
573	msg += str(e)
574	msg += '\n '
575	msg += 'I will use maximum discharge %.2f m^3/s ' % Q
576	msg += 'for culvert "%s"' % self.label
577
578	if self.log_filename is not None:
579	log_to_file(self.log_filename, msg)
580	else:
581	# User culvert routine
582	Q, barrel_velocity, culvert_outlet_depth =\
583	self.culvert_routine(inlet.depth,
584	outlet.depth,
585	inlet.specific_energy,
586	delta_total_energy,
587	g,
588	culvert_length=self.length,
589	culvert_width=self.width,
590	culvert_height=self.height,
591	culvert_type=self.culvert_type,
592	manning=self.manning,
593	sum_loss=self.sum_loss,
594	log_filename=self.log_filename)
595
596
597
598	# Adjust discharge for multiple barrels
599	Q *= self.number_of_barrels
600
601
602	Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
603
604	self.inlet.rate = -Q
605	self.outlet.rate = Q
606
607	# Log timeseries to file
608	try:
609	fid = open(self.timeseries_filename, 'a')
610	except:
611	pass
612	else:
613	fid.write('%.2f, %.2f\n' %(time, Q))
614	fid.close()
615
616
617	# OBSOLETE (Except for momentum jet in Culvert_flow_energy)
618	class Culvert_flow_rating:
619	"""Culvert flow - transfer water from one hole to another
620
621
622	Input: Two points, pipe_size (either diameter or width, height),
623	mannings_rougness,
624	inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
625	top-down_blockage_factor and bottom_up_blockage_factor
626
627	"""
628
629	def __init__(self,
630	domain,
631	culvert_description_filename=None,
632	end_point0=None,
633	end_point1=None,
634	enquiry_point0=None,
635	enquiry_point1=None,
636	update_interval=None,
637	log_file=False,
638	discharge_hydrograph=False,
639	verbose=False):
640
641
642
643	label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
644
645
646	# Store culvert information
647	self.label = label
648	self.description = description
649	self.culvert_type = type
650	self.rating_curve = ensure_numeric(rating_curve)
651	self.number_of_barrels = number_of_barrels
652
653	if label is None: label = 'culvert_flow'
654	label += '_' + str(id(self))
655	self.label = label
656
657	# File for storing discharge_hydrograph
658	if discharge_hydrograph is True:
659	self.timeseries_filename = label + '_timeseries.csv'
660	fid = open(self.timeseries_filename, 'w')
661	fid.write('time, discharge\n')
662	fid.close()
663
664	# Log file for storing general textual output
665	if log_file is True:
666	self.log_filename = label + '.log'
667	log_to_file(self.log_filename, self.label)
668	log_to_file(self.log_filename, description)
669	log_to_file(self.log_filename, self.culvert_type)
670
671
672	# Create the fundamental culvert polygons from POLYGON
673	#if self.culvert_type == 'circle':
674	# # Redefine width and height for use with create_culvert_polygons
675	# width = height = diameter
676
677	P = create_culvert_polygons(end_point0,
678	end_point1,
679	width=width,
680	height=height,
681	number_of_barrels=number_of_barrels)
682
683	# Select enquiry points
684	if enquiry_point0 is None:
685	enquiry_point0 = P['enquiry_point0']
686
687	if enquiry_point1 is None:
688	enquiry_point1 = P['enquiry_point1']
689
690	if verbose is True:
691	pass
692	#plot_polygons([[end_point0, end_point1],
693	# P['exchange_polygon0'],
694	# P['exchange_polygon1'],
695	# [enquiry_point0, 1.005*enquiry_point0],
696	# [enquiry_point1, 1.005*enquiry_point1]],
697	# figname='culvert_polygon_output')
698
699
700
701	self.enquiry_points = [enquiry_point0, enquiry_point1]
702
703	self.enquiry_indices = []
704	for point in self.enquiry_points:
705	# Find nearest centroid
706	N = len(domain)
707	points = domain.get_centroid_coordinates(absolute=True)
708
709	# Calculate indices in exchange area for this forcing term
710
711	triangle_id = min_dist = sys.maxint
712	for k in range(N):
713	x, y = points[k,:] # Centroid
714
715	c = point
716	distance = (x-c[0])2+(y-c[1])2
717	if distance < min_dist:
718	min_dist = distance
719	triangle_id = k
720
721
722	if triangle_id < sys.maxint:
723	msg = 'found triangle with centroid (%f, %f)'\
724	%tuple(points[triangle_id, :])
725	msg += ' for point (%f, %f)' %tuple(point)
726
727	self.enquiry_indices.append(triangle_id)
728	else:
729	msg = 'Triangle not found for point (%f, %f)' %point
730	raise Exception, msg
731
732
733
734	# Check that all polygons lie within the mesh.
735	bounding_polygon = domain.get_boundary_polygon()
736	for key in P.keys():
737	if key in ['exchange_polygon0',
738	'exchange_polygon1']:
739	for point in list(P[key]) + self.enquiry_points:
740	msg = 'Point %s in polygon %s for culvert %s did not'\
741	%(str(point), key, self.label)
742	msg += 'fall within the domain boundary.'
743	assert is_inside_polygon(point, bounding_polygon), msg
744
745
746	# Create inflow object at each end of the culvert.
747	self.openings = []
748	self.openings.append(Inflow(domain,
749	polygon=P['exchange_polygon0']))
750
751	self.openings.append(Inflow(domain,
752	polygon=P['exchange_polygon1']))
753
754
755
756	dq = domain.quantities
757	for i, opening in enumerate(self.openings):
758	elevation = dq['elevation'].get_values(location='centroids',
759	indices=[self.enquiry_indices[i]])
760	opening.elevation = elevation
761	opening.stage = elevation # Simple assumption that culvert is dry initially
762
763	# Assume two openings for now: Referred to as 0 and 1
764	assert len(self.openings) == 2
765
766	# Determine pipe direction
767	self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation
768	if delta_z > 0.0:
769	self.inlet = self.openings[0]
770	self.outlet = self.openings[1]
771	else:
772	self.outlet = self.openings[0]
773	self.inlet = self.openings[1]
774
775
776	# Store basic geometry
777	self.end_points = [end_point0, end_point1]
778	self.vector = P['vector']
779	self.length = P['length']; assert self.length > 0.0
780	if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
781	msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length)
782	msg += ' does not match distance between specified'
783	msg += ' end points (%.2f m)' %self.length
784	print msg
785
786	self.verbose = verbose
787	self.last_update = 0.0 # For use with update_interval
788	self.last_time = 0.0
789	self.update_interval = update_interval
790
791
792	# Print something
793	if hasattr(self, 'log_filename'):
794	s = 'Culvert Effective Length = %.2f m' %(self.length)
795	log_to_file(self.log_filename, s)
796
797	s = 'Culvert Direction is %s\n' %str(self.vector)
798	log_to_file(self.log_filename, s)
799
800
801
802
803
804	def __call__(self, domain):
805
806	# Time stuff
807	time = domain.get_time()
808
809
810	update = False
811	if self.update_interval is None:
812	update = True
813	delta_t = domain.timestep # Next timestep has been computed in domain.py
814	else:
815	if time - self.last_update > self.update_interval or time == 0.0:
816	update = True
817	delta_t = self.update_interval
818
819	s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
820	if hasattr(self, 'log_filename'):
821	log_to_file(self.log_filename, s)
822
823
824	if update is True:
825	self.last_update = time
826
827	dq = domain.quantities
828
829	# Get average water depths at each opening
830	openings = self.openings # There are two Opening [0] and [1]
831	for i, opening in enumerate(openings):
832
833	# Compute mean values of selected quantitites in the
834	# enquiry area in front of the culvert
835
836	stage = dq['stage'].get_values(location='centroids',
837	indices=[self.enquiry_indices[i]])
838
839	# Store current average stage and depth with each opening object
840	opening.depth = stage - opening.elevation
841	opening.stage = stage
842
843
844
845	################# End of the FOR loop ################################################
846
847	# We now need to deal with each opening individually
848
849	# Determine flow direction based on total energy difference
850
851	delta_w = self.inlet.stage - self.outlet.stage
852
853	if hasattr(self, 'log_filename'):
854	s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time,
855	self.inlet.stage,
856	self.outlet.stage)
857	log_to_file(self.log_filename, s)
858
859
860	if self.inlet.depth <= 0.01:
861	Q = 0.0
862	else:
863	# Calculate discharge for one barrel and set inlet.rate and outlet.rate
864
865	try:
866	Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1])
867	except Below_interval, e:
868	Q = self.rating_curve[0,1]
869	msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
870	msg += str(e)
871	msg += '\n I will use minimum discharge %.2f m^3/s for culvert "%s"'\
872	%(Q, self.label)
873	if hasattr(self, 'log_filename'):
874	log_to_file(self.log_filename, msg)
875	except Above_interval, e:
876	Q = self.rating_curve[-1,1]
877	msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
878	msg += str(e)
879	msg += '\n I will use maximum discharge %.2f m^3/s for culvert "%s"'\
880	%(Q, self.label)
881	if hasattr(self, 'log_filename'):
882	log_to_file(self.log_filename, msg)
883
884
885
886
887	# Adjust discharge for multiple barrels
888	Q *= self.number_of_barrels
889
890
891	# Adjust Q downwards depending on available water at inlet
892	stage = self.inlet.get_quantity_values(quantity_name='stage')
893	elevation = self.inlet.get_quantity_values(quantity_name='elevation')
894	depth = stage-elevation
895
896
897	V = 0
898	for i, d in enumerate(depth):
899	V += d * domain.areas[i]
900
901	#Vsimple = mean(depth)*self.inlet.exchange_area # Current volume in exchange area
902	#print 'Q', Q, 'dt', delta_t, 'Qdt', Qdelta_t, 'V', V, 'Vsimple', Vsimple
903
904	dt = delta_t
905	if Q*dt > V:
906
907	Q_reduced = 0.9*V/dt # Reduce with safety factor
908
909	msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
910	msg += 'greater than current volume (V) at inlet.\n'
911	msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
912
913	#print msg
914
915	if self.verbose is True:
916	print msg
917	if hasattr(self, 'log_filename'):
918	log_to_file(self.log_filename, msg)
919
920	Q = Q_reduced
921
922	self.inlet.rate = -Q
923	self.outlet.rate = Q
924
925	# Log timeseries to file
926	try:
927	fid = open(self.timeseries_filename, 'a')
928	except:
929	pass
930	else:
931	fid.write('%.2f, %.2f\n' %(time, Q))
932	fid.close()
933
934	# Store value of time
935	self.last_time = time
936
937
938
939	# Execute flow term for each opening
940	# This is where Inflow objects are evaluated using the last rate that has been calculated
941	#
942	# This will take place at every internal timestep and update the domain
943	for opening in self.openings:
944	opening(domain)
945
946
947
948
949
950
951	class Culvert_flow_energy:
952	"""Culvert flow - transfer water from one hole to another
953
954	Using Momentum as Calculated by Culvert Flow !!
955	Could be Several Methods Investigated to do This !!!
956
957	2008_May_08
958	To Ole:
959	OK so here we need to get the Polygon Creating code to create
960	polygons for the culvert Based on
961	the two input Points (X0,Y0) and (X1,Y1) - need to be passed
962	to create polygon
963
964	The two centers are now passed on to create_polygon.
965
966
967	Input: Two points, pipe_size (either diameter or width, height),
968	mannings_rougness,
969	inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
970	top-down_blockage_factor and bottom_up_blockage_factor
971
972
973	And the Delta H enquiry should be change from Openings in line 412
974	to the enquiry Polygons infront of the culvert
975	At the moment this script uses only Depth, later we can change it to
976	Energy...
977
978	Once we have Delta H can calculate a Flow Rate and from Flow Rate
979	an Outlet Velocity
980	The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
981
982	Invert levels are optional. If left out they will default to the
983	elevation at the opening.
984
985	"""
986
987	def __init__(self,
988	domain,
989	label=None,
990	description=None,
991	end_point0=None,
992	end_point1=None,
993	width=None,
994	height=None,
995	diameter=None,
996	manning=None, # Mannings Roughness for Culvert
997	invert_level0=None, # Invert level at opening 0
998	invert_level1=None, # Invert level at opening 1
999	loss_exit=None,
1000	loss_entry=None,
1001	loss_bend=None,
1002	loss_special=None,
1003	blockage_topdwn=None,
1004	blockage_bottup=None,
1005	culvert_routine=None,
1006	number_of_barrels=1,
1007	enquiry_point0=None,
1008	enquiry_point1=None,
1009	update_interval=None,
1010	verbose=False):
1011
1012	# Input check
1013	if diameter is not None:
1014	self.culvert_type = 'circle'
1015	self.diameter = diameter
1016	if height is not None or width is not None:
1017	msg = 'Either diameter or width&height must be specified, '
1018	msg += 'but not both.'
1019	raise Exception, msg
1020	else:
1021	if height is not None:
1022	if width is None:
1023	self.culvert_type = 'square'
1024	width = height
1025	else:
1026	self.culvert_type = 'rectangle'
1027	elif width is not None:
1028	if height is None:
1029	self.culvert_type = 'square'
1030	height = width
1031	else:
1032	msg = 'Either diameter or width&height must be specified.'
1033	raise Exception, msg
1034
1035	if height == width:
1036	self.culvert_type = 'square'
1037
1038	self.height = height
1039	self.width = width
1040
1041
1042	assert self.culvert_type in ['circle', 'square', 'rectangle']
1043
1044	assert number_of_barrels >= 1
1045	self.number_of_barrels = number_of_barrels
1046
1047
1048	# Set defaults
1049	if manning is None: manning = 0.012 # Default roughness for pipe
1050	if loss_exit is None: loss_exit = 1.00
1051	if loss_entry is None: loss_entry = 0.50
1052	if loss_bend is None: loss_bend=0.00
1053	if loss_special is None: loss_special=0.00
1054	if blockage_topdwn is None: blockage_topdwn=0.00
1055	if blockage_bottup is None: blockage_bottup=0.00
1056	if culvert_routine is None:
1057	culvert_routine=boyd_generalised_culvert_model
1058
1059	if label is None: label = 'culvert_flow'
1060	label += '_' + str(id(self))
1061	self.label = label
1062
1063	# File for storing culvert quantities
1064	self.timeseries_filename = label + '_timeseries.csv'
1065	fid = open(self.timeseries_filename, 'w')
1066	fid.write('time, E0, E1, Velocity, Discharge\n')
1067	fid.close()
1068
1069	# Log file for storing general textual output
1070	self.log_filename = label + '.log'
1071	log_to_file(self.log_filename, self.label)
1072	log_to_file(self.log_filename, description)
1073	log_to_file(self.log_filename, self.culvert_type)
1074
1075
1076	# Create the fundamental culvert polygons from POLYGON
1077	if self.culvert_type == 'circle':
1078	# Redefine width and height for use with create_culvert_polygons
1079	width = height = diameter
1080
1081	P = create_culvert_polygons(end_point0,
1082	end_point1,
1083	width=width,
1084	height=height,
1085	number_of_barrels=number_of_barrels)
1086
1087	# Select enquiry points
1088	if enquiry_point0 is None:
1089	enquiry_point0 = P['enquiry_point0']
1090
1091	if enquiry_point1 is None:
1092	enquiry_point1 = P['enquiry_point1']
1093
1094	if verbose is True:
1095	pass
1096	#plot_polygons([[end_point0, end_point1],
1097	# P['exchange_polygon0'],
1098	# P['exchange_polygon1'],
1099	# [enquiry_point0, 1.005*enquiry_point0],
1100	# [enquiry_point1, 1.005*enquiry_point1]],
1101	# figname='culvert_polygon_output')
1102
1103
1104	self.enquiry_points = [enquiry_point0, enquiry_point1]
1105
1106
1107	self.enquiry_indices = []
1108	for point in self.enquiry_points:
1109	# Find nearest centroid
1110	N = len(domain)
1111	points = domain.get_centroid_coordinates(absolute=True)
1112
1113	# Calculate indices in exchange area for this forcing term
1114
1115	triangle_id = min_dist = sys.maxint
1116	for k in range(N):
1117	x, y = points[k,:] # Centroid
1118
1119	c = point
1120	distance = (x-c[0])2+(y-c[1])2
1121	if distance < min_dist:
1122	min_dist = distance
1123	triangle_id = k
1124
1125
1126	if triangle_id < sys.maxint:
1127	msg = 'found triangle with centroid (%f, %f)'\
1128	%tuple(points[triangle_id, :])
1129	msg += ' for point (%f, %f)' %tuple(point)
1130
1131	self.enquiry_indices.append(triangle_id)
1132	else:
1133	msg = 'Triangle not found for point (%f, %f)' %point
1134	raise Exception, msg
1135
1136
1137
1138
1139
1140
1141	# Check that all polygons lie within the mesh.
1142	bounding_polygon = domain.get_boundary_polygon()
1143	for key in P.keys():
1144	if key in ['exchange_polygon0',
1145	'exchange_polygon1']:
1146	for point in P[key]:
1147
1148	msg = 'Point %s in polygon %s for culvert %s did not'\
1149	%(str(point), key, self.label)
1150	msg += 'fall within the domain boundary.'
1151	assert is_inside_polygon(point, bounding_polygon), msg
1152
1153
1154	# Create inflow object at each end of the culvert.
1155	self.openings = []
1156	self.openings.append(Inflow(domain,
1157	polygon=P['exchange_polygon0']))
1158
1159	self.openings.append(Inflow(domain,
1160	polygon=P['exchange_polygon1']))
1161
1162
1163	# Assume two openings for now: Referred to as 0 and 1
1164	assert len(self.openings) == 2
1165
1166	# Store basic geometry
1167	self.end_points = [end_point0, end_point1]
1168	self.invert_levels = [invert_level0, invert_level1]
1169	#self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
1170	#self.enquiry_polylines = [P['enquiry_polygon0'][:2],
1171	# P['enquiry_polygon1'][:2]]
1172	self.vector = P['vector']
1173	self.length = P['length']; assert self.length > 0.0
1174	self.verbose = verbose
1175	self.last_time = 0.0
1176	self.last_update = 0.0 # For use with update_interval
1177	self.update_interval = update_interval
1178
1179
1180	# Store hydraulic parameters
1181	self.manning = manning
1182	self.loss_exit = loss_exit
1183	self.loss_entry = loss_entry
1184	self.loss_bend = loss_bend
1185	self.loss_special = loss_special
1186	self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
1187	self.blockage_topdwn = blockage_topdwn
1188	self.blockage_bottup = blockage_bottup
1189
1190	# Store culvert routine
1191	self.culvert_routine = culvert_routine
1192
1193
1194	# Create objects to update momentum (a bit crude at this stage)
1195
1196
1197	xmom0 = General_forcing(domain, 'xmomentum',
1198	polygon=P['exchange_polygon0'])
1199
1200	xmom1 = General_forcing(domain, 'xmomentum',
1201	polygon=P['exchange_polygon1'])
1202
1203	ymom0 = General_forcing(domain, 'ymomentum',
1204	polygon=P['exchange_polygon0'])
1205
1206	ymom1 = General_forcing(domain, 'ymomentum',
1207	polygon=P['exchange_polygon1'])
1208
1209	self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
1210
1211
1212	# Print something
1213	s = 'Culvert Effective Length = %.2f m' %(self.length)
1214	log_to_file(self.log_filename, s)
1215
1216	s = 'Culvert Direction is %s\n' %str(self.vector)
1217	log_to_file(self.log_filename, s)
1218
1219
1220	def __call__(self, domain):
1221
1222	log_filename = self.log_filename
1223
1224	# Time stuff
1225	time = domain.get_time()
1226
1227	# Short hand
1228	dq = domain.quantities
1229
1230
1231	update = False
1232	if self.update_interval is None:
1233	update = True
1234	delta_t = domain.timestep # Next timestep has been computed in domain.py
1235	else:
1236	if time - self.last_update > self.update_interval or time == 0.0:
1237	update = True
1238	delta_t = self.update_interval
1239
1240	s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
1241	if hasattr(self, 'log_filename'):
1242	log_to_file(log_filename, s)
1243
1244
1245	if update is True:
1246	self.last_update = time
1247
1248	msg = 'Time did not advance'
1249	if time > 0.0: assert delta_t > 0.0, msg
1250
1251
1252	# Get average water depths at each opening
1253	openings = self.openings # There are two Opening [0] and [1]
1254	for i, opening in enumerate(openings):
1255
1256	# Compute mean values of selected quantitites in the
1257	# exchange area in front of the culvert
1258
1259	stage = opening.get_quantity_values(quantity_name='stage')
1260	w = mean(stage) # Average stage
1261
1262	# Use invert level instead of elevation if specified
1263	invert_level = self.invert_levels[i]
1264	if invert_level is not None:
1265	z = invert_level
1266	else:
1267	elevation = opening.get_quantity_values(quantity_name='elevation')
1268	z = mean(elevation) # Average elevation
1269
1270	# Estimated depth above the culvert inlet
1271	d = w - z # Used for calculations involving depth
1272	if d < 0.0:
1273	# This is possible since w and z are taken at different locations
1274	#msg = 'D < 0.0: %f' %d
1275	#raise msg
1276	d = 0.0
1277
1278
1279	# Ratio of depth to culvert height.
1280	# If ratio > 1 then culvert is running full
1281	if self.culvert_type == 'circle':
1282	ratio = d/self.diameter
1283	else:
1284	ratio = d/self.height
1285	opening.ratio = ratio
1286
1287
1288	# Average measures of energy in front of this opening
1289
1290	id = [self.enquiry_indices[i]]
1291	stage = dq['stage'].get_values(location='centroids',
1292	indices=id)
1293	elevation = dq['elevation'].get_values(location='centroids',
1294	indices=id)
1295	xmomentum = dq['xmomentum'].get_values(location='centroids',
1296	indices=id)
1297	ymomentum = dq['xmomentum'].get_values(location='centroids',
1298	indices=id)
1299	depth = stage-elevation
1300	if depth > 0.0:
1301	u = xmomentum/(depth + velocity_protection/depth)
1302	v = ymomentum/(depth + velocity_protection/depth)
1303	else:
1304	u = v = 0.0
1305
1306
1307	opening.total_energy = 0.5(uu + v*v)/g + stage
1308	#print 'Et = %.3f m' %opening.total_energy
1309
1310	# Store current average stage and depth with each opening object
1311	opening.depth = d
1312	opening.depth_trigger = d # Use this for now
1313	opening.stage = w
1314	opening.elevation = z
1315
1316
1317	################# End of the FOR loop ################################################
1318
1319	# We now need to deal with each opening individually
1320
1321	# Determine flow direction based on total energy difference
1322	delta_Et = openings[0].total_energy - openings[1].total_energy
1323
1324	if delta_Et > 0:
1325	#print 'Flow U/S ---> D/S'
1326	inlet = openings[0]
1327	outlet = openings[1]
1328
1329	inlet.momentum = self.opening_momentum[0]
1330	outlet.momentum = self.opening_momentum[1]
1331
1332	else:
1333	#print 'Flow D/S ---> U/S'
1334	inlet = openings[1]
1335	outlet = openings[0]
1336
1337	inlet.momentum = self.opening_momentum[1]
1338	outlet.momentum = self.opening_momentum[0]
1339
1340	delta_Et = -delta_Et
1341
1342	self.inlet = inlet
1343	self.outlet = outlet
1344
1345	msg = 'Total energy difference is negative'
1346	assert delta_Et > 0.0, msg
1347
1348	delta_h = inlet.stage - outlet.stage
1349	delta_z = inlet.elevation - outlet.elevation
1350	culvert_slope = (delta_z/self.length)
1351
1352	if culvert_slope < 0.0:
1353	# Adverse gradient - flow is running uphill
1354	# Flow will be purely controlled by uphill outlet face
1355	if self.verbose is True:
1356	print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
1357
1358
1359	s = 'Delta total energy = %.3f' %(delta_Et)
1360	log_to_file(log_filename, s)
1361
1362
1363	# Calculate discharge for one barrel and set inlet.rate and outlet.rate
1364	Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
1365
1366	# Adjust discharge for multiple barrels
1367	Q *= self.number_of_barrels
1368
1369	# Compute barrel momentum
1370	barrel_momentum = barrel_velocity*culvert_outlet_depth
1371
1372	s = 'Barrel velocity = %f' %barrel_velocity
1373	log_to_file(log_filename, s)
1374
1375	# Compute momentum vector at outlet
1376	outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
1377
1378	s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
1379	log_to_file(log_filename, s)
1380
1381	# Log timeseries to file
1382	fid = open(self.timeseries_filename, 'a')
1383	fid.write('%f, %f, %f, %f, %f\n'\
1384	%(time,
1385	openings[0].total_energy,
1386	openings[1].total_energy,
1387	barrel_velocity,
1388	Q))
1389	fid.close()
1390
1391	# Update momentum
1392
1393	if delta_t > 0.0:
1394	xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
1395	xmomentum_rate /= delta_t
1396
1397	ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
1398	ymomentum_rate /= delta_t
1399
1400	s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
1401	log_to_file(log_filename, s)
1402	else:
1403	xmomentum_rate = ymomentum_rate = 0.0
1404
1405
1406	# Set momentum rates for outlet jet
1407	outlet.momentum[0].rate = xmomentum_rate
1408	outlet.momentum[1].rate = ymomentum_rate
1409
1410	# Remember this value for next step (IMPORTANT)
1411	outlet.momentum[0].value = outlet_mom_x
1412	outlet.momentum[1].value = outlet_mom_y
1413
1414	if int(domain.time*100) % 100 == 0:
1415	s = 'T=%.5f, Culvert Discharge = %.3f f'\
1416	%(time, Q)
1417	s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\
1418	%(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
1419	s += ' Momentum rate: (%.4f, %.4f)'\
1420	%(xmomentum_rate, ymomentum_rate)
1421	s+='Outlet Vel= %.3f'\
1422	%(barrel_velocity)
1423	log_to_file(log_filename, s)
1424
1425	# Store value of time
1426	self.last_time = time
1427
1428
1429
1430	# Execute flow term for each opening
1431	# This is where Inflow objects are evaluated and update the domain
1432	for opening in self.openings:
1433	opening(domain)
1434
1435	# Execute momentum terms
1436	# This is where Inflow objects are evaluated and update the domain
1437	self.outlet.momentum[0](domain)
1438	self.outlet.momentum[1](domain)
1439
1440
1441
1442	Culvert_flow = Culvert_flow_general

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

Download in other formats: