source: trunk/anuga_work/development/sudi/sw_1d/avalanche/A_velocity/shallow_water_domain_avalanche.py @ 7914

Last change on this file since 7914 was 7914, checked in by mungkasi, 14 years ago

Modifying debris avalanche codes using s-h-v reconstruction.

File size: 38.5 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
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] <= epsilon:
611            uh_C[i] = 0.0
612            u_C[i] = 0.0
613            w_C[i] = z_C[i]
614        else:
615            u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
616
617    for name in ['stage', 'height', 'velocity']:
618        Q = domain.quantities[name]
619        if domain.order == 1:
620            Q.extrapolate_first_order()
621        elif domain.order == 2:
622            Q.extrapolate_second_order()
623        else:
624            raise 'Unknown order'
625
626    stage_V = domain.quantities['stage'].vertex_values
627    bed_V   = domain.quantities['elevation'].vertex_values
628    h_V     = domain.quantities['height'].vertex_values
629    u_V     = domain.quantities['velocity'].vertex_values               
630    xmom_V  = domain.quantities['xmomentum'].vertex_values
631
632    bed_V[:] = stage_V - h_V
633    xmom_V[:] = u_V * h_V
634   
635    return
636
637
638
639def distribute_to_vertices_and_edges_shm(domain):
640    # shm stands for STAGE, HEIGHT, MOMENTUM
641    """Distribution from centroids to vertices specific to the
642    shallow water wave
643    equation.
644
645    It will ensure that h (w-z) is always non-negative even in the
646    presence of steep bed-slopes by taking a weighted average between shallow
647    and deep cases.
648
649    In addition, all conserved quantities get distributed as per either a
650    constant (order==1) or a piecewise linear function (order==2).
651
652    FIXME: more explanation about removal of artificial variability etc
653
654    Precondition:
655      All quantities defined at centroids and bed elevation defined at
656      vertices.
657
658    Postcondition
659      Conserved quantities defined at vertices
660
661    """
662
663    #from config import optimised_gradient_limiter
664
665    #Remove very thin layers of water
666    #protect_against_infinitesimal_and_negative_heights(domain) 
667
668    import sys
669    from Numeric import Float
670    from numpy import zeros
671    from config import epsilon, h0
672
673    N = domain.number_of_elements
674
675    #Shortcuts
676    Stage = domain.quantities['stage']
677    Xmom = domain.quantities['xmomentum']
678    Bed = domain.quantities['elevation']
679    Height = domain.quantities['height']
680    Velocity = domain.quantities['velocity']
681
682    #Arrays   
683    w_C   = Stage.centroid_values   
684    uh_C  = Xmom.centroid_values   
685    z_C   = Bed.centroid_values
686    h_C   = Height.centroid_values
687    u_C   = Velocity.centroid_values
688
689    for i in range(N):
690        h_C[i] = w_C[i] - z_C[i]                                               
691        if h_C[i] <= 0.01:  #epsilon :
692            #h_C[i] = 0.0
693            #w_C[i] = z_C[i]
694            uh_C[i] = 0.0
695            u_C[i] = 0.0
696           
697        if h_C[i] < epsilon:
698            h_C[i] = 0.0
699            w_C[i] = z_C[i]
700        else:
701            u_C[i]  = uh_C[i]/h_C[i] #+  h0/h_C[i])
702
703    for name in ['stage', 'height', 'xmomentum']:
704        Q = domain.quantities[name]
705        if domain.order == 1:
706            Q.extrapolate_first_order()
707        elif domain.order == 2:
708            Q.extrapolate_second_order()
709        else:
710            raise 'Unknown order'
711
712    stage_V = domain.quantities['stage'].vertex_values
713    bed_V   = domain.quantities['elevation'].vertex_values
714    h_V     = domain.quantities['height'].vertex_values
715    u_V     = domain.quantities['velocity'].vertex_values               
716    xmom_V  = domain.quantities['xmomentum'].vertex_values
717
718    bed_V[:] = stage_V - h_V
719
720    for i in range(N):                                         
721        if h_V[i] <= 0.0:
722            h_V[i] = 0.0
723            stage_V[i] = bed_V[i]
724            xmom_V[i] = 0.0
725   
726    u_V[:]  = xmom_V/(h_V + h0/h_V)   
727   
728    return
729
730
731
732
733   
734#
735def protect_against_infinitesimal_and_negative_heights(domain):
736    """Protect against infinitesimal heights and associated high velocities
737    """
738
739    #Shortcuts
740    wc = domain.quantities['stage'].centroid_values
741    zc = domain.quantities['elevation'].centroid_values
742    xmomc = domain.quantities['xmomentum'].centroid_values
743    hc = wc - zc  #Water depths at centroids
744
745    zv = domain.quantities['elevation'].vertex_values
746    wv = domain.quantities['stage'].vertex_values
747    hv = wv-zv
748    xmomv = domain.quantities['xmomentum'].vertex_values
749    #remove the above two lines and corresponding code below
750
751    #Update
752    for k in range(domain.number_of_elements):
753
754        if hc[k] < domain.minimum_allowed_height:
755            #Control stage
756            if hc[k] < domain.epsilon:
757                wc[k] = zc[k] # Contain 'lost mass' error
758                wv[k,0] = zv[k,0]
759                wv[k,1] = zv[k,1]
760               
761            xmomc[k] = 0.0
762           
763   
764def h_limiter(domain):
765    """Limit slopes for each volume to eliminate artificial variance
766    introduced by e.g. second order extrapolator
767
768    limit on h = w-z
769
770    This limiter depends on two quantities (w,z) so it resides within
771    this module rather than within quantity.py
772    """
773
774    from Numeric import Float
775    from numpy import zeros   
776
777    N = domain.number_of_elements
778    beta_h = domain.beta_h
779
780    #Shortcuts
781    wc = domain.quantities['stage'].centroid_values
782    zc = domain.quantities['elevation'].centroid_values
783    hc = wc - zc
784
785    wv = domain.quantities['stage'].vertex_values
786    zv = domain.quantities['elevation'].vertex_values
787    hv = wv-zv
788
789    hvbar = zeros(hv.shape, Float) #h-limited values
790
791    #Find min and max of this and neighbour's centroid values
792    hmax = zeros(hc.shape, Float)
793    hmin = zeros(hc.shape, Float)
794
795    for k in range(N):
796        hmax[k] = hmin[k] = hc[k]
797        for i in range(2):   
798            n = domain.neighbours[k,i]
799            if n >= 0:
800                hn = hc[n] #Neighbour's centroid value
801
802                hmin[k] = min(hmin[k], hn)
803                hmax[k] = max(hmax[k], hn)
804
805
806    #Diffences between centroids and maxima/minima
807    dhmax = hmax - hc
808    dhmin = hmin - hc
809
810    #Deltas between vertex and centroid values
811    dh = zeros(hv.shape, Float)
812    for i in range(2):
813        dh[:,i] = hv[:,i] - hc
814
815    #Phi limiter
816    for k in range(N):
817
818        #Find the gradient limiter (phi) across vertices
819        phi = 1.0
820        for i in range(2):
821            r = 1.0
822            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
823            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
824
825            phi = min( min(r*beta_h, 1), phi )
826
827        #Then update using phi limiter
828        for i in range(2):
829            hvbar[k,i] = hc[k] + phi*dh[k,i]
830
831    return hvbar
832
833def balance_deep_and_shallow(domain):
834    """Compute linear combination between stage as computed by
835    gradient-limiters limiting using w, and stage computed by
836    gradient-limiters limiting using h (h-limiter).
837    The former takes precedence when heights are large compared to the
838    bed slope while the latter takes precedence when heights are
839    relatively small.  Anything in between is computed as a balanced
840    linear combination in order to avoid numerical disturbances which
841    would otherwise appear as a result of hard switching between
842    modes.
843
844    The h-limiter is always applied irrespective of the order.
845    """
846
847    #Shortcuts
848    wc = domain.quantities['stage'].centroid_values
849    zc = domain.quantities['elevation'].centroid_values
850    hc = wc - zc
851
852    wv = domain.quantities['stage'].vertex_values
853    zv = domain.quantities['elevation'].vertex_values
854    hv = wv-zv
855
856    #Limit h
857    hvbar = h_limiter(domain)
858
859    for k in range(domain.number_of_elements):
860        #Compute maximal variation in bed elevation
861        #  This quantitiy is
862        #    dz = max_i abs(z_i - z_c)
863        #  and it is independent of dimension
864        #  In the 1d case zc = (z0+z1)/2
865        #  In the 2d case zc = (z0+z1+z2)/3
866
867        dz = max(abs(zv[k,0]-zc[k]),
868                 abs(zv[k,1]-zc[k]))#,
869#                 abs(zv[k,2]-zc[k]))
870
871
872        hmin = min( hv[k,:] )
873
874        #Create alpha in [0,1], where alpha==0 means using the h-limited
875        #stage and alpha==1 means using the w-limited stage as
876        #computed by the gradient limiter (both 1st or 2nd order)
877
878        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
879        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
880
881        if dz > 0.0:
882            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
883        else:
884            #Flat bed
885            alpha = 1.0
886
887        alpha  = 0.0
888        #Let
889        #
890        #  wvi be the w-limited stage (wvi = zvi + hvi)
891        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
892        #
893        #
894        #where i=0,1,2 denotes the vertex ids
895        #
896        #Weighted balance between w-limited and h-limited stage is
897        #
898        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
899        #
900        #It follows that the updated wvi is
901        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
902        #
903        # Momentum is balanced between constant and limited
904
905
906        #for i in range(3):
907        #    wv[k,i] = zv[k,i] + hvbar[k,i]
908
909        #return
910
911        if alpha < 1:
912
913            #for i in range(3):
914            for i in range(2):
915                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
916
917            #Momentums at centroids
918            xmomc = domain.quantities['xmomentum'].centroid_values
919
920
921            #Momentums at vertices
922            xmomv = domain.quantities['xmomentum'].vertex_values
923
924
925            # Update momentum as a linear combination of
926            # xmomc and ymomc (shallow) and momentum
927            # from extrapolator xmomv and ymomv (deep).
928            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
929
930
931
932
933
934#########################
935#Standard forcing terms:
936#
937def gravity(domain):
938    """Apply gravitational pull in the presence of bed slope
939    """
940
941    from util import gradient
942    from Numeric import Float
943    from numpy import zeros, array, sum
944
945    xmom = domain.quantities['xmomentum'].explicit_update
946    stage = domain.quantities['stage'].explicit_update
947
948    Stage = domain.quantities['stage']
949    Elevation = domain.quantities['elevation']
950    h = Stage.vertex_values - Elevation.vertex_values
951    b = Elevation.vertex_values
952    w = Stage.vertex_values
953
954    x = domain.get_vertex_coordinates()
955    g = domain.g
956
957    for k in range(domain.number_of_elements):
958        avg_h = sum( h[k,:] )/2
959
960        #Compute bed slope
961        x0, x1 = x[k,:]
962        b0, b1 = b[k,:]
963
964        w0, w1 = w[k,:]
965        wx = gradient(x0, x1, w0, w1)
966        bx = gradient(x0, x1, b0, b1)
967       
968        #Update momentum (explicit update is reset to source values)
969        xmom[k] += -g*bx*avg_h
970
971
972def gravity_F1(domain):
973    """Apply gravitational pull in the presence of bed slope
974    """
975
976    from util import gradient
977    from Numeric import Float
978    from numpy import zeros, array, sum   
979    from parameters import F1#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
980
981    xmom = domain.quantities['xmomentum'].explicit_update
982    stage = domain.quantities['stage'].explicit_update
983
984    Stage = domain.quantities['stage']
985    Elevation = domain.quantities['elevation']
986    h = Stage.vertex_values - Elevation.vertex_values
987    b = Elevation.vertex_values
988    w = Stage.vertex_values
989
990    x = domain.get_vertex_coordinates()
991    g = domain.g
992
993    for k in range(domain.number_of_elements):
994        avg_h = sum( h[k,:] )/2
995
996        #Compute bed slope
997        x0, x1 = x[k,:]
998        b0, b1 = b[k,:]
999
1000        w0, w1 = w[k,:]
1001        wx = gradient(x0, x1, w0, w1)
1002        bx = gradient(x0, x1, b0, b1)
1003       
1004        #Update momentum (explicit update is reset to source values)
1005        xmom[k] += -g*bx*avg_h + avg_h*F1#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1006
1007
1008
1009def gravity_F2(domain):
1010    """Apply gravitational pull in the presence of bed slope
1011    """
1012
1013    from util import gradient
1014    from Numeric import Float
1015    from numpy import zeros, array, sum   
1016    from parameters import F2#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1017
1018    xmom = domain.quantities['xmomentum'].explicit_update
1019    stage = domain.quantities['stage'].explicit_update
1020
1021    Stage = domain.quantities['stage']
1022    Elevation = domain.quantities['elevation']
1023    h = Stage.vertex_values - Elevation.vertex_values
1024    b = Elevation.vertex_values
1025    w = Stage.vertex_values
1026
1027    x = domain.get_vertex_coordinates()
1028    g = domain.g
1029
1030    for k in range(domain.number_of_elements):
1031        avg_h = sum( h[k,:] )/2
1032
1033        #Compute bed slope
1034        x0, x1 = x[k,:]
1035        b0, b1 = b[k,:]
1036
1037        w0, w1 = w[k,:]
1038        wx = gradient(x0, x1, w0, w1)
1039        bx = gradient(x0, x1, b0, b1)
1040       
1041        #Update momentum (explicit update is reset to source values)
1042        xmom[k] += -g*bx*avg_h + avg_h*F2#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1043
1044
1045
1046 
1047 
1048def manning_friction(domain):
1049    """Apply (Manning) friction to water momentum
1050    """
1051
1052    from math import sqrt
1053
1054    w = domain.quantities['stage'].centroid_values
1055    z = domain.quantities['elevation'].centroid_values
1056    h = w-z
1057
1058    uh = domain.quantities['xmomentum'].centroid_values
1059    eta = domain.quantities['friction'].centroid_values
1060
1061    xmom_update = domain.quantities['xmomentum'].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 * uh[k]
1071                S /= h[k]**(7.0/3)
1072
1073                #Update momentum
1074                xmom_update[k] += S*uh[k]
1075
1076def linear_friction(domain):
1077    """Apply linear friction to water momentum
1078
1079    Assumes quantity: 'linear_friction' to be present
1080    """
1081
1082    from math import sqrt
1083
1084    w = domain.quantities['stage'].centroid_values
1085    z = domain.quantities['elevation'].centroid_values
1086    h = w-z
1087
1088    uh = domain.quantities['xmomentum'].centroid_values
1089    tau = domain.quantities['linear_friction'].centroid_values
1090
1091    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1092
1093    N = domain.number_of_elements
1094    eps = domain.minimum_allowed_height
1095    g = domain.g #Not necessary? Why was this added?
1096
1097    for k in range(N):
1098        if tau[k] >= eps:
1099            if h[k] >= eps:
1100                S = -tau[k]/h[k]
1101
1102                #Update momentum
1103                xmom_update[k] += S*uh[k]
1104
1105
1106
1107def check_forcefield(f):
1108    """Check that f is either
1109    1: a callable object f(t,x,y), where x and y are vectors
1110       and that it returns an array or a list of same length
1111       as x and y
1112    2: a scalar
1113    """
1114
1115    from Numeric import Float
1116    from numpy import ones, array   
1117
1118
1119    if callable(f):
1120        N = 2
1121        x = ones(2, Float)
1122       
1123        try:
1124            q = f(1.0, x=x)
1125        except Exception, e:
1126            msg = 'Function %s could not be executed:\n%s' %(f, e)
1127            #FIXME: Reconsider this semantics
1128            raise msg
1129
1130        try:
1131            q = array(q).astype(Float)
1132        except:
1133            msg = 'Return value from vector function %s could ' %f
1134            msg += 'not be converted into a Numeric array of floats.\n'
1135            msg += 'Specified function should return either list or array.'
1136            raise msg
1137
1138        #Is this really what we want?
1139        msg = 'Return vector from function %s ' %f
1140        msg += 'must have same lenght as input vectors'
1141        assert len(q) == N, msg
1142
1143    else:
1144        try:
1145            f = float(f)
1146        except:
1147            msg = 'Force field %s must be either a scalar' %f
1148            msg += ' or a vector function'
1149            raise msg
1150    return f
1151
1152class Wind_stress:
1153    """Apply wind stress to water momentum in terms of
1154    wind speed [m/s] and wind direction [degrees]
1155    """
1156
1157    def __init__(self, *args, **kwargs):
1158        """Initialise windfield from wind speed s [m/s]
1159        and wind direction phi [degrees]
1160
1161        Inputs v and phi can be either scalars or Python functions, e.g.
1162
1163        W = Wind_stress(10, 178)
1164
1165        #FIXME - 'normal' degrees are assumed for now, i.e. the
1166        vector (1,0) has zero degrees.
1167        We may need to convert from 'compass' degrees later on and also
1168        map from True north to grid north.
1169
1170        Arguments can also be Python functions of t,x,y as in
1171
1172        def speed(t,x,y):
1173            ...
1174            return s
1175
1176        def angle(t,x,y):
1177            ...
1178            return phi
1179
1180        where x and y are vectors.
1181
1182        and then pass the functions in
1183
1184        W = Wind_stress(speed, angle)
1185
1186        The instantiated object W can be appended to the list of
1187        forcing_terms as in
1188
1189        Alternatively, one vector valued function for (speed, angle)
1190        can be applied, providing both quantities simultaneously.
1191        As in
1192        W = Wind_stress(F), where returns (speed, angle) for each t.
1193
1194        domain.forcing_terms.append(W)
1195        """
1196
1197        from config import rho_a, rho_w, eta_w
1198        from Numeric import Float
1199        from numpy import array       
1200
1201        if len(args) == 2:
1202            s = args[0]
1203            phi = args[1]
1204        elif len(args) == 1:
1205            #Assume vector function returning (s, phi)(t,x,y)
1206            vector_function = args[0]
1207            s = lambda t,x: vector_function(t,x=x)[0]
1208            phi = lambda t,x: vector_function(t,x=x)[1]
1209        else:
1210           #Assume info is in 2 keyword arguments
1211
1212           if len(kwargs) == 2:
1213               s = kwargs['s']
1214               phi = kwargs['phi']
1215           else:
1216               raise 'Assumes two keyword arguments: s=..., phi=....'
1217
1218        print 'phi', phi
1219        self.speed = check_forcefield(s)
1220        self.phi = check_forcefield(phi)
1221
1222        self.const = eta_w*rho_a/rho_w
1223
1224
1225    def __call__(self, domain):
1226        """Evaluate windfield based on values found in domain
1227        """
1228
1229        from math import pi, cos, sin, sqrt
1230        from Numeric import Float, ArrayType
1231        from numpy import ones
1232
1233        xmom_update = domain.quantities['xmomentum'].explicit_update
1234
1235        N = domain.number_of_elements
1236        t = domain.time
1237
1238        if callable(self.speed):
1239            xc = domain.get_centroid_coordinates()
1240            s_vec = self.speed(t, xc)
1241        else:
1242            #Assume s is a scalar
1243
1244            try:
1245                s_vec = self.speed * ones(N, Float)
1246            except:
1247                msg = 'Speed must be either callable or a scalar: %s' %self.s
1248                raise msg
1249
1250
1251        if callable(self.phi):
1252            xc = domain.get_centroid_coordinates()
1253            phi_vec = self.phi(t, xc)
1254        else:
1255            #Assume phi is a scalar
1256
1257            try:
1258                phi_vec = self.phi * ones(N, Float)
1259            except:
1260                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1261                raise msg
1262
1263        #assign_windfield_values(xmom_update, ymom_update,
1264        #                        s_vec, phi_vec, self.const)
1265        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
1266
1267
1268#def assign_windfield_values(xmom_update, ymom_update,
1269#                            s_vec, phi_vec, const):
1270def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
1271    """Python version of assigning wind field to update vectors.
1272    A c version also exists (for speed)
1273    """
1274    from math import pi, cos, sin, sqrt
1275
1276    N = len(s_vec)
1277    for k in range(N):
1278        s = s_vec[k]
1279        phi = phi_vec[k]
1280
1281        #Convert to radians
1282        phi = phi*pi/180
1283
1284        #Compute velocity vector (u, v)
1285        u = s*cos(phi)
1286        v = s*sin(phi)
1287
1288        #Compute wind stress
1289        S = const * u
1290        xmom_update[k] += S*u
Note: See TracBrowser for help on using the repository browser.