source: anuga_work/development/anuga_1d/shallow_water_domain_suggestion1.py @ 5827

Last change on this file since 5827 was 5827, checked in by sudi, 16 years ago

Adding new versions of shallow_water_domain

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