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

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

Moving 2010 project

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