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

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

Moving 2010 project

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