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

Last change on this file since 7837 was 7837, checked in by mungkasi, 12 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

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