Context Navigation

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

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

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

Download in other formats: