1 | import os |
---|
2 | from math import sqrt |
---|
3 | from shallow_water_h import * |
---|
4 | from Numeric import zeros, Float |
---|
5 | from analytic_dam_sudi import AnalyticDam |
---|
6 | |
---|
7 | h0=5.0 |
---|
8 | h1=10.0 |
---|
9 | |
---|
10 | analytical_sol=AnalyticDam(h0,h1) |
---|
11 | |
---|
12 | """ |
---|
13 | def newLinePlot(title='Simple Plot'): |
---|
14 | import Gnuplot |
---|
15 | gg=Gnuplot.Gnuplot(persist=0) |
---|
16 | gg.title(title) |
---|
17 | gg('set data style linespoints') |
---|
18 | gg.xlabel('x') |
---|
19 | gg.ylabel('y') |
---|
20 | return gg |
---|
21 | |
---|
22 | def linePlot(gg, x1, y1, x2, y2): |
---|
23 | import Gnuplot |
---|
24 | plot1=Gnuplot.PlotItems.Data(x1.flat, y1.flat, with="linespoints") |
---|
25 | plot2=Gnuplot.PlotItems.Data(x2.flat, y2.flat, with="lines 3") |
---|
26 | gg.plot(plot1, plot2) |
---|
27 | """ |
---|
28 | |
---|
29 | |
---|
30 | print "TEST 1D-SOLUTION I" |
---|
31 | |
---|
32 | L=2000.0 |
---|
33 | N=100 |
---|
34 | |
---|
35 | cell_len=L/N |
---|
36 | |
---|
37 | points=zeros(N+1, Float) |
---|
38 | for i in range(N+1): |
---|
39 | points[i]=i*cell_len |
---|
40 | |
---|
41 | domain=Domain(points) |
---|
42 | |
---|
43 | domain.default_order = 1 |
---|
44 | domain.default_time_order = 1 |
---|
45 | #domain.cfl = 1.0 |
---|
46 | #domain.limiter = "minmod" |
---|
47 | |
---|
48 | |
---|
49 | |
---|
50 | |
---|
51 | |
---|
52 | |
---|
53 | |
---|
54 | def height(x): |
---|
55 | y=zeros(len(x), Float) |
---|
56 | for i in range (len(x)): |
---|
57 | if x[i]<=L/4.0: |
---|
58 | y[i]=0.0 #h0 |
---|
59 | elif x[i]<=3*L/4.0: |
---|
60 | y[i]=h1 |
---|
61 | else: |
---|
62 | y[i]=h0 |
---|
63 | return y |
---|
64 | |
---|
65 | domain.set_quantity('height', height) |
---|
66 | domain.order=domain.default_order |
---|
67 | print "domain order", domain.order |
---|
68 | |
---|
69 | domain.set_boundary({'exterior':Reflective_boundary(domain)}) |
---|
70 | |
---|
71 | X=domain.vertices |
---|
72 | C=domain.centroids |
---|
73 | #plot1x=newLinePlot("Height") |
---|
74 | #plot2x=newLinePlot("Momentum") |
---|
75 | |
---|
76 | |
---|
77 | import time |
---|
78 | t0=time.time() |
---|
79 | yieldstep=30.0 |
---|
80 | finaltime=30.0 |
---|
81 | print "integral", domain.quantities['height'].get_integral() |
---|
82 | for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime): |
---|
83 | domain.write_time() |
---|
84 | print "integral", domain.quantities['height'].get_integral() |
---|
85 | if t>0.0: |
---|
86 | HeightQ=domain.quantities['height'].vertex_values |
---|
87 | MomentumQ=domain.quantities['xmomentum'].vertex_values |
---|
88 | h, uh=analytical_sol(X.flat, domain.time) |
---|
89 | #linePlot(plot1x, X, HeightQ, X, h) |
---|
90 | #linePlot(plot2x, X, MomentumQ, X, uh) |
---|
91 | #print "press return" |
---|
92 | #pass |
---|
93 | |
---|
94 | from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot |
---|
95 | #print 'Test1' |
---|
96 | hold(False) |
---|
97 | #print 'test 2' |
---|
98 | plot1 = subplot(211) |
---|
99 | #print 'test 3' |
---|
100 | |
---|
101 | plot(X,h,X,HeightQ) |
---|
102 | #print 'Test4' |
---|
103 | plot1.set_ylim([0,11]) |
---|
104 | xlabel('Position') |
---|
105 | ylabel('Stage') |
---|
106 | #legend(('Analytical Solution', 'Numerical Solution'), |
---|
107 | # 'lower right', shadow=False) |
---|
108 | plot2 = subplot(212) |
---|
109 | plot(X,uh,X,MomentumQ) |
---|
110 | #plot2.set_ylim([-5,35]) |
---|
111 | legend(('Analytical Solution', 'Numerical Solution'), |
---|
112 | 'lower right', shadow=False) |
---|
113 | |
---|
114 | xlabel('Position') |
---|
115 | ylabel('Xmomentum') |
---|
116 | |
---|
117 | file = "dam_h_" |
---|
118 | #file += str(number_of_cells[i]) |
---|
119 | file += ".eps" |
---|
120 | #savefig(file) |
---|
121 | show() |
---|
122 | |
---|
123 | print 'That took %.2f seconds'%(time.time()-t0) |
---|