Changeset 6167


Ignore:
Timestamp:
Jan 14, 2009, 4:10:50 PM (10 years ago)
Author:
wilson
Message:

update 14/01/09

File:
1 edited

Legend:

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

    r6038 r6167  
    11import os
    22from math import sqrt, pi
    3 from shallow_water_vel_domain import *
     3from channel_domain_Ab import *
    44from Numeric import allclose, array, zeros, ones, Float, take, sqrt
    55from config import g, epsilon
    66
    77
    8 h1 = 10.0
     8h1 = 5.0
    99h0 = 0.0
    1010
     
    6565#    g.plot(plot1,plot2)
    6666
    67 
     67h2=5.0
     68k=1
    6869
    6970print "TEST 1D-SOLUTION III -- DRY BED"
     
    7374    for i in range(len(x)):
    7475        if x[i]<=L/4.0:
    75             y[i] = h0
     76            y[i] = h0*width([x[i]])
    7677        elif x[i]<=3*L/4.0:
    77             y[i] = h1
     78            y[i] = h2*width([x[i]])
    7879        else:
    79             y[i] = h0
     80            y[i] = h0*width([x[i]])
    8081    return y
    8182
     83def width(x):
     84    return k
    8285
     86 
    8387import time
    8488
     
    8690yieldstep = finaltime
    8791L = 2000.0     # Length of channel (m)
    88 number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000]
    89 h_error = zeros(len(number_of_cells),Float)
    90 uh_error = zeros(len(number_of_cells),Float)
     92number_of_cells = [200]
    9193k = 0
    92 for i in range(len(number_of_cells)):
    93     N = int(number_of_cells[i])
    94     print "Evaluating domain with %d cells" %N
    95     cell_len = L/N # Origin = 0.0
    96     points = zeros(N+1,Float)
    97     for j in range(N+1):
    98         points[j] = j*cell_len
     94widths = [1]
     95heights= []
     96for i in range(len(widths)):
     97    k=widths[i]
     98    for i in range(len(number_of_cells)):
     99        N = int(number_of_cells[i])
     100        print "Evaluating domain with %d cells" %N
     101        cell_len = L/N # Origin = 0.0
     102        points = zeros(N+1,Float)
     103        for j in range(N+1):
     104            points[j] = j*cell_len
    99105       
    100     domain = Domain(points)
     106        domain = Domain(points)
    101107   
    102     domain.set_quantity('stage', stage)
    103     domain.set_boundary({'exterior': Reflective_boundary(domain)})
    104     domain.order = 2
    105     domain.set_timestepping_method('rk2')
    106     domain.set_CFL(1.0)
    107     domain.set_limiter("vanleer")
    108     #domain.h0=0.0001
     108        domain.set_quantity('area', stage)
     109        domain.set_quantity('width',width)
     110        print "width in cell 1",domain.quantities['width'].vertex_values[1]
     111        domain.set_boundary({'exterior': Reflective_boundary(domain)})
     112        domain.order = 2
     113        domain.set_timestepping_method('rk2')
     114        domain.set_CFL(1.0)
     115        domain.set_limiter("vanleer")
     116        #domain.h0=0.0001
     117 
     118        t0 = time.time()
    109119
    110     t0 = time.time()
     120        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     121            domain.write_time()
    111122
    112     for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    113         domain.write_time()
    114 
    115     N = float(N)
    116     StageC = domain.quantities['stage'].centroid_values
    117     XmomC = domain.quantities['xmomentum'].centroid_values
    118     C = domain.centroids
    119     h, uh, u = analytical_sol(C,domain.time)
    120     h_error[k] = 1.0/(N)*sum(abs(h-StageC))
    121     uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
    122     print "h_error %.10f" %(h_error[k])
    123     print "uh_error %.10f"% (uh_error[k])
    124     k = k+1
    125     print 'That took %.2f seconds' %(time.time()-t0)
    126     X = domain.vertices
    127     StageQ = domain.quantities['stage'].vertex_values
    128     XmomQ = domain.quantities['xmomentum'].vertex_values
    129     velQ = domain.quantities['velocity'].vertex_values
     123        N = float(N)
     124        HeightC = domain.quantities['height'].centroid_values
     125        DischargeC = domain.quantities['discharge'].centroid_values
     126        C = domain.centroids
     127        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))
     130        #print "h_error %.10f" %(h_error[k])
     131        #print "uh_error %.10f"% (uh_error[k])
     132        k = k+1
     133        print 'That took %.2f seconds' %(time.time()-t0)
     134        X = domain.vertices
     135        heights.append(domain.quantities['height'].vertex_values)
     136        VelocityQ = domain.quantities['velocity'].vertex_values
     137        #stage = domain.quantities['stage'].vertex_values
     138        h, uh, u = analytical_sol(X.flat,domain.time)
     139        x = X.flat
    130140   
    131     h, uh, u = analytical_sol(X.flat,domain.time)
    132     x = X.flat
     141        #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
     142from pylab import *
     143import pylab as p
     144import matplotlib.axes3d as p3
     145print 'test 2'
     146hold(False)
     147print 'test 3'
     148plot1 = subplot(211)
     149print 'test 4'
     150plot(x,h,x,heights[0].flat)
     151print 'test 5'
     152plot1.set_ylim([-1,11])
     153xlabel('Position')
     154ylabel('Stage')
     155legend(('Analytical Solution', 'Numerical Solution'),
     156           'upper right', shadow=True)
     157plot2 = subplot(212)
     158plot(x,u,x,VelocityQ.flat)
     159plot2.set_ylim([-35,35])
    133160   
    134     from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
    135     print 'test 2'
    136     hold(False)
    137     print 'test 3'
    138     plot1 = subplot(211)
    139     print 'test 4'
    140     plot(x,h,x,StageQ.flat)
    141     print 'test 5'
    142     plot1.set_ylim([-1,11])
    143     xlabel('Position')
    144     ylabel('Stage')
    145     legend(('Analytical Solution', 'Numerical Solution'),
    146            'upper right', shadow=True)
    147     plot2 = subplot(212)
    148     plot(x,u,x,velQ.flat)
    149     plot2.set_ylim([-35,35])
    150    
    151     xlabel('Position')
    152     ylabel('Velocity')
     161xlabel('Position')
     162ylabel('Velocity')
    153163   
    154     file = "dry_bed_"
    155     file += str(number_of_cells[i])
    156     file += ".eps"
    157     #savefig(file)
    158     show()
     164file = "dry_bed_"
     165file += str(number_of_cells[i])
     166file += ".eps"
     167#savefig(file)
     168show()
    159169   
    160170   
    161 print "Error in height", h_error
    162 print "Error in xmom", uh_error
     171#print "Error in height", h_error
     172#print "Error in xmom", uh_error
Note: See TracChangeset for help on using the changeset viewer.