Context Navigation

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

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

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

Download in other formats: