Context Navigation

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

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

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

Download in other formats: