source: trunk/anuga_work/development/gareth/tests/parabolic/parabolaplot.py @ 9358

Last change on this file since 9358 was 9022, checked in by davies, 11 years ago

Removing most significant velocity artefacts from dev_audusse

File size: 2.6 KB
Line 
1"""View results of runup.py
2"""
3#---------------
4# Import Modules
5#---------------
6#import anuga
7#import struct
8import numpy
9#import scipy
10from matplotlib import pyplot as pyplot
11#import util
12from bal_and import plot_utils as util
13
14p=util.get_output('parabola_v2.sww', 0.001)
15p2=util.get_centroids(p, velocity_extrapolation=True)
16
17# Define some 'lines' along which to plot
18v=(p2.y==p2.y[0])
19
20#--------------------
21# Make plot animation
22#--------------------
23pyplot.close() #If the plot is open, there will be problems
24pyplot.ion()
25
26if False:
27    line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,(p2.stage[:,v].min(),p2.stage[:,v].max() ) )
28    for i in range(700,p2.xmom.shape[0]):
29        line.set_xdata(p2.x[v])
30        line.set_ydata(p2.stage[i,v])
31        pyplot.draw()
32        #pyplot.plot( (0,1),(0,0), 'r' )
33        #pyplot.plot(p2.x[v],p2.elev[v])
34        pyplot.title(str(i)+'/200') # : velocity does not converge to zero' )
35        pyplot.xlabel('x')
36        pyplot.ylabel('stage (m/s)')
37
38if False:
39    pyplot.clf()
40 
41    line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,(p2.xvel[:,v].min(),p2.xvel[:,v].max() ), 'ro' )
42    for i in range(700,p2.xmom.shape[0]):
43        line.set_xdata(p2.x[v])
44        line.set_ydata(p2.xvel[i,v])
45        pyplot.draw()
46        pyplot.plot( (p2.x[v].min(),p2.x[v].max()),(0,0), 'ro' )
47        #pyplot.plot(x[v],elev[v])
48        pyplot.title(str(i)+'/200') # : velocity does not converge to zero' )
49        pyplot.xlabel('x')
50        pyplot.ylabel('xvel (m/s)')
51
52D0=4.0
53L=10.0
54A=2.0
55g=9.8
56#t=time[30]
57
58omega=numpy.sqrt(2*D0*g)/L
59T= 2*numpy.pi/omega
60ppp= (abs(p2.x-4.*L/2.)).argmin()
61ppp2= (abs(p2.x-4.*L/4.)).argmin()
62
63w = D0 + 2*A*D0/(L**2)*numpy.cos(omega*p2.time)*( (p2.x[ppp]-4.*L/2.) -A/2.*numpy.cos(omega*p2.time))
64w2 = D0 + 2*A*D0/(L**2)*numpy.cos(omega*p2.time)*( (p2.x[ppp2]-4.*L/2.) -A/2.*numpy.cos(omega*p2.time))
65w2 = w2*(w2>p2.elev[ppp2])+p2.elev[ppp2]*(w2<=p2.elev[ppp2])
66
67pyplot.clf()
68pyplot.plot(p2.time,w, color='blue')
69pyplot.plot(p2.time,p2.stage[:,ppp], color='green')
70pyplot.savefig('Stage_centre_v2.png')
71#       pyplot.savefig('runup_x_velocities.png')
72pyplot.clf()
73pyplot.plot(p2.time,w2, color='blue')
74pyplot.plot(p2.time,p2.stage[:,ppp2], color='green')
75pyplot.savefig('Stage_centre_v3.png')
76
77
78pltind=numpy.argmin(abs(p2.time-T*3))
79
80pyplot.clf()
81pyplot.plot(p2.x[v],p2.xvel[pltind,v])
82pyplot.title('Velocity near to time=3T')
83pyplot.savefig('Vel_3T_v2.png')
84
85pltind=numpy.argmin(abs(p2.time-T*3.25))
86
87pyplot.clf()
88pyplot.plot(p2.x[v],p2.xvel[pltind,v])
89pyplot.title('Velocity near to time=3.25T')
90pyplot.savefig('Vel_3_5T_v2.png')
91#pyplot.clf()
92#pyplot.close()
Note: See TracBrowser for help on using the repository browser.