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

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

Moving 2010 project

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