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

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

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

File size: 11.7 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_johns 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_JOHNS)
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  #Note: Tp=15.0*60.0
137yieldstep = finaltime = waktu
138t0 = time.time()
139counter=1
140
141shorelines_numerical_johns = zeros(int(WAKTU/waktu), Float)
142shorelines_analytical = zeros(int(WAKTU/waktu), Float)
143time_instants = zeros(int(WAKTU/waktu), Float)
144print "the initial time_instants=", time_instants
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_johns_", 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_johns[counter-1] = Ver.flat[s+1]
202    #print "shorelines_numerical_johns=",shorelines_numerical_johns[counter-1]
203
204
205   
206    #print "Now the ANALYTIC"
207    pos_shore, vel_shore = shore(domain.time*sqrt(g*h_0)/L) #dimensionless
208    pos_shore = pos_shore*L
209    vel_shore = vel_shore*sqrt(g*h_0)
210
211    shorelines_analytical[counter-1] = pos_shore
212
213    #The following is for calculating the error at centroids
214    for i in range(n_C):
215        x = Cen[i]
216        if x < pos_shore:
217            sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless
218            AnalitW_C[i] = sta*h_0
219            AnalitU_C[i] = vel                            #It needs dimensionalisation
220            AnalitZ_C[i] = (h_0/L)*x - h_0
221            AnalitH_C[i] = AnalitW_C[i] - AnalitZ_C[i]
222            AnalitP_C[i] = AnalitH_C[i]*AnalitU_C[i]  #It needs dimensionalisation
223        else:
224            AnalitW_C[i] = (h_0/L)*x - h_0                 
225            AnalitU_C[i] = 0.0
226            AnalitZ_C[i] = (h_0/L)*x - h_0
227            AnalitH_C[i] = 0.0
228            AnalitP_C[i] = 0.0
229    AnalitU_C = AnalitU_C*sqrt(g*h_0) #This is the dimensionalisation
230    AnalitP_C = AnalitP_C*sqrt(g*h_0) #This is the dimensionalisation
231
232    error_W = (1.0/n_C)*sum(abs(StageC-AnalitW_C))
233    error_P = (1.0/n_C)*sum(abs(MomC-AnalitP_C))
234    error_U = (1.0/n_C)*sum(abs(VelC-AnalitU_C))
235    table = zeros((1,4),Float)
236    for r in range(1):
237        for c in range(4):
238            if c==0:
239                table[r][c] = domain.time
240            elif c==1:
241                table[r][c] = error_W
242            elif c==2:
243                table[r][c] = error_P
244            else:
245                table[r][c] = error_U
246    outname = "%s%04i%s%f%s" %("error_johns_", counter, "_", domain.time, ".csv")               
247    outfile = open(outname, 'w')
248    writer = csv.writer(outfile)
249    for row in table:
250        writer.writerow(row)       
251    outfile.close()
252
253
254    #The following is for ploting the quantities at vertex values
255    for i in range(n_V):
256        vector_x = Ver[i]
257        for k in range(2):
258            x = vector_x[k]
259            if x < pos_shore:
260                sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless
261                AnalitW_V[i,k] = sta*h_0
262                AnalitU_V[i,k] = vel                            #It needs dimensionalisation
263                AnalitZ_V[i,k] = (h_0/L)*x - h_0
264                AnalitH_V[i,k] = AnalitW_V[i,k] - AnalitZ_V[i,k]
265                AnalitP_V[i,k] = AnalitH_V[i,k]*AnalitU_V[i,k]  #It needs dimensionalisation
266            else:
267                AnalitW_V[i,k] = (h_0/L)*x - h_0                 
268                AnalitU_V[i,k] = 0.0
269                AnalitZ_V[i,k] = (h_0/L)*x - h_0
270                AnalitH_V[i,k] = 0.0
271                AnalitP_V[i,k] = 0.0
272    AnalitU_V = AnalitU_V*sqrt(g*h_0) #This is the dimensionalisation
273    AnalitP_V = AnalitP_V*sqrt(g*h_0) #This is the dimensionalisation
274
275
276    table = zeros((len(Ver.flat),6),Float)
277    for r in range(len(Ver.flat)):
278        for c in range(6):
279            if c==0:
280                table[r][c] = Ver.flat[r]
281            elif c==1:
282                table[r][c] = AnalitW_V.flat[r]
283            elif c==2:
284                table[r][c] = AnalitP_V.flat[r]
285            elif c==3:
286                table[r][c] = AnalitZ_V.flat[r]
287            elif c==4:
288                table[r][c] = AnalitH_V.flat[r]               
289            else:
290                table[r][c] = AnalitU_V.flat[r]
291
292    outname = "%s%04i%s%f%s" %("analytical_", counter, "_", domain.time, ".csv")               
293    outfile = open(outname, 'w')
294    writer = csv.writer(outfile)
295    for row in table:
296        writer.writerow(row)       
297    outfile.close()
298
299    #put this
300    from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
301    hold(False)
302    clf()
303
304    plot1 = subplot(311)
305    plot(Ver/1e+4,AnalitW_V,'b-',  Ver/1e+4,StageV,'g-',  Ver/1e+4,ElevationV,'k-')
306    #plot(Ver,StageV, Ver,ElevationV)
307    #plot(Ver/L,StageV/h_0,  Ver/L,ElevationV/h_0)
308    #xlabel('Position')
309    ylabel('Stage')
310    #plot1.set_xlim([0.0,1.2])
311    plot1.set_ylim([-6.0,6.0])#([-9.0e-3,9.0e-3])
312    #legend(('Analytical Solution', 'Numerical Solution', 'Discretized Bed'),
313    #       'upper right', shadow=False)
314   
315    plot2 = subplot(312)
316    plot(Ver/1e+4,AnalitP_V,'b-',  Ver/1e+4,MomV,'g-')
317    #plot(Ver/L, VelV/sqrt(g*h_0))
318    #xlabel('Position')
319    ylabel('Momentum')
320    #plot2.set_xlim([0.0,1.2])
321    #legend(('Analytical Solution','Numerical Solution'),
322    #       'upper right', shadow=False)       
323
324    plot3 = subplot(313)
325    plot(Ver/1e+4,AnalitU_V,'b-',  Ver/1e+4,VelV,'g-')
326    #plot(Ver/L, VelV/sqrt(g*h_0))
327    xlabel('Position / 10,000')
328    ylabel('Velocity')
329    #plot2.set_xlim([0.0,1.2])
330    legend(('Analytical Solution','Numerical Solution'),
331           'upper center', shadow=False)
332   
333    filename = "%s%04i%s%f%s" %("numerical_johns_", counter,"_", domain.time, ".png")
334    savefig(filename)
335    #show()
336    #raw_input("Press ENTER to continue")
337    #put this
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.