Changeset 7781


Ignore:
Timestamp:
Jun 5, 2010, 7:01:28 PM (15 years ago)
Author:
steve
Message:

cleaning up shallow water domain

File:
1 copied

Legend:

Unmodified
Added
Removed

  • anuga_work/development/pipeflow/sww_domain.py

    r7780 r7781  
    4545import numpy
    4646
    47 from domain import *
    48 Generic_Domain = Domain #Rename
     47from generic_domain import *
     48from sww_boundary_conditions import *
     49from sww_forcing_terms import *
    4950
    5051#Shallow water domain
    … …  
    6970        self.h0 = h0
    7071
    71         #forcing terms not included in 1d domain ?WHy?
     72        #forcing terms
    7273        self.forcing_terms.append(gravity)
    7374        #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
     75       
    8376       
    8477        #Stored output
    … …  
    131124       
    132125
    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
    147126
    148127    def check_integrity(self):
    … …  
    155134        assert self.conserved_quantities[1] == 'xmomentum', msg
    156135
    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)
     136
    161137
    162138    def compute_fluxes(self):
    163139        #Call correct module function
    164140        #(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)
     141
     142        import sys
     143
     144
     145        timestep = numpy.float(sys.maxint)
     146
     147        stage    = self.quantities['stage']
     148        xmom     = self.quantities['xmomentum']
     149        bed      = self.quantities['elevation']
     150        height   = self.quantities['height']
     151        velocity = self.quantities['velocity']
     152
     153
     154        from sww_comp_flux_ext import compute_fluxes_ext
     155
     156        self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity)
     157
    170158
    171159    def distribute_to_vertices_and_edges(self):
    … …  
    210198
    211199#=============== End of Shallow Water Domain ===============================
    212 
    213 #Rotation of momentum vector
    214 def 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 
    254 def 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 #
    343 def 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 #
    450 def 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
    512 def 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
    603 def 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 # ###################################
    623 def 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 
    706200
    707201
    … …  
    1019513
    1020514
    1021 ###############################################
    1022 #Boundaries - specific to the shallow water wave equation
    1023 class 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 
    1070 class 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 #
    1095 def 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  
    1137 def 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 
    1169 def 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 
    1203 def 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 
    1250 class 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):
    1369 def 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
     515
Note: See TracChangeset for help on using the changeset viewer.