source: trunk/anuga_work/development/sudi/sw_1d/avalanche/B_momentum/shallow_water_domain_avalanche.py @ 7915

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