source: anuga_work/development/anuga_1d/stead_state_varying_width.py @ 6453

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

Added Padarn's (2008/2009 summer student) files

File size: 3.8 KB
Line 
1import os
2from math import sqrt, pow, pi
3from channel_domain import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6
7
8print "Radial Dam Break Test"
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 = [200]
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
82finaltime2 = 2.0
83yieldstep2 = finaltime2
84
85# Length of channel (m)
86L2 = 100.0   
87# Define the number of cells
88number_of_cells2 = [2000]
89
90# Define cells for finite volume and their size
91N2 = int(number_of_cells2[0])     
92print "Evaluating domain with %d cells" %N2
93cell_len2 = L2/N2 # Origin = 0.0
94points2 = zeros(N2+1,Float)
95
96# Define the centroid points
97for j in range(N2+1):
98    points2[j] = j*cell_len2
99
100# Create domain with centroid points as defined above       
101domain2 = Domain(points2)
102
103# Set initial values of quantities - default to zero
104domain2.set_quantity('area', initial_area)
105domain2.set_quantity('width',width)
106
107
108# Set boundry type, order, timestepping method and limiter
109domain2.set_boundary({'exterior':Reflective_boundary(domain2)})
110domain2.order = 2
111domain2.set_timestepping_method('rk2')
112domain2.set_CFL(1.0)
113domain2.set_limiter("vanleer")
114#domain.h0=0.0001
115
116# Start timer
117t02 = time.time()
118
119for t in domain2.evolve(yieldstep = yieldstep2, finaltime = finaltime2):
120    domain2.write_time()
121
122N2 = float(N2)
123HeightC2 = domain2.quantities['height'].centroid_values
124DischargeC2 = domain2.quantities['discharge'].centroid_values
125C2 = domain2.centroids
126print 'That took %.2f seconds' %(time.time()-t02)
127X2 = domain2.vertices
128HeightQ2 = domain2.quantities['height'].vertex_values
129VelocityQ2 = domain2.quantities['velocity'].vertex_values
130x2 = X2.flat
131z2 = domain2.quantities['elevation'].vertex_values.flat
132stage2=HeightQ2.flat+z2
133
134################# REPEAT ABOVE
135
136from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
137
138hold(False)
139plot1 = subplot(211)
140
141plot(x,stage,'ro',x2,stage2)
142 
143plot1.set_ylim([-1,11])
144xlabel('Position')
145ylabel('Stage')
146legend(('Analytical Solution', 'Numerical Solution'),
147           'upper right', shadow=True)
148plot2 = subplot(212)
149plot(x,VelocityQ.flat)
150plot2.set_ylim([-10,10])
151   
152xlabel('Position')
153ylabel('Velocity')
154   
155show()
156
157
158
Note: See TracBrowser for help on using the repository browser.