source: trunk/anuga_work/development/2010-projects/anuga_1d/channel/var_depth_only.py @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 12 years ago

Moving 2010 project

File size: 3.7 KB
Line 
1import os
2import random
3from math import sqrt, pow, pi
4from channel_domain_Ab import *
5from Numeric import allclose, array, zeros, ones, Float, take, sqrt
6from config import g, epsilon
7
8
9print "Variable Width Only Test"
10
11# Define functions for initial quantities
12def initial_area(x):
13    y = zeros(len(x),Float)
14    for i in range(len(x)):
15        y[i]=(10-randomarray[i])
16
17   
18    return y
19
20def bed(x):
21
22    y = zeros(len(x),Float)
23    for i in range(len(x)):
24        y[i]=randomarray[i]
25        randomarray[i]=random.normalvariate(3,1)
26    return y
27
28def stage(x):
29
30    y = zeros(len(x),Float)
31    for i in range(len(x)):
32        y[i]=3+x[i]/2000
33    return y
34 
35
36
37def width(x):
38    return 1
39 
40import time
41
42# Set final time and yield time for simulation
43finaltime = 1
44yieldstep = 0.1
45
46# Length of channel (m)
47L = 1000.0   
48# Define the number of cells
49number_of_cells = [50]
50
51# Define cells for finite volume and their size
52N = int(number_of_cells[0])     
53print "Evaluating domain with %d cells" %N
54cell_len = L/N # Origin = 0.0
55points = zeros(N+1,Float)
56
57# Define the centroid points
58for j in range(N+1):
59    points[j] = j*cell_len
60
61# Create domain with centroid points as defined above       
62domain = Domain(points)
63
64# Define random array for width
65
66randomarray=zeros(len(points),Float)
67for j in range(N+1):
68    randomarray[j]=random.normalvariate(3,1)
69
70
71# Set initial values of quantities - default to zero
72domain.set_quantity('stage',stage)
73#domain.set_quantity('area', initial_area)
74domain.set_quantity('elevation',bed)
75domain.set_quantity('width',width)
76domain.setstageflag = True
77# Set boundry type, order, timestepping method and limiter
78domain.set_boundary({'exterior':Reflective_boundary(domain)})
79domain.order = 2
80domain.set_timestepping_method('rk2')
81domain.set_CFL(1.0)
82domain.set_limiter("vanleer")
83#domain.h0=0.0001
84
85# Start timer
86t0 = time.time()
87i=1
88
89print domain.quantities['elevation'].vertex_values
90
91for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
92    domain.write_time()
93
94    N = float(N)
95    HeightC = domain.quantities['height'].centroid_values
96    DischargeC = domain.quantities['discharge'].centroid_values
97    C = domain.centroids
98    print 'That took %.2f seconds' %(time.time()-t0)
99    X = domain.vertices
100    HeightQ = domain.quantities['height'].vertex_values
101    VelocityQ = domain.quantities['velocity'].vertex_values
102    x = X.flat
103    z = domain.quantities['elevation'].vertex_values.flat
104    stage=HeightQ.flat+z
105
106
107
108
109    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,savefig
110    import pylab as p
111#hold(False)
112   
113    plot1=subplot(211)
114
115    plot(x,stage,x,z)
116 
117    plot1.set_ylim([-1,11])
118    xlabel('Position')
119    ylabel('Stage')
120    legend(('Solution', 'Bed Height'),
121                'upper right', shadow=True)
122
123    plot2=subplot(212)
124    plot(x,VelocityQ.flat)
125    plot2.set_ylim([-5,5])
126    xlabel('Position') 
127    ylabel('Velocity')
128
129    savefig(str(i)+'.png')
130    i+=1
131    plot1.clear()
132
133print domain.quantities['elevation'].vertex_values
134
135N = float(N)
136HeightC = domain.quantities['height'].centroid_values
137DischargeC = domain.quantities['discharge'].centroid_values
138C = domain.centroids
139print 'That took %.2f seconds' %(time.time()-t0)
140X = domain.vertices
141HeightQ = domain.quantities['height'].vertex_values
142VelocityQ = domain.quantities['velocity'].vertex_values
143x = X.flat
144z = domain.quantities['elevation'].vertex_values.flat
145stage=HeightQ.flat+z
146
147
148
149
150from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
151
152#hold(False)
153   
154plot1=subplot(211)
155
156plot(x,stage,x,z)
157 
158plot1.set_ylim([-1,11])
159xlabel('Position')
160ylabel('Stage')
161legend(('Solution', 'Bed Height'),
162           'upper right', shadow=True)
163
164plot2=subplot(212)
165plot(x,VelocityQ.flat)
166plot2.set_ylim([-5,5])
167xlabel('Position')
168ylabel('Velocity')
169   
170show()
171
Note: See TracBrowser for help on using the repository browser.