Ignore:
Timestamp:
Jun 15, 2010, 5:34:24 PM (12 years ago)
Author:
steve
Message:

Adding in a few new files

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py

    r7839 r7842  
    33the shallow water wave equation.
    44
    5 This module contains a specialisation of class Domain from module domain.py
    6 consisting of methods specific to the Shallow Water Wave Equation
     5This module contains a specialisation of class Generic_domain from module
     6generic_domain.py
     7consisting of methods specific to Channel flow using the Shallow Water Wave Equation
    78
    89This particular modification of the Domain class implements the ability to
     
    1415
    1516where
    16 ------------!!!! NOTE THIS NEEDS UPDATING !!!!------------------
    17 U = [A, Q]
    18 E = [Q, Q^2/A + gh^2/2]
     17
     18U = [A, Q] = [b*h, u*b*h]
     19E = [Q, Q^2/A + g*b*h^2/2]
    1920S represents source terms forcing the system
    20 (e.g. gravity, friction, wind stress, ...)
     21(e.g. gravity, boundary_stree, friction, wind stress, ...)
     22gravity = -g*b*h*z_x
     23boundary_stress = 1/2*g*b_x*h^2
    2124
    2225and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
     
    2528
    2629symbol    variable name    explanation
     30A         area             Wetted area = b*h
     31Q         discharge        flux of water = u*b*h
    2732x         x                horizontal distance from origin [m]
    2833z         elevation        elevation of bed on which flow is modelled [m]
     
    3136u                          speed in the x direction [m/s]
    3237uh        xmomentum        momentum in the x direction [m^2/s]
    33 
     38b         width            width of channel
    3439eta                        mannings friction coefficient [to appear]
    3540nu                         wind stress coefficient [to appear]
    3641
    37 The conserved quantities are w, uh
     42The conserved quantities are A, Q
    3843--------------------------------------------------------------------------
    3944For details see e.g.
     
    4853
    4954
    50 from domain import *
    51 Generic_Domain = Domain #Rename
     55from anuga_1d.generic.generic_domain import *
     56import numpy
     57
    5258
    5359#Shallow water domain
    54 class Domain(Generic_Domain):
     60class Domain(Generic_domain):
    5561
    5662    def __init__(self, coordinates, boundary = None, tagged_elements = None):
     
    5965        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','stage']
    6066        other_quantities = ['friction']
    61         Generic_Domain.__init__(self,
     67        Generic_domain.__init__(self,
    6268                                coordinates = coordinates,
    6369                                boundary = boundary,
     
    6773                                tagged_elements = tagged_elements)
    6874       
    69         from config import minimum_allowed_height, g, h0
     75        from anuga_1d.config import minimum_allowed_height, g, h0
    7076        self.minimum_allowed_height = minimum_allowed_height
    7177        self.g = g
     
    7480        self.discontinousb = False
    7581
    76      
     82
     83        # forcing terms gravity and boundary stress are included in the flux calculation
    7784        #self.forcing_terms.append(gravity)
    7885        #self.forcing_terms.append(boundary_stress)
     
    8895       
    8996        #Reduction operation for get_vertex_values
    90         from util import mean
     97        from anuga_1d.generic.util import mean
    9198        self.reduction = mean
    9299        #self.reduction = min  #Looks better near steep slopes
     
    94101        self.set_quantities_to_be_stored(['area','discharge'])
    95102
    96         self.__doc__ = 'channel_domain_Ab'
     103        self.__doc__ = 'channel_domain'
    97104
    98105        self.check_integrity()
     
    118125        msg = 'Fifth evolved quantity must be "velocity"'
    119126        assert self.evolved_quantities[4] == 'velocity', msg
    120         msg = 'Fifth evolved quantity must be "width"'
     127        msg = 'Sixth evolved quantity must be "width"'
    121128        assert self.evolved_quantities[5] == 'width', msg
    122 
    123         Generic_Domain.check_integrity(self)
     129        msg = 'Seventh evolved quantity must be "stage"'
     130        assert self.evolved_quantities[6] == 'stage', msg
     131
     132        Generic_domain.check_integrity(self)
    124133
    125134    def compute_fluxes(self):
     
    140149#-----------------------------------
    141150def compute_fluxes_channel(domain):
    142     from Numeric import zeros, Float
    143151    import sys
    144 
    145        
    146152    timestep = float(sys.maxint)
    147153
    148     area    = domain.quantities['area']
    149     discharge     = domain.quantities['discharge']
    150     bed      = domain.quantities['elevation']
    151     height   = domain.quantities['height']
    152     velocity = domain.quantities['velocity']
    153     width    = domain.quantities['width']
    154 
    155 
    156     from channel_domain_ext import compute_fluxes_channel_ext
     154    area       = domain.quantities['area']
     155    discharge  = domain.quantities['discharge']
     156    bed        = domain.quantities['elevation']
     157    height     = domain.quantities['height']
     158    velocity   = domain.quantities['velocity']
     159    width      = domain.quantities['width']
     160
     161
     162    from anuga_1d.channel.channel_domain_ext import compute_fluxes_channel_ext
    157163    domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,velocity,width)
    158164
     
    168174
    169175    import sys
    170     from Numeric import zeros, Float
    171     from config import epsilon, h0
    172    ##  linearb(domain)
    173 
    174 
     176    from anuga_1d.config import epsilon, h0
    175177
    176178   
     
    178180
    179181    #Shortcuts
    180     Area = domain.quantities['area']
    181     Discharge = domain.quantities['discharge']
    182     Bed = domain.quantities['elevation']
    183     Height = domain.quantities['height']
    184     Velocity = domain.quantities['velocity']
    185     Width = domain.quantities['width']
    186     Stage = domain.quantities['stage']
     182    area      = domain.quantities['area']
     183    discharge = domain.quantities['discharge']
     184    bed      = domain.quantities['elevation']
     185    height    = domain.quantities['height']
     186    velocity = domain.quantities['velocity']
     187    width    = domain.quantities['width']
     188    stage    = domain.quantities['stage']
    187189
    188190    #Arrays   
    189     a_C   = Area.centroid_values   
    190     d_C  = Discharge.centroid_values   
    191     z_C   = Bed.centroid_values
    192     h_C   = Height.centroid_values
    193     u_C   = Velocity.centroid_values
    194     b_C   = Width.centroid_values
    195     w_C   = Stage.centroid_values
     191    a_C   = area.centroid_values
     192    d_C   = discharge.centroid_values
     193    z_C   = bed.centroid_values
     194    h_C   = height.centroid_values
     195    u_C   = velocity.centroid_values
     196    b_C   = width.centroid_values
     197    w_C   = stage.centroid_values
    196198
    197199    if domain.setstageflag:
    198         for i in range(len(a_C)):
    199             a_C[i]=(w_C[i]-z_C[i])*b_C[i]
    200      
     200        a_C[:,] = (w_C-z_C)*b_C
    201201        domain.setstageflag = False
    202202
     
    205205       
    206206
    207     h0 = 1.0e-12       
    208     #print id(h_C)
    209     for i in range(N):
    210                                                        
    211         if a_C[i] <= h0:
    212             a_C[i] = 0.0
    213             h_C[i] = 0.0
    214             d_C[i] = 0.0
    215             u_C[i] = 0.0
    216             w_C[i] = z_C[i]
    217            
    218            
    219 
    220         else:
    221 
    222             if b_C[i]<=h0:
    223                 a_C[i] = 0.0
    224                 h_C[i] = 0.0
    225                 d_C[i] = 0.0
    226                 u_C[i] = 0.0
    227                 w_C[i] = z_C[i]
    228 
    229             else:
    230                 h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])
    231                 w_C[i] = h_C[i]+z_C[i]
    232                 u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i])
     207    # FIXME (SGR): Replace with numpy.where to speed up code
     208    h0 = 1.0e-12
     209
     210    h_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C/(b_C+h0/b_C), 0.0 )
     211    u_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C/(a_C+h0/a_C), 0.0 )
     212
     213    a_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C, 0.0 )
     214    b_C[:] = numpy.where( (a_C>h0) | (b_C>h0), b_C, 0.0 )
     215    d_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C, 0.0 )
     216
     217    w_C[:] = h_C + z_C
     218   
     219#    for i in range(N):
     220#
     221#        if a_C[i] <= h0:
     222#           a_C[i] = 0.0
     223#           h_C[i] = 0.0
     224#            d_C[i] = 0.0
     225#            u_C[i] = 0.0
     226#            w_C[i] = z_C[i]
     227#        else:
     228#
     229#            if b_C[i]<=h0:
     230#                a_C[i] = 0.0
     231#                h_C[i] = 0.0
     232#                d_C[i] = 0.0
     233#                u_C[i] = 0.0
     234#                w_C[i] = z_C[i]
     235#
     236#            else:
     237#                h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])
     238#                w_C[i] = h_C[i]+z_C[i]
     239#                u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i])
    233240               
    234    
    235    
    236        
     241
     242    bed.extrapolate_second_order()
     243   
    237244    for name in ['velocity','stage']:
    238245        Q = domain.quantities[name]
     
    243250        else:
    244251            raise 'Unknown order'
    245     a_V  = domain.quantities['area'].vertex_values
    246     w_V  = domain.quantities['stage'].vertex_values                     
    247     z_V  = domain.quantities['elevation'].vertex_values
    248     h_V  = domain.quantities['height'].vertex_values
    249     u_V  = domain.quantities['velocity'].vertex_values         
    250     d_V = domain.quantities['discharge'].vertex_values
    251     b_V = domain.quantities['width'].vertex_values
     252
     253
     254    a_V  = area.vertex_values
     255    w_V  = stage.vertex_values
     256    z_V  = bed.vertex_values
     257    h_V  = height.vertex_values
     258    u_V  = velocity.vertex_values
     259    d_V  = discharge.vertex_values
     260    b_V  = width.vertex_values
    252261
    253262   
    254    
    255     for i in range(len(h_C)):
    256         for j in range(2):
    257 ##             if b_V[i,j] < h0 :
    258 ##                 a_V[i,j]=0
    259 ##                 h_V[i,j]=0
    260 ##                 d_V[i,j]=0
    261 ##                 u_V[i,j]=0
    262 ##             else:
    263                
    264             h_V[i,j] = w_V[i,j]-z_V[i,j]
    265             if h_V[i,j]<h0:
    266                 h_V[i,j]=0
    267                 w_V[i,j]=z_V[i,j]
    268             a_V[i,j] = b_V[i,j]*h_V[i,j]
    269             d_V[i,j]=u_V[i,j]*a_V[i,j]
    270            
    271                
    272                
    273    
    274 
     263    #FIXME (SGR): replace with numpy.where
     264
     265
     266    h_V[:] = w_V-z_V
     267
     268    w_V[:] = numpy.where(h_V > h0, w_V, z_V)
     269    h_V[:] = numpy.where(h_V > h0, h_V, 0.0)
     270    a_V[:] = b_V*h_V
     271    d_V[:] = u_V*a_V
     272
     273#    for i in range(len(h_C)):
     274#        for j in range(2):
     275#            h_V[i,j] = w_V[i,j]-z_V[i,j]
     276#            if h_V[i,j]<h0:
     277#                h_V[i,j]=0
     278#                w_V[i,j]=z_V[i,j]
     279#            a_V[i,j] = b_V[i,j]*h_V[i,j]
     280#            d_V[i,j]=u_V[i,j]*a_V[i,j]
    275281   
    276282    return
     
    278284
    279285#--------------------------------------------------------
    280 #Boundaries - specific to the shallow_water_vel_domain
     286#Boundaries - specific to the channel_domain
    281287#--------------------------------------------------------
    282288class Reflective_boundary(Boundary):
     
    298304        #Handy shorthands
    299305        self.normals  = domain.normals
    300         self.area    = domain.quantities['area'].vertex_values
     306        self.area     = domain.quantities['area'].vertex_values
    301307        self.discharge     = domain.quantities['discharge'].vertex_values
    302308        self.bed      = domain.quantities['elevation'].vertex_values
     
    306312        self.stage    = domain.quantities['stage'].vertex_values
    307313
    308         from Numeric import zeros, Float
    309         #self.conserved_quantities = zeros(3, Float)
    310         self.evolved_quantities = zeros(7, Float)
     314        self.evolved_quantities = numpy.zeros(7, numpy.float)
    311315
    312316    def __repr__(self):
     
    320324
    321325        q = self.evolved_quantities
    322         q[0] = self.area[vol_id, edge_id]
     326        q[0] =  self.area[vol_id, edge_id]
    323327        q[1] = -self.discharge[vol_id, edge_id]
    324         q[2] = self.bed[vol_id, edge_id]
    325         q[3] = self.height[vol_id, edge_id]
     328        q[2] =  self.bed[vol_id, edge_id]
     329        q[3] =  self.height[vol_id, edge_id]
    326330        q[4] = -self.velocity[vol_id, edge_id]
    327         q[5] = self.width[vol_id,edge_id]
    328         q[6] = self.stage[vol_id,edge_id]
    329 
    330         #print "In Reflective q ",q
    331 
     331        q[5] =  self.width[vol_id,edge_id]
     332        q[6] =  self.stage[vol_id,edge_id]
    332333
    333334        return q
     
    348349            raise msg
    349350
    350         from Numeric import array, Float
    351         self.evolved_quantities=array(evolved_quantities).astype(Float)
     351        assert len(evolved_quantities) == 7
     352
     353        self.evolved_quantities=numpy.array(evolved_quantities,numpy.float)
    352354
    353355    def __repr__(self):
     
    370372
    371373
    372     Area     = domain.quantities['area']
     374    Area      = domain.quantities['area']
    373375    Discharge  = domain.quantities['discharge']
    374376    Elevation = domain.quantities['elevation']
     
    402404def boundary_stress(domain):
    403405
    404 
    405406    from util import gradient
    406407    from Numeric import zeros, Float, array, sum
    407 
    408 
    409408
    410409    Area     = domain.quantities['area']
     
    415414
    416415    discharge_ud  = Discharge.explicit_update
    417 
    418 
    419        
     416 
    420417    h = Height.vertex_values
    421418    b = Width.vertex_values
     
    451448
    452449    uh = domain.quantities['xmomentum'].centroid_values
    453     #vh = domain.quantities['ymomentum'].centroid_values
    454450    eta = domain.quantities['friction'].centroid_values
    455451
    456452    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
    457     #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
    458453
    459454    N = domain.number_of_elements
     
    502497
    503498
    504 def check_forcefield(f):
    505     """Check that f is either
    506     1: a callable object f(t,x,y), where x and y are vectors
    507        and that it returns an array or a list of same length
    508        as x and y
    509     2: a scalar
    510     """
    511 
    512     from Numeric import ones, Float, array
    513 
    514 
    515     if callable(f):
    516         #N = 3
    517         N = 2
    518         #x = ones(3, Float)
    519         #y = ones(3, Float)
    520         x = ones(2, Float)
    521         #y = ones(2, Float)
    522        
    523         try:
    524             #q = f(1.0, x=x, y=y)
    525             q = f(1.0, x=x)
    526         except Exception, e:
    527             msg = 'Function %s could not be executed:\n%s' %(f, e)
    528             #FIXME: Reconsider this semantics
    529             raise msg
    530 
    531         try:
    532             q = array(q).astype(Float)
    533         except:
    534             msg = 'Return value from vector function %s could ' %f
    535             msg += 'not be converted into a Numeric array of floats.\n'
    536             msg += 'Specified function should return either list or array.'
    537             raise msg
    538 
    539         #Is this really what we want?
    540         msg = 'Return vector from function %s ' %f
    541         msg += 'must have same lenght as input vectors'
    542         assert len(q) == N, msg
    543 
    544     else:
    545         try:
    546             f = float(f)
    547         except:
    548             msg = 'Force field %s must be either a scalar' %f
    549             msg += ' or a vector function'
    550             raise msg
    551     return f
    552499
    553500def linearb(domain):
Note: See TracChangeset for help on using the changeset viewer.