Context Navigation

culvert_class.py @ 7184

Last change on this file since 7184 was 7035, checked in by rwilson, 16 years ago
Back-merge from Numeric trunk.
File size: 53.9 KB

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

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

Context Navigation

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

Download in other formats: