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

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

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