source: trunk/anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_2_padarn.py @ 8204

Last change on this file since 8204 was 8204, checked in by steve, 13 years ago

Getting the 1d tsunami model working

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