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