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

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