source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/johns/run-program.py

Last change on this file was 7933, checked in by mungkasi, 14 years ago

Renaming a module "bisect.py" to "bisect_function.py" since it crashes with "bisect" in pylab.
File size: 11.4 KB
Line 
1from scipy import sin, cos, sqrt, linspace, pi, dot
2from Numeric import Float
3from numpy import zeros,array
4from gaussPivot import *
5from analytical_prescription import *
6from parameter import *
7import os, time, csv, pprint
8from domain_johns import *
9from config import g, epsilon
10from rootsearch import *
11from bisect_function import *
12
13#Analytical computations#################################################################
14def root_g(a,b,t):
15    dx = 0.01
16    def g(u):
17        return u + 2.0*A*pi/T*sin(2.0*pi/T*(t+u))
18    while 1:
19        x1,x2 = rootsearch(g,a,b,dx)
20        if x1 != None:
21            a = x2
22            root = bisect(g,x1,x2,1)
23        else:
24            break
25    return root
26def shore(t):
27    a = -1.0#-0.2#-1.0
28    b = 1.0#0.2#1.0
29    #dx = 0.01
30    u = root_g(a,b,t)
31    xi = -0.5*u*u + A*cos(2.0*pi/T*(t+u))
32    position = 1.0 + xi
33    return position, u
34
35
36
37
38#Numerical computations###################################################################
39def newtonRaphson2(f,q,tol=1.0e-15): ##1.0e-9 may be too large
40    for i in range(30):                                                                                                                                   
41        h = 1.0e-4      ##1.0e-4 may be too large.
42        n = len(q)
43        jac = zeros((n,n),Float)
44        if 1.0+q[0]-x<0.0:
45            temp1 = 1.0+q[0]-x
46            q[0] = q[0]-temp1
47            q[1] = v
48            return q
49        f0 = f(q)
50        for i in range(n):
51            temp = q[i]
52            q[i] = temp + h
53            f1 = f(q)
54            q[i] = temp
55            jac[:,i] = (f1 - f0)/h
56        if sqrt(dot(f0,f0)/len(q)) < tol: return q
57        dq = gaussPivot(jac,-f0)
58        q = q + dq
59        if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q
60    print 'Too many iterations'
61
62def elevation(X):
63    N = len(X)
64    z = zeros(N,Float)
65    for i in range(N):
66        z[i] = (h_0/L)*X[i] - h_0
67    return z
68
69def height(X):
70    N = len(X)
71    z = zeros(N,Float)
72    for i in range(N):
73        z[i] = h_0 - (h_0/L)*X[i]
74    return z
75
76def velocity(X):
77    N = len(X)
78    return zeros(N,Float)
79   
80boundary = { (0,0): 'left',(N-1,1): 'right'}
81
82domain = Domain(points,boundary)
83domain.order = 2
84domain.set_timestepping_method('rk2')
85domain.set_CFL(1.0)
86domain.beta = 1.0
87domain.set_limiter("minmod")
88
89
90def f_CG(t):
91    timing = t*sqrt(g*h_0)/L
92    w, u = prescribe(0.0,timing)
93    wO = w*h_0
94    uO = u*sqrt(g*h_0)
95    zO = -h_0
96    hO = wO - zO
97    pO = uO * hO
98    #[    'stage', 'xmomentum', 'elevation', 'height', 'velocity']
99    return [wO,  pO,   zO,      hO,      uO]
100def f_JOHNS(t):
101    timing = t*sqrt(g*h_0)/L
102    w, u = prescribe_at_O_JOHNS(timing)
103    wO = w*h_0
104    uO = u*sqrt(g*h_0)
105    zO = -h_0
106    hO = wO - zO
107    pO = uO * hO
108    #[    'stage', 'xmomentum', 'elevation', 'height', 'velocity']
109    return [wO,  pO,   zO,      hO,      uO]
110
111T1 = Time_boundary(domain,f_JOHNS)
112D2 = Dirichlet_boundary([50.5,  0.0,  50.5,  0.0,  0.0])
113domain.set_boundary({'left':T1,'right':D2})
114
115domain.set_quantity('height',height)
116domain.set_quantity('elevation',elevation)
117domain.set_quantity('velocity',velocity)
118
119
120Ver = domain.vertices
121n_V = len(Ver)
122AnalitW_V = zeros((n_V,2), Float)
123AnalitP_V = zeros((n_V,2), Float)
124AnalitZ_V = zeros((n_V,2), Float)
125AnalitH_V = zeros((n_V,2), Float)
126AnalitU_V = zeros((n_V,2), Float)
127
128Cen = domain.centroids
129n_C = len(Cen)
130AnalitW_C = zeros(n_C, Float)
131AnalitP_C = zeros(n_C, Float)
132AnalitZ_C = zeros(n_C, Float)
133AnalitH_C = zeros(n_C, Float)
134AnalitU_C = zeros(n_C, Float)
135
136waktu = 10.0     #3.0*60.0
137WAKTU = 10.0#12690.0  #Note: Tp=15.0*60.0
138yieldstep = finaltime = waktu
139t0 = time.time()
140counter=1
141
142shorelines_numerical_johns = zeros(int(WAKTU/waktu), Float)
143shorelines_analytical = zeros(int(WAKTU/waktu), Float)
144time_instants = zeros(int(WAKTU/waktu), Float)
145print "the initial time_instants=", time_instants
146
147
148while finaltime < WAKTU+0.1:
149    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
150        domain.write_time()
151    time_instants[counter-1] = domain.time
152   
153    Stage = domain.quantities['stage']
154    Momentum = domain.quantities['xmomentum']
155    Elevation = domain.quantities['elevation']
156    Height = domain.quantities['height']
157    Velocity = domain.quantities['velocity']
158
159    StageV = Stage.vertex_values
160    MomV = Momentum.vertex_values
161    ElevationV = Elevation.vertex_values
162    HeightV = Height.vertex_values
163    VelV = Velocity.vertex_values
164
165    StageC = Stage.centroid_values
166    MomC = Momentum.centroid_values
167    ElevationC = Elevation.centroid_values
168    HeightC = Height.centroid_values
169    VelC = Velocity.centroid_values
170
171    table = zeros((len(Ver.flat),6),Float)
172    for r in range(len(Ver.flat)):
173        for c in range(6):
174            if c==0:
175                table[r][c] = Ver.flat[r]
176            elif c==1:
177                table[r][c] = StageV.flat[r]
178            elif c==2:
179                table[r][c] = MomV.flat[r]
180            elif c==3:
181                table[r][c] = ElevationV.flat[r]
182            elif c==4:
183                table[r][c] = HeightV.flat[r]               
184            else:
185                table[r][c] = VelV.flat[r]
186
187    outname = "%s%04i%s%f%s" %("numerical_johns_", counter, "_", domain.time, ".csv")               
188    outfile = open(outname, 'w')
189    writer = csv.writer(outfile)
190    for row in table:
191        writer.writerow(row)       
192    outfile.close()   
193
194
195    for s in range(2*n_V):
196        heiL = HeightV.flat[s]
197        momR = MomV.flat[s+1]
198        if heiL >= 1e-6:
199            if abs(momR)==0.0: #<1e-15:
200                break
201    #print "s+1=",s+1
202    shorelines_numerical_johns[counter-1] = Ver.flat[s+1]
203    #print "shorelines_numerical_johns=",shorelines_numerical_johns[counter-1]
204
205
206   
207    #print "Now the ANALYTIC"
208    pos_shore, vel_shore = shore(domain.time*sqrt(g*h_0)/L) #dimensionless
209    pos_shore = pos_shore*L
210    vel_shore = vel_shore*sqrt(g*h_0)
211
212    shorelines_analytical[counter-1] = pos_shore
213
214    #The following is for calculating the error at centroids
215    for i in range(n_C):
216        x = Cen[i]
217        if x < pos_shore:
218            sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless
219            AnalitW_C[i] = sta*h_0
220            AnalitU_C[i] = vel                            #It needs dimensionalisation
221            AnalitZ_C[i] = (h_0/L)*x - h_0
222            AnalitH_C[i] = AnalitW_C[i] - AnalitZ_C[i]
223            AnalitP_C[i] = AnalitH_C[i]*AnalitU_C[i]  #It needs dimensionalisation
224        else:
225            AnalitW_C[i] = (h_0/L)*x - h_0                 
226            AnalitU_C[i] = 0.0
227            AnalitZ_C[i] = (h_0/L)*x - h_0
228            AnalitH_C[i] = 0.0
229            AnalitP_C[i] = 0.0
230    AnalitU_C = AnalitU_C*sqrt(g*h_0) #This is the dimensionalisation
231    AnalitP_C = AnalitP_C*sqrt(g*h_0) #This is the dimensionalisation
232
233    error_W = (1.0/n_C)*sum(abs(StageC-AnalitW_C))
234    error_P = (1.0/n_C)*sum(abs(MomC-AnalitP_C))
235    error_U = (1.0/n_C)*sum(abs(VelC-AnalitU_C))
236    table = zeros((1,4),Float)
237    for r in range(1):
238        for c in range(4):
239            if c==0:
240                table[r][c] = domain.time
241            elif c==1:
242                table[r][c] = error_W
243            elif c==2:
244                table[r][c] = error_P
245            else:
246                table[r][c] = error_U
247    outname = "%s%04i%s%f%s" %("error_johns_", counter, "_", domain.time, ".csv")               
248    outfile = open(outname, 'w')
249    writer = csv.writer(outfile)
250    for row in table:
251        writer.writerow(row)       
252    outfile.close()
253
254
255    #The following is for ploting the quantities at vertex values
256    for i in range(n_V):
257        vector_x = Ver[i]
258        for k in range(2):
259            x = vector_x[k]
260            if x < pos_shore:
261                sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless
262                AnalitW_V[i,k] = sta*h_0
263                AnalitU_V[i,k] = vel                            #It needs dimensionalisation
264                AnalitZ_V[i,k] = (h_0/L)*x - h_0
265                AnalitH_V[i,k] = AnalitW_V[i,k] - AnalitZ_V[i,k]
266                AnalitP_V[i,k] = AnalitH_V[i,k]*AnalitU_V[i,k]  #It needs dimensionalisation
267            else:
268                AnalitW_V[i,k] = (h_0/L)*x - h_0                 
269                AnalitU_V[i,k] = 0.0
270                AnalitZ_V[i,k] = (h_0/L)*x - h_0
271                AnalitH_V[i,k] = 0.0
272                AnalitP_V[i,k] = 0.0
273    AnalitU_V = AnalitU_V*sqrt(g*h_0) #This is the dimensionalisation
274    AnalitP_V = AnalitP_V*sqrt(g*h_0) #This is the dimensionalisation
275
276
277    table = zeros((len(Ver.flat),6),Float)
278    for r in range(len(Ver.flat)):
279        for c in range(6):
280            if c==0:
281                table[r][c] = Ver.flat[r]
282            elif c==1:
283                table[r][c] = AnalitW_V.flat[r]
284            elif c==2:
285                table[r][c] = AnalitP_V.flat[r]
286            elif c==3:
287                table[r][c] = AnalitZ_V.flat[r]
288            elif c==4:
289                table[r][c] = AnalitH_V.flat[r]               
290            else:
291                table[r][c] = AnalitU_V.flat[r]
292
293    outname = "%s%04i%s%f%s" %("analytical_", counter, "_", domain.time, ".csv")               
294    outfile = open(outname, 'w')
295    writer = csv.writer(outfile)
296    for row in table:
297        writer.writerow(row)       
298    outfile.close()
299
300    #put this
301    from pylab import clf,plot,title,xlabel,ylabel,legend,show,hold,subplot,ion#,savefig
302    hold(False)
303    clf()
304
305    plot1 = subplot(311)
306    plot(Ver.flat/1e+4,AnalitW_V.flat,'b-',  Ver.flat/1e+4,StageV.flat,'g-',  Ver.flat/1e+4,ElevationV.flat,'k-')
307    #plot(Ver,StageV, Ver,ElevationV)
308    #plot(Ver/L,StageV/h_0,  Ver/L,ElevationV/h_0)
309    #xlabel('Position')
310    ylabel('Stage')
311    #plot1.set_xlim([0.0,1.2])
312    plot1.set_ylim([-6.0,6.0])#([-9.0e-3,9.0e-3])
313    #legend(('Analytical Solution', 'Numerical Solution', 'Discretized Bed'),
314    #       'upper right', shadow=False)
315   
316    plot2 = subplot(312)
317    plot(Ver.flat/1e+4,AnalitP_V.flat,'b-',  Ver.flat/1e+4,MomV.flat,'g-')
318    #plot(Ver/L, VelV/sqrt(g*h_0))
319    #xlabel('Position')
320    ylabel('Momentum')
321    #plot2.set_xlim([0.0,1.2])
322    #legend(('Analytical Solution','Numerical Solution'),
323    #       'upper right', shadow=False)       
324
325    plot3 = subplot(313)
326    plot(Ver.flat/1e+4,AnalitU_V.flat,'b-',  Ver.flat/1e+4,VelV.flat,'g-')
327    #plot(Ver/L, VelV/sqrt(g*h_0))
328    xlabel('Position / 10,000')
329    ylabel('Velocity')
330    #plot2.set_xlim([0.0,1.2])
331    legend(('Analytical Solution','Numerical Solution'),
332           'upper center', shadow=False)
333   
334    #filename = "%s%04i%s%f%s" %("numerical_johns_", counter,"_", domain.time, ".png")
335    #savefig(filename)
336    show()
337    #raw_input("Press ENTER to continue")
338    counter = counter+1
339    finaltime = finaltime + waktu
340
341
342
343table = zeros((int(WAKTU/waktu), 2),Float)
344for r in range(int(WAKTU/waktu)):
345    for c in range(2):
346        if c==0:
347            table[r][c] = time_instants[r]               
348        else:
349            table[r][c] = shorelines_numerical_johns[r]
350
351outname = "%s" %("shore_numerical_johns.csv")               
352outfile = open(outname, 'w')
353writer = csv.writer(outfile)
354for row in table:
355    writer.writerow(row)       
356outfile.close()
357
358
359table = zeros((int(WAKTU/waktu), 2),Float)
360for r in range(int(WAKTU/waktu)):
361    for c in range(2):
362        if c==0:
363            table[r][c] = time_instants[r]               
364        else:
365            table[r][c] = shorelines_analytical[r]
366
367outname = "%s" %("shore_analytical.csv")               
368outfile = open(outname, 'w')
369writer = csv.writer(outfile)
370for row in table:
371    writer.writerow(row)       
372outfile.close()
373
Note: See TracBrowser for help on using the repository browser.