Changeset 5536


Ignore:
Timestamp:
Jul 18, 2008, 4:37:03 PM (11 years ago)
Author:
steve
Message:

Added unit test files for quantity and shallow water

Location:
anuga_work/development/anuga_1d
Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/anuga_1d/domain.py

    r5535 r5536  
    3434        self.beta = 1.0
    3535        self.limiter = "minmod_kurganov"
    36         self.split = False
    3736        self.wet_nodes = zeros((N,2), Int) # should this be here
    3837
     
    125124        #FIXME: remove later - maybe OK, though....
    126125        for name in self.conserved_quantities:
    127             self.quantities[name] = Conserved_quantity(self)
     126            self.quantities[name] = Quantity(self)
    128127        for name in self.other_quantities:
    129128            self.quantities[name] = Quantity(self)
  • anuga_work/development/anuga_1d/dry_dam_sudi.py

    r5535 r5536  
    11import os
    22from math import sqrt, pi
    3 from shallow_water_1d import *
     3from shallow_water_domain import *
    44from Numeric import allclose, array, zeros, ones, Float, take, sqrt
    55from config import g, epsilon
  • anuga_work/development/anuga_1d/quantity.py

    r5535 r5536  
    2222        #from domain import Domain
    2323        #from domain_order2 import Domain
    24         from domain_t2 import Domain
     24        from domain import Domain
    2525        from Numeric import array, zeros, Float
    2626
     
    5656        #Intialise centroid values
    5757        self.interpolate()
     58
     59
     60        from Numeric import zeros, Float
     61       
     62        #Allocate space for boundary values
     63        #L = len(domain.boundary)
     64        self.boundary_values = zeros(2, Float) #assumes no parrellism
     65
     66        #Allocate space for updates of conserved quantities by
     67        #flux calculations and forcing functions
     68       
     69        N = domain.number_of_elements
     70        self.explicit_update = zeros(N, Float )
     71        self.semi_implicit_update = zeros(N, Float )
     72
     73        self.gradients = zeros(N, Float)
     74        self.qmax = zeros(self.centroid_values.shape, Float)
     75        self.qmin = zeros(self.centroid_values.shape, Float)
     76
     77        self.beta = domain.beta       
    5878
    5979    #Methods for operator overloading
     
    381401        return integral
    382402
    383 class Conserved_quantity(Quantity):
    384     """Class conserved quantity adds to Quantity:
    385 
    386     storage and method for updating, and
    387     methods for extrapolation from centropid to vertices inluding
    388     gradients and limiters
    389     """
    390 
    391     def __init__(self, domain, vertex_values=None):
    392         Quantity.__init__(self, domain, vertex_values)
    393 
    394         from Numeric import zeros, Float
    395        
    396         #Allocate space for boundary values
    397         #L = len(domain.boundary)
    398         self.boundary_values = zeros(2, Float) #assumes no parrellism
    399 
    400         #Allocate space for updates of conserved quantities by
    401         #flux calculations and forcing functions
    402        
    403         N = domain.number_of_elements
    404         self.explicit_update = zeros(N, Float )
    405         self.semi_implicit_update = zeros(N, Float )
    406 
    407         self.gradients = zeros(N, Float)
    408         self.qmax = zeros(self.centroid_values.shape, Float)
    409         self.qmin = zeros(self.centroid_values.shape, Float)
    410 
    411         self.beta = domain.beta
    412403
    413404    def update(self, timestep):
     
    863854                else:
    864855                    qv[k,i] = qc[k]
     856   
     857
     858class Conserved_quantity(Quantity):
     859    """Class conserved quantity adds to Quantity:
     860
     861    storage and method for updating, and
     862    methods for extrapolation from centropid to vertices inluding
     863    gradients and limiters
     864    """
     865
     866    def __init__(self, domain, vertex_values=None):
     867        Quantity.__init__(self, domain, vertex_values)
     868
     869        print "Use Quantity instead of Conserved_quantity"
     870
    865871                   
    866872
     
    884890if __name__ == "__main__":
    885891    #from domain import Domain
    886     from shallow_water_1d import Domain     
     892    from shallow_water_domain import Domain     
    887893    from Numeric import arange
    888894   
     
    892898    D1 = Domain(points1)
    893899
    894     Q1 = Conserved_quantity(D1, vertex_values)
     900    Q1 = Quantity(D1, vertex_values)
    895901
    896902    print Q1.vertex_values
     
    941947    D2 = Domain(points2)
    942948
    943     Q2 = Conserved_quantity(D2)
     949    Q2 = Quantity(D2)
    944950    Q2.set_values(fun,'vertices')
    945951    Xc = Q2.domain.vertices
  • anuga_work/development/anuga_1d/shallow_water_domain.py

    r5535 r5536  
    4545#from domain import *
    4646#from domain_order2 import *
    47 from domain_t2 import *
     47from domain import *
    4848Generic_Domain = Domain #Rename
    4949
     
    388388
    389389
    390 def flux_function_split(normal, ql, qr, zl, zr):
    391     from config import g, epsilon
    392     from math import sqrt
    393     from Numeric import array
    394 
    395     #print 'ql',ql
    396 
    397     #Align momentums with x-axis
    398     #q_left  = rotate(ql, normal, direction = 1)
    399     #q_right = rotate(qr, normal, direction = 1)
    400     q_left = ql
    401     q_left[1] = q_left[1]*normal
    402     q_right = qr
    403     q_right[1] = q_right[1]*normal
    404 
    405     #z = (zl+zr)/2 #Take average of field values
    406     z = 0.5*(zl+zr) #Take average of field values
    407 
    408     w_left  = q_left[0]   #w=h+z
    409     h_left  = w_left-z
    410     uh_left = q_left[1]
    411 
    412     if h_left < epsilon:
    413         u_left = 0.0  #Could have been negative
    414         h_left = 0.0
    415     else:
    416         u_left  = uh_left/h_left
    417 
    418 
    419     w_right  = q_right[0]  #w=h+z
    420     h_right  = w_right-z
    421     uh_right = q_right[1]
    422 
    423 
    424     if h_right < epsilon:
    425         u_right = 0.0  #Could have been negative
    426         h_right = 0.0
    427     else:
    428         u_right  = uh_right/h_right
    429 
    430     #vh_left  = q_left[2]
    431     #vh_right = q_right[2]
    432 
    433     #soundspeed_left  = sqrt(g*h_left)
    434     #soundspeed_right = sqrt(g*h_right)
    435 
    436     #Maximal wave speed
    437     #s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
    438     s_max = max(u_left, u_right, 0)
    439    
    440     #Minimal wave speed
    441     #s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
    442     s_min = min(u_left, u_right, 0)
    443    
    444     #Flux computation
    445 
    446     #flux_left  = array([u_left*h_left,
    447     #                    u_left*uh_left + 0.5*g*h_left*h_left])
    448     #flux_right = array([u_right*h_right,
    449     #                    u_right*uh_right + 0.5*g*h_right*h_right])
    450     flux_left  = array([u_left*h_left,
    451                         u_left*uh_left])# + 0.5*g*h_left*h_left])
    452     flux_right = array([u_right*h_right,
    453                         u_right*uh_right])# + 0.5*g*h_right*h_right])
    454    
    455     denom = s_max-s_min
    456     if denom == 0.0:
    457         edgeflux = array([0.0, 0.0])
    458         max_speed = 0.0
    459     else:
    460         edgeflux = (s_max*flux_left - s_min*flux_right)/denom
    461         edgeflux += s_max*s_min*(q_right-q_left)/denom
    462        
    463         edgeflux[1] = edgeflux[1]*normal
    464 
    465         max_speed = max(abs(s_max), abs(s_min))
    466 
    467     return edgeflux, max_speed   
     390 
    468391
    469392def compute_timestep(domain):
     
    518441            normal = domain.normals[k, i]
    519442
    520             if domain.split == False:
    521                 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
    522             elif domain.split == True:
    523                 edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
     443            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
     444
    524445            #Update optimal_timestep
    525446            try:
     
    612533                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
    613534                qr[0] = stage[n, m]
    614                 qr[1] =  xmom[n, m]
     535                qr[1] = xmom[n, m]
    615536                zr = bed[n, m]
    616537
     
    625546           
    626547
    627             if domain.split == False:
    628                 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
    629             elif domain.split == True:
    630                 edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
     548            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
     549
    631550            #print 'edgeflux', edgeflux
    632551
     
    646565
    647566        Stage.explicit_update[k] = flux[0]
    648         Xmom.explicit_update[k] = flux[1]
     567        Xmom.explicit_update[k] =  flux[1]
    649568        #Ymom.explicit_update[k] = flux[2]
    650569        #print "flux cell",k,flux[0]
     
    1054973       
    1055974        #Update momentum (explicit update is reset to source values)
    1056         if domain.split == False:
    1057             xmom[k] += -g*bx*avg_h
    1058             #xmom[k] = -g*bx*avg_h
    1059             #stage[k] = 0.0
    1060         elif domain.split == True:
    1061             xmom[k] += -g*wx*avg_h
    1062             #xmom[k] = -g*wx*avg_h
    1063         #ymom[k] += -g*zy*avg_h
     975        xmom[k] += -g*bx*avg_h
     976        #xmom[k] = -g*bx*avg_h
     977        #stage[k] = 0.0
     978 
    1064979 
    1065980def manning_friction(domain):
Note: See TracChangeset for help on using the changeset viewer.