Changeset 6453


Ignore:
Timestamp:
Mar 4, 2009, 3:30:29 PM (15 years ago)
Author:
steve
Message:

Added Padarn's (2008/2009 summer student) files

Location:
anuga_work/development/anuga_1d
Files:
19 added
5 edited

Legend:

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

    r6042 r6453  
    66consisting of methods specific to the Shallow Water Wave Equation
    77
     8This particular modification of the Domain class implements the ability to
     9vary the width of the 1D channel that the water flows in. As a result the
     10conserved variables are different than previous implementations and so are the
     11equations.
    812
    913U_t + E_x = S
    1014
    1115where
    12 
    13 U = [w, uh]
    14 E = [uh, u^2h + gh^2/2]
     16------------!!!! NOTE THIS NEEDS UPDATING !!!!------------------
     17U = [A, Q]
     18E = [Q, Q^2/A + gh^2/2]
    1519S represents source terms forcing the system
    1620(e.g. gravity, friction, wind stress, ...)
     
    3236
    3337The conserved quantities are w, uh
    34 
     38--------------------------------------------------------------------------
    3539For details see e.g.
    3640Christopher Zoppou and Stephen Roberts,
     
    3943
    4044
    41 John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
    42 Geoscience Australia, 2006
     45John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou,
     46Padarn Wilson, Geoscience Australia, 2008
    4347"""
    4448
  • anuga_work/development/anuga_1d/dam_padarn.py

    r6167 r6453  
    1010
    1111def analytical_sol(C,t):
    12    
     12    if t==0:
     13        t=0.0001
    1314    #t  = 0.0     # time (s)
    1415    # gravity (m/s^2)
     
    9091yieldstep = finaltime
    9192L = 2000.0     # Length of channel (m)
    92 number_of_cells = [200]
     93number_of_cells = [810]
    9394k = 0
    94 widths = [1]
     95widths = [1,2,5]
    9596heights= []
     97velocities = []
     98
    9699for i in range(len(widths)):
    97100    k=widths[i]
     
    118121        t0 = time.time()
    119122
    120         for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     123        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):           
    121124            domain.write_time()
    122125
     
    126129        C = domain.centroids
    127130        h, uh, u = analytical_sol(C,domain.time)
    128         #h_error[k] = 1.0/(N)*sum(abs(h-StageC))
    129         #uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
     131        h_error = 1.0/(N)*sum(abs(h-HeightC))
     132        u_error = 1.0/(N)*sum(abs(uh-DischargeC))
    130133        #print "h_error %.10f" %(h_error[k])
    131134        #print "uh_error %.10f"% (uh_error[k])
     
    134137        X = domain.vertices
    135138        heights.append(domain.quantities['height'].vertex_values)
    136         VelocityQ = domain.quantities['velocity'].vertex_values
     139        velocities.append( domain.quantities['velocity'].vertex_values)
    137140        #stage = domain.quantities['stage'].vertex_values
    138141        h, uh, u = analytical_sol(X.flat,domain.time)
    139142        x = X.flat
    140    
     143           
     144        print "Error in height", h_error
     145        print "Error in xmom", u_error
    141146        #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
    142147from pylab import *
     
    144149import matplotlib.axes3d as p3
    145150print 'test 2'
    146 hold(False)
     151#hold(False)
    147152print 'test 3'
    148153plot1 = subplot(211)
    149154print 'test 4'
    150 plot(x,h,x,heights[0].flat)
     155plot(x,heights[0].flat,x,h)
    151156print 'test 5'
    152157plot1.set_ylim([-1,11])
    153158xlabel('Position')
    154159ylabel('Stage')
    155 legend(('Analytical Solution', 'Numerical Solution'),
    156            'upper right', shadow=True)
     160legend(('Numerical Solution','Analytic Soltuion'), 'upper right', shadow=True)
    157161plot2 = subplot(212)
    158 plot(x,u,x,VelocityQ.flat)
     162
     163
     164plot(x,velocities[0].flat,x,u)
    159165plot2.set_ylim([-35,35])
    160    
    161166xlabel('Position')
    162167ylabel('Velocity')
     168
     169print heights[0].flat-heights[1].flat
    163170   
    164171file = "dry_bed_"
     
    168175show()
    169176   
    170    
    171 #print "Error in height", h_error
    172 #print "Error in xmom", uh_error
     177
  • anuga_work/development/anuga_1d/domain.py

    r6042 r6453  
    916916        # Update ghosts
    917917        self.update_ghosts()
    918 
     918       
    919919        # Initial update of vertex and edge values
    920         self.distribute_to_vertices_and_edges()
    921 
     920        self.distribute_to_vertices_and_edges() 
     921       
    922922        # Update extrema if necessary (for reporting)
    923923        self.update_extrema()
  • anuga_work/development/anuga_1d/dry_dam_vel.py

    r6452 r6453  
    7878for j in range(N+1):
    7979    points[j] = j*cell_len
    80 
    81 boundary = {}
    82 boundary[0,0] = 'left'
    83 boundary[N-1,1] = 'right'
    84 
    85 domain = Domain(points, boundary = boundary)
    86 
     80       
     81domain = Domain(points)
     82   
    8783domain.set_quantity('stage', stage)
    88 
    89 Br =  Reflective_boundary(domain)
    90 domain.set_boundary({'left':  Br, 'right': Br})
     84domain.set_boundary({'exterior': Reflective_boundary(domain)})
    9185domain.order = 2
    9286domain.set_timestepping_method('euler')
  • anuga_work/development/anuga_1d/quantity.py

    r5858 r6453  
    135135            #Use function specific method
    136136            self.set_function_values(X, location)
     137         
    137138        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
    138139            if location == 'centroids':
     
    163164       
    164165        if location == 'centroids':
    165             
     166         
    166167            P = self.domain.centroids
    167168            self.set_values(f(P), location)
    168169        else:
    169170            #Vertices
     171           
    170172            P = self.domain.get_vertices()
    171173           
    172174            for i in range(2):
    173                 self.vertex_values[:,i] = f(P[:,i])       
    174 
     175               
     176                self.vertex_values[:,i] = f(P[:,i])
     177               
    175178    def set_array_values(self, values, location='vertices'):
    176179        """Set values for quantity
Note: See TracChangeset for help on using the changeset viewer.