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

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

Moving 2010 project

File size: 2.7 KB
RevLine 
[7884]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 width(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]=8.0
33    return y
34 
35
36
37 
38import time
39
40# Set final time and yield time for simulation
41finaltime = 10.0
42yieldstep = finaltime
43
44# Length of channel (m)
45L = 1000.0   
46# Define the number of cells
47number_of_cells = [50]
48
49# Define cells for finite volume and their size
50N = int(number_of_cells[0])     
51print "Evaluating domain with %d cells" %N
52cell_len = L/N # Origin = 0.0
53points = zeros(N+1,Float)
54
55# Define the centroid points
56for j in range(N+1):
57    points[j] = j*cell_len
58
59# Create domain with centroid points as defined above       
60domain = Domain(points)
61
62# Define random array for width
63
64randomarray=zeros(len(points),Float)
65for j in range(N+1):
66    randomarray[j]=random.normalvariate(3,1)
67
68
69# Set initial values of quantities - default to zero
70domain.set_quantity('stage',stage)
71#domain.set_quantity('area', initial_area)
72domain.set_quantity('width',width)
73domain.setstageflag = True
74# Set boundry type, order, timestepping method and limiter
75domain.set_boundary({'exterior':Reflective_boundary(domain)})
76domain.order = 2
77domain.set_timestepping_method('rk2')
78domain.set_CFL(1.0)
79domain.set_limiter("vanleer")
80#domain.h0=0.0001
81
82# Start timer
83t0 = time.time()
84
85print domain.quantities['elevation'].vertex_values
86
87for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
88    domain.write_time()
89
90print domain.quantities['elevation'].vertex_values
91
92N = float(N)
93HeightC = domain.quantities['height'].centroid_values
94DischargeC = domain.quantities['discharge'].centroid_values
95C = domain.centroids
96print 'That took %.2f seconds' %(time.time()-t0)
97X = domain.vertices
98HeightQ = domain.quantities['height'].vertex_values
99VelocityQ = domain.quantities['velocity'].vertex_values
100x = X.flat
101z = domain.quantities['width'].vertex_values.flat
102stage=HeightQ.flat+0
103
104
105
106from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
107
108#hold(False)
109   
110plot1=subplot(211)
111
112plot(x,stage,x,z)
113 
114plot1.set_ylim([-1,11])
115xlabel('Position')
116ylabel('Stage')
117legend(('Solution', 'Width'),
118            'upper right', shadow=True)
119
120plot2=subplot(212)
121plot(x,VelocityQ.flat)
122plot2.set_ylim([-5,5])
123xlabel('Position')
124ylabel('Velocity')
125   
126show()
127
Note: See TracBrowser for help on using the repository browser.