# source:anuga_work/development/anuga_1d/var_depth_only.py@7818

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

Fixed pylab so that it would work on bogong, by flattening the arrays for plotting

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'),
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'),