source: trunk/anuga_work/development/sudi/sw_1d/numerical_entropy_production/sww_domain_shv.py @ 8041

Last change on this file since 8041 was 8041, checked in by mungkasi, 14 years ago

Adding numerical_entropy_production folder containing codes on some smoothness indicators of solutions.
File size: 27.4 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
40John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
41Geoscience Australia, 2006
42
43Sudi Mungkasi, ANU, 2010
44"""
45
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        conserved_quantities = ['stage', 'xmomentum'] #['height', 'xmomentum']
55        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
56        other_quantities = ['friction']
57        Generic_Domain.__init__(self, coordinates, boundary,
58                                conserved_quantities, evolved_quantities, other_quantities,
59                                tagged_elements)
60       
61        from config import minimum_allowed_height, g, h0
62        self.minimum_allowed_height = minimum_allowed_height
63        self.g = g
64        self.h0 = h0
65
66        self.forcing_terms.append(gravity_F2)
67        #self.forcing_terms.append(manning_friction)
68
69        #Realtime visualisation
70        self.visualiser = None
71        self.visualise  = False
72        self.visualise_color_stage = False
73        self.visualise_stage_range = 1.0
74        self.visualise_timer = True
75        self.visualise_range_z = None
76       
77        #Stored output
78        self.store = True
79        self.format = 'sww'
80        self.smooth = True
81
82        #Evolve parametrs
83        self.cfl = 1.0
84       
85        #Reduction operation for get_vertex_values
86        from util import mean
87        self.reduction = mean
88        #self.reduction = min  #Looks better near steep slopes
89
90        self.quantities_to_be_stored = ['stage','xmomentum']
91
92        self.__doc__ = 'sww_domain_shv'
93        self.check_integrity()
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        if q is None:
111            self.quantities_to_be_stored = []           
112            self.store = False
113            return
114
115        if isinstance(q, basestring):
116            q = [q] # Turn argument into a list
117
118        #Check correcness   
119        for quantity_name in q:
120            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
121            assert quantity_name in self.conserved_quantities, msg
122       
123        self.quantities_to_be_stored = q
124       
125
126    def initialise_visualiser(self,scale_z=1.0,rect=None):
127        #Realtime visualisation
128        if self.visualiser is None:
129            from realtime_visualisation_new import Visualiser
130            self.visualiser = Visualiser(self,scale_z,rect)
131            self.visualiser.setup['elevation']=True
132            self.visualiser.updating['stage']=True
133        self.visualise = True
134        if self.visualise_color_stage == True:
135            self.visualiser.coloring['stage'] = True
136            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
137        print 'initialise visualiser'
138        print self.visualiser.setup
139        print self.visualiser.updating
140
141    def check_integrity(self):
142        Generic_Domain.check_integrity(self)
143        #Check that we are solving the shallow water wave equation
144
145        msg = 'First conserved quantity must be "stage"'
146        assert self.conserved_quantities[0] == 'stage', msg
147        msg = 'Second conserved quantity must be "xmomentum"'
148        assert self.conserved_quantities[1] == 'xmomentum', msg
149               
150        msg = 'First evolved quantity must be "stage"'
151        assert self.evolved_quantities[0] == 'stage', msg
152        msg = 'Second evolved quantity must be "xmomentum"'
153        assert self.evolved_quantities[1] == 'xmomentum', msg
154        msg = 'Third evolved quantity must be "elevation"'
155        assert self.evolved_quantities[2] == 'elevation', msg
156        msg = 'Fourth evolved quantity must be "height"'
157        assert self.evolved_quantities[3] == 'height', msg
158        msg = 'Fifth evolved quantity must be "velocity"'
159        assert self.evolved_quantities[4] == 'velocity', msg
160
161    def extrapolate_second_order_sw(self):
162        #Call correct module function
163        #(either from this module or C-extension)
164        extrapolate_second_order_sw(self)
165
166    def compute_fluxes(self):
167        #Call correct module function(either from this module or C-extension)
168        #compute_fluxes_C_wellbalanced(self)
169        #compute_fluxes_quantity_and_entropy(self)
170        compute_fluxes_C_include_entropy(self)
171
172    def compute_timestep(self):
173        #Call correct module function
174        compute_timestep(self)
175
176    def distribute_to_vertices_and_edges(self):
177        #Call correct module function(either from this module or C-extension)
178        distribute_to_vertices_and_edges_shv(self)
179        #distribute_to_vertices_and_edges_shm(self)
180       
181    def evolve(self, yieldstep = None, finaltime = None, duration = None, skip_initial_step = False):
182        #Call basic machinery from parent class
183        for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration, skip_initial_step):
184            yield(t)
185
186    def initialise_storage(self):
187        """Create and initialise self.writer object for storing data.
188        Also, save x and bed elevation
189        """
190
191        import data_manager
192
193        #Initialise writer
194        self.writer = data_manager.get_dataobject(self, mode = 'w')
195
196        #Store vertices and connectivity
197        self.writer.store_connectivity()
198
199
200    def store_timestep(self, name):
201        """Store named quantity and time.
202
203        Precondition:
204           self.write has been initialised
205        """
206        self.writer.store_timestep(name)
207
208
209#=============== End of Shallow Water Domain ===============================
210
211# Compute flux definition
212def flux_function_quantity_and_entropy(normal, ql, qr, zl, zr):
213    """Compute fluxes between volumes for the shallow water wave equation
214    cast in terms of w = h+z using the 'central scheme' as described in
215
216    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
217    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
218    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
219
220    The implemented formula is given in equation (3.15) on page 714
221
222    Conserved quantities w, uh, are stored as elements 0 and 1
223    in the numerical vectors ql an qr.
224
225    Bed elevations zl and zr.
226    """
227
228    from config import g, epsilon, h0
229    from math import sqrt
230    from numpy import array
231
232    q_left = ql
233    q_left[1] = q_left[1]*normal
234    q_right = qr
235    q_right[1] = q_right[1]*normal
236
237    z = 0.5*(zl+zr) #Take average of field values
238
239    w_left  = q_left[0]   #w=h+z
240    h_left  = w_left-z
241    uh_left = q_left[1]
242    z_left = zl
243    if h_left < epsilon:
244        u_left = 0.0  #Could have been negative
245        h_left = 0.0
246    else:
247        u_left  = uh_left/(h_left +  h0/h_left)
248    uh_left = u_left*h_left
249
250
251    w_right  = q_right[0]  #w=h+z
252    h_right  = w_right-z
253    uh_right = q_right[1]
254    z_right = zr
255    if h_right < epsilon:
256        u_right = 0.0  #Could have been negative
257        h_right = 0.0
258    else:
259        u_right  = uh_right/(h_right + h0/h_right)
260    uh_right = u_right*h_right
261
262    soundspeed_left  = sqrt(g*h_left)
263    soundspeed_right = sqrt(g*h_right)
264
265    #Maximal wave speed
266    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
267   
268    #Minimal wave speed
269    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
270   
271    #Flux computation
272    flux_left  = array([u_left*h_left,
273                        u_left*uh_left + 0.5*g*h_left*h_left])
274    flux_right = array([u_right*h_right,
275                        u_right*uh_right + 0.5*g*h_right*h_right])
276
277    entropy_left = 0.5*h_left*u_left**2.0 + 0.5*g*h_left**2.0 + g*h_left*z_left
278    entropy_right= 0.5*h_right*u_right**2.0 + 0.5*g*h_right**2.0 + g*h_right*z_right
279   
280    entropy_flux_left  = 0.5*h_left*u_left**3.0 + g*h_left**2.0*u_left + g*h_left*z_left*u_left
281    entropy_flux_right = 0.5*h_right*u_right**3.0 + g*h_right**2.0*u_right + g*h_right*z_right*u_right
282   
283
284    denom = s_max-s_min
285    if denom == 0.0:
286        edgeflux = array([0.0, 0.0])
287        max_speed = 0.0
288        entropy_edgeflux = 0.0
289    else:
290        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
291        edgeflux += s_max*s_min*(q_right-q_left)/denom
292       
293        edgeflux[1] = edgeflux[1]*normal
294
295        max_speed = max(abs(s_max), abs(s_min))
296
297        entropy_edgeflux = (s_max*entropy_flux_left - s_min*entropy_flux_right)/denom
298        entropy_edgeflux += s_max*s_min*(entropy_right-entropy_left)/denom
299
300    return edgeflux, max_speed, entropy_edgeflux
301
302def compute_fluxes_quantity_and_entropy(domain):
303    """Compute all fluxes and the timestep suitable for all volumes
304    in domain.
305
306    Compute total flux for each conserved quantity using "flux_function"
307
308    Fluxes across each edge are scaled by edgelengths and summed up
309    Resulting flux is then scaled by area and stored in
310    explicit_update for each of the three conserved quantities
311    stage, xmomentum and ymomentum
312
313    The maximal allowable speed computed by the flux_function for each volume
314    is converted to a timestep that must not be exceeded. The minimum of
315    those is computed as the next overall timestep.
316
317    Post conditions:
318      domain.explicit_update is reset to computed flux values
319      domain.timestep is set to the largest step satisfying all volumes.
320    """
321
322    import sys
323    from numpy import zeros
324    from config import g
325
326    N = domain.number_of_elements
327
328    #Shortcuts
329    Stage = domain.quantities['stage']
330    Xmom = domain.quantities['xmomentum']
331    Bed = domain.quantities['elevation']
332
333    #Arrays
334    stage = Stage.vertex_values
335    xmom =  Xmom.vertex_values
336    bed =   Bed.vertex_values
337   
338    stage_bdry = Stage.boundary_values
339    xmom_bdry =  Xmom.boundary_values
340   
341    flux = zeros(2, float) #Work array for summing up fluxes
342    ql = zeros(2, float)
343    qr = zeros(2, float)
344
345    #Loop
346    timestep = float(sys.maxint)
347    enter = True
348    for k in range(N):
349        flux[:] = 0.0  #Reset work array
350        entropy_flux = 0.0 #Reset too
351        for i in range(2):
352            #Quantities inside volume facing neighbour i
353            ql = [stage[k, i], xmom[k, i]]
354            zl = bed[k, i]
355
356            #Quantities at neighbour on nearest face
357            n = domain.neighbours[k,i]
358            if n < 0:
359                m = -n-1 #Convert negative flag to index
360                qr[0] = stage_bdry[m]
361                qr[1] = xmom_bdry[m]
362                zr = zl #Extend bed elevation to boundary
363            else:
364                m = domain.neighbour_vertices[k,i]
365                qr[0] = stage[n, m]
366                qr[1] = xmom[n, m]
367                zr = bed[n, m]
368
369            #Outward pointing normal vector
370            normal = domain.normals[k, i]
371       
372            #Flux computation using provided function
373            edgeflux, max_speed, entropy_edgeflux = flux_function_quantity_and_entropy(normal, ql, qr, zl, zr)
374
375            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
376            flux -= edgeflux
377            entropy_flux -= entropy_edgeflux           
378            #Update optimal_timestep
379            try:
380                timestep = min(timestep, domain.CFL*0.5*domain.areas[k]/max_speed)
381            except ZeroDivisionError:
382                pass
383
384        #Normalise by area and store for when all conserved quantities get updated
385        flux /= domain.areas[k]
386        entropy_flux /= domain.areas[k]       
387
388        Stage.explicit_update[k] = flux[0]
389        Xmom.explicit_update[k] =  flux[1]
390        domain.entropy_explicit_update[k] = entropy_flux       
391
392    domain.flux_timestep = timestep
393
394
395# ###################################
396def compute_fluxes_C_include_entropy(domain):
397    import sys
398    from numpy import zeros   
399   
400    N = domain.number_of_elements
401    timestep = float(sys.maxint)
402    epsilon = domain.epsilon
403    g = domain.g
404
405    cfl = domain.CFL
406    neighbours = domain.neighbours
407    neighbour_vertices = domain.neighbour_vertices
408    normals = domain.normals
409    areas = domain.areas
410   
411    Stage    = domain.quantities['stage']
412    Xmom     = domain.quantities['xmomentum']
413    Bed      = domain.quantities['elevation']
414    Height   = domain.quantities['height']
415    Velocity = domain.quantities['velocity']
416
417
418    stage_boundary_values = Stage.boundary_values
419    xmom_boundary_values  = Xmom.boundary_values
420    bed_boundary_values   = Bed.boundary_values
421    height_boundary_values= Height.boundary_values
422    vel_boundary_values   = Velocity.boundary_values
423   
424    stage_explicit_update = Stage.explicit_update
425    xmom_explicit_update  = Xmom.explicit_update
426    bed_explicit_values   = Bed.explicit_update
427    height_explicit_values= Height.explicit_update
428    vel_explicit_values   = Velocity.explicit_update
429
430    entropy_explicit_update = domain.entropy_explicit_update
431   
432    max_speed_array = domain.max_speed_array
433   
434    domain.distribute_to_vertices_and_edges()
435    domain.update_boundary()
436   
437    stage_V = Stage.vertex_values
438    xmom_V  = Xmom.vertex_values
439    bed_V   = Bed.vertex_values
440    height_V= Height.vertex_values
441    vel_V   = Velocity.vertex_values
442
443    number_of_elements = len(stage_V)
444
445    from comp_flux_ext_include_entropy import compute_fluxes_ext_wellbalanced
446    domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
447                                  epsilon,
448                                  g,
449                                  cfl,
450                                                           
451                                  neighbours,
452                                  neighbour_vertices,
453                                  normals,
454                                  areas,
455                                                           
456                                  stage_V,
457                                  xmom_V,
458                                  bed_V,
459                                  height_V,
460                                  vel_V,
461                                                           
462                                  stage_boundary_values,
463                                  xmom_boundary_values,
464                                  bed_boundary_values,
465                                  height_boundary_values,
466                                  vel_boundary_values,
467                                                           
468                                  stage_explicit_update,
469                                  xmom_explicit_update,
470                                  bed_explicit_values,
471                                  height_explicit_values,
472                                  vel_explicit_values,
473
474                                  entropy_explicit_update,
475                                                           
476                                  number_of_elements,
477                                  max_speed_array)
478
479
480
481# ###################################
482def compute_fluxes_C_wellbalanced(domain):
483    import sys
484    from numpy import zeros   
485   
486    N = domain.number_of_elements
487    timestep = float(sys.maxint)
488    epsilon = domain.epsilon
489    g = domain.g
490    neighbours = domain.neighbours
491    neighbour_vertices = domain.neighbour_vertices
492    normals = domain.normals
493    areas = domain.areas
494   
495    Stage    = domain.quantities['stage']
496    Xmom     = domain.quantities['xmomentum']
497    Bed      = domain.quantities['elevation']
498    Height   = domain.quantities['height']
499    Velocity = domain.quantities['velocity']
500
501
502    stage_boundary_values = Stage.boundary_values
503    xmom_boundary_values  = Xmom.boundary_values
504    bed_boundary_values   = Bed.boundary_values
505    height_boundary_values= Height.boundary_values
506    vel_boundary_values   = Velocity.boundary_values
507   
508    stage_explicit_update = Stage.explicit_update
509    xmom_explicit_update  = Xmom.explicit_update
510    bed_explicit_values   = Bed.explicit_update
511    height_explicit_values= Height.explicit_update
512    vel_explicit_values   = Velocity.explicit_update
513   
514    max_speed_array = domain.max_speed_array
515   
516    domain.distribute_to_vertices_and_edges()
517    domain.update_boundary()
518   
519    stage_V = Stage.vertex_values
520    xmom_V  = Xmom.vertex_values
521    bed_V   = Bed.vertex_values
522    height_V= Height.vertex_values
523    vel_V   = Velocity.vertex_values
524
525    number_of_elements = len(stage_V)
526
527    from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced
528    domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
529                                  epsilon,
530                                  g,
531                                                           
532                                  neighbours,
533                                  neighbour_vertices,
534                                  normals,
535                                  areas,
536                                                           
537                                  stage_V,
538                                  xmom_V,
539                                  bed_V,
540                                  height_V,
541                                  vel_V,
542                                                           
543                                  stage_boundary_values,
544                                  xmom_boundary_values,
545                                  bed_boundary_values,
546                                  height_boundary_values,
547                                  vel_boundary_values,
548                                                           
549                                  stage_explicit_update,
550                                  xmom_explicit_update,
551                                  bed_explicit_values,
552                                  height_explicit_values,
553                                  vel_explicit_values,
554                                                           
555                                  number_of_elements,
556                                  max_speed_array)
557
558
559# ###################################
560# Module functions for gradient limiting (distribute_to_vertices_and_edges)
561def distribute_to_vertices_and_edges_shv(domain):
562    """Distribution from centroids to vertices specific to the
563    shallow water wave
564    equation.
565
566    It will ensure that h (w-z) is always non-negative even in the
567    presence of steep bed-slopes by taking a weighted average between shallow
568    and deep cases.
569
570    In addition, all conserved quantities get distributed as per either a
571    constant (order==1) or a piecewise linear function (order==2).
572
573    FIXME: more explanation about removal of artificial variability etc
574
575    Precondition:
576      All quantities defined at centroids and bed elevation defined at
577      vertices.
578
579    Postcondition
580      Conserved quantities defined at vertices
581
582    """
583
584    #from config import optimised_gradient_limiter
585
586    #Remove very thin layers of water
587    #protect_against_infinitesimal_and_negative_heights(domain) 
588
589    import sys
590    from numpy import zeros   
591    from config import epsilon, h0
592
593    N = domain.number_of_elements
594
595    #Shortcuts
596    Stage = domain.quantities['stage']
597    Xmom = domain.quantities['xmomentum']
598    Bed = domain.quantities['elevation']
599    Height = domain.quantities['height']
600    Velocity = domain.quantities['velocity']
601
602    #Arrays   
603    w_C   = Stage.centroid_values   
604    uh_C  = Xmom.centroid_values   
605    z_C   = Bed.centroid_values
606    h_C   = Height.centroid_values
607    u_C   = Velocity.centroid_values
608
609    for i in range(N):
610        h_C[i] = w_C[i] - z_C[i]                                               
611        if h_C[i] <= epsilon:
612            uh_C[i] = 0.0
613            u_C[i] = 0.0
614            #w_C[i] = z_C[i]
615        else:
616            u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
617
618    for name in ['stage', 'height', 'velocity']:
619        Q = domain.quantities[name]
620        if domain.order == 1:
621            Q.extrapolate_first_order()
622        elif domain.order == 2:
623            Q.extrapolate_second_order()
624        else:
625            raise 'Unknown order'
626
627    stage_V = domain.quantities['stage'].vertex_values
628    bed_V   = domain.quantities['elevation'].vertex_values
629    h_V     = domain.quantities['height'].vertex_values
630    u_V     = domain.quantities['velocity'].vertex_values               
631    xmom_V  = domain.quantities['xmomentum'].vertex_values
632
633    bed_V[:] = stage_V - h_V
634    xmom_V[:] = u_V * h_V
635   
636    return
637
638def distribute_to_vertices_and_edges_shm(domain):
639    # shm stands for STAGE, HEIGHT, MOMENTUM
640    """Distribution from centroids to vertices specific to the
641    shallow water wave
642    equation.
643
644    It will ensure that h (w-z) is always non-negative even in the
645    presence of steep bed-slopes by taking a weighted average between shallow
646    and deep cases.
647
648    In addition, all conserved quantities get distributed as per either a
649    constant (order==1) or a piecewise linear function (order==2).
650
651    FIXME: more explanation about removal of artificial variability etc
652
653    Precondition:
654      All quantities defined at centroids and bed elevation defined at
655      vertices.
656
657    Postcondition
658      Conserved quantities defined at vertices
659
660    """
661
662    #from config import optimised_gradient_limiter
663
664    #Remove very thin layers of water
665    #protect_against_infinitesimal_and_negative_heights(domain) 
666
667    import sys
668    from numpy import array, zeros
669    from config import epsilon, h0
670
671    N = domain.number_of_elements
672
673    #Shortcuts
674    Stage = domain.quantities['stage']
675    Xmom = domain.quantities['xmomentum']
676    Bed = domain.quantities['elevation']
677    Height = domain.quantities['height']
678    Velocity = domain.quantities['velocity']
679
680    #Arrays   
681    w_C   = Stage.centroid_values   
682    uh_C  = Xmom.centroid_values   
683    z_C   = Bed.centroid_values
684    h_C   = Height.centroid_values
685    u_C   = Velocity.centroid_values
686
687   
688    for i in range(N):
689        h_C[i] = w_C[i] - z_C[i]                                               
690        if h_C[i] <= epsilon:
691            uh_C[i] = 0.0
692            u_C[i] = 0.0
693        else:
694            u_C[i] = uh_C[i]/(h_C[i] + h0/h_C[i])
695
696    for name in ['stage', 'height', 'xmomentum']:
697        Q = domain.quantities[name]
698        if domain.order == 1:
699            Q.extrapolate_first_order()
700        elif domain.order == 2:
701            Q.extrapolate_second_order()
702        else:
703            raise 'Unknown order'
704
705    stage_V = domain.quantities['stage'].vertex_values
706    bed_V   = domain.quantities['elevation'].vertex_values
707    h_V     = domain.quantities['height'].vertex_values
708    u_V     = domain.quantities['velocity'].vertex_values               
709    xmom_V  = domain.quantities['xmomentum'].vertex_values
710
711    bed_V[:] = stage_V - h_V
712
713    for i in range(N):                                         
714        if min(h_V[i]) <= 0.0:
715            h_V[i] = array([0.0, 0.0])
716            stage_V[i] = bed_V[i]
717            xmom_V[i] = array([0.0, 0.0])
718   
719    u_V[:]  = xmom_V/(h_V + h0/h_V)   
720   
721    return
722
723
724   
725#
726def protect_against_infinitesimal_and_negative_heights(domain):
727    """Protect against infinitesimal heights and associated high velocities
728    """
729
730    #Shortcuts
731    wc = domain.quantities['stage'].centroid_values
732    zc = domain.quantities['elevation'].centroid_values
733    xmomc = domain.quantities['xmomentum'].centroid_values
734    hc = wc - zc  #Water depths at centroids
735
736    zv = domain.quantities['elevation'].vertex_values
737    wv = domain.quantities['stage'].vertex_values
738    hv = wv-zv
739    xmomv = domain.quantities['xmomentum'].vertex_values
740    #remove the above two lines and corresponding code below
741
742    #Update
743    for k in range(domain.number_of_elements):
744
745        if hc[k] < domain.minimum_allowed_height:
746            #Control stage
747            if hc[k] < domain.epsilon:
748                wc[k] = zc[k] # Contain 'lost mass' error
749                wv[k,0] = zv[k,0]
750                wv[k,1] = zv[k,1]
751               
752            xmomc[k] = 0.0
753           
754   
755
756
757
758#########################
759#Standard forcing terms:
760#
761def gravity(domain):
762    """Apply gravitational pull in the presence of bed slope
763    """
764
765    from util import gradient
766    from numpy import zeros, array, sum   
767
768    xmom = domain.quantities['xmomentum'].explicit_update
769    stage = domain.quantities['stage'].explicit_update
770
771    Stage = domain.quantities['stage']
772    Elevation = domain.quantities['elevation']
773    h = Stage.vertex_values - Elevation.vertex_values
774    b = Elevation.vertex_values
775    w = Stage.vertex_values
776
777    x = domain.get_vertex_coordinates()
778    g = domain.g
779
780    for k in range(domain.number_of_elements):
781        avg_h = sum( h[k,:] )/2
782
783        #Compute bed slope
784        x0, x1 = x[k,:]
785        b0, b1 = b[k,:]
786
787        w0, w1 = w[k,:]
788        wx = gradient(x0, x1, w0, w1)
789        bx = gradient(x0, x1, b0, b1)
790       
791        #Update momentum (explicit update is reset to source values)
792        xmom[k] += -g*bx*avg_h
793
794def gravity_F2(domain):
795    """Apply gravitational pull in the presence of bed slope
796    """
797
798    from util import gradient
799    from numpy import zeros, array, sum   
800    from parameters import F2
801
802    xmom = domain.quantities['xmomentum'].explicit_update
803    stage = domain.quantities['stage'].explicit_update
804
805    Stage = domain.quantities['stage']
806    Elevation = domain.quantities['elevation']
807    h = Stage.vertex_values - Elevation.vertex_values
808    b = Elevation.vertex_values
809    w = Stage.vertex_values
810
811    x = domain.get_vertex_coordinates()
812    g = domain.g
813
814    for k in range(domain.number_of_elements):
815        avg_h = sum( h[k,:] )/2
816
817        #Compute bed slope
818        x0, x1 = x[k,:]
819        b0, b1 = b[k,:]
820
821        w0, w1 = w[k,:]
822        wx = gradient(x0, x1, w0, w1)
823        bx = gradient(x0, x1, b0, b1)
824       
825        #Update momentum (explicit update is reset to source values)
826        xmom[k] += -g*bx*avg_h + avg_h*F2
827
828 
829def manning_friction(domain):
830    """Apply (Manning) friction to water momentum
831    """
832
833    from math import sqrt
834
835    w = domain.quantities['stage'].centroid_values
836    z = domain.quantities['elevation'].centroid_values
837    h = w-z
838
839    uh = domain.quantities['xmomentum'].centroid_values
840    eta = domain.quantities['friction'].centroid_values
841
842    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
843
844    N = domain.number_of_elements
845    eps = domain.minimum_allowed_height
846    g = domain.g
847
848    for k in range(N):
849        if eta[k] >= eps:
850            if h[k] >= eps:
851                S = -g * eta[k]**2 * uh[k]
852                S /= h[k]**(7.0/3)
853
854                #Update momentum
855                xmom_update[k] += S*uh[k]
Note: See TracBrowser for help on using the repository browser.