source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/run_parabolic_canal.py

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

Channel flow reconstruction fixed up (width was being cutoff)

File size: 3.2 KB
Line 
1import os
2from math import sqrt, pi, sin, cos
3import numpy
4import time
5
6
7from anuga_1d.sww.sww_domain import *
8from anuga_1d.config import g, epsilon
9from anuga_1d.base.generic_mesh import uniform_mesh
10
11#===============================================================================
12# setup problem
13#===============================================================================
14
15
16z_infty = 10.0       ## max equilibrium water depth at lowest point.
17L_x = 2500.0         ## width of channel
18A0 = 0.5*L_x                  ## determines amplitudes of oscillations
19omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
20
21def analytic_canal(x,t):
22    u = numpy.zeros_like(x)    ## water velocity
23    h = numpy.zeros_like(x)    ## water depth
24
25    ## Define Basin Bathymetry
26    z = numpy.zeros_like(x) ## elevation of basin
27    w = numpy.zeros_like(x)   ## elevation of water surface
28
29    z[:] = z_infty*(x**2/L_x**2)
30    u[:] = -A0*omega*sin(omega*t)
31    w[:] = numpy.maximum(z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x-0.5*A0/(L_x)*cos(omega*t)),z)
32    h[:] = numpy.maximum(w-z, 0.0)
33
34    T = 2.0*pi/omega
35
36    return u,h,w,z, T
37
38
39def stage(x):
40    t=0.0
41    u,h,w,z,T = analytic_canal(x,t)
42
43    return w
44
45def elevation(x):
46    t=0.0
47    u,h,w,z,T = analytic_canal(x,t)
48
49    return z
50
51
52def height(x):
53    t=0.0
54    u,h,w,z,T = analytic_canal(x,t)
55
56    return h
57
58
59#===============================================================================
60finaltime = 500.0
61yieldstep = 10.0
62
63
64
65N = 100
66print "Evaluating domain with %d cells" %N
67domain = Domain(*uniform_mesh(N, x_0 = -2.0*L_x, x_1 = 2.0*L_x))
68
69domain.set_spatial_order(2)
70domain.set_timestepping_method('rk2')
71domain.set_CFL(1.0)
72#domain.set_limiter("minmod_kurganov")
73domain.set_limiter("vanleer")
74
75domain.set_beta(1.0)
76
77domain.set_quantity('stage', stage)
78domain.set_quantity('elevation',elevation)
79
80domain.quantities['elevation'].extrapolate_second_order()
81domain.quantities['stage'].extrapolate_second_order()
82
83Br = Reflective_boundary(domain)
84
85domain.set_boundary({'left': Br, 'right' : Br})
86
87#domain.h0=0.0001
88
89t0 = time.time()
90
91for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
92    domain.write_time()
93
94print 'That took %.2f seconds' %(time.time()-t0)
95
96
97N = float(N)
98HeightC    = domain.quantities['height'].centroid_values
99StageC     = domain.quantities['stage'].centroid_values
100BedC       = domain.quantities['elevation'].centroid_values
101C          = domain.centroids
102
103HeightV    = domain.quantities['height'].vertex_values
104StageV     = domain.quantities['stage'].vertex_values
105BedV       = domain.quantities['elevation'].vertex_values
106VelocityV  = domain.quantities['velocity'].vertex_values
107X          = domain.vertices
108
109
110import pylab
111
112pylab.hold(False)
113
114plot1 = pylab.subplot(311)
115
116pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat)
117
118plot1.set_ylim([-1,40])
119
120pylab.xlabel('Position')
121pylab.ylabel('Stage')
122pylab.legend(('Analytical Solution', 'Numerical Solution'),
123           'upper right', shadow=True)
124
125
126plot2 = pylab.subplot(312)
127
128pylab.plot(X.flat,HeightV.flat)
129
130plot2.set_ylim([-1,10])
131
132pylab.xlabel('Position')
133pylab.ylabel('Height')
134
135plot3 = pylab.subplot(313)
136
137pylab.plot(X.flat,VelocityV.flat)
138plot3.set_ylim([-15,15])
139
140pylab.xlabel('Position')
141pylab.ylabel('Velocity')
142
143pylab.show()
144
145
146
147
148
149
Note: See TracBrowser for help on using the repository browser.