source: development/pyvolution-1d/shallow_water_1d.py @ 3322

Last change on this file since 3322 was 3293, checked in by jakeman, 19 years ago

Updating lots of changes
File size: 40.4 KB
Line 
1"""Class Domain -
21D interval 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 = S
10
11where
12
13U = [w, uh]
14E = [uh, u^2h + gh^2/2]
15S represents source terms forcing the system
16(e.g. gravity, friction, wind stress, ...)
17
18and _t, _x, _y denote the derivative with respect to t, x and y respectively.
19
20The quantities are
21
22symbol    variable name    explanation
23x         x                horizontal distance from origin [m]
24z         elevation        elevation of bed on which flow is modelled [m]
25h         height           water height above z [m]
26w         stage            absolute water level, w = z+h [m]
27u                          speed in the x direction [m/s]
28uh        xmomentum        momentum in the x direction [m^2/s]
29
30eta                        mannings friction coefficient [to appear]
31nu                         wind stress coefficient [to appear]
32
33The conserved quantities are w, uh
34
35For details see e.g.
36Christopher Zoppou and Stephen Roberts,
37Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
38Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
39
40
41John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
42Geoscience Australia, 2006
43"""
44
45from domain import *
46Generic_Domain = Domain #Rename
47
48#Shallow water domain
49class Domain(Generic_Domain):
50
51    def __init__(self, coordinates, boundary = None, tagged_elements = None,
52                 geo_reference = None):
53
54        conserved_quantities = ['stage', 'xmomentum']
55        other_quantities = ['elevation', 'friction']
56        Generic_Domain.__init__(self, coordinates, boundary,
57                                conserved_quantities, other_quantities,
58                                tagged_elements, geo_reference)
59       
60        from config import minimum_allowed_height, g
61        self.minimum_allowed_height = minimum_allowed_height
62        self.g = g
63
64        #forcing terms not included in 1d domain ?WHy?
65        #self.forcing_terms.append(gravity)
66        #self.forcing_terms.append(manning_friction)
67        #print "\nI have Removed forcing terms line 64 1dsw"
68
69        #Realtime visualisation
70        self.visualiser = None
71        self.visualise  = False
72        self.visualise_color_stage = False
73        self.visualise_stage_range = 1.0
74        self.visualise_timer = True
75        self.visualise_range_z = None
76       
77        #Stored output
78        self.store = True
79        self.format = 'sww'
80        self.smooth = True
81
82        #Reduction operation for get_vertex_values
83        from util import mean
84        self.reduction = mean
85        #self.reduction = min  #Looks better near steep slopes
86
87        self.quantities_to_be_stored = ['stage','xmomentum']
88
89    def set_quantities_to_be_stored(self, q):
90        """Specify which quantities will be stored in the sww file.
91
92        q must be either:
93          - the name of a quantity
94          - a list of quantity names
95          - None
96
97        In the two first cases, the named quantities will be stored at each
98        yieldstep
99        (This is in addition to the quantities elevation and friction) 
100        If q is None, storage will be switched off altogether.
101        """
102
103
104        if q is None:
105            self.quantities_to_be_stored = []           
106            self.store = False
107            return
108
109        if isinstance(q, basestring):
110            q = [q] # Turn argument into a list
111
112        #Check correcness   
113        for quantity_name in q:
114            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
115            assert quantity_name in self.conserved_quantities, msg
116       
117        self.quantities_to_be_stored = q
118       
119
120    def initialise_visualiser(self,scale_z=1.0,rect=None):
121        #Realtime visualisation
122        if self.visualiser is None:
123            from realtime_visualisation_new import Visualiser
124            self.visualiser = Visualiser(self,scale_z,rect)
125            self.visualiser.setup['elevation']=True
126            self.visualiser.updating['stage']=True
127        self.visualise = True
128        if self.visualise_color_stage == True:
129            self.visualiser.coloring['stage'] = True
130            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
131        print 'initialise visualiser'
132        print self.visualiser.setup
133        print self.visualiser.updating
134
135    def check_integrity(self):
136        Generic_Domain.check_integrity(self)
137        #Check that we are solving the shallow water wave equation
138
139        msg = 'First conserved quantity must be "stage"'
140        assert self.conserved_quantities[0] == 'stage', msg
141        msg = 'Second conserved quantity must be "xmomentum"'
142        assert self.conserved_quantities[1] == 'xmomentum', msg
143
144    def extrapolate_second_order_sw(self):
145        #Call correct module function
146        #(either from this module or C-extension)
147        extrapolate_second_order_sw(self)
148
149    def compute_fluxes(self):
150        #Call correct module function
151        #(either from this module or C-extension)
152        compute_fluxes(self)
153
154    def distribute_to_vertices_and_edges(self):
155        #Call correct module function
156        #(either from this module or C-extension)
157        distribute_to_vertices_and_edges(self)
158
159    def evolve(self, yieldstep = None, finaltime = None,
160               skip_initial_step = False):
161        """Specialisation of basic evolve method from parent class
162        """
163
164        #Call check integrity here rather than from user scripts
165        #self.check_integrity()
166
167        msg = 'Parameter beta_h must be in the interval [0, 1)'
168        assert 0 <= self.beta_h < 1.0, msg
169        msg = 'Parameter beta_w must be in the interval [0, 1)'
170        assert 0 <= self.beta_w < 1.0, msg
171
172
173        #Initial update of vertex and edge values before any storage
174        #and or visualisation
175        self.distribute_to_vertices_and_edges()
176
177        #Initialise real time viz if requested
178        #if self.visualise is True and self.time == 0.0:
179        #    if self.visualiser is None:
180        #        self.initialise_visualiser()
181        #
182        #    self.visualiser.update_timer()
183        #    self.visualiser.setup_all()
184
185        #Store model data, e.g. for visualisation
186        #if self.store is True and self.time == 0.0:
187        #    self.initialise_storage()
188        #    #print 'Storing results in ' + self.writer.filename
189        #else:
190        #    pass
191        #    #print 'Results will not be stored.'
192        #    #print 'To store results set domain.store = True'
193        #    #FIXME: Diagnostic output should be controlled by
194        #    # a 'verbose' flag living in domain (or in a parent class)
195
196        #Call basic machinery from parent class
197        for t in Generic_Domain.evolve(self, yieldstep, finaltime,
198                                       skip_initial_step):
199            #Real time viz
200        #    if self.visualise is True:
201        #        self.visualiser.update_all()
202        #        self.visualiser.update_timer()
203
204
205            #Store model data, e.g. for subsequent visualisation
206        #    if self.store is True:
207        #        self.store_timestep(self.quantities_to_be_stored)
208
209            #FIXME: Could maybe be taken from specified list
210            #of 'store every step' quantities
211
212            #Pass control on to outer loop for more specific actions
213            yield(t)
214
215    def initialise_storage(self):
216        """Create and initialise self.writer object for storing data.
217        Also, save x and bed elevation
218        """
219
220        import data_manager
221
222        #Initialise writer
223        self.writer = data_manager.get_dataobject(self, mode = 'w')
224
225        #Store vertices and connectivity
226        self.writer.store_connectivity()
227
228
229    def store_timestep(self, name):
230        """Store named quantity and time.
231
232        Precondition:
233           self.write has been initialised
234        """
235        self.writer.store_timestep(name)
236
237
238#=============== End of Shallow Water Domain ===============================
239
240#Rotation of momentum vector
241def rotate(q, normal, direction = 1):
242    """Rotate the momentum component q (q[1], q[2])
243    from x,y coordinates to coordinates based on normal vector.
244
245    If direction is negative the rotation is inverted.
246
247    Input vector is preserved
248
249    This function is specific to the shallow water wave equation
250    """
251
252    from Numeric import zeros, Float
253
254    assert len(q) == 3,\
255           'Vector of conserved quantities must have length 3'\
256           'for 2D shallow water equation'
257
258    try:
259        l = len(normal)
260    except:
261        raise 'Normal vector must be an Numeric array'
262
263    assert l == 2, 'Normal vector must have 2 components'
264
265
266    n1 = normal[0]
267    n2 = normal[1]
268
269    r = zeros(len(q), Float) #Rotated quantities
270    r[0] = q[0]              #First quantity, height, is not rotated
271
272    if direction == -1:
273        n2 = -n2
274
275
276    r[1] =  n1*q[1] + n2*q[2]
277    r[2] = -n2*q[1] + n1*q[2]
278
279    return r
280
281
282# Flux computation
283#def flux_function(normal, ql, qr, zl, zr):
284def flux_function(normal, ql, qr, zl, zr):
285    """Compute fluxes between volumes for the shallow water wave equation
286    cast in terms of w = h+z using the 'central scheme' as described in
287
288    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
289    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
290    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
291
292    The implemented formula is given in equation (3.15) on page 714
293
294    Conserved quantities w, uh, are stored as elements 0 and 1
295    in the numerical vectors ql an qr.
296
297    Bed elevations zl and zr.
298    """
299
300    from config import g, epsilon
301    from math import sqrt
302    from Numeric import array
303
304    #print 'ql',ql
305
306    #Align momentums with x-axis
307    #q_left  = rotate(ql, normal, direction = 1)
308    #q_right = rotate(qr, normal, direction = 1)
309    q_left = ql
310    q_left[1] = q_left[1]*normal
311    q_right = qr
312    q_right[1] = q_right[1]*normal
313
314    #z = (zl+zr)/2 #Take average of field values
315    z = 0.5*(zl+zr) #Take average of field values
316
317    w_left  = q_left[0]   #w=h+z
318    h_left  = w_left-z
319    uh_left = q_left[1]
320
321    if h_left < epsilon:
322        u_left = 0.0  #Could have been negative
323        h_left = 0.0
324    else:
325        u_left  = uh_left/h_left
326
327
328    w_right  = q_right[0]  #w=h+z
329    h_right  = w_right-z
330    uh_right = q_right[1]
331
332
333    if h_right < epsilon:
334        u_right = 0.0  #Could have been negative
335        h_right = 0.0
336    else:
337        u_right  = uh_right/h_right
338
339    #vh_left  = q_left[2]
340    #vh_right = q_right[2]
341
342    soundspeed_left  = sqrt(g*h_left)
343    soundspeed_right = sqrt(g*h_right)
344
345    #Maximal wave speed
346    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
347   
348    #Minimal wave speed
349    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
350   
351    #Flux computation
352
353    #flux_left  = array([u_left*h_left,
354    #                    u_left*uh_left + 0.5*g*h_left**2])
355    #flux_right = array([u_right*h_right,
356    #                    u_right*uh_right + 0.5*g*h_right**2])
357    flux_left  = array([u_left*h_left,
358                        u_left*uh_left + 0.5*g*h_left*h_left])
359    flux_right = array([u_right*h_right,
360                        u_right*uh_right + 0.5*g*h_right*h_right])
361
362    denom = s_max-s_min
363    if denom == 0.0:
364        edgeflux = array([0.0, 0.0])
365        max_speed = 0.0
366    else:
367        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
368        edgeflux += s_max*s_min*(q_right-q_left)/denom
369       
370        edgeflux[1] = edgeflux[1]*normal
371
372        max_speed = max(abs(s_max), abs(s_min))
373
374    return edgeflux, max_speed
375
376def compute_fluxes(domain):
377    """Compute all fluxes and the timestep suitable for all volumes
378    in domain.
379
380    Compute total flux for each conserved quantity using "flux_function"
381
382    Fluxes across each edge are scaled by edgelengths and summed up
383    Resulting flux is then scaled by area and stored in
384    explicit_update for each of the three conserved quantities
385    stage, xmomentum and ymomentum
386
387    The maximal allowable speed computed by the flux_function for each volume
388    is converted to a timestep that must not be exceeded. The minimum of
389    those is computed as the next overall timestep.
390
391    Post conditions:
392      domain.explicit_update is reset to computed flux values
393      domain.timestep is set to the largest step satisfying all volumes.
394    """
395
396    import sys
397    from Numeric import zeros, Float
398
399    N = domain.number_of_elements
400
401    #Shortcuts
402    Stage = domain.quantities['stage']
403    Xmom = domain.quantities['xmomentum']
404#    Ymom = domain.quantities['ymomentum']
405    Bed = domain.quantities['elevation']
406
407    #Arrays
408    #stage = Stage.edge_values
409    #xmom =  Xmom.edge_values
410 #   ymom =  Ymom.edge_values
411    #bed =   Bed.edge_values
412   
413    stage = Stage.vertex_values
414    xmom =  Xmom.vertex_values
415    bed =   Bed.vertex_values
416   
417    #print 'stage edge values', stage
418    #print 'xmom edge values', xmom
419    #print 'bed values', bed
420
421    stage_bdry = Stage.boundary_values
422    xmom_bdry =  Xmom.boundary_values
423    #print 'stage_bdry',stage_bdry
424    #print 'xmom_bdry', xmom_bdry
425#    ymom_bdry =  Ymom.boundary_values
426
427#    flux = zeros(3, Float) #Work array for summing up fluxes
428    flux = zeros(2, Float) #Work array for summing up fluxes
429    ql = zeros(2, Float)
430    qr = zeros(2, Float)
431
432    #Loop
433    timestep = float(sys.maxint)
434    for k in range(N):
435
436        flux[:] = 0.  #Reset work array
437        #for i in range(3):
438        for i in range(2):
439            #Quantities inside volume facing neighbour i
440            #ql[0] = stage[k, i]
441            #ql[1] = xmom[k, i]
442            ql = [stage[k, i], xmom[k, i]]
443            zl = bed[k, i]
444
445            #Quantities at neighbour on nearest face
446            n = domain.neighbours[k,i]
447            if n < 0:
448                m = -n-1 #Convert negative flag to index
449                qr[0] = stage_bdry[m]
450                qr[1] = xmom_bdry[m]
451                zr = zl #Extend bed elevation to boundary
452            else:
453                #m = domain.neighbour_edges[k,i]
454                m = domain.neighbour_vertices[k,i]
455                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
456                qr[0] = stage[n, m]
457                qr[1] =  xmom[n, m]
458                zr = bed[n, m]
459
460
461            #Outward pointing normal vector
462            normal = domain.normals[k, i]
463       
464            #Flux computation using provided function
465            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
466            #print 'ql',ql
467            #print 'qr',qr
468           
469            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
470            #print 'edgeflux', edgeflux
471
472            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
473            # flux = edgefluxleft - edgefluxright
474            flux -= edgeflux #* domain.edgelengths[k,i]
475
476            #Update optimal_timestep
477            try:
478                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
479                #timestep = 0.01
480                timestep = min(timestep, 0.5*domain.areas[k]/max_speed)
481                if timestep < 0.00001:
482                    #print 'max_speed', max_speed
483                    s = domain.quantities['stage']
484                    s = s.centroid_values
485                    xm = domain.quantities['xmomentum']
486                    xm = xm.centroid_values
487                    #print 'h', s[k]
488                    #print 'xm', xm[k]
489                    #print 'u', xm[k]/s[k]
490                    #break
491                #timestep = 0.01
492                #print 'areas', domain.areas[k]
493                #print "timestep", timestep
494            except ZeroDivisionError:
495                pass
496
497        #Normalise by area and store for when all conserved
498        #quantities get updated
499        #flux /= domain.areas[k]
500        # ADD ABOVE LINE AGAIN
501        Stage.explicit_update[k] = flux[0]
502        Xmom.explicit_update[k] = flux[1]
503        #Ymom.explicit_update[k] = flux[2]
504
505
506    domain.timestep = timestep
507
508#see comments in the corresponding method in shallow_water_ext.c
509def extrapolate_second_order_sw_c(domain):
510    """Wrapper calling C version of extrapolate_second_order_sw
511    """
512    import sys
513    from Numeric import zeros, Float
514
515    N = domain.number_of_elements
516
517    #Shortcuts
518    Stage = domain.quantities['stage']
519    Xmom = domain.quantities['xmomentum']
520#    Ymom = domain.quantities['ymomentum']
521    from shallow_water_ext import extrapolate_second_order_sw
522    extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
523# cant find this in domain                      domain.number_of_boundaries,
524                                domain.centroid_coordinates,
525                                Stage.centroid_values,
526                                Xmom.centroid_values,
527#                                Ymom.centroid_values,
528                                domain.vertex_coordinates,
529                                Stage.vertex_values,
530                                Xmom.vertex_values)#,
531#                                Ymom.vertex_values)
532
533def compute_fluxes_c(domain):
534    """Wrapper calling C version of compute fluxes
535    """
536
537    import sys
538    from Numeric import zeros, Float
539
540    N = domain.number_of_elements
541
542    #Shortcuts
543    Stage = domain.quantities['stage']
544    Xmom = domain.quantities['xmomentum']
545    #Ymom = domain.quantities['ymomentum']
546    Bed = domain.quantities['elevation']
547
548    timestep = float(sys.maxint)
549    from shallow_water_ext import compute_fluxes
550    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
551                                     domain.neighbours,
552                                     #domain.neighbour_edges,
553                                     domain.neighbour_vertices,
554                                     domain.normals,
555                                     #domain.edgelengths,
556                                     #domain.radii,
557                                     domain.areas,
558                                     #Stage.edge_values,
559                                     #Xmom.edge_values,
560                                     #Ymom.edge_values,
561                                     #Bed.edge_values,
562                                     Stage.vertex_values,
563                                     Xmom.vertex_values,
564                                     Bed.vertex_values,
565                                     Stage.boundary_values,
566                                     Xmom.boundary_values,
567                                     #Ymom.boundary_values,
568                                     Stage.explicit_update,
569                                     Xmom.explicit_update,
570                                     #Ymom.explicit_update,
571                                     domain.already_computed_flux)
572
573
574####################################
575
576# Module functions for gradient limiting (distribute_to_vertices_and_edges)
577
578def distribute_to_vertices_and_edges(domain):
579    """Distribution from centroids to vertices specific to the
580    shallow water wave
581    equation.
582
583    It will ensure that h (w-z) is always non-negative even in the
584    presence of steep bed-slopes by taking a weighted average between shallow
585    and deep cases.
586
587    In addition, all conserved quantities get distributed as per either a
588    constant (order==1) or a piecewise linear function (order==2).
589
590    FIXME: more explanation about removal of artificial variability etc
591
592    Precondition:
593      All quantities defined at centroids and bed elevation defined at
594      vertices.
595
596    Postcondition
597      Conserved quantities defined at vertices
598
599    """
600
601    from config import optimised_gradient_limiter
602
603    #Remove very thin layers of water
604    #protect_against_infinitesimal_and_negative_heights(domain)
605
606    #Extrapolate all conserved quantities
607    #if optimised_gradient_limiter:
608    #    #MH090605 if second order,
609    #    #perform the extrapolation and limiting on
610    #    #all of the conserved quantities
611
612    #    if (domain.order == 1):
613    #        for name in domain.conserved_quantities:
614    #            Q = domain.quantities[name]
615    #            Q.extrapolate_first_order()
616    #    elif domain.order == 2:
617    #        domain.extrapolate_second_order_sw()
618    #    else:
619    #        raise 'Unknown order'
620    #else:
621        #old code:
622    for name in domain.conserved_quantities:
623        Q = domain.quantities[name]
624        if domain.order == 1:
625            Q.extrapolate_first_order()
626        elif domain.order == 2:
627            Q.extrapolate_second_order()
628            Q.limit()
629        else:
630            raise 'Unknown order'
631
632    #Take bed elevation into account when water heights are small
633    #balance_deep_and_shallow(domain)
634
635    #Compute edge values by interpolation
636    #for name in domain.conserved_quantities:
637    #    Q = domain.quantities[name]
638    #    Q.interpolate_from_vertices_to_edges()
639
640
641def protect_against_infinitesimal_and_negative_heights(domain):
642    """Protect against infinitesimal heights and associated high velocities
643    """
644
645    #Shortcuts
646    wc = domain.quantities['stage'].centroid_values
647    zc = domain.quantities['elevation'].centroid_values
648    xmomc = domain.quantities['xmomentum'].centroid_values
649#    ymomc = domain.quantities['ymomentum'].centroid_values
650    hc = wc - zc  #Water depths at centroids
651
652    #Update
653    for k in range(domain.number_of_elements):
654
655        if hc[k] < domain.minimum_allowed_height:
656            #Control stage
657            if hc[k] < domain.epsilon:
658                wc[k] = zc[k] # Contain 'lost mass' error
659
660            #Control momentum
661#            xmomc[k] = ymomc[k] = 0.0
662            xmomc[k] = 0.0
663
664def protect_against_infinitesimal_and_negative_heights_c(domain):
665    """Protect against infinitesimal heights and associated high velocities
666    """
667
668    #Shortcuts
669    wc = domain.quantities['stage'].centroid_values
670    zc = domain.quantities['elevation'].centroid_values
671    xmomc = domain.quantities['xmomentum'].centroid_values
672#    ymomc = domain.quantities['ymomentum'].centroid_values
673
674    from shallow_water_ext import protect
675
676    protect(domain.minimum_allowed_height, domain.epsilon,
677            wc, zc, xmomc)#, ymomc)
678
679def h_limiter(domain):
680    """Limit slopes for each volume to eliminate artificial variance
681    introduced by e.g. second order extrapolator
682
683    limit on h = w-z
684
685    This limiter depends on two quantities (w,z) so it resides within
686    this module rather than within quantity.py
687    """
688
689    from Numeric import zeros, Float
690
691    N = domain.number_of_elements
692    beta_h = domain.beta_h
693
694    #Shortcuts
695    wc = domain.quantities['stage'].centroid_values
696    zc = domain.quantities['elevation'].centroid_values
697    hc = wc - zc
698
699    wv = domain.quantities['stage'].vertex_values
700    zv = domain.quantities['elevation'].vertex_values
701    hv = wv-zv
702
703    hvbar = zeros(hv.shape, Float) #h-limited values
704
705    #Find min and max of this and neighbour's centroid values
706    hmax = zeros(hc.shape, Float)
707    hmin = zeros(hc.shape, Float)
708
709    for k in range(N):
710        hmax[k] = hmin[k] = hc[k]
711        #for i in range(3):
712        for i in range(2):   
713            n = domain.neighbours[k,i]
714            if n >= 0:
715                hn = hc[n] #Neighbour's centroid value
716
717                hmin[k] = min(hmin[k], hn)
718                hmax[k] = max(hmax[k], hn)
719
720
721    #Diffences between centroids and maxima/minima
722    dhmax = hmax - hc
723    dhmin = hmin - hc
724
725    #Deltas between vertex and centroid values
726    dh = zeros(hv.shape, Float)
727    #for i in range(3):
728    for i in range(2):
729        dh[:,i] = hv[:,i] - hc
730
731    #Phi limiter
732    for k in range(N):
733
734        #Find the gradient limiter (phi) across vertices
735        phi = 1.0
736        #for i in range(3):
737        for i in range(2):
738            r = 1.0
739            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
740            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
741
742            phi = min( min(r*beta_h, 1), phi )
743
744        #Then update using phi limiter
745        #for i in range(3):
746        for i in range(2):
747            hvbar[k,i] = hc[k] + phi*dh[k,i]
748
749    return hvbar
750
751
752
753def h_limiter_c(domain):
754    """Limit slopes for each volume to eliminate artificial variance
755    introduced by e.g. second order extrapolator
756
757    limit on h = w-z
758
759    This limiter depends on two quantities (w,z) so it resides within
760    this module rather than within quantity.py
761
762    Wrapper for c-extension
763    """
764
765    from Numeric import zeros, Float
766
767    N = domain.number_of_elements
768    beta_h = domain.beta_h
769
770    #Shortcuts
771    wc = domain.quantities['stage'].centroid_values
772    zc = domain.quantities['elevation'].centroid_values
773    hc = wc - zc
774
775    wv = domain.quantities['stage'].vertex_values
776    zv = domain.quantities['elevation'].vertex_values
777    hv = wv - zv
778
779    #Call C-extension
780    from shallow_water_ext import h_limiter_sw as h_limiter
781    hvbar = h_limiter(domain, hc, hv)
782
783    return hvbar
784
785def balance_deep_and_shallow(domain):
786    """Compute linear combination between stage as computed by
787    gradient-limiters limiting using w, and stage computed by
788    gradient-limiters limiting using h (h-limiter).
789    The former takes precedence when heights are large compared to the
790    bed slope while the latter takes precedence when heights are
791    relatively small.  Anything in between is computed as a balanced
792    linear combination in order to avoid numerical disturbances which
793    would otherwise appear as a result of hard switching between
794    modes.
795
796    The h-limiter is always applied irrespective of the order.
797    """
798
799    #Shortcuts
800    wc = domain.quantities['stage'].centroid_values
801    zc = domain.quantities['elevation'].centroid_values
802    hc = wc - zc
803
804    wv = domain.quantities['stage'].vertex_values
805    zv = domain.quantities['elevation'].vertex_values
806    hv = wv-zv
807
808    #Limit h
809    hvbar = h_limiter(domain)
810
811    for k in range(domain.number_of_elements):
812        #Compute maximal variation in bed elevation
813        #  This quantitiy is
814        #    dz = max_i abs(z_i - z_c)
815        #  and it is independent of dimension
816        #  In the 1d case zc = (z0+z1)/2
817        #  In the 2d case zc = (z0+z1+z2)/3
818
819        dz = max(abs(zv[k,0]-zc[k]),
820                 abs(zv[k,1]-zc[k]))#,
821#                 abs(zv[k,2]-zc[k]))
822
823
824        hmin = min( hv[k,:] )
825
826        #Create alpha in [0,1], where alpha==0 means using the h-limited
827        #stage and alpha==1 means using the w-limited stage as
828        #computed by the gradient limiter (both 1st or 2nd order)
829
830        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
831        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
832
833        if dz > 0.0:
834            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
835        else:
836            #Flat bed
837            alpha = 1.0
838
839        #Let
840        #
841        #  wvi be the w-limited stage (wvi = zvi + hvi)
842        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
843        #
844        #
845        #where i=0,1,2 denotes the vertex ids
846        #
847        #Weighted balance between w-limited and h-limited stage is
848        #
849        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
850        #
851        #It follows that the updated wvi is
852        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
853        #
854        # Momentum is balanced between constant and limited
855
856
857        #for i in range(3):
858        #    wv[k,i] = zv[k,i] + hvbar[k,i]
859
860        #return
861
862        if alpha < 1:
863
864            #for i in range(3):
865            for i in range(2):
866                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
867
868            #Momentums at centroids
869            xmomc = domain.quantities['xmomentum'].centroid_values
870#            ymomc = domain.quantities['ymomentum'].centroid_values
871
872            #Momentums at vertices
873            xmomv = domain.quantities['xmomentum'].vertex_values
874#           ymomv = domain.quantities['ymomentum'].vertex_values
875
876            # Update momentum as a linear combination of
877            # xmomc and ymomc (shallow) and momentum
878            # from extrapolator xmomv and ymomv (deep).
879            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
880#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
881
882
883def balance_deep_and_shallow_c(domain):
884    """Wrapper for C implementation
885    """
886
887    #Shortcuts
888    wc = domain.quantities['stage'].centroid_values
889    zc = domain.quantities['elevation'].centroid_values
890    hc = wc - zc
891
892    wv = domain.quantities['stage'].vertex_values
893    zv = domain.quantities['elevation'].vertex_values
894    hv = wv - zv
895
896    #Momentums at centroids
897    xmomc = domain.quantities['xmomentum'].centroid_values
898 #   ymomc = domain.quantities['ymomentum'].centroid_values
899
900    #Momentums at vertices
901    xmomv = domain.quantities['xmomentum'].vertex_values
902#    ymomv = domain.quantities['ymomentum'].vertex_values
903
904    #Limit h
905    hvbar = h_limiter(domain)
906
907    #This is how one would make a first order h_limited value
908    #as in the old balancer (pre 17 Feb 2005):
909    #from Numeric import zeros, Float
910    #hvbar = zeros( (len(hc), 3), Float)
911    #for i in range(3):
912    #    hvbar[:,i] = hc[:]
913
914    from shallow_water_ext import balance_deep_and_shallow
915    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
916    #                         xmomc, ymomc, xmomv, ymomv)
917                             xmomc, xmomv)
918
919
920###############################################
921#Boundaries - specific to the shallow water wave equation
922class Reflective_boundary(Boundary):
923    """Reflective boundary returns same conserved quantities as
924    those present in its neighbour volume but reflected.
925
926    This class is specific to the shallow water equation as it
927    works with the momentum quantities assumed to be the second
928    and third conserved quantities.
929    """
930
931    def __init__(self, domain = None):
932        Boundary.__init__(self)
933
934        if domain is None:
935            msg = 'Domain must be specified for reflective boundary'
936            raise msg
937
938        #Handy shorthands
939        #self.stage   = domain.quantities['stage'].edge_values
940        #self.xmom    = domain.quantities['xmomentum'].edge_values
941        #self.ymom    = domain.quantities['ymomentum'].edge_values
942        #self.normals = domain.normals
943        self.stage   = domain.quantities['stage'].vertex_values
944        self.xmom    = domain.quantities['xmomentum'].vertex_values
945
946        from Numeric import zeros, Float
947        #self.conserved_quantities = zeros(3, Float)
948        self.conserved_quantities = zeros(2, Float)
949
950    def __repr__(self):
951        return 'Reflective_boundary'
952
953
954    def evaluate(self, vol_id, edge_id):
955        """Reflective boundaries reverses the outward momentum
956        of the volume they serve.
957        """
958
959        q = self.conserved_quantities
960        q[0] = self.stage[vol_id, edge_id]
961        q[1] = self.xmom[vol_id, edge_id]
962        #q[2] = self.ymom[vol_id, edge_id]
963
964        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
965        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
966
967
968        #r = rotate(q, normal, direction = 1)
969        #r[1] = -r[1]
970        #q = rotate(r, normal, direction = -1)
971        r = q
972        r[1] = -q[1]
973        q = r
974        #For start interval there is no outward momentum so do not need to
975        #reverse direction in this case
976
977        return q
978
979
980#########################
981#Standard forcing terms:
982#
983def gravity(domain):
984    """Apply gravitational pull in the presence of bed slope
985    """
986
987    from util import gradient
988    from Numeric import zeros, Float, array, sum
989
990    xmom = domain.quantities['xmomentum'].explicit_update
991#    ymom = domain.quantities['ymomentum'].explicit_update
992
993    Stage = domain.quantities['stage']
994    Elevation = domain.quantities['elevation']
995    #h = Stage.edge_values - Elevation.edge_values
996    h = Stage.vertex_values - Elevation.vertex_values
997    v = Elevation.vertex_values
998
999    x = domain.get_vertex_coordinates()
1000    g = domain.g
1001
1002    for k in range(domain.number_of_elements):
1003#        avg_h = sum( h[k,:] )/3
1004        avg_h = sum( h[k,:] )/2
1005
1006        #Compute bed slope
1007        #x0, y0, x1, y1, x2, y2 = x[k,:]
1008        x0, x1 = x[k,:]
1009        #z0, z1, z2 = v[k,:]
1010        z0, z1 = v[k,:]
1011
1012        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
1013        zx = gradient(x0, x1, z0, z1)
1014       
1015        #Update momentum
1016        xmom[k] += -g*zx*avg_h
1017#        ymom[k] += -g*zy*avg_h
1018
1019
1020def graity_c(domain):
1021    """Wrapper calling C version
1022    """
1023
1024    xmom = domain.quantities['xmomentum'].explicit_update
1025#   ymom = domain.quantities['ymomentum'].explicit_update
1026
1027    Stage = domain.quantities['stage']
1028    Elevation = domain.quantities['elevation']
1029    h = Stage.edge_values - Elevation.edge_values
1030    v = Elevation.vertex_values
1031
1032    x = domain.get_vertex_coordinates()
1033    g = domain.g
1034
1035
1036    from shallow_water_ext import gravity
1037    gravity(g, h, v, x, xmom, ymom)
1038
1039
1040def manning_friction(domain):
1041    """Apply (Manning) friction to water momentum
1042    """
1043
1044    from math import sqrt
1045
1046    w = domain.quantities['stage'].centroid_values
1047    z = domain.quantities['elevation'].centroid_values
1048    h = w-z
1049
1050    uh = domain.quantities['xmomentum'].centroid_values
1051    #vh = domain.quantities['ymomentum'].centroid_values
1052    eta = domain.quantities['friction'].centroid_values
1053
1054    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1055    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1056
1057    N = domain.number_of_elements
1058    eps = domain.minimum_allowed_height
1059    g = domain.g
1060
1061    for k in range(N):
1062        if eta[k] >= eps:
1063            if h[k] >= eps:
1064                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1065                S = -g * eta[k]**2 * uh[k]
1066                S /= h[k]**(7.0/3)
1067
1068                #Update momentum
1069                xmom_update[k] += S*uh[k]
1070                #ymom_update[k] += S*vh[k]
1071
1072
1073def manning_friction_c(domain):
1074    """Wrapper for c version
1075    """
1076
1077    xmom = domain.quantities['xmomentum']
1078#    ymom = domain.quantities['ymomentum']
1079
1080    w = domain.quantities['stage'].centroid_values
1081    z = domain.quantities['elevation'].centroid_values
1082
1083    uh = xmom.centroid_values
1084#    vh = ymom.centroid_values
1085    eta = domain.quantities['friction'].centroid_values
1086
1087    xmom_update = xmom.semi_implicit_update
1088#    ymom_update = ymom.semi_implicit_update
1089
1090    N = domain.number_of_elements
1091    eps = domain.minimum_allowed_height
1092    g = domain.g
1093
1094    from shallow_water_ext import manning_friction
1095#    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1096    manning_friction(g, eps, w, z, uh, eta, xmom_update)
1097
1098def linear_friction(domain):
1099    """Apply linear friction to water momentum
1100
1101    Assumes quantity: 'linear_friction' to be present
1102    """
1103
1104    from math import sqrt
1105
1106    w = domain.quantities['stage'].centroid_values
1107    z = domain.quantities['elevation'].centroid_values
1108    h = w-z
1109
1110    uh = domain.quantities['xmomentum'].centroid_values
1111#    vh = domain.quantities['ymomentum'].centroid_values
1112    tau = domain.quantities['linear_friction'].centroid_values
1113
1114    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1115#    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1116
1117    N = domain.number_of_elements
1118    eps = domain.minimum_allowed_height
1119    g = domain.g #Not necessary? Why was this added?
1120
1121    for k in range(N):
1122        if tau[k] >= eps:
1123            if h[k] >= eps:
1124                S = -tau[k]/h[k]
1125
1126                #Update momentum
1127                xmom_update[k] += S*uh[k]
1128 #              ymom_update[k] += S*vh[k]
1129
1130
1131
1132def check_forcefield(f):
1133    """Check that f is either
1134    1: a callable object f(t,x,y), where x and y are vectors
1135       and that it returns an array or a list of same length
1136       as x and y
1137    2: a scalar
1138    """
1139
1140    from Numeric import ones, Float, array
1141
1142
1143    if callable(f):
1144        #N = 3
1145        N = 2
1146        #x = ones(3, Float)
1147        #y = ones(3, Float)
1148        x = ones(2, Float)
1149        #y = ones(2, Float)
1150       
1151        try:
1152            #q = f(1.0, x=x, y=y)
1153            q = f(1.0, x=x)
1154        except Exception, e:
1155            msg = 'Function %s could not be executed:\n%s' %(f, e)
1156            #FIXME: Reconsider this semantics
1157            raise msg
1158
1159        try:
1160            q = array(q).astype(Float)
1161        except:
1162            msg = 'Return value from vector function %s could ' %f
1163            msg += 'not be converted into a Numeric array of floats.\n'
1164            msg += 'Specified function should return either list or array.'
1165            raise msg
1166
1167        #Is this really what we want?
1168        msg = 'Return vector from function %s ' %f
1169        msg += 'must have same lenght as input vectors'
1170        assert len(q) == N, msg
1171
1172    else:
1173        try:
1174            f = float(f)
1175        except:
1176            msg = 'Force field %s must be either a scalar' %f
1177            msg += ' or a vector function'
1178            raise msg
1179    return f
1180
1181class Wind_stress:
1182    """Apply wind stress to water momentum in terms of
1183    wind speed [m/s] and wind direction [degrees]
1184    """
1185
1186    def __init__(self, *args, **kwargs):
1187        """Initialise windfield from wind speed s [m/s]
1188        and wind direction phi [degrees]
1189
1190        Inputs v and phi can be either scalars or Python functions, e.g.
1191
1192        W = Wind_stress(10, 178)
1193
1194        #FIXME - 'normal' degrees are assumed for now, i.e. the
1195        vector (1,0) has zero degrees.
1196        We may need to convert from 'compass' degrees later on and also
1197        map from True north to grid north.
1198
1199        Arguments can also be Python functions of t,x,y as in
1200
1201        def speed(t,x,y):
1202            ...
1203            return s
1204
1205        def angle(t,x,y):
1206            ...
1207            return phi
1208
1209        where x and y are vectors.
1210
1211        and then pass the functions in
1212
1213        W = Wind_stress(speed, angle)
1214
1215        The instantiated object W can be appended to the list of
1216        forcing_terms as in
1217
1218        Alternatively, one vector valued function for (speed, angle)
1219        can be applied, providing both quantities simultaneously.
1220        As in
1221        W = Wind_stress(F), where returns (speed, angle) for each t.
1222
1223        domain.forcing_terms.append(W)
1224        """
1225
1226        from config import rho_a, rho_w, eta_w
1227        from Numeric import array, Float
1228
1229        if len(args) == 2:
1230            s = args[0]
1231            phi = args[1]
1232        elif len(args) == 1:
1233            #Assume vector function returning (s, phi)(t,x,y)
1234            vector_function = args[0]
1235            #s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1236            #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1237            s = lambda t,x: vector_function(t,x=x)[0]
1238            phi = lambda t,x: vector_function(t,x=x)[1]
1239        else:
1240           #Assume info is in 2 keyword arguments
1241
1242           if len(kwargs) == 2:
1243               s = kwargs['s']
1244               phi = kwargs['phi']
1245           else:
1246               raise 'Assumes two keyword arguments: s=..., phi=....'
1247
1248        print 'phi', phi
1249        self.speed = check_forcefield(s)
1250        self.phi = check_forcefield(phi)
1251
1252        self.const = eta_w*rho_a/rho_w
1253
1254
1255    def __call__(self, domain):
1256        """Evaluate windfield based on values found in domain
1257        """
1258
1259        from math import pi, cos, sin, sqrt
1260        from Numeric import Float, ones, ArrayType
1261
1262        xmom_update = domain.quantities['xmomentum'].explicit_update
1263        #ymom_update = domain.quantities['ymomentum'].explicit_update
1264
1265        N = domain.number_of_elements
1266        t = domain.time
1267
1268        if callable(self.speed):
1269            xc = domain.get_centroid_coordinates()
1270            #s_vec = self.speed(t, xc[:,0], xc[:,1])
1271            s_vec = self.speed(t, xc)
1272        else:
1273            #Assume s is a scalar
1274
1275            try:
1276                s_vec = self.speed * ones(N, Float)
1277            except:
1278                msg = 'Speed must be either callable or a scalar: %s' %self.s
1279                raise msg
1280
1281
1282        if callable(self.phi):
1283            xc = domain.get_centroid_coordinates()
1284            #phi_vec = self.phi(t, xc[:,0], xc[:,1])
1285            phi_vec = self.phi(t, xc)
1286        else:
1287            #Assume phi is a scalar
1288
1289            try:
1290                phi_vec = self.phi * ones(N, Float)
1291            except:
1292                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1293                raise msg
1294
1295        #assign_windfield_values(xmom_update, ymom_update,
1296        #                        s_vec, phi_vec, self.const)
1297        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
1298
1299
1300#def assign_windfield_values(xmom_update, ymom_update,
1301#                            s_vec, phi_vec, const):
1302def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
1303    """Python version of assigning wind field to update vectors.
1304    A c version also exists (for speed)
1305    """
1306    from math import pi, cos, sin, sqrt
1307
1308    N = len(s_vec)
1309    for k in range(N):
1310        s = s_vec[k]
1311        phi = phi_vec[k]
1312
1313        #Convert to radians
1314        phi = phi*pi/180
1315
1316        #Compute velocity vector (u, v)
1317        u = s*cos(phi)
1318        v = s*sin(phi)
1319
1320        #Compute wind stress
1321        #S = const * sqrt(u**2 + v**2)
1322        S = const * u
1323        xmom_update[k] += S*u
1324        #ymom_update[k] += S*v
Note: See TracBrowser for help on using the repository browser.