source: anuga_work/development/flow_1d/channel_flow/channel_flow_1_padarn.py @ 7833

Last change on this file since 7833 was 7833, checked in by steve, 12 years ago
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.