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

Last change on this file since 3365 was 3362, checked in by jakeman, 19 years ago

New files dam2.py and dry_dam.py further test numerical solution
against analytical solution of
Stoker 57. Limiting is still not working correctly.
File size: 40.6 KB
Line 
1u"""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    enter = True
435    for k in range(N):
436
437        flux[:] = 0.  #Reset work array
438        #for i in range(3):
439        for i in range(2):
440            #Quantities inside volume facing neighbour i
441            #ql[0] = stage[k, i]
442            #ql[1] = xmom[k, i]
443            ql = [stage[k, i], xmom[k, i]]
444            zl = bed[k, i]
445
446            #Quantities at neighbour on nearest face
447            n = domain.neighbours[k,i]
448            if n < 0:
449                m = -n-1 #Convert negative flag to index
450                qr[0] = stage_bdry[m]
451                qr[1] = xmom_bdry[m]
452                zr = zl #Extend bed elevation to boundary
453            else:
454                #m = domain.neighbour_edges[k,i]
455                m = domain.neighbour_vertices[k,i]
456                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
457                qr[0] = stage[n, m]
458                qr[1] =  xmom[n, m]
459                zr = bed[n, m]
460
461
462            #Outward pointing normal vector
463            normal = domain.normals[k, i]
464       
465            #Flux computation using provided function
466            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
467            #print 'ql',ql
468            #print 'qr',qr
469           
470            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
471            #print 'edgeflux', edgeflux
472
473            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
474            # flux = edgefluxleft - edgefluxright
475            flux -= edgeflux #* domain.edgelengths[k,i]
476            #Update optimal_timestep
477            try:
478                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
479                #timestep = 0.01
480           
481                timestep = min(timestep, 0.5*domain.areas[k]/max_speed)
482                if (timestep < 1e-6) & (enter == True):
483                    #print "domain.order", domain.order
484                    #domain.write_time()
485                    print "cell number", k
486                    print "timestep", timestep
487                    print 'max_speed', max_speed
488                    s = domain.quantities['stage']
489                    s = s.centroid_values
490                    xm = domain.quantities['xmomentum']
491                    xm = xm.centroid_values
492                    print 'h', s[k]
493                    print 'xm', xm[k]
494                    print 'u', xm[k]/s[k]
495                    enter = False
496           
497            except ZeroDivisionError:
498                pass
499
500        #Normalise by area and store for when all conserved
501        #quantities get updated
502        flux /= domain.areas[k]
503        # ADD ABOVE LINE AGAIN
504        Stage.explicit_update[k] = flux[0]
505        Xmom.explicit_update[k] = flux[1]
506        #Ymom.explicit_update[k] = flux[2]
507
508
509    domain.timestep = timestep
510
511#see comments in the corresponding method in shallow_water_ext.c
512def extrapolate_second_order_sw_c(domain):
513    """Wrapper calling C version of extrapolate_second_order_sw
514    """
515    import sys
516    from Numeric import zeros, Float
517
518    N = domain.number_of_elements
519
520    #Shortcuts
521    Stage = domain.quantities['stage']
522    Xmom = domain.quantities['xmomentum']
523#    Ymom = domain.quantities['ymomentum']
524    from shallow_water_ext import extrapolate_second_order_sw
525    extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
526# cant find this in domain                      domain.number_of_boundaries,
527                                domain.centroid_coordinates,
528                                Stage.centroid_values,
529                                Xmom.centroid_values,
530#                                Ymom.centroid_values,
531                                domain.vertex_coordinates,
532                                Stage.vertex_values,
533                                Xmom.vertex_values)#,
534#                                Ymom.vertex_values)
535
536def compute_fluxes_c(domain):
537    """Wrapper calling C version of compute fluxes
538    """
539
540    import sys
541    from Numeric import zeros, Float
542
543    N = domain.number_of_elements
544
545    #Shortcuts
546    Stage = domain.quantities['stage']
547    Xmom = domain.quantities['xmomentum']
548    #Ymom = domain.quantities['ymomentum']
549    Bed = domain.quantities['elevation']
550
551    timestep = float(sys.maxint)
552    from shallow_water_ext import compute_fluxes
553    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
554                                     domain.neighbours,
555                                     #domain.neighbour_edges,
556                                     domain.neighbour_vertices,
557                                     domain.normals,
558                                     #domain.edgelengths,
559                                     #domain.radii,
560                                     domain.areas,
561                                     #Stage.edge_values,
562                                     #Xmom.edge_values,
563                                     #Ymom.edge_values,
564                                     #Bed.edge_values,
565                                     Stage.vertex_values,
566                                     Xmom.vertex_values,
567                                     Bed.vertex_values,
568                                     Stage.boundary_values,
569                                     Xmom.boundary_values,
570                                     #Ymom.boundary_values,
571                                     Stage.explicit_update,
572                                     Xmom.explicit_update,
573                                     #Ymom.explicit_update,
574                                     domain.already_computed_flux)
575
576
577####################################
578
579# Module functions for gradient limiting (distribute_to_vertices_and_edges)
580
581def distribute_to_vertices_and_edges(domain):
582    """Distribution from centroids to vertices specific to the
583    shallow water wave
584    equation.
585
586    It will ensure that h (w-z) is always non-negative even in the
587    presence of steep bed-slopes by taking a weighted average between shallow
588    and deep cases.
589
590    In addition, all conserved quantities get distributed as per either a
591    constant (order==1) or a piecewise linear function (order==2).
592
593    FIXME: more explanation about removal of artificial variability etc
594
595    Precondition:
596      All quantities defined at centroids and bed elevation defined at
597      vertices.
598
599    Postcondition
600      Conserved quantities defined at vertices
601
602    """
603
604    #from config import optimised_gradient_limiter
605
606    #Remove very thin layers of water
607    protect_against_infinitesimal_and_negative_heights(domain)
608
609    #Extrapolate all conserved quantities
610    #if optimised_gradient_limiter:
611    #    #MH090605 if second order,
612    #    #perform the extrapolation and limiting on
613    #    #all of the conserved quantities
614
615    #    if (domain.order == 1):
616    #        for name in domain.conserved_quantities:
617    #            Q = domain.quantities[name]
618    #            Q.extrapolate_first_order()
619    #    elif domain.order == 2:
620    #        domain.extrapolate_second_order_sw()
621    #    else:
622    #        raise 'Unknown order'
623    #else:
624        #old code:
625    for name in domain.conserved_quantities:
626        Q = domain.quantities[name]
627        if domain.order == 1:
628            Q.extrapolate_first_order()
629        elif domain.order == 2:
630            Q.extrapolate_second_order()
631            Q.limit()
632        else:
633            raise 'Unknown order'
634
635    #Take bed elevation into account when water heights are small
636    balance_deep_and_shallow(domain)
637
638    #Compute edge values by interpolation
639    #for name in domain.conserved_quantities:
640    #    Q = domain.quantities[name]
641    #    Q.interpolate_from_vertices_to_edges()
642
643
644def protect_against_infinitesimal_and_negative_heights(domain):
645    """Protect against infinitesimal heights and associated high velocities
646    """
647
648    #Shortcuts
649    wc = domain.quantities['stage'].centroid_values
650    zc = domain.quantities['elevation'].centroid_values
651    xmomc = domain.quantities['xmomentum'].centroid_values
652#    ymomc = domain.quantities['ymomentum'].centroid_values
653    hc = wc - zc  #Water depths at centroids
654
655    #Update
656    for k in range(domain.number_of_elements):
657
658        if hc[k] < domain.minimum_allowed_height:
659            #Control stage
660            if hc[k] < domain.epsilon:
661                wc[k] = zc[k] # Contain 'lost mass' error
662
663            #Control momentum
664            #xmomc[k] = ymomc[k] = 0.0
665            xmomc[k] = 0.0
666
667def protect_against_infinitesimal_and_negative_heights_c(domain):
668    """Protect against infinitesimal heights and associated high velocities
669    """
670
671    #Shortcuts
672    wc = domain.quantities['stage'].centroid_values
673    zc = domain.quantities['elevation'].centroid_values
674    xmomc = domain.quantities['xmomentum'].centroid_values
675#    ymomc = domain.quantities['ymomentum'].centroid_values
676
677    from shallow_water_ext import protect
678
679    protect(domain.minimum_allowed_height, domain.epsilon,
680            wc, zc, xmomc)#, ymomc)
681
682def h_limiter(domain):
683    """Limit slopes for each volume to eliminate artificial variance
684    introduced by e.g. second order extrapolator
685
686    limit on h = w-z
687
688    This limiter depends on two quantities (w,z) so it resides within
689    this module rather than within quantity.py
690    """
691
692    from Numeric import zeros, Float
693
694    N = domain.number_of_elements
695    beta_h = domain.beta_h
696
697    #Shortcuts
698    wc = domain.quantities['stage'].centroid_values
699    zc = domain.quantities['elevation'].centroid_values
700    hc = wc - zc
701
702    wv = domain.quantities['stage'].vertex_values
703    zv = domain.quantities['elevation'].vertex_values
704    hv = wv-zv
705
706    hvbar = zeros(hv.shape, Float) #h-limited values
707
708    #Find min and max of this and neighbour's centroid values
709    hmax = zeros(hc.shape, Float)
710    hmin = zeros(hc.shape, Float)
711
712    for k in range(N):
713        hmax[k] = hmin[k] = hc[k]
714        #for i in range(3):
715        for i in range(2):   
716            n = domain.neighbours[k,i]
717            if n >= 0:
718                hn = hc[n] #Neighbour's centroid value
719
720                hmin[k] = min(hmin[k], hn)
721                hmax[k] = max(hmax[k], hn)
722
723
724    #Diffences between centroids and maxima/minima
725    dhmax = hmax - hc
726    dhmin = hmin - hc
727
728    #Deltas between vertex and centroid values
729    dh = zeros(hv.shape, Float)
730    #for i in range(3):
731    for i in range(2):
732        dh[:,i] = hv[:,i] - hc
733
734    #Phi limiter
735    for k in range(N):
736
737        #Find the gradient limiter (phi) across vertices
738        phi = 1.0
739        #for i in range(3):
740        for i in range(2):
741            r = 1.0
742            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
743            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
744
745            phi = min( min(r*beta_h, 1), phi )
746
747        #Then update using phi limiter
748        #for i in range(3):
749        for i in range(2):
750            hvbar[k,i] = hc[k] + phi*dh[k,i]
751
752    return hvbar
753
754
755
756def h_limiter_c(domain):
757    """Limit slopes for each volume to eliminate artificial variance
758    introduced by e.g. second order extrapolator
759
760    limit on h = w-z
761
762    This limiter depends on two quantities (w,z) so it resides within
763    this module rather than within quantity.py
764
765    Wrapper for c-extension
766    """
767
768    from Numeric import zeros, Float
769
770    N = domain.number_of_elements
771    beta_h = domain.beta_h
772
773    #Shortcuts
774    wc = domain.quantities['stage'].centroid_values
775    zc = domain.quantities['elevation'].centroid_values
776    hc = wc - zc
777
778    wv = domain.quantities['stage'].vertex_values
779    zv = domain.quantities['elevation'].vertex_values
780    hv = wv - zv
781
782    #Call C-extension
783    from shallow_water_ext import h_limiter_sw as h_limiter
784    hvbar = h_limiter(domain, hc, hv)
785
786    return hvbar
787
788def balance_deep_and_shallow(domain):
789    """Compute linear combination between stage as computed by
790    gradient-limiters limiting using w, and stage computed by
791    gradient-limiters limiting using h (h-limiter).
792    The former takes precedence when heights are large compared to the
793    bed slope while the latter takes precedence when heights are
794    relatively small.  Anything in between is computed as a balanced
795    linear combination in order to avoid numerical disturbances which
796    would otherwise appear as a result of hard switching between
797    modes.
798
799    The h-limiter is always applied irrespective of the order.
800    """
801
802    #Shortcuts
803    wc = domain.quantities['stage'].centroid_values
804    zc = domain.quantities['elevation'].centroid_values
805    hc = wc - zc
806
807    wv = domain.quantities['stage'].vertex_values
808    zv = domain.quantities['elevation'].vertex_values
809    hv = wv-zv
810
811    #Limit h
812    hvbar = h_limiter(domain)
813
814    for k in range(domain.number_of_elements):
815        #Compute maximal variation in bed elevation
816        #  This quantitiy is
817        #    dz = max_i abs(z_i - z_c)
818        #  and it is independent of dimension
819        #  In the 1d case zc = (z0+z1)/2
820        #  In the 2d case zc = (z0+z1+z2)/3
821
822        dz = max(abs(zv[k,0]-zc[k]),
823                 abs(zv[k,1]-zc[k]))#,
824#                 abs(zv[k,2]-zc[k]))
825
826
827        hmin = min( hv[k,:] )
828
829        #Create alpha in [0,1], where alpha==0 means using the h-limited
830        #stage and alpha==1 means using the w-limited stage as
831        #computed by the gradient limiter (both 1st or 2nd order)
832
833        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
834        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
835
836        if dz > 0.0:
837            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
838        else:
839            #Flat bed
840            alpha = 1.0
841
842        #Let
843        #
844        #  wvi be the w-limited stage (wvi = zvi + hvi)
845        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
846        #
847        #
848        #where i=0,1,2 denotes the vertex ids
849        #
850        #Weighted balance between w-limited and h-limited stage is
851        #
852        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
853        #
854        #It follows that the updated wvi is
855        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
856        #
857        # Momentum is balanced between constant and limited
858
859
860        #for i in range(3):
861        #    wv[k,i] = zv[k,i] + hvbar[k,i]
862
863        #return
864
865        if alpha < 1:
866
867            #for i in range(3):
868            for i in range(2):
869                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
870
871            #Momentums at centroids
872            xmomc = domain.quantities['xmomentum'].centroid_values
873#            ymomc = domain.quantities['ymomentum'].centroid_values
874
875            #Momentums at vertices
876            xmomv = domain.quantities['xmomentum'].vertex_values
877#           ymomv = domain.quantities['ymomentum'].vertex_values
878
879            # Update momentum as a linear combination of
880            # xmomc and ymomc (shallow) and momentum
881            # from extrapolator xmomv and ymomv (deep).
882            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
883#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
884
885
886def balance_deep_and_shallow_c(domain):
887    """Wrapper for C implementation
888    """
889
890    #Shortcuts
891    wc = domain.quantities['stage'].centroid_values
892    zc = domain.quantities['elevation'].centroid_values
893    hc = wc - zc
894
895    wv = domain.quantities['stage'].vertex_values
896    zv = domain.quantities['elevation'].vertex_values
897    hv = wv - zv
898
899    #Momentums at centroids
900    xmomc = domain.quantities['xmomentum'].centroid_values
901 #   ymomc = domain.quantities['ymomentum'].centroid_values
902
903    #Momentums at vertices
904    xmomv = domain.quantities['xmomentum'].vertex_values
905#    ymomv = domain.quantities['ymomentum'].vertex_values
906
907    #Limit h
908    hvbar = h_limiter(domain)
909
910    #This is how one would make a first order h_limited value
911    #as in the old balancer (pre 17 Feb 2005):
912    #from Numeric import zeros, Float
913    #hvbar = zeros( (len(hc), 3), Float)
914    #for i in range(3):
915    #    hvbar[:,i] = hc[:]
916
917    from shallow_water_ext import balance_deep_and_shallow
918    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
919    #                         xmomc, ymomc, xmomv, ymomv)
920                             xmomc, xmomv)
921
922
923###############################################
924#Boundaries - specific to the shallow water wave equation
925class Reflective_boundary(Boundary):
926    """Reflective boundary returns same conserved quantities as
927    those present in its neighbour volume but reflected.
928
929    This class is specific to the shallow water equation as it
930    works with the momentum quantities assumed to be the second
931    and third conserved quantities.
932    """
933
934    def __init__(self, domain = None):
935        Boundary.__init__(self)
936
937        if domain is None:
938            msg = 'Domain must be specified for reflective boundary'
939            raise msg
940
941        #Handy shorthands
942        #self.stage   = domain.quantities['stage'].edge_values
943        #self.xmom    = domain.quantities['xmomentum'].edge_values
944        #self.ymom    = domain.quantities['ymomentum'].edge_values
945        #self.normals = domain.normals
946        self.stage   = domain.quantities['stage'].vertex_values
947        self.xmom    = domain.quantities['xmomentum'].vertex_values
948
949        from Numeric import zeros, Float
950        #self.conserved_quantities = zeros(3, Float)
951        self.conserved_quantities = zeros(2, Float)
952
953    def __repr__(self):
954        return 'Reflective_boundary'
955
956
957    def evaluate(self, vol_id, edge_id):
958        """Reflective boundaries reverses the outward momentum
959        of the volume they serve.
960        """
961
962        q = self.conserved_quantities
963        q[0] = self.stage[vol_id, edge_id]
964        q[1] = self.xmom[vol_id, edge_id]
965        #q[2] = self.ymom[vol_id, edge_id]
966
967        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
968        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
969
970
971        #r = rotate(q, normal, direction = 1)
972        #r[1] = -r[1]
973        #q = rotate(r, normal, direction = -1)
974        r = q
975        r[1] = -q[1]
976        q = r
977        #For start interval there is no outward momentum so do not need to
978        #reverse direction in this case
979
980        return q
981
982
983#########################
984#Standard forcing terms:
985#
986def gravity(domain):
987    """Apply gravitational pull in the presence of bed slope
988    """
989
990    from util import gradient
991    from Numeric import zeros, Float, array, sum
992
993    xmom = domain.quantities['xmomentum'].explicit_update
994#    ymom = domain.quantities['ymomentum'].explicit_update
995
996    Stage = domain.quantities['stage']
997    Elevation = domain.quantities['elevation']
998    #h = Stage.edge_values - Elevation.edge_values
999    h = Stage.vertex_values - Elevation.vertex_values
1000    v = Elevation.vertex_values
1001
1002    x = domain.get_vertex_coordinates()
1003    g = domain.g
1004
1005    for k in range(domain.number_of_elements):
1006#        avg_h = sum( h[k,:] )/3
1007        avg_h = sum( h[k,:] )/2
1008
1009        #Compute bed slope
1010        #x0, y0, x1, y1, x2, y2 = x[k,:]
1011        x0, x1 = x[k,:]
1012        #z0, z1, z2 = v[k,:]
1013        z0, z1 = v[k,:]
1014
1015        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
1016        zx = gradient(x0, x1, z0, z1)
1017       
1018        #Update momentum
1019        xmom[k] += -g*zx*avg_h
1020#        ymom[k] += -g*zy*avg_h
1021
1022
1023def graity_c(domain):
1024    """Wrapper calling C version
1025    """
1026
1027    xmom = domain.quantities['xmomentum'].explicit_update
1028#   ymom = domain.quantities['ymomentum'].explicit_update
1029
1030    Stage = domain.quantities['stage']
1031    Elevation = domain.quantities['elevation']
1032    h = Stage.edge_values - Elevation.edge_values
1033    v = Elevation.vertex_values
1034
1035    x = domain.get_vertex_coordinates()
1036    g = domain.g
1037
1038
1039    from shallow_water_ext import gravity
1040    gravity(g, h, v, x, xmom, ymom)
1041
1042
1043def manning_friction(domain):
1044    """Apply (Manning) friction to water momentum
1045    """
1046
1047    from math import sqrt
1048
1049    w = domain.quantities['stage'].centroid_values
1050    z = domain.quantities['elevation'].centroid_values
1051    h = w-z
1052
1053    uh = domain.quantities['xmomentum'].centroid_values
1054    #vh = domain.quantities['ymomentum'].centroid_values
1055    eta = domain.quantities['friction'].centroid_values
1056
1057    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1058    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1059
1060    N = domain.number_of_elements
1061    eps = domain.minimum_allowed_height
1062    g = domain.g
1063
1064    for k in range(N):
1065        if eta[k] >= eps:
1066            if h[k] >= eps:
1067                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1068                S = -g * eta[k]**2 * uh[k]
1069                S /= h[k]**(7.0/3)
1070
1071                #Update momentum
1072                xmom_update[k] += S*uh[k]
1073                #ymom_update[k] += S*vh[k]
1074
1075
1076def manning_friction_c(domain):
1077    """Wrapper for c version
1078    """
1079
1080    xmom = domain.quantities['xmomentum']
1081#    ymom = domain.quantities['ymomentum']
1082
1083    w = domain.quantities['stage'].centroid_values
1084    z = domain.quantities['elevation'].centroid_values
1085
1086    uh = xmom.centroid_values
1087#    vh = ymom.centroid_values
1088    eta = domain.quantities['friction'].centroid_values
1089
1090    xmom_update = xmom.semi_implicit_update
1091#    ymom_update = ymom.semi_implicit_update
1092
1093    N = domain.number_of_elements
1094    eps = domain.minimum_allowed_height
1095    g = domain.g
1096
1097    from shallow_water_ext import manning_friction
1098#    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1099    manning_friction(g, eps, w, z, uh, eta, xmom_update)
1100
1101def linear_friction(domain):
1102    """Apply linear friction to water momentum
1103
1104    Assumes quantity: 'linear_friction' to be present
1105    """
1106
1107    from math import sqrt
1108
1109    w = domain.quantities['stage'].centroid_values
1110    z = domain.quantities['elevation'].centroid_values
1111    h = w-z
1112
1113    uh = domain.quantities['xmomentum'].centroid_values
1114#    vh = domain.quantities['ymomentum'].centroid_values
1115    tau = domain.quantities['linear_friction'].centroid_values
1116
1117    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1118#    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1119
1120    N = domain.number_of_elements
1121    eps = domain.minimum_allowed_height
1122    g = domain.g #Not necessary? Why was this added?
1123
1124    for k in range(N):
1125        if tau[k] >= eps:
1126            if h[k] >= eps:
1127                S = -tau[k]/h[k]
1128
1129                #Update momentum
1130                xmom_update[k] += S*uh[k]
1131 #              ymom_update[k] += S*vh[k]
1132
1133
1134
1135def check_forcefield(f):
1136    """Check that f is either
1137    1: a callable object f(t,x,y), where x and y are vectors
1138       and that it returns an array or a list of same length
1139       as x and y
1140    2: a scalar
1141    """
1142
1143    from Numeric import ones, Float, array
1144
1145
1146    if callable(f):
1147        #N = 3
1148        N = 2
1149        #x = ones(3, Float)
1150        #y = ones(3, Float)
1151        x = ones(2, Float)
1152        #y = ones(2, Float)
1153       
1154        try:
1155            #q = f(1.0, x=x, y=y)
1156            q = f(1.0, x=x)
1157        except Exception, e:
1158            msg = 'Function %s could not be executed:\n%s' %(f, e)
1159            #FIXME: Reconsider this semantics
1160            raise msg
1161
1162        try:
1163            q = array(q).astype(Float)
1164        except:
1165            msg = 'Return value from vector function %s could ' %f
1166            msg += 'not be converted into a Numeric array of floats.\n'
1167            msg += 'Specified function should return either list or array.'
1168            raise msg
1169
1170        #Is this really what we want?
1171        msg = 'Return vector from function %s ' %f
1172        msg += 'must have same lenght as input vectors'
1173        assert len(q) == N, msg
1174
1175    else:
1176        try:
1177            f = float(f)
1178        except:
1179            msg = 'Force field %s must be either a scalar' %f
1180            msg += ' or a vector function'
1181            raise msg
1182    return f
1183
1184class Wind_stress:
1185    """Apply wind stress to water momentum in terms of
1186    wind speed [m/s] and wind direction [degrees]
1187    """
1188
1189    def __init__(self, *args, **kwargs):
1190        """Initialise windfield from wind speed s [m/s]
1191        and wind direction phi [degrees]
1192
1193        Inputs v and phi can be either scalars or Python functions, e.g.
1194
1195        W = Wind_stress(10, 178)
1196
1197        #FIXME - 'normal' degrees are assumed for now, i.e. the
1198        vector (1,0) has zero degrees.
1199        We may need to convert from 'compass' degrees later on and also
1200        map from True north to grid north.
1201
1202        Arguments can also be Python functions of t,x,y as in
1203
1204        def speed(t,x,y):
1205            ...
1206            return s
1207
1208        def angle(t,x,y):
1209            ...
1210            return phi
1211
1212        where x and y are vectors.
1213
1214        and then pass the functions in
1215
1216        W = Wind_stress(speed, angle)
1217
1218        The instantiated object W can be appended to the list of
1219        forcing_terms as in
1220
1221        Alternatively, one vector valued function for (speed, angle)
1222        can be applied, providing both quantities simultaneously.
1223        As in
1224        W = Wind_stress(F), where returns (speed, angle) for each t.
1225
1226        domain.forcing_terms.append(W)
1227        """
1228
1229        from config import rho_a, rho_w, eta_w
1230        from Numeric import array, Float
1231
1232        if len(args) == 2:
1233            s = args[0]
1234            phi = args[1]
1235        elif len(args) == 1:
1236            #Assume vector function returning (s, phi)(t,x,y)
1237            vector_function = args[0]
1238            #s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1239            #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1240            s = lambda t,x: vector_function(t,x=x)[0]
1241            phi = lambda t,x: vector_function(t,x=x)[1]
1242        else:
1243           #Assume info is in 2 keyword arguments
1244
1245           if len(kwargs) == 2:
1246               s = kwargs['s']
1247               phi = kwargs['phi']
1248           else:
1249               raise 'Assumes two keyword arguments: s=..., phi=....'
1250
1251        print 'phi', phi
1252        self.speed = check_forcefield(s)
1253        self.phi = check_forcefield(phi)
1254
1255        self.const = eta_w*rho_a/rho_w
1256
1257
1258    def __call__(self, domain):
1259        """Evaluate windfield based on values found in domain
1260        """
1261
1262        from math import pi, cos, sin, sqrt
1263        from Numeric import Float, ones, ArrayType
1264
1265        xmom_update = domain.quantities['xmomentum'].explicit_update
1266        #ymom_update = domain.quantities['ymomentum'].explicit_update
1267
1268        N = domain.number_of_elements
1269        t = domain.time
1270
1271        if callable(self.speed):
1272            xc = domain.get_centroid_coordinates()
1273            #s_vec = self.speed(t, xc[:,0], xc[:,1])
1274            s_vec = self.speed(t, xc)
1275        else:
1276            #Assume s is a scalar
1277
1278            try:
1279                s_vec = self.speed * ones(N, Float)
1280            except:
1281                msg = 'Speed must be either callable or a scalar: %s' %self.s
1282                raise msg
1283
1284
1285        if callable(self.phi):
1286            xc = domain.get_centroid_coordinates()
1287            #phi_vec = self.phi(t, xc[:,0], xc[:,1])
1288            phi_vec = self.phi(t, xc)
1289        else:
1290            #Assume phi is a scalar
1291
1292            try:
1293                phi_vec = self.phi * ones(N, Float)
1294            except:
1295                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1296                raise msg
1297
1298        #assign_windfield_values(xmom_update, ymom_update,
1299        #                        s_vec, phi_vec, self.const)
1300        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
1301
1302
1303#def assign_windfield_values(xmom_update, ymom_update,
1304#                            s_vec, phi_vec, const):
1305def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
1306    """Python version of assigning wind field to update vectors.
1307    A c version also exists (for speed)
1308    """
1309    from math import pi, cos, sin, sqrt
1310
1311    N = len(s_vec)
1312    for k in range(N):
1313        s = s_vec[k]
1314        phi = phi_vec[k]
1315
1316        #Convert to radians
1317        phi = phi*pi/180
1318
1319        #Compute velocity vector (u, v)
1320        u = s*cos(phi)
1321        v = s*sin(phi)
1322
1323        #Compute wind stress
1324        #S = const * sqrt(u**2 + v**2)
1325        S = const * u
1326        xmom_update[k] += S*u
1327        #ymom_update[k] += S*v
Note: See TracBrowser for help on using the repository browser.