source: anuga_work/development/anuga_1d/anuga_1d_suggestion1_2_3/parabolic_cannal_sudi.py @ 7576

Last change on this file since 7576 was 7576, checked in by sudi, 14 years ago
File size: 10.2 KB
Line 
1import os
2from math import sqrt, pi, sin, cos
3from shallow_water_domain_suggestion3 import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6import time
7
8def analytic_cannal(C,t):
9    N = len(C)
10    u = zeros(N,Float)              ## water velocity
11    h = zeros(N,Float)              ## water depth
12    x = C
13    g = 9.81
14    ## Define Basin Bathymetry
15    z_b = zeros(N,Float)            ## elevation of basin
16    w = zeros(N,Float)              ## elevation of water surface
17    D0 = 10.0                       ## max equilibrium water depth at lowest point.
18    L_x = 2500.0                    ## width of channel
19    A0 = 0.5*L_x                    ## determines amplitudes of oscillations
20    omega = sqrt(2*g*D0)/L_x        ## angular frequency of osccilation
21    for i in range(N):
22        z_b[i] = D0*(x[i]**2/L_x**2)
23        u[i] = -A0*omega*sin(omega*t)
24        w[i] = D0+2*A0*D0/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
25        if w[i] <= z_b[i] :
26            u[i] = 0.0
27            w[i] = z_b[i]
28    h = w - z_b
29    T = 2.0*pi/omega
30    return u,h,w,z_b, T
31
32
33L_x = 2500.0                        # Length of channel (m)
34N = 400                               # Number of computational cells
35cell_len = 4*L_x/N                  # Origin = 0.0
36
37points = zeros(N+1,Float)
38for i in range(N+1):
39        points[i] = -2*L_x +i*cell_len
40
41domain = Domain(points)
42domain.order = 2
43domain.set_timestepping_method('rk2')
44domain.set_CFL(1.0)
45domain.beta = 1.0
46domain.set_limiter("minmod")
47
48def stage(x):
49    D0 = 10.0                       ## max equilibrium water depth at lowest point.
50    L_x = 2500.0                    ## width of channel
51    A0 = 0.5*L_x                    ## determines amplitudes of oscillations
52    omega = sqrt(2*g*D0)/L_x        ## angular frequency of osccilation
53    t=0.0
54    y = zeros(len(x),Float)
55    for i in range(len(x)):
56        y[i] = D0+2*A0*D0/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
57        #y[i] = 12.0
58    return y
59
60def elevation(x):
61    N = len(x)
62    D0 = 10.0
63    z = zeros(N,Float)
64    L_x = 2500.0
65    A0 = 0.5*L_x
66    omega = sqrt(2*g*D0)/L_x
67    for i in range(N):
68        z[i] = D0*(x[i]**2/L_x**2)
69    return z
70
71def height(x):
72    D0 = 10.0                       ## max equilibrium water depth at lowest point.
73    L_x = 2500.0                    ## width of channel
74    A0 = 0.5*L_x                    ## determines amplitudes of oscillations
75    omega = sqrt(2*g*D0)/L_x        ## angular frequency of osccilation
76    t=0.0
77    y = zeros(len(x),Float)
78    for i in range(len(x)):
79        y[i] = max(D0+2*A0*D0/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t)) - D0*(x[i]**2/L_x**2),  0.0)
80    return y
81
82
83domain.set_quantity('stage', stage)
84domain.set_quantity('elevation',elevation)
85#domain.set_quantity('height',height)
86domain.set_boundary({'exterior': Reflective_boundary(domain)})
87
88C = domain.centroids
89X = domain.vertices
90u,h,w,z_b,T = analytic_cannal(X.flat,domain.time)
91print 'T = ',T
92
93yieldstep = finaltime = T/16 #3.0*T/2.0
94StageQ = domain.quantities['stage'].vertex_values
95XmomQ = domain.quantities['xmomentum'].vertex_values
96
97import time
98while finaltime < T + 1.0:
99    yieldstep = finaltime
100    t0 = time.time()
101    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
102        domain.write_time()
103        #print "integral", domain.quantities['stage'].get_integral()
104        if t>0.0:
105            print "t=",t
106            uC,hC,wC,z_bC,TC = analytic_cannal(C,t)
107            for k in range(len(hC)):
108                if hC[k] < 0.0:
109                    hC[k] = 0.0
110                    wC[k] = z_bC[k]
111               
112            VelC = domain.quantities['velocity'].centroid_values
113            StageC = domain.quantities['stage'].centroid_values
114            ElevC = domain.quantities['elevation'].centroid_values
115            HeightC = domain.quantities['height'].centroid_values
116            XmomC = domain.quantities['xmomentum'].centroid_values
117
118           
119            #for i in range(domain.number_of_elements):
120            #    if HeightC[i] < 0.0 :
121            #        VelC[i] = 0.0
122               
123            #HeightC = domain.quantities['height'].centroid_values
124            #print 'Difference=',HeightC-hC
125            #print 'Difference=',StageC-wC
126       
127
128            error_h = (1.0/N)*sum(abs(hC-HeightC))
129            error_uh = (1.0/N)*sum(abs(uC*hC - XmomC))
130            error_u = (1.0/N)*sum(abs(uC - VelC))
131            #print 'len(hC)=',len(hC)
132            #print 'N=',N
133            print 'Height error measured at centroids = ', error_h
134            print 'Momentum error measured at centroids = ', error_uh
135            print 'Velocity error measured at centroids = ', error_u
136       
137            StageV = domain.quantities['stage'].vertex_values
138            ElevV = domain.quantities['elevation'].vertex_values
139            HeightV = domain.quantities['height'].vertex_values
140            XmomV = domain.quantities['xmomentum'].vertex_values
141            VelV = domain.quantities['velocity'].vertex_values
142       
143            u,h,w,z_b,T = analytic_cannal(X.flat,t)
144            for k in range(len(h)):
145                if h[k] < 0.0:
146                    h[k] = 0.0
147                    w[k] = z_b[k]
148       
149            from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
150            hold(False)
151       
152            plot1 = subplot(311)
153            plot(X,w,  X,StageV,  X,z_b)
154            #plot1.set_xlim([-6000,6000])
155            #xlabel('Position')
156            ylabel('Stage')
157            legend(('Analytical Solution', 'Numerical Solution', 'Canal Bed'),
158                   'upper center', shadow=False)
159       
160            plot2 = subplot(312)
161            plot(C,uC*hC, C,XmomC) #(X,u*h,X,XmomV)
162            #plot2.set_ylim([-1,25])
163            #xlabel('Position')
164            ylabel('Momentum')
165            legend(('Analytical Solution', 'Numerical Solution'),
166                   'upper left', shadow=False)
167
168       
169
170            #h_V = zeros(len(VelV),Float)
171            #print 'len(h_V)=',len(h_V)
172            #print 'len(VelV)=', len(VelV)
173            #print 'h_V=',h_V[:]
174            #h_V[:] = HeightV
175       
176            #for i in range(domain.number_of_elements):
177            #    for j in range(2):
178            #        if HeightV[i,j] <= 0.0:
179            #            VelV[i,j] = 0.0
180                   
181            #print 'len(h_V)=',len(h_V)
182            #print 'le(VelV)=', len(VelV)
183            #print 'h_V=',h_V[:]
184       
185            plot3 = subplot(313)
186            plot(C,uC, C,VelC) #(X,u, X,VelV)
187            xlabel('Position')
188            ylabel('Velocity')
189            legend(('Analytical Solution', 'Numerical Solution'),
190                   'upper right', shadow=False)
191            #print 'domain.order=',domain.order
192            #print 'domain.limiter=',domain.limiter
193            #print 'domain.time=',domain.time
194            #print 'domain.timestepping = ',domain.get_timestepping_method()
195
196            #file = "parabolic_new_"
197            #file += str(t)
198            #file += ".eps"
199            #savefig(file)
200            #show()
201    finaltime=finaltime + T/16#10.0 #T/20.0
202raw_input("Press the return key!")
203
204
205
206"""
207domain.set_quantity('stage', stage)
208domain.set_quantity('elevation',elevation)
209domain.set_boundary({'exterior': Reflective_boundary(domain)})
210
211C = domain.centroids
212X = domain.vertices
213u,h,w,z_b,T = analytic_cannal(X.flat,domain.time)
214print 'T = ',T
215
216StageQ = domain.quantities['stage'].vertex_values
217XmomQ = domain.quantities['xmomentum'].vertex_values
218
219print 'domain.order=',domain.order
220print 'domain.limiter=',domain.limiter
221print 'domain.time=',domain.time
222print 'domain.timestepping = ',domain.get_timestepping_method()
223
224yieldstep = finaltime = T
225t0 = time.time()
226while finaltime < T+1.0:
227    #yieldstep = finaltime
228    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
229        domain.write_time()
230        print "integral", domain.quantities['stage'].get_integral()
231        if t>= 0.0:       
232            uC,hC,wC,z_bC,TC = analytic_cannal(C,t)
233               
234            VelC = domain.quantities['velocity'].centroid_values
235            StageC = domain.quantities['stage'].centroid_values
236            ElevC = domain.quantities['elevation'].centroid_values
237            HeightC = domain.quantities['height'].centroid_values
238            XmomC = domain.quantities['xmomentum'].centroid_values
239           
240
241            error_h = cell_len*sum(abs(hC-HeightC))
242            error_uh = cell_len*sum(abs(uC*hC - XmomC))
243            error_u = cell_len*sum(abs(uC - VelC))
244
245            print 'Height error measured at centroids = ', error_h
246            print 'Momentum error measured at centroids = ', error_uh
247            print 'Velocity error measured at centroids = ', error_u
248       
249            StageV = domain.quantities['stage'].vertex_values
250            ElevV = domain.quantities['elevation'].vertex_values
251            HeightV = domain.quantities['height'].vertex_values
252            XmomV = domain.quantities['xmomentum'].vertex_values
253            VelV = domain.quantities['velocity'].vertex_values
254       
255            u,h,w,z_b,T = analytic_cannal(X.flat,t)
256       
257            from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
258            hold(False)
259       
260            plot1 = subplot(311)
261            plot(X,w, X,StageV,  X,ElevV)
262            #plot1.set_xlim([-6000,6000])
263            #xlabel('Position')
264            ylabel('Stage')
265            legend(('Analytical Solution', 'Numerical Solution', 'Canal Bed'),
266                   'upper center', shadow=False)
267       
268            plot2 = subplot(312)
269            plot(C,uC*hC, C,XmomC)
270            #plot2.set_ylim([-1,25])
271            #xlabel('Position')
272            ylabel('Momentum')
273            #legend(('Analytical Solution', 'Numerical Solution'),
274            #       'upper right', shadow=False)
275       
276            plot3 = subplot(313)
277            plot(C,uC, C,VelC)
278            #xlabel('Position')
279            ylabel('Velocity')
280
281            #plot4 = subplot(414)
282            #plot(X,h, X,HeightV)
283            #xlabel('Position')
284            #ylabel('Height')
285            #file = "parabolic_2nd_vanleer_"
286            #file += str(t)
287            #file += ".png"
288            #savefig(file)
289            #show()
290    finaltime=finaltime + T/16.0
291raw_input("Press the return key!")
292"""
Note: See TracBrowser for help on using the repository browser.