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

Last change on this file since 5587 was 5587, checked in by steve, 16 years ago

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