source: anuga_work/development/anuga_1d/shallow_water_domain_steve.py @ 7818

Last change on this file since 7818 was 6453, checked in by steve, 15 years ago

Added Padarn's (2008/2009 summer student) files

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