source: anuga_work/development/pipeflow/shallow_water_domain.py @ 7777

Last change on this file since 7777 was 7777, checked in by steve, 14 years ago

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