source: anuga_validation/convergence_study/animatesww2d_alt.py @ 7444

Last change on this file since 7444 was 5162, checked in by steve, 17 years ago

Updated some methods for quantity. Looks like we can use old
limiting system with larger values of beta.

File size: 3.9 KB
RevLine 
[5162]1## def animatesww2d_alt(sww_filename, movie_filename, range):
2##     """Plot cross section of model output
3##     """
[4695]4
[5162]5##     # Read SWW file 
6##     from Scientific.IO.NetCDF import NetCDFFile
7##     fid = NetCDFFile(sww_filename, 'r')
[4695]8   
[5162]9##     x = fid.variables['x']
10##     y = fid.variables['y']
11##     volumes = fid.variables['volumes']
[4695]12   
[5162]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']
[4695]18   
[4700]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
[4703]35    # number of increments in x
[5162]36    m = 1000 
[4700]37    max_x = 100000.
38    x = 0
39    points = []
40    x_points = []
41    for i in range(m):
42        x += max_x/m
[4703]43        points.append([x,0.])
[4700]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)
[4703]69    max_stages = zeros((n,m), Float)
[4700]70    max_stage = None
[4703]71    min_stages = zeros((n,m), Float)
72    min_stage = 100.
[4709]73    fid = open('stage.csv', 'w')
[4700]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
[4709]88            s = '%.2f,' %(w)
89            fid.write(s)
[4700]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
[4703]97            if w < min_stage: min_stage = w
98            max_stages[i,k] = max_stage
99            min_stages[i,k] = min_stage
[4709]100        s = '\n'
101        fid.write(s)
[4700]102
[4703]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]
[4838]106    from pylab import plot, xlabel, ylabel, savefig, close, hold, axis, title,show
[4700]107    hold(False)
[4703]108    figs = []
109    j = 1
110    #max_vec = zeros(100, Float)
111
[4838]112    #if sys.platform == 'win32':
113    try:
[4709]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)
[4838]124            show()
125    except:
126        pass
127   
[4700]128
129    close('all')
[4703]130
[4715]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)
[4709]139
[4703]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.