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

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