source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/cg/domain_cg.py @ 7922

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

Modifying codes so that arrays are compatible with numpy instead of Numeric.

File size: 30.7 KB
Line 
1"""Class Domain -
21D interval domains for finite-volume computations of
3the shallow water wave equation.
4
5This module contains a specialisation of class Domain from module domain.py
6consisting of methods specific to the Shallow Water Wave Equation
7
8
9U_t + E_x = S
10
11where
12
13U = [w, uh]
14E = [uh, u^2h + gh^2/2]
15S represents source terms forcing the system
16(e.g. gravity, friction, wind stress, ...)
17
18and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
19
20The quantities are
21
22symbol    variable name    explanation
23x         x                horizontal distance from origin [m]
24z         elevation        elevation of bed on which flow is modelled [m]
25h         height           water height above z [m]
26w         stage            absolute water level, w = z+h [m]
27u                          speed in the x direction [m/s]
28uh        xmomentum        momentum in the x direction [m^2/s]
29
30eta                        mannings friction coefficient [to appear]
31nu                         wind stress coefficient [to appear]
32
33The conserved quantities are w, uh
34
35For details see e.g.
36Christopher Zoppou and Stephen Roberts,
37Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
38Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
39
40John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
41Geoscience Australia, 2006
42"""
43
44
45from domain import *
46Generic_Domain = Domain #Rename
47
48#Shallow water domain
49class Domain(Generic_Domain):
50
51    def __init__(self, coordinates, boundary = None, tagged_elements = None):
52        conserved_quantities = ['stage', 'xmomentum'] #['height', 'xmomentum']
53        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
54        other_quantities = ['friction']
55        Generic_Domain.__init__(self, coordinates, boundary,
56                                conserved_quantities, evolved_quantities, other_quantities,
57                                tagged_elements)
58       
59        from config import minimum_allowed_height, g, h0
60        self.minimum_allowed_height = minimum_allowed_height
61        self.g = g
62        self.h0 = h0
63
64        self.forcing_terms.append(gravity)
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__ = 'domain_cg'
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
166        #(either from this module or C-extension)
167        compute_fluxes_C_wellbalanced(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
175        #(either from this module or C-extension)
176        distribute_to_vertices_and_edges(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 
478
479# ###################################
480# Module functions for gradient limiting (distribute_to_vertices_and_edges)
481
482def distribute_to_vertices_and_edges(domain):
483    """Distribution from centroids to vertices specific to the
484    shallow water wave
485    equation.
486
487    It will ensure that h (w-z) is always non-negative even in the
488    presence of steep bed-slopes by taking a weighted average between shallow
489    and deep cases.
490
491    In addition, all conserved quantities get distributed as per either a
492    constant (order==1) or a piecewise linear function (order==2).
493
494    FIXME: more explanation about removal of artificial variability etc
495
496    Precondition:
497      All quantities defined at centroids and bed elevation defined at
498      vertices.
499
500    Postcondition
501      Conserved quantities defined at vertices
502
503    """
504
505    #from config import optimised_gradient_limiter
506
507    #Remove very thin layers of water
508    #protect_against_infinitesimal_and_negative_heights(domain) 
509
510    import sys
511    from Numeric import Float
512    from numpy import zeros
513    from config import epsilon, h0, h_min
514
515    N = domain.number_of_elements
516
517    #Shortcuts
518    Stage = domain.quantities['stage']
519    Xmom = domain.quantities['xmomentum']
520    Bed = domain.quantities['elevation']
521    Height = domain.quantities['height']
522    Velocity = domain.quantities['velocity']
523
524    #Arrays   
525    w_C   = Stage.centroid_values   
526    uh_C  = Xmom.centroid_values   
527    z_C   = Bed.centroid_values
528    h_C   = Height.centroid_values
529    u_C   = Velocity.centroid_values
530
531    for i in range(N):
532        h_C[i] = w_C[i] - z_C[i]                                               
533        if h_C[i] <= 0.01:  #epsilon :
534            #h_C[i] = 0.0
535            #w_C[i] = z_C[i]
536            uh_C[i] = 0.0
537            u_C[i] = 0.0
538           
539        if h_C[i] < epsilon: #h_min:
540            h_C[i] = 0.0
541            w_C[i] = z_C[i]
542        else:
543            u_C[i]  = uh_C[i]/h_C[i] #+  h0/h_C[i])
544
545    for name in ['stage', 'height', 'velocity']:
546        Q = domain.quantities[name]
547        if domain.order == 1:
548            Q.extrapolate_first_order()
549        elif domain.order == 2:
550            Q.extrapolate_second_order()
551        else:
552            raise 'Unknown order'
553
554    stage_V = domain.quantities['stage'].vertex_values
555    bed_V   = domain.quantities['elevation'].vertex_values
556    h_V     = domain.quantities['height'].vertex_values
557    u_V     = domain.quantities['velocity'].vertex_values               
558    xmom_V  = domain.quantities['xmomentum'].vertex_values
559
560    bed_V[:] = stage_V - h_V
561    xmom_V[:] = u_V * h_V
562   
563    return
564
565   
566#
567def protect_against_infinitesimal_and_negative_heights(domain):
568    """Protect against infinitesimal heights and associated high velocities
569    """
570
571    #Shortcuts
572    wc = domain.quantities['stage'].centroid_values
573    zc = domain.quantities['elevation'].centroid_values
574    xmomc = domain.quantities['xmomentum'].centroid_values
575    hc = wc - zc  #Water depths at centroids
576
577    zv = domain.quantities['elevation'].vertex_values
578    wv = domain.quantities['stage'].vertex_values
579    hv = wv-zv
580    xmomv = domain.quantities['xmomentum'].vertex_values
581    #remove the above two lines and corresponding code below
582
583    #Update
584    for k in range(domain.number_of_elements):
585
586        if hc[k] < domain.minimum_allowed_height:
587            #Control stage
588            if hc[k] < domain.epsilon:
589                wc[k] = zc[k] # Contain 'lost mass' error
590                wv[k,0] = zv[k,0]
591                wv[k,1] = zv[k,1]
592               
593            xmomc[k] = 0.0
594           
595   
596def h_limiter(domain):
597    """Limit slopes for each volume to eliminate artificial variance
598    introduced by e.g. second order extrapolator
599
600    limit on h = w-z
601
602    This limiter depends on two quantities (w,z) so it resides within
603    this module rather than within quantity.py
604    """
605
606    from Numeric import Float
607    from numpy import zeros
608
609    N = domain.number_of_elements
610    beta_h = domain.beta_h
611
612    #Shortcuts
613    wc = domain.quantities['stage'].centroid_values
614    zc = domain.quantities['elevation'].centroid_values
615    hc = wc - zc
616
617    wv = domain.quantities['stage'].vertex_values
618    zv = domain.quantities['elevation'].vertex_values
619    hv = wv-zv
620
621    hvbar = zeros(hv.shape, Float) #h-limited values
622
623    #Find min and max of this and neighbour's centroid values
624    hmax = zeros(hc.shape, Float)
625    hmin = zeros(hc.shape, Float)
626
627    for k in range(N):
628        hmax[k] = hmin[k] = hc[k]
629        for i in range(2):   
630            n = domain.neighbours[k,i]
631            if n >= 0:
632                hn = hc[n] #Neighbour's centroid value
633
634                hmin[k] = min(hmin[k], hn)
635                hmax[k] = max(hmax[k], hn)
636
637
638    #Diffences between centroids and maxima/minima
639    dhmax = hmax - hc
640    dhmin = hmin - hc
641
642    #Deltas between vertex and centroid values
643    dh = zeros(hv.shape, Float)
644    for i in range(2):
645        dh[:,i] = hv[:,i] - hc
646
647    #Phi limiter
648    for k in range(N):
649
650        #Find the gradient limiter (phi) across vertices
651        phi = 1.0
652        for i in range(2):
653            r = 1.0
654            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
655            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
656
657            phi = min( min(r*beta_h, 1), phi )
658
659        #Then update using phi limiter
660        for i in range(2):
661            hvbar[k,i] = hc[k] + phi*dh[k,i]
662
663    return hvbar
664
665def balance_deep_and_shallow(domain):
666    """Compute linear combination between stage as computed by
667    gradient-limiters limiting using w, and stage computed by
668    gradient-limiters limiting using h (h-limiter).
669    The former takes precedence when heights are large compared to the
670    bed slope while the latter takes precedence when heights are
671    relatively small.  Anything in between is computed as a balanced
672    linear combination in order to avoid numerical disturbances which
673    would otherwise appear as a result of hard switching between
674    modes.
675
676    The h-limiter is always applied irrespective of the order.
677    """
678
679    #Shortcuts
680    wc = domain.quantities['stage'].centroid_values
681    zc = domain.quantities['elevation'].centroid_values
682    hc = wc - zc
683
684    wv = domain.quantities['stage'].vertex_values
685    zv = domain.quantities['elevation'].vertex_values
686    hv = wv-zv
687
688    #Limit h
689    hvbar = h_limiter(domain)
690
691    for k in range(domain.number_of_elements):
692        #Compute maximal variation in bed elevation
693        #  This quantitiy is
694        #    dz = max_i abs(z_i - z_c)
695        #  and it is independent of dimension
696        #  In the 1d case zc = (z0+z1)/2
697        #  In the 2d case zc = (z0+z1+z2)/3
698
699        dz = max(abs(zv[k,0]-zc[k]),
700                 abs(zv[k,1]-zc[k]))#,
701#                 abs(zv[k,2]-zc[k]))
702
703
704        hmin = min( hv[k,:] )
705
706        #Create alpha in [0,1], where alpha==0 means using the h-limited
707        #stage and alpha==1 means using the w-limited stage as
708        #computed by the gradient limiter (both 1st or 2nd order)
709
710        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
711        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
712
713        if dz > 0.0:
714            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
715        else:
716            #Flat bed
717            alpha = 1.0
718
719        alpha  = 0.0
720        #Let
721        #
722        #  wvi be the w-limited stage (wvi = zvi + hvi)
723        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
724        #
725        #
726        #where i=0,1,2 denotes the vertex ids
727        #
728        #Weighted balance between w-limited and h-limited stage is
729        #
730        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
731        #
732        #It follows that the updated wvi is
733        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
734        #
735        # Momentum is balanced between constant and limited
736
737
738        #for i in range(3):
739        #    wv[k,i] = zv[k,i] + hvbar[k,i]
740
741        #return
742
743        if alpha < 1:
744
745            #for i in range(3):
746            for i in range(2):
747                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
748
749            #Momentums at centroids
750            xmomc = domain.quantities['xmomentum'].centroid_values
751
752
753            #Momentums at vertices
754            xmomv = domain.quantities['xmomentum'].vertex_values
755
756
757            # Update momentum as a linear combination of
758            # xmomc and ymomc (shallow) and momentum
759            # from extrapolator xmomv and ymomv (deep).
760            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
761
762
763
764
765
766#########################
767#Standard forcing terms:
768#
769def gravity(domain):
770    """Apply gravitational pull in the presence of bed slope
771    """
772
773    from util import gradient
774    from Numeric import Float
775    from numpy import zeros, array, sum
776
777    xmom = domain.quantities['xmomentum'].explicit_update
778    stage = domain.quantities['stage'].explicit_update
779
780    Stage = domain.quantities['stage']
781    Elevation = domain.quantities['elevation']
782    h = Stage.vertex_values - Elevation.vertex_values
783    b = Elevation.vertex_values
784    w = Stage.vertex_values
785
786    x = domain.get_vertex_coordinates()
787    g = domain.g
788
789    for k in range(domain.number_of_elements):
790        avg_h = sum( h[k,:] )/2
791
792        #Compute bed slope
793        x0, x1 = x[k,:]
794        b0, b1 = b[k,:]
795
796        w0, w1 = w[k,:]
797        wx = gradient(x0, x1, w0, w1)
798        bx = gradient(x0, x1, b0, b1)
799       
800        #Update momentum (explicit update is reset to source values)
801        xmom[k] += -g*bx*avg_h
802
803 
804 
805def manning_friction(domain):
806    """Apply (Manning) friction to water momentum
807    """
808
809    from math import sqrt
810
811    w = domain.quantities['stage'].centroid_values
812    z = domain.quantities['elevation'].centroid_values
813    h = w-z
814
815    uh = domain.quantities['xmomentum'].centroid_values
816    eta = domain.quantities['friction'].centroid_values
817
818    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
819
820    N = domain.number_of_elements
821    eps = domain.minimum_allowed_height
822    g = domain.g
823
824    for k in range(N):
825        if eta[k] >= eps:
826            if h[k] >= eps:
827                S = -g * eta[k]**2 * uh[k]
828                S /= h[k]**(7.0/3)
829
830                #Update momentum
831                xmom_update[k] += S*uh[k]
832
833def linear_friction(domain):
834    """Apply linear friction to water momentum
835
836    Assumes quantity: 'linear_friction' to be present
837    """
838
839    from math import sqrt
840
841    w = domain.quantities['stage'].centroid_values
842    z = domain.quantities['elevation'].centroid_values
843    h = w-z
844
845    uh = domain.quantities['xmomentum'].centroid_values
846    tau = domain.quantities['linear_friction'].centroid_values
847
848    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
849
850    N = domain.number_of_elements
851    eps = domain.minimum_allowed_height
852    g = domain.g #Not necessary? Why was this added?
853
854    for k in range(N):
855        if tau[k] >= eps:
856            if h[k] >= eps:
857                S = -tau[k]/h[k]
858
859                #Update momentum
860                xmom_update[k] += S*uh[k]
861 #              ymom_update[k] += S*vh[k]
862
863
864
865def check_forcefield(f):
866    """Check that f is either
867    1: a callable object f(t,x,y), where x and y are vectors
868       and that it returns an array or a list of same length
869       as x and y
870    2: a scalar
871    """
872
873    from Numeric import Float
874    from numpy import ones, array   
875
876
877    if callable(f):
878        N = 2
879        x = ones(2, Float)
880       
881        try:
882            q = f(1.0, x=x)
883        except Exception, e:
884            msg = 'Function %s could not be executed:\n%s' %(f, e)
885            #FIXME: Reconsider this semantics
886            raise msg
887
888        try:
889            q = array(q).astype(Float)
890        except:
891            msg = 'Return value from vector function %s could ' %f
892            msg += 'not be converted into a Numeric array of floats.\n'
893            msg += 'Specified function should return either list or array.'
894            raise msg
895
896        #Is this really what we want?
897        msg = 'Return vector from function %s ' %f
898        msg += 'must have same lenght as input vectors'
899        assert len(q) == N, msg
900
901    else:
902        try:
903            f = float(f)
904        except:
905            msg = 'Force field %s must be either a scalar' %f
906            msg += ' or a vector function'
907            raise msg
908    return f
909
910class Wind_stress:
911    """Apply wind stress to water momentum in terms of
912    wind speed [m/s] and wind direction [degrees]
913    """
914
915    def __init__(self, *args, **kwargs):
916        """Initialise windfield from wind speed s [m/s]
917        and wind direction phi [degrees]
918
919        Inputs v and phi can be either scalars or Python functions, e.g.
920
921        W = Wind_stress(10, 178)
922
923        #FIXME - 'normal' degrees are assumed for now, i.e. the
924        vector (1,0) has zero degrees.
925        We may need to convert from 'compass' degrees later on and also
926        map from True north to grid north.
927
928        Arguments can also be Python functions of t,x,y as in
929
930        def speed(t,x,y):
931            ...
932            return s
933
934        def angle(t,x,y):
935            ...
936            return phi
937
938        where x and y are vectors.
939
940        and then pass the functions in
941
942        W = Wind_stress(speed, angle)
943
944        The instantiated object W can be appended to the list of
945        forcing_terms as in
946
947        Alternatively, one vector valued function for (speed, angle)
948        can be applied, providing both quantities simultaneously.
949        As in
950        W = Wind_stress(F), where returns (speed, angle) for each t.
951
952        domain.forcing_terms.append(W)
953        """
954
955        from config import rho_a, rho_w, eta_w
956        from Numeric import Float
957        from numpy import array
958
959        if len(args) == 2:
960            s = args[0]
961            phi = args[1]
962        elif len(args) == 1:
963            #Assume vector function returning (s, phi)(t,x,y)
964            vector_function = args[0]
965            s = lambda t,x: vector_function(t,x=x)[0]
966            phi = lambda t,x: vector_function(t,x=x)[1]
967        else:
968           #Assume info is in 2 keyword arguments
969
970           if len(kwargs) == 2:
971               s = kwargs['s']
972               phi = kwargs['phi']
973           else:
974               raise 'Assumes two keyword arguments: s=..., phi=....'
975
976        print 'phi', phi
977        self.speed = check_forcefield(s)
978        self.phi = check_forcefield(phi)
979
980        self.const = eta_w*rho_a/rho_w
981
982
983    def __call__(self, domain):
984        """Evaluate windfield based on values found in domain
985        """
986
987        from math import pi, cos, sin, sqrt
988        from Numeric import Float
989        from numpy import ones, ArrayType       
990
991        xmom_update = domain.quantities['xmomentum'].explicit_update
992
993        N = domain.number_of_elements
994        t = domain.time
995
996        if callable(self.speed):
997            xc = domain.get_centroid_coordinates()
998            s_vec = self.speed(t, xc)
999        else:
1000            #Assume s is a scalar
1001
1002            try:
1003                s_vec = self.speed * ones(N, Float)
1004            except:
1005                msg = 'Speed must be either callable or a scalar: %s' %self.s
1006                raise msg
1007
1008
1009        if callable(self.phi):
1010            xc = domain.get_centroid_coordinates()
1011            phi_vec = self.phi(t, xc)
1012        else:
1013            #Assume phi is a scalar
1014
1015            try:
1016                phi_vec = self.phi * ones(N, Float)
1017            except:
1018                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1019                raise msg
1020
1021        #assign_windfield_values(xmom_update, ymom_update,
1022        #                        s_vec, phi_vec, self.const)
1023        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
1024
1025
1026#def assign_windfield_values(xmom_update, ymom_update,
1027#                            s_vec, phi_vec, const):
1028def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
1029    """Python version of assigning wind field to update vectors.
1030    A c version also exists (for speed)
1031    """
1032    from math import pi, cos, sin, sqrt
1033
1034    N = len(s_vec)
1035    for k in range(N):
1036        s = s_vec[k]
1037        phi = phi_vec[k]
1038
1039        #Convert to radians
1040        phi = phi*pi/180
1041
1042        #Compute velocity vector (u, v)
1043        u = s*cos(phi)
1044        v = s*sin(phi)
1045
1046        #Compute wind stress
1047        S = const * u
1048        xmom_update[k] += S*u
Note: See TracBrowser for help on using the repository browser.