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

Last change on this file since 2705 was 2705, checked in by jakeman, 18 years ago

domain with working neighbour structure

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