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

Last change on this file since 7922 was 7922, checked in by mungkasi, 13 years ago

Modifying codes so that arrays are compatible with numpy instead of Numeric.
File size: 11.3 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 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,savefig,show,hold,subplot,ion
302    hold(False)
303    clf()
304
305    plot1 = subplot(311)
306    plot(Ver/1e+4,AnalitW_V,'b-',  Ver/1e+4,StageV,'g-',  Ver/1e+4,ElevationV,'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/1e+4,AnalitP_V,'b-',  Ver/1e+4,MomV,'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/1e+4,AnalitU_V,'b-',  Ver/1e+4,VelV,'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    #put this
339    counter = counter+1
340    finaltime = finaltime + waktu
341
342
343
344table = zeros((int(WAKTU/waktu), 2),Float)
345for r in range(int(WAKTU/waktu)):
346    for c in range(2):
347        if c==0:
348            table[r][c] = time_instants[r]               
349        else:
350            table[r][c] = shorelines_numerical_johns[r]
351
352outname = "%s" %("shore_numerical_johns.csv")               
353outfile = open(outname, 'w')
354writer = csv.writer(outfile)
355for row in table:
356    writer.writerow(row)       
357outfile.close()
358
359
360table = zeros((int(WAKTU/waktu), 2),Float)
361for r in range(int(WAKTU/waktu)):
362    for c in range(2):
363        if c==0:
364            table[r][c] = time_instants[r]               
365        else:
366            table[r][c] = shorelines_analytical[r]
367
368outname = "%s" %("shore_analytical.csv")               
369outfile = open(outname, 'w')
370writer = csv.writer(outfile)
371for row in table:
372    writer.writerow(row)       
373outfile.close()
374
Note: See TracBrowser for help on using the repository browser.