Context Navigation

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

Last change on this file since 7967 was 7967, checked in by steve, 13 years ago
Commits from session with Nariman
Property svn:keywords set to `LastChangedDate LastChangedRevision LastChangedBy HeadURL Id`
File size: 47.5 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: 7967 $)
72	ModifiedBy:
73	$Author: steve $
74	$Date: 2010-08-22 23:15:04 +0000 (Sun, 22 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 distribute_to_vertices_and_edges(self):
536	""" Call correct module function """
537	if self.use_edge_limiter:
538	distribute_using_edge_limiter(self)
539	else:
540	distribute_using_vertex_limiter(self)
541
542
543
544	def evolve(self,
545	yieldstep=None,
546	finaltime=None,
547	duration=None,
548	skip_initial_step=False):
549	"""Specialisation of basic evolve method from parent class.
550
551	Evolve the model by 1 step.
552	"""
553
554	# Call check integrity here rather than from user scripts
555	# self.check_integrity()
556
557	msg = 'Attribute self.beta_w must be in the interval [0, 2]'
558	assert 0 <= self.beta_w <= 2.0, msg
559
560	# Initial update of vertex and edge values before any STORAGE
561	# and or visualisation.
562	# This is done again in the initialisation of the Generic_Domain
563	# evolve loop but we do it here to ensure the values are ok for storage.
564	self.distribute_to_vertices_and_edges()
565
566	if self.store is True and self.time == 0.0:
567	self.initialise_storage()
568
569	# Call basic machinery from parent class
570	for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
571	finaltime=finaltime, duration=duration,
572	skip_initial_step=skip_initial_step):
573	# Store model data, e.g. for subsequent visualisation
574	if self.store is True:
575	self.store_timestep()
576
577	# Pass control on to outer loop for more specific actions
578	yield(t)
579
580
581	def initialise_storage(self):
582	"""Create and initialise self.writer object for storing data.
583	Also, save x,y and bed elevation
584	"""
585
586	# Initialise writer
587	self.writer = SWW_file(self)
588
589	# Store vertices and connectivity
590	self.writer.store_connectivity()
591
592
593	def store_timestep(self):
594	"""Store time dependent quantities and time.
595
596	Precondition:
597	self.writer has been initialised
598	"""
599
600	self.writer.store_timestep()
601
602
603	def timestepping_statistics(self,
604	track_speeds=False,
605	triangle_id=None):
606	"""Return string with time stepping statistics for printing or logging
607
608	Optional boolean keyword track_speeds decides whether to report
609	location of smallest timestep as well as a histogram and percentile
610	report.
611	"""
612
613	from anuga.config import epsilon, g
614
615	# Call basic machinery from parent class
616	msg = Generic_Domain.timestepping_statistics(self, track_speeds,
617	triangle_id)
618
619	if track_speeds is True:
620	# qwidth determines the text field used for quantities
621	qwidth = self.qwidth
622
623	# Selected triangle
624	k = self.k
625
626	# Report some derived quantities at vertices, edges and centroid
627	# specific to the shallow water wave equation
628	z = self.quantities['elevation']
629	w = self.quantities['stage']
630
631	Vw = w.get_values(location='vertices', indices=[k])[0]
632	Ew = w.get_values(location='edges', indices=[k])[0]
633	Cw = w.get_values(location='centroids', indices=[k])
634
635	Vz = z.get_values(location='vertices', indices=[k])[0]
636	Ez = z.get_values(location='edges', indices=[k])[0]
637	Cz = z.get_values(location='centroids', indices=[k])
638
639	name = 'depth'
640	Vh = Vw-Vz
641	Eh = Ew-Ez
642	Ch = Cw-Cz
643
644	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\
645	% (name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
646
647	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\
648	% (name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
649
650	message += ' %s: centroid_value = %.4f\n'\
651	% (name.ljust(qwidth), Ch[0])
652
653	msg += message
654
655	uh = self.quantities['xmomentum']
656	vh = self.quantities['ymomentum']
657
658	Vuh = uh.get_values(location='vertices', indices=[k])[0]
659	Euh = uh.get_values(location='edges', indices=[k])[0]
660	Cuh = uh.get_values(location='centroids', indices=[k])
661
662	Vvh = vh.get_values(location='vertices', indices=[k])[0]
663	Evh = vh.get_values(location='edges', indices=[k])[0]
664	Cvh = vh.get_values(location='centroids', indices=[k])
665
666	# Speeds in each direction
667	Vu = Vuh/(Vh + epsilon)
668	Eu = Euh/(Eh + epsilon)
669	Cu = Cuh/(Ch + epsilon)
670	name = 'U'
671	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n' \
672	% (name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
673
674	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n' \
675	% (name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
676
677	message += ' %s: centroid_value = %.4f\n' \
678	% (name.ljust(qwidth), Cu[0])
679
680	msg += message
681
682	Vv = Vvh/(Vh + epsilon)
683	Ev = Evh/(Eh + epsilon)
684	Cv = Cvh/(Ch + epsilon)
685	name = 'V'
686	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n' \
687	% (name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
688
689	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n' \
690	% (name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
691
692	message += ' %s: centroid_value = %.4f\n'\
693	%(name.ljust(qwidth), Cv[0])
694
695	msg += message
696
697	# Froude number in each direction
698	name = 'Froude (x)'
699	Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
700	Efx = Eu/(num.sqrt(g*Eh) + epsilon)
701	Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
702
703	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\
704	% (name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
705
706	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\
707	% (name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
708
709	message += ' %s: centroid_value = %.4f\n'\
710	% (name.ljust(qwidth), Cfx[0])
711
712	msg += message
713
714	name = 'Froude (y)'
715	Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
716	Efy = Ev/(num.sqrt(g*Eh) + epsilon)
717	Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
718
719	message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\
720	% (name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
721
722	message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\
723	% (name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
724
725	message += ' %s: centroid_value = %.4f\n'\
726	% (name.ljust(qwidth), Cfy[0])
727
728	msg += message
729
730	return msg
731
732
733
734	def compute_boundary_flows(self):
735	"""Compute boundary flows at current timestep.
736
737	Quantities computed are:
738	Total inflow across boundary
739	Total outflow across boundary
740	Flow across each tagged boundary segment
741	"""
742
743	# Run through boundary array and compute for each segment
744	# the normal momentum ((uh, vh) dot normal) times segment length.
745	# Based on sign accumulate this into boundary_inflow and
746	# boundary_outflow.
747
748	# Compute flows along boundary
749
750	uh = self.get_quantity('xmomentum').get_values(location='edges')
751	vh = self.get_quantity('ymomentum').get_values(location='edges')
752
753	# Loop through edges that lie on the boundary and calculate
754	# flows
755	boundary_flows = {}
756	total_boundary_inflow = 0.0
757	total_boundary_outflow = 0.0
758	for vol_id, edge_id in self.boundary:
759	# Compute normal flow across edge. Since normal vector points
760	# away from triangle, a positive sign means that water
761	# flows out from this triangle.
762
763	momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
764	normal = self.mesh.get_normal(vol_id, edge_id)
765	length = self.mesh.get_edgelength(vol_id, edge_id)
766	normal_flow = num.dot(momentum, normal)*length
767
768	# Reverse sign so that + is taken to mean inflow
769	# and - means outflow. This is more intuitive.
770	edge_flow = -normal_flow
771
772	# Tally up inflows and outflows separately
773	if edge_flow > 0:
774	# Flow is inflow
775	total_boundary_inflow += edge_flow
776	else:
777	# Flow is outflow
778	total_boundary_outflow += edge_flow
779
780	# Tally up flows by boundary tag
781	tag = self.boundary[(vol_id, edge_id)]
782
783	if tag not in boundary_flows:
784	boundary_flows[tag] = 0.0
785	boundary_flows[tag] += edge_flow
786
787
788	return boundary_flows, total_boundary_inflow, total_boundary_outflow
789
790
791	def compute_forcing_flows(self):
792	"""Compute flows in and out of domain due to forcing terms.
793
794	Quantities computed are:
795
796
797	Total inflow through forcing terms
798	Total outflow through forcing terms
799	Current total volume in domain
800
801	"""
802
803	#FIXME(Ole): We need to separate what part of explicit_update was
804	# due to the normal flux calculations and what is due to forcing terms.
805
806	pass
807
808
809	def compute_total_volume(self):
810	"""Compute total volume (m^3) of water in entire domain
811	"""
812
813	area = self.mesh.get_areas()
814
815	stage = self.get_quantity('stage').get_values(location='centroids')
816	elevation = \
817	self.get_quantity('elevation').get_values(location='centroids')
818	depth = stage-elevation
819
820	return num.sum(depth*area)
821
822
823	def volumetric_balance_statistics(self):
824	"""Create volumetric balance report suitable for printing or logging.
825	"""
826
827	(boundary_flows, total_boundary_inflow,
828	total_boundary_outflow) = self.compute_boundary_flows()
829
830	message = '---------------------------\n'
831	message += 'Volumetric balance report:\n'
832	message += '--------------------------\n'
833	message += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
834	message += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow
835	message += 'Net boundary flow by tags [m^3/s]\n'
836	for tag in boundary_flows:
837	message += ' %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
838
839	message += 'Total net boundary flow [m^3/s]: %.2f\n' % \
840	(total_boundary_inflow + total_boundary_outflow)
841	message += 'Total volume in domain [m^3]: %.2f\n' % \
842	self.compute_total_volume()
843
844	# The go through explicit forcing update and record the rate of change
845	# for stage and
846	# record into forcing_inflow and forcing_outflow. Finally compute
847	# integral of depth to obtain total volume of domain.
848
849	# FIXME(Ole): This part is not yet done.
850
851	return message
852
853	################################################################################
854	# End of class Shallow Water Domain
855	################################################################################
856
857	#-----------------
858	# Flux computation
859	#-----------------
860
861	## @brief Compute fluxes and timestep suitable for all volumes in domain.
862	# @param domain The domain to calculate fluxes for.
863	def compute_fluxes(domain):
864	"""Compute fluxes and timestep suitable for all volumes in domain.
865
866	Compute total flux for each conserved quantity using "flux_function"
867
868	Fluxes across each edge are scaled by edgelengths and summed up
869	Resulting flux is then scaled by area and stored in
870	explicit_update for each of the three conserved quantities
871	stage, xmomentum and ymomentum
872
873	The maximal allowable speed computed by the flux_function for each volume
874	is converted to a timestep that must not be exceeded. The minimum of
875	those is computed as the next overall timestep.
876
877	Post conditions:
878	domain.explicit_update is reset to computed flux values
879	domain.timestep is set to the largest step satisfying all volumes.
880
881	This wrapper calls the underlying C version of compute fluxes
882	"""
883
884	import sys
885	from shallow_water_ext import compute_fluxes_ext_central \
886	as compute_fluxes_ext
887
888	# Shortcuts
889	Stage = domain.quantities['stage']
890	Xmom = domain.quantities['xmomentum']
891	Ymom = domain.quantities['ymomentum']
892	Bed = domain.quantities['elevation']
893
894	timestep = float(sys.maxint)
895
896	flux_timestep = compute_fluxes_ext(timestep,
897	domain.epsilon,
898	domain.H0,
899	domain.g,
900	domain.neighbours,
901	domain.neighbour_edges,
902	domain.normals,
903	domain.edgelengths,
904	domain.radii,
905	domain.areas,
906	domain.tri_full_flag,
907	Stage.edge_values,
908	Xmom.edge_values,
909	Ymom.edge_values,
910	Bed.edge_values,
911	Stage.boundary_values,
912	Xmom.boundary_values,
913	Ymom.boundary_values,
914	Stage.explicit_update,
915	Xmom.explicit_update,
916	Ymom.explicit_update,
917	domain.already_computed_flux,
918	domain.max_speed,
919	int(domain.optimise_dry_cells))
920
921	domain.flux_timestep = flux_timestep
922
923	################################################################################
924	# Module functions for gradient limiting
925	################################################################################
926
927	##
928	# @brief Wrapper for C version of extrapolate_second_order_sw.
929	# @param domain The domain to operate on.
930	# @note MH090605 The following method belongs to the shallow_water domain class
931	# see comments in the corresponding method in shallow_water_ext.c
932	def extrapolate_second_order_sw(domain):
933	"""Wrapper calling C version of extrapolate_second_order_sw"""
934
935	from shallow_water_ext import extrapolate_second_order_sw as extrapol2
936
937	# Shortcuts
938	Stage = domain.quantities['stage']
939	Xmom = domain.quantities['xmomentum']
940	Ymom = domain.quantities['ymomentum']
941	Elevation = domain.quantities['elevation']
942
943	extrapol2(domain,
944	domain.surrogate_neighbours,
945	domain.number_of_boundaries,
946	domain.centroid_coordinates,
947	Stage.centroid_values,
948	Xmom.centroid_values,
949	Ymom.centroid_values,
950	Elevation.centroid_values,
951	domain.vertex_coordinates,
952	Stage.vertex_values,
953	Xmom.vertex_values,
954	Ymom.vertex_values,
955	Elevation.vertex_values,
956	int(domain.optimise_dry_cells))
957
958	##
959	# @brief Distribution from centroids to vertices specific to the SWW eqn.
960	# @param domain The domain to operate on.
961	def distribute_using_vertex_limiter(domain):
962	"""Distribution from centroids to vertices specific to the SWW equation.
963
964	It will ensure that h (w-z) is always non-negative even in the
965	presence of steep bed-slopes by taking a weighted average between shallow
966	and deep cases.
967
968	In addition, all conserved quantities get distributed as per either a
969	constant (order==1) or a piecewise linear function (order==2).
970
971	FIXME: more explanation about removal of artificial variability etc
972
973	Precondition:
974	All quantities defined at centroids and bed elevation defined at
975	vertices.
976
977	Postcondition
978	Conserved quantities defined at vertices
979	"""
980
981	# Remove very thin layers of water
982	protect_against_infinitesimal_and_negative_heights(domain)
983
984	# Extrapolate all conserved quantities
985	if domain.optimised_gradient_limiter:
986	# MH090605 if second order,
987	# perform the extrapolation and limiting on
988	# all of the conserved quantities
989
990	if (domain._order_ == 1):
991	for name in domain.conserved_quantities:
992	Q = domain.quantities[name]
993	Q.extrapolate_first_order()
994	elif domain._order_ == 2:
995	domain.extrapolate_second_order_sw()
996	else:
997	raise Exception('Unknown order')
998	else:
999	# Old code:
1000	for name in domain.conserved_quantities:
1001	Q = domain.quantities[name]
1002
1003	if domain._order_ == 1:
1004	Q.extrapolate_first_order()
1005	elif domain._order_ == 2:
1006	Q.extrapolate_second_order_and_limit_by_vertex()
1007	else:
1008	raise Exception('Unknown order')
1009
1010	# Take bed elevation into account when water heights are small
1011	balance_deep_and_shallow(domain)
1012
1013	# Compute edge values by interpolation
1014	for name in domain.conserved_quantities:
1015	Q = domain.quantities[name]
1016	Q.interpolate_from_vertices_to_edges()
1017
1018	##
1019	# @brief Distribution from centroids to edges specific to the SWW eqn.
1020	# @param domain The domain to operate on.
1021	def distribute_using_edge_limiter(domain):
1022	"""Distribution from centroids to edges specific to the SWW eqn.
1023
1024	It will ensure that h (w-z) is always non-negative even in the
1025	presence of steep bed-slopes by taking a weighted average between shallow
1026	and deep cases.
1027
1028	In addition, all conserved quantities get distributed as per either a
1029	constant (order==1) or a piecewise linear function (order==2).
1030
1031
1032	Precondition:
1033	All quantities defined at centroids and bed elevation defined at
1034	vertices.
1035
1036	Postcondition
1037	Conserved quantities defined at vertices
1038	"""
1039
1040	# Remove very thin layers of water
1041	protect_against_infinitesimal_and_negative_heights(domain)
1042
1043	for name in domain.conserved_quantities:
1044	Q = domain.quantities[name]
1045	if domain._order_ == 1:
1046	Q.extrapolate_first_order()
1047	elif domain._order_ == 2:
1048	Q.extrapolate_second_order_and_limit_by_edge()
1049	else:
1050	raise Exception('Unknown order')
1051
1052	balance_deep_and_shallow(domain)
1053
1054	# Compute edge values by interpolation
1055	for name in domain.conserved_quantities:
1056	Q = domain.quantities[name]
1057	Q.interpolate_from_vertices_to_edges()
1058
1059	##
1060	# @brief Protect against infinitesimal heights and associated high velocities.
1061	# @param domain The domain to operate on.
1062	def protect_against_infinitesimal_and_negative_heights(domain):
1063	"""Protect against infinitesimal heights and associated high velocities"""
1064
1065	from shallow_water_ext import protect
1066
1067	# Shortcuts
1068	wc = domain.quantities['stage'].centroid_values
1069	zc = domain.quantities['elevation'].centroid_values
1070	xmomc = domain.quantities['xmomentum'].centroid_values
1071	ymomc = domain.quantities['ymomentum'].centroid_values
1072
1073	protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1074	domain.epsilon, wc, zc, xmomc, ymomc)
1075
1076	##
1077	# @brief Wrapper for C function balance_deep_and_shallow_c().
1078	# @param domain The domain to operate on.
1079	def balance_deep_and_shallow(domain):
1080	"""Compute linear combination between stage as computed by
1081	gradient-limiters limiting using w, and stage computed by
1082	gradient-limiters limiting using h (h-limiter).
1083	The former takes precedence when heights are large compared to the
1084	bed slope while the latter takes precedence when heights are
1085	relatively small. Anything in between is computed as a balanced
1086	linear combination in order to avoid numerical disturbances which
1087	would otherwise appear as a result of hard switching between
1088	modes.
1089
1090	Wrapper for C implementation
1091	"""
1092
1093	from shallow_water_ext import balance_deep_and_shallow \
1094	as balance_deep_and_shallow_c
1095
1096	# Shortcuts
1097	wc = domain.quantities['stage'].centroid_values
1098	zc = domain.quantities['elevation'].centroid_values
1099	wv = domain.quantities['stage'].vertex_values
1100	zv = domain.quantities['elevation'].vertex_values
1101
1102	# Momentums at centroids
1103	xmomc = domain.quantities['xmomentum'].centroid_values
1104	ymomc = domain.quantities['ymomentum'].centroid_values
1105
1106	# Momentums at vertices
1107	xmomv = domain.quantities['xmomentum'].vertex_values
1108	ymomv = domain.quantities['ymomentum'].vertex_values
1109
1110	balance_deep_and_shallow_c(domain,
1111	wc, zc, wv, zv, wc,
1112	xmomc, ymomc, xmomv, ymomv)
1113
1114
1115
1116	################################################################################
1117	# Standard forcing terms
1118	################################################################################
1119
1120	##
1121	# @brief Apply gravitational pull in the presence of bed slope.
1122	# @param domain The domain to apply gravity to.
1123	# @note Wrapper for C function gravity_c().
1124	def gravity(domain):
1125	"""Apply gravitational pull in the presence of bed slope
1126	Wrapper calls underlying C implementation
1127	"""
1128
1129	from shallow_water_ext import gravity as gravity_c
1130
1131	xmom_update = domain.quantities['xmomentum'].explicit_update
1132	ymom_update = domain.quantities['ymomentum'].explicit_update
1133
1134	stage = domain.quantities['stage']
1135	elevation = domain.quantities['elevation']
1136
1137	#FIXME SR Should avoid allocating memory!
1138	height = stage.centroid_values - elevation.centroid_values
1139	elevation = elevation.vertex_values
1140
1141	point = domain.get_vertex_coordinates()
1142
1143	gravity_c(domain.g, height, elevation, point, xmom_update, ymom_update)
1144
1145	##
1146	# @brief Apply friction to a surface (implicit).
1147	# @param domain The domain to apply Manning friction to.
1148	# @note Wrapper for C function manning_friction_c().
1149	def manning_friction_implicit(domain):
1150	"""Apply (Manning) friction to water momentum
1151	Wrapper for c version
1152	"""
1153
1154	from shallow_water_ext import manning_friction_old
1155	from shallow_water_ext import manning_friction_new
1156
1157	xmom = domain.quantities['xmomentum']
1158	ymom = domain.quantities['ymomentum']
1159
1160	x = domain.get_vertex_coordinates()
1161
1162	w = domain.quantities['stage'].centroid_values
1163	z = domain.quantities['elevation'].vertex_values
1164
1165	uh = xmom.centroid_values
1166	vh = ymom.centroid_values
1167	eta = domain.quantities['friction'].centroid_values
1168
1169	xmom_update = xmom.semi_implicit_update
1170	ymom_update = ymom.semi_implicit_update
1171
1172	eps = domain.minimum_allowed_height
1173	g = domain.g
1174
1175	if domain.use_new_mannings:
1176	manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, \
1177	ymom_update)
1178	else:
1179	manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, \
1180	ymom_update)
1181
1182
1183	##
1184	# @brief Apply friction to a surface (explicit).
1185	# @param domain The domain to apply Manning friction to.
1186	# @note Wrapper for C function manning_friction_c().
1187	def manning_friction_explicit(domain):
1188	"""Apply (Manning) friction to water momentum
1189	Wrapper for c version
1190	"""
1191
1192	from shallow_water_ext import manning_friction_old
1193	from shallow_water_ext import manning_friction_new
1194
1195	xmom = domain.quantities['xmomentum']
1196	ymom = domain.quantities['ymomentum']
1197
1198	x = domain.get_vertex_coordinates()
1199
1200	w = domain.quantities['stage'].centroid_values
1201	z = domain.quantities['elevation'].vertex_values
1202
1203	uh = xmom.centroid_values
1204	vh = ymom.centroid_values
1205	eta = domain.quantities['friction'].centroid_values
1206
1207	xmom_update = xmom.explicit_update
1208	ymom_update = ymom.explicit_update
1209
1210	eps = domain.minimum_allowed_height
1211
1212	if domain.use_new_mannings:
1213	manning_friction_new(domain.g, eps, x, w, uh, vh, z, eta, xmom_update, \
1214	ymom_update)
1215	else:
1216	manning_friction_old(domain.g, eps, w, uh, vh, z, eta, xmom_update, \
1217	ymom_update)
1218
1219
1220
1221	# FIXME (Ole): This was implemented for use with one of the analytical solutions
1222	##
1223	# @brief Apply linear friction to a surface.
1224	# @param domain The domain to apply Manning friction to.
1225	# @note Is this still used (30 Oct 2007)?
1226	def linear_friction(domain):
1227	"""Apply linear friction to water momentum
1228
1229	Assumes quantity: 'linear_friction' to be present
1230	"""
1231
1232	w = domain.quantities['stage'].centroid_values
1233	z = domain.quantities['elevation'].centroid_values
1234	h = w-z
1235
1236	uh = domain.quantities['xmomentum'].centroid_values
1237	vh = domain.quantities['ymomentum'].centroid_values
1238	tau = domain.quantities['linear_friction'].centroid_values
1239
1240	xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1241	ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1242
1243	num_tris = len(domain)
1244	eps = domain.minimum_allowed_height
1245
1246	for k in range(num_tris):
1247	if tau[k] >= eps:
1248	if h[k] >= eps:
1249	S = -tau[k]/h[k]
1250
1251	#Update momentum
1252	xmom_update[k] += S*uh[k]
1253	ymom_update[k] += S*vh[k]
1254
1255	def depth_dependent_friction(domain, default_friction,
1256	surface_roughness_data,
1257	verbose=False):
1258	"""Returns an array of friction values for each wet element adjusted for
1259	depth.
1260
1261	Inputs:
1262	domain - computational domain object
1263	default_friction - depth independent bottom friction
1264	surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values
1265	for each friction region.
1266
1267	Outputs:
1268	wet_friction - Array that can be used directly to update friction as
1269	follows:
1270	domain.set_quantity('friction', wet_friction)
1271
1272
1273
1274	"""
1275
1276	default_n0 = 0 # James - this was missing, don't know what it should be
1277
1278	# Create a temp array to store updated depth dependent
1279	# friction for wet elements
1280	# EHR this is outwardly inneficient but not obvious how to avoid
1281	# recreating each call??????
1282
1283	wet_friction = num.zeros(len(domain), num.float)
1284	wet_friction[:] = default_n0 # Initially assign default_n0 to all array so
1285	# sure have no zeros values
1286
1287	# create depth instance for this timestep
1288	depth = domain.create_quantity_from_expression('stage - elevation')
1289	# Recompute depth as vector
1290	d_vals = depth.get_values(location='centroids')
1291
1292	# rebuild the 'friction' values adjusted for depth at this instant
1293	# loop for each wet element in domain
1294
1295	for i in domain.get_wet_elements():
1296	# Get roughness data for each element
1297	d1 = float(surface_roughness_data[i, 1])
1298	n1 = float(surface_roughness_data[i, 2])
1299	d2 = float(surface_roughness_data[i, 3])
1300	n2 = float(surface_roughness_data[i, 4])
1301
1302
1303	# Recompute friction values from depth for this element
1304
1305	if d_vals[i] <= d1:
1306	ddf = n1
1307	elif d_vals[i] >= d2:
1308	ddf = n2
1309	else:
1310	ddf = n1 + ((n2-n1)/(d2-d1))*(d_vals[i]-d1)
1311
1312	# check sanity of result
1313	if (ddf < 0.010 or \
1314	ddf > 9999.0) :
1315	log.critical('>>>> WARNING: computed depth_dependent friction '
1316	'out of range, ddf%f, n1=%f, n2=%f'
1317	% (ddf, n1, n2))
1318
1319	# update depth dependent friction for that wet element
1320	wet_friction[i] = ddf
1321
1322	# EHR add code to show range of 'friction across domain at this instant as
1323	# sanity check?????????
1324
1325	if verbose :
1326	# return array of domain nvals
1327	nvals = domain.get_quantity('friction').get_values(location='centroids')
1328	n_min = min(nvals)
1329	n_max = max(nvals)
1330
1331	log.critical(' ++++ calculate_depth_dependent_friction - '
1332	'Updated friction - range %7.3f to %7.3f'
1333	% (n_min, n_max))
1334
1335	return wet_friction
1336
1337
1338
1339	################################################################################
1340	# Initialise module
1341	################################################################################
1342
1343	def _raise_compile_exception():
1344	""" Raise exception if compiler not available. """
1345	msg = 'C implementations could not be accessed by %s.\n ' % __file__
1346	msg += 'Make sure compile_all.py has been run as described in '
1347	msg += 'the ANUGA installation guide.'
1348	raise Exception(msg)
1349
1350	from anuga.utilities import compile
1351	if not compile.can_use_C_extension('shallow_water_ext.c'):
1352	_raise_compile_exception()
1353
1354	if __name__ == "__main__":
1355	pass

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

Download in other formats: