source: anuga_work/development/sudi/sw_1d/periodic_waves/johns/domain_johns.py @ 7837

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

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

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