source: anuga_work/development/anuga_1d/shallow_water_domain.py @ 6038

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