Ignore:
Timestamp:
Apr 13, 2006, 4:34:22 PM (18 years ago)
Author:
jakeman
Message:

Adding new files

File:
1 edited

Legend:

Unmodified
Added
Removed

  • development/pyvolution-1d/shallow_water_1d.py

    r2705 r2716  
    4949class Domain(Generic_Domain):
    5050
    51     def __init__(self, coordinates, boundary = None, tagged_elements = None):
     51    def __init__(self, coordinates, boundary = None, tagged_elements = None,
     52                 geo_reference = None):
    5253
    5354        conserved_quantities = ['stage', 'xmomentum']
    … …  
    5556        Generic_Domain.__init__(self, coordinates, boundary,
    5657                                conserved_quantities, other_quantities,
    57                                 tagged_elements)
     58                                tagged_elements, geo_reference)
    5859       
    5960        from config import minimum_allowed_height, g
    … …  
    281282# Flux computation
    282283#def flux_function(normal, ql, qr, zl, zr):
    283 def flux_function(ql, qr, zl, zr):
     284def flux_function(normal, ql, qr, zl, zr):
    284285    """Compute fluxes between volumes for the shallow water wave equation
    285286    cast in terms of w = h+z using the 'central scheme' as described in
    … …  
    305306    #q_right = rotate(qr, normal, direction = 1)
    306307    q_left = ql
     308    q_left[1] = q_left[1]*normal
    307309    q_right = qr
     310    q_right[1] = q_right[1]*normal
    308311
    309312    z = (zl+zr)/2 #Take average of field values
    … …  
    338341
    339342    #Maximal wave speed
    340     #s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
    341     s_max = max(u_left+soundspeed_left, u_right+soundspeed_right,0)
     343    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
    342344   
    343345    #Minimal wave speed
    344     #s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
    345     s_min = min(u_left-soundspeed_left, u_right-soundspeed_right,0)
     346    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
    346347   
    347348    #Flux computation
    … …  
    363364    denom = s_max-s_min
    364365    if denom == 0.0:
    365         edgeflux = array([0.0, 0.0, 0.0])
     366        #edgeflux = array([0.0, 0.0, 0.0])
    366367        edgeflux = array([0.0, 0.0])
    367368        max_speed = 0.0
    … …  
    369370        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
    370371        edgeflux += s_max*s_min*(q_right-q_left)/denom
     372       
     373        edgeflux[1] = -edgeflux[1]*normal
    371374
    372375        #edgeflux = rotate(edgeflux, normal, direction=-1)
    … …  
    407410
    408411    #Arrays
    409     stage = Stage.edge_values
    410     xmom =  Xmom.edge_values
     412    #stage = Stage.edge_values
     413    #xmom =  Xmom.edge_values
    411414 #   ymom =  Ymom.edge_values
    412     bed =   Bed.edge_values
     415    #bed =   Bed.edge_values
     416   
     417    stage = Stage.vertex_values
     418    xmom =  Xmom.vertex_values
     419    bed =   Bed.vertex_values
     420   
     421    print 'stage edge values'
     422    print stage
     423    print 'xmom edge values'
     424    print xmom
     425    print 'bed values'
     426    print bed
    413427
    414428    stage_bdry = Stage.boundary_values
    415429    xmom_bdry =  Xmom.boundary_values
     430    print 'stage_bdry'
     431    print stage_bdry
     432    print 'xmom_bdry'
     433    print xmom_bdry
    416434#    ymom_bdry =  Ymom.boundary_values
    417435
    … …  
    439457                zr = zl #Extend bed elevation to boundary
    440458            else:
    441                 m = domain.neighbour_edges[k,i]
     459                #m = domain.neighbour_edges[k,i]
     460                m = domain.neighbour_vertices[k,i]
    442461                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
    443462                qr = [stage[n, m], xmom[n, m]]
    … …  
    446465
    447466            #Outward pointing normal vector
    448             normal = domain.normals[k, 2*i:2*i+2]
     467            #normal = domain.normals[k, 2*i:2*i+2]
    449468            #CHECK THIS LINE [k, 2*i:2*i+1]
    450469
    451470            #Flux computation using provided function
    452471            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
    453             edgeflux, max_speed = flux_function(ql, qr, zl, zr)
    454             flux -= edgeflux * domain.edgelengths[k,i]
     472            print 'ql'
     473            print ql
     474            print 'qr'
     475            print qr
     476           
     477            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
     478            print 'edgeflux'
     479            print edgeflux
     480            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
     481            # flux = edgefluxleft - edgefluxright
     482            flux -= edgeflux #* domain.edgelengths[k,i]
    455483
    456484            #Update optimal_timestep
    457485            try:
    458                 timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
     486                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
     487                timestep = min(timestep, 0.5/max_speed)
     488                print 'check time step in compute fluxes is ok'
    459489            except ZeroDivisionError:
    460490                pass
    … …  
    514544    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
    515545                                     domain.neighbours,
    516                                      domain.neighbour_edges,
     546                                     #domain.neighbour_edges,
     547                                     domain.neighbour_vertices,
    517548                                     domain.normals,
    518                                      domain.edgelengths,
    519                                      domain.radii,
     549                                     #domain.edgelengths,
     550                                     #domain.radii,
    520551                                     domain.areas,
    521                                      Stage.edge_values,
    522                                      Xmom.edge_values,
     552                                     #Stage.edge_values,
     553                                     #Xmom.edge_values,
    523554                                     #Ymom.edge_values,
    524                                      Bed.edge_values,
     555                                     #Bed.edge_values,
     556                                     Stage.vertex_values,
     557                                     Xmom.vertex_values,
     558                                     Bed.vertex_values,
    525559                                     Stage.boundary_values,
    526560                                     Xmom.boundary_values,
    … …  
    568602        #MH090605 if second order,
    569603        #perform the extrapolation and limiting on
    570         #all of the conserved quantitie
     604        #all of the conserved quantities
    571605
    572606        if (domain.order == 1):
    … …  
    595629
    596630    #Compute edge values by interpolation
    597     for name in domain.conserved_quantities:
    598         Q = domain.quantities[name]
    599         Q.interpolate_from_vertices_to_edges()
     631    #for name in domain.conserved_quantities:
     632    #    Q = domain.quantities[name]
     633    #    Q.interpolate_from_vertices_to_edges()
    600634
    601635
    … …  
    898932
    899933        #Handy shorthands
    900         self.stage   = domain.quantities['stage'].edge_values
    901         self.xmom    = domain.quantities['xmomentum'].edge_values
     934        #self.stage   = domain.quantities['stage'].edge_values
     935        #self.xmom    = domain.quantities['xmomentum'].edge_values
    902936        #self.ymom    = domain.quantities['ymomentum'].edge_values
    903937        #self.normals = domain.normals
     938        self.stage   = domain.quantities['stage'].vertex_values
     939        self.xmom    = domain.quantities['xmomentum'].vertex_values
    904940
    905941        from Numeric import zeros, Float
    … …  
    931967        r[1] = -q[1]
    932968        q = r
     969        #For start interval there is no outward momentum so do not need to
     970        #reverse direction in this case
    933971
    934972        return q
    … …  
    950988    Stage = domain.quantities['stage']
    951989    Elevation = domain.quantities['elevation']
    952     h = Stage.edge_values - Elevation.edge_values
     990    #h = Stage.edge_values - Elevation.edge_values
     991    h = Stage.vertex_values - Elevation.vertex_values
    953992    v = Elevation.vertex_values
    954993
    … …  
    10051044
    10061045    uh = domain.quantities['xmomentum'].centroid_values
    1007     vh = domain.quantities['ymomentum'].centroid_values
     1046    #vh = domain.quantities['ymomentum'].centroid_values
    10081047    eta = domain.quantities['friction'].centroid_values
    10091048
    10101049    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
    1011     ymom_update = domain.quantities['ymomentum'].semi_implicit_update
     1050    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
    10121051
    10131052    N = domain.number_of_elements
    … …  
    10181057        if eta[k] >= eps:
    10191058            if h[k] >= eps:
    1020                 S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
     1059                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
     1060                S = -g * eta[k]**2 * uh[k]
    10211061                S /= h[k]**(7.0/3)
    10221062
    10231063                #Update momentum
    10241064                xmom_update[k] += S*uh[k]
    1025                 ymom_update[k] += S*vh[k]
     1065                #ymom_update[k] += S*vh[k]
    10261066
    10271067
Note: See TracChangeset for help on using the changeset viewer.