Last change on this file since 7818 was 6453, checked in by steve, 15 years ago

File size: 4.4 KB
Line
1import os
2from math import sqrt, pi
3from channel_domain import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6
7
8h1 = 5.0
9h0 = 0.0
10
11## def analytical_sol(C,t):
12
13##     #t  = 0.0     # time (s)
14##     # gravity (m/s^2)
15##     #h1 = 10.0    # depth upstream (m)
16##     #h0 = 0.0     # depth downstream (m)
17##     L = 2000.0   # length of stream/domain (m)
18##     n = len(C)    # number of cells
19
20##     u = zeros(n,Float)
21##     h = zeros(n,Float)
22##     x = C-3*L/4.0
23
24
25##     for i in range(n):
26##         # Calculate Analytical Solution at time t > 0
27##         u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t)
28##         h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
29##         u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1))
30##         h3_ = 1.0/(9.0*g)*((x[i]+L/2.0)/t+2*sqrt(g*h1))*((x[i]+L/2.0)/t+2*sqrt(g*h1))
31
32##         if ( x[i] <= -1*L/2.0+2*(-sqrt(g*h1)*t)):
33##             u[i] = 0.0
34##             h[i] = h0
35##         elif ( x[i] <= -1*L/2.0-(-sqrt(g*h1)*t)):
36##             u[i] = u3_
37##             h[i] = h3_
38
39##         elif ( x[i] <= -t*sqrt(g*h1) ):
40##             u[i] = 0.0
41##             h[i] = h1
42##         elif ( x[i] <= 2.0*t*sqrt(g*h1) ):
43##             u[i] = u3
44##             h[i] = h3
45##         else:
46##             u[i] = 0.0
47##             h[i] = h0
48
49##     return h , u*h, u
50
51#def newLinePlot(title='Simple Plot'):
52#   import Gnuplot
53#    gg = Gnuplot.Gnuplot(persist=0)
54#    gg.terminal(postscript)
55#    gg.title(title)
56#    gg('set data style linespoints')
57#    gg.xlabel('x')
58#    gg.ylabel('y')
59#    return gg
60
61#def linePlot(gg,x1,y1,x2,y2):
62#    import Gnuplot
63#    plot1 = Gnuplot.PlotItems.Data(x1.flat,y1.flat,with="linespoints")
64#    plot2 = Gnuplot.PlotItems.Data(x2.flat,y2.flat, with="lines 3")
65#    g.plot(plot1,plot2)
66
67h2=10.0
68
69print "TEST 1D-SOLUTION III -- DRY BED"
70
71def stage(x):
72    return 5
73
74def width(x):
75    return 1
76
77def bot(x):
78    return x/500
79
80import time
81
82finaltime = 20.0
83yieldstep = finaltime
84L = 2000.0     # Length of channel (m)
85number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000]
86h_error = zeros(len(number_of_cells),Float)
87uh_error = zeros(len(number_of_cells),Float)
88k = 0
89for i in range(len(number_of_cells)):
90    N = int(number_of_cells[i])
91    print "Evaluating domain with %d cells" %N
92    cell_len = L/N # Origin = 0.0
93    points = zeros(N+1,Float)
94    for j in range(N+1):
95        points[j] = j*cell_len
96
97    domain = Domain(points)
98
99    domain.set_quantity('area', stage)
100    domain.set_quantity('width',width)
101    domain.set_quantity('elevation',bot)
102    domain.set_boundary({'exterior': Reflective_boundary(domain)})
103    domain.order = 2
104    domain.set_timestepping_method('rk2')
105    domain.set_CFL(1.0)
106    domain.set_limiter("vanleer")
107    #domain.h0=0.0001
108
109    t0 = time.time()
110
111    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
112        domain.write_time()
113
114    N = float(N)
115    HeightC = domain.quantities['height'].centroid_values
116    DischargeC = domain.quantities['discharge'].centroid_values
117    C = domain.centroids
118    #h, uh, u = analytical_sol(C,domain.time)
119    #h_error[k] = 1.0/(N)*sum(abs(h-StageC))
120    #uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
121    #print "h_error %.10f" %(h_error[k])
122    #print "uh_error %.10f"% (uh_error[k])
123    k = k+1
124    print 'That took %.2f seconds' %(time.time()-t0)
125    X = domain.vertices
126    HeightQ = domain.quantities['height'].vertex_values
127    VelocityQ = domain.quantities['velocity'].vertex_values
128    #stage = domain.quantities['stage'].vertex_values
129    #h, uh, u = analytical_sol(X.flat,domain.time)
130    x = X.flat
131    z =bot(x)
132    w=HeightQ.flat+z
133
134    #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
135    from pylab import *
136    import pylab as p
137    import matplotlib.axes3d as p3
138    print 'test 2'
139    hold(False)
140    print 'test 3'
141    plot1 = subplot(211)
142    print 'test 4'
143    plot(x,z,x,w)
144    print 'test 5'
145    plot1.set_ylim([-1,11])
146    xlabel('Position')
147    ylabel('Stage')
148    #legend(( 'Numerical Solution'),
150    plot2 = subplot(212)
151    plot(x,VelocityQ.flat)
152    plot2.set_ylim([-2,2])
153
154    xlabel('Position')
155    ylabel('Velocity')
156
157    file = "dry_bed_"
158    file += str(number_of_cells[i])
159    file += ".eps"
160    #savefig(file)
161    show()
162
163
164#print "Error in height", h_error
165#print "Error in xmom", uh_error
Note: See TracBrowser for help on using the repository browser.