source: anuga_core/source/obsolete_code/shallow_water_domain_kinetic.py @ 6011

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

Completed committing changeset:4707 and changeset:4708

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