source: anuga_work/development/2010-projects/anuga_1d/channel/test_padarn.py @ 7839

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

Changing name of 1d projects so that it will be easy to moveto the trunk

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.