source: anuga_work/development/anuga_1d/test_padarn.py @ 7818

Last change on this file since 7818 was 6453, checked in by steve, 15 years ago

Added Padarn's (2008/2009 summer student) files

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