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

Last change on this file since 3158 was 2791, checked in by jakeman, 19 years ago

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