source: anuga_work/development/euler/euler.py @ 3970

Last change on this file since 3970 was 3514, checked in by duncan, 19 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 51.4 KB
Line 
1"""Class Domain -
22D triangular domains for finite-volume computations of
3the Euler 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 + G_y = 0
10
11where
12
13U = [rho, mx, my, E]
14u = mx/rho
15v = my/rho
16
17E = [u rho, u mx + p, u my, u(E + p)]
18G = [v rho, v mx, v my + p, v(E + p)]
19
20and _t, _x, _y denote the derivative with respect to t, x and y respectively.
21
22The quantities are
23
24symbol    variable name    explanation
25x         x                horizontal distance from origin [m]
26y         y                vertical distance from origin [m]
27rho       density          density [kg/m^2]
28u                          speed in the x direction [m/s]
29v                          speed in the y direction [m/s]
30mx        xmomentum        momentum density in the x direction [kg/m/s^2]
31my        ymomentum        momentum density in the y direction [kg/m/s^2]
32E         energy           energy density [J/m^2]
33gamma     gamma            gamma law gas (1.4 for air)
34p         pressure         pressure [N/m
35
36
37The conserved quantities are rho, mx, my, E
38
39Gamma law relates pressure with conserved quantities
40
41p = (gamma-1)[E - rho/2(u^2 + v^2)]
42
43
44Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
45Geoscience Australia, 2004
46"""
47
48#Subversion keywords:
49#
50#$LastChangedDate: 2005-10-14 12:02:18 +1000 (Fri, 14 Oct 2005) $
51#$LastChangedRevision: 1922 $
52#$LastChangedBy: duncan $
53
54
55from domain import *
56from region import *
57
58Generic_domain = Domain #Rename
59
60#Euler domain
61class Domain(Generic_domain):
62
63    def __init__(self, coordinates, vertices, boundary = None,
64                 tagged_elements = None, geo_reference = None,
65                 use_inscribed_circle=False):
66
67        conserved_quantities = ['density', 'xmomentum', 'ymomentum', 'energy']
68        other_quantities = ['pressure']
69        Generic_domain.__init__(self, coordinates, vertices, boundary,
70                                conserved_quantities, other_quantities,
71                                tagged_elements, geo_reference, use_inscribed_circle)
72
73        from euler_config import minimum_allowed_density
74        self.minimum_allowed_density = minimum_allowed_density
75
76        #Realtime visualisation
77        self.visualiser = None
78        self.visualise  = False
79        self.visualise_color_stage = False
80        self.visualise_stage_range = 1.0
81        self.visualise_timer = True
82        self.visualise_range_z = None
83
84        #Stored output
85        self.store = False
86        self.format = 'sww'
87        self.smooth = False
88
89        #Reduction operation for get_vertex_values
90        from anuga.pyvolution.util import mean
91        self.reduction = mean
92        #self.reduction = min  #Looks better near steep slopes
93
94        self.quantities_to_be_stored = ['density', 'xmomentum', 'ymomentum', 'energy', 'pressure']
95
96
97        #Establish shortcuts to relevant quantities (for efficiency)
98        self.rho = self.quantities['density']
99        self.mx = self.quantities['xmomentum']
100        self.my = self.quantities['ymomentum']
101        self.E = self.quantities['energy']
102
103
104    def initialise_visualiser(self,scale_z=1.0,rect=None):
105        #Realtime visualisation
106        if self.visualiser is None:
107            from realtime_visualisation_new import Visualiser
108            self.visualiser = Visualiser(self,scale_z,rect)
109        self.visualise = True
110
111    def check_integrity(self):
112        Generic_domain.check_integrity(self)
113
114        #Check that we are solving the Euler equation
115
116        msg = 'First conserved quantity must be "density"'
117        assert self.conserved_quantities[0] == 'stage', msg
118        msg = 'Second conserved quantity must be "xmomentum"'
119        assert self.conserved_quantities[1] == 'xmomentum', msg
120        msg = 'Third conserved quantity must be "ymomentum"'
121        assert self.conserved_quantities[2] == 'ymomentum', msg
122        msg = 'Fourth conserved quantity must be "energy"'
123        assert self.conserved_quantities[2] == 'ymomentum', msg
124
125    def extrapolate_second_order_sw(self):
126        #Call correct module function
127        #(either from this module or C-extension)
128        extrapolate_second_order_sw(self)
129
130    def compute_fluxes(self):
131        #Call correct module function
132        #(either from this module or C-extension)
133        compute_fluxes(self)
134
135    def distribute_to_vertices_and_edges(self):
136        #Call correct module function
137        #(either from this module or C-extension)
138        distribute_to_vertices_and_edges(self)
139
140
141    def evolve(self,
142               yieldstep = None,
143               finaltime = None,
144               duration = None,               
145               skip_initial_step = False):
146       
147        """Specialisation of basic evolve method from parent class
148        """
149
150        #Call check integrity here rather than from user scripts
151        #self.check_integrity()
152
153        msg = 'Parameter beta must be in the interval [0, 1['
154        assert 0 <= self.beta < 1.0, msg
155
156        #Initial update of vertex and edge values before any storage
157        #and or visualisation
158        self.distribute_to_vertices_and_edges()
159
160
161        #Initialise real time viz if requested
162        if self.visualise is True and self.time == 0.0:
163
164            import realtime_visualisation_new as visualise
165            self.initialise_visualiser()
166            self.visualiser.setup_all()
167            self.visualiser.update_timer()
168
169        #Store model data, e.g. for visualisation
170        if self.store is True and self.time == 0.0:
171            self.initialise_storage()
172        else:
173            pass
174
175        #Call basic machinery from parent class
176        #for t in Generic_domain.evolve(self, yieldstep, finaltime,
177        #                               skip_initial_step):
178
179        for t in Generic_domain.evolve(self,
180                                       yieldstep=yieldstep,
181                                       finaltime=finaltime,
182                                       duration=duration,
183                                       skip_initial_step=skip_initial_step):           
184            #Real time viz
185            if self.visualise is True:
186                self.visualiser.update_all()
187                self.visualiser.update_timer()
188
189            #Store model data, e.g. for subsequent visualisation
190            if self.store is True:
191                self.store_timestep(self.quantities_to_be_stored)
192
193            #FIXME: Could maybe be taken from specified list
194            #of 'store every step' quantities
195
196            #Pass control on to outer loop for more specific actions
197            yield(t)
198
199
200    def initialise_storage(self):
201        """
202        Create and initialise self.writer object for storing data.
203        Also, save x,y and bed elevation
204        """
205
206        import anuga.pyvolution.data_manager
207
208        #Initialise writer
209        self.writer = data_manager.get_dataobject(self, mode = 'w')
210
211        #Store vertices and connectivity
212        self.writer.store_connectivity()
213
214
215    def store_timestep(self, name):
216        """Store named quantity and time.
217
218        Precondition:
219           self.write has been initialised
220        """
221        self.writer.store_timestep(name)
222
223
224####################MH 090605 new extrapolation function belonging to domain class
225def extrapolate_second_order_sw(domain):
226    """extrapolate conserved quntities to the vertices of the triangles
227    Python version to be written after the C version
228    """
229    msg = 'Method extrapolate_second_order_sw should be implemented in C'
230    raise msg
231####################MH 090605 ###########################################
232
233    def initialise_visualiser(self,scale_z=1.0,rect=None):
234        #Realtime visualisation
235        if self.visualiser is None:
236            from realtime_visualisation_new import Visualiser
237            from Numeric import array, Float
238
239            if rect is None:
240                rect = array(self.xy_extent, Float)
241
242            self.visualiser = Visualiser(self,scale_z,rect)
243        self.visualise = True
244
245
246#Rotation of momentum vector
247def rotate(q, normal, direction = 1):
248    """Rotate the momentum component q (q[1], q[2])
249    from x,y coordinates to coordinates based on normal vector.
250
251    If direction is negative the rotation is inverted.
252
253    Input vector is preserved
254
255    This function is specific to the shallow water wave equation
256    """
257
258    from Numeric import zeros, Float
259
260    assert len(q) == 3,\
261           'Vector of conserved quantities must have length 3'\
262           'for 2D shallow water equation'
263
264    try:
265        l = len(normal)
266    except:
267        raise 'Normal vector must be an Numeric array'
268
269    assert l == 2, 'Normal vector must have 2 components'
270
271
272    n1 = normal[0]
273    n2 = normal[1]
274
275    r = zeros(len(q), Float) #Rotated quantities
276    r[0] = q[0]              #First quantity, height, is not rotated
277
278    if direction == -1:
279        n2 = -n2
280
281
282    r[1] =  n1*q[1] + n2*q[2]
283    r[2] = -n2*q[1] + n1*q[2]
284
285    return r
286
287
288
289####################################
290# Flux computation
291def flux_function(normal, ql, qr, zl, zr):
292    """Compute fluxes between volumes for the shallow water wave equation
293    cast in terms of w = h+z using the 'central scheme' as described in
294
295    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
296    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
297    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
298
299    The implemented formula is given in equation (3.15) on page 714
300
301    Conserved quantities w, uh, vh are stored as elements 0, 1 and 2
302    in the numerical vectors ql an qr.
303
304    Bed elevations zl and zr.
305    """
306
307    from anuga.pyvolution.config import g, epsilon
308    from math import sqrt
309    from Numeric import array
310
311    #Align momentums with x-axis
312    q_left  = rotate(ql, normal, direction = 1)
313    q_right = rotate(qr, normal, direction = 1)
314
315    z = (zl+zr)/2 #Take average of field values
316
317    w_left  = q_left[0]   #w=h+z
318    h_left  = w_left-z
319    uh_left = q_left[1]
320
321    if h_left < epsilon:
322        u_left = 0.0  #Could have been negative
323        h_left = 0.0
324    else:
325        u_left  = uh_left/h_left
326
327
328    w_right  = q_right[0]  #w=h+z
329    h_right  = w_right-z
330    uh_right = q_right[1]
331
332
333    if h_right < epsilon:
334        u_right = 0.0  #Could have been negative
335        h_right = 0.0
336    else:
337        u_right  = uh_right/h_right
338
339    vh_left  = q_left[2]
340    vh_right = q_right[2]
341
342    soundspeed_left  = sqrt(g*h_left)
343    soundspeed_right = sqrt(g*h_right)
344
345    #Maximal wave speed
346    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
347
348    #Minimal wave speed
349    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
350
351    #Flux computation
352
353    #FIXME(Ole): Why is it again that we don't
354    #use uh_left and uh_right directly in the first entries?
355    flux_left  = array([u_left*h_left,
356                        u_left*uh_left + 0.5*g*h_left**2,
357                        u_left*vh_left])
358    flux_right = array([u_right*h_right,
359                        u_right*uh_right + 0.5*g*h_right**2,
360                        u_right*vh_right])
361
362    denom = s_max-s_min
363    if denom == 0.0:
364        edgeflux = array([0.0, 0.0, 0.0])
365        max_speed = 0.0
366    else:
367        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
368        edgeflux += s_max*s_min*(q_right-q_left)/denom
369
370        edgeflux = rotate(edgeflux, normal, direction=-1)
371        max_speed = max(abs(s_max), abs(s_min))
372
373    return edgeflux, max_speed
374
375
376def compute_fluxes(domain):
377    """Compute all fluxes and the timestep suitable for all volumes
378    in domain.
379
380    Compute total flux for each conserved quantity using "flux_function"
381
382    Fluxes across each edge are scaled by edgelengths and summed up
383    Resulting flux is then scaled by area and stored in
384    explicit_update for each of the three conserved quantities
385    stage, xmomentum and ymomentum
386
387    The maximal allowable speed computed by the flux_function for each volume
388    is converted to a timestep that must not be exceeded. The minimum of
389    those is computed as the next overall timestep.
390
391    Post conditions:
392      domain.explicit_update is reset to computed flux values
393      domain.timestep is set to the largest step satisfying all volumes.
394    """
395
396    import sys
397    from Numeric import zeros, Float
398
399    N = domain.number_of_elements
400
401    #Shortcuts
402    Stage = domain.quantities['stage']
403    Xmom = domain.quantities['xmomentum']
404    Ymom = domain.quantities['ymomentum']
405    Bed = domain.quantities['elevation']
406
407    #Arrays
408    stage = Stage.edge_values
409    xmom =  Xmom.edge_values
410    ymom =  Ymom.edge_values
411    bed =   Bed.edge_values
412
413    stage_bdry = Stage.boundary_values
414    xmom_bdry =  Xmom.boundary_values
415    ymom_bdry =  Ymom.boundary_values
416
417    flux = zeros(3, Float) #Work array for summing up fluxes
418
419
420    #Loop
421    timestep = float(sys.maxint)
422    for k in range(N):
423
424        flux[:] = 0.  #Reset work array
425        for i in range(3):
426            #Quantities inside volume facing neighbour i
427            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
428            zl = bed[k, i]
429
430            #Quantities at neighbour on nearest face
431            n = domain.neighbours[k,i]
432            if n < 0:
433                m = -n-1 #Convert negative flag to index
434                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
435                zr = zl #Extend bed elevation to boundary
436            else:
437                m = domain.neighbour_edges[k,i]
438                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
439                zr = bed[n, m]
440
441
442            #Outward pointing normal vector
443            normal = domain.normals[k, 2*i:2*i+2]
444
445            #Flux computation using provided function
446            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
447            flux -= edgeflux * domain.edgelengths[k,i]
448
449            #Update optimal_timestep
450            try:
451                timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
452            except ZeroDivisionError:
453                pass
454
455        #Normalise by area and store for when all conserved
456        #quantities get updated
457        flux /= domain.areas[k]
458        Stage.explicit_update[k] = flux[0]
459        Xmom.explicit_update[k] = flux[1]
460        Ymom.explicit_update[k] = flux[2]
461
462
463    domain.timestep = timestep
464
465#MH090605 The following method belongs to the shallow_water domain class
466#see comments in the corresponding method in shallow_water_ext.c
467def extrapolate_second_order_sw_c(domain):
468    """Wrapper calling C version of extrapolate_second_order_sw
469    """
470    import sys
471    from Numeric import zeros, Float
472
473    N = domain.number_of_elements
474
475    #Shortcuts
476    Stage = domain.quantities['stage']
477    Xmom = domain.quantities['xmomentum']
478    Ymom = domain.quantities['ymomentum']
479    from shallow_water_ext import extrapolate_second_order_sw
480    extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
481                                domain.number_of_boundaries,
482                                domain.centroid_coordinates,
483                                Stage.centroid_values,
484                                Xmom.centroid_values,
485                                Ymom.centroid_values,
486                                domain.vertex_coordinates,
487                                Stage.vertex_values,
488                                Xmom.vertex_values,
489                                Ymom.vertex_values)
490
491def compute_fluxes_c(domain):
492    """Wrapper calling C version of compute fluxes
493    """
494
495    import sys
496    from Numeric import zeros, Float
497
498    N = domain.number_of_elements
499
500    #Shortcuts
501    Stage = domain.quantities['stage']
502    Xmom = domain.quantities['xmomentum']
503    Ymom = domain.quantities['ymomentum']
504    Bed = domain.quantities['elevation']
505
506    timestep = float(sys.maxint)
507    from shallow_water_ext import compute_fluxes
508    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
509                                     domain.neighbours,
510                                     domain.neighbour_edges,
511                                     domain.normals,
512                                     domain.edgelengths,
513                                     domain.radii,
514                                     domain.areas,
515                                     Stage.edge_values,
516                                     Xmom.edge_values,
517                                     Ymom.edge_values,
518                                     Bed.edge_values,
519                                     Stage.boundary_values,
520                                     Xmom.boundary_values,
521                                     Ymom.boundary_values,
522                                     Stage.explicit_update,
523                                     Xmom.explicit_update,
524                                     Ymom.explicit_update,
525                                     domain.already_computed_flux)
526
527
528####################################
529# Module functions for gradient limiting (distribute_to_vertices_and_edges)
530
531def distribute_to_vertices_and_edges(domain):
532    """Distribution from centroids to vertices specific to the
533    shallow water wave
534    equation.
535
536    It will ensure that h (w-z) is always non-negative even in the
537    presence of steep bed-slopes by taking a weighted average between shallow
538    and deep cases.
539
540    In addition, all conserved quantities get distributed as per either a
541    constant (order==1) or a piecewise linear function (order==2).
542
543    FIXME: more explanation about removal of artificial variability etc
544
545    Precondition:
546      All quantities defined at centroids and bed elevation defined at
547      vertices.
548
549    Postcondition
550      Conserved quantities defined at vertices
551
552    """
553
554    from config import optimised_gradient_limiter
555
556    #Remove very thin layers of water
557    protect_against_infinitesimal_and_negative_heights(domain)
558
559    #Extrapolate all conserved quantities
560    if optimised_gradient_limiter:
561        #MH090605 if second order,
562        #perform the extrapolation and limiting on
563        #all of the conserved quantitie
564
565        if (domain.order == 1):
566            for name in domain.conserved_quantities:
567                Q = domain.quantities[name]
568                Q.extrapolate_first_order()
569        elif domain.order == 2:
570            domain.extrapolate_second_order_sw()
571        else:
572            raise 'Unknown order'
573    else:
574        #old code:
575        for name in domain.conserved_quantities:
576            Q = domain.quantities[name]
577            if domain.order == 1:
578                Q.extrapolate_first_order()
579            elif domain.order == 2:
580                Q.extrapolate_second_order()
581                Q.limit()
582            else:
583                raise 'Unknown order'
584
585
586    #Take bed elevation into account when water heights are small
587    balance_deep_and_shallow(domain)
588
589    #Compute edge values by interpolation
590    for name in domain.conserved_quantities:
591        Q = domain.quantities[name]
592        Q.interpolate_from_vertices_to_edges()
593
594
595def protect_against_infinitesimal_and_negative_heights(domain):
596    """Protect against infinitesimal heights and associated high velocities
597    """
598
599    #Shortcuts
600    wc = domain.quantities['stage'].centroid_values
601    zc = domain.quantities['elevation'].centroid_values
602    xmomc = domain.quantities['xmomentum'].centroid_values
603    ymomc = domain.quantities['ymomentum'].centroid_values
604    hc = wc - zc  #Water depths at centroids
605
606    #Update
607    for k in range(domain.number_of_elements):
608
609        if hc[k] < domain.minimum_allowed_height:
610            #Control stage
611            if hc[k] < domain.epsilon:
612                wc[k] = zc[k] # Contain 'lost mass' error
613
614            #Control momentum
615            xmomc[k] = ymomc[k] = 0.0
616
617
618def protect_against_infinitesimal_and_negative_heights_c(domain):
619    """Protect against infinitesimal heights and associated high velocities
620    """
621
622    #Shortcuts
623    wc = domain.quantities['stage'].centroid_values
624    zc = domain.quantities['elevation'].centroid_values
625    xmomc = domain.quantities['xmomentum'].centroid_values
626    ymomc = domain.quantities['ymomentum'].centroid_values
627
628    from shallow_water_ext import protect
629
630    protect(domain.minimum_allowed_height, domain.epsilon,
631            wc, zc, xmomc, ymomc)
632
633
634
635def h_limiter(domain):
636    """Limit slopes for each volume to eliminate artificial variance
637    introduced by e.g. second order extrapolator
638
639    limit on h = w-z
640
641    This limiter depends on two quantities (w,z) so it resides within
642    this module rather than within quantity.py
643    """
644
645    from Numeric import zeros, Float
646
647    N = domain.number_of_elements
648    beta_h = domain.beta_h
649
650    #Shortcuts
651    wc = domain.quantities['stage'].centroid_values
652    zc = domain.quantities['elevation'].centroid_values
653    hc = wc - zc
654
655    wv = domain.quantities['stage'].vertex_values
656    zv = domain.quantities['elevation'].vertex_values
657    hv = wv-zv
658
659    hvbar = zeros(hv.shape, Float) #h-limited values
660
661    #Find min and max of this and neighbour's centroid values
662    hmax = zeros(hc.shape, Float)
663    hmin = zeros(hc.shape, Float)
664
665    for k in range(N):
666        hmax[k] = hmin[k] = hc[k]
667        for i in range(3):
668            n = domain.neighbours[k,i]
669            if n >= 0:
670                hn = hc[n] #Neighbour's centroid value
671
672                hmin[k] = min(hmin[k], hn)
673                hmax[k] = max(hmax[k], hn)
674
675
676    #Diffences between centroids and maxima/minima
677    dhmax = hmax - hc
678    dhmin = hmin - hc
679
680    #Deltas between vertex and centroid values
681    dh = zeros(hv.shape, Float)
682    for i in range(3):
683        dh[:,i] = hv[:,i] - hc
684
685    #Phi limiter
686    for k in range(N):
687
688        #Find the gradient limiter (phi) across vertices
689        phi = 1.0
690        for i in range(3):
691            r = 1.0
692            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
693            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
694
695            phi = min( min(r*beta_h, 1), phi )
696
697        #Then update using phi limiter
698        for i in range(3):
699            hvbar[k,i] = hc[k] + phi*dh[k,i]
700
701    return hvbar
702
703
704
705def h_limiter_c(domain):
706    """Limit slopes for each volume to eliminate artificial variance
707    introduced by e.g. second order extrapolator
708
709    limit on h = w-z
710
711    This limiter depends on two quantities (w,z) so it resides within
712    this module rather than within quantity.py
713
714    Wrapper for c-extension
715    """
716
717    from Numeric import zeros, Float
718
719    N = domain.number_of_elements
720    beta_h = domain.beta_h
721
722    #Shortcuts
723    wc = domain.quantities['stage'].centroid_values
724    zc = domain.quantities['elevation'].centroid_values
725    hc = wc - zc
726
727    wv = domain.quantities['stage'].vertex_values
728    zv = domain.quantities['elevation'].vertex_values
729    hv = wv - zv
730
731    #Call C-extension
732    from shallow_water_ext import h_limiter_sw as h_limiter
733    hvbar = h_limiter(domain, hc, hv)
734
735    return hvbar
736
737
738def balance_deep_and_shallow(domain):
739    """Compute linear combination between stage as computed by
740    gradient-limiters limiting using w, and stage computed by
741    gradient-limiters limiting using h (h-limiter).
742    The former takes precedence when heights are large compared to the
743    bed slope while the latter takes precedence when heights are
744    relatively small.  Anything in between is computed as a balanced
745    linear combination in order to avoid numerical disturbances which
746    would otherwise appear as a result of hard switching between
747    modes.
748
749    The h-limiter is always applied irrespective of the order.
750    """
751
752    #Shortcuts
753    wc = domain.quantities['stage'].centroid_values
754    zc = domain.quantities['elevation'].centroid_values
755    hc = wc - zc
756
757    wv = domain.quantities['stage'].vertex_values
758    zv = domain.quantities['elevation'].vertex_values
759    hv = wv-zv
760
761    #Limit h
762    hvbar = h_limiter(domain)
763
764    for k in range(domain.number_of_elements):
765        #Compute maximal variation in bed elevation
766        #  This quantitiy is
767        #    dz = max_i abs(z_i - z_c)
768        #  and it is independent of dimension
769        #  In the 1d case zc = (z0+z1)/2
770        #  In the 2d case zc = (z0+z1+z2)/3
771
772        dz = max(abs(zv[k,0]-zc[k]),
773                 abs(zv[k,1]-zc[k]),
774                 abs(zv[k,2]-zc[k]))
775
776
777        hmin = min( hv[k,:] )
778
779        #Create alpha in [0,1], where alpha==0 means using the h-limited
780        #stage and alpha==1 means using the w-limited stage as
781        #computed by the gradient limiter (both 1st or 2nd order)
782
783        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
784        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
785
786        if dz > 0.0:
787            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
788        else:
789            #Flat bed
790            alpha = 1.0
791
792        #Let
793        #
794        #  wvi be the w-limited stage (wvi = zvi + hvi)
795        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
796        #
797        #
798        #where i=0,1,2 denotes the vertex ids
799        #
800        #Weighted balance between w-limited and h-limited stage is
801        #
802        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
803        #
804        #It follows that the updated wvi is
805        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
806        #
807        # Momentum is balanced between constant and limited
808
809
810        #for i in range(3):
811        #    wv[k,i] = zv[k,i] + hvbar[k,i]
812
813        #return
814
815        if alpha < 1:
816
817            for i in range(3):
818                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
819
820            #Momentums at centroids
821            xmomc = domain.quantities['xmomentum'].centroid_values
822            ymomc = domain.quantities['ymomentum'].centroid_values
823
824            #Momentums at vertices
825            xmomv = domain.quantities['xmomentum'].vertex_values
826            ymomv = domain.quantities['ymomentum'].vertex_values
827
828            # Update momentum as a linear combination of
829            # xmomc and ymomc (shallow) and momentum
830            # from extrapolator xmomv and ymomv (deep).
831            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
832            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
833
834
835def balance_deep_and_shallow_c(domain):
836    """Wrapper for C implementation
837    """
838
839    #Shortcuts
840    wc = domain.quantities['stage'].centroid_values
841    zc = domain.quantities['elevation'].centroid_values
842    hc = wc - zc
843
844    wv = domain.quantities['stage'].vertex_values
845    zv = domain.quantities['elevation'].vertex_values
846    hv = wv - zv
847
848    #Momentums at centroids
849    xmomc = domain.quantities['xmomentum'].centroid_values
850    ymomc = domain.quantities['ymomentum'].centroid_values
851
852    #Momentums at vertices
853    xmomv = domain.quantities['xmomentum'].vertex_values
854    ymomv = domain.quantities['ymomentum'].vertex_values
855
856    #Limit h
857    hvbar = h_limiter(domain)
858
859    #This is how one would make a first order h_limited value
860    #as in the old balancer (pre 17 Feb 2005):
861    #from Numeric import zeros, Float
862    #hvbar = zeros( (len(hc), 3), Float)
863    #for i in range(3):
864    #    hvbar[:,i] = hc[:]
865
866    from shallow_water_ext import balance_deep_and_shallow
867    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
868                             xmomc, ymomc, xmomv, ymomv)
869
870
871
872
873###############################################
874#Boundaries - specific to the shallow water wave equation
875class Reflective_boundary(Boundary):
876    """Reflective boundary returns same conserved quantities as
877    those present in its neighbour volume but reflected.
878
879    This class is specific to the shallow water equation as it
880    works with the momentum quantities assumed to be the second
881    and third conserved quantities.
882    """
883
884    def __init__(self, domain = None):
885        Boundary.__init__(self)
886
887        if domain is None:
888            msg = 'Domain must be specified for reflective boundary'
889            raise msg
890
891        #Handy shorthands
892        self.stage   = domain.quantities['stage'].edge_values
893        self.xmom    = domain.quantities['xmomentum'].edge_values
894        self.ymom    = domain.quantities['ymomentum'].edge_values
895        self.normals = domain.normals
896
897        from Numeric import zeros, Float
898        self.conserved_quantities = zeros(3, Float)
899
900    def __repr__(self):
901        return 'Reflective_boundary'
902
903
904    def evaluate(self, vol_id, edge_id):
905        """Reflective boundaries reverses the outward momentum
906        of the volume they serve.
907        """
908
909        q = self.conserved_quantities
910        q[0] = self.stage[vol_id, edge_id]
911        q[1] = self.xmom[vol_id, edge_id]
912        q[2] = self.ymom[vol_id, edge_id]
913
914        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
915
916
917        r = rotate(q, normal, direction = 1)
918        r[1] = -r[1]
919        q = rotate(r, normal, direction = -1)
920
921        return q
922
923
924
925class Transmissive_Momentum_Set_Stage_boundary(Boundary):
926    """Returns same momentum conserved quantities as
927    those present in its neighbour volume. Sets stage
928
929    Underlying domain must be specified when boundary is instantiated
930    """
931
932    def __init__(self, domain = None, function=None):
933        Boundary.__init__(self)
934
935        if domain is None:
936            msg = 'Domain must be specified for this type boundary'
937            raise msg
938
939        if function is None:
940            msg = 'Function must be specified for this type boundary'
941            raise msg
942
943        self.domain   = domain
944        self.function = function
945
946    def __repr__(self):
947        return 'Transmissive_Momentum_Set_Stage_boundary(%s)' %self.domain
948
949    def evaluate(self, vol_id, edge_id):
950        """Transmissive Momentum Set Stage boundaries return the edge momentum
951        values of the volume they serve.
952        """
953
954        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
955        value = self.function(self.domain.time)
956        q[0] = value[0]
957        return q
958
959
960        #FIXME: Consider this (taken from File_boundary) to allow
961        #spatial variation
962        #if vol_id is not None and edge_id is not None:
963        #    i = self.boundary_indices[ vol_id, edge_id ]
964        #    return self.F(t, point_id = i)
965        #else:
966        #    return self.F(t)
967
968
969
970class Dirichlet_Discharge_boundary(Boundary):
971    """Sets stage (h0) - FIXME: Should we use w0?
972    Sets momentum (wh0) in the inward normal direction.
973
974    Underlying domain must be specified when boundary is instantiated
975    """
976
977    def __init__(self, domain = None, h0=None, wh0=None):
978        Boundary.__init__(self)
979
980        if domain is None:
981            msg = 'Domain must be specified for this type boundary'
982            raise msg
983
984        if h0 is None:
985            h0 = 0.0
986
987        if wh0 is None:
988            wh0 = 0.0
989
990        self.domain   = domain
991        self.h0  = h0
992        self.wh0 = wh0
993
994    def __repr__(self):
995        return 'Dirichlet_Discharge_boundary(%s)' %self.domain
996
997    def evaluate(self, vol_id, edge_id):
998        """Set discharge in the (inward) normal direction
999        """
1000
1001        normal = self.domain.get_normal(vol_id,edge_id)
1002        q = [self.h0, -self.wh0*normal[0], -self.wh0*normal[1]]
1003        return q
1004
1005
1006        #FIXME: Consider this (taken from File_boundary) to allow
1007        #spatial variation
1008        #if vol_id is not None and edge_id is not None:
1009        #    i = self.boundary_indices[ vol_id, edge_id ]
1010        #    return self.F(t, point_id = i)
1011        #else:
1012        #    return self.F(t)
1013
1014
1015
1016#class Spatio_temporal_boundary(Boundary):
1017#    """The spatio-temporal boundary, reads values for the conserved
1018#    quantities from an sww NetCDF file, and returns interpolated values
1019#    at the midpoints of each associated boundaty segment.
1020#    Time dependency is interpolated linearly as in util.File_function.#
1021#
1022#    Example:
1023#    Bf = Spatio_temporal_boundary('source_file.sww', domain)
1024#
1025#    """
1026Spatio_temporal_boundary = File_boundary
1027
1028
1029
1030
1031#########################
1032#Standard forcing terms:
1033#
1034def gravity(domain):
1035    """Apply gravitational pull in the presence of bed slope
1036    """
1037
1038    from anuga.pyvolution.util import gradient
1039    from Numeric import zeros, Float, array, sum
1040
1041    xmom = domain.quantities['xmomentum'].explicit_update
1042    ymom = domain.quantities['ymomentum'].explicit_update
1043
1044    Stage = domain.quantities['stage']
1045    Elevation = domain.quantities['elevation']
1046    h = Stage.edge_values - Elevation.edge_values
1047    v = Elevation.vertex_values
1048
1049    x = domain.get_vertex_coordinates()
1050    g = domain.g
1051
1052    for k in range(domain.number_of_elements):
1053        avg_h = sum( h[k,:] )/3
1054
1055        #Compute bed slope
1056        x0, y0, x1, y1, x2, y2 = x[k,:]
1057        z0, z1, z2 = v[k,:]
1058
1059        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
1060
1061        #Update momentum
1062        xmom[k] += -g*zx*avg_h
1063        ymom[k] += -g*zy*avg_h
1064
1065
1066def gravity_c(domain):
1067    """Wrapper calling C version
1068    """
1069
1070    xmom = domain.quantities['xmomentum'].explicit_update
1071    ymom = domain.quantities['ymomentum'].explicit_update
1072
1073    Stage = domain.quantities['stage']
1074    Elevation = domain.quantities['elevation']
1075    h = Stage.edge_values - Elevation.edge_values
1076    v = Elevation.vertex_values
1077
1078    x = domain.get_vertex_coordinates()
1079    g = domain.g
1080
1081
1082    from shallow_water_ext import gravity
1083    gravity(g, h, v, x, xmom, ymom)
1084
1085
1086def manning_friction(domain):
1087    """Apply (Manning) friction to water momentum
1088    """
1089
1090    from math import sqrt
1091
1092    w = domain.quantities['stage'].centroid_values
1093    z = domain.quantities['elevation'].centroid_values
1094    h = w-z
1095
1096    uh = domain.quantities['xmomentum'].centroid_values
1097    vh = domain.quantities['ymomentum'].centroid_values
1098    eta = domain.quantities['friction'].centroid_values
1099
1100    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1101    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1102
1103    N = domain.number_of_elements
1104    eps = domain.minimum_allowed_height
1105    g = domain.g
1106
1107    for k in range(N):
1108        if eta[k] >= eps:
1109            if h[k] >= eps:
1110                S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1111                S /= h[k]**(7.0/3)
1112
1113                #Update momentum
1114                xmom_update[k] += S*uh[k]
1115                ymom_update[k] += S*vh[k]
1116
1117
1118def manning_friction_c(domain):
1119    """Wrapper for c version
1120    """
1121
1122
1123    xmom = domain.quantities['xmomentum']
1124    ymom = domain.quantities['ymomentum']
1125
1126    w = domain.quantities['stage'].centroid_values
1127    z = domain.quantities['elevation'].centroid_values
1128
1129    uh = xmom.centroid_values
1130    vh = ymom.centroid_values
1131    eta = domain.quantities['friction'].centroid_values
1132
1133    xmom_update = xmom.semi_implicit_update
1134    ymom_update = ymom.semi_implicit_update
1135
1136    N = domain.number_of_elements
1137    eps = domain.minimum_allowed_height
1138    g = domain.g
1139
1140    from shallow_water_ext import manning_friction
1141    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1142
1143
1144def linear_friction(domain):
1145    """Apply linear friction to water momentum
1146
1147    Assumes quantity: 'linear_friction' to be present
1148    """
1149
1150    from math import sqrt
1151
1152    w = domain.quantities['stage'].centroid_values
1153    z = domain.quantities['elevation'].centroid_values
1154    h = w-z
1155
1156    uh = domain.quantities['xmomentum'].centroid_values
1157    vh = domain.quantities['ymomentum'].centroid_values
1158    tau = domain.quantities['linear_friction'].centroid_values
1159
1160    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1161    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1162
1163    N = domain.number_of_elements
1164    eps = domain.minimum_allowed_height
1165    g = domain.g #Not necessary? Why was this added?
1166
1167    for k in range(N):
1168        if tau[k] >= eps:
1169            if h[k] >= eps:
1170                S = -tau[k]/h[k]
1171
1172                #Update momentum
1173                xmom_update[k] += S*uh[k]
1174                ymom_update[k] += S*vh[k]
1175
1176
1177
1178def check_forcefield(f):
1179    """Check that f is either
1180    1: a callable object f(t,x,y), where x and y are vectors
1181       and that it returns an array or a list of same length
1182       as x and y
1183    2: a scalar
1184    """
1185
1186    from Numeric import ones, Float, array
1187
1188
1189    if callable(f):
1190        N = 3
1191        x = ones(3, Float)
1192        y = ones(3, Float)
1193        try:
1194            q = f(1.0, x=x, y=y)
1195        except Exception, e:
1196            msg = 'Function %s could not be executed:\n%s' %(f, e)
1197            #FIXME: Reconsider this semantics
1198            raise msg
1199
1200        try:
1201            q = array(q).astype(Float)
1202        except:
1203            msg = 'Return value from vector function %s could ' %f
1204            msg += 'not be converted into a Numeric array of floats.\n'
1205            msg += 'Specified function should return either list or array.'
1206            raise msg
1207
1208        #Is this really what we want?
1209        msg = 'Return vector from function %s ' %f
1210        msg += 'must have same lenght as input vectors'
1211        assert len(q) == N, msg
1212
1213    else:
1214        try:
1215            f = float(f)
1216        except:
1217            msg = 'Force field %s must be either a scalar' %f
1218            msg += ' or a vector function'
1219            raise msg
1220    return f
1221
1222
1223class Wind_stress:
1224    """Apply wind stress to water momentum in terms of
1225    wind speed [m/s] and wind direction [degrees]
1226    """
1227
1228    def __init__(self, *args, **kwargs):
1229        """Initialise windfield from wind speed s [m/s]
1230        and wind direction phi [degrees]
1231
1232        Inputs v and phi can be either scalars or Python functions, e.g.
1233
1234        W = Wind_stress(10, 178)
1235
1236        #FIXME - 'normal' degrees are assumed for now, i.e. the
1237        vector (1,0) has zero degrees.
1238        We may need to convert from 'compass' degrees later on and also
1239        map from True north to grid north.
1240
1241        Arguments can also be Python functions of t,x,y as in
1242
1243        def speed(t,x,y):
1244            ...
1245            return s
1246
1247        def angle(t,x,y):
1248            ...
1249            return phi
1250
1251        where x and y are vectors.
1252
1253        and then pass the functions in
1254
1255        W = Wind_stress(speed, angle)
1256
1257        The instantiated object W can be appended to the list of
1258        forcing_terms as in
1259
1260        Alternatively, one vector valued function for (speed, angle)
1261        can be applied, providing both quantities simultaneously.
1262        As in
1263        W = Wind_stress(F), where returns (speed, angle) for each t.
1264
1265        domain.forcing_terms.append(W)
1266        """
1267
1268        from config import rho_a, rho_w, eta_w
1269        from Numeric import array, Float
1270
1271        if len(args) == 2:
1272            s = args[0]
1273            phi = args[1]
1274        elif len(args) == 1:
1275            #Assume vector function returning (s, phi)(t,x,y)
1276            vector_function = args[0]
1277            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1278            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1279        else:
1280           #Assume info is in 2 keyword arguments
1281
1282           if len(kwargs) == 2:
1283               s = kwargs['s']
1284               phi = kwargs['phi']
1285           else:
1286               raise 'Assumes two keyword arguments: s=..., phi=....'
1287
1288        self.speed = check_forcefield(s)
1289        self.phi = check_forcefield(phi)
1290
1291        self.const = eta_w*rho_a/rho_w
1292
1293
1294    def __call__(self, domain):
1295        """Evaluate windfield based on values found in domain
1296        """
1297
1298        from math import pi, cos, sin, sqrt
1299        from Numeric import Float, ones, ArrayType
1300
1301        xmom_update = domain.quantities['xmomentum'].explicit_update
1302        ymom_update = domain.quantities['ymomentum'].explicit_update
1303
1304        N = domain.number_of_elements
1305        t = domain.time
1306
1307        if callable(self.speed):
1308            xc = domain.get_centroid_coordinates()
1309            s_vec = self.speed(t, xc[:,0], xc[:,1])
1310        else:
1311            #Assume s is a scalar
1312
1313            try:
1314                s_vec = self.speed * ones(N, Float)
1315            except:
1316                msg = 'Speed must be either callable or a scalar: %s' %self.s
1317                raise msg
1318
1319
1320        if callable(self.phi):
1321            xc = domain.get_centroid_coordinates()
1322            phi_vec = self.phi(t, xc[:,0], xc[:,1])
1323        else:
1324            #Assume phi is a scalar
1325
1326            try:
1327                phi_vec = self.phi * ones(N, Float)
1328            except:
1329                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1330                raise msg
1331
1332        assign_windfield_values(xmom_update, ymom_update,
1333                                s_vec, phi_vec, self.const)
1334
1335
1336def assign_windfield_values(xmom_update, ymom_update,
1337                            s_vec, phi_vec, const):
1338    """Python version of assigning wind field to update vectors.
1339    A c version also exists (for speed)
1340    """
1341    from math import pi, cos, sin, sqrt
1342
1343    N = len(s_vec)
1344    for k in range(N):
1345        s = s_vec[k]
1346        phi = phi_vec[k]
1347
1348        #Convert to radians
1349        phi = phi*pi/180
1350
1351        #Compute velocity vector (u, v)
1352        u = s*cos(phi)
1353        v = s*sin(phi)
1354
1355        #Compute wind stress
1356        S = const * sqrt(u**2 + v**2)
1357        xmom_update[k] += S*u
1358        ymom_update[k] += S*v
1359
1360
1361
1362##############################
1363#OBSOLETE STUFF
1364
1365def balance_deep_and_shallow_old(domain):
1366    """Compute linear combination between stage as computed by
1367    gradient-limiters and stage computed as constant height above bed.
1368    The former takes precedence when heights are large compared to the
1369    bed slope while the latter takes precedence when heights are
1370    relatively small.  Anything in between is computed as a balanced
1371    linear combination in order to avoid numerical disturbances which
1372    would otherwise appear as a result of hard switching between
1373    modes.
1374    """
1375
1376    #OBSOLETE
1377
1378    #Shortcuts
1379    wc = domain.quantities['stage'].centroid_values
1380    zc = domain.quantities['elevation'].centroid_values
1381    hc = wc - zc
1382
1383    wv = domain.quantities['stage'].vertex_values
1384    zv = domain.quantities['elevation'].vertex_values
1385    hv = wv-zv
1386
1387
1388    #Computed linear combination between constant stages and and
1389    #stages parallel to the bed elevation.
1390    for k in range(domain.number_of_elements):
1391        #Compute maximal variation in bed elevation
1392        #  This quantitiy is
1393        #    dz = max_i abs(z_i - z_c)
1394        #  and it is independent of dimension
1395        #  In the 1d case zc = (z0+z1)/2
1396        #  In the 2d case zc = (z0+z1+z2)/3
1397
1398        dz = max(abs(zv[k,0]-zc[k]),
1399                 abs(zv[k,1]-zc[k]),
1400                 abs(zv[k,2]-zc[k]))
1401
1402
1403        hmin = min( hv[k,:] )
1404
1405        #Create alpha in [0,1], where alpha==0 means using shallow
1406        #first order scheme and alpha==1 means using the stage w as
1407        #computed by the gradient limiter (1st or 2nd order)
1408        #
1409        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
1410        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
1411
1412        if dz > 0.0:
1413            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
1414        else:
1415            #Flat bed
1416            alpha = 1.0
1417
1418        #Weighted balance between stage parallel to bed elevation
1419        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
1420        #order gradient limiter
1421        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
1422        #
1423        #It follows that the updated wvi is
1424        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
1425        #  zvi + hc + alpha*(hvi - hc)
1426        #
1427        #Note that hvi = zc+hc-zvi in the first order case (constant).
1428
1429        if alpha < 1:
1430            for i in range(3):
1431                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
1432
1433
1434            #Momentums at centroids
1435            xmomc = domain.quantities['xmomentum'].centroid_values
1436            ymomc = domain.quantities['ymomentum'].centroid_values
1437
1438            #Momentums at vertices
1439            xmomv = domain.quantities['xmomentum'].vertex_values
1440            ymomv = domain.quantities['ymomentum'].vertex_values
1441
1442            # Update momentum as a linear combination of
1443            # xmomc and ymomc (shallow) and momentum
1444            # from extrapolator xmomv and ymomv (deep).
1445            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
1446            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
1447
1448
1449
1450
1451
1452###########################
1453###########################
1454#Geometries
1455
1456
1457#FIXME: Rethink this way of creating values.
1458
1459
1460class Weir:
1461    """Set a bathymetry for weir with a hole and a downstream gutter
1462    x,y are assumed to be in the unit square
1463    """
1464
1465    def __init__(self, stage):
1466        self.inflow_stage = stage
1467
1468    def __call__(self, x, y):
1469        from Numeric import zeros, Float
1470        from math import sqrt
1471
1472        N = len(x)
1473        assert N == len(y)
1474
1475        z = zeros(N, Float)
1476        for i in range(N):
1477            z[i] = -x[i]/2  #General slope
1478
1479            #Flattish bit to the left
1480            if x[i] < 0.3:
1481                z[i] = -x[i]/10
1482
1483            #Weir
1484            if x[i] >= 0.3 and x[i] < 0.4:
1485                z[i] = -x[i]+0.9
1486
1487            #Dip
1488            x0 = 0.6
1489            #depth = -1.3
1490            depth = -1.0
1491            #plateaux = -0.9
1492            plateaux = -0.6
1493            if y[i] < 0.7:
1494                if x[i] > x0 and x[i] < 0.9:
1495                    z[i] = depth
1496
1497                #RHS plateaux
1498                if x[i] >= 0.9:
1499                    z[i] = plateaux
1500
1501
1502            elif y[i] >= 0.7 and y[i] < 1.5:
1503                #Restrict and deepen
1504                if x[i] >= x0 and x[i] < 0.8:
1505                    z[i] = depth-(y[i]/3-0.3)
1506                    #z[i] = depth-y[i]/5
1507                    #z[i] = depth
1508                elif x[i] >= 0.8:
1509                    #RHS plateaux
1510                    z[i] = plateaux
1511
1512            elif y[i] >= 1.5:
1513                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
1514                    #Widen up and stay at constant depth
1515                    z[i] = depth-1.5/5
1516                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
1517                    #RHS plateaux
1518                    z[i] = plateaux
1519
1520
1521            #Hole in weir (slightly higher than inflow condition)
1522            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
1523                z[i] = -x[i]+self.inflow_stage + 0.02
1524
1525            #Channel behind weir
1526            x0 = 0.5
1527            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
1528                z[i] = -x[i]+self.inflow_stage + 0.02
1529
1530            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
1531                #Flatten it out towards the end
1532                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
1533
1534            #Hole to the east
1535            x0 = 1.1; y0 = 0.35
1536            #if x[i] < -0.2 and y < 0.5:
1537            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
1538                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
1539
1540            #Tiny channel draining hole
1541            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
1542                z[i] = -0.9 #North south
1543
1544            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
1545                z[i] = -1.0 + (x[i]-0.9)/3 #East west
1546
1547
1548
1549            #Stuff not in use
1550
1551            #Upward slope at inlet to the north west
1552            #if x[i] < 0.0: # and y[i] > 0.5:
1553            #    #z[i] = -y[i]+0.5  #-x[i]/2
1554            #    z[i] = x[i]/4 - y[i]**2 + 0.5
1555
1556            #Hole to the west
1557            #x0 = -0.4; y0 = 0.35 # center
1558            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
1559            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
1560
1561
1562
1563
1564
1565        return z/2
1566
1567class Weir_simple:
1568    """Set a bathymetry for weir with a hole and a downstream gutter
1569    x,y are assumed to be in the unit square
1570    """
1571
1572    def __init__(self, stage):
1573        self.inflow_stage = stage
1574
1575    def __call__(self, x, y):
1576        from Numeric import zeros, Float
1577
1578        N = len(x)
1579        assert N == len(y)
1580
1581        z = zeros(N, Float)
1582        for i in range(N):
1583            z[i] = -x[i]  #General slope
1584
1585            #Flat bit to the left
1586            if x[i] < 0.3:
1587                z[i] = -x[i]/10  #General slope
1588
1589            #Weir
1590            if x[i] > 0.3 and x[i] < 0.4:
1591                z[i] = -x[i]+0.9
1592
1593            #Dip
1594            if x[i] > 0.6 and x[i] < 0.9:
1595                z[i] = -x[i]-0.5  #-y[i]/5
1596
1597            #Hole in weir (slightly higher than inflow condition)
1598            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
1599                z[i] = -x[i]+self.inflow_stage + 0.05
1600
1601
1602        return z/2
1603
1604
1605
1606class Constant_stage:
1607    """Set an initial condition with constant stage
1608    """
1609    def __init__(self, s):
1610        self.s = s
1611
1612    def __call__(self, x, y):
1613        return self.s
1614
1615class Constant_height:
1616    """Set an initial condition with constant water height, e.g
1617    stage s = z+h
1618    """
1619
1620    def __init__(self, W, h):
1621        self.W = W
1622        self.h = h
1623
1624    def __call__(self, x, y):
1625        if self.W is None:
1626            from Numeric import ones, Float
1627            return self.h*ones(len(x), Float)
1628        else:
1629            return self.W(x,y) + self.h
1630
1631
1632
1633
1634def compute_fluxes_python(domain):
1635    """Compute all fluxes and the timestep suitable for all volumes
1636    in domain.
1637
1638    Compute total flux for each conserved quantity using "flux_function"
1639
1640    Fluxes across each edge are scaled by edgelengths and summed up
1641    Resulting flux is then scaled by area and stored in
1642    explicit_update for each of the three conserved quantities
1643    stage, xmomentum and ymomentum
1644
1645    The maximal allowable speed computed by the flux_function for each volume
1646    is converted to a timestep that must not be exceeded. The minimum of
1647    those is computed as the next overall timestep.
1648
1649    Post conditions:
1650      domain.explicit_update is reset to computed flux values
1651      domain.timestep is set to the largest step satisfying all volumes.
1652    """
1653
1654    import sys
1655    from Numeric import zeros, Float
1656
1657    N = domain.number_of_elements
1658
1659    #Shortcuts
1660    Stage = domain.quantities['stage']
1661    Xmom = domain.quantities['xmomentum']
1662    Ymom = domain.quantities['ymomentum']
1663    Bed = domain.quantities['elevation']
1664
1665    #Arrays
1666    stage = Stage.edge_values
1667    xmom =  Xmom.edge_values
1668    ymom =  Ymom.edge_values
1669    bed =   Bed.edge_values
1670
1671    stage_bdry = Stage.boundary_values
1672    xmom_bdry =  Xmom.boundary_values
1673    ymom_bdry =  Ymom.boundary_values
1674
1675    flux = zeros((N,3), Float) #Work array for summing up fluxes
1676
1677    #Loop
1678    timestep = float(sys.maxint)
1679    for k in range(N):
1680
1681        for i in range(3):
1682            #Quantities inside volume facing neighbour i
1683            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
1684            zl = bed[k, i]
1685
1686            #Quantities at neighbour on nearest face
1687            n = domain.neighbours[k,i]
1688            if n < 0:
1689                m = -n-1 #Convert negative flag to index
1690                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
1691                zr = zl #Extend bed elevation to boundary
1692            else:
1693                m = domain.neighbour_edges[k,i]
1694                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
1695                zr = bed[n, m]
1696
1697
1698            #Outward pointing normal vector
1699            normal = domain.normals[k, 2*i:2*i+2]
1700
1701            #Flux computation using provided function
1702            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
1703
1704            flux[k,:] = edgeflux
1705
1706    return flux
1707
1708
1709
1710
1711
1712
1713
1714##############################################
1715#Initialise module
1716
1717
1718from anuga.utilities import compile
1719if compile.can_use_C_extension('shallow_water_ext.c'):
1720    #Replace python version with c implementations
1721
1722    from shallow_water_ext import rotate, assign_windfield_values
1723    compute_fluxes = compute_fluxes_c
1724    extrapolate_second_order_sw=extrapolate_second_order_sw_c
1725    gravity = gravity_c
1726    manning_friction = manning_friction_c
1727    h_limiter = h_limiter_c
1728    balance_deep_and_shallow = balance_deep_and_shallow_c
1729    protect_against_infinitesimal_and_negative_heights = protect_against_infinitesimal_and_negative_heights_c
1730
1731
1732    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c #(like MH's)
1733
1734
1735
1736#Optimisation with psyco
1737from config import use_psyco
1738if use_psyco:
1739    try:
1740        import psyco
1741    except:
1742        import os
1743        if os.name == 'posix' and os.uname()[4] == 'x86_64':
1744            pass
1745            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1746        else:
1747            msg = 'WARNING: psyco (speedup) could not import'+\
1748                  ', you may want to consider installing it'
1749            print msg
1750    else:
1751        psyco.bind(Domain.distribute_to_vertices_and_edges)
1752        psyco.bind(Domain.compute_fluxes)
1753
1754if __name__ == "__main__":
1755    pass
Note: See TracBrowser for help on using the repository browser.