source: trunk/anuga_validation/convergence_study/plot_timeseries.py @ 7877

Last change on this file since 7877 was 5306, checked in by steve, 16 years ago

Working version of wave.py for Ole to work on.

File size: 3.9 KB
Line 
1## def animatesww2d_alt(sww_filename, movie_filename, range):
2##     """Plot cross section of model output
3##     """
4
5##     # Read SWW file 
6##     from Scientific.IO.NetCDF import NetCDFFile
7##     fid = NetCDFFile(sww_filename, 'r')
8   
9##     x = fid.variables['x']
10##     y = fid.variables['y']
11##     volumes = fid.variables['volumes']
12   
13##     elevation = fid.variables['elevation']
14##     time = fid.variables['time']
15##     stage = fid.variables['stage']
16##     xmomentum = fid.variables['xmomentum']
17##     ymomentum = fid.variables['ymomentum']
18   
19def animatesww2d(swwfile):
20    """Read in sww file and plot cross section of model output
21    """
22
23    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
24
25    try:
26        fid = open(swwfile)
27    except Exception, e:
28        msg = 'File "%s" could not be opened: Error="%s"'\
29              %(swwfile, e)
30        raise msg
31
32    print 'swwfile', swwfile
33
34    # interpolation points are for y = 0
35    # number of increments in x
36    m = 1000 
37    max_x = 100000.
38    x = 0
39    points = []
40    x_points = []
41    for i in range(m):
42        x += max_x/m
43        points.append([x,0.])
44        x_points.append(x)
45           
46    time_thinning = 1
47    use_cache = True
48    verbose = True
49
50    from anuga.abstract_2d_finite_volumes.util import file_function
51   
52    f = file_function(swwfile,
53                      quantities = sww_quantity,
54                      interpolation_points = points,
55                      time_thinning = time_thinning,
56                      verbose = verbose,
57                      use_cache = use_cache)
58
59    n = len(f.get_time()) # number of time steps
60
61    from Numeric import zeros, Float
62    from math import sqrt
63    stage = zeros((n,m), Float)
64    elevation = zeros((n,m), Float)
65    xmom = zeros((n,m), Float)
66    ymom = zeros((n,m), Float)
67    momenta = zeros((n,m), Float)
68    speed = zeros((n,m), Float)
69    max_stages = zeros((n,m), Float)
70    max_stage = None
71    min_stages = zeros((n,m), Float)
72    min_stage = 100.
73    fid = open('stage.csv', 'w')
74    for k, location in enumerate(x_points):
75        for i, t in enumerate(f.get_time()):
76            w = f(t, point_id = k)[0]
77            z = f(t, point_id = k)[1]
78            uh = f(t, point_id = k)[2]
79            vh = f(t, point_id = k)[3]
80            depth = w-z     
81            m = sqrt(uh*uh + vh*vh)
82            if depth < 0.001:
83                vel = 0.0
84            else:
85                vel = m / (depth + 1.e-6/depth)
86
87            stage[i,k] = w
88            s = '%.2f,' %(w)
89            fid.write(s)
90            elevation[i,k] = z
91            xmom[i,k] = uh
92            ymom[i,k] = vh
93            momenta[i,k] = m
94            speed[i,k] = vel
95
96            if w > max_stage: max_stage = w
97            if w < min_stage: min_stage = w
98            max_stages[i,k] = max_stage
99            min_stages[i,k] = min_stage
100        s = '\n'
101        fid.write(s)
102
103    min_stage = min(min(min_stages))
104    max_stage = max(max(max_stages))
105    stage_axis = [0, max_x, min_stage-1, max_stage+1]
106    from pylab import plot, xlabel, ylabel, savefig, close, hold, axis, title,show
107    hold(False)
108    figs = []
109    j = 1
110    #max_vec = zeros(100, Float)
111
112    #if sys.platform == 'win32':
113    try:
114        for i in range(0,n,20):
115            #max_vec[i] = max(stage[i,:])
116            plot(x_points,stage[i,:]) #x_points,max_vec,'+')
117            xlabel('x (m)')
118            ylabel('stage (m)')
119            axis(stage_axis)
120            name = 'time_%i' %i
121            savefig(name)
122            title(name)
123            figs.append(name)
124            show()
125    except:
126        pass
127   
128
129    close('all')
130
131if __name__ == '__main__':
132    import sys
133    swwfile = sys.argv[1]
134   
135    animatesww2d(swwfile)
136   
137#swwfile = join(getenv('INUNDATIONHOME'),'data','validation','convergence_study','5000','myexample2.sww')
138#animatesww2d(swwfile)
139
140# once png files have been created, they can then be dragged into
141# Windows Movie Maker and exported to a mpeg.
Note: See TracBrowser for help on using the repository browser.