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

Last change on this file since 7837 was 7837, checked in by mungkasi, 12 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

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