source: anuga_work/development/shallow_water_1d_old/shallow_water_domain.py @ 5533

Last change on this file since 5533 was 4959, checked in by steve, 16 years ago

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

File size: 27.7 KB
Line 
1"""Class Domain -
22D triangular domains for finite-volume computations of
3the shallow water wave equation.
4
5
6$Description:
7This module contains a specialisation of class Domain from module domain.py
8consisting of methods specific to the Shallow Water Wave Equation
9
10
11U_t + E_x + G_y = S
12
13where
14
15U = [w, uh, vh]
16E = [uh, u^2h + gh^2/2, uvh]
17G = [vh, uvh, v^2h + gh^2/2]
18S represents source terms forcing the system
19(e.g. gravity, friction, wind stress, ...)
20
21and _t, _x, _y denote the derivative with respect to t, x and y respectively.
22
23The quantities are
24
25symbol    variable name    explanation
26x         x                horizontal distance from origin [m]
27y         y                vertical distance from origin [m]
28z         elevation        elevation of bed on which flow is modelled [m]
29h         height           water height above z [m]
30w         stage            absolute water level, w = z+h [m]
31u                          speed in the x direction [m/s]
32v                          speed in the y direction [m/s]
33uh        xmomentum        momentum in the x direction [m^2/s]
34vh        ymomentum        momentum in the y direction [m^2/s]
35
36eta                        mannings friction coefficient [to appear]
37nu                         wind stress coefficient [to appear]
38
39The conserved quantities are w, uh, vh
40
41$References
42Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
43Christopher Zoppou and Stephen Roberts,
44Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
45
46Hydrodynamic modelling of coastal inundation.
47Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
48In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
49Modelling and Simulation. Modelling and Simulation Society of Australia and
50New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.
51http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
52
53
54Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
55Geoscience Australia, 2004
56
57
58$Author: ole $
59$Revision: 3642 $
60$Date: 2006-09-21 11:30:59 +1000 (Thu, 21 Sep 2006) $
61$LastChangedDate: 2006-09-21 11:30:59 +1000 (Thu, 21 Sep 2006) $
62$LastChangedRevision: 3642 $
63$LastChangedBy: ole $
64$HeadURL: https://datamining/svn/ga/anuga_core/source/anuga/shallow_water/shallow_water_domain.py $
65"""
66
67#Subversion keywords:
68#
69#$LastChangedDate: 2006-09-21 11:30:59 +1000 (Thu, 21 Sep 2006) $
70#$LastChangedRevision: 3642 $
71#$LastChangedBy: ole $
72
73
74from domain_t2 import *
75Generic_Domain = Domain #Rename
76
77#Shallow water domain
78class Domain(Generic_Domain):
79
80    def __init__(self,
81                 coordinates=None,
82                 vertices=None,
83                 boundary=None,
84                 tagged_elements=None,
85                 geo_reference=None):
86
87
88        conserved_quantities = ['stage', 'xmomentum']
89        other_quantities = ['elevation', 'friction']
90        Generic_Domain.__init__(self,
91                                coordinates,
92                                boundary,
93                                conserved_quantities,
94                                other_quantities,
95                                tagged_elements,
96                                geo_reference)
97
98
99        from config import minimum_allowed_height, maximum_allowed_speed, g
100        self.minimum_allowed_height = minimum_allowed_height
101        self.maximum_allowed_speed = maximum_allowed_speed
102        self.g = g
103
104
105        #self.forcing_terms.append(manning_friction)
106        self.forcing_terms.append(gravity)
107
108    def set_quantities_to_be_stored(self, q):
109        """Specify which quantities will be stored in the sww file.
110
111        q must be either:
112          - the name of a quantity
113          - a list of quantity names
114          - None
115
116        In the two first cases, the named quantities will be stored at each yieldstep
117        (This is in addition to the quantities elevation and friction)
118        If q is None, storage will be switched off altogether.
119        """
120
121
122        if q is None:
123            self.quantities_to_be_stored = []
124            self.store = False
125            return
126
127        if isinstance(q, basestring):
128            q = [q] # Turn argument into a list
129
130        #Check correcness
131        for quantity_name in q:
132            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
133            assert quantity_name in self.conserved_quantities, msg
134
135        self.quantities_to_be_stored = q
136
137
138    def check_integrity(self):
139        Generic_Domain.check_integrity(self)
140
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    def extrapolate_second_order_sw(self):
149        #Call correct module function
150        #(either from this module or C-extension)
151        extrapolate_second_order_sw(self)
152
153    def compute_timestep(self):
154        #Call correct module function
155        compute_timestep(self)
156
157    def compute_fluxes(self):
158        #Call correct module function
159        #(either from this module or C-extension)
160        compute_fluxes(self)
161
162    def distribute_to_vertices_and_edges(self):
163        #Call correct module function
164        #(either from this module or C-extension)
165        distribute_to_vertices_and_edges(self)
166
167    def evolve(self,
168               yieldstep = None,
169               finaltime = None,
170               duration = None,
171               skip_initial_step = False):
172        """Specialisation of basic evolve method from parent class
173        """
174
175        #Call check integrity here rather than from user scripts
176        #self.check_integrity()
177
178        msg = 'Parameter beta_h must be in the interval [0, 1['
179        assert 0 <= self.beta_h < 1.0, msg
180        msg = 'Parameter beta_w must be in the interval [0, 1['
181        assert 0 <= self.beta_w < 1.0, msg
182
183
184        #Initial update of vertex and edge values before any storage
185        #and or visualisation
186        self.distribute_to_vertices_and_edges()
187
188        #Call basic machinery from parent class
189        for t in Generic_Domain.evolve(self,
190                                       yieldstep=yieldstep,
191                                       finaltime=finaltime,
192                                       skip_initial_step=skip_initial_step):
193
194
195            #Pass control on to outer loop for more specific actions
196            yield(t)
197
198
199#=============== End of Shallow Water Domain ===============================
200
201
202####################################
203# Flux computation
204def flux_function(normal, ql, qr, zl, zr):
205    """Compute fluxes between volumes for the shallow water wave equation
206    cast in terms of w = h+z using the 'central scheme' as described in
207
208    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
209    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
210    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
211
212    The implemented formula is given in equation (3.15) on page 714
213
214    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
215    in the numerical vectors ql and qr.
216
217    Bed elevations zl and zr.
218    """
219
220    from config import g, epsilon
221    from math import sqrt
222    from Numeric import array, zeros, Float
223
224    #Align momentums with x-axis
225    #q_left  = rotate(ql, normal, direction = 1)
226    #q_right = rotate(qr, normal, direction = 1)
227    q_left = zeros(len(ql),Float)
228    q_right = zeros(len(qr),Float)
229   
230    q_left[:] = ql[:]
231    q_right[:] = qr[:]
232   
233    q_left[1] = q_left[1]*normal
234    q_right = qr
235    q_right[1] = q_right[1]*normal
236
237    z = (zl+zr)/2.0 #Take average of field values
238
239    w_left  = q_left[0]   #w=h+z
240    h_left  = w_left-z
241    uh_left = q_left[1]
242
243    if h_left < epsilon:
244        u_left = 0.0  #Could have been negative
245        h_left = 0.0
246    else:
247        u_left  = uh_left/h_left
248
249
250    w_right  = q_right[0]  #w=h+z
251    h_right  = w_right-z
252    uh_right = q_right[1]
253
254
255    if h_right < epsilon:
256        u_right = 0.0  #Could have been negative
257        h_right = 0.0
258    else:
259        u_right  = uh_right/h_right
260
261    soundspeed_left  = sqrt(g*h_left)
262    soundspeed_right = sqrt(g*h_right)
263
264    #Maximal wave speed
265    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
266
267    #Minimal wave speed
268    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
269
270    #Flux computation
271
272    #FIXME(Ole): Why is it again that we don't
273    #use uh_left and uh_right directly in the first entries?
274    flux_left  = array([u_left*h_left,
275                        u_left*uh_left + 0.5*g*h_left**2])
276    flux_right = array([u_right*h_right,
277                        u_right*uh_right + 0.5*g*h_right**2])
278
279    denom = s_max-s_min
280    if denom == 0.0:
281        edgeflux = array([0.0, 0.0])
282        max_speed = 0.0
283    else:
284        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
285        edgeflux += s_max*s_min*(q_right-q_left)/denom
286
287        #edgeflux = rotate(edgeflux, normal, direction=-1)
288        edgeflux[1] = edgeflux[1]*normal
289        max_speed = max(abs(s_max), abs(s_min))
290
291    return edgeflux, max_speed
292
293def compute_timestep(domain):
294    import sys
295    from Numeric import zeros, Float
296
297    N = domain.number_of_elements
298
299    #Shortcuts
300    Stage = domain.quantities['stage']
301    Xmom = domain.quantities['xmomentum']
302    Bed = domain.quantities['elevation']
303   
304    stage = Stage.vertex_values
305    xmom =  Xmom.vertex_values
306    bed =   Bed.vertex_values
307
308    stage_bdry = Stage.boundary_values
309    xmom_bdry =  Xmom.boundary_values
310
311    flux = zeros(2, Float) #Work array for summing up fluxes
312    ql = zeros(2, Float)
313    qr = zeros(2, Float)
314
315    #Loop
316    timestep = float(sys.maxint)
317    enter = True
318    for k in range(N):
319
320        flux[:] = 0.  #Reset work array
321        for i in range(2):
322            #Quantities inside volume facing neighbour i
323            ql = [stage[k, i], xmom[k, i]]
324            zl = bed[k, i]
325
326            #Quantities at neighbour on nearest face
327            n = domain.neighbours[k,i]
328            if n < 0:
329                m = -n-1 #Convert negative flag to index
330                qr[0] = stage_bdry[m]
331                qr[1] = xmom_bdry[m]
332                zr = zl #Extend bed elevation to boundary
333            else:
334                #m = domain.neighbour_edges[k,i]
335                m = domain.neighbour_vertices[k,i]
336                qr[0] = stage[n, m]
337                qr[1] =  xmom[n, m]
338                zr = bed[n, m]
339
340
341            #Outward pointing normal vector
342            normal = domain.normals[k, i]
343
344            if domain.split == False:
345                edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
346            elif domain.split == True:
347                edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
348            #Update optimal_timestep
349            try:
350                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
351            except ZeroDivisionError:
352                pass
353
354    domain.timestep = timestep
355
356def compute_fluxes(domain):
357    """Compute all fluxes and the timestep suitable for all volumes
358    in domain.
359
360    Compute total flux for each conserved quantity using "flux_function"
361
362    Fluxes across each edge are scaled by edgelengths and summed up
363    Resulting flux is then scaled by area and stored in
364    explicit_update for each of the three conserved quantities
365    stage, xmomentum and ymomentum
366
367    The maximal allowable speed computed by the flux_function for each volume
368    is converted to a timestep that must not be exceeded. The minimum of
369    those is computed as the next overall timestep.
370
371    Post conditions:
372      domain.explicit_update is reset to computed flux values
373      domain.timestep is set to the largest step satisfying all volumes.
374    """
375
376    import sys
377    from Numeric import zeros, Float
378
379    N = domain.number_of_elements
380
381    #Shortcuts
382    Stage = domain.quantities['stage']
383    Xmom = domain.quantities['xmomentum']
384    Bed = domain.quantities['elevation']
385
386    stage = Stage.vertex_values
387    xmom =  Xmom.vertex_values
388    bed =   Bed.vertex_values
389
390    #Arrays
391    stage_bdry = Stage.boundary_values
392    xmom_bdry =  Xmom.boundary_values
393
394    flux = zeros(2, Float) #Work array for summing up fluxes
395
396
397    #Loop
398    timestep = float(sys.maxint)
399    for k in range(N):
400
401        flux[:] = 0.  #Reset work array
402        for i in range(2):
403            #Quantities inside volume facing neighbour i
404            ql = [stage[k, i], xmom[k, i]]
405            zl = bed[k, i]
406
407            #Quantities at neighbour on nearest face
408            n = domain.neighbours[k,i]
409            if n < 0:
410                m = -n-1 #Convert negative flag to index
411                qr = [stage_bdry[m], xmom_bdry[m]]
412                zr = zl #Extend bed elevation to boundary
413            else:
414                m = domain.neighbour_vertices[k,i]
415                qr = [stage[n, m], xmom[n, m]]
416                zr = bed[n, m]
417
418
419            #Outward pointing normal vector
420            normal = domain.normals[k, i]
421
422            #Flux computation using provided function
423            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
424            flux -= edgeflux
425
426            #Update optimal_timestep on full cells
427            #if  domain.tri_full_flag[k] == 1:
428            try:
429                timestep = min(timestep, 0.5*domain.areas[k]/max_speed)
430            except ZeroDivisionError:
431                pass
432
433        #Normalise by area and store for when all conserved
434        #quantities get updated
435        flux /= domain.areas[k]
436        Stage.explicit_update[k] = flux[0]
437        Xmom.explicit_update[k] = flux[1]
438
439    domain.timestep = timestep
440
441####################################
442# Module functions for gradient limiting (distribute_to_vertices_and_edges)
443
444def distribute_to_vertices_and_edges(domain):
445    """Distribution from centroids to vertices specific to the
446    shallow water wave
447    equation.
448
449    It will ensure that h (w-z) is always non-negative even in the
450    presence of steep bed-slopes by taking a weighted average between shallow
451    and deep cases.
452
453    In addition, all conserved quantities get distributed as per either a
454    constant (order==1) or a piecewise linear function (order==2).
455
456    FIXME: more explanation about removal of artificial variability etc
457
458    Precondition:
459      All quantities defined at centroids and bed elevation defined at
460      vertices.
461
462    Postcondition
463      Conserved quantities defined at vertices
464
465    """
466
467    #Remove very thin layers of water
468    protect_against_infinitesimal_and_negative_heights(domain)
469
470    #Extrapolate all conserved quantities
471    #old code:
472    for name in domain.conserved_quantities:
473        Q = domain.quantities[name]
474        if domain.order == 1:
475            Q.extrapolate_first_order()
476        elif domain.order == 2:
477            Q.extrapolate_second_order()
478        else:
479            raise 'Unknown order'
480
481
482    #Take bed elevation into account when water heights are small
483    balance_deep_and_shallow(domain)
484
485def protect_against_infinitesimal_and_negative_heights(domain):
486    """Protect against infinitesimal heights and associated high velocities
487    """
488
489    #Shortcuts
490    wc = domain.quantities['stage'].centroid_values
491    zc = domain.quantities['elevation'].centroid_values
492    xmomc = domain.quantities['xmomentum'].centroid_values
493    hc = wc - zc  #Water depths at centroids
494
495    #Update
496    #FIXME: Modify accroditg to c-version - or discard altogether.
497    for k in range(domain.number_of_elements):
498
499        if hc[k] < domain.minimum_allowed_height:
500            #Control stage
501            if hc[k] < domain.epsilon:
502                wc[k] = zc[k] # Contain 'lost mass' error
503
504            #Control momentum
505            xmomc[k] = 0.0
506
507
508def h_limiter(domain):
509    """Limit slopes for each volume to eliminate artificial variance
510    introduced by e.g. second order extrapolator
511
512    limit on h = w-z
513
514    This limiter depends on two quantities (w,z) so it resides within
515    this module rather than within quantity.py
516    """
517
518    from Numeric import zeros, Float
519
520    N = domain.number_of_elements
521    beta_h = domain.beta_h
522
523    #Shortcuts
524    wc = domain.quantities['stage'].centroid_values
525    zc = domain.quantities['elevation'].centroid_values
526    hc = wc - zc
527
528    wv = domain.quantities['stage'].vertex_values
529    zv = domain.quantities['elevation'].vertex_values
530    hv = wv-zv
531
532    hvbar = zeros(hv.shape, Float) #h-limited values
533
534    #Find min and max of this and neighbour's centroid values
535    hmax = zeros(hc.shape, Float)
536    hmin = zeros(hc.shape, Float)
537
538    for k in range(N):
539        hmax[k] = hmin[k] = hc[k]
540        for i in range(2):
541            n = domain.neighbours[k,i]
542            if n >= 0:
543                hn = hc[n] #Neighbour's centroid value
544
545                hmin[k] = min(hmin[k], hn)
546                hmax[k] = max(hmax[k], hn)
547
548
549    #Diffences between centroids and maxima/minima
550    dhmax = hmax - hc
551    dhmin = hmin - hc
552
553    #Deltas between vertex and centroid values
554    dh = zeros(hv.shape, Float)
555    for i in range(2):
556        dh[:,i] = hv[:,i] - hc
557
558    #Phi limiter
559    for k in range(N):
560
561        #Find the gradient limiter (phi) across vertices
562        phi = 1.0
563        for i in range(2):
564            r = 1.0
565            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
566            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
567
568            phi = min( min(r*beta_h, 1.0), phi )
569
570        #Then update using phi limiter
571        for i in range(2):
572            hvbar[k,i] = hc[k] + phi*dh[k,i]
573
574    return hvbar
575
576def balance_deep_and_shallow(domain):
577    """Compute linear combination between stage as computed by
578    gradient-limiters limiting using w, and stage computed by
579    gradient-limiters limiting using h (h-limiter).
580    The former takes precedence when heights are large compared to the
581    bed slope while the latter takes precedence when heights are
582    relatively small.  Anything in between is computed as a balanced
583    linear combination in order to avoid numerical disturbances which
584    would otherwise appear as a result of hard switching between
585    modes.
586
587    The h-limiter is always applied irrespective of the order.
588    """
589
590    #Shortcuts
591    wc = domain.quantities['stage'].centroid_values
592    zc = domain.quantities['elevation'].centroid_values
593    hc = wc - zc
594
595    wv = domain.quantities['stage'].vertex_values
596    zv = domain.quantities['elevation'].vertex_values
597    hv = wv-zv
598
599    #Limit h
600    hvbar = h_limiter(domain)
601
602    for k in range(domain.number_of_elements):
603        #Compute maximal variation in bed elevation
604        #  This quantitiy is
605        #    dz = max_i abs(z_i - z_c)
606        #  and it is independent of dimension
607        #  In the 1d case zc = (z0+z1)/2
608        #  In the 2d case zc = (z0+z1+z2)/3
609
610        dz = max(abs(zv[k,0]-zc[k]),
611                 abs(zv[k,1]-zc[k]))
612
613
614        hmin = min( hv[k,:] )
615
616        #Create alpha in [0,1], where alpha==0 means using the h-limited
617        #stage and alpha==1 means using the w-limited stage as
618        #computed by the gradient limiter (both 1st or 2nd order)
619
620        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
621        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
622
623        if dz > 0.0:
624            alpha = max( min( 2.0*hmin/dz, 1.0), 0.0 )
625        else:
626            #Flat bed
627            alpha = 1.0
628
629        #Let
630        #
631        #  wvi be the w-limited stage (wvi = zvi + hvi)
632        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
633        #
634        #
635        #where i=0,1,2 denotes the vertex ids
636        #
637        #Weighted balance between w-limited and h-limited stage is
638        #
639        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
640        #
641        #It follows that the updated wvi is
642        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
643        #
644        # Momentum is balanced between constant and limited
645
646
647        #for i in range(3):
648        #    wv[k,i] = zv[k,i] + hvbar[k,i]
649
650        #return
651
652        if alpha < 1:
653
654            for i in range(2):
655                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
656
657            #Momentums at centroids
658            xmomc = domain.quantities['xmomentum'].centroid_values
659
660            #Momentums at vertices
661            xmomv = domain.quantities['xmomentum'].vertex_values
662
663            # Update momentum as a linear combination of
664            # xmomc and ymomc (shallow) and momentum
665            # from extrapolator xmomv and ymomv (deep).
666            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
667
668
669###############################################
670#Boundaries - specific to the shallow water wave equation
671class Reflective_boundary(Boundary):
672    """Reflective boundary returns same conserved quantities as
673    those present in its neighbour volume but reflected.
674
675    This class is specific to the shallow water equation as it
676    works with the momentum quantities assumed to be the second
677    and third conserved quantities.
678    """
679
680    def __init__(self, domain = None):
681        Boundary.__init__(self)
682
683        if domain is None:
684            msg = 'Domain must be specified for reflective boundary'
685            raise msg
686
687        #Handy shorthands
688        self.stage   = domain.quantities['stage'].vertex_values
689        self.xmom    = domain.quantities['xmomentum'].vertex_values
690        self.normals = domain.normals
691
692        from Numeric import zeros, Float
693        self.conserved_quantities = zeros(2, Float)
694
695    def __repr__(self):
696        return 'Reflective_boundary'
697
698
699    def evaluate(self, vol_id, edge_id):
700        """Reflective boundaries reverses the outward momentum
701        of the volume they serve.
702        """
703
704        q = self.conserved_quantities
705        q[0] = self.stage[vol_id, edge_id]
706        q[1] = self.xmom[vol_id, edge_id]
707
708
709        #r = rotate(q, normal, direction = 1)
710        #r[1] = -r[1]
711        #q = rotate(r, normal, direction = -1)
712        r = q
713        r[1] = -q[1]
714        q = r
715
716        return q
717
718
719#########################
720#Standard forcing terms:
721#
722def gravity(domain):
723    """Apply gravitational pull in the presence of bed slope
724    """
725
726    from Numeric import zeros, Float, array, sum
727    from util import gradient
728
729    xmom = domain.quantities['xmomentum'].explicit_update
730
731    Stage = domain.quantities['stage']
732    Elevation = domain.quantities['elevation']
733    h = Stage.vertex_values - Elevation.vertex_values
734    v = Elevation.vertex_values
735
736    x = domain.get_vertex_coordinates()
737    g = domain.g
738
739    for k in range(domain.number_of_elements):
740        avg_h = sum( h[k,:] )/2.0
741
742        #Compute bed slope
743        x0, x1 = x[k,:]
744        z0, z1 = v[k,:]
745
746        zx = gradient(x0, x1, z0, z1)
747
748        #Update momentumw
749        xmom[k] += -g*zx*avg_h
750
751def manning_friction(domain):
752    """Apply (Manning) friction to water momentum
753    (Python version)
754    """
755
756    from math import sqrt
757
758    w = domain.quantities['stage'].centroid_values
759    z = domain.quantities['elevation'].centroid_values
760    h = w-z
761
762    uh = domain.quantities['xmomentum'].centroid_values
763    eta = domain.quantities['friction'].centroid_values
764
765    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
766
767    N = domain.number_of_elements
768    eps = domain.minimum_allowed_height
769    g = domain.g
770
771    for k in range(N):
772        if eta[k] >= eps:
773            if h[k] >= eps:
774                S = -g * eta[k]**2 * sqrt(uh[k]**2)
775                S /= h[k]**(7.0/3.0)
776
777                #Update momentum
778                xmom_update[k] += S*uh[k]
779
780def check_forcefield(f):
781    """Check that f is either
782    1: a callable object f(t,x,y), where x and y are vectors
783       and that it returns an array or a list of same length
784       as x and y
785    2: a scalar
786    """
787
788    from Numeric import ones, Float, array
789
790
791    if callable(f):
792        N = 2
793        x = ones(2, Float)
794        y = ones(2, Float)
795        try:
796            q = f(1.0, x=x, y=y)
797        except Exception, e:
798            msg = 'Function %s could not be executed:\n%s' %(f, e)
799            #FIXME: Reconsider this semantics
800            raise msg
801
802        try:
803            q = array(q).astype(Float)
804        except:
805            msg = 'Return value from vector function %s could ' %f
806            msg += 'not be converted into a Numeric array of floats.\n'
807            msg += 'Specified function should return either list or array.'
808            raise msg
809
810        #Is this really what we want?
811        msg = 'Return vector from function %s ' %f
812        msg += 'must have same lenght as input vectors'
813        assert len(q) == N, msg
814
815    else:
816        try:
817            f = float(f)
818        except:
819            msg = 'Force field %s must be either a scalar' %f
820            msg += ' or a vector function'
821            raise Exception(msg)
822    return f
823
824
825###########################
826###########################
827#Geometries
828
829class Constant_stage:
830    """Set an initial condition with constant stage
831    """
832    def __init__(self, s):
833        self.s = s
834
835    def __call__(self, x, y):
836        return self.s
837
838class Constant_height:
839    """Set an initial condition with constant water height, e.g
840    stage s = z+h
841    """
842
843    def __init__(self, W, h):
844        self.W = W
845        self.h = h
846
847    def __call__(self, x, y):
848        if self.W is None:
849            from Numeric import ones, Float
850            return self.h*ones(len(x), Float)
851        else:
852            return self.W(x,y) + self.h
853
854
855
856"""
857def compute_fluxes_python(domain):
858    Compute all fluxes and the timestep suitable for all volumes
859    in domain.
860
861    Compute total flux for each conserved quantity using "flux_function"
862
863    Fluxes across each edge are scaled by edgelengths and summed up
864    Resulting flux is then scaled by area and stored in
865    explicit_update for each of the three conserved quantities
866    stage, xmomentum and ymomentum
867
868    The maximal allowable speed computed by the flux_function for each volume
869    is converted to a timestep that must not be exceeded. The minimum of
870    those is computed as the next overall timestep.
871
872    Post conditions:
873      domain.explicit_update is reset to computed flux values
874      domain.timestep is set to the largest step satisfying all volumes.
875   
876
877    import sys
878    from Numeric import zeros, Float
879
880    N = domain.number_of_elements
881
882    #Shortcuts
883    Stage = domain.quantities['stage']
884    Xmom = domain.quantities['xmomentum']
885    Bed = domain.quantities['elevation']
886
887    #Arrays
888    stage = Stage.edge_values
889    xmom =  Xmom.edge_values
890    bed =   Bed.edge_values
891
892    stage_bdry = Stage.boundary_values
893    xmom_bdry =  Xmom.boundary_values
894
895    flux = zeros((N,2), Float) #Work array for summing up fluxes
896
897    #Loop
898    timestep = float(sys.maxint)
899    for k in range(N):
900
901        for i in range(2):
902            #Quantities inside volume facing neighbour i
903            ql = [stage[k, i], xmom[k, i]]
904            zl = bed[k, i]
905
906            #Quantities at neighbour on nearest face
907            n = domain.neighbours[k,i]
908            if n < 0:
909                m = -n-1 #Convert negative flag to index
910                qr = [stage_bdry[m], xmom_bdry[m]]
911                zr = zl #Extend bed elevation to boundary
912            else:
913                m = domain.neighbour_edges[k,i]
914                qr = [stage[n, m], xmom[n, m]]
915                zr = bed[n, m]
916
917
918            #Outward pointing normal vector
919            normal = domain.normals[k, 2*i:2*i+2]
920
921            #Flux computation using provided function
922            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
923
924            flux[k,:] = edgeflux
925
926    return flux
927"""
928
929
930
931
932
933
934##############################################
935#Initialise module
936
937
938#Optimisation with psyco
939from config import use_psyco
940if use_psyco:
941    try:
942        import psyco
943    except:
944        import os
945        if os.name == 'posix' and os.uname()[4] == 'x86_64':
946            pass
947            #Psyco isn't supported on 64 bit systems, but it doesn't matter
948        else:
949            msg = 'WARNING: psyco (speedup) could not import'+\
950                  ', you may want to consider installing it'
951            print msg
952    else:
953        psyco.bind(Domain.distribute_to_vertices_and_edges)
954        psyco.bind(Domain.compute_fluxes)
955
956if __name__ == "__main__":
957    pass
958
959# Profiling stuff
960import profile
961profiler = profile.Profile()
Note: See TracBrowser for help on using the repository browser.