source: anuga_work/development/shallow_water_1d/shallow_water_h.py @ 4959

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

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

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