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

Last change on this file since 7578 was 4494, checked in by steve, 17 years ago

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