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

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

Moving 2010 project

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