source: development/momentum_sink/scripts/plot_velocity_CS.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 3.1 KB
Line 
1"""Read in sww file, interpolate at specified locations and
2    Cross section at specified T
3
4"""
5from os import sep
6import Numeric
7import project
8from anuga.pyvolution.util import file_function
9from create_buildings import  WidtH
10# from run_friction import fric
11from pylab import *
12
13fric = 33   
14#swwfile = project.newoutputname + '.sww'
15#swwfile = project.outputname
16
17#swwfile_B = project.outputdir + sep  + 'Buildings_3790.sww'
18swwfile_B = project.outputdir + sep  + 'buildings_BR=17_1106.sww'
19#swwfile_B = project.outputdir + sep  + 'friction_n=10_3135.sww'
20swwfile_F = project.outputdir + sep  + 'friction_n=Dbl_12_25.sww'
21
22gauge_depth = Numeric.arrayrange(41.5, 5*WidtH, 25) # Random special
23#gauge_depth = Numeric.arrayrange(0, 1.5*WidtH, 25)
24gauge_breadth = 100
25gauge_locations = []
26
27for GD in gauge_depth:
28    gauge_location = [GD,gauge_breadth]
29    gauge_locations.append(gauge_location)
30
31
32#Read model output
33quantities = ['stage', 'xmomentum', 'ymomentum']
34
35f_B = file_function(swwfile_B,
36                  quantities = quantities,
37                  interpolation_points = gauge_locations,
38                  verbose = True,
39                  use_cache = True)
40
41f_F = file_function(swwfile_F,
42                  quantities = quantities,
43                  interpolation_points = gauge_locations,
44                  verbose = True,
45                  use_cache = True)
46
47
48T = Numeric.arrayrange(700,1000,10)
49#T = [ 20, 50, 100, 150, 200 ]
50
51for t in T:
52    model_time = []
53    stages_B = []
54    stages_F = []
55    elevations = []
56    momenta = []
57    velocityB = []
58    velocityF=[]
59    for i, GL in enumerate(gauge_depth):
60
61        # collects data from specified time values in sww file
62        wB = f_B(t, point_id = i)[0]
63        uhB = f_B(t, point_id = i)[1]
64        vhB = f_B(t, point_id = i)[2]
65
66        wF = f_F(t, point_id = i)[0]       
67        uhF = f_B(t, point_id = i)[1]
68        vhF = f_B(t, point_id = i)[2]
69        #m = sqrt(uh*uh + vh*vh)   #Absolute momentum
70        velB = sqrt(uhB*uhB + vhB*vhB) / (wB + 1.e-30) #Absolute velocity
71        velF = sqrt(uhF*uhF + vhF*vhF) / (wF + 1.e-30) #Absolute velocity                   
72        #print velB, "building velocity"
73        #print velF, "friction velocity"
74
75        stages_B.append(wB)
76        stages_F.append(wF)
77        #momenta.append(mB)
78        velocityB.append(velB)
79        velocityF.append(velF)
80       
81    #print len(stages), "stages"
82    diff = abs(Numeric.array(stages_B) - Numeric.array(stages_F))
83    #print diff, "difference"
84
85
86    ion()
87    hold(True)
88    subplot(211)
89   
90    plot(gauge_depth, velocityB, '-b',gauge_depth,velocityF,'-r')
91    title('Velocity_%d_ vs gauge length (red: fric, blue: buildings)' %t)
92    xlabel('gauge length(m)')
93    ylabel('water velocities (m/s)')
94    #savefig('Comp_Time%d_B_vs_F' %t)
95    hold(False)
96#savefig('Final_comp_n=%s' %fric)
97subplot(212)
98subplots_adjust(hspace = 0.4)
99plot(gauge_depth, diff, 'g')
100title('Final Differences in water levels')
101xlabel('gauge length(m)')
102ylabel('water differences (m)')
103savefig('Final_Water_velocities_n=%s' %fric)
104print "finished"
105
106
Note: See TracBrowser for help on using the repository browser.