source: anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_1_padarn.py @ 7839

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

Changing name of 1d projects so that it will be easy to moveto the trunk

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