Last change on this file since 7818 was 6453, checked in by steve, 14 years ago

File size: 4.1 KB
Line
1import os
2from math import sqrt, pow, pi
3from channel_domain_Ab import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6
7
9
10# Define functions for initial quantities
11
12def initial_area(x):
13    y=zeros(len(x),Float)
14    for i in range (len(x)):
15         if x[i]<=40:
16            y[i]=10*(x[i]*pi*2)
17         else:
18            y[i]=1*(x[i]*pi*2)
19    return y
20
21
22def width(x):
23    return 2*pi*x
24
25import time
26
27# Set final time and yield time for simulation
28finaltime = 2.0
29yieldstep = finaltime
30
31# Length of channel (m)
32L = 100.0
33# Define the number of cells
34number_of_cells = [100]
35
36# Define cells for finite volume and their size
37N = int(number_of_cells[0])
38print "Evaluating domain with %d cells" %N
39cell_len = L/N # Origin = 0.0
40points = zeros(N+1,Float)
41
42# Define the centroid points
43for j in range(N+1):
44    points[j] = j*cell_len
45
46# Create domain with centroid points as defined above
47domain = Domain(points)
48
49# Set initial values of quantities - default to zero
50domain.set_quantity('area', initial_area)
51domain.set_quantity('width',width)
52
53
54# Set boundry type, order, timestepping method and limiter
55domain.set_boundary({'exterior':Reflective_boundary(domain)})
56domain.order = 2
57domain.set_timestepping_method('rk2')
58domain.set_CFL(1.0)
59domain.set_limiter("vanleer")
60#domain.h0=0.0001
61
62# Start timer
63t0 = time.time()
64
65for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
66    domain.write_time()
67
68N = float(N)
69HeightC = domain.quantities['height'].centroid_values
70DischargeC = domain.quantities['discharge'].centroid_values
71C = domain.centroids
72print 'That took %.2f seconds' %(time.time()-t0)
73X = domain.vertices
74HeightQ = domain.quantities['height'].vertex_values
75VelocityQ = domain.quantities['velocity'].vertex_values
76x = X.flat
77z = domain.quantities['elevation'].vertex_values.flat
78stage=HeightQ.flat+z
79
80## ################# REPEAT ABOVE WITH MORE CELLS
81## # Set final time and yield time for simulation
82## finaltime2 = 2.0
83## yieldstep2 = finaltime2
84
85## # Length of channel (m)
86## L2 = 100.0
87## # Define the number of cells
88## number_of_cells2 = [100]
89
90## # Define cells for finite volume and their size
91## N2 = int(number_of_cells2[0])
92## print "Evaluating domain with %d cells" %N2
93## cell_len2 = L2/N2 # Origin = 0.0
94## points2 = zeros(N2+1,Float)
95
96## # Define the centroid points
97## for j in range(N2+1):
98##     points2[j] = j*cell_len2
99
100## # Create domain with centroid points as defined above
101## domain2 = Domain(points2)
102
103## # Set initial values of quantities - default to zero
104## domain2.set_quantity('area', initial_area)
105## domain2.set_quantity('width',width)
106
107
108## # Set boundry type, order, timestepping method and limiter
109## domain2.set_boundary({'exterior':Reflective_boundary(domain2)})
110## domain2.order = 2
111## domain2.set_timestepping_method('rk2')
112## domain2.set_CFL(1.0)
113## domain2.set_limiter("vanleer")
114## #domain.h0=0.0001
115
116## # Start timer
117## t02 = time.time()
118
119## for t in domain2.evolve(yieldstep = yieldstep2, finaltime = finaltime2):
120##     domain2.write_time()
121
122## N2 = float(N2)
123## HeightC2 = domain2.quantities['height'].centroid_values
124## DischargeC2 = domain2.quantities['discharge'].centroid_values
125## C2 = domain2.centroids
126## print 'That took %.2f seconds' %(time.time()-t02)
127## X2 = domain2.vertices
128## HeightQ2 = domain2.quantities['height'].vertex_values
129## VelocityQ2 = domain2.quantities['velocity'].vertex_values
130## x2 = X2.flat
131## z2 = domain2.quantities['elevation'].vertex_values.flat
132## stage2=HeightQ2.flat+z2
133
134## ################# REPEAT ABOVE
135
136from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
137import pickle
138f=open('highresdam.txt','r')
139## Data=[[x,stage],[x,VelocityQ.flat]]
140## pickle.dump(Data,f)
141## f.close()
143
144
145hold(False)
146plot1 = subplot(211)
147
148plot(x2,stage2,x,stage,'.')
149
150plot1.set_ylim([-1,11])
151xlabel('Position')
152ylabel('Stage')
153legend(('Analytical Solution', 'Numerical Solution'),