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

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

Putting 1d stuff in one folder

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