Last change on this file since 6167 was 6167, checked in by wilson, 15 years ago

update 14/01/09

File size: 4.5 KB
Line
1import os
2from math import sqrt, pi
3from channel_domain_Ab 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
11def 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=5.0
68k=1
69
70print "TEST 1D-SOLUTION III -- DRY BED"
71
72def stage(x):
73    y = zeros(len(x),Float)
74    for i in range(len(x)):
75        if x[i]<=L/4.0:
76            y[i] = h0*width([x[i]])
77        elif x[i]<=3*L/4.0:
78            y[i] = h2*width([x[i]])
79        else:
80            y[i] = h0*width([x[i]])
81    return y
82
83def width(x):
84    return k
85
86
87import time
88
89finaltime = 10.0
90yieldstep = finaltime
91L = 2000.0     # Length of channel (m)
92number_of_cells = [200]
93k = 0
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
105
106        domain = Domain(points)
107
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()
119
120        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
121            domain.write_time()
122
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
140
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'),
157plot2 = subplot(212)
158plot(x,u,x,VelocityQ.flat)
159plot2.set_ylim([-35,35])
160
161xlabel('Position')
162ylabel('Velocity')
163
164file = "dry_bed_"
165file += str(number_of_cells[i])
166file += ".eps"
167#savefig(file)
168show()
169
170
171#print "Error in height", h_error
172#print "Error in xmom", uh_error
Note: See TracBrowser for help on using the repository browser.