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

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

Working version of comp_flux_ext

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