source: trunk/anuga_work/development/2010-projects/anuga_1d/channel/dam_padarn.py @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 12 years ago

Moving 2010 project

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