source: development/pyvolution-1d/shallow_water_minmod.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 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: 40.8 KB
Line 
1"""Class Domain -
21D interval domains for finite-volume computations of
3the shallow water wave equation.
4
5This module contains a specialisation of class Domain from module domain.py
6consisting of methods specific to the Shallow Water Wave Equation
7
8
9U_t + E_x = S
10
11where
12
13U = [w, uh]
14E = [uh, u^2h + gh^2/2]
15S represents source terms forcing the system
16(e.g. gravity, friction, wind stress, ...)
17
18and _t, _x, _y denote the derivative with respect to t, x and y respectively.
19
20The quantities are
21
22symbol    variable name    explanation
23x         x                horizontal distance from origin [m]
24z         elevation        elevation of bed on which flow is modelled [m]
25h         height           water height above z [m]
26w         stage            absolute water level, w = z+h [m]
27u                          speed in the x direction [m/s]
28uh        xmomentum        momentum in the x direction [m^2/s]
29
30eta                        mannings friction coefficient [to appear]
31nu                         wind stress coefficient [to appear]
32
33The conserved quantities are w, uh
34
35For details see e.g.
36Christopher Zoppou and Stephen Roberts,
37Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
38Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
39
40
41John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
42Geoscience Australia, 2006
43"""
44
45from domain import *
46Generic_Domain = Domain #Rename
47
48#Shallow water domain
49class Domain(Generic_Domain):
50
51    def __init__(self, coordinates, boundary = None, tagged_elements = None,
52                 geo_reference = None):
53
54        conserved_quantities = ['stage', 'xmomentum']
55        other_quantities = ['elevation', 'friction']
56        Generic_Domain.__init__(self, coordinates, boundary,
57                                conserved_quantities, other_quantities,
58                                tagged_elements, geo_reference)
59       
60        from config import minimum_allowed_height, g
61        self.minimum_allowed_height = minimum_allowed_height
62        self.g = g
63
64        #forcing terms not included in 1d domain ?WHy?
65        self#.forcing_terms.append(gravity)
66        #self.forcing_terms.append(manning_friction)
67        #print "\nI have Removed forcing terms line 64 1dsw"
68
69        #Realtime visualisation
70        self.visualiser = None
71        self.visualise  = False
72        self.visualise_color_stage = False
73        self.visualise_stage_range = 1.0
74        self.visualise_timer = True
75        self.visualise_range_z = None
76       
77        #Stored output
78        self.store = True
79        self.format = 'sww'
80        self.smooth = True
81
82        #Evolve parametrs
83        self.cfl = 1.0
84       
85        #Reduction operation for get_vertex_values
86        from anuga.pyvolution.util import mean
87        self.reduction = mean
88        #self.reduction = min  #Looks better near steep slopes
89
90        self.quantities_to_be_stored = ['stage','xmomentum']
91
92    def set_quantities_to_be_stored(self, q):
93        """Specify which quantities will be stored in the sww file.
94
95        q must be either:
96          - the name of a quantity
97          - a list of quantity names
98          - None
99
100        In the two first cases, the named quantities will be stored at each
101        yieldstep
102        (This is in addition to the quantities elevation and friction) 
103        If q is None, storage will be switched off altogether.
104        """
105
106
107        if q is None:
108            self.quantities_to_be_stored = []           
109            self.store = False
110            return
111
112        if isinstance(q, basestring):
113            q = [q] # Turn argument into a list
114
115        #Check correcness   
116        for quantity_name in q:
117            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
118            assert quantity_name in self.conserved_quantities, msg
119       
120        self.quantities_to_be_stored = q
121       
122
123    def initialise_visualiser(self,scale_z=1.0,rect=None):
124        #Realtime visualisation
125        if self.visualiser is None:
126            from realtime_visualisation_new import Visualiser
127            self.visualiser = Visualiser(self,scale_z,rect)
128            self.visualiser.setup['elevation']=True
129            self.visualiser.updating['stage']=True
130        self.visualise = True
131        if self.visualise_color_stage == True:
132            self.visualiser.coloring['stage'] = True
133            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
134        print 'initialise visualiser'
135        print self.visualiser.setup
136        print self.visualiser.updating
137
138    def check_integrity(self):
139        Generic_Domain.check_integrity(self)
140        #Check that we are solving the shallow water wave equation
141
142        msg = 'First conserved quantity must be "stage"'
143        assert self.conserved_quantities[0] == 'stage', msg
144        msg = 'Second conserved quantity must be "xmomentum"'
145        assert self.conserved_quantities[1] == 'xmomentum', msg
146
147    def extrapolate_second_order_sw(self):
148        #Call correct module function
149        #(either from this module or C-extension)
150        extrapolate_second_order_sw(self)
151
152    def compute_fluxes(self):
153        #Call correct module function
154        #(either from this module or C-extension)
155        compute_fluxes(self)
156
157    def distribute_to_vertices_and_edges(self):
158        #Call correct module function
159        #(either from this module or C-extension)
160        distribute_to_vertices_and_edges(self)
161
162    def evolve(self, yieldstep = None, finaltime = None,
163               skip_initial_step = False):
164        """Specialisation of basic evolve method from parent class
165        """
166
167        #Call check integrity here rather than from user scripts
168        #self.check_integrity()
169
170        #msg = 'Parameter beta_h must be in the interval [0, 1)'
171        #assert 0 <= self.beta_h < 1.0, msg
172        #msg = 'Parameter beta_w must be in the interval [0, 1)'
173        #assert 0 <= self.beta_w < 1.0, msg
174
175
176        #Initial update of vertex and edge values before any storage
177        #and or visualisation
178        self.distribute_to_vertices_and_edges()
179       
180        #Initialise real time viz if requested
181        #if self.visualise is True and self.time == 0.0:
182        #    if self.visualiser is None:
183        #        self.initialise_visualiser()
184        #
185        #    self.visualiser.update_timer()
186        #    self.visualiser.setup_all()
187
188        #Store model data, e.g. for visualisation
189        #if self.store is True and self.time == 0.0:
190        #    self.initialise_storage()
191        #    #print 'Storing results in ' + self.writer.filename
192        #else:
193        #    pass
194        #    #print 'Results will not be stored.'
195        #    #print 'To store results set domain.store = True'
196        #    #FIXME: Diagnostic output should be controlled by
197        #    # a 'verbose' flag living in domain (or in a parent class)
198
199        #Call basic machinery from parent class
200        for t in Generic_Domain.evolve(self, yieldstep, finaltime,
201                                       skip_initial_step):
202            #Real time viz
203        #    if self.visualise is True:
204        #        self.visualiser.update_all()
205        #        self.visualiser.update_timer()
206
207
208            #Store model data, e.g. for subsequent visualisation
209        #    if self.store is True:
210        #        self.store_timestep(self.quantities_to_be_stored)
211
212            #FIXME: Could maybe be taken from specified list
213            #of 'store every step' quantities
214
215            #Pass control on to outer loop for more specific actions
216            yield(t)
217
218    def initialise_storage(self):
219        """Create and initialise self.writer object for storing data.
220        Also, save x and bed elevation
221        """
222
223        import anuga.pyvolution.data_manager
224
225        #Initialise writer
226        self.writer = data_manager.get_dataobject(self, mode = 'w')
227
228        #Store vertices and connectivity
229        self.writer.store_connectivity()
230
231
232    def store_timestep(self, name):
233        """Store named quantity and time.
234
235        Precondition:
236           self.write has been initialised
237        """
238        self.writer.store_timestep(name)
239
240
241#=============== End of Shallow Water Domain ===============================
242
243#Rotation of momentum vector
244def rotate(q, normal, direction = 1):
245    """Rotate the momentum component q (q[1], q[2])
246    from x,y coordinates to coordinates based on normal vector.
247
248    If direction is negative the rotation is inverted.
249
250    Input vector is preserved
251
252    This function is specific to the shallow water wave equation
253    """
254
255    from Numeric import zeros, Float
256
257    assert len(q) == 3,\
258           'Vector of conserved quantities must have length 3'\
259           'for 2D shallow water equation'
260
261    try:
262        l = len(normal)
263    except:
264        raise 'Normal vector must be an Numeric array'
265
266    assert l == 2, 'Normal vector must have 2 components'
267
268
269    n1 = normal[0]
270    n2 = normal[1]
271
272    r = zeros(len(q), Float) #Rotated quantities
273    r[0] = q[0]              #First quantity, height, is not rotated
274
275    if direction == -1:
276        n2 = -n2
277
278
279    r[1] =  n1*q[1] + n2*q[2]
280    r[2] = -n2*q[1] + n1*q[2]
281
282    return r
283
284
285# Flux computation
286#def flux_function(normal, ql, qr, zl, zr):
287def flux_function(normal, ql, qr, zl, zr):
288    """Compute fluxes between volumes for the shallow water wave equation
289    cast in terms of w = h+z using the 'central scheme' as described in
290
291    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
292    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
293    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
294
295    The implemented formula is given in equation (3.15) on page 714
296
297    Conserved quantities w, uh, are stored as elements 0 and 1
298    in the numerical vectors ql an qr.
299
300    Bed elevations zl and zr.
301    """
302
303    from config import g, epsilon
304    from math import sqrt
305    from Numeric import array
306
307    #print 'ql',ql
308
309    #Align momentums with x-axis
310    #q_left  = rotate(ql, normal, direction = 1)
311    #q_right = rotate(qr, normal, direction = 1)
312    q_left = ql
313    q_left[1] = q_left[1]*normal
314    q_right = qr
315    q_right[1] = q_right[1]*normal
316
317    #z = (zl+zr)/2 #Take average of field values
318    z = 0.5*(zl+zr) #Take average of field values
319
320    w_left  = q_left[0]   #w=h+z
321    h_left  = w_left-z
322    uh_left = q_left[1]
323
324    if h_left < epsilon:
325        u_left = 0.0  #Could have been negative
326        h_left = 0.0
327    else:
328        u_left  = uh_left/h_left
329
330
331    w_right  = q_right[0]  #w=h+z
332    h_right  = w_right-z
333    uh_right = q_right[1]
334
335
336    if h_right < epsilon:
337        u_right = 0.0  #Could have been negative
338        h_right = 0.0
339    else:
340        u_right  = uh_right/h_right
341
342    #vh_left  = q_left[2]
343    #vh_right = q_right[2]
344
345    soundspeed_left  = sqrt(g*h_left)
346    soundspeed_right = sqrt(g*h_right)
347
348    #Maximal wave speed
349    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
350   
351    #Minimal wave speed
352    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
353   
354    #Flux computation
355
356    #flux_left  = array([u_left*h_left,
357    #                    u_left*uh_left + 0.5*g*h_left**2])
358    #flux_right = array([u_right*h_right,
359    #                    u_right*uh_right + 0.5*g*h_right**2])
360    flux_left  = array([u_left*h_left,
361                        u_left*uh_left + 0.5*g*h_left*h_left])
362    flux_right = array([u_right*h_right,
363                        u_right*uh_right + 0.5*g*h_right*h_right])
364
365    denom = s_max-s_min
366    if denom == 0.0:
367        edgeflux = array([0.0, 0.0])
368        max_speed = 0.0
369    else:
370        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
371        edgeflux += s_max*s_min*(q_right-q_left)/denom
372       
373        edgeflux[1] = edgeflux[1]*normal
374
375        max_speed = max(abs(s_max), abs(s_min))
376
377    return edgeflux, max_speed
378
379def compute_fluxes(domain):
380    """Compute all fluxes and the timestep suitable for all volumes
381    in domain.
382
383    Compute total flux for each conserved quantity using "flux_function"
384
385    Fluxes across each edge are scaled by edgelengths and summed up
386    Resulting flux is then scaled by area and stored in
387    explicit_update for each of the three conserved quantities
388    stage, xmomentum and ymomentum
389
390    The maximal allowable speed computed by the flux_function for each volume
391    is converted to a timestep that must not be exceeded. The minimum of
392    those is computed as the next overall timestep.
393
394    Post conditions:
395      domain.explicit_update is reset to computed flux values
396      domain.timestep is set to the largest step satisfying all volumes.
397    """
398
399    import sys
400    from Numeric import zeros, Float
401
402    N = domain.number_of_elements
403
404    #Shortcuts
405    Stage = domain.quantities['stage']
406    Xmom = domain.quantities['xmomentum']
407#    Ymom = domain.quantities['ymomentum']
408    Bed = domain.quantities['elevation']
409
410    #Arrays
411    #stage = Stage.edge_values
412    #xmom =  Xmom.edge_values
413 #   ymom =  Ymom.edge_values
414    #bed =   Bed.edge_values
415   
416    stage = Stage.vertex_values
417    xmom =  Xmom.vertex_values
418    bed =   Bed.vertex_values
419   
420    #print 'stage edge values', stage
421    #print 'xmom edge values', xmom
422    #print 'bed values', bed
423
424    stage_bdry = Stage.boundary_values
425    xmom_bdry =  Xmom.boundary_values
426    #print 'stage_bdry',stage_bdry
427    #print 'xmom_bdry', xmom_bdry
428#    ymom_bdry =  Ymom.boundary_values
429
430#    flux = zeros(3, Float) #Work array for summing up fluxes
431    flux = zeros(2, Float) #Work array for summing up fluxes
432    ql = zeros(2, Float)
433    qr = zeros(2, Float)
434
435    #Loop
436    timestep = float(sys.maxint)
437    enter = True
438    for k in range(N):
439
440        flux[:] = 0.  #Reset work array
441        #for i in range(3):
442        for i in range(2):
443            #Quantities inside volume facing neighbour i
444            #ql[0] = stage[k, i]
445            #ql[1] = xmom[k, i]
446            ql = [stage[k, i], xmom[k, i]]
447            zl = bed[k, i]
448
449            #Quantities at neighbour on nearest face
450            n = domain.neighbours[k,i]
451            if n < 0:
452                m = -n-1 #Convert negative flag to index
453                qr[0] = stage_bdry[m]
454                qr[1] = xmom_bdry[m]
455                zr = zl #Extend bed elevation to boundary
456            else:
457                #m = domain.neighbour_edges[k,i]
458                m = domain.neighbour_vertices[k,i]
459                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
460                qr[0] = stage[n, m]
461                qr[1] =  xmom[n, m]
462                zr = bed[n, m]
463
464
465            #Outward pointing normal vector
466            normal = domain.normals[k, i]
467       
468            #Flux computation using provided function
469            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
470            #print 'ql',ql
471            #print 'qr',qr
472           
473            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
474            #print 'edgeflux', edgeflux
475
476            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
477            # flux = edgefluxleft - edgefluxright
478            flux -= edgeflux #* domain.edgelengths[k,i]
479            #Update optimal_timestep
480            try:
481                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
482                #timestep = 0.01
483           
484                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
485                if (timestep < 1e-6) & (enter == True):
486                    #print "domain.order", domain.order
487                    #domain.write_time()
488                    print "cell number", k
489                    print "timestep", timestep
490                    print 'max_speed', max_speed
491                    s = domain.quantities['stage']
492                    s = s.centroid_values
493                    xm = domain.quantities['xmomentum']
494                    xm = xm.centroid_values
495                    print 'h', s[k]
496                    print 'xm', xm[k]
497                    print 'u', xm[k]/s[k]
498                    enter = False
499           
500            except ZeroDivisionError:
501                pass
502
503        #Normalise by area and store for when all conserved
504        #quantities get updated
505        flux /= domain.areas[k]
506        # ADD ABOVE LINE AGAIN
507        Stage.explicit_update[k] = flux[0]
508        Xmom.explicit_update[k] = flux[1]
509        #Ymom.explicit_update[k] = flux[2]
510
511
512    domain.timestep = timestep
513
514## #see comments in the corresponding method in shallow_water_ext.c
515## def extrapolate_second_order_sw_c(domain):
516##     """Wrapper calling C version of extrapolate_second_order_sw
517##     """
518##     import sys
519##     from Numeric import zeros, Float
520
521##     N = domain.number_of_elements
522
523##     #Shortcuts
524##     Stage = domain.quantities['stage']
525##     Xmom = domain.quantities['xmomentum']
526## #    Ymom = domain.quantities['ymomentum']
527##     from shallow_water_ext import extrapolate_second_order_sw
528##     extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
529## # cant find this in domain                      domain.number_of_boundaries,
530##                                 domain.centroid_coordinates,
531##                                 Stage.centroid_values,
532##                                 Xmom.centroid_values,
533## #                                Ymom.centroid_values,
534##                                 domain.vertex_coordinates,
535##                                 Stage.vertex_values,
536##                                 Xmom.vertex_values)#,
537## #                                Ymom.vertex_values)
538
539## def compute_fluxes_c(domain):
540##     """Wrapper calling C version of compute fluxes
541##     """
542
543##     import sys
544##     from Numeric import zeros, Float
545
546##     N = domain.number_of_elements
547
548##     #Shortcuts
549##     Stage = domain.quantities['stage']
550##     Xmom = domain.quantities['xmomentum']
551##     #Ymom = domain.quantities['ymomentum']
552##     Bed = domain.quantities['elevation']
553
554##     timestep = float(sys.maxint)
555##     from shallow_water_ext import compute_fluxes
556##     domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
557##                                      domain.neighbours,
558##                                      #domain.neighbour_edges,
559##                                      domain.neighbour_vertices,
560##                                      domain.normals,
561##                                      #domain.edgelengths,
562##                                      #domain.radii,
563##                                      domain.areas,
564##                                      #Stage.edge_values,
565##                                      #Xmom.edge_values,
566##                                      #Ymom.edge_values,
567##                                      #Bed.edge_values,
568##                                      Stage.vertex_values,
569##                                      Xmom.vertex_values,
570##                                      Bed.vertex_values,
571##                                      Stage.boundary_values,
572##                                      Xmom.boundary_values,
573##                                      #Ymom.boundary_values,
574##                                      Stage.explicit_update,
575##                                      Xmom.explicit_update,
576##                                      #Ymom.explicit_update,
577##                                      domain.already_computed_flux)
578
579
580####################################
581
582# Module functions for gradient limiting (distribute_to_vertices_and_edges)
583
584def distribute_to_vertices_and_edges(domain):
585    """Distribution from centroids to vertices specific to the
586    shallow water wave
587    equation.
588
589    It will ensure that h (w-z) is always non-negative even in the
590    presence of steep bed-slopes by taking a weighted average between shallow
591    and deep cases.
592
593    In addition, all conserved quantities get distributed as per either a
594    constant (order==1) or a piecewise linear function (order==2).
595
596    FIXME: more explanation about removal of artificial variability etc
597
598    Precondition:
599      All quantities defined at centroids and bed elevation defined at
600      vertices.
601
602    Postcondition
603      Conserved quantities defined at vertices
604
605    """
606
607    #from config import optimised_gradient_limiter
608
609    #Remove very thin layers of water
610    ##protect_against_infinitesimal_and_negative_heights(domain)
611
612    #Extrapolate all conserved quantities
613    #if optimised_gradient_limiter:
614    #    #MH090605 if second order,
615    #    #perform the extrapolation and limiting on
616    #    #all of the conserved quantities
617
618    #    if (domain.order == 1):
619    #        for name in domain.conserved_quantities:
620    #            Q = domain.quantities[name]
621    #            Q.extrapolate_first_order()
622    #    elif domain.order == 2:
623    #        domain.extrapolate_second_order_sw()
624    #    else:
625    #        raise 'Unknown order'
626    #else:
627        #old code:
628    for name in domain.conserved_quantities:
629        Q = domain.quantities[name]
630        if domain.order == 1:
631            Q.extrapolate_first_order()
632        elif domain.order == 2:
633            Q.extrapolate_second_order()
634            #Q.limit()
635        else:
636            raise 'Unknown order'
637
638    #Take bed elevation into account when water heights are small
639    #balance_deep_and_shallow(domain)
640
641    #Compute edge values by interpolation
642    #for name in domain.conserved_quantities:
643    #    Q = domain.quantities[name]
644    #    Q.interpolate_from_vertices_to_edges()
645
646
647def protect_against_infinitesimal_and_negative_heights(domain):
648    """Protect against infinitesimal heights and associated high velocities
649    """
650
651    #Shortcuts
652    wc = domain.quantities['stage'].centroid_values
653    zc = domain.quantities['elevation'].centroid_values
654    xmomc = domain.quantities['xmomentum'].centroid_values
655#    ymomc = domain.quantities['ymomentum'].centroid_values
656    hc = wc - zc  #Water depths at centroids
657
658    #Update
659    for k in range(domain.number_of_elements):
660
661        if hc[k] < domain.minimum_allowed_height:
662            #Control stage
663            if hc[k] < domain.epsilon:
664                wc[k] = zc[k] # Contain 'lost mass' error
665
666            #Control momentum
667            #xmomc[k] = ymomc[k] = 0.0
668            xmomc[k] = 0.0
669
670def protect_against_infinitesimal_and_negative_heights_c(domain):
671    """Protect against infinitesimal heights and associated high velocities
672    """
673
674    #Shortcuts
675    wc = domain.quantities['stage'].centroid_values
676    zc = domain.quantities['elevation'].centroid_values
677    xmomc = domain.quantities['xmomentum'].centroid_values
678#    ymomc = domain.quantities['ymomentum'].centroid_values
679
680    from shallow_water_ext import protect
681
682    protect(domain.minimum_allowed_height, domain.epsilon,
683            wc, zc, xmomc)#, ymomc)
684
685def h_limiter(domain):
686    """Limit slopes for each volume to eliminate artificial variance
687    introduced by e.g. second order extrapolator
688
689    limit on h = w-z
690
691    This limiter depends on two quantities (w,z) so it resides within
692    this module rather than within quantity.py
693    """
694
695    from Numeric import zeros, Float
696
697    N = domain.number_of_elements
698    beta_h = domain.beta_h
699
700    #Shortcuts
701    wc = domain.quantities['stage'].centroid_values
702    zc = domain.quantities['elevation'].centroid_values
703    hc = wc - zc
704
705    wv = domain.quantities['stage'].vertex_values
706    zv = domain.quantities['elevation'].vertex_values
707    hv = wv-zv
708
709    hvbar = zeros(hv.shape, Float) #h-limited values
710
711    #Find min and max of this and neighbour's centroid values
712    hmax = zeros(hc.shape, Float)
713    hmin = zeros(hc.shape, Float)
714
715    for k in range(N):
716        hmax[k] = hmin[k] = hc[k]
717        #for i in range(3):
718        for i in range(2):   
719            n = domain.neighbours[k,i]
720            if n >= 0:
721                hn = hc[n] #Neighbour's centroid value
722
723                hmin[k] = min(hmin[k], hn)
724                hmax[k] = max(hmax[k], hn)
725
726
727    #Diffences between centroids and maxima/minima
728    dhmax = hmax - hc
729    dhmin = hmin - hc
730
731    #Deltas between vertex and centroid values
732    dh = zeros(hv.shape, Float)
733    #for i in range(3):
734    for i in range(2):
735        dh[:,i] = hv[:,i] - hc
736
737    #Phi limiter
738    for k in range(N):
739
740        #Find the gradient limiter (phi) across vertices
741        phi = 1.0
742        #for i in range(3):
743        for i in range(2):
744            r = 1.0
745            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
746            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
747
748            phi = min( min(r*beta_h, 1), phi )
749
750        #Then update using phi limiter
751        #for i in range(3):
752        for i in range(2):
753            hvbar[k,i] = hc[k] + phi*dh[k,i]
754
755    return hvbar
756
757
758
759def h_limiter_c(domain):
760    """Limit slopes for each volume to eliminate artificial variance
761    introduced by e.g. second order extrapolator
762
763    limit on h = w-z
764
765    This limiter depends on two quantities (w,z) so it resides within
766    this module rather than within quantity.py
767
768    Wrapper for c-extension
769    """
770
771    from Numeric import zeros, Float
772
773    N = domain.number_of_elements
774    beta_h = domain.beta_h
775
776    #Shortcuts
777    wc = domain.quantities['stage'].centroid_values
778    zc = domain.quantities['elevation'].centroid_values
779    hc = wc - zc
780
781    wv = domain.quantities['stage'].vertex_values
782    zv = domain.quantities['elevation'].vertex_values
783    hv = wv - zv
784
785    #Call C-extension
786    from shallow_water_ext import h_limiter_sw as h_limiter
787    hvbar = h_limiter(domain, hc, hv)
788
789    return hvbar
790
791def balance_deep_and_shallow(domain):
792    """Compute linear combination between stage as computed by
793    gradient-limiters limiting using w, and stage computed by
794    gradient-limiters limiting using h (h-limiter).
795    The former takes precedence when heights are large compared to the
796    bed slope while the latter takes precedence when heights are
797    relatively small.  Anything in between is computed as a balanced
798    linear combination in order to avoid numerical disturbances which
799    would otherwise appear as a result of hard switching between
800    modes.
801
802    The h-limiter is always applied irrespective of the order.
803    """
804
805    #Shortcuts
806    wc = domain.quantities['stage'].centroid_values
807    zc = domain.quantities['elevation'].centroid_values
808    hc = wc - zc
809
810    wv = domain.quantities['stage'].vertex_values
811    zv = domain.quantities['elevation'].vertex_values
812    hv = wv-zv
813
814    #Limit h
815    hvbar = h_limiter(domain)
816
817    for k in range(domain.number_of_elements):
818        #Compute maximal variation in bed elevation
819        #  This quantitiy is
820        #    dz = max_i abs(z_i - z_c)
821        #  and it is independent of dimension
822        #  In the 1d case zc = (z0+z1)/2
823        #  In the 2d case zc = (z0+z1+z2)/3
824
825        dz = max(abs(zv[k,0]-zc[k]),
826                 abs(zv[k,1]-zc[k]))#,
827#                 abs(zv[k,2]-zc[k]))
828
829
830        hmin = min( hv[k,:] )
831
832        #Create alpha in [0,1], where alpha==0 means using the h-limited
833        #stage and alpha==1 means using the w-limited stage as
834        #computed by the gradient limiter (both 1st or 2nd order)
835
836        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
837        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
838
839        if dz > 0.0:
840            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
841        else:
842            #Flat bed
843            alpha = 1.0
844
845        #Let
846        #
847        #  wvi be the w-limited stage (wvi = zvi + hvi)
848        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
849        #
850        #
851        #where i=0,1,2 denotes the vertex ids
852        #
853        #Weighted balance between w-limited and h-limited stage is
854        #
855        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
856        #
857        #It follows that the updated wvi is
858        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
859        #
860        # Momentum is balanced between constant and limited
861
862
863        #for i in range(3):
864        #    wv[k,i] = zv[k,i] + hvbar[k,i]
865
866        #return
867
868        if alpha < 1:
869
870            #for i in range(3):
871            for i in range(2):
872                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
873
874            #Momentums at centroids
875            xmomc = domain.quantities['xmomentum'].centroid_values
876#            ymomc = domain.quantities['ymomentum'].centroid_values
877
878            #Momentums at vertices
879            xmomv = domain.quantities['xmomentum'].vertex_values
880#           ymomv = domain.quantities['ymomentum'].vertex_values
881
882            # Update momentum as a linear combination of
883            # xmomc and ymomc (shallow) and momentum
884            # from extrapolator xmomv and ymomv (deep).
885            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
886#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
887
888
889def balance_deep_and_shallow_c(domain):
890    """Wrapper for C implementation
891    """
892
893    #Shortcuts
894    wc = domain.quantities['stage'].centroid_values
895    zc = domain.quantities['elevation'].centroid_values
896    hc = wc - zc
897
898    wv = domain.quantities['stage'].vertex_values
899    zv = domain.quantities['elevation'].vertex_values
900    hv = wv - zv
901
902    #Momentums at centroids
903    xmomc = domain.quantities['xmomentum'].centroid_values
904 #   ymomc = domain.quantities['ymomentum'].centroid_values
905
906    #Momentums at vertices
907    xmomv = domain.quantities['xmomentum'].vertex_values
908#    ymomv = domain.quantities['ymomentum'].vertex_values
909
910    #Limit h
911    hvbar = h_limiter(domain)
912
913    #This is how one would make a first order h_limited value
914    #as in the old balancer (pre 17 Feb 2005):
915    #from Numeric import zeros, Float
916    #hvbar = zeros( (len(hc), 3), Float)
917    #for i in range(3):
918    #    hvbar[:,i] = hc[:]
919
920    from shallow_water_ext import balance_deep_and_shallow
921    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
922    #                         xmomc, ymomc, xmomv, ymomv)
923                             xmomc, xmomv)
924
925
926###############################################
927#Boundaries - specific to the shallow water wave equation
928class Reflective_boundary(Boundary):
929    """Reflective boundary returns same conserved quantities as
930    those present in its neighbour volume but reflected.
931
932    This class is specific to the shallow water equation as it
933    works with the momentum quantities assumed to be the second
934    and third conserved quantities.
935    """
936
937    def __init__(self, domain = None):
938        Boundary.__init__(self)
939
940        if domain is None:
941            msg = 'Domain must be specified for reflective boundary'
942            raise msg
943
944        #Handy shorthands
945        #self.stage   = domain.quantities['stage'].edge_values
946        #self.xmom    = domain.quantities['xmomentum'].edge_values
947        #self.ymom    = domain.quantities['ymomentum'].edge_values
948        #self.normals = domain.normals
949        self.stage   = domain.quantities['stage'].vertex_values
950        self.xmom    = domain.quantities['xmomentum'].vertex_values
951
952        from Numeric import zeros, Float
953        #self.conserved_quantities = zeros(3, Float)
954        self.conserved_quantities = zeros(2, Float)
955
956    def __repr__(self):
957        return 'Reflective_boundary'
958
959
960    def evaluate(self, vol_id, edge_id):
961        """Reflective boundaries reverses the outward momentum
962        of the volume they serve.
963        """
964
965        q = self.conserved_quantities
966        q[0] = self.stage[vol_id, edge_id]
967        q[1] = self.xmom[vol_id, edge_id]
968        #q[2] = self.ymom[vol_id, edge_id]
969
970        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
971        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
972
973
974        #r = rotate(q, normal, direction = 1)
975        #r[1] = -r[1]
976        #q = rotate(r, normal, direction = -1)
977        r = q
978        r[1] = -q[1]
979        q = r
980        #For start interval there is no outward momentum so do not need to
981        #reverse direction in this case
982
983        return q
984
985
986#########################
987#Standard forcing terms:
988#
989def gravity(domain):
990    """Apply gravitational pull in the presence of bed slope
991    """
992
993    from anuga.pyvolution.util import gradient
994    from Numeric import zeros, Float, array, sum
995
996    xmom = domain.quantities['xmomentum'].explicit_update
997#    ymom = domain.quantities['ymomentum'].explicit_update
998
999    Stage = domain.quantities['stage']
1000    Elevation = domain.quantities['elevation']
1001    #h = Stage.edge_values - Elevation.edge_values
1002    h = Stage.vertex_values - Elevation.vertex_values
1003    v = Elevation.vertex_values
1004
1005    x = domain.get_vertex_coordinates()
1006    g = domain.g
1007
1008    for k in range(domain.number_of_elements):
1009#        avg_h = sum( h[k,:] )/3
1010        avg_h = sum( h[k,:] )/2
1011
1012        #Compute bed slope
1013        #x0, y0, x1, y1, x2, y2 = x[k,:]
1014        x0, x1 = x[k,:]
1015        #z0, z1, z2 = v[k,:]
1016        z0, z1 = v[k,:]
1017
1018        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
1019        zx = gradient(x0, x1, z0, z1)
1020       
1021        #Update momentum
1022        xmom[k] += -g*zx*avg_h
1023#        ymom[k] += -g*zy*avg_h
1024
1025
1026def graity_c(domain):
1027    """Wrapper calling C version
1028    """
1029
1030    xmom = domain.quantities['xmomentum'].explicit_update
1031#   ymom = domain.quantities['ymomentum'].explicit_update
1032
1033    Stage = domain.quantities['stage']
1034    Elevation = domain.quantities['elevation']
1035    h = Stage.edge_values - Elevation.edge_values
1036    v = Elevation.vertex_values
1037
1038    x = domain.get_vertex_coordinates()
1039    g = domain.g
1040
1041
1042    from shallow_water_ext import gravity
1043    gravity(g, h, v, x, xmom, ymom)
1044
1045
1046def manning_friction(domain):
1047    """Apply (Manning) friction to water momentum
1048    """
1049
1050    from math import sqrt
1051
1052    w = domain.quantities['stage'].centroid_values
1053    z = domain.quantities['elevation'].centroid_values
1054    h = w-z
1055
1056    uh = domain.quantities['xmomentum'].centroid_values
1057    #vh = domain.quantities['ymomentum'].centroid_values
1058    eta = domain.quantities['friction'].centroid_values
1059
1060    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1061    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1062
1063    N = domain.number_of_elements
1064    eps = domain.minimum_allowed_height
1065    g = domain.g
1066
1067    for k in range(N):
1068        if eta[k] >= eps:
1069            if h[k] >= eps:
1070                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1071                S = -g * eta[k]**2 * uh[k]
1072                S /= h[k]**(7.0/3)
1073
1074                #Update momentum
1075                xmom_update[k] += S*uh[k]
1076                #ymom_update[k] += S*vh[k]
1077
1078
1079def manning_friction_c(domain):
1080    """Wrapper for c version
1081    """
1082
1083    xmom = domain.quantities['xmomentum']
1084#    ymom = domain.quantities['ymomentum']
1085
1086    w = domain.quantities['stage'].centroid_values
1087    z = domain.quantities['elevation'].centroid_values
1088
1089    uh = xmom.centroid_values
1090#    vh = ymom.centroid_values
1091    eta = domain.quantities['friction'].centroid_values
1092
1093    xmom_update = xmom.semi_implicit_update
1094#    ymom_update = ymom.semi_implicit_update
1095
1096    N = domain.number_of_elements
1097    eps = domain.minimum_allowed_height
1098    g = domain.g
1099
1100    from shallow_water_ext import manning_friction
1101#    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
1102    manning_friction(g, eps, w, z, uh, eta, xmom_update)
1103
1104def linear_friction(domain):
1105    """Apply linear friction to water momentum
1106
1107    Assumes quantity: 'linear_friction' to be present
1108    """
1109
1110    from math import sqrt
1111
1112    w = domain.quantities['stage'].centroid_values
1113    z = domain.quantities['elevation'].centroid_values
1114    h = w-z
1115
1116    uh = domain.quantities['xmomentum'].centroid_values
1117#    vh = domain.quantities['ymomentum'].centroid_values
1118    tau = domain.quantities['linear_friction'].centroid_values
1119
1120    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1121#    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1122
1123    N = domain.number_of_elements
1124    eps = domain.minimum_allowed_height
1125    g = domain.g #Not necessary? Why was this added?
1126
1127    for k in range(N):
1128        if tau[k] >= eps:
1129            if h[k] >= eps:
1130                S = -tau[k]/h[k]
1131
1132                #Update momentum
1133                xmom_update[k] += S*uh[k]
1134 #              ymom_update[k] += S*vh[k]
1135
1136
1137
1138def check_forcefield(f):
1139    """Check that f is either
1140    1: a callable object f(t,x,y), where x and y are vectors
1141       and that it returns an array or a list of same length
1142       as x and y
1143    2: a scalar
1144    """
1145
1146    from Numeric import ones, Float, array
1147
1148
1149    if callable(f):
1150        #N = 3
1151        N = 2
1152        #x = ones(3, Float)
1153        #y = ones(3, Float)
1154        x = ones(2, Float)
1155        #y = ones(2, Float)
1156       
1157        try:
1158            #q = f(1.0, x=x, y=y)
1159            q = f(1.0, x=x)
1160        except Exception, e:
1161            msg = 'Function %s could not be executed:\n%s' %(f, e)
1162            #FIXME: Reconsider this semantics
1163            raise msg
1164
1165        try:
1166            q = array(q).astype(Float)
1167        except:
1168            msg = 'Return value from vector function %s could ' %f
1169            msg += 'not be converted into a Numeric array of floats.\n'
1170            msg += 'Specified function should return either list or array.'
1171            raise msg
1172
1173        #Is this really what we want?
1174        msg = 'Return vector from function %s ' %f
1175        msg += 'must have same lenght as input vectors'
1176        assert len(q) == N, msg
1177
1178    else:
1179        try:
1180            f = float(f)
1181        except:
1182            msg = 'Force field %s must be either a scalar' %f
1183            msg += ' or a vector function'
1184            raise msg
1185    return f
1186
1187class Wind_stress:
1188    """Apply wind stress to water momentum in terms of
1189    wind speed [m/s] and wind direction [degrees]
1190    """
1191
1192    def __init__(self, *args, **kwargs):
1193        """Initialise windfield from wind speed s [m/s]
1194        and wind direction phi [degrees]
1195
1196        Inputs v and phi can be either scalars or Python functions, e.g.
1197
1198        W = Wind_stress(10, 178)
1199
1200        #FIXME - 'normal' degrees are assumed for now, i.e. the
1201        vector (1,0) has zero degrees.
1202        We may need to convert from 'compass' degrees later on and also
1203        map from True north to grid north.
1204
1205        Arguments can also be Python functions of t,x,y as in
1206
1207        def speed(t,x,y):
1208            ...
1209            return s
1210
1211        def angle(t,x,y):
1212            ...
1213            return phi
1214
1215        where x and y are vectors.
1216
1217        and then pass the functions in
1218
1219        W = Wind_stress(speed, angle)
1220
1221        The instantiated object W can be appended to the list of
1222        forcing_terms as in
1223
1224        Alternatively, one vector valued function for (speed, angle)
1225        can be applied, providing both quantities simultaneously.
1226        As in
1227        W = Wind_stress(F), where returns (speed, angle) for each t.
1228
1229        domain.forcing_terms.append(W)
1230        """
1231
1232        from config import rho_a, rho_w, eta_w
1233        from Numeric import array, Float
1234
1235        if len(args) == 2:
1236            s = args[0]
1237            phi = args[1]
1238        elif len(args) == 1:
1239            #Assume vector function returning (s, phi)(t,x,y)
1240            vector_function = args[0]
1241            #s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1242            #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1243            s = lambda t,x: vector_function(t,x=x)[0]
1244            phi = lambda t,x: vector_function(t,x=x)[1]
1245        else:
1246           #Assume info is in 2 keyword arguments
1247
1248           if len(kwargs) == 2:
1249               s = kwargs['s']
1250               phi = kwargs['phi']
1251           else:
1252               raise 'Assumes two keyword arguments: s=..., phi=....'
1253
1254        print 'phi', phi
1255        self.speed = check_forcefield(s)
1256        self.phi = check_forcefield(phi)
1257
1258        self.const = eta_w*rho_a/rho_w
1259
1260
1261    def __call__(self, domain):
1262        """Evaluate windfield based on values found in domain
1263        """
1264
1265        from math import pi, cos, sin, sqrt
1266        from Numeric import Float, ones, ArrayType
1267
1268        xmom_update = domain.quantities['xmomentum'].explicit_update
1269        #ymom_update = domain.quantities['ymomentum'].explicit_update
1270
1271        N = domain.number_of_elements
1272        t = domain.time
1273
1274        if callable(self.speed):
1275            xc = domain.get_centroid_coordinates()
1276            #s_vec = self.speed(t, xc[:,0], xc[:,1])
1277            s_vec = self.speed(t, xc)
1278        else:
1279            #Assume s is a scalar
1280
1281            try:
1282                s_vec = self.speed * ones(N, Float)
1283            except:
1284                msg = 'Speed must be either callable or a scalar: %s' %self.s
1285                raise msg
1286
1287
1288        if callable(self.phi):
1289            xc = domain.get_centroid_coordinates()
1290            #phi_vec = self.phi(t, xc[:,0], xc[:,1])
1291            phi_vec = self.phi(t, xc)
1292        else:
1293            #Assume phi is a scalar
1294
1295            try:
1296                phi_vec = self.phi * ones(N, Float)
1297            except:
1298                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1299                raise msg
1300
1301        #assign_windfield_values(xmom_update, ymom_update,
1302        #                        s_vec, phi_vec, self.const)
1303        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
1304
1305
1306#def assign_windfield_values(xmom_update, ymom_update,
1307#                            s_vec, phi_vec, const):
1308def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
1309    """Python version of assigning wind field to update vectors.
1310    A c version also exists (for speed)
1311    """
1312    from math import pi, cos, sin, sqrt
1313
1314    N = len(s_vec)
1315    for k in range(N):
1316        s = s_vec[k]
1317        phi = phi_vec[k]
1318
1319        #Convert to radians
1320        phi = phi*pi/180
1321
1322        #Compute velocity vector (u, v)
1323        u = s*cos(phi)
1324        v = s*sin(phi)
1325
1326        #Compute wind stress
1327        #S = const * sqrt(u**2 + v**2)
1328        S = const * u
1329        xmom_update[k] += S*u
1330        #ymom_update[k] += S*v
Note: See TracBrowser for help on using the repository browser.