Context Navigation

source: trunk/anuga_core/source/anuga/shallow_water/forcing.py @ 7967

Last change on this file since 7967 was 7967, checked in by steve, 13 years ago
Commits from session with Nariman
File size: 28.7 KB

Line
1	"""
2	Environmental forcing functions, such as wind and rainfall.
3
4	Constraints: See GPL license in the user guide
5	Version: 1.0 ($Revision: 7731 $)
6	ModifiedBy:
7	$Author: hudson $
8	$Date: 2010-05-18 14:54:05 +1000 (Tue, 18 May 2010) $
9	"""
10
11
12	from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
13	from anuga.utilities.numerical_tools import ensure_numeric
14	from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
15	Modeltime_too_late
16	from anuga.geometry.polygon import is_inside_polygon, inside_polygon, \
17	polygon_area
18	from types import IntType, FloatType
19	from anuga.geospatial_data.geospatial_data import ensure_geospatial
20
21	from warnings import warn
22	import numpy as num
23
24
25
26	def check_forcefield(f):
27	"""Check that force object is as expected.
28
29	Check that f is either:
30	1: a callable object f(t,x,y), where x and y are vectors
31	and that it returns an array or a list of same length
32	as x and y
33	2: a scalar
34	"""
35
36	if callable(f):
37	N = 3
38	x = num.ones(3, num.float)
39	y = num.ones(3, num.float)
40	try:
41	q = f(1.0, x=x, y=y)
42	except Exception, e:
43	msg = 'Function %s could not be executed:\n%s' %(f, e)
44	# FIXME: Reconsider this semantics
45	raise Exception, msg
46
47	try:
48	q = num.array(q, num.float)
49	except:
50	msg = ('Return value from vector function %s could not '
51	'be converted into a numeric array of floats.\nSpecified '
52	'function should return either list or array.' % f)
53	raise Exception, msg
54
55	# Is this really what we want?
56	# info is "(func name, filename, defining line)"
57	func_info = (f.func_name, f.func_code.co_filename,
58	f.func_code.co_firstlineno)
59	func_msg = 'Function %s (defined in %s, line %d)' % func_info
60	try:
61	result_len = len(q)
62	except:
63	msg = '%s must return vector' % func_msg
64	self.fail(msg)
65	msg = '%s must return vector of length %d' % (func_msg, N)
66	assert result_len == N, msg
67	else:
68	try:
69	f = float(f)
70	except:
71	msg = ('Force field %s must be a scalar value coercible to float.'
72	% str(f))
73	raise Exception, msg
74
75	return f
76
77
78
79	class Wind_stress:
80	"""Apply wind stress to water momentum in terms of
81	wind speed [m/s] and wind direction [degrees]
82	"""
83	def __init__(self, args, *kwargs):
84	"""Initialise windfield from wind speed s [m/s]
85	and wind direction phi [degrees]
86
87	Inputs v and phi can be either scalars or Python functions, e.g.
88
89	W = Wind_stress(10, 178)
90
91	#FIXME - 'normal' degrees are assumed for now, i.e. the
92	vector (1,0) has zero degrees.
93	We may need to convert from 'compass' degrees later on and also
94	map from True north to grid north.
95
96	Arguments can also be Python functions of t,x,y as in
97
98	def speed(t,x,y):
99	...
100	return s
101
102	def angle(t,x,y):
103	...
104	return phi
105
106	where x and y are vectors.
107
108	and then pass the functions in
109
110	W = Wind_stress(speed, angle)
111
112	The instantiated object W can be appended to the list of
113	forcing_terms as in
114
115	Alternatively, one vector valued function for (speed, angle)
116	can be applied, providing both quantities simultaneously.
117	As in
118	W = Wind_stress(F), where returns (speed, angle) for each t.
119
120	domain.forcing_terms.append(W)
121	"""
122
123	from anuga.config import rho_a, rho_w, eta_w
124
125	if len(args) == 2:
126	s = args[0]
127	phi = args[1]
128	elif len(args) == 1:
129	# Assume vector function returning (s, phi)(t,x,y)
130	vector_function = args[0]
131	s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
132	phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
133	else:
134	# Assume info is in 2 keyword arguments
135	if len(kwargs) == 2:
136	s = kwargs['s']
137	phi = kwargs['phi']
138	else:
139	raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
140
141	self.speed = check_forcefield(s)
142	self.phi = check_forcefield(phi)
143
144	self.const = eta_w*rho_a/rho_w
145
146	##
147	# @brief 'execute' this class instance.
148	# @param domain
149	def __call__(self, domain):
150	"""Evaluate windfield based on values found in domain"""
151
152	from math import pi, cos, sin, sqrt
153
154	xmom_update = domain.quantities['xmomentum'].explicit_update
155	ymom_update = domain.quantities['ymomentum'].explicit_update
156
157	N = len(domain) # number_of_triangles
158	t = domain.time
159
160	if callable(self.speed):
161	xc = domain.get_centroid_coordinates()
162	s_vec = self.speed(t, xc[:,0], xc[:,1])
163	else:
164	# Assume s is a scalar
165	try:
166	s_vec = self.speed * num.ones(N, num.float)
167	except:
168	msg = 'Speed must be either callable or a scalar: %s' %self.s
169	raise msg
170
171	if callable(self.phi):
172	xc = domain.get_centroid_coordinates()
173	phi_vec = self.phi(t, xc[:,0], xc[:,1])
174	else:
175	# Assume phi is a scalar
176
177	try:
178	phi_vec = self.phi * num.ones(N, num.float)
179	except:
180	msg = 'Angle must be either callable or a scalar: %s' %self.phi
181	raise msg
182
183	assign_windfield_values(xmom_update, ymom_update,
184	s_vec, phi_vec, self.const)
185
186
187	##
188	# @brief Assign wind field values
189	# @param xmom_update
190	# @param ymom_update
191	# @param s_vec
192	# @param phi_vec
193	# @param const
194	def assign_windfield_values(xmom_update, ymom_update,
195	s_vec, phi_vec, const):
196	"""Python version of assigning wind field to update vectors.
197	A C version also exists (for speed)
198	"""
199
200	from math import pi, cos, sin, sqrt
201
202	N = len(s_vec)
203	for k in range(N):
204	s = s_vec[k]
205	phi = phi_vec[k]
206
207	# Convert to radians
208	phi = phi*pi/180
209
210	# Compute velocity vector (u, v)
211	u = s*cos(phi)
212	v = s*sin(phi)
213
214	# Compute wind stress
215	S = const * sqrt(u2 + v2)
216	xmom_update[k] += S*u
217	ymom_update[k] += S*v
218
219
220	##
221	# @brief A class for a general explicit forcing term.
222	class General_forcing:
223	"""General explicit forcing term for update of quantity
224
225	This is used by Inflow and Rainfall for instance
226
227
228	General_forcing(quantity_name, rate, center, radius, polygon)
229
230	domain: ANUGA computational domain
231	quantity_name: Name of quantity to update.
232	It must be a known conserved quantity.
233
234	rate [?/s]: Total rate of change over the specified area.
235	This parameter can be either a constant or a
236	function of time. Positive values indicate increases,
237	negative values indicate decreases.
238	Rate can be None at initialisation but must be specified
239	before forcing term is applied (i.e. simulation has started).
240
241	center [m]: Coordinates at center of flow point
242	radius [m]: Size of circular area
243	polygon: Arbitrary polygon
244	default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
245	Admissible types: None, constant number or function of t
246
247
248	Either center, radius or polygon can be specified but not both.
249	If neither are specified the entire domain gets updated.
250	All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.
251
252	Inflow or Rainfall for examples of use
253	"""
254
255
256	# FIXME (AnyOne) : Add various methods to allow spatial variations
257
258	##
259	# @brief Create an instance of this forcing term.
260	# @param domain
261	# @param quantity_name
262	# @param rate
263	# @param center
264	# @param radius
265	# @param polygon
266	# @param default_rate
267	# @param verbose
268	def __init__(self,
269	domain,
270	quantity_name,
271	rate=0.0,
272	center=None,
273	radius=None,
274	polygon=None,
275	default_rate=None,
276	verbose=False):
277
278	from math import pi, cos, sin
279
280	if center is None:
281	msg = 'I got radius but no center.'
282	assert radius is None, msg
283
284	if radius is None:
285	msg += 'I got center but no radius.'
286	assert center is None, msg
287
288	self.domain = domain
289	self.quantity_name = quantity_name
290	self.rate = rate
291	self.center = ensure_numeric(center)
292	self.radius = radius
293	self.polygon = polygon
294	self.verbose = verbose
295	self.value = 0.0 # Can be used to remember value at
296	# previous timestep in order to obtain rate
297
298	# Get boundary (in absolute coordinates)
299	bounding_polygon = domain.get_boundary_polygon()
300
301	# Update area if applicable
302	if center is not None and radius is not None:
303	assert len(center) == 2
304	msg = 'Polygon cannot be specified when center and radius are'
305	assert polygon is None, msg
306
307	# Check that circle center lies within the mesh.
308	msg = 'Center %s specified for forcing term did not' % str(center)
309	msg += 'fall within the domain boundary.'
310	assert is_inside_polygon(center, bounding_polygon), msg
311
312	# Check that circle periphery lies within the mesh.
313	N = 100
314	periphery_points = []
315	for i in range(N):
316	theta = 2pii/100
317
318	x = center[0] + radius*cos(theta)
319	y = center[1] + radius*sin(theta)
320
321	periphery_points.append([x,y])
322
323	for point in periphery_points:
324	msg = 'Point %s on periphery for forcing term' % str(point)
325	msg += ' did not fall within the domain boundary.'
326	assert is_inside_polygon(point, bounding_polygon), msg
327
328	if polygon is not None:
329	# Check that polygon lies within the mesh.
330	for point in self.polygon:
331	msg = 'Point %s in polygon for forcing term' % str(point)
332	msg += ' did not fall within the domain boundary.'
333	assert is_inside_polygon(point, bounding_polygon), msg
334
335	# Pointer to update vector
336	self.update = domain.quantities[self.quantity_name].explicit_update
337
338	# Determine indices in flow area
339	N = len(domain)
340	points = domain.get_centroid_coordinates(absolute=True)
341
342	# Calculate indices in exchange area for this forcing term
343	self.exchange_indices = None
344	if self.center is not None and self.radius is not None:
345	# Inlet is circular
346	inlet_region = 'center=%s, radius=%s' % (self.center, self.radius)
347
348	self.exchange_indices = []
349	for k in range(N):
350	x, y = points[k,:] # Centroid
351
352	c = self.center
353	if ((x-c[0])2+(y-c[1])2) < self.radius**2:
354	self.exchange_indices.append(k)
355
356	if self.polygon is not None:
357	# Inlet is polygon
358	self.exchange_indices = inside_polygon(points, self.polygon)
359
360	if self.exchange_indices is None:
361	self.exchange_area = polygon_area(bounding_polygon)
362	else:
363	if len(self.exchange_indices) == 0:
364	msg = 'No triangles have been identified in '
365	msg += 'specified region: %s' % inlet_region
366	raise Exception, msg
367
368	# Compute exchange area as the sum of areas of triangles identified
369	# by circle or polygon
370	self.exchange_area = 0.0
371	for i in self.exchange_indices:
372	self.exchange_area += domain.areas[i]
373
374
375	msg = 'Exchange area in forcing term'
376	msg += ' has area = %f' %self.exchange_area
377	assert self.exchange_area > 0.0
378
379
380
381
382	# Check and store default_rate
383	msg = ('Keyword argument default_rate must be either None '
384	'or a function of time.\nI got %s.' % str(default_rate))
385	assert (default_rate is None or
386	type(default_rate) in [IntType, FloatType] or
387	callable(default_rate)), msg
388
389	if default_rate is not None:
390	# If it is a constant, make it a function
391	if not callable(default_rate):
392	tmp = default_rate
393	default_rate = lambda t: tmp
394
395	# Check that default_rate is a function of one argument
396	try:
397	default_rate(0.0)
398	except:
399	raise Exception, msg
400
401	self.default_rate = default_rate
402	self.default_rate_invoked = False # Flag
403
404	##
405	# @brief Execute this instance.
406	# @param domain
407	def __call__(self, domain):
408	"""Apply inflow function at time specified in domain, update stage"""
409
410	# Call virtual method allowing local modifications
411	t = domain.get_time()
412	try:
413	rate = self.update_rate(t)
414	except Modeltime_too_early, e:
415	raise Modeltime_too_early, e
416	except Modeltime_too_late, e:
417	if self.default_rate is None:
418	msg = '%s: ANUGA is trying to run longer than specified data.\n' %str(e)
419	msg += 'You can specify keyword argument default_rate in the '
420	msg += 'forcing function to tell it what to do in the absence of time data.'
421	raise Modeltime_too_late, msg
422	else:
423	# Pass control to default rate function
424	rate = self.default_rate(t)
425
426	if self.default_rate_invoked is False:
427	# Issue warning the first time
428	msg = ('%s\n'
429	'Instead I will use the default rate: %s\n'
430	'Note: Further warnings will be supressed'
431	% (str(e), str(self.default_rate)))
432	warn(msg)
433
434	# FIXME (Ole): Replace this crude flag with
435	# Python's ability to print warnings only once.
436	# See http://docs.python.org/lib/warning-filter.html
437	self.default_rate_invoked = True
438
439	if rate is None:
440	msg = ('Attribute rate must be specified in General_forcing '
441	'or its descendants before attempting to call it')
442	raise Exception, msg
443
444	# Now rate is a number
445	if self.verbose is True:
446	log.critical('Rate of %s at time = %.2f = %f'
447	% (self.quantity_name, domain.get_time(), rate))
448
449	if self.exchange_indices is None:
450	self.update[:] += rate
451	else:
452	# Brute force assignment of restricted rate
453	for k in self.exchange_indices:
454	self.update[k] += rate
455
456	##
457	# @brief Update the internal rate.
458	# @param t A callable or scalar used to set the rate.
459	# @return The new rate.
460	def update_rate(self, t):
461	"""Virtual method allowing local modifications by writing an
462	overriding version in descendant
463	"""
464
465	if callable(self.rate):
466	rate = self.rate(t)
467	else:
468	rate = self.rate
469
470	return rate
471
472	##
473	# @brief Get values for the specified quantity.
474	# @param quantity_name Name of the quantity of interest.
475	# @return The value(s) of the quantity.
476	# @note If 'quantity_name' is None, use self.quantity_name.
477	def get_quantity_values(self, quantity_name=None):
478	"""Return values for specified quantity restricted to opening
479
480	Optionally a quantity name can be specified if values from another
481	quantity is sought
482	"""
483
484	if quantity_name is None:
485	quantity_name = self.quantity_name
486
487	q = self.domain.quantities[quantity_name]
488	return q.get_values(location='centroids',
489	indices=self.exchange_indices)
490
491	##
492	# @brief Set value for the specified quantity.
493	# @param val The value object used to set value.
494	# @param quantity_name Name of the quantity of interest.
495	# @note If 'quantity_name' is None, use self.quantity_name.
496	def set_quantity_values(self, val, quantity_name=None):
497	"""Set values for specified quantity restricted to opening
498
499	Optionally a quantity name can be specified if values from another
500	quantity is sought
501	"""
502
503	if quantity_name is None:
504	quantity_name = self.quantity_name
505
506	q = self.domain.quantities[self.quantity_name]
507	q.set_values(val,
508	location='centroids',
509	indices=self.exchange_indices)
510
511
512	##
513	# @brief A class for rainfall forcing function.
514	# @note Inherits from General_forcing.
515	class Rainfall(General_forcing):
516	"""Class Rainfall - general 'rain over entire domain' forcing term.
517
518	Used for implementing Rainfall over the entire domain.
519
520	Current Limited to only One Gauge..
521
522	Need to add Spatial Varying Capability
523	(This module came from copying and amending the Inflow Code)
524
525	Rainfall(rain)
526
527	domain
528	rain [mm/s]: Total rain rate over the specified domain.
529	NOTE: Raingauge Data needs to reflect the time step.
530	IE: if Gauge is mm read at a time step, then the input
531	here is as mm/(timeStep) so 10mm in 5minutes becomes
532	10/(5x60) = 0.0333mm/s.
533
534	This parameter can be either a constant or a
535	function of time. Positive values indicate inflow,
536	negative values indicate outflow.
537	(and be used for Infiltration - Write Seperate Module)
538	The specified flow will be divided by the area of
539	the inflow region and then applied to update the
540	stage quantity.
541
542	polygon: Specifies a polygon to restrict the rainfall.
543
544	Examples
545	How to put them in a run File...
546
547	#------------------------------------------------------------------------
548	# Setup specialised forcing terms
549	#------------------------------------------------------------------------
550	# This is the new element implemented by Ole and Rudy to allow direct
551	# input of Rainfall in mm/s
552
553	catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms'))
554	# Note need path to File in String.
555	# Else assumed in same directory
556
557	domain.forcing_terms.append(catchmentrainfall)
558	"""
559
560	##
561	# @brief Create an instance of the class.
562	# @param domain Domain of interest.
563	# @param rate Total rain rate over the specified domain (mm/s).
564	# @param center
565	# @param radius
566	# @param polygon Polygon to restrict rainfall.
567	# @param default_rate
568	# @param verbose True if this instance is to be verbose.
569	def __init__(self,
570	domain,
571	rate=0.0,
572	center=None,
573	radius=None,
574	polygon=None,
575	default_rate=None,
576	verbose=False):
577
578	# Converting mm/s to m/s to apply in ANUGA)
579	if callable(rate):
580	rain = lambda t: rate(t)/1000.0
581	else:
582	rain = rate/1000.0
583
584	if default_rate is not None:
585	if callable(default_rate):
586	default_rain = lambda t: default_rate(t)/1000.0
587	else:
588	default_rain = default_rate/1000.0
589	else:
590	default_rain = None
591
592
593
594	General_forcing.__init__(self,
595	domain,
596	'stage',
597	rate=rain,
598	center=center,
599	radius=radius,
600	polygon=polygon,
601	default_rate=default_rain,
602	verbose=verbose)
603
604
605	##
606	# @brief A class for inflow (rain and drain) forcing function.
607	# @note Inherits from General_forcing.
608	class Inflow(General_forcing):
609	"""Class Inflow - general 'rain and drain' forcing term.
610
611	Useful for implementing flows in and out of the domain.
612
613	Inflow(flow, center, radius, polygon)
614
615	domain
616	rate [m^3/s]: Total flow rate over the specified area.
617	This parameter can be either a constant or a
618	function of time. Positive values indicate inflow,
619	negative values indicate outflow.
620	The specified flow will be divided by the area of
621	the inflow region and then applied to update stage.
622	center [m]: Coordinates at center of flow point
623	radius [m]: Size of circular area
624	polygon: Arbitrary polygon.
625
626	Either center, radius or polygon must be specified
627
628	Examples
629
630	# Constant drain at 0.003 m^3/s.
631	# The outflow area is 0.07*2pi=0.0154 m^2
632	# This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
633	#
634	Inflow((0.7, 0.4), 0.07, -0.003)
635
636
637	# Tap turning up to a maximum inflow of 0.0142 m^3/s.
638	# The inflow area is 0.03*2pi = 0.00283 m^2
639	# This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s
640	# over the specified area
641	Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
642
643
644	#------------------------------------------------------------------------
645	# Setup specialised forcing terms
646	#------------------------------------------------------------------------
647	# This is the new element implemented by Ole to allow direct input
648	# of Inflow in m^3/s
649
650	hydrograph = Inflow(center=(320, 300), radius=10,
651	rate=file_function('Q/QPMF_Rot_Sub13.tms'))
652
653	domain.forcing_terms.append(hydrograph)
654	"""
655
656	##
657	# @brief Create an instance of the class.
658	# @param domain Domain of interest.
659	# @param rate Total rain rate over the specified domain (mm/s).
660	# @param center
661	# @param radius
662	# @param polygon Polygon to restrict rainfall.
663	# @param default_rate
664	# @param verbose True if this instance is to be verbose.
665	def __init__(self,
666	domain,
667	rate=0.0,
668	center=None,
669	radius=None,
670	polygon=None,
671	default_rate=None,
672	verbose=False):
673	# Create object first to make area is available
674	General_forcing.__init__(self,
675	domain,
676	'stage',
677	rate=rate,
678	center=center,
679	radius=radius,
680	polygon=polygon,
681	default_rate=default_rate,
682	verbose=verbose)
683
684	##
685	# @brief Update the instance rate.
686	# @param t New rate object.
687	def update_rate(self, t):
688	"""Virtual method allowing local modifications by writing an
689	overriding version in descendant
690
691	This one converts m^3/s to m/s which can be added directly
692	to 'stage' in ANUGA
693	"""
694
695	if callable(self.rate):
696	_rate = self.rate(t)/self.exchange_area
697	else:
698	_rate = self.rate/self.exchange_area
699
700	return _rate
701
702
703	##
704	# @brief A class for creating cross sections.
705	# @note Inherits from General_forcing.
706	class Cross_section:
707	"""Class Cross_section - a class to setup a cross section from
708	which you can then calculate flow and energy through cross section
709
710
711	Cross_section(domain, polyline)
712
713	domain:
714	polyline: Representation of desired cross section - it may contain
715	multiple sections allowing for complex shapes. Assume
716	absolute UTM coordinates.
717	Format [[x0, y0], [x1, y1], ...]
718	verbose:
719	"""
720
721	##
722	# @brief Create an instance of the class.
723	# @param domain Domain of interest.
724	# @param polyline Polyline defining cross section
725	# @param verbose True if this instance is to be verbose.
726	def __init__(self,
727	domain,
728	polyline=None,
729	verbose=False):
730
731	self.domain = domain
732	self.polyline = polyline
733	self.verbose = verbose
734
735	# Find all intersections and associated triangles.
736	self.segments = self.domain.get_intersecting_segments(self.polyline,
737	use_cache=True,
738	verbose=self.verbose)
739
740	# Get midpoints
741	self.midpoints = segment_midpoints(self.segments)
742
743	# Make midpoints Geospatial instances
744	self.midpoints = ensure_geospatial(self.midpoints, self.domain.geo_reference)
745
746	##
747	# @brief set verbose mode
748	def set_verbose(self,verbose=True):
749	"""Set verbose mode true or flase
750	"""
751
752	self.verbose=verbose
753
754	##
755	# @brief calculate current flow through cross section
756	def get_flow_through_cross_section(self):
757	""" Output: Total flow [m^3/s] across cross section.
758	"""
759
760	# Get interpolated values
761	xmomentum = self.domain.get_quantity('xmomentum')
762	ymomentum = self.domain.get_quantity('ymomentum')
763
764	uh = xmomentum.get_values(interpolation_points=self.midpoints,
765	use_cache=True)
766	vh = ymomentum.get_values(interpolation_points=self.midpoints,
767	use_cache=True)
768
769	# Compute and sum flows across each segment
770	total_flow = 0
771	for i in range(len(uh)):
772	# Inner product of momentum vector with segment normal [m^2/s]
773	normal = self.segments[i].normal
774	normal_momentum = uh[i]normal[0] + vh[i]normal[1]
775
776	# Flow across this segment [m^3/s]
777	segment_flow = normal_momentum*self.segments[i].length
778
779	# Accumulate
780	total_flow += segment_flow
781
782	return total_flow
783
784
785	##
786	# @brief calculate current energy flow through cross section
787	def get_energy_through_cross_section(self, kind='total'):
788	"""Obtain average energy head [m] across specified cross section.
789
790	Output:
791	E: Average energy [m] across given segments for all stored times.
792
793	The average velocity is computed for each triangle intersected by
794	the polyline and averaged weighted by segment lengths.
795
796	The typical usage of this function would be to get average energy of
797	flow in a channel, and the polyline would then be a cross section
798	perpendicular to the flow.
799
800	#FIXME (Ole) - need name for this energy reflecting that its dimension
801	is [m].
802	"""
803
804	from anuga.config import g, epsilon, velocity_protection as h0
805
806	# Get interpolated values
807	stage = self.domain.get_quantity('stage')
808	elevation = self.domain.get_quantity('elevation')
809	xmomentum = self.domain.get_quantity('xmomentum')
810	ymomentum = self.domain.get_quantity('ymomentum')
811
812	w = stage.get_values(interpolation_points=self.midpoints, use_cache=True)
813	z = elevation.get_values(interpolation_points=self.midpoints, use_cache=True)
814	uh = xmomentum.get_values(interpolation_points=self.midpoints,
815	use_cache=True)
816	vh = ymomentum.get_values(interpolation_points=self.midpoints,
817	use_cache=True)
818	h = w-z # Depth
819
820	# Compute total length of polyline for use with weighted averages
821	total_line_length = 0.0
822	for segment in self.segments:
823	total_line_length += segment.length
824
825	# Compute and sum flows across each segment
826	average_energy = 0.0
827	for i in range(len(w)):
828	# Average velocity across this segment
829	if h[i] > epsilon:
830	# Use protection against degenerate velocities
831	u = uh[i]/(h[i] + h0/h[i])
832	v = vh[i]/(h[i] + h0/h[i])
833	else:
834	u = v = 0.0
835
836	speed_squared = uu + vv
837	kinetic_energy = 0.5*speed_squared/g
838
839	if kind == 'specific':
840	segment_energy = h[i] + kinetic_energy
841	elif kind == 'total':
842	segment_energy = w[i] + kinetic_energy
843	else:
844	msg = 'Energy kind must be either "specific" or "total".'
845	msg += ' I got %s' %kind
846
847	# Add to weighted average
848	weigth = self.segments[i].length/total_line_length
849	average_energy += segment_energy*weigth
850
851	return average_energy
852

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

Download in other formats: