source: anuga_work/development/anuga_1d/channel_flow_2_padarn.py @ 7818

Last change on this file since 7818 was 6453, checked in by steve, 15 years ago

Added Padarn's (2008/2009 summer student) files

File size: 2.5 KB
Line 
1import os
2from math import sqrt, pow, pi
3from channel_domain import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6
7
8print "Channel Flow 1 Padarn Test"
9
10# Define functions for initial quantities
11def initial_area(x):
12    return 1.4691*width(x)
13
14def width(x):
15    x1=(x/1000)*(x/1000)
16    x2=x1*(x/1000)
17    x3=x2*(x/1000)
18    return 10-64*(x1-2*x2+x3)
19
20def bed(x):
21    y = zeros(len(x),Float)
22    for i in range(len(x)):
23        if x[i]<525 and x[i]>475:
24            y[i]=0
25        else:
26            y[i]=0
27    return y
28 
29
30def initial_discharge(x):
31    return 20
32 
33import time
34
35# Set final time and yield time for simulation
36finaltime = 10.0
37yieldstep = finaltime
38
39# Length of channel (m)
40L = 1000.0   
41# Define the number of cells
42number_of_cells = [200]
43
44# Define cells for finite volume and their size
45N = int(number_of_cells[0])     
46print "Evaluating domain with %d cells" %N
47cell_len = L/N # Origin = 0.0
48points = zeros(N+1,Float)
49
50# Define the centroid points
51for j in range(N+1):
52    points[j] = j*cell_len
53
54# Create domain with centroid points as defined above       
55domain = Domain(points)
56
57# Set initial values of quantities - default to zero
58domain.set_quantity('area', initial_area)
59domain.set_quantity('width',width)
60domain.set_quantity('elevation',bed)
61domain.set_quantity('discharge',initial_discharge)
62
63# Set boundry type, order, timestepping method and limiter
64domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])})
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
74for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
75    domain.write_time()
76
77N = float(N)
78HeightC = domain.quantities['height'].centroid_values
79DischargeC = domain.quantities['discharge'].centroid_values
80C = domain.centroids
81print 'That took %.2f seconds' %(time.time()-t0)
82X = domain.vertices
83HeightQ = domain.quantities['height'].vertex_values
84VelocityQ = domain.quantities['velocity'].vertex_values
85x = X.flat
86z = domain.quantities['elevation'].vertex_values.flat
87stage=HeightQ.flat+z
88
89from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
90
91hold(False)
92   
93plot1 = subplot(211)
94
95plot(x,z,x,stage)
96 
97plot1.set_ylim([-1,11])
98xlabel('Position')
99ylabel('Stage')
100legend(('Analytical Solution', 'Numerical Solution'),
101           'upper right', shadow=True)
102plot2 = subplot(212)
103plot(x,VelocityQ.flat)
104plot2.set_ylim([-10,10])
105   
106xlabel('Position')
107ylabel('Velocity')
108   
109show()
110   
Note: See TracBrowser for help on using the repository browser.