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

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

Updated timestepping
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
46from domain import *
47Generic_Domain = Domain #Rename
48
49#Shallow water domain
50class Domain(Generic_Domain):
51
52    def __init__(self, coordinates, boundary = None, tagged_elements = None):
53
54        conserved_quantities = ['stage', 'xmomentum']
55        other_quantities = ['elevation', 'friction']
56        Generic_Domain.__init__(self, coordinates, boundary,
57                                conserved_quantities, other_quantities,
58                                tagged_elements)
59       
60        from config import minimum_allowed_height, g, h0
61        self.minimum_allowed_height = minimum_allowed_height
62        self.g = g
63        self.h0 = h0
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        self.__doc__ = 'shallow_water_domain'
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, duration = 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,duration,
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, h0
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 +  h0/h_left)
336
337
338    uh_left = u_left*h_left
339
340
341    w_right  = q_right[0]  #w=h+z
342    h_right  = w_right-z
343    uh_right = q_right[1]
344
345    if h_right < epsilon:
346        u_right = 0.0  #Could have been negative
347        h_right = 0.0
348    else:
349        u_right  = uh_right/(h_right + h0/h_right)
350
351    uh_right = u_right*h_right
352
353    #vh_left  = q_left[2]
354    #vh_right = q_right[2]
355
356    #print h_right
357    #print u_right
358    #print h_left
359    #print u_right
360
361    soundspeed_left  = sqrt(g*h_left)
362    soundspeed_right = sqrt(g*h_right)
363
364    #Maximal wave speed
365    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
366   
367    #Minimal wave speed
368    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
369   
370    #Flux computation
371
372    #flux_left  = array([u_left*h_left,
373    #                    u_left*uh_left + 0.5*g*h_left**2])
374    #flux_right = array([u_right*h_right,
375    #                    u_right*uh_right + 0.5*g*h_right**2])
376    flux_left  = array([u_left*h_left,
377                        u_left*uh_left + 0.5*g*h_left*h_left])
378    flux_right = array([u_right*h_right,
379                        u_right*uh_right + 0.5*g*h_right*h_right])
380
381    denom = s_max-s_min
382    if denom == 0.0:
383        edgeflux = array([0.0, 0.0])
384        max_speed = 0.0
385    else:
386        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
387        edgeflux += s_max*s_min*(q_right-q_left)/denom
388       
389        edgeflux[1] = edgeflux[1]*normal
390
391        max_speed = max(abs(s_max), abs(s_min))
392
393    return edgeflux, max_speed
394
395
396 
397
398def compute_timestep(domain):
399    import sys
400    from Numeric import zeros, Float
401
402    N = domain.number_of_elements
403
404    #Shortcuts
405    Stage = domain.quantities['stage']
406    Xmom = domain.quantities['xmomentum']
407    Bed = domain.quantities['elevation']
408   
409    stage = Stage.vertex_values
410    xmom =  Xmom.vertex_values
411    bed =   Bed.vertex_values
412
413    stage_bdry = Stage.boundary_values
414    xmom_bdry =  Xmom.boundary_values
415
416    flux = zeros(2, Float) #Work array for summing up fluxes
417    ql = zeros(2, Float)
418    qr = zeros(2, Float)
419
420    #Loop
421    timestep = float(sys.maxint)
422    enter = True
423    for k in range(N):
424
425        flux[:] = 0.  #Reset work array
426        for i in range(2):
427            #Quantities inside volume facing neighbour i
428            ql = [stage[k, i], xmom[k, i]]
429            zl = bed[k, i]
430
431            #Quantities at neighbour on nearest face
432            n = domain.neighbours[k,i]
433            if n < 0:
434                m = -n-1 #Convert negative flag to index
435                qr[0] = stage_bdry[m]
436                qr[1] = xmom_bdry[m]
437                zr = zl #Extend bed elevation to boundary
438            else:
439                #m = domain.neighbour_edges[k,i]
440                m = domain.neighbour_vertices[k,i]
441                qr[0] = stage[n, m]
442                qr[1] =  xmom[n, m]
443                zr = bed[n, m]
444
445
446            #Outward pointing normal vector
447            normal = domain.normals[k, i]
448
449            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
450
451            #Update optimal_timestep
452            try:
453                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
454            except ZeroDivisionError:
455                pass
456
457    domain.timestep = timestep
458
459def compute_fluxes(domain):
460    """Compute all fluxes and the timestep suitable for all volumes
461    in domain.
462
463    Compute total flux for each conserved quantity using "flux_function"
464
465    Fluxes across each edge are scaled by edgelengths and summed up
466    Resulting flux is then scaled by area and stored in
467    explicit_update for each of the three conserved quantities
468    stage, xmomentum and ymomentum
469
470    The maximal allowable speed computed by the flux_function for each volume
471    is converted to a timestep that must not be exceeded. The minimum of
472    those is computed as the next overall timestep.
473
474    Post conditions:
475      domain.explicit_update is reset to computed flux values
476      domain.timestep is set to the largest step satisfying all volumes.
477    """
478
479    import sys
480    from Numeric import zeros, Float
481
482    N = domain.number_of_elements
483
484    #Shortcuts
485    Stage = domain.quantities['stage']
486    Xmom = domain.quantities['xmomentum']
487#    Ymom = domain.quantities['ymomentum']
488    Bed = domain.quantities['elevation']
489
490    #Arrays
491    #stage = Stage.edge_values
492    #xmom =  Xmom.edge_values
493 #   ymom =  Ymom.edge_values
494    #bed =   Bed.edge_values
495   
496    stage = Stage.vertex_values
497    xmom =  Xmom.vertex_values
498    bed =   Bed.vertex_values
499   
500    #print 'stage edge values', stage
501    #print 'xmom edge values', xmom
502    #print 'bed values', bed
503
504    stage_bdry = Stage.boundary_values
505    xmom_bdry =  Xmom.boundary_values
506    #print 'stage_bdry',stage_bdry
507    #print 'xmom_bdry', xmom_bdry
508#    ymom_bdry =  Ymom.boundary_values
509
510#    flux = zeros(3, Float) #Work array for summing up fluxes
511    flux = zeros(2, Float) #Work array for summing up fluxes
512    ql = zeros(2, Float)
513    qr = zeros(2, Float)
514
515    #Loop
516    timestep = float(sys.maxint)
517    enter = True
518    for k in range(N):
519
520        flux[:] = 0.  #Reset work array
521        #for i in range(3):
522        for i in range(2):
523            #Quantities inside volume facing neighbour i
524            #ql[0] = stage[k, i]
525            #ql[1] = xmom[k, i]
526            ql = [stage[k, i], xmom[k, i]]
527            zl = bed[k, i]
528
529            #Quantities at neighbour on nearest face
530            n = domain.neighbours[k,i]
531            if n < 0:
532                m = -n-1 #Convert negative flag to index
533                qr[0] = stage_bdry[m]
534                qr[1] = xmom_bdry[m]
535                zr = zl #Extend bed elevation to boundary
536            else:
537                #m = domain.neighbour_edges[k,i]
538                m = domain.neighbour_vertices[k,i]
539                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
540                qr[0] = stage[n, m]
541                qr[1] = xmom[n, m]
542                zr = bed[n, m]
543
544
545            #Outward pointing normal vector
546            normal = domain.normals[k, i]
547       
548            #Flux computation using provided function
549            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
550            #print 'ql',ql
551            #print 'qr',qr
552           
553
554            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
555
556            #print 'edgeflux', edgeflux
557
558            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
559            # flux = edgefluxleft - edgefluxright
560            flux -= edgeflux #* domain.edgelengths[k,i]
561            #Update optimal_timestep
562            try:
563                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
564                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
565            except ZeroDivisionError:
566                pass
567
568        #Normalise by area and store for when all conserved
569        #quantities get updated
570        flux /= domain.areas[k]
571
572        Stage.explicit_update[k] = flux[0]
573        Xmom.explicit_update[k] =  flux[1]
574        #Ymom.explicit_update[k] = flux[2]
575        #print "flux cell",k,flux[0]
576
577    domain.flux_timestep = timestep
578    #print domain.quantities['stage'].centroid_values
579
580
581# Compute flux definition
582def compute_fluxes_C_long(domain):
583        from Numeric import zeros, Float
584        import sys
585
586       
587        timestep = float(sys.maxint)
588        #print 'timestep=',timestep
589        #print 'The type of timestep is',type(timestep)
590       
591        epsilon = domain.epsilon
592        #print 'epsilon=',epsilon
593        #print 'The type of epsilon is',type(epsilon)
594       
595        g = domain.g
596        #print 'g=',g
597        #print 'The type of g is',type(g)
598       
599        neighbours = domain.neighbours
600        #print 'neighbours=',neighbours
601        #print 'The type of neighbours is',type(neighbours)
602       
603        neighbour_vertices = domain.neighbour_vertices
604        #print 'neighbour_vertices=',neighbour_vertices
605        #print 'The type of neighbour_vertices is',type(neighbour_vertices)
606       
607        normals = domain.normals
608        #print 'normals=',normals
609        #print 'The type of normals is',type(normals)
610       
611        areas = domain.areas
612        #print 'areas=',areas
613        #print 'The type of areas is',type(areas)
614       
615        stage_edge_values = domain.quantities['stage'].vertex_values
616        #print 'stage_edge_values=',stage_edge_values
617        #print 'The type of stage_edge_values is',type(stage_edge_values)
618       
619        xmom_edge_values = domain.quantities['xmomentum'].vertex_values
620        #print 'xmom_edge_values=',xmom_edge_values
621        #print 'The type of xmom_edge_values is',type(xmom_edge_values)
622       
623        bed_edge_values = domain.quantities['elevation'].vertex_values
624        #print 'bed_edge_values=',bed_edge_values
625        #print 'The type of bed_edge_values is',type(bed_edge_values)
626       
627        stage_boundary_values = domain.quantities['stage'].boundary_values
628        #print 'stage_boundary_values=',stage_boundary_values
629        #print 'The type of stage_boundary_values is',type(stage_boundary_values)
630       
631        xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
632        #print 'xmom_boundary_values=',xmom_boundary_values
633        #print 'The type of xmom_boundary_values is',type(xmom_boundary_values)
634       
635        stage_explicit_update = domain.quantities['stage'].explicit_update
636        #print 'stage_explicit_update=',stage_explicit_update
637        #print 'The type of stage_explicit_update is',type(stage_explicit_update)
638       
639        xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
640        #print 'xmom_explicit_update=',xmom_explicit_update
641        #print 'The type of xmom_explicit_update is',type(xmom_explicit_update)
642       
643        number_of_elements = len(stage_edge_values)
644        #print 'number_of_elements=',number_of_elements
645        #print 'The type of number_of_elements is',type(number_of_elements)
646       
647        max_speed_array = domain.max_speed_array
648        #print 'max_speed_array=',max_speed_array
649        #print 'The type of max_speed_array is',type(max_speed_array)
650       
651       
652        from comp_flux_ext import compute_fluxes_ext
653       
654        domain.flux_timestep = compute_fluxes_ext(timestep,
655                                  epsilon,
656                                  g,
657                                  neighbours,
658                                  neighbour_vertices,
659                                  normals,
660                                  areas, 
661                                  stage_edge_values,
662                                  xmom_edge_values,
663                                  bed_edge_values,
664                                  stage_boundary_values,
665                                  xmom_boundary_values,
666                                  stage_explicit_update,
667                                  xmom_explicit_update,
668                                  number_of_elements,
669                                  max_speed_array)
670
671
672# Compute flux definition
673def compute_fluxes_C_short(domain):
674        from Numeric import zeros, Float
675        import sys
676
677       
678        timestep = float(sys.maxint)
679
680        stage = domain.quantities['stage']
681        xmom = domain.quantities['xmomentum']
682        bed = domain.quantities['elevation']
683
684
685        from comp_flux_ext import compute_fluxes_ext_short
686
687        domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
688
689
690
691####################################
692
693
694
695
696
697
698# Module functions for gradient limiting (distribute_to_vertices_and_edges)
699
700def distribute_to_vertices_and_edges(domain):
701    """Distribution from centroids to vertices specific to the
702    shallow water wave
703    equation.
704
705    It will ensure that h (w-z) is always non-negative even in the
706    presence of steep bed-slopes by taking a weighted average between shallow
707    and deep cases.
708
709    In addition, all conserved quantities get distributed as per either a
710    constant (order==1) or a piecewise linear function (order==2).
711
712    FIXME: more explanation about removal of artificial variability etc
713
714    Precondition:
715      All quantities defined at centroids and bed elevation defined at
716      vertices.
717
718    Postcondition
719      Conserved quantities defined at vertices
720
721    """
722
723    #from config import optimised_gradient_limiter
724
725    #Remove very thin layers of water
726    protect_against_infinitesimal_and_negative_heights(domain)
727   
728
729    #Extrapolate all conserved quantities
730    #if optimised_gradient_limiter:
731    #    #MH090605 if second order,
732    #    #perform the extrapolation and limiting on
733    #    #all of the conserved quantities
734
735    #    if (domain.order == 1):
736    #        for name in domain.conserved_quantities:
737    #            Q = domain.quantities[name]
738    #            Q.extrapolate_first_order()
739    #    elif domain.order == 2:
740    #        domain.extrapolate_second_order_sw()
741    #    else:
742    #        raise 'Unknown order'
743    #else:
744        #old code:
745
746    for name in domain.conserved_quantities:
747        Q = domain.quantities[name]
748        if domain.order == 1:
749            Q.extrapolate_first_order()
750        elif domain.order == 2:
751            #print "add extrapolate second order to shallow water"
752            #if name != 'height':
753            Q.extrapolate_second_order()
754            #Q.limit()
755        else:
756            raise 'Unknown order'
757
758    #Take bed elevation into account when water heights are small
759    #balance_deep_and_shallow(domain)
760    #protect_against_infinitesimal_and_negative_heights(domain)
761
762#Compute edge values by interpolation
763    #for name in domain.conserved_quantities:
764    #    Q = domain.quantities[name]
765    #    Q.interpolate_from_vertices_to_edges()   
766   
767def protect_against_infinitesimal_and_negative_heights(domain):
768    """Protect against infinitesimal heights and associated high velocities
769    """
770
771    #Shortcuts
772    wc = domain.quantities['stage'].centroid_values
773    zc = domain.quantities['elevation'].centroid_values
774    xmomc = domain.quantities['xmomentum'].centroid_values
775#    ymomc = domain.quantities['ymomentum'].centroid_values
776    hc = wc - zc  #Water depths at centroids
777
778    zv = domain.quantities['elevation'].vertex_values
779    wv = domain.quantities['stage'].vertex_values
780    hv = wv-zv
781    xmomv = domain.quantities['xmomentum'].vertex_values
782    #remove the above two lines and corresponding code below
783
784    #Update
785    for k in range(domain.number_of_elements):
786
787        if hc[k] < domain.minimum_allowed_height:
788            #Control stage
789            if hc[k] < domain.epsilon:
790                wc[k] = zc[k] # Contain 'lost mass' error
791                wv[k,0] = zv[k,0]
792                wv[k,1] = zv[k,1]
793               
794            xmomc[k] = 0.0
795           
796        #N = domain.number_of_elements
797        #if (k == 0) | (k==N-1):
798        #    wc[k] = zc[k] # Contain 'lost mass' error
799        #    wv[k,0] = zv[k,0]
800        #    wv[k,1] = zv[k,1]
801   
802def h_limiter(domain):
803    """Limit slopes for each volume to eliminate artificial variance
804    introduced by e.g. second order extrapolator
805
806    limit on h = w-z
807
808    This limiter depends on two quantities (w,z) so it resides within
809    this module rather than within quantity.py
810    """
811
812    from Numeric import zeros, Float
813
814    N = domain.number_of_elements
815    beta_h = domain.beta_h
816
817    #Shortcuts
818    wc = domain.quantities['stage'].centroid_values
819    zc = domain.quantities['elevation'].centroid_values
820    hc = wc - zc
821
822    wv = domain.quantities['stage'].vertex_values
823    zv = domain.quantities['elevation'].vertex_values
824    hv = wv-zv
825
826    hvbar = zeros(hv.shape, Float) #h-limited values
827
828    #Find min and max of this and neighbour's centroid values
829    hmax = zeros(hc.shape, Float)
830    hmin = zeros(hc.shape, Float)
831
832    for k in range(N):
833        hmax[k] = hmin[k] = hc[k]
834        #for i in range(3):
835        for i in range(2):   
836            n = domain.neighbours[k,i]
837            if n >= 0:
838                hn = hc[n] #Neighbour's centroid value
839
840                hmin[k] = min(hmin[k], hn)
841                hmax[k] = max(hmax[k], hn)
842
843
844    #Diffences between centroids and maxima/minima
845    dhmax = hmax - hc
846    dhmin = hmin - hc
847
848    #Deltas between vertex and centroid values
849    dh = zeros(hv.shape, Float)
850    #for i in range(3):
851    for i in range(2):
852        dh[:,i] = hv[:,i] - hc
853
854    #Phi limiter
855    for k in range(N):
856
857        #Find the gradient limiter (phi) across vertices
858        phi = 1.0
859        #for i in range(3):
860        for i in range(2):
861            r = 1.0
862            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
863            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
864
865            phi = min( min(r*beta_h, 1), phi )
866
867        #Then update using phi limiter
868        #for i in range(3):
869        for i in range(2):
870            hvbar[k,i] = hc[k] + phi*dh[k,i]
871
872    return hvbar
873
874def balance_deep_and_shallow(domain):
875    """Compute linear combination between stage as computed by
876    gradient-limiters limiting using w, and stage computed by
877    gradient-limiters limiting using h (h-limiter).
878    The former takes precedence when heights are large compared to the
879    bed slope while the latter takes precedence when heights are
880    relatively small.  Anything in between is computed as a balanced
881    linear combination in order to avoid numerical disturbances which
882    would otherwise appear as a result of hard switching between
883    modes.
884
885    The h-limiter is always applied irrespective of the order.
886    """
887
888    #Shortcuts
889    wc = domain.quantities['stage'].centroid_values
890    zc = domain.quantities['elevation'].centroid_values
891    hc = wc - zc
892
893    wv = domain.quantities['stage'].vertex_values
894    zv = domain.quantities['elevation'].vertex_values
895    hv = wv-zv
896
897    #Limit h
898    hvbar = h_limiter(domain)
899
900    for k in range(domain.number_of_elements):
901        #Compute maximal variation in bed elevation
902        #  This quantitiy is
903        #    dz = max_i abs(z_i - z_c)
904        #  and it is independent of dimension
905        #  In the 1d case zc = (z0+z1)/2
906        #  In the 2d case zc = (z0+z1+z2)/3
907
908        dz = max(abs(zv[k,0]-zc[k]),
909                 abs(zv[k,1]-zc[k]))#,
910#                 abs(zv[k,2]-zc[k]))
911
912
913        hmin = min( hv[k,:] )
914
915        #Create alpha in [0,1], where alpha==0 means using the h-limited
916        #stage and alpha==1 means using the w-limited stage as
917        #computed by the gradient limiter (both 1st or 2nd order)
918
919        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
920        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
921
922        if dz > 0.0:
923            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
924        else:
925            #Flat bed
926            alpha = 1.0
927
928        alpha  = 0.0
929        #Let
930        #
931        #  wvi be the w-limited stage (wvi = zvi + hvi)
932        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
933        #
934        #
935        #where i=0,1,2 denotes the vertex ids
936        #
937        #Weighted balance between w-limited and h-limited stage is
938        #
939        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
940        #
941        #It follows that the updated wvi is
942        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
943        #
944        # Momentum is balanced between constant and limited
945
946
947        #for i in range(3):
948        #    wv[k,i] = zv[k,i] + hvbar[k,i]
949
950        #return
951
952        if alpha < 1:
953
954            #for i in range(3):
955            for i in range(2):
956                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
957
958            #Momentums at centroids
959            xmomc = domain.quantities['xmomentum'].centroid_values
960#            ymomc = domain.quantities['ymomentum'].centroid_values
961
962            #Momentums at vertices
963            xmomv = domain.quantities['xmomentum'].vertex_values
964#           ymomv = domain.quantities['ymomentum'].vertex_values
965
966            # Update momentum as a linear combination of
967            # xmomc and ymomc (shallow) and momentum
968            # from extrapolator xmomv and ymomv (deep).
969            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
970#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
971
972
973###############################################
974#Boundaries - specific to the shallow water wave equation
975class Reflective_boundary(Boundary):
976    """Reflective boundary returns same conserved quantities as
977    those present in its neighbour volume but reflected.
978
979    This class is specific to the shallow water equation as it
980    works with the momentum quantities assumed to be the second
981    and third conserved quantities.
982    """
983
984    def __init__(self, domain = None):
985        Boundary.__init__(self)
986
987        if domain is None:
988            msg = 'Domain must be specified for reflective boundary'
989            raise msg
990
991        #Handy shorthands
992        #self.stage   = domain.quantities['stage'].edge_values
993        #self.xmom    = domain.quantities['xmomentum'].edge_values
994        #self.ymom    = domain.quantities['ymomentum'].edge_values
995        self.normals = domain.normals
996        self.stage   = domain.quantities['stage'].vertex_values
997        self.xmom    = domain.quantities['xmomentum'].vertex_values
998
999        from Numeric import zeros, Float
1000        #self.conserved_quantities = zeros(3, Float)
1001        self.conserved_quantities = zeros(2, Float)
1002
1003    def __repr__(self):
1004        return 'Reflective_boundary'
1005
1006
1007    def evaluate(self, vol_id, edge_id):
1008        """Reflective boundaries reverses the outward momentum
1009        of the volume they serve.
1010        """
1011
1012        q = self.conserved_quantities
1013        q[0] = self.stage[vol_id, edge_id]
1014        q[1] = self.xmom[vol_id, edge_id]
1015        #q[2] = self.ymom[vol_id, edge_id]
1016        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
1017        #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
1018        normal = self.normals[vol_id,edge_id]
1019
1020        #r = rotate(q, normal, direction = 1)
1021        #r[1] = -r[1]
1022        #q = rotate(r, normal, direction = -1)
1023        r = q
1024        r[1] = normal*r[1]
1025        r[1] = -r[1]
1026        r[1] = normal*r[1]
1027        q = r
1028        #For start interval there is no outward momentum so do not need to
1029        #reverse direction in this case
1030
1031        return q
1032
1033class Dirichlet_boundary(Boundary):
1034    """Dirichlet boundary returns constant values for the
1035    conserved quantities
1036    """
1037
1038
1039    def __init__(self, conserved_quantities=None):
1040        Boundary.__init__(self)
1041
1042        if conserved_quantities is None:
1043            msg = 'Must specify one value for each conserved quantity'
1044            raise msg
1045
1046        from Numeric import array, Float
1047        self.conserved_quantities=array(conserved_quantities).astype(Float)
1048
1049    def __repr__(self):
1050        return 'Dirichlet boundary (%s)' %self.conserved_quantities
1051
1052    def evaluate(self, vol_id=None, edge_id=None):
1053        return self.conserved_quantities
1054
1055
1056#########################
1057#Standard forcing terms:
1058#
1059def gravity(domain):
1060    """Apply gravitational pull in the presence of bed slope
1061    """
1062
1063    from util import gradient
1064    from Numeric import zeros, Float, array, sum
1065
1066    xmom = domain.quantities['xmomentum'].explicit_update
1067    stage = domain.quantities['stage'].explicit_update
1068#    ymom = domain.quantities['ymomentum'].explicit_update
1069
1070    Stage = domain.quantities['stage']
1071    Elevation = domain.quantities['elevation']
1072    #h = Stage.edge_values - Elevation.edge_values
1073    h = Stage.vertex_values - Elevation.vertex_values
1074    b = Elevation.vertex_values
1075    w = Stage.vertex_values
1076
1077    x = domain.get_vertex_coordinates()
1078    g = domain.g
1079
1080    for k in range(domain.number_of_elements):
1081#        avg_h = sum( h[k,:] )/3
1082        avg_h = sum( h[k,:] )/2
1083
1084        #Compute bed slope
1085        #x0, y0, x1, y1, x2, y2 = x[k,:]
1086        x0, x1 = x[k,:]
1087        #z0, z1, z2 = v[k,:]
1088        b0, b1 = b[k,:]
1089
1090        w0, w1 = w[k,:]
1091        wx = gradient(x0, x1, w0, w1)
1092
1093        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
1094        bx = gradient(x0, x1, b0, b1)
1095       
1096        #Update momentum (explicit update is reset to source values)
1097        xmom[k] += -g*bx*avg_h
1098        #xmom[k] = -g*bx*avg_h
1099        #stage[k] = 0.0
1100 
1101 
1102def manning_friction(domain):
1103    """Apply (Manning) friction to water momentum
1104    """
1105
1106    from math import sqrt
1107
1108    w = domain.quantities['stage'].centroid_values
1109    z = domain.quantities['elevation'].centroid_values
1110    h = w-z
1111
1112    uh = domain.quantities['xmomentum'].centroid_values
1113    #vh = domain.quantities['ymomentum'].centroid_values
1114    eta = domain.quantities['friction'].centroid_values
1115
1116    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1117    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1118
1119    N = domain.number_of_elements
1120    eps = domain.minimum_allowed_height
1121    g = domain.g
1122
1123    for k in range(N):
1124        if eta[k] >= eps:
1125            if h[k] >= eps:
1126                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
1127                S = -g * eta[k]**2 * uh[k]
1128                S /= h[k]**(7.0/3)
1129
1130                #Update momentum
1131                xmom_update[k] += S*uh[k]
1132                #ymom_update[k] += S*vh[k]
1133
1134def linear_friction(domain):
1135    """Apply linear friction to water momentum
1136
1137    Assumes quantity: 'linear_friction' to be present
1138    """
1139
1140    from math import sqrt
1141
1142    w = domain.quantities['stage'].centroid_values
1143    z = domain.quantities['elevation'].centroid_values
1144    h = w-z
1145
1146    uh = domain.quantities['xmomentum'].centroid_values
1147#    vh = domain.quantities['ymomentum'].centroid_values
1148    tau = domain.quantities['linear_friction'].centroid_values
1149
1150    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1151#    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1152
1153    N = domain.number_of_elements
1154    eps = domain.minimum_allowed_height
1155    g = domain.g #Not necessary? Why was this added?
1156
1157    for k in range(N):
1158        if tau[k] >= eps:
1159            if h[k] >= eps:
1160                S = -tau[k]/h[k]
1161
1162                #Update momentum
1163                xmom_update[k] += S*uh[k]
1164 #              ymom_update[k] += S*vh[k]
1165
1166
1167
1168def check_forcefield(f):
1169    """Check that f is either
1170    1: a callable object f(t,x,y), where x and y are vectors
1171       and that it returns an array or a list of same length
1172       as x and y
1173    2: a scalar
1174    """
1175
1176    from Numeric import ones, Float, array
1177
1178
1179    if callable(f):
1180        #N = 3
1181        N = 2
1182        #x = ones(3, Float)
1183        #y = ones(3, Float)
1184        x = ones(2, Float)
1185        #y = ones(2, Float)
1186       
1187        try:
1188            #q = f(1.0, x=x, y=y)
1189            q = f(1.0, x=x)
1190        except Exception, e:
1191            msg = 'Function %s could not be executed:\n%s' %(f, e)
1192            #FIXME: Reconsider this semantics
1193            raise msg
1194
1195        try:
1196            q = array(q).astype(Float)
1197        except:
1198            msg = 'Return value from vector function %s could ' %f
1199            msg += 'not be converted into a Numeric array of floats.\n'
1200            msg += 'Specified function should return either list or array.'
1201            raise msg
1202
1203        #Is this really what we want?
1204        msg = 'Return vector from function %s ' %f
1205        msg += 'must have same lenght as input vectors'
1206        assert len(q) == N, msg
1207
1208    else:
1209        try:
1210            f = float(f)
1211        except:
1212            msg = 'Force field %s must be either a scalar' %f
1213            msg += ' or a vector function'
1214            raise msg
1215    return f
1216
1217class Wind_stress:
1218    """Apply wind stress to water momentum in terms of
1219    wind speed [m/s] and wind direction [degrees]
1220    """
1221
1222    def __init__(self, *args, **kwargs):
1223        """Initialise windfield from wind speed s [m/s]
1224        and wind direction phi [degrees]
1225
1226        Inputs v and phi can be either scalars or Python functions, e.g.
1227
1228        W = Wind_stress(10, 178)
1229
1230        #FIXME - 'normal' degrees are assumed for now, i.e. the
1231        vector (1,0) has zero degrees.
1232        We may need to convert from 'compass' degrees later on and also
1233        map from True north to grid north.
1234
1235        Arguments can also be Python functions of t,x,y as in
1236
1237        def speed(t,x,y):
1238            ...
1239            return s
1240
1241        def angle(t,x,y):
1242            ...
1243            return phi
1244
1245        where x and y are vectors.
1246
1247        and then pass the functions in
1248
1249        W = Wind_stress(speed, angle)
1250
1251        The instantiated object W can be appended to the list of
1252        forcing_terms as in
1253
1254        Alternatively, one vector valued function for (speed, angle)
1255        can be applied, providing both quantities simultaneously.
1256        As in
1257        W = Wind_stress(F), where returns (speed, angle) for each t.
1258
1259        domain.forcing_terms.append(W)
1260        """
1261
1262        from config import rho_a, rho_w, eta_w
1263        from Numeric import array, Float
1264
1265        if len(args) == 2:
1266            s = args[0]
1267            phi = args[1]
1268        elif len(args) == 1:
1269            #Assume vector function returning (s, phi)(t,x,y)
1270            vector_function = args[0]
1271            #s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1272            #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1273            s = lambda t,x: vector_function(t,x=x)[0]
1274            phi = lambda t,x: vector_function(t,x=x)[1]
1275        else:
1276           #Assume info is in 2 keyword arguments
1277
1278           if len(kwargs) == 2:
1279               s = kwargs['s']
1280               phi = kwargs['phi']
1281           else:
1282               raise 'Assumes two keyword arguments: s=..., phi=....'
1283
1284        print 'phi', phi
1285        self.speed = check_forcefield(s)
1286        self.phi = check_forcefield(phi)
1287
1288        self.const = eta_w*rho_a/rho_w
1289
1290
1291    def __call__(self, domain):
1292        """Evaluate windfield based on values found in domain
1293        """
1294
1295        from math import pi, cos, sin, sqrt
1296        from Numeric import Float, ones, ArrayType
1297
1298        xmom_update = domain.quantities['xmomentum'].explicit_update
1299        #ymom_update = domain.quantities['ymomentum'].explicit_update
1300
1301        N = domain.number_of_elements
1302        t = domain.time
1303
1304        if callable(self.speed):
1305            xc = domain.get_centroid_coordinates()
1306            #s_vec = self.speed(t, xc[:,0], xc[:,1])
1307            s_vec = self.speed(t, xc)
1308        else:
1309            #Assume s is a scalar
1310
1311            try:
1312                s_vec = self.speed * ones(N, Float)
1313            except:
1314                msg = 'Speed must be either callable or a scalar: %s' %self.s
1315                raise msg
1316
1317
1318        if callable(self.phi):
1319            xc = domain.get_centroid_coordinates()
1320            #phi_vec = self.phi(t, xc[:,0], xc[:,1])
1321            phi_vec = self.phi(t, xc)
1322        else:
1323            #Assume phi is a scalar
1324
1325            try:
1326                phi_vec = self.phi * ones(N, Float)
1327            except:
1328                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1329                raise msg
1330
1331        #assign_windfield_values(xmom_update, ymom_update,
1332        #                        s_vec, phi_vec, self.const)
1333        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)
1334
1335
1336#def assign_windfield_values(xmom_update, ymom_update,
1337#                            s_vec, phi_vec, const):
1338def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
1339    """Python version of assigning wind field to update vectors.
1340    A c version also exists (for speed)
1341    """
1342    from math import pi, cos, sin, sqrt
1343
1344    N = len(s_vec)
1345    for k in range(N):
1346        s = s_vec[k]
1347        phi = phi_vec[k]
1348
1349        #Convert to radians
1350        phi = phi*pi/180
1351
1352        #Compute velocity vector (u, v)
1353        u = s*cos(phi)
1354        v = s*sin(phi)
1355
1356        #Compute wind stress
1357        #S = const * sqrt(u**2 + v**2)
1358        S = const * u
1359        xmom_update[k] += S*u
1360        #ymom_update[k] += S*v
Note: See TracBrowser for help on using the repository browser.