source: anuga_work/development/anuga_1d/anuga_1d_suggestion1_2_3/shallow_water_domain_suggestion1.py @ 7576

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