source: anuga_work/development/anuga_1d/shallow_water_h.py @ 7818

Last change on this file since 7818 was 5535, checked in by steve, 16 years ago

Changed name of shallow_water_1d to anuga_1d

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