Context Navigation

source: trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 7939

Last change on this file since 7939 was 7939, checked in by steve, 13 years ago
Adding in culvert routines into structures directory
Property svn:keywords set to `LastChangedDate LastChangedRevision LastChangedBy HeadURL Id`
File size: 47.7 KB

Line
1	"""
2	Finite-volume computations of the shallow water wave equation.
3
4	Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
5	computations of the shallow water wave equation.
6
7
8	Author: Ole Nielsen, Ole.Nielsen@ga.gov.au
9	Stephen Roberts, Stephen.Roberts@anu.edu.au
10	Duncan Gray, Duncan.Gray@ga.gov.au
11
12	CreationDate: 2004
13
14	Description:
15	This module contains a specialisation of class Domain from
16	module domain.py consisting of methods specific to the
17	Shallow Water Wave Equation
18
19	U_t + E_x + G_y = S
20
21	where
22
23	U = [w, uh, vh]
24	E = [uh, u^2h + gh^2/2, uvh]
25	G = [vh, uvh, v^2h + gh^2/2]
26	S represents source terms forcing the system
27	(e.g. gravity, friction, wind stress, ...)
28
29	and _t, _x, _y denote the derivative with respect to t, x and y
30	respectively.
31
32
33	The quantities are
34
35	symbol variable name explanation
36	x x horizontal distance from origin [m]
37	y y vertical distance from origin [m]
38	z elevation elevation of bed on which flow is modelled [m]
39	h height water height above z [m]
40	w stage absolute water level, w = z+h [m]
41	u speed in the x direction [m/s]
42	v speed in the y direction [m/s]
43	uh xmomentum momentum in the x direction [m^2/s]
44	vh ymomentum momentum in the y direction [m^2/s]
45
46	eta mannings friction coefficient [to appear]
47	nu wind stress coefficient [to appear]
48
49	The conserved quantities are w, uh, vh
50
51	Reference:
52	Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
53	Christopher Zoppou and Stephen Roberts,
54	Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
55
56	Hydrodynamic modelling of coastal inundation.
57	Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
58	In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
59	Modelling and Simulation. Modelling and Simulation Society of Australia and
60	New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.
61	http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
62
63
64	SeeAlso:
65	TRAC administration of ANUGA (User Manuals etc) at
66	https://datamining.anu.edu.au/anuga and Subversion repository at
67	$HeadURL: https://datamining.anu.edu.au/svn/anuga/trunk/anuga_core/source/
68	anuga/shallow_water/shallow_water_domain.py $
69
70	Constraints: See GPL license in the user guide
71	Version: 1.0 ($Revision: 7939 $)
72	ModifiedBy:
73	$Author: steve $
74	$Date: 2010-08-13 05:35:39 +0000 (Fri, 13 Aug 2010) $
75	"""
76
77
78	import numpy as num
79
80	from anuga.abstract_2d_finite_volumes.generic_domain \
81	import Generic_Domain
82
83	from anuga.shallow_water.forcing import Cross_section
84	from anuga.utilities.numerical_tools import mean
85	from anuga.file.sww import SWW_file
86
87	import anuga.utilities.log as log
88
89	import types
90
91	class Domain(Generic_Domain):
92	""" Class for a shallow water domain."""
93	def __init__(self,
94	coordinates=None,
95	vertices=None,
96	boundary=None,
97	tagged_elements=None,
98	geo_reference=None,
99	use_inscribed_circle=False,
100	mesh_filename=None,
101	use_cache=False,
102	verbose=False,
103	conserved_quantities = None,
104	evolved_quantities = None,
105	other_quantities = None,
106	full_send_dict=None,
107	ghost_recv_dict=None,
108	processor=0,
109	numproc=1,
110	number_of_full_nodes=None,
111	number_of_full_triangles=None):
112	"""
113	Instantiate a shallow water domain.
114	coordinates - vertex locations for the mesh
115	vertices - vertex indices for the mesh
116	boundary - boundaries of the mesh
117	# @param tagged_elements
118	# @param geo_reference
119	# @param use_inscribed_circle
120	# @param mesh_filename
121	# @param use_cache
122	# @param verbose
123	# @param evolved_quantities
124	# @param full_send_dict
125	# @param ghost_recv_dict
126	# @param processor
127	# @param numproc
128	# @param number_of_full_nodes
129	# @param number_of_full_triangles
130	"""
131
132	# Define quantities for the shallow_water domain
133	if conserved_quantities == None:
134	conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
135
136	if evolved_quantities == None:
137	evolved_quantities = ['stage', 'xmomentum', 'ymomentum']
138
139	if other_quantities == None:
140	other_quantities = ['elevation', 'friction']
141
142	Generic_Domain.__init__(self,
143	coordinates,
144	vertices,
145	boundary,
146	conserved_quantities,
147	evolved_quantities,
148	other_quantities,
149	tagged_elements,
150	geo_reference,
151	use_inscribed_circle,
152	mesh_filename,
153	use_cache,
154	verbose,
155	full_send_dict,
156	ghost_recv_dict,
157	processor,
158	numproc,
159	number_of_full_nodes=number_of_full_nodes,
160	number_of_full_triangles=number_of_full_triangles)
161
162	self.set_defaults()
163
164
165	self.forcing_terms.append(manning_friction_implicit)
166	self.forcing_terms.append(gravity)
167
168
169	self.fractional_step_operators = []
170
171	# Stored output
172	self.store = True
173	self.set_store_vertices_uniquely(False)
174	self.quantities_to_be_stored = {'elevation': 1,
175	'stage': 2,
176	'xmomentum': 2,
177	'ymomentum': 2}
178
179
180	def set_defaults(self):
181	"""Set the default values in this routine. That way we can inherit class
182	and just over redefine the defaults for the new class
183	"""
184
185	from anuga.config import minimum_storable_height
186	from anuga.config import minimum_allowed_height, maximum_allowed_speed
187	from anuga.config import g, beta_w, beta_w_dry, \
188	beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
189	from anuga.config import alpha_balance
190	from anuga.config import optimise_dry_cells
191	from anuga.config import optimised_gradient_limiter
192	from anuga.config import use_edge_limiter
193	from anuga.config import use_centroid_velocities
194
195	self.set_minimum_allowed_height(minimum_allowed_height)
196	self.maximum_allowed_speed = maximum_allowed_speed
197
198	self.g = g
199	self.beta_w = beta_w
200	self.beta_w_dry = beta_w_dry
201	self.beta_uh = beta_uh
202	self.beta_uh_dry = beta_uh_dry
203	self.beta_vh = beta_vh
204	self.beta_vh_dry = beta_vh_dry
205	self.alpha_balance = alpha_balance
206
207	self.tight_slope_limiters = tight_slope_limiters
208	self.optimise_dry_cells = optimise_dry_cells
209
210
211	self.set_new_mannings_function(False)
212
213	self.minimum_storable_height = minimum_storable_height
214
215	# Limiters
216	self.use_edge_limiter = use_edge_limiter
217	self.optimised_gradient_limiter = optimised_gradient_limiter
218	self.use_centroid_velocities = use_centroid_velocities
219
220
221	def set_new_mannings_function(self, flag=True):
222	"""Cludge to allow unit test to pass, but to
223	also introduce new mannings friction function
224	which takes into account the slope of the bed.
225	The flag is tested in the python wrapper
226	mannings_friction_implicit
227	"""
228	if flag:
229	self.use_new_mannings = True
230	else:
231	self.use_new_mannings = False
232
233
234	def set_use_edge_limiter(self, flag=True):
235	"""Cludge to allow unit test to pass, but to
236	also introduce new edge limiting. The flag is
237	tested in distribute_to_vertices_and_edges
238	"""
239	if flag:
240	self.use_edge_limiter = True
241	else:
242	self.use_edge_limiter = False
243
244
245	def set_all_betas(self, beta):
246	"""Shorthand to assign one constant value [0,1] to all limiters.
247	0 Corresponds to first order, where as larger values make use of
248	the second order scheme.
249	"""
250
251	self.beta_w = beta
252	self.beta_w_dry = beta
253	self.quantities['stage'].beta = beta
254
255	self.beta_uh = beta
256	self.beta_uh_dry = beta
257	self.quantities['xmomentum'].beta = beta
258
259	self.beta_vh = beta
260	self.beta_vh_dry = beta
261	self.quantities['ymomentum'].beta = beta
262
263
264	def set_store_vertices_uniquely(self, flag, reduction=None):
265	"""Decide whether vertex values should be stored uniquely as
266	computed in the model (True) or whether they should be reduced to one
267	value per vertex using self.reduction (False).
268	"""
269
270	# FIXME (Ole): how about using the word "continuous vertex values" or
271	# "continuous stage surface"
272	self.smooth = not flag
273
274	# Reduction operation for get_vertex_values
275	if reduction is None:
276	self.reduction = mean
277	#self.reduction = min #Looks better near steep slopes
278
279	##
280	# @brief Set the minimum depth that will be written to an SWW file.
281	# @param minimum_storable_height The minimum stored height (in m).
282	def set_minimum_storable_height(self, minimum_storable_height):
283	"""Set the minimum depth that will be recognised when writing
284	to an sww file. This is useful for removing thin water layers
285	that seems to be caused by friction creep.
286
287	The minimum allowed sww depth is in meters.
288	"""
289
290	self.minimum_storable_height = minimum_storable_height
291
292	##
293	# @brief
294	# @param minimum_allowed_height
295	def set_minimum_allowed_height(self, minimum_allowed_height):
296	"""Set minimum depth that will be recognised in the numerical scheme.
297
298	The minimum allowed depth is in meters.
299
300	The parameter H0 (Minimal height for flux computation)
301	is also set by this function
302	"""
303
304	#FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
305
306	#FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
307	#and deal with them adaptively similarly to how we used to use 1 order
308	#steps to recover.
309
310	self.minimum_allowed_height = minimum_allowed_height
311	self.H0 = minimum_allowed_height
312
313	##
314	# @brief
315	# @param maximum_allowed_speed
316	def set_maximum_allowed_speed(self, maximum_allowed_speed):
317	"""Set the maximum particle speed that is allowed in water
318	shallower than minimum_allowed_height. This is useful for
319	controlling speeds in very thin layers of water and at the same time
320	allow some movement avoiding pooling of water.
321	"""
322
323	self.maximum_allowed_speed = maximum_allowed_speed
324
325	##
326	# @brief
327	# @param points_file_block_line_size
328	def set_points_file_block_line_size(self, points_file_block_line_size):
329	"""Set the minimum depth that will be recognised when writing
330	to an sww file. This is useful for removing thin water layers
331	that seems to be caused by friction creep.
332
333	The minimum allowed sww depth is in meters.
334	"""
335	self.points_file_block_line_size = points_file_block_line_size
336
337
338	# FIXME: Probably obsolete in its curren form
339	##
340	# @brief Set the quantities that will be written to an SWW file.
341	# @param q The quantities to be written.
342	# @note Param 'q' may be None, single quantity or list of quantity strings.
343	# @note If 'q' is None, no quantities will be stored in the SWW file.
344	def set_quantities_to_be_stored(self, q):
345	"""Specify which quantities will be stored in the sww file
346
347	q must be either:
348	- a dictionary with quantity names
349	- a list of quantity names (for backwards compatibility)
350	- None
351
352	The format of the dictionary is as follows
353
354	quantity_name: flag where flag must be either 1 or 2.
355	If flag is 1, the quantity is considered static and will be
356	stored once at the beginning of the simulation in a 1D array.
357
358	If flag is 2, the quantity is considered time dependent and
359	it will be stored at each yieldstep by appending it to the
360	appropriate 2D array in the sww file.
361
362	If q is None, storage will be switched off altogether.
363
364	Once the simulation has started and thw sww file opened,
365	this function will have no effect.
366
367	The format, where q is a list of names is for backwards compatibility
368	only.
369	It will take the specified quantities to be time dependent and assume
370	'elevation' to be static regardless.
371	"""
372
373	if q is None:
374	self.quantities_to_be_stored = {}
375	self.store = False
376	return
377
378	# Check correctness
379	for quantity_name in q:
380	msg = ('Quantity %s is not a valid conserved quantity'
381	% quantity_name)
382	assert quantity_name in self.quantities, msg
383
384	assert type(q) == types.DictType
385	self.quantities_to_be_stored = q
386
387	##
388	# @brief
389	# @param indices
390	def get_wet_elements(self, indices=None):
391	"""Return indices for elements where h > minimum_allowed_height
392
393	Optional argument:
394	indices is the set of element ids that the operation applies to.
395
396	Usage:
397	indices = get_wet_elements()
398
399	Note, centroid values are used for this operation
400	"""
401
402	# Water depth below which it is considered to be 0 in the model
403	# FIXME (Ole): Allow this to be specified as a keyword argument as well
404	from anuga.config import minimum_allowed_height
405
406	elevation = self.get_quantity('elevation').\
407	get_values(location='centroids', indices=indices)
408	stage = self.get_quantity('stage').\
409	get_values(location='centroids', indices=indices)
410	depth = stage - elevation
411
412	# Select indices for which depth > 0
413	wet_indices = num.compress(depth > minimum_allowed_height,
414	num.arange(len(depth)))
415	return wet_indices
416
417	##
418	# @brief
419	# @param indices
420	def get_maximum_inundation_elevation(self, indices=None):
421	"""Return highest elevation where h > 0
422
423	Optional argument:
424	indices is the set of element ids that the operation applies to.
425
426	Usage:
427	q = get_maximum_inundation_elevation()
428
429	Note, centroid values are used for this operation
430	"""
431
432	wet_elements = self.get_wet_elements(indices)
433	return self.get_quantity('elevation').\
434	get_maximum_value(indices=wet_elements)
435
436	##
437	# @brief
438	# @param indices
439	def get_maximum_inundation_location(self, indices=None):
440	"""Return location of highest elevation where h > 0
441
442	Optional argument:
443	indices is the set of element ids that the operation applies to.
444
445	Usage:
446	q = get_maximum_inundation_location()
447
448	Note, centroid values are used for this operation
449	"""
450
451	wet_elements = self.get_wet_elements(indices)
452	return self.get_quantity('elevation').\
453	get_maximum_location(indices=wet_elements)
454
455
456	def get_flow_through_cross_section(self, polyline, verbose=False):
457	"""Get the total flow through an arbitrary poly line.
458
459	This is a run-time equivalent of the function with same name
460	in sww_interrogate.py
461
462	Input:
463	polyline: Representation of desired cross section - it may contain
464	multiple sections allowing for complex shapes. Assume
465	absolute UTM coordinates.
466	Format [[x0, y0], [x1, y1], ...]
467
468	Output:
469	Q: Total flow [m^3/s] across given segments.
470	"""
471
472
473	cross_section = Cross_section(self, polyline, verbose)
474
475	return cross_section.get_flow_through_cross_section()
476
477
478	def get_energy_through_cross_section(self, polyline,
479	kind='total',
480	verbose=False):
481	"""Obtain average energy head [m] across specified cross section.
482
483	Inputs:
484	polyline: Representation of desired cross section - it may contain
485	multiple sections allowing for complex shapes. Assume
486	absolute UTM coordinates.
487	Format [[x0, y0], [x1, y1], ...]
488	kind: Select which energy to compute.
489	Options are 'specific' and 'total' (default)
490
491	Output:
492	E: Average energy [m] across given segments for all stored times.
493
494	The average velocity is computed for each triangle intersected by
495	the polyline and averaged weighted by segment lengths.
496
497	The typical usage of this function would be to get average energy of
498	flow in a channel, and the polyline would then be a cross section
499	perpendicular to the flow.
500
501	#FIXME (Ole) - need name for this energy reflecting that its dimension
502	is [m].
503	"""
504
505
506
507	cross_section = Cross_section(self, polyline, verbose)
508
509	return cross_section.get_energy_through_cross_section(kind)
510
511
512	def check_integrity(self):
513	""" Run integrity checks on shallow water domain. """
514	Generic_Domain.check_integrity(self)
515
516	#Check that we are solving the shallow water wave equation
517	msg = 'First conserved quantity must be "stage"'
518	assert self.conserved_quantities[0] == 'stage', msg
519	msg = 'Second conserved quantity must be "xmomentum"'
520	assert self.conserved_quantities[1] == 'xmomentum', msg
521	msg = 'Third conserved quantity must be "ymomentum"'
522	assert self.conserved_quantities[2] == 'ymomentum', msg
523
524	def extrapolate_second_order_sw(self):
525	"""Call correct module function
526	(either from this module or C-extension)"""
527	extrapolate_second_order_sw(self)
528
529	def compute_fluxes(self):
530	"""Call correct module function
531	(either from this module or C-extension)"""
532	compute_fluxes(self)
533
534
535	def apply_fractional_steps(self):
536	"""Loop throughand apply all extra fractional steps"""
537	for operator in self.fractional_step_operators:
538	operator.apply()
539
540	def distribute_to_vertices_and_edges(self):
541	""" Call correct module function """
542	if self.use_edge_limiter:
543	distribute_using_edge_limiter(self)
544	else:
545	distribute_using_vertex_limiter(self)
546
547
548
549	def evolve(self,
550	yieldstep=None,
551	finaltime=None,
552	duration=None,
553	skip_initial_step=False):
554	"""Specialisation of basic evolve method from parent class.
555
556	Evolve the model by 1 step.
557	"""
558
559	# Call check integrity here rather than from user scripts
560	# self.check_integrity()
561
562	msg = 'Attribute self.beta_w must be in the interval [0, 2]'
563	assert 0 <= self.beta_w <= 2.0, msg
564
565	# Initial update of vertex and edge values before any STORAGE
566	# and or visualisation.
567	# This is done again in the initialisation of the Generic_Domain
568	# evolve loop but we do it here to ensure the values are ok for storage.
569	self.distribute_to_vertices_and_edges()
570
571	if self.store is True and self.time == 0.0:
572	self.initialise_storage()
573
574	# Call basic machinery from parent class
575	for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
576	finaltime=finaltime, duration=duration,
577	skip_initial_step=skip_initial_step):
578	# Store model data, e.g. for subsequent visualisation
579	if self.store is True:
580	self.store_timestep()
581
582	# Pass control on to outer loop for more specific actions
583	yield(t)
584
585
586	def initialise_storage(self):
587	"""Create and initialise self.writer object for storing data.
588	Also, save x,y and bed elevation
589	"""
590
591	# Initialise writer
592	self.writer = SWW_file(self)
593
594	# Store vertices and connectivity
595	self.writer.store_connectivity()
596
597
598	def store_timestep(self):
599	"""Store time dependent quantities and time.
600
601	Precondition:
602	self.writer has been initialised
603	"""
604
605	self.writer.store_timestep()
606
607
608	def timestepping_statistics(self,
609	track_speeds=False,
610	triangle_id=None):
611	"""Return string with time stepping statistics for printing or logging
612
613	Optional boolean keyword track_speeds decides whether to report
614	location of smallest timestep as well as a histogram and percentile
615	report.
616	"""
617
618	from anuga.config import epsilon, g
619
620	# Call basic machinery from parent class
621	msg = Generic_Domain.timestepping_statistics(self, track_speeds,
622	triangle_id)
623
624	if track_speeds is True:
625	# qwidth determines the text field used for quantities
626	qwidth = self.qwidth
627
628	# Selected triangle
629	k = self.k
630
631	# Report some derived quantities at vertices, edges and centroid
632	# specific to the shallow water wave equation
633	z = self.quantities['elevation']
634	w = self.quantities['stage']
635
636	Vw = w.get_values(location='vertices', indices=[k])[0]
637	Ew = w.get_values(location='edges', indices=[k])[0]
638	Cw = w.get_values(location='centroids', indices=[k])
639
640	Vz = z.get_values(location='vertices', indices=[k])[0]
641	Ez = z.get_values(location='edges', indices=[k])[0]
642	Cz = z.get_values(location='centroids', indices=[k])
643
644	name = 'depth'
645	Vh = Vw-Vz
646	Eh = Ew-Ez
647	Ch = Cw-Cz
648
649	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\
650	% (name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
651
652	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\
653	% (name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
654
655	message += ' %s: centroid_value = %.4f\n'\
656	% (name.ljust(qwidth), Ch[0])
657
658	msg += message
659
660	uh = self.quantities['xmomentum']
661	vh = self.quantities['ymomentum']
662
663	Vuh = uh.get_values(location='vertices', indices=[k])[0]
664	Euh = uh.get_values(location='edges', indices=[k])[0]
665	Cuh = uh.get_values(location='centroids', indices=[k])
666
667	Vvh = vh.get_values(location='vertices', indices=[k])[0]
668	Evh = vh.get_values(location='edges', indices=[k])[0]
669	Cvh = vh.get_values(location='centroids', indices=[k])
670
671	# Speeds in each direction
672	Vu = Vuh/(Vh + epsilon)
673	Eu = Euh/(Eh + epsilon)
674	Cu = Cuh/(Ch + epsilon)
675	name = 'U'
676	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n' \
677	% (name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
678
679	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n' \
680	% (name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
681
682	message += ' %s: centroid_value = %.4f\n' \
683	% (name.ljust(qwidth), Cu[0])
684
685	msg += message
686
687	Vv = Vvh/(Vh + epsilon)
688	Ev = Evh/(Eh + epsilon)
689	Cv = Cvh/(Ch + epsilon)
690	name = 'V'
691	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n' \
692	% (name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
693
694	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n' \
695	% (name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
696
697	message += ' %s: centroid_value = %.4f\n'\
698	%(name.ljust(qwidth), Cv[0])
699
700	msg += message
701
702	# Froude number in each direction
703	name = 'Froude (x)'
704	Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
705	Efx = Eu/(num.sqrt(g*Eh) + epsilon)
706	Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
707
708	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\
709	% (name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
710
711	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\
712	% (name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
713
714	message += ' %s: centroid_value = %.4f\n'\
715	% (name.ljust(qwidth), Cfx[0])
716
717	msg += message
718
719	name = 'Froude (y)'
720	Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
721	Efy = Ev/(num.sqrt(g*Eh) + epsilon)
722	Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
723
724	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\
725	% (name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
726
727	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\
728	% (name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
729
730	message += ' %s: centroid_value = %.4f\n'\
731	% (name.ljust(qwidth), Cfy[0])
732
733	msg += message
734
735	return msg
736
737
738
739	def compute_boundary_flows(self):
740	"""Compute boundary flows at current timestep.
741
742	Quantities computed are:
743	Total inflow across boundary
744	Total outflow across boundary
745	Flow across each tagged boundary segment
746	"""
747
748	# Run through boundary array and compute for each segment
749	# the normal momentum ((uh, vh) dot normal) times segment length.
750	# Based on sign accumulate this into boundary_inflow and
751	# boundary_outflow.
752
753	# Compute flows along boundary
754
755	uh = self.get_quantity('xmomentum').get_values(location='edges')
756	vh = self.get_quantity('ymomentum').get_values(location='edges')
757
758	# Loop through edges that lie on the boundary and calculate
759	# flows
760	boundary_flows = {}
761	total_boundary_inflow = 0.0
762	total_boundary_outflow = 0.0
763	for vol_id, edge_id in self.boundary:
764	# Compute normal flow across edge. Since normal vector points
765	# away from triangle, a positive sign means that water
766	# flows out from this triangle.
767
768	momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
769	normal = self.mesh.get_normal(vol_id, edge_id)
770	length = self.mesh.get_edgelength(vol_id, edge_id)
771	normal_flow = num.dot(momentum, normal)*length
772
773	# Reverse sign so that + is taken to mean inflow
774	# and - means outflow. This is more intuitive.
775	edge_flow = -normal_flow
776
777	# Tally up inflows and outflows separately
778	if edge_flow > 0:
779	# Flow is inflow
780	total_boundary_inflow += edge_flow
781	else:
782	# Flow is outflow
783	total_boundary_outflow += edge_flow
784
785	# Tally up flows by boundary tag
786	tag = self.boundary[(vol_id, edge_id)]
787
788	if tag not in boundary_flows:
789	boundary_flows[tag] = 0.0
790	boundary_flows[tag] += edge_flow
791
792
793	return boundary_flows, total_boundary_inflow, total_boundary_outflow
794
795
796	def compute_forcing_flows(self):
797	"""Compute flows in and out of domain due to forcing terms.
798
799	Quantities computed are:
800
801
802	Total inflow through forcing terms
803	Total outflow through forcing terms
804	Current total volume in domain
805
806	"""
807
808	#FIXME(Ole): We need to separate what part of explicit_update was
809	# due to the normal flux calculations and what is due to forcing terms.
810
811	pass
812
813
814	def compute_total_volume(self):
815	"""Compute total volume (m^3) of water in entire domain
816	"""
817
818	area = self.mesh.get_areas()
819
820	stage = self.get_quantity('stage').get_values(location='centroids')
821	elevation = \
822	self.get_quantity('elevation').get_values(location='centroids')
823	depth = stage-elevation
824
825	return num.sum(depth*area)
826
827
828	def volumetric_balance_statistics(self):
829	"""Create volumetric balance report suitable for printing or logging.
830	"""
831
832	(boundary_flows, total_boundary_inflow,
833	total_boundary_outflow) = self.compute_boundary_flows()
834
835	message = '---------------------------\n'
836	message += 'Volumetric balance report:\n'
837	message += '--------------------------\n'
838	message += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
839	message += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow
840	message += 'Net boundary flow by tags [m^3/s]\n'
841	for tag in boundary_flows:
842	message += ' %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
843
844	message += 'Total net boundary flow [m^3/s]: %.2f\n' % \
845	(total_boundary_inflow + total_boundary_outflow)
846	message += 'Total volume in domain [m^3]: %.2f\n' % \
847	self.compute_total_volume()
848
849	# The go through explicit forcing update and record the rate of change
850	# for stage and
851	# record into forcing_inflow and forcing_outflow. Finally compute
852	# integral of depth to obtain total volume of domain.
853
854	# FIXME(Ole): This part is not yet done.
855
856	return message
857
858	################################################################################
859	# End of class Shallow Water Domain
860	################################################################################
861
862	#-----------------
863	# Flux computation
864	#-----------------
865
866	## @brief Compute fluxes and timestep suitable for all volumes in domain.
867	# @param domain The domain to calculate fluxes for.
868	def compute_fluxes(domain):
869	"""Compute fluxes and timestep suitable for all volumes in domain.
870
871	Compute total flux for each conserved quantity using "flux_function"
872
873	Fluxes across each edge are scaled by edgelengths and summed up
874	Resulting flux is then scaled by area and stored in
875	explicit_update for each of the three conserved quantities
876	stage, xmomentum and ymomentum
877
878	The maximal allowable speed computed by the flux_function for each volume
879	is converted to a timestep that must not be exceeded. The minimum of
880	those is computed as the next overall timestep.
881
882	Post conditions:
883	domain.explicit_update is reset to computed flux values
884	domain.timestep is set to the largest step satisfying all volumes.
885
886	This wrapper calls the underlying C version of compute fluxes
887	"""
888
889	import sys
890	from shallow_water_ext import compute_fluxes_ext_central \
891	as compute_fluxes_ext
892
893	# Shortcuts
894	Stage = domain.quantities['stage']
895	Xmom = domain.quantities['xmomentum']
896	Ymom = domain.quantities['ymomentum']
897	Bed = domain.quantities['elevation']
898
899	timestep = float(sys.maxint)
900
901	flux_timestep = compute_fluxes_ext(timestep,
902	domain.epsilon,
903	domain.H0,
904	domain.g,
905	domain.neighbours,
906	domain.neighbour_edges,
907	domain.normals,
908	domain.edgelengths,
909	domain.radii,
910	domain.areas,
911	domain.tri_full_flag,
912	Stage.edge_values,
913	Xmom.edge_values,
914	Ymom.edge_values,
915	Bed.edge_values,
916	Stage.boundary_values,
917	Xmom.boundary_values,
918	Ymom.boundary_values,
919	Stage.explicit_update,
920	Xmom.explicit_update,
921	Ymom.explicit_update,
922	domain.already_computed_flux,
923	domain.max_speed,
924	int(domain.optimise_dry_cells))
925
926	domain.flux_timestep = flux_timestep
927
928	################################################################################
929	# Module functions for gradient limiting
930	################################################################################
931
932	##
933	# @brief Wrapper for C version of extrapolate_second_order_sw.
934	# @param domain The domain to operate on.
935	# @note MH090605 The following method belongs to the shallow_water domain class
936	# see comments in the corresponding method in shallow_water_ext.c
937	def extrapolate_second_order_sw(domain):
938	"""Wrapper calling C version of extrapolate_second_order_sw"""
939
940	from shallow_water_ext import extrapolate_second_order_sw as extrapol2
941
942	# Shortcuts
943	Stage = domain.quantities['stage']
944	Xmom = domain.quantities['xmomentum']
945	Ymom = domain.quantities['ymomentum']
946	Elevation = domain.quantities['elevation']
947
948	extrapol2(domain,
949	domain.surrogate_neighbours,
950	domain.number_of_boundaries,
951	domain.centroid_coordinates,
952	Stage.centroid_values,
953	Xmom.centroid_values,
954	Ymom.centroid_values,
955	Elevation.centroid_values,
956	domain.vertex_coordinates,
957	Stage.vertex_values,
958	Xmom.vertex_values,
959	Ymom.vertex_values,
960	Elevation.vertex_values,
961	int(domain.optimise_dry_cells))
962
963	##
964	# @brief Distribution from centroids to vertices specific to the SWW eqn.
965	# @param domain The domain to operate on.
966	def distribute_using_vertex_limiter(domain):
967	"""Distribution from centroids to vertices specific to the SWW equation.
968
969	It will ensure that h (w-z) is always non-negative even in the
970	presence of steep bed-slopes by taking a weighted average between shallow
971	and deep cases.
972
973	In addition, all conserved quantities get distributed as per either a
974	constant (order==1) or a piecewise linear function (order==2).
975
976	FIXME: more explanation about removal of artificial variability etc
977
978	Precondition:
979	All quantities defined at centroids and bed elevation defined at
980	vertices.
981
982	Postcondition
983	Conserved quantities defined at vertices
984	"""
985
986	# Remove very thin layers of water
987	protect_against_infinitesimal_and_negative_heights(domain)
988
989	# Extrapolate all conserved quantities
990	if domain.optimised_gradient_limiter:
991	# MH090605 if second order,
992	# perform the extrapolation and limiting on
993	# all of the conserved quantities
994
995	if (domain._order_ == 1):
996	for name in domain.conserved_quantities:
997	Q = domain.quantities[name]
998	Q.extrapolate_first_order()
999	elif domain._order_ == 2:
1000	domain.extrapolate_second_order_sw()
1001	else:
1002	raise Exception('Unknown order')
1003	else:
1004	# Old code:
1005	for name in domain.conserved_quantities:
1006	Q = domain.quantities[name]
1007
1008	if domain._order_ == 1:
1009	Q.extrapolate_first_order()
1010	elif domain._order_ == 2:
1011	Q.extrapolate_second_order_and_limit_by_vertex()
1012	else:
1013	raise Exception('Unknown order')
1014
1015	# Take bed elevation into account when water heights are small
1016	balance_deep_and_shallow(domain)
1017
1018	# Compute edge values by interpolation
1019	for name in domain.conserved_quantities:
1020	Q = domain.quantities[name]
1021	Q.interpolate_from_vertices_to_edges()
1022
1023	##
1024	# @brief Distribution from centroids to edges specific to the SWW eqn.
1025	# @param domain The domain to operate on.
1026	def distribute_using_edge_limiter(domain):
1027	"""Distribution from centroids to edges specific to the SWW eqn.
1028
1029	It will ensure that h (w-z) is always non-negative even in the
1030	presence of steep bed-slopes by taking a weighted average between shallow
1031	and deep cases.
1032
1033	In addition, all conserved quantities get distributed as per either a
1034	constant (order==1) or a piecewise linear function (order==2).
1035
1036
1037	Precondition:
1038	All quantities defined at centroids and bed elevation defined at
1039	vertices.
1040
1041	Postcondition
1042	Conserved quantities defined at vertices
1043	"""
1044
1045	# Remove very thin layers of water
1046	protect_against_infinitesimal_and_negative_heights(domain)
1047
1048	for name in domain.conserved_quantities:
1049	Q = domain.quantities[name]
1050	if domain._order_ == 1:
1051	Q.extrapolate_first_order()
1052	elif domain._order_ == 2:
1053	Q.extrapolate_second_order_and_limit_by_edge()
1054	else:
1055	raise Exception('Unknown order')
1056
1057	balance_deep_and_shallow(domain)
1058
1059	# Compute edge values by interpolation
1060	for name in domain.conserved_quantities:
1061	Q = domain.quantities[name]
1062	Q.interpolate_from_vertices_to_edges()
1063
1064	##
1065	# @brief Protect against infinitesimal heights and associated high velocities.
1066	# @param domain The domain to operate on.
1067	def protect_against_infinitesimal_and_negative_heights(domain):
1068	"""Protect against infinitesimal heights and associated high velocities"""
1069
1070	from shallow_water_ext import protect
1071
1072	# Shortcuts
1073	wc = domain.quantities['stage'].centroid_values
1074	zc = domain.quantities['elevation'].centroid_values
1075	xmomc = domain.quantities['xmomentum'].centroid_values
1076	ymomc = domain.quantities['ymomentum'].centroid_values
1077
1078	protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1079	domain.epsilon, wc, zc, xmomc, ymomc)
1080
1081	##
1082	# @brief Wrapper for C function balance_deep_and_shallow_c().
1083	# @param domain The domain to operate on.
1084	def balance_deep_and_shallow(domain):
1085	"""Compute linear combination between stage as computed by
1086	gradient-limiters limiting using w, and stage computed by
1087	gradient-limiters limiting using h (h-limiter).
1088	The former takes precedence when heights are large compared to the
1089	bed slope while the latter takes precedence when heights are
1090	relatively small. Anything in between is computed as a balanced
1091	linear combination in order to avoid numerical disturbances which
1092	would otherwise appear as a result of hard switching between
1093	modes.
1094
1095	Wrapper for C implementation
1096	"""
1097
1098	from shallow_water_ext import balance_deep_and_shallow \
1099	as balance_deep_and_shallow_c
1100
1101	# Shortcuts
1102	wc = domain.quantities['stage'].centroid_values
1103	zc = domain.quantities['elevation'].centroid_values
1104	wv = domain.quantities['stage'].vertex_values
1105	zv = domain.quantities['elevation'].vertex_values
1106
1107	# Momentums at centroids
1108	xmomc = domain.quantities['xmomentum'].centroid_values
1109	ymomc = domain.quantities['ymomentum'].centroid_values
1110
1111	# Momentums at vertices
1112	xmomv = domain.quantities['xmomentum'].vertex_values
1113	ymomv = domain.quantities['ymomentum'].vertex_values
1114
1115	balance_deep_and_shallow_c(domain,
1116	wc, zc, wv, zv, wc,
1117	xmomc, ymomc, xmomv, ymomv)
1118
1119
1120
1121	################################################################################
1122	# Standard forcing terms
1123	################################################################################
1124
1125	##
1126	# @brief Apply gravitational pull in the presence of bed slope.
1127	# @param domain The domain to apply gravity to.
1128	# @note Wrapper for C function gravity_c().
1129	def gravity(domain):
1130	"""Apply gravitational pull in the presence of bed slope
1131	Wrapper calls underlying C implementation
1132	"""
1133
1134	from shallow_water_ext import gravity as gravity_c
1135
1136	xmom_update = domain.quantities['xmomentum'].explicit_update
1137	ymom_update = domain.quantities['ymomentum'].explicit_update
1138
1139	stage = domain.quantities['stage']
1140	elevation = domain.quantities['elevation']
1141
1142	#FIXME SR Should avoid allocating memory!
1143	height = stage.centroid_values - elevation.centroid_values
1144	elevation = elevation.vertex_values
1145
1146	point = domain.get_vertex_coordinates()
1147
1148	gravity_c(domain.g, height, elevation, point, xmom_update, ymom_update)
1149
1150	##
1151	# @brief Apply friction to a surface (implicit).
1152	# @param domain The domain to apply Manning friction to.
1153	# @note Wrapper for C function manning_friction_c().
1154	def manning_friction_implicit(domain):
1155	"""Apply (Manning) friction to water momentum
1156	Wrapper for c version
1157	"""
1158
1159	from shallow_water_ext import manning_friction_old
1160	from shallow_water_ext import manning_friction_new
1161
1162	xmom = domain.quantities['xmomentum']
1163	ymom = domain.quantities['ymomentum']
1164
1165	x = domain.get_vertex_coordinates()
1166
1167	w = domain.quantities['stage'].centroid_values
1168	z = domain.quantities['elevation'].vertex_values
1169
1170	uh = xmom.centroid_values
1171	vh = ymom.centroid_values
1172	eta = domain.quantities['friction'].centroid_values
1173
1174	xmom_update = xmom.semi_implicit_update
1175	ymom_update = ymom.semi_implicit_update
1176
1177	eps = domain.minimum_allowed_height
1178	g = domain.g
1179
1180	if domain.use_new_mannings:
1181	manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, \
1182	ymom_update)
1183	else:
1184	manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, \
1185	ymom_update)
1186
1187
1188	##
1189	# @brief Apply friction to a surface (explicit).
1190	# @param domain The domain to apply Manning friction to.
1191	# @note Wrapper for C function manning_friction_c().
1192	def manning_friction_explicit(domain):
1193	"""Apply (Manning) friction to water momentum
1194	Wrapper for c version
1195	"""
1196
1197	from shallow_water_ext import manning_friction_old
1198	from shallow_water_ext import manning_friction_new
1199
1200	xmom = domain.quantities['xmomentum']
1201	ymom = domain.quantities['ymomentum']
1202
1203	x = domain.get_vertex_coordinates()
1204
1205	w = domain.quantities['stage'].centroid_values
1206	z = domain.quantities['elevation'].vertex_values
1207
1208	uh = xmom.centroid_values
1209	vh = ymom.centroid_values
1210	eta = domain.quantities['friction'].centroid_values
1211
1212	xmom_update = xmom.explicit_update
1213	ymom_update = ymom.explicit_update
1214
1215	eps = domain.minimum_allowed_height
1216
1217	if domain.use_new_mannings:
1218	manning_friction_new(domain.g, eps, x, w, uh, vh, z, eta, xmom_update, \
1219	ymom_update)
1220	else:
1221	manning_friction_old(domain.g, eps, w, uh, vh, z, eta, xmom_update, \
1222	ymom_update)
1223
1224
1225
1226	# FIXME (Ole): This was implemented for use with one of the analytical solutions
1227	##
1228	# @brief Apply linear friction to a surface.
1229	# @param domain The domain to apply Manning friction to.
1230	# @note Is this still used (30 Oct 2007)?
1231	def linear_friction(domain):
1232	"""Apply linear friction to water momentum
1233
1234	Assumes quantity: 'linear_friction' to be present
1235	"""
1236
1237	w = domain.quantities['stage'].centroid_values
1238	z = domain.quantities['elevation'].centroid_values
1239	h = w-z
1240
1241	uh = domain.quantities['xmomentum'].centroid_values
1242	vh = domain.quantities['ymomentum'].centroid_values
1243	tau = domain.quantities['linear_friction'].centroid_values
1244
1245	xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1246	ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1247
1248	num_tris = len(domain)
1249	eps = domain.minimum_allowed_height
1250
1251	for k in range(num_tris):
1252	if tau[k] >= eps:
1253	if h[k] >= eps:
1254	S = -tau[k]/h[k]
1255
1256	#Update momentum
1257	xmom_update[k] += S*uh[k]
1258	ymom_update[k] += S*vh[k]
1259
1260	def depth_dependent_friction(domain, default_friction,
1261	surface_roughness_data,
1262	verbose=False):
1263	"""Returns an array of friction values for each wet element adjusted for
1264	depth.
1265
1266	Inputs:
1267	domain - computational domain object
1268	default_friction - depth independent bottom friction
1269	surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values
1270	for each friction region.
1271
1272	Outputs:
1273	wet_friction - Array that can be used directly to update friction as
1274	follows:
1275	domain.set_quantity('friction', wet_friction)
1276
1277
1278
1279	"""
1280
1281	default_n0 = 0 # James - this was missing, don't know what it should be
1282
1283	# Create a temp array to store updated depth dependent
1284	# friction for wet elements
1285	# EHR this is outwardly inneficient but not obvious how to avoid
1286	# recreating each call??????
1287
1288	wet_friction = num.zeros(len(domain), num.float)
1289	wet_friction[:] = default_n0 # Initially assign default_n0 to all array so
1290	# sure have no zeros values
1291
1292	# create depth instance for this timestep
1293	depth = domain.create_quantity_from_expression('stage - elevation')
1294	# Recompute depth as vector
1295	d_vals = depth.get_values(location='centroids')
1296
1297	# rebuild the 'friction' values adjusted for depth at this instant
1298	# loop for each wet element in domain
1299
1300	for i in domain.get_wet_elements():
1301	# Get roughness data for each element
1302	d1 = float(surface_roughness_data[i, 1])
1303	n1 = float(surface_roughness_data[i, 2])
1304	d2 = float(surface_roughness_data[i, 3])
1305	n2 = float(surface_roughness_data[i, 4])
1306
1307
1308	# Recompute friction values from depth for this element
1309
1310	if d_vals[i] <= d1:
1311	ddf = n1
1312	elif d_vals[i] >= d2:
1313	ddf = n2
1314	else:
1315	ddf = n1 + ((n2-n1)/(d2-d1))*(d_vals[i]-d1)
1316
1317	# check sanity of result
1318	if (ddf < 0.010 or \
1319	ddf > 9999.0) :
1320	log.critical('>>>> WARNING: computed depth_dependent friction '
1321	'out of range, ddf%f, n1=%f, n2=%f'
1322	% (ddf, n1, n2))
1323
1324	# update depth dependent friction for that wet element
1325	wet_friction[i] = ddf
1326
1327	# EHR add code to show range of 'friction across domain at this instant as
1328	# sanity check?????????
1329
1330	if verbose :
1331	# return array of domain nvals
1332	nvals = domain.get_quantity('friction').get_values(location='centroids')
1333	n_min = min(nvals)
1334	n_max = max(nvals)
1335
1336	log.critical(' ++++ calculate_depth_dependent_friction - '
1337	'Updated friction - range %7.3f to %7.3f'
1338	% (n_min, n_max))
1339
1340	return wet_friction
1341
1342
1343
1344	################################################################################
1345	# Initialise module
1346	################################################################################
1347
1348	def _raise_compile_exception():
1349	""" Raise exception if compiler not available. """
1350	msg = 'C implementations could not be accessed by %s.\n ' % __file__
1351	msg += 'Make sure compile_all.py has been run as described in '
1352	msg += 'the ANUGA installation guide.'
1353	raise Exception(msg)
1354
1355	from anuga.utilities import compile
1356	if not compile.can_use_C_extension('shallow_water_ext.c'):
1357	_raise_compile_exception()
1358
1359	if __name__ == "__main__":
1360	pass

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

Download in other formats: