source: anuga_work/development/sudi/sw_1d/periodic_waves/johns/run-plot.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: 7.3 KB
Line 
1#################Opening analytical solution################
2infile = open('analytical_1230_12300.000000.csv', 'r')
3table = []
4for row in csv.reader(infile):
5    table.append(row)
6infile.close()
7for r in range(len(table)):
8    for c in range(len(table[0])):
9        table[r][c] = float(table[r][c])
10
11position = zeros(len(table),Float)
12analitW12300 = zeros(len(table),Float)
13analitP12300 = zeros(len(table),Float)
14analitZ12300 = zeros(len(table),Float)
15analitH12300 = zeros(len(table),Float)
16analitU12300 = zeros(len(table),Float)
17for r in range(len(table)):
18    for c in range(len(table[0])):
19        if c==0:
20            position[r] = table[r][c]
21        elif c==1:
22            analitW12300[r] = table[r][c]
23        elif c==2:
24            analitP12300[r] = table[r][c]
25        elif c==3:
26            analitZ12300[r] = table[r][c]
27        elif c==4:
28            analitH12300[r] = table[r][c]
29        elif c==5:
30            analitU12300[r] = table[r][c]           
31
32#################Opening numerical_cg solution################
33infile = open('numerical_cg_1230_12300.000000.csv', 'r')
34table = []
35for row in csv.reader(infile):
36    table.append(row)
37infile.close()
38for r in range(len(table)):
39    for c in range(len(table[0])):
40        table[r][c] = float(table[r][c])
41
42numeric_cgW12300 = zeros(len(table),Float)
43numeric_cgP12300 = zeros(len(table),Float)
44numeric_cgZ12300 = zeros(len(table),Float)
45numeric_cgH12300 = zeros(len(table),Float)
46numeric_cgU12300 = zeros(len(table),Float)
47
48for r in range(len(table)):
49    for c in range(len(table[0])):
50        if c==1:
51            numeric_cgW12300[r] = table[r][c]
52        elif c==2:
53            numeric_cgP12300[r] = table[r][c]
54        elif c==3:
55            numeric_cgZ12300[r] = table[r][c]
56        elif c==4:
57            numeric_cgH12300[r] = table[r][c]
58        elif c==5:
59            numeric_cgU12300[r] = table[r][c]
60           
61
62
63#################Opening numerical_johns solution################
64infile = open('numerical_johns_1230_12300.000000.csv', 'r')
65table = []
66for row in csv.reader(infile):
67    table.append(row)
68infile.close()
69for r in range(len(table)):
70    for c in range(len(table[0])):
71        table[r][c] = float(table[r][c])
72
73numeric_johnsW12300 = zeros(len(table),Float)
74numeric_johnsP12300 = zeros(len(table),Float)
75numeric_johnsZ12300 = zeros(len(table),Float)
76numeric_johnsH12300 = zeros(len(table),Float)
77numeric_johnsU12300 = zeros(len(table),Float)
78
79for r in range(len(table)):
80    for c in range(len(table[0])):
81        if c==1:
82            numeric_johnsW12300[r] = table[r][c]
83        elif c==2:
84            numeric_johnsP12300[r] = table[r][c]
85        elif c==3:
86            numeric_johnsZ12300[r] = table[r][c]
87        elif c==4:
88            numeric_johnsH12300[r] = table[r][c]
89        elif c==5:
90            numeric_johnsU12300[r] = table[r][c]
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111#################Opening analytical solution################
112infile = open('analytical_1260_12600.000000.csv', 'r')
113table = []
114for row in csv.reader(infile):
115    table.append(row)
116infile.close()
117for r in range(len(table)):
118    for c in range(len(table[0])):
119        table[r][c] = float(table[r][c])
120
121position = zeros(len(table),Float)
122analitW12600 = zeros(len(table),Float)
123analitP12600 = zeros(len(table),Float)
124analitZ12600 = zeros(len(table),Float)
125analitH12600 = zeros(len(table),Float)
126analitU12600 = zeros(len(table),Float)
127for r in range(len(table)):
128    for c in range(len(table[0])):
129        if c==0:
130            position[r] = table[r][c]
131        elif c==1:
132            analitW12600[r] = table[r][c]
133        elif c==2:
134            analitP12600[r] = table[r][c]
135        elif c==3:
136            analitZ12600[r] = table[r][c]
137        elif c==4:
138            analitH12600[r] = table[r][c]
139        elif c==5:
140            analitU12600[r] = table[r][c]           
141
142#################Opening numerical_cg solution################
143infile = open('numerical_cg_1260_12600.000000.csv', 'r')
144table = []
145for row in csv.reader(infile):
146    table.append(row)
147infile.close()
148for r in range(len(table)):
149    for c in range(len(table[0])):
150        table[r][c] = float(table[r][c])
151
152numeric_cgW12600 = zeros(len(table),Float)
153numeric_cgP12600 = zeros(len(table),Float)
154numeric_cgZ12600 = zeros(len(table),Float)
155numeric_cgH12600 = zeros(len(table),Float)
156numeric_cgU12600 = zeros(len(table),Float)
157
158for r in range(len(table)):
159    for c in range(len(table[0])):
160        if c==1:
161            numeric_cgW12600[r] = table[r][c]
162        elif c==2:
163            numeric_cgP12600[r] = table[r][c]
164        elif c==3:
165            numeric_cgZ12600[r] = table[r][c]
166        elif c==4:
167            numeric_cgH12600[r] = table[r][c]
168        elif c==5:
169            numeric_cgU12600[r] = table[r][c]
170           
171
172
173#################Opening numerical_johns solution################
174infile = open('numerical_johns_1260_12600.000000.csv', 'r')
175table = []
176for row in csv.reader(infile):
177    table.append(row)
178infile.close()
179for r in range(len(table)):
180    for c in range(len(table[0])):
181        table[r][c] = float(table[r][c])
182
183numeric_johnsW12600 = zeros(len(table),Float)
184numeric_johnsP12600 = zeros(len(table),Float)
185numeric_johnsZ12600 = zeros(len(table),Float)
186numeric_johnsH12600 = zeros(len(table),Float)
187numeric_johnsU12600 = zeros(len(table),Float)
188
189for r in range(len(table)):
190    for c in range(len(table[0])):
191        if c==1:
192            numeric_johnsW12600[r] = table[r][c]
193        elif c==2:
194            numeric_johnsP12600[r] = table[r][c]
195        elif c==3:
196            numeric_johnsZ12600[r] = table[r][c]
197        elif c==4:
198            numeric_johnsH12600[r] = table[r][c]
199        elif c==5:
200            numeric_johnsU12600[r] = table[r][c]
201
202
203
204#################PLOTTING THE RESULTS##############################################
205from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
206hold(False)
207clf()
208
209plot1 = subplot(311)
210plot(#position*1e-4,analitW12300, position*1e-4,numeric_cgW12300, position*1e-4,numeric_johnsW12300,
211     position*1e-4,analitW12600, position*1e-4,numeric_cgW12600, position*1e-4,numeric_johnsW12600,
212     position*1e-4,analitZ12300,'k-')
213#xlabel('Position')
214ylabel('Stage')
215#plot1.set_xlim([0.0,1.2])
216#plot1.set_ylim([-5.0,3.0])#([-9.0e-3,9.0e-3])
217#legend(('Analytical Solution', 'C-G Solution', 'Johns Solution'),
218#       'lower left', shadow=False)
219
220plot2 = subplot(312)
221plot(#position*1e-4,analitP12300, position*1e-4,numeric_cgP12300, position*1e-4,numeric_johnsP12300)
222     position*1e-4,analitP12600, position*1e-4,numeric_cgP12600, position*1e-4,numeric_johnsP12600)
223#xlabel('Position')
224ylabel('Momentum')
225#plot1.set_xlim([0.0,1.2])
226#plot1.set_ylim([-4.0,3.0])#([-9.0e-3,9.0e-3])
227#legend(('Analytical Solution', 'C-G Solution', 'Johns Solution'),
228#       'lower left', shadow=False)
229
230
231plot3 = subplot(313)
232plot(#position*1e-4,analitU12300, position*1e-4,numeric_cgU12300, position*1e-4,numeric_johnsU12300)
233     position*1e-4,analitU12600, position*1e-4,numeric_cgU12600, position*1e-4,numeric_johnsU12600)
234xlabel('Position/10,000')
235ylabel('Velocity')
236#plot1.set_xlim([0.0,1.2])
237#plot3.set_ylim([-4.0,2.0])#([-9.0e-3,9.0e-3])
238legend(('Exact Solution', 'C-G Solution', 'Johns Solution'),
239       'lower left', shadow=False)
240
241filename = "%s" %("12600.eps")
242#savefig(filename)
243show()
244
245
Note: See TracBrowser for help on using the repository browser.