source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/cg/run-program.py @ 7924

Last change on this file since 7924 was 7924, checked in by mungkasi, 14 years ago
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_cg 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_CG)
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  #20.0*Tp  #Note: Tp=15.0*60.0
138yieldstep = finaltime = waktu
139t0 = time.time()
140counter=1
141
142shorelines_numerical_cg = zeros(int(WAKTU/waktu), Float)
143shorelines_analytical = zeros(int(WAKTU/waktu), Float)
144time_instants = zeros(int(WAKTU/waktu), Float)
145
146
147while finaltime < WAKTU+0.1:
148    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
149        domain.write_time()
150    time_instants[counter-1] = domain.time
151   
152    Stage = domain.quantities['stage']
153    Momentum = domain.quantities['xmomentum']
154    Elevation = domain.quantities['elevation']
155    Height = domain.quantities['height']
156    Velocity = domain.quantities['velocity']
157
158    StageV = Stage.vertex_values
159    MomV = Momentum.vertex_values
160    ElevationV = Elevation.vertex_values
161    HeightV = Height.vertex_values
162    VelV = Velocity.vertex_values
163   
164    StageC = Stage.centroid_values
165    MomC = Momentum.centroid_values
166    ElevationC = Elevation.centroid_values
167    HeightC = Height.centroid_values
168    VelC = Velocity.centroid_values
169
170    table = zeros((len(Ver.flat),6),Float)
171    for r in range(len(Ver.flat)):
172        for c in range(6):
173            if c==0:
174                table[r][c] = Ver.flat[r]
175            elif c==1:
176                table[r][c] = StageV.flat[r]
177            elif c==2:
178                table[r][c] = MomV.flat[r]
179            elif c==3:
180                table[r][c] = ElevationV.flat[r]
181            elif c==4:
182                table[r][c] = HeightV.flat[r]               
183            else:
184                table[r][c] = VelV.flat[r]
185
186    outname = "%s%04i%s%f%s" %("numerical_cg_", counter, "_", domain.time, ".csv")               
187    outfile = open(outname, 'w')
188    writer = csv.writer(outfile)
189    for row in table:
190        writer.writerow(row)       
191    outfile.close()   
192   
193
194    for s in range(2*n_V):
195        heiL = HeightV.flat[s]
196        momR = MomV.flat[s+1]
197        if heiL >= 1e-6:
198            if abs(momR)==0.0: #<1e-15:
199                break
200    #print "s+1=",s+1
201    shorelines_numerical_cg[counter-1] = Ver.flat[s+1]
202    #print "shorelines_numerical_cg=",shorelines_numerical_cg[counter-1]
203
204
205    #print "Now the ANALYTIC"
206    pos_shore, vel_shore = shore(domain.time*sqrt(g*h_0)/L) #dimensionless
207    pos_shore = pos_shore*L
208    vel_shore = vel_shore*sqrt(g*h_0)
209
210    shorelines_analytical[counter-1] = pos_shore
211
212    #The following is for calculating the error at centroids
213    for i in range(n_C):
214        x = Cen[i]
215        if x < pos_shore:
216            sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless
217            AnalitW_C[i] = sta*h_0
218            AnalitU_C[i] = vel                            #It needs dimensionalisation
219            AnalitZ_C[i] = (h_0/L)*x - h_0
220            AnalitH_C[i] = AnalitW_C[i] - AnalitZ_C[i]
221            AnalitP_C[i] = AnalitH_C[i]*AnalitU_C[i]  #It needs dimensionalisation
222        else:
223            AnalitW_C[i] = (h_0/L)*x - h_0                 
224            AnalitU_C[i] = 0.0
225            AnalitZ_C[i] = (h_0/L)*x - h_0
226            AnalitH_C[i] = 0.0
227            AnalitP_C[i] = 0.0
228    AnalitU_C = AnalitU_C*sqrt(g*h_0) #This is the dimensionalisation
229    AnalitP_C = AnalitP_C*sqrt(g*h_0) #This is the dimensionalisation
230
231    error_W = (1.0/n_C)*sum(abs(StageC-AnalitW_C))
232    error_P = (1.0/n_C)*sum(abs(MomC-AnalitP_C))
233    error_U = (1.0/n_C)*sum(abs(VelC-AnalitU_C))
234    table = zeros((1,4),Float)
235    for r in range(1):
236        for c in range(4):
237            if c==0:
238                table[r][c] = domain.time
239            elif c==1:
240                table[r][c] = error_W
241            elif c==2:
242                table[r][c] = error_P
243            else:
244                table[r][c] = error_U
245    outname = "%s%04i%s%f%s" %("error_cg_", counter, "_", domain.time, ".csv")               
246    outfile = open(outname, 'w')
247    writer = csv.writer(outfile)
248    for row in table:
249        writer.writerow(row)       
250    outfile.close()
251
252
253    #The following is for ploting the quantities at vertex values
254    for i in range(n_V):
255        vector_x = Ver[i]
256        for k in range(2):
257            x = vector_x[k]
258            if x < pos_shore:
259                sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless
260                AnalitW_V[i,k] = sta*h_0
261                AnalitU_V[i,k] = vel                            #It needs dimensionalisation
262                AnalitZ_V[i,k] = (h_0/L)*x - h_0
263                AnalitH_V[i,k] = AnalitW_V[i,k] - AnalitZ_V[i,k]
264                AnalitP_V[i,k] = AnalitH_V[i,k]*AnalitU_V[i,k]  #It needs dimensionalisation
265            else:
266                AnalitW_V[i,k] = (h_0/L)*x - h_0                 
267                AnalitU_V[i,k] = 0.0
268                AnalitZ_V[i,k] = (h_0/L)*x - h_0
269                AnalitH_V[i,k] = 0.0
270                AnalitP_V[i,k] = 0.0
271    AnalitU_V = AnalitU_V*sqrt(g*h_0) #This is the dimensionalisation
272    AnalitP_V = AnalitP_V*sqrt(g*h_0) #This is the dimensionalisation
273
274
275    table = zeros((len(Ver.flat),6),Float)
276    for r in range(len(Ver.flat)):
277        for c in range(6):
278            if c==0:
279                table[r][c] = Ver.flat[r]
280            elif c==1:
281                table[r][c] = AnalitW_V.flat[r]
282            elif c==2:
283                table[r][c] = AnalitP_V.flat[r]
284            elif c==3:
285                table[r][c] = AnalitZ_V.flat[r]
286            elif c==4:
287                table[r][c] = AnalitH_V.flat[r]               
288            else:
289                table[r][c] = AnalitU_V.flat[r]
290
291    outname = "%s%04i%s%f%s" %("analytical_", counter, "_", domain.time, ".csv")               
292    outfile = open(outname, 'w')
293    writer = csv.writer(outfile)
294    for row in table:
295        writer.writerow(row)       
296    outfile.close()
297
298   
299    from pylab import clf,plot,title,xlabel,ylabel,legend,show,hold,subplot,ion#,savefig
300    hold(False)
301    clf()
302
303    plot1 = subplot(311)
304    plot(Ver.flat/1e+4,AnalitW_V.flat,'b-',  Ver.flat/1e+4,StageV.flat,'g-',  Ver.flat/1e+4,ElevationV.flat,'k-')
305    #plot(Ver,StageV, Ver,ElevationV)
306    #plot(Ver/L,StageV/h_0,  Ver/L,ElevationV/h_0)
307    #xlabel('Position')
308    ylabel('Stage')
309    #plot1.set_xlim([0.0,1.2])
310    plot1.set_ylim([-6.0,6.0])#([-9.0e-3,9.0e-3])
311    #legend(('Analytical Solution', 'Numerical Solution', 'Discretized Bed'),
312    #       'upper right', shadow=False)
313   
314    plot2 = subplot(312)
315    plot(Ver.flat/1e+4,AnalitP_V.flat,'b-',  Ver.flat/1e+4,MomV.flat,'g-')
316    #plot(Ver/L, VelV/sqrt(g*h_0))
317    #xlabel('Position')
318    ylabel('Momentum')
319    #plot2.set_xlim([0.0,1.2])
320    #legend(('Analytical Solution','Numerical Solution'),
321    #       'upper right', shadow=False)       
322
323    plot3 = subplot(313)
324    plot(Ver.flat/1e+4,AnalitU_V.flat,'b-',  Ver.flat/1e+4,VelV.flat,'g-')
325    #plot(Ver/L, VelV/sqrt(g*h_0))
326    xlabel('Position / 10,000')
327    ylabel('Velocity')
328    #plot2.set_xlim([0.0,1.2])
329    legend(('Analytical Solution','Numerical Solution'),
330           'upper center', shadow=False)
331   
332    #filename = "%s%04i%s%f%s" %("numerical_cg_", counter,"_", domain.time, ".png")
333    #savefig(filename)
334    show()
335    #raw_input("Press ENTER to continue")
336   
337    counter = counter+1
338    finaltime = finaltime + waktu
339   
340
341table = zeros((int(WAKTU/waktu), 2),Float)
342for r in range(int(WAKTU/waktu)):
343    for c in range(2):
344        if c==0:
345            table[r][c] = time_instants[r]               
346        else:
347            table[r][c] = shorelines_numerical_cg[r]
348
349outname = "%s" %("shore_numerical_cg.csv")               
350outfile = open(outname, 'w')
351writer = csv.writer(outfile)
352for row in table:
353    writer.writerow(row)       
354outfile.close()
355
356
357table = zeros((int(WAKTU/waktu), 2),Float)
358for r in range(int(WAKTU/waktu)):
359    for c in range(2):
360        if c==0:
361            table[r][c] = time_instants[r]               
362        else:
363            table[r][c] = shorelines_analytical[r]
364
365outname = "%s" %("shore_analytical.csv")               
366outfile = open(outname, 'w')
367writer = csv.writer(outfile)
368for row in table:
369    writer.writerow(row)       
370outfile.close()
371
Note: See TracBrowser for help on using the repository browser.