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

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

Moving 2010 project

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