source: anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 3926

Last change on this file since 3926 was 3926, checked in by ole, 17 years ago

First step towards keeping track of full nodes and triangles in
parallel domains.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 63.8 KB
RevLine 
[3804]1"""Class Domain -
22D triangular domains for finite-volume computations of
3the shallow water wave equation.
4
5
6$Description:
7This module contains a specialisation of class Domain from module domain.py
8consisting of methods specific to the Shallow Water Wave Equation
9
10
11U_t + E_x + G_y = S
12
13where
14
15U = [w, uh, vh]
16E = [uh, u^2h + gh^2/2, uvh]
17G = [vh, uvh, v^2h + gh^2/2]
18S represents source terms forcing the system
19(e.g. gravity, friction, wind stress, ...)
20
21and _t, _x, _y denote the derivative with respect to t, x and y respectively.
22
23The quantities are
24
25symbol    variable name    explanation
26x         x                horizontal distance from origin [m]
27y         y                vertical distance from origin [m]
28z         elevation        elevation of bed on which flow is modelled [m]
29h         height           water height above z [m]
30w         stage            absolute water level, w = z+h [m]
31u                          speed in the x direction [m/s]
32v                          speed in the y direction [m/s]
33uh        xmomentum        momentum in the x direction [m^2/s]
34vh        ymomentum        momentum in the y direction [m^2/s]
35
36eta                        mannings friction coefficient [to appear]
37nu                         wind stress coefficient [to appear]
38
39The conserved quantities are w, uh, vh
40
41$References
42Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
43Christopher Zoppou and Stephen Roberts,
44Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
45
46Hydrodynamic modelling of coastal inundation.
47Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
48In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
49Modelling and Simulation. Modelling and Simulation Society of Australia and
50New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.
51http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
52
53
54Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
55Geoscience Australia, 2004
56
57
58$Author: ole $
59$Revision: 3926 $
60$Date: 2006-11-06 23:57:58 +0000 (Mon, 06 Nov 2006) $
61$LastChangedDate: 2006-11-06 23:57:58 +0000 (Mon, 06 Nov 2006) $
62$LastChangedRevision: 3926 $
63$LastChangedBy: ole $
64$HeadURL: anuga_core/source/anuga/shallow_water/shallow_water_domain.py $
65"""
66
67#Subversion keywords:
68#
69#$LastChangedDate: 2006-11-06 23:57:58 +0000 (Mon, 06 Nov 2006) $
70#$LastChangedRevision: 3926 $
71#$LastChangedBy: ole $
72
73from Numeric import zeros, ones, Float, array, sum, size
74from Numeric import compress, arange
75
76
77from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain
78from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
79     import Boundary
80from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
81     import File_boundary
82from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
83     import Dirichlet_boundary
84from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
85     import Time_boundary
86from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
87     import Transmissive_boundary
88
89from anuga.utilities.numerical_tools import gradient, mean
90from anuga.config import minimum_storable_height
91from anuga.config import minimum_allowed_height, maximum_allowed_speed
92from anuga.config import g, beta_h, beta_w, beta_w_dry,\
93     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry
[3876]94from anuga.config import alpha_balance
[3804]95
96
97#Shallow water domain
98class Domain(Generic_Domain):
99
100    def __init__(self,
101                 coordinates=None,
102                 vertices=None,
103                 boundary=None,
104                 tagged_elements=None,
105                 geo_reference=None,
106                 use_inscribed_circle=False,
107                 mesh_filename=None,
108                 use_cache=False,
109                 verbose=False,
110                 full_send_dict=None,
111                 ghost_recv_dict=None,
112                 processor=0,
[3926]113                 numproc=1,
114                 number_of_full_nodes=0,
115                 number_of_full_triangles=0):
[3804]116
117
118        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
119        other_quantities = ['elevation', 'friction']
120        Generic_Domain.__init__(self,
121                                coordinates,
122                                vertices,
123                                boundary,
124                                conserved_quantities,
125                                other_quantities,
126                                tagged_elements,
127                                geo_reference,
128                                use_inscribed_circle,
129                                mesh_filename,
130                                use_cache,
131                                verbose,
132                                full_send_dict,
133                                ghost_recv_dict,
134                                processor,
[3926]135                                numproc,
136                                number_of_full_nodes=number_of_full_nodes,
137                                number_of_full_triangles=number_of_full_triangles) 
[3804]138
139        self.minimum_allowed_height = minimum_allowed_height
140        self.maximum_allowed_speed = maximum_allowed_speed
141        self.g = g
142        self.beta_w      = beta_w
143        self.beta_w_dry  = beta_w_dry
144        self.beta_uh     = beta_uh
145        self.beta_uh_dry = beta_uh_dry
146        self.beta_vh     = beta_vh
147        self.beta_vh_dry = beta_vh_dry
148        self.beta_h      = beta_h
[3876]149        self.alpha_balance = alpha_balance
[3804]150
151        self.flux_function = flux_function_central
152        #self.flux_function = flux_function_kinetic
153       
154        self.forcing_terms.append(manning_friction)
155        self.forcing_terms.append(gravity)
156
157        #Realtime visualisation
158        self.visualiser = None
159        self.visualise  = False
160        self.visualise_color_stage = False
161        self.visualise_stage_range = 1.0
162        self.visualise_timer = True
163        self.visualise_range_z = None
164
165        #Stored output
166        self.store = True
167        self.format = 'sww'
168        self.set_store_vertices_uniquely(False)
169        self.minimum_storable_height = minimum_storable_height
170        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
171               
172
[3848]173    def set_all_limiters(self, beta):
[3847]174        """Shorthand to assign one constant value [0,1[ to all limiters.
175        0 Corresponds to first order, where as larger values make use of
176        the second order scheme.
177        """
[3804]178
[3847]179        self.beta_w      = beta
180        self.beta_w_dry  = beta
181        self.beta_uh     = beta
182        self.beta_uh_dry = beta
183        self.beta_vh     = beta
184        self.beta_vh_dry = beta
185        self.beta_h      = beta
186       
187
[3804]188    def set_store_vertices_uniquely(self, flag, reduction=None):
189        """Decide whether vertex values should be stored uniquely as
190        computed in the model or whether they should be reduced to one
191        value per vertex using self.reduction.
192        """
193        self.smooth = not flag
194
195        #Reduction operation for get_vertex_values
196        if reduction is None:
197            self.reduction = mean
198            #self.reduction = min  #Looks better near steep slopes
199
200
201    def set_minimum_storable_height(self, minimum_storable_height):
202        """
203        Set the minimum depth that will be recognised when writing
204        to an sww file. This is useful for removing thin water layers
205        that seems to be caused by friction creep.
206
207        The minimum allowed sww depth is in meters.
208        """
209        self.minimum_storable_height = minimum_storable_height
210       
211
212    def set_maximum_allowed_speed(self, maximum_allowed_speed):
213        """
214        Set the maximum particle speed that is allowed in water
215        shallower than minimum_allowed_height. This is useful for
216        controlling speeds in very thin layers of water and at the same time
217        allow some movement avoiding pooling of water.
218
219        """
220        self.maximum_allowed_speed = maximum_allowed_speed
221
222
223    def set_quantities_to_be_stored(self, q):
224        """Specify which quantities will be stored in the sww file.
225
226        q must be either:
227          - the name of a quantity
228          - a list of quantity names
229          - None
230
231        In the two first cases, the named quantities will be stored at
232        each yieldstep (This is in addition to the quantities elevation
233        and friction)
234       
235        If q is None, storage will be switched off altogether.
236        """
237
238
239        if q is None:
240            self.quantities_to_be_stored = []
241            self.store = False
242            return
243
244        if isinstance(q, basestring):
245            q = [q] # Turn argument into a list
246
247        #Check correcness
248        for quantity_name in q:
249            msg = 'Quantity %s is not a valid conserved quantity'\
250                  %quantity_name
251           
252            assert quantity_name in self.conserved_quantities, msg
253
254        self.quantities_to_be_stored = q
255
256
257
258    def get_wet_elements(self, indices=None):
259        """Return indices for elements where h > minimum_allowed_height
260
261        Optional argument:
262            indices is the set of element ids that the operation applies to.
263
264        Usage:
265            indices = get_wet_elements()
266
267        Note, centroid values are used for this operation           
268        """
269
270        # Water depth below which it is considered to be 0 in the model
271        # FIXME (Ole): Allow this to be specified as a keyword argument as well
272        from anuga.config import minimum_allowed_height
273
274       
275        elevation = self.get_quantity('elevation').\
276                    get_values(location='centroids', indices=indices)
277        stage = self.get_quantity('stage').\
278                get_values(location='centroids', indices=indices)       
279        depth = stage - elevation
280
281        # Select indices for which depth > 0
282        wet_indices = compress(depth > minimum_allowed_height,
283                               arange(len(depth)))
284        return wet_indices
285
286
287    def get_maximum_inundation_elevation(self, indices=None):
288        """Return highest elevation where h > 0
289
290        Optional argument:
291            indices is the set of element ids that the operation applies to.
292
293        Usage:
294            q = get_maximum_inundation_elevation()
295
296        Note, centroid values are used for this operation           
297        """
298
299        wet_elements = self.get_wet_elements(indices)
300        return self.get_quantity('elevation').\
301               get_maximum_value(indices=wet_elements)
302
303
304    def get_maximum_inundation_location(self, indices=None):
305        """Return highest elevation where h > 0
306
307        Optional argument:
308            indices is the set of element ids that the operation applies to.
309
310        Usage:
311            q = get_maximum_inundation_elevation()
312
313        Note, centroid values are used for this operation           
314        """
315
316        wet_elements = self.get_wet_elements(indices)
317        return self.get_quantity('elevation').\
318               get_maximum_location(indices=wet_elements)   
319
320
321
322
323    def initialise_visualiser(self,scale_z=1.0,rect=None):
324        #Realtime visualisation
325        if self.visualiser is None:
326            from realtime_visualisation_new import Visualiser
327            self.visualiser = Visualiser(self,scale_z,rect)
328            self.visualiser.setup['elevation']=True
329            self.visualiser.updating['stage']=True
330        self.visualise = True
331        if self.visualise_color_stage == True:
332            self.visualiser.coloring['stage'] = True
333            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
334        print 'initialise visualiser'
335        print self.visualiser.setup
336        print self.visualiser.updating
337
338    def check_integrity(self):
339        Generic_Domain.check_integrity(self)
340
341        #Check that we are solving the shallow water wave equation
342
343        msg = 'First conserved quantity must be "stage"'
344        assert self.conserved_quantities[0] == 'stage', msg
345        msg = 'Second conserved quantity must be "xmomentum"'
346        assert self.conserved_quantities[1] == 'xmomentum', msg
347        msg = 'Third conserved quantity must be "ymomentum"'
348        assert self.conserved_quantities[2] == 'ymomentum', msg
349
350    def extrapolate_second_order_sw(self):
351        #Call correct module function
352        #(either from this module or C-extension)
353        extrapolate_second_order_sw(self)
354
355    def compute_fluxes(self):
356        #Call correct module function
357        #(either from this module or C-extension)
358        compute_fluxes(self)
359
360    def distribute_to_vertices_and_edges(self):
361        #Call correct module function
362        #(either from this module or C-extension)
363        distribute_to_vertices_and_edges(self)
364
365
366    #FIXME: Under construction
367#     def set_defaults(self):
368#         """Set default values for uninitialised quantities.
369#         This is specific to the shallow water wave equation
370#         Defaults for 'elevation', 'friction', 'xmomentum' and 'ymomentum'
371#         are 0.0. Default for 'stage' is whatever the value of 'elevation'.
372#         """
373
374#         for name in self.other_quantities + self.conserved_quantities:
375#             print name
376#             print self.quantities.keys()
377#             if not self.quantities.has_key(name):
378#                 if name == 'stage':
379
380#                     if self.quantities.has_key('elevation'):
381#                         z = self.quantities['elevation'].vertex_values
382#                         self.set_quantity(name, z)
383#                     else:
384#                         self.set_quantity(name, 0.0)
385#                 else:
386#                     self.set_quantity(name, 0.0)
387
388
389
390#         #Lift negative heights up
391#         #z = self.quantities['elevation'].vertex_values
392#         #w = self.quantities['stage'].vertex_values
393
394#         #h = w-z
395
396#         #for k in range(h.shape[0]):
397#         #    for i in range(3):
398#         #        if h[k, i] < 0.0:
399#         #            w[k, i] = z[k, i]
400
401
402#         #self.quantities['stage'].interpolate()
403
404
405    def evolve(self,
406               yieldstep = None,
407               finaltime = None,
408               duration = None,
409               skip_initial_step = False):
410        """Specialisation of basic evolve method from parent class
411        """
412
413        #Call check integrity here rather than from user scripts
414        #self.check_integrity()
415
416        msg = 'Parameter beta_h must be in the interval [0, 1['
417        assert 0 <= self.beta_h <= 1.0, msg
418        msg = 'Parameter beta_w must be in the interval [0, 1['
419        assert 0 <= self.beta_w <= 1.0, msg
420
421
422        #Initial update of vertex and edge values before any storage
423        #and or visualisation
424        self.distribute_to_vertices_and_edges()
425
426        #Initialise real time viz if requested
427        if self.visualise is True and self.time == 0.0:
428            if self.visualiser is None:
429                self.initialise_visualiser()
430
431            self.visualiser.update_timer()
432            self.visualiser.setup_all()
433
434        #Store model data, e.g. for visualisation
435        if self.store is True and self.time == 0.0:
436            self.initialise_storage()
437            #print 'Storing results in ' + self.writer.filename
438        else:
439            pass
440            #print 'Results will not be stored.'
441            #print 'To store results set domain.store = True'
442            #FIXME: Diagnostic output should be controlled by
443            # a 'verbose' flag living in domain (or in a parent class)
444
445        #Call basic machinery from parent class
446        for t in Generic_Domain.evolve(self,
447                                       yieldstep=yieldstep,
448                                       finaltime=finaltime,
449                                       duration=duration,
450                                       skip_initial_step=skip_initial_step):
451            #Real time viz
452            if self.visualise is True:
453                self.visualiser.update_all()
454                self.visualiser.update_timer()
455
456
457            #Store model data, e.g. for subsequent visualisation
458            if self.store is True:
459                self.store_timestep(self.quantities_to_be_stored)
460
461            #FIXME: Could maybe be taken from specified list
462            #of 'store every step' quantities
463
464            #Pass control on to outer loop for more specific actions
465            yield(t)
466
467    def initialise_storage(self):
468        """Create and initialise self.writer object for storing data.
469        Also, save x,y and bed elevation
470        """
471
472        from anuga.shallow_water.data_manager import get_dataobject
473
474        #Initialise writer
475        self.writer = get_dataobject(self, mode = 'w')
476
477        #Store vertices and connectivity
478        self.writer.store_connectivity()
479
480
481    def store_timestep(self, name):
482        """Store named quantity and time.
483
484        Precondition:
485           self.write has been initialised
486        """
487        self.writer.store_timestep(name)
488
489
490#=============== End of Shallow Water Domain ===============================
491
492
493
494#Rotation of momentum vector
495def rotate(q, normal, direction = 1):
496    """Rotate the momentum component q (q[1], q[2])
497    from x,y coordinates to coordinates based on normal vector.
498
499    If direction is negative the rotation is inverted.
500
501    Input vector is preserved
502
503    This function is specific to the shallow water wave equation
504    """
505
506    assert len(q) == 3,\
507           'Vector of conserved quantities must have length 3'\
508           'for 2D shallow water equation'
509
510    try:
511        l = len(normal)
512    except:
513        raise 'Normal vector must be an Numeric array'
514
515    assert l == 2, 'Normal vector must have 2 components'
516
517
518    n1 = normal[0]
519    n2 = normal[1]
520
521    r = zeros(len(q), Float) #Rotated quantities
522    r[0] = q[0]              #First quantity, height, is not rotated
523
524    if direction == -1:
525        n2 = -n2
526
527
528    r[1] =  n1*q[1] + n2*q[2]
529    r[2] = -n2*q[1] + n1*q[2]
530
531    return r
532
533
534
535####################################
536# Flux computation
537def flux_function_central(normal, ql, qr, zl, zr):
538    """Compute fluxes between volumes for the shallow water wave equation
539    cast in terms of w = h+z using the 'central scheme' as described in
540
541    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
542    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
543    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
544
545    The implemented formula is given in equation (3.15) on page 714
546
547    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
548    in the numerical vectors ql and qr.
549
550    Bed elevations zl and zr.
551    """
552
553    from anuga.config import g, epsilon
554    from math import sqrt
555
556    #Align momentums with x-axis
557    q_left  = rotate(ql, normal, direction = 1)
558    q_right = rotate(qr, normal, direction = 1)
559
560    z = (zl+zr)/2 #Take average of field values
561
562    w_left  = q_left[0]   #w=h+z
563    h_left  = w_left-z
564    uh_left = q_left[1]
565
566    if h_left < epsilon:
567        u_left = 0.0  #Could have been negative
568        h_left = 0.0
569    else:
570        u_left  = uh_left/h_left
571
572
573    w_right  = q_right[0]  #w=h+z
574    h_right  = w_right-z
575    uh_right = q_right[1]
576
577
578    if h_right < epsilon:
579        u_right = 0.0  #Could have been negative
580        h_right = 0.0
581    else:
582        u_right  = uh_right/h_right
583
584    vh_left  = q_left[2]
585    vh_right = q_right[2]
586
587    soundspeed_left  = sqrt(g*h_left)
588    soundspeed_right = sqrt(g*h_right)
589
590    #Maximal wave speed
591    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
592
593    #Minimal wave speed
594    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
595
596    #Flux computation
597
598    #FIXME(Ole): Why is it again that we don't
599    #use uh_left and uh_right directly in the first entries?
600    flux_left  = array([u_left*h_left,
601                        u_left*uh_left + 0.5*g*h_left**2,
602                        u_left*vh_left])
603    flux_right = array([u_right*h_right,
604                        u_right*uh_right + 0.5*g*h_right**2,
605                        u_right*vh_right])
606
607    denom = s_max-s_min
608    if denom == 0.0:
609        edgeflux = array([0.0, 0.0, 0.0])
610        max_speed = 0.0
611    else:
612        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
613        edgeflux += s_max*s_min*(q_right-q_left)/denom
614
615        edgeflux = rotate(edgeflux, normal, direction=-1)
616        max_speed = max(abs(s_max), abs(s_min))
617
618    return edgeflux, max_speed
619
620def erfcc(x):
621
622    from math import fabs, exp
623
624    z=fabs(x)
625    t=1.0/(1.0+0.5*z)
626    result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
627         t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
628         t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
629    if x < 0.0:
630        result = 2.0-result
631
632    return result
633
634def flux_function_kinetic(normal, ql, qr, zl, zr):
635    """Compute fluxes between volumes for the shallow water wave equation
636    cast in terms of w = h+z using the 'central scheme' as described in
637
638    Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
639
640
641    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
642    in the numerical vectors ql an qr.
643
644    Bed elevations zl and zr.
645    """
646
647    from anuga.config import g, epsilon
648    from math import sqrt
649    from Numeric import array
650
651    #Align momentums with x-axis
652    q_left  = rotate(ql, normal, direction = 1)
653    q_right = rotate(qr, normal, direction = 1)
654
655    z = (zl+zr)/2 #Take average of field values
656
657    w_left  = q_left[0]   #w=h+z
658    h_left  = w_left-z
659    uh_left = q_left[1]
660
661    if h_left < epsilon:
662        u_left = 0.0  #Could have been negative
663        h_left = 0.0
664    else:
665        u_left  = uh_left/h_left
666
667
668    w_right  = q_right[0]  #w=h+z
669    h_right  = w_right-z
670    uh_right = q_right[1]
671
672
673    if h_right < epsilon:
674        u_right = 0.0  #Could have been negative
675        h_right = 0.0
676    else:
677        u_right  = uh_right/h_right
678
679    vh_left  = q_left[2]
680    vh_right = q_right[2]
681
682    soundspeed_left  = sqrt(g*h_left)
683    soundspeed_right = sqrt(g*h_right)
684
685    #Maximal wave speed
686    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
687
688    #Minimal wave speed
689    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
690
691    #Flux computation
692
693    F_left  = 0.0
694    F_right = 0.0
695    from math import sqrt, pi, exp
696    if h_left > 0.0:
697        F_left = u_left/sqrt(g*h_left)
698    if h_right > 0.0:
699        F_right = u_right/sqrt(g*h_right)
700
701    edgeflux = array([0.0, 0.0, 0.0])
702
703    edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \
704          h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
705          h_right*u_right/2.0*erfcc(F_right) -  \
706          h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
707
708    edgeflux[1] = (h_left*u_left**2 + g/2.0*h_left**2)/2.0*erfcc(-F_left) + \
709          u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
710          (h_right*u_right**2 + g/2.0*h_right**2)/2.0*erfcc(F_right) -  \
711          u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
712
713    edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \
714          vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \
715          vh_right*u_right/2.0*erfcc(F_right) - \
716          vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2))
717
718
719    edgeflux = rotate(edgeflux, normal, direction=-1)
720    max_speed = max(abs(s_max), abs(s_min))
721
722    return edgeflux, max_speed
723
724
725
726def compute_fluxes(domain):
727    """Compute all fluxes and the timestep suitable for all volumes
728    in domain.
729
730    Compute total flux for each conserved quantity using "flux_function"
731
732    Fluxes across each edge are scaled by edgelengths and summed up
733    Resulting flux is then scaled by area and stored in
734    explicit_update for each of the three conserved quantities
735    stage, xmomentum and ymomentum
736
737    The maximal allowable speed computed by the flux_function for each volume
738    is converted to a timestep that must not be exceeded. The minimum of
739    those is computed as the next overall timestep.
740
741    Post conditions:
742      domain.explicit_update is reset to computed flux values
743      domain.timestep is set to the largest step satisfying all volumes.
744    """
745
746    import sys
747
748    N = domain.number_of_elements
749
750    #Shortcuts
751    Stage = domain.quantities['stage']
752    Xmom = domain.quantities['xmomentum']
753    Ymom = domain.quantities['ymomentum']
754    Bed = domain.quantities['elevation']
755
756    #Arrays
757    stage = Stage.edge_values
758    xmom =  Xmom.edge_values
759    ymom =  Ymom.edge_values
760    bed =   Bed.edge_values
761
762    stage_bdry = Stage.boundary_values
763    xmom_bdry =  Xmom.boundary_values
764    ymom_bdry =  Ymom.boundary_values
765
766    flux = zeros(3, Float) #Work array for summing up fluxes
767
768
769    #Loop
770    timestep = float(sys.maxint)
771    for k in range(N):
772
773        flux[:] = 0.  #Reset work array
774        for i in range(3):
775            #Quantities inside volume facing neighbour i
776            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
777            zl = bed[k, i]
778
779            #Quantities at neighbour on nearest face
780            n = domain.neighbours[k,i]
781            if n < 0:
782                m = -n-1 #Convert negative flag to index
783                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
784                zr = zl #Extend bed elevation to boundary
785            else:
786                m = domain.neighbour_edges[k,i]
787                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
788                zr = bed[n, m]
789
790
791            #Outward pointing normal vector
792            normal = domain.normals[k, 2*i:2*i+2]
793
794            #Flux computation using provided function
795            edgeflux, max_speed = domain.flux_function(normal, ql, qr, zl, zr)
796            flux -= edgeflux * domain.edgelengths[k,i]
797
798            #Update optimal_timestep on full cells
799            if  domain.tri_full_flag[k] == 1:
800                try:
801                    timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
802                except ZeroDivisionError:
803                    pass
804
805        #Normalise by area and store for when all conserved
806        #quantities get updated
807        flux /= domain.areas[k]
808        Stage.explicit_update[k] = flux[0]
809        Xmom.explicit_update[k] = flux[1]
810        Ymom.explicit_update[k] = flux[2]
811
812
813    domain.timestep = timestep
814
815#MH090605 The following method belongs to the shallow_water domain class
816#see comments in the corresponding method in shallow_water_ext.c
817def extrapolate_second_order_sw_c(domain):
818    """Wrapper calling C version of extrapolate_second_order_sw
819    """
820    import sys
821
822    N = domain.number_of_elements
823
824    #Shortcuts
825    Stage = domain.quantities['stage']
826    Xmom = domain.quantities['xmomentum']
827    Ymom = domain.quantities['ymomentum']
828    Elevation = domain.quantities['elevation']
829    from shallow_water_ext import extrapolate_second_order_sw
830    extrapolate_second_order_sw(domain,
831                                domain.surrogate_neighbours,
832                                domain.number_of_boundaries,
833                                domain.centroid_coordinates,
834                                Stage.centroid_values,
835                                Xmom.centroid_values,
836                                Ymom.centroid_values,
837                                Elevation.centroid_values,
838                                domain.vertex_coordinates,
839                                Stage.vertex_values,
840                                Xmom.vertex_values,
841                                Ymom.vertex_values,
842                                Elevation.vertex_values)
843
844def compute_fluxes_c(domain):
845    """Wrapper calling C version of compute fluxes
846    """
847
848    import sys
849
850    N = domain.number_of_elements
851
852    #Shortcuts
853    Stage = domain.quantities['stage']
854    Xmom = domain.quantities['xmomentum']
855    Ymom = domain.quantities['ymomentum']
856    Bed = domain.quantities['elevation']
857
858    timestep = float(sys.maxint)
859    from shallow_water_ext import\
860         compute_fluxes_ext_central as compute_fluxes_ext
861   
862    domain.timestep = compute_fluxes_ext(timestep, domain.epsilon, domain.g,
863                                     domain.neighbours,
864                                     domain.neighbour_edges,
865                                     domain.normals,
866                                     domain.edgelengths,
867                                     domain.radii,
868                                     domain.areas,
869                                     domain.tri_full_flag,
870                                     Stage.edge_values,
871                                     Xmom.edge_values,
872                                     Ymom.edge_values,
873                                     Bed.edge_values,
874                                     Stage.boundary_values,
875                                     Xmom.boundary_values,
876                                     Ymom.boundary_values,
877                                     Stage.explicit_update,
878                                     Xmom.explicit_update,
879                                     Ymom.explicit_update,
880                                     domain.already_computed_flux)
881
882
883####################################
884# Module functions for gradient limiting (distribute_to_vertices_and_edges)
885
886def distribute_to_vertices_and_edges(domain):
887    """Distribution from centroids to vertices specific to the
888    shallow water wave
889    equation.
890
891    It will ensure that h (w-z) is always non-negative even in the
892    presence of steep bed-slopes by taking a weighted average between shallow
893    and deep cases.
894
895    In addition, all conserved quantities get distributed as per either a
896    constant (order==1) or a piecewise linear function (order==2).
897
898    FIXME: more explanation about removal of artificial variability etc
899
900    Precondition:
901      All quantities defined at centroids and bed elevation defined at
902      vertices.
903
904    Postcondition
905      Conserved quantities defined at vertices
906
907    """
908
909    from anuga.config import optimised_gradient_limiter
910
911    #Remove very thin layers of water
912    protect_against_infinitesimal_and_negative_heights(domain)
913
914    #Extrapolate all conserved quantities
915    if optimised_gradient_limiter:
916        #MH090605 if second order,
917        #perform the extrapolation and limiting on
918        #all of the conserved quantitie
919
920        if (domain._order_ == 1):
921            for name in domain.conserved_quantities:
922                Q = domain.quantities[name]
923                Q.extrapolate_first_order()
924        elif domain._order_ == 2:
925            domain.extrapolate_second_order_sw()
926        else:
927            raise 'Unknown order'
928    else:
929        #old code:
930        for name in domain.conserved_quantities:
931            Q = domain.quantities[name]
932            if domain._order_ == 1:
933                Q.extrapolate_first_order()
934            elif domain._order_ == 2:
935                Q.extrapolate_second_order()
936                Q.limit()
937            else:
938                raise 'Unknown order'
939
940
941    #Take bed elevation into account when water heights are small
942    balance_deep_and_shallow(domain)
943
944    #Compute edge values by interpolation
945    for name in domain.conserved_quantities:
946        Q = domain.quantities[name]
947        Q.interpolate_from_vertices_to_edges()
948
949
950def protect_against_infinitesimal_and_negative_heights(domain):
951    """Protect against infinitesimal heights and associated high velocities
952    """
953
954    #Shortcuts
955    wc = domain.quantities['stage'].centroid_values
956    zc = domain.quantities['elevation'].centroid_values
957    xmomc = domain.quantities['xmomentum'].centroid_values
958    ymomc = domain.quantities['ymomentum'].centroid_values
959    hc = wc - zc  #Water depths at centroids
960
961    #Update
962    #FIXME: Modify accroditg to c-version - or discard altogether.
963    for k in range(domain.number_of_elements):
964
965        if hc[k] < domain.minimum_allowed_height:
966            #Control stage
967            if hc[k] < domain.epsilon:
968                wc[k] = zc[k] # Contain 'lost mass' error
969
970            #Control momentum
971            xmomc[k] = ymomc[k] = 0.0
972
973
974def protect_against_infinitesimal_and_negative_heights_c(domain):
975    """Protect against infinitesimal heights and associated high velocities
976    """
977
978    #Shortcuts
979    wc = domain.quantities['stage'].centroid_values
980    zc = domain.quantities['elevation'].centroid_values
981    xmomc = domain.quantities['xmomentum'].centroid_values
982    ymomc = domain.quantities['ymomentum'].centroid_values
983
984    from shallow_water_ext import protect
985
986    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
987            domain.epsilon, wc, zc, xmomc, ymomc)
988
989
990
991def h_limiter(domain):
992    """Limit slopes for each volume to eliminate artificial variance
993    introduced by e.g. second order extrapolator
994
995    limit on h = w-z
996
997    This limiter depends on two quantities (w,z) so it resides within
998    this module rather than within quantity.py
999    """
1000
1001    N = domain.number_of_elements
1002    beta_h = domain.beta_h
1003
1004    #Shortcuts
1005    wc = domain.quantities['stage'].centroid_values
1006    zc = domain.quantities['elevation'].centroid_values
1007    hc = wc - zc
1008
1009    wv = domain.quantities['stage'].vertex_values
1010    zv = domain.quantities['elevation'].vertex_values
1011    hv = wv-zv
1012
1013    hvbar = zeros(hv.shape, Float) #h-limited values
1014
1015    #Find min and max of this and neighbour's centroid values
1016    hmax = zeros(hc.shape, Float)
1017    hmin = zeros(hc.shape, Float)
1018
1019    for k in range(N):
1020        hmax[k] = hmin[k] = hc[k]
1021        for i in range(3):
1022            n = domain.neighbours[k,i]
1023            if n >= 0:
1024                hn = hc[n] #Neighbour's centroid value
1025
1026                hmin[k] = min(hmin[k], hn)
1027                hmax[k] = max(hmax[k], hn)
1028
1029
1030    #Diffences between centroids and maxima/minima
1031    dhmax = hmax - hc
1032    dhmin = hmin - hc
1033
1034    #Deltas between vertex and centroid values
1035    dh = zeros(hv.shape, Float)
1036    for i in range(3):
1037        dh[:,i] = hv[:,i] - hc
1038
1039    #Phi limiter
1040    for k in range(N):
1041
1042        #Find the gradient limiter (phi) across vertices
1043        phi = 1.0
1044        for i in range(3):
1045            r = 1.0
1046            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
1047            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
1048
1049            phi = min( min(r*beta_h, 1), phi )
1050
1051        #Then update using phi limiter
1052        for i in range(3):
1053            hvbar[k,i] = hc[k] + phi*dh[k,i]
1054
1055    return hvbar
1056
1057
1058
1059def h_limiter_c(domain):
1060    """Limit slopes for each volume to eliminate artificial variance
1061    introduced by e.g. second order extrapolator
1062
1063    limit on h = w-z
1064
1065    This limiter depends on two quantities (w,z) so it resides within
1066    this module rather than within quantity.py
1067
1068    Wrapper for c-extension
1069    """
1070
1071    N = domain.number_of_elements
1072    beta_h = domain.beta_h
1073
1074    #Shortcuts
1075    wc = domain.quantities['stage'].centroid_values
1076    zc = domain.quantities['elevation'].centroid_values
1077    hc = wc - zc
1078
1079    wv = domain.quantities['stage'].vertex_values
1080    zv = domain.quantities['elevation'].vertex_values
1081    hv = wv - zv
1082
1083    #Call C-extension
1084    from shallow_water_ext import h_limiter_sw as h_limiter
1085    hvbar = h_limiter(domain, hc, hv)
1086
1087    return hvbar
1088
1089
1090def balance_deep_and_shallow(domain):
1091    """Compute linear combination between stage as computed by
1092    gradient-limiters limiting using w, and stage computed by
1093    gradient-limiters limiting using h (h-limiter).
1094    The former takes precedence when heights are large compared to the
1095    bed slope while the latter takes precedence when heights are
1096    relatively small.  Anything in between is computed as a balanced
1097    linear combination in order to avoid numerical disturbances which
1098    would otherwise appear as a result of hard switching between
1099    modes.
1100
1101    The h-limiter is always applied irrespective of the order.
1102    """
1103
1104    #Shortcuts
1105    wc = domain.quantities['stage'].centroid_values
1106    zc = domain.quantities['elevation'].centroid_values
1107    hc = wc - zc
1108
1109    wv = domain.quantities['stage'].vertex_values
1110    zv = domain.quantities['elevation'].vertex_values
1111    hv = wv-zv
1112
1113    #Limit h
1114    hvbar = h_limiter(domain)
1115
1116    for k in range(domain.number_of_elements):
1117        #Compute maximal variation in bed elevation
1118        #  This quantitiy is
1119        #    dz = max_i abs(z_i - z_c)
1120        #  and it is independent of dimension
1121        #  In the 1d case zc = (z0+z1)/2
1122        #  In the 2d case zc = (z0+z1+z2)/3
1123
1124        dz = max(abs(zv[k,0]-zc[k]),
1125                 abs(zv[k,1]-zc[k]),
1126                 abs(zv[k,2]-zc[k]))
1127
1128
1129        hmin = min( hv[k,:] )
1130
1131        #Create alpha in [0,1], where alpha==0 means using the h-limited
1132        #stage and alpha==1 means using the w-limited stage as
1133        #computed by the gradient limiter (both 1st or 2nd order)
1134
1135        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
1136        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
1137
1138        if dz > 0.0:
1139            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
1140        else:
1141            #Flat bed
1142            alpha = 1.0
1143
1144        #Let
1145        #
1146        #  wvi be the w-limited stage (wvi = zvi + hvi)
1147        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
1148        #
1149        #
1150        #where i=0,1,2 denotes the vertex ids
1151        #
1152        #Weighted balance between w-limited and h-limited stage is
1153        #
1154        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
1155        #
1156        #It follows that the updated wvi is
1157        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
1158        #
1159        # Momentum is balanced between constant and limited
1160
1161
1162        #for i in range(3):
1163        #    wv[k,i] = zv[k,i] + hvbar[k,i]
1164
1165        #return
1166
1167        if alpha < 1:
1168
1169            for i in range(3):
1170                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
1171
1172            #Momentums at centroids
1173            xmomc = domain.quantities['xmomentum'].centroid_values
1174            ymomc = domain.quantities['ymomentum'].centroid_values
1175
1176            #Momentums at vertices
1177            xmomv = domain.quantities['xmomentum'].vertex_values
1178            ymomv = domain.quantities['ymomentum'].vertex_values
1179
1180            # Update momentum as a linear combination of
1181            # xmomc and ymomc (shallow) and momentum
1182            # from extrapolator xmomv and ymomv (deep).
1183            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
1184            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
1185
1186
1187def balance_deep_and_shallow_c(domain):
1188    """Wrapper for C implementation
1189    """
1190
1191    #Shortcuts
1192    wc = domain.quantities['stage'].centroid_values
1193    zc = domain.quantities['elevation'].centroid_values
1194    hc = wc - zc
1195
1196    wv = domain.quantities['stage'].vertex_values
1197    zv = domain.quantities['elevation'].vertex_values
1198    hv = wv - zv
1199
1200    #Momentums at centroids
1201    xmomc = domain.quantities['xmomentum'].centroid_values
1202    ymomc = domain.quantities['ymomentum'].centroid_values
1203
1204    #Momentums at vertices
1205    xmomv = domain.quantities['xmomentum'].vertex_values
1206    ymomv = domain.quantities['ymomentum'].vertex_values
1207
1208    #Limit h
1209    hvbar = h_limiter(domain)
1210
1211    #This is how one would make a first order h_limited value
1212    #as in the old balancer (pre 17 Feb 2005):
1213    #from Numeric import zeros, Float
1214    #hvbar = zeros( (len(hc), 3), Float)
1215    #for i in range(3):
1216    #    hvbar[:,i] = hc[:]
1217
1218    from shallow_water_ext import balance_deep_and_shallow
[3876]1219    balance_deep_and_shallow(domain, wc, zc, hc, wv, zv, hv, hvbar,
[3804]1220                             xmomc, ymomc, xmomv, ymomv)
1221
1222
1223
1224
1225###############################################
1226#Boundaries - specific to the shallow water wave equation
1227class Reflective_boundary(Boundary):
1228    """Reflective boundary returns same conserved quantities as
1229    those present in its neighbour volume but reflected.
1230
1231    This class is specific to the shallow water equation as it
1232    works with the momentum quantities assumed to be the second
1233    and third conserved quantities.
1234    """
1235
1236    def __init__(self, domain = None):
1237        Boundary.__init__(self)
1238
1239        if domain is None:
1240            msg = 'Domain must be specified for reflective boundary'
1241            raise msg
1242
1243        #Handy shorthands
1244        self.stage   = domain.quantities['stage'].edge_values
1245        self.xmom    = domain.quantities['xmomentum'].edge_values
1246        self.ymom    = domain.quantities['ymomentum'].edge_values
1247        self.normals = domain.normals
1248
1249        self.conserved_quantities = zeros(3, Float)
1250
1251    def __repr__(self):
1252        return 'Reflective_boundary'
1253
1254
1255    def evaluate(self, vol_id, edge_id):
1256        """Reflective boundaries reverses the outward momentum
1257        of the volume they serve.
1258        """
1259
1260        q = self.conserved_quantities
1261        q[0] = self.stage[vol_id, edge_id]
1262        q[1] = self.xmom[vol_id, edge_id]
1263        q[2] = self.ymom[vol_id, edge_id]
1264
1265        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
1266
1267
1268        r = rotate(q, normal, direction = 1)
1269        r[1] = -r[1]
1270        q = rotate(r, normal, direction = -1)
1271
1272        return q
1273
1274
1275
1276class Transmissive_Momentum_Set_Stage_boundary(Boundary):
1277    """Returns same momentum conserved quantities as
1278    those present in its neighbour volume.
1279    Sets stage by specifying a function f of time which may either be a
1280    vector function or a scalar function
1281
1282    Example:
1283
1284    def waveform(t):
1285        return sea_level + normalized_amplitude/cosh(t-25)**2
1286
1287    Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
1288   
1289
1290    Underlying domain must be specified when boundary is instantiated
1291    """
1292
1293    def __init__(self, domain = None, function=None):
1294        Boundary.__init__(self)
1295
1296        if domain is None:
1297            msg = 'Domain must be specified for this type boundary'
1298            raise msg
1299
1300        if function is None:
1301            msg = 'Function must be specified for this type boundary'
1302            raise msg
1303
1304        self.domain   = domain
1305        self.function = function
1306
1307    def __repr__(self):
1308        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
1309
1310    def evaluate(self, vol_id, edge_id):
1311        """Transmissive Momentum Set Stage boundaries return the edge momentum
1312        values of the volume they serve.
1313        """
1314
1315        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
1316        value = self.function(self.domain.time)
1317
1318        try:
1319            x = float(value)
1320        except:   
1321            x = float(value[0])
1322           
1323        q[0] = x
1324        return q
1325
1326
1327        #FIXME: Consider this (taken from File_boundary) to allow
1328        #spatial variation
1329        #if vol_id is not None and edge_id is not None:
1330        #    i = self.boundary_indices[ vol_id, edge_id ]
1331        #    return self.F(t, point_id = i)
1332        #else:
1333        #    return self.F(t)
1334
1335
1336
1337class Dirichlet_Discharge_boundary(Boundary):
1338    """
1339    Sets stage (stage0)
1340    Sets momentum (wh0) in the inward normal direction.
1341
1342    Underlying domain must be specified when boundary is instantiated
1343    """
1344
1345    def __init__(self, domain = None, stage0=None, wh0=None):
1346        Boundary.__init__(self)
1347
1348        if domain is None:
1349            msg = 'Domain must be specified for this type boundary'
1350            raise msg
1351
1352        if stage0 is None:
1353            raise 'set stage'
1354
1355        if wh0 is None:
1356            wh0 = 0.0
1357
1358        self.domain   = domain
1359        self.stage0  = stage0
1360        self.wh0 = wh0
1361
1362    def __repr__(self):
1363        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
1364
1365    def evaluate(self, vol_id, edge_id):
1366        """Set discharge in the (inward) normal direction
1367        """
1368
1369        normal = self.domain.get_normal(vol_id,edge_id)
1370        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
1371        return q
1372
1373
1374        #FIXME: Consider this (taken from File_boundary) to allow
1375        #spatial variation
1376        #if vol_id is not None and edge_id is not None:
1377        #    i = self.boundary_indices[ vol_id, edge_id ]
1378        #    return self.F(t, point_id = i)
1379        #else:
1380        #    return self.F(t)
1381
1382
1383
1384
1385
1386#########################
1387#Standard forcing terms:
1388#
1389def gravity(domain):
1390    """Apply gravitational pull in the presence of bed slope
1391    """
1392
1393    xmom = domain.quantities['xmomentum'].explicit_update
1394    ymom = domain.quantities['ymomentum'].explicit_update
1395
1396    Stage = domain.quantities['stage']
1397    Elevation = domain.quantities['elevation']
1398    h = Stage.edge_values - Elevation.edge_values
1399    v = Elevation.vertex_values
1400
1401    x = domain.get_vertex_coordinates()
1402    g = domain.g
1403
1404    for k in range(domain.number_of_elements):
1405        avg_h = sum( h[k,:] )/3
1406
1407        #Compute bed slope
1408        x0, y0, x1, y1, x2, y2 = x[k,:]
1409        z0, z1, z2 = v[k,:]
1410
1411        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
1412
1413        #Update momentum
1414        xmom[k] += -g*zx*avg_h
1415        ymom[k] += -g*zy*avg_h
1416
1417
1418def gravity_c(domain):
1419    """Wrapper calling C version
1420    """
1421
1422    xmom = domain.quantities['xmomentum'].explicit_update
1423    ymom = domain.quantities['ymomentum'].explicit_update
1424
1425    Stage = domain.quantities['stage']
1426    Elevation = domain.quantities['elevation']
1427    h = Stage.edge_values - Elevation.edge_values
1428    v = Elevation.vertex_values
1429
1430    x = domain.get_vertex_coordinates()
1431    g = domain.g
1432
1433
1434    from shallow_water_ext import gravity
1435    gravity(g, h, v, x, xmom, ymom)
1436
1437
1438def manning_friction(domain):
1439    """Apply (Manning) friction to water momentum
1440    (Python version)
1441    """
1442
1443    from math import sqrt
1444
1445    w = domain.quantities['stage'].centroid_values
1446    z = domain.quantities['elevation'].centroid_values
1447    h = w-z
1448
1449    uh = domain.quantities['xmomentum'].centroid_values
1450    vh = domain.quantities['ymomentum'].centroid_values
1451    eta = domain.quantities['friction'].centroid_values
1452
1453    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1454    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1455
1456    N = domain.number_of_elements
1457    eps = domain.minimum_allowed_height
1458    g = domain.g
1459
1460    for k in range(N):
1461        if eta[k] >= eps:
1462            if h[k] >= eps:
1463                S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1464                S /= h[k]**(7.0/3)
1465
1466                #Update momentum
1467                xmom_update[k] += S*uh[k]
1468                ymom_update[k] += S*vh[k]
1469
1470
1471def manning_friction_implicit_c(domain):
1472    """Wrapper for c version
1473    """
1474
1475
1476    #print 'Implicit friction'
1477
1478    xmom = domain.quantities['xmomentum']
1479    ymom = domain.quantities['ymomentum']
1480
1481    w = domain.quantities['stage'].centroid_values
1482    z = domain.quantities['elevation'].centroid_values
1483
1484    uh = xmom.centroid_values
1485    vh = ymom.centroid_values
1486    eta = domain.quantities['friction'].centroid_values
1487
1488    xmom_update = xmom.semi_implicit_update
1489    ymom_update = ymom.semi_implicit_update
1490
1491    N = domain.number_of_elements
1492    eps = domain.minimum_allowed_height
1493    g = domain.g
1494
1495    from shallow_water_ext import manning_friction
1496    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1497
1498
1499def manning_friction_explicit_c(domain):
1500    """Wrapper for c version
1501    """
1502
1503    #print 'Explicit friction'
1504
1505    xmom = domain.quantities['xmomentum']
1506    ymom = domain.quantities['ymomentum']
1507
1508    w = domain.quantities['stage'].centroid_values
1509    z = domain.quantities['elevation'].centroid_values
1510
1511    uh = xmom.centroid_values
1512    vh = ymom.centroid_values
1513    eta = domain.quantities['friction'].centroid_values
1514
1515    xmom_update = xmom.explicit_update
1516    ymom_update = ymom.explicit_update
1517
1518    N = domain.number_of_elements
1519    eps = domain.minimum_allowed_height
1520    g = domain.g
1521
1522    from shallow_water_ext import manning_friction
1523    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1524
1525
1526def linear_friction(domain):
1527    """Apply linear friction to water momentum
1528
1529    Assumes quantity: 'linear_friction' to be present
1530    """
1531
1532    from math import sqrt
1533
1534    w = domain.quantities['stage'].centroid_values
1535    z = domain.quantities['elevation'].centroid_values
1536    h = w-z
1537
1538    uh = domain.quantities['xmomentum'].centroid_values
1539    vh = domain.quantities['ymomentum'].centroid_values
1540    tau = domain.quantities['linear_friction'].centroid_values
1541
1542    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1543    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1544
1545    N = domain.number_of_elements
1546    eps = domain.minimum_allowed_height
1547    g = domain.g #Not necessary? Why was this added?
1548
1549    for k in range(N):
1550        if tau[k] >= eps:
1551            if h[k] >= eps:
1552                S = -tau[k]/h[k]
1553
1554                #Update momentum
1555                xmom_update[k] += S*uh[k]
1556                ymom_update[k] += S*vh[k]
1557
1558
1559
1560def check_forcefield(f):
1561    """Check that f is either
1562    1: a callable object f(t,x,y), where x and y are vectors
1563       and that it returns an array or a list of same length
1564       as x and y
1565    2: a scalar
1566    """
1567
1568    if callable(f):
1569        N = 3
1570        x = ones(3, Float)
1571        y = ones(3, Float)
1572        try:
1573            q = f(1.0, x=x, y=y)
1574        except Exception, e:
1575            msg = 'Function %s could not be executed:\n%s' %(f, e)
1576            #FIXME: Reconsider this semantics
1577            raise msg
1578
1579        try:
1580            q = array(q).astype(Float)
1581        except:
1582            msg = 'Return value from vector function %s could ' %f
1583            msg += 'not be converted into a Numeric array of floats.\n'
1584            msg += 'Specified function should return either list or array.'
1585            raise msg
1586
1587        #Is this really what we want?
1588        msg = 'Return vector from function %s ' %f
1589        msg += 'must have same lenght as input vectors'
1590        assert len(q) == N, msg
1591
1592    else:
1593        try:
1594            f = float(f)
1595        except:
1596            msg = 'Force field %s must be either a scalar' %f
1597            msg += ' or a vector function'
1598            raise Exception(msg)
1599    return f
1600
1601
1602class Wind_stress:
1603    """Apply wind stress to water momentum in terms of
1604    wind speed [m/s] and wind direction [degrees]
1605    """
1606
1607    def __init__(self, *args, **kwargs):
1608        """Initialise windfield from wind speed s [m/s]
1609        and wind direction phi [degrees]
1610
1611        Inputs v and phi can be either scalars or Python functions, e.g.
1612
1613        W = Wind_stress(10, 178)
1614
1615        #FIXME - 'normal' degrees are assumed for now, i.e. the
1616        vector (1,0) has zero degrees.
1617        We may need to convert from 'compass' degrees later on and also
1618        map from True north to grid north.
1619
1620        Arguments can also be Python functions of t,x,y as in
1621
1622        def speed(t,x,y):
1623            ...
1624            return s
1625
1626        def angle(t,x,y):
1627            ...
1628            return phi
1629
1630        where x and y are vectors.
1631
1632        and then pass the functions in
1633
1634        W = Wind_stress(speed, angle)
1635
1636        The instantiated object W can be appended to the list of
1637        forcing_terms as in
1638
1639        Alternatively, one vector valued function for (speed, angle)
1640        can be applied, providing both quantities simultaneously.
1641        As in
1642        W = Wind_stress(F), where returns (speed, angle) for each t.
1643
1644        domain.forcing_terms.append(W)
1645        """
1646
1647        from anuga.config import rho_a, rho_w, eta_w
1648        from Numeric import array, Float
1649
1650        if len(args) == 2:
1651            s = args[0]
1652            phi = args[1]
1653        elif len(args) == 1:
1654            #Assume vector function returning (s, phi)(t,x,y)
1655            vector_function = args[0]
1656            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1657            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1658        else:
1659           #Assume info is in 2 keyword arguments
1660
1661           if len(kwargs) == 2:
1662               s = kwargs['s']
1663               phi = kwargs['phi']
1664           else:
1665               raise 'Assumes two keyword arguments: s=..., phi=....'
1666
1667        self.speed = check_forcefield(s)
1668        self.phi = check_forcefield(phi)
1669
1670        self.const = eta_w*rho_a/rho_w
1671
1672
1673    def __call__(self, domain):
1674        """Evaluate windfield based on values found in domain
1675        """
1676
1677        from math import pi, cos, sin, sqrt
1678        from Numeric import Float, ones, ArrayType
1679
1680        xmom_update = domain.quantities['xmomentum'].explicit_update
1681        ymom_update = domain.quantities['ymomentum'].explicit_update
1682
1683        N = domain.number_of_elements
1684        t = domain.time
1685
1686        if callable(self.speed):
1687            xc = domain.get_centroid_coordinates()
1688            s_vec = self.speed(t, xc[:,0], xc[:,1])
1689        else:
1690            #Assume s is a scalar
1691
1692            try:
1693                s_vec = self.speed * ones(N, Float)
1694            except:
1695                msg = 'Speed must be either callable or a scalar: %s' %self.s
1696                raise msg
1697
1698
1699        if callable(self.phi):
1700            xc = domain.get_centroid_coordinates()
1701            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1702        else:
1703            #Assume phi is a scalar
1704
1705            try:
1706                phi_vec = self.phi * ones(N, Float)
1707            except:
1708                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1709                raise msg
1710
1711        assign_windfield_values(xmom_update, ymom_update,
1712                                s_vec, phi_vec, self.const)
1713
1714
1715def assign_windfield_values(xmom_update, ymom_update,
1716                            s_vec, phi_vec, const):
1717    """Python version of assigning wind field to update vectors.
1718    A c version also exists (for speed)
1719    """
1720    from math import pi, cos, sin, sqrt
1721
1722    N = len(s_vec)
1723    for k in range(N):
1724        s = s_vec[k]
1725        phi = phi_vec[k]
1726
1727        #Convert to radians
1728        phi = phi*pi/180
1729
1730        #Compute velocity vector (u, v)
1731        u = s*cos(phi)
1732        v = s*sin(phi)
1733
1734        #Compute wind stress
1735        S = const * sqrt(u**2 + v**2)
1736        xmom_update[k] += S*u
1737        ymom_update[k] += S*v
1738
1739
1740
1741##############################
1742#OBSOLETE STUFF
1743
1744def balance_deep_and_shallow_old(domain):
1745    """Compute linear combination between stage as computed by
1746    gradient-limiters and stage computed as constant height above bed.
1747    The former takes precedence when heights are large compared to the
1748    bed slope while the latter takes precedence when heights are
1749    relatively small.  Anything in between is computed as a balanced
1750    linear combination in order to avoid numerical disturbances which
1751    would otherwise appear as a result of hard switching between
1752    modes.
1753    """
1754
1755    #OBSOLETE
1756
1757    #Shortcuts
1758    wc = domain.quantities['stage'].centroid_values
1759    zc = domain.quantities['elevation'].centroid_values
1760    hc = wc - zc
1761
1762    wv = domain.quantities['stage'].vertex_values
1763    zv = domain.quantities['elevation'].vertex_values
1764    hv = wv-zv
1765
1766
1767    #Computed linear combination between constant stages and and
1768    #stages parallel to the bed elevation.
1769    for k in range(domain.number_of_elements):
1770        #Compute maximal variation in bed elevation
1771        #  This quantitiy is
1772        #    dz = max_i abs(z_i - z_c)
1773        #  and it is independent of dimension
1774        #  In the 1d case zc = (z0+z1)/2
1775        #  In the 2d case zc = (z0+z1+z2)/3
1776
1777        dz = max(abs(zv[k,0]-zc[k]),
1778                 abs(zv[k,1]-zc[k]),
1779                 abs(zv[k,2]-zc[k]))
1780
1781
1782        hmin = min( hv[k,:] )
1783
1784        #Create alpha in [0,1], where alpha==0 means using shallow
1785        #first order scheme and alpha==1 means using the stage w as
1786        #computed by the gradient limiter (1st or 2nd order)
1787        #
1788        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
1789        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
1790
1791        if dz > 0.0:
1792            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
1793        else:
1794            #Flat bed
1795            alpha = 1.0
1796
1797        #Weighted balance between stage parallel to bed elevation
1798        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
1799        #order gradient limiter
1800        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
1801        #
1802        #It follows that the updated wvi is
1803        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
1804        #  zvi + hc + alpha*(hvi - hc)
1805        #
1806        #Note that hvi = zc+hc-zvi in the first order case (constant).
1807
1808        if alpha < 1:
1809            for i in range(3):
1810                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
1811
1812
1813            #Momentums at centroids
1814            xmomc = domain.quantities['xmomentum'].centroid_values
1815            ymomc = domain.quantities['ymomentum'].centroid_values
1816
1817            #Momentums at vertices
1818            xmomv = domain.quantities['xmomentum'].vertex_values
1819            ymomv = domain.quantities['ymomentum'].vertex_values
1820
1821            # Update momentum as a linear combination of
1822            # xmomc and ymomc (shallow) and momentum
1823            # from extrapolator xmomv and ymomv (deep).
1824            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
1825            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
1826
1827
1828
1829
1830
1831###########################
1832###########################
1833#Geometries
1834
1835
1836#FIXME: Rethink this way of creating values.
1837
1838
1839class Weir:
1840    """Set a bathymetry for weir with a hole and a downstream gutter
1841    x,y are assumed to be in the unit square
1842    """
1843
1844    def __init__(self, stage):
1845        self.inflow_stage = stage
1846
1847    def __call__(self, x, y):
1848        from Numeric import zeros, Float
1849        from math import sqrt
1850
1851        N = len(x)
1852        assert N == len(y)
1853
1854        z = zeros(N, Float)
1855        for i in range(N):
1856            z[i] = -x[i]/2  #General slope
1857
1858            #Flattish bit to the left
1859            if x[i] < 0.3:
1860                z[i] = -x[i]/10
1861
1862            #Weir
1863            if x[i] >= 0.3 and x[i] < 0.4:
1864                z[i] = -x[i]+0.9
1865
1866            #Dip
1867            x0 = 0.6
1868            #depth = -1.3
1869            depth = -1.0
1870            #plateaux = -0.9
1871            plateaux = -0.6
1872            if y[i] < 0.7:
1873                if x[i] > x0 and x[i] < 0.9:
1874                    z[i] = depth
1875
1876                #RHS plateaux
1877                if x[i] >= 0.9:
1878                    z[i] = plateaux
1879
1880
1881            elif y[i] >= 0.7 and y[i] < 1.5:
1882                #Restrict and deepen
1883                if x[i] >= x0 and x[i] < 0.8:
1884                    z[i] = depth-(y[i]/3-0.3)
1885                    #z[i] = depth-y[i]/5
1886                    #z[i] = depth
1887                elif x[i] >= 0.8:
1888                    #RHS plateaux
1889                    z[i] = plateaux
1890
1891            elif y[i] >= 1.5:
1892                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
1893                    #Widen up and stay at constant depth
1894                    z[i] = depth-1.5/5
1895                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
1896                    #RHS plateaux
1897                    z[i] = plateaux
1898
1899
1900            #Hole in weir (slightly higher than inflow condition)
1901            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
1902                z[i] = -x[i]+self.inflow_stage + 0.02
1903
1904            #Channel behind weir
1905            x0 = 0.5
1906            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
1907                z[i] = -x[i]+self.inflow_stage + 0.02
1908
1909            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
1910                #Flatten it out towards the end
1911                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
1912
1913            #Hole to the east
1914            x0 = 1.1; y0 = 0.35
1915            #if x[i] < -0.2 and y < 0.5:
1916            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
1917                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
1918
1919            #Tiny channel draining hole
1920            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
1921                z[i] = -0.9 #North south
1922
1923            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
1924                z[i] = -1.0 + (x[i]-0.9)/3 #East west
1925
1926
1927
1928            #Stuff not in use
1929
1930            #Upward slope at inlet to the north west
1931            #if x[i] < 0.0: # and y[i] > 0.5:
1932            #    #z[i] = -y[i]+0.5  #-x[i]/2
1933            #    z[i] = x[i]/4 - y[i]**2 + 0.5
1934
1935            #Hole to the west
1936            #x0 = -0.4; y0 = 0.35 # center
1937            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
1938            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
1939
1940
1941
1942
1943
1944        return z/2
1945
1946class Weir_simple:
1947    """Set a bathymetry for weir with a hole and a downstream gutter
1948    x,y are assumed to be in the unit square
1949    """
1950
1951    def __init__(self, stage):
1952        self.inflow_stage = stage
1953
1954    def __call__(self, x, y):
1955        from Numeric import zeros, Float
1956
1957        N = len(x)
1958        assert N == len(y)
1959
1960        z = zeros(N, Float)
1961        for i in range(N):
1962            z[i] = -x[i]  #General slope
1963
1964            #Flat bit to the left
1965            if x[i] < 0.3:
1966                z[i] = -x[i]/10  #General slope
1967
1968            #Weir
1969            if x[i] > 0.3 and x[i] < 0.4:
1970                z[i] = -x[i]+0.9
1971
1972            #Dip
1973            if x[i] > 0.6 and x[i] < 0.9:
1974                z[i] = -x[i]-0.5  #-y[i]/5
1975
1976            #Hole in weir (slightly higher than inflow condition)
1977            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
1978                z[i] = -x[i]+self.inflow_stage + 0.05
1979
1980
1981        return z/2
1982
1983
1984
1985class Constant_stage:
1986    """Set an initial condition with constant stage
1987    """
1988    def __init__(self, s):
1989        self.s = s
1990
1991    def __call__(self, x, y):
1992        return self.s
1993
1994class Constant_height:
1995    """Set an initial condition with constant water height, e.g
1996    stage s = z+h
1997    """
1998
1999    def __init__(self, W, h):
2000        self.W = W
2001        self.h = h
2002
2003    def __call__(self, x, y):
2004        if self.W is None:
2005            from Numeric import ones, Float
2006            return self.h*ones(len(x), Float)
2007        else:
2008            return self.W(x,y) + self.h
2009
2010
2011
2012
2013def compute_fluxes_python(domain):
2014    """Compute all fluxes and the timestep suitable for all volumes
2015    in domain.
2016
2017    Compute total flux for each conserved quantity using "flux_function"
2018
2019    Fluxes across each edge are scaled by edgelengths and summed up
2020    Resulting flux is then scaled by area and stored in
2021    explicit_update for each of the three conserved quantities
2022    stage, xmomentum and ymomentum
2023
2024    The maximal allowable speed computed by the flux_function for each volume
2025    is converted to a timestep that must not be exceeded. The minimum of
2026    those is computed as the next overall timestep.
2027
2028    Post conditions:
2029      domain.explicit_update is reset to computed flux values
2030      domain.timestep is set to the largest step satisfying all volumes.
2031    """
2032
2033    import sys
2034    from Numeric import zeros, Float
2035
2036    N = domain.number_of_elements
2037
2038    #Shortcuts
2039    Stage = domain.quantities['stage']
2040    Xmom = domain.quantities['xmomentum']
2041    Ymom = domain.quantities['ymomentum']
2042    Bed = domain.quantities['elevation']
2043
2044    #Arrays
2045    stage = Stage.edge_values
2046    xmom =  Xmom.edge_values
2047    ymom =  Ymom.edge_values
2048    bed =   Bed.edge_values
2049
2050    stage_bdry = Stage.boundary_values
2051    xmom_bdry =  Xmom.boundary_values
2052    ymom_bdry =  Ymom.boundary_values
2053
2054    flux = zeros((N,3), Float) #Work array for summing up fluxes
2055
2056    #Loop
2057    timestep = float(sys.maxint)
2058    for k in range(N):
2059
2060        for i in range(3):
2061            #Quantities inside volume facing neighbour i
2062            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
2063            zl = bed[k, i]
2064
2065            #Quantities at neighbour on nearest face
2066            n = domain.neighbours[k,i]
2067            if n < 0:
2068                m = -n-1 #Convert negative flag to index
2069                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
2070                zr = zl #Extend bed elevation to boundary
2071            else:
2072                m = domain.neighbour_edges[k,i]
2073                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
2074                zr = bed[n, m]
2075
2076
2077            #Outward pointing normal vector
2078            normal = domain.normals[k, 2*i:2*i+2]
2079
2080            #Flux computation using provided function
2081            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
2082
2083            flux[k,:] = edgeflux
2084
2085    return flux
2086
2087
2088
2089
2090
2091
2092
2093##############################################
2094#Initialise module
2095
2096
2097from anuga.utilities import compile
2098if compile.can_use_C_extension('shallow_water_ext.c'):
2099    #Replace python version with c implementations
2100
2101    from shallow_water_ext import rotate, assign_windfield_values
2102    compute_fluxes = compute_fluxes_c
2103    extrapolate_second_order_sw=extrapolate_second_order_sw_c
2104    gravity = gravity_c
2105    manning_friction = manning_friction_implicit_c
2106    h_limiter = h_limiter_c
2107    balance_deep_and_shallow = balance_deep_and_shallow_c
2108    protect_against_infinitesimal_and_negative_heights =\
2109                    protect_against_infinitesimal_and_negative_heights_c
2110
2111
2112    #distribute_to_vertices_and_edges =\
2113    #              distribute_to_vertices_and_edges_c #(like MH's)
2114
2115
2116
2117#Optimisation with psyco
2118from anuga.config import use_psyco
2119if use_psyco:
2120    try:
2121        import psyco
2122    except:
2123        import os
2124        if os.name == 'posix' and os.uname()[4] == 'x86_64':
2125            pass
2126            #Psyco isn't supported on 64 bit systems, but it doesn't matter
2127        else:
2128            msg = 'WARNING: psyco (speedup) could not import'+\
2129                  ', you may want to consider installing it'
2130            print msg
2131    else:
2132        psyco.bind(Domain.distribute_to_vertices_and_edges)
2133        psyco.bind(Domain.compute_fluxes)
2134
2135if __name__ == "__main__":
2136    pass
2137
2138# Profiling stuff
2139import profile
2140profiler = profile.Profile()
Note: See TracBrowser for help on using the repository browser.