1 | import os |
---|
2 | from math import sqrt |
---|
3 | #from shallow_water_h import * |
---|
4 | from shallow_water_domain_suggestion1 import * |
---|
5 | from Numeric import zeros, Float |
---|
6 | from analytic_dam_sudi import AnalyticDam |
---|
7 | from scipy import linspace |
---|
8 | |
---|
9 | h0=5.0 |
---|
10 | h1=10.0 |
---|
11 | |
---|
12 | analytical_sol=AnalyticDam(h0,h1) |
---|
13 | |
---|
14 | """ |
---|
15 | def newLinePlot(title='Simple Plot'): |
---|
16 | import Gnuplot |
---|
17 | gg=Gnuplot.Gnuplot(persist=0) |
---|
18 | gg.title(title) |
---|
19 | gg('set data style linespoints') |
---|
20 | gg.xlabel('x') |
---|
21 | gg.ylabel('y') |
---|
22 | return gg |
---|
23 | |
---|
24 | def linePlot(gg, x1, y1, x2, y2): |
---|
25 | import Gnuplot |
---|
26 | plot1=Gnuplot.PlotItems.Data(x1.flat, y1.flat, with="linespoints") |
---|
27 | plot2=Gnuplot.PlotItems.Data(x2.flat, y2.flat, with="lines 3") |
---|
28 | gg.plot(plot1, plot2) |
---|
29 | """ |
---|
30 | |
---|
31 | |
---|
32 | print "TEST 1D-SOLUTION I" |
---|
33 | |
---|
34 | L=2000.0 |
---|
35 | N=8000 |
---|
36 | |
---|
37 | cell_len=L/N |
---|
38 | |
---|
39 | points=zeros(N+1, Float) |
---|
40 | for i in range(N+1): |
---|
41 | points[i]=i*cell_len |
---|
42 | |
---|
43 | domain=Domain(points) |
---|
44 | |
---|
45 | domain.order = 2 |
---|
46 | domain.set_timestepping_method('rk2') |
---|
47 | domain.cfl = 1.0 |
---|
48 | domain.limiter = "minmod" |
---|
49 | |
---|
50 | |
---|
51 | |
---|
52 | def height(x): |
---|
53 | y=zeros(len(x), Float) |
---|
54 | for i in range (len(x)): |
---|
55 | if x[i]<=L/4.0: |
---|
56 | y[i]=0.0 #h0 |
---|
57 | elif x[i]<=3*L/4.0: |
---|
58 | y[i]=h1 |
---|
59 | else: |
---|
60 | y[i]=h0 |
---|
61 | return y |
---|
62 | |
---|
63 | def nol(x): |
---|
64 | y=zeros(len(x), Float) |
---|
65 | for i in range (len(x)): |
---|
66 | y[i]=0.0 |
---|
67 | return y |
---|
68 | |
---|
69 | domain.set_quantity('stage',height) |
---|
70 | print "domain order=", domain.order |
---|
71 | print "domain limiter=", domain.limiter |
---|
72 | |
---|
73 | domain.set_boundary({'exterior':Reflective_boundary(domain)}) |
---|
74 | |
---|
75 | X=domain.vertices |
---|
76 | C=domain.centroids |
---|
77 | #plot1x=newLinePlot("Height") |
---|
78 | #plot2x=newLinePlot("Momentum") |
---|
79 | |
---|
80 | P=linspace(0,2000,1000) |
---|
81 | H=height(P) |
---|
82 | Enol=nol(P) |
---|
83 | |
---|
84 | |
---|
85 | import time |
---|
86 | yieldstep=finaltime=0.01 |
---|
87 | t0=time.time() |
---|
88 | |
---|
89 | """ |
---|
90 | for i in range(10): |
---|
91 | L=2000.0 |
---|
92 | N=400 |
---|
93 | cell_len=L/N |
---|
94 | points=zeros(N+1, Float) |
---|
95 | for i in range(N+1): |
---|
96 | points[i]=i*cell_len |
---|
97 | domain=Domain(points) |
---|
98 | domain.order = 2 |
---|
99 | domain.set_timestepping_method('rk2') |
---|
100 | domain.cfl = 1.0 |
---|
101 | domain.limiter = "vanleer" |
---|
102 | domain.set_quantity('stage',height) |
---|
103 | domain.set_boundary({'exterior':Reflective_boundary(domain)}) |
---|
104 | #print "integral", domain.quantities['stage'].get_integral() |
---|
105 | for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime): |
---|
106 | domain.write_time() |
---|
107 | #print "integral", domain.quantities['stage'].get_integral() |
---|
108 | #print '===================================================================' |
---|
109 | print 'The average time is %.2f seconds'%((time.time()-t0)/10.0) |
---|
110 | """ |
---|
111 | while finaltime<0.011: |
---|
112 | i = 0 |
---|
113 | yieldstep=finaltime |
---|
114 | for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime): |
---|
115 | domain.write_time() |
---|
116 | #if t>=0.0: |
---|
117 | print "t=",t |
---|
118 | N = float(N) |
---|
119 | StageC = domain.quantities['stage'].centroid_values |
---|
120 | XmomC = domain.quantities['xmomentum'].centroid_values |
---|
121 | VelC = domain.quantities['velocity'].centroid_values |
---|
122 | C = domain.centroids |
---|
123 | |
---|
124 | """ |
---|
125 | hC, uhC, uC = analytical_sol(C,domain.time) |
---|
126 | h_error = cell_len*sum(abs(hC-StageC)) |
---|
127 | uh_error = cell_len*sum(abs(uhC-XmomC)) |
---|
128 | u_error = cell_len*sum(abs(uC-VelC)) |
---|
129 | print "h_error %.10f" %(h_error) |
---|
130 | print "uh_error %.10f"%(uh_error) |
---|
131 | print "u_error %.10f" %(u_error) |
---|
132 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
133 | """ |
---|
134 | |
---|
135 | X = domain.vertices |
---|
136 | StageQ = domain.quantities['stage'].vertex_values |
---|
137 | XmomQ = domain.quantities['xmomentum'].vertex_values |
---|
138 | VelQ = domain.quantities['velocity'].vertex_values |
---|
139 | |
---|
140 | #h, uh, u = analytical_sol(X.flat, domain.time) |
---|
141 | #linePlot(plot1x, X, HeightQ, X, h) |
---|
142 | #linePlot(plot2x, X, MomentumQ, X, uh) |
---|
143 | #print "press return" |
---|
144 | #pass |
---|
145 | |
---|
146 | from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot |
---|
147 | hold(False) |
---|
148 | |
---|
149 | plot1 = subplot(311) |
---|
150 | plot(P,H, X,StageQ) |
---|
151 | plot1.set_ylim([-1,11]) |
---|
152 | #plot1.set_xlim([-100.0,2000.0]) |
---|
153 | #xlabel('Position') |
---|
154 | ylabel('Stage') |
---|
155 | #legend(('Analytical Solution', 'Numerical Solution'), |
---|
156 | # 'lower right', shadow=False) |
---|
157 | |
---|
158 | plot2 = subplot(312) |
---|
159 | plot(P,Enol, X,XmomQ) |
---|
160 | #plot2.set_ylim([-5,35]) |
---|
161 | #plot2.set_xlim([-100.0,2000.0]) |
---|
162 | #legend(('Analytical Solution', 'Numerical Solution'), |
---|
163 | # 'lower right', shadow=False) |
---|
164 | #xlabel('Position') |
---|
165 | ylabel('Momentum') |
---|
166 | |
---|
167 | plot3 = subplot(313) |
---|
168 | plot(P,Enol, X,VelQ) |
---|
169 | #plot2.set_ylim([-5,35]) |
---|
170 | #plot3.set_xlim([-100.0,2000.0]) |
---|
171 | legend(('Penyelesaian Analitis', 'Penyelesaian Numeris'), |
---|
172 | 'lower right', shadow=False) |
---|
173 | xlabel('Position') |
---|
174 | ylabel('Kecepatan') |
---|
175 | |
---|
176 | #file = "dam_h_" |
---|
177 | #file += str(number_of_cells[i]) |
---|
178 | #file += ".eps" |
---|
179 | #savefig(file) |
---|
180 | #show() |
---|
181 | |
---|
182 | #file = "dam_" |
---|
183 | #file += str(finaltime) |
---|
184 | #file += ".png" |
---|
185 | #savefig(file) |
---|
186 | #show() |
---|
187 | filename = "%s%04i%s" %("dambreak", i, ".eps") |
---|
188 | savefig(filename) |
---|
189 | i = i + 1 |
---|
190 | print "The domain.limiter is", domain.limiter |
---|
191 | print 'That took %.2f seconds'%(time.time()-t0) |
---|
192 | print '=============================================================================' |
---|
193 | finaltime = finaltime + 0.25 |
---|
194 | |
---|