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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.