source: anuga_work/development/shallow_water_1d/parabolic_cannal_h.py @ 4959

Last change on this file since 4959 was 4959, checked in by steve, 16 years ago

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

File size: 4.9 KB
Line 
1import os
2from math import sqrt, pi, sin, cos
3from shallow_water_h import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6
7def newLinePlot(title='Simple Plot'):
8    import Gnuplot
9    g = Gnuplot.Gnuplot(persist=0)
10    g.title(title)
11    g('set data style linespoints') 
12    g.xlabel('x')
13    g.ylabel('y')
14    return g
15
16def linePlot(g,x1,y1,x2,y2,x3,y3):
17    import Gnuplot
18    plot1 = Gnuplot.PlotItems.Data(x1.flat,y1.flat,with="lines 2")
19    plot2 = Gnuplot.PlotItems.Data(x2.flat,y2.flat,with="lines 3")
20    plot3 = Gnuplot.PlotItems.Data(x3.flat,y3.flat,with="linespoints 1")
21    g.plot(plot1,plot2,plot3)
22
23
24def analytic_cannal(C,t):
25    N = len(C)
26    u = zeros(N,Float)    ## water velocity
27    h = zeros(N,Float)    ## water depth
28    x = C
29    g = 9.81
30
31
32    ## Define Basin Bathymetry
33    z_b = zeros(N,Float) ## elevation of basin
34    z = zeros(N,Float)   ## elevation of water surface
35    z_infty = 10.0       ## max equilibrium water depth at lowest point.
36    L_x = 2500.0         ## width of channel
37
38    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
39    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
40
41    for i in range(N):
42        z_b[i] = z_infty*(x[i]**2/L_x**2) ## or A0*cos(omega*t)\pmL_x
43        u[i] = -A0*omega*sin(omega*t)
44        z[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
45
46    h = z-z_b
47    return u,h,z,z_b
48
49
50#plot2 = newLinePlot("Momentum") 
51
52L_x = 2500.0     # Length of channel (m)
53N = 100          # Number of compuational cells
54cell_len = 4*L_x/N   # Origin = 0.0
55
56points = zeros(N+1,Float)
57for i in range(N+1):
58        points[i] = -2*L_x +i*cell_len
59
60domain = Domain(points)
61
62domain.default_order = 2
63domain.default_time_order = 2
64domain.cfl = 1.0
65domain.beta = 1.0
66domain.flux_first = False
67
68
69
70def stage(x):
71
72    z_infty = 10.0       ## max equilibrium water depth at lowest point.
73    L_x = 2500.0         ## width of channel
74    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
75    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
76    t=0.0
77    y = zeros(len(x),Float)
78    for i in range(len(x)):
79        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))
80    return y
81
82def elevation(x):
83    N = len(x)
84    z_infty = 10.0
85    z = zeros(N,Float)
86    L_x = 2500.0
87    A0 = 0.5*L_x
88    omega = sqrt(2*g*z_infty)/L_x
89    for i in range(N):
90        z[i] = z_infty*(x[i]**2/L_x**2)
91    return z
92
93def height(x):
94    z_infty = 10.0       ## max equilibrium water depth at lowest point.
95    L_x = 2500.0         ## width of channel
96    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
97    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
98    t=0.0
99    y = zeros(len(x),Float)
100    for i in range(len(x)):
101        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)
102        #y[i] = 12.0-z_infty*(x[i]**2/L_x**2)
103    return y
104
105   
106#domain.set_quantity('stage', stage)
107domain.set_quantity('elevation',elevation)
108domain.set_quantity('height',height)
109#print "Centroid Heights"
110#print domain.quantities['height'].centroid_values
111#print "Vertex Heights"
112#print domain.quantities['height'].vertex_values
113#print "Centroids Stage"
114#print domain.quantities['height'].centroid_values+domain.quantities['elevation'].centroid_values
115#raw_input()
116domain.set_boundary({'exterior': Reflective_boundary(domain)})
117 
118import time
119t0 = time.time()
120#yieldstep = 50.0
121yieldstep = 1122.0
122#yieldstep = 1.0
123finaltime = 1122.0
124
125plot1 = newLinePlot("Stage")
126plot2 = newLinePlot("Xmomentum")
127C = domain.centroids
128X = domain.vertices
129HeightQ = domain.quantities['height'].vertex_values
130ElevationQ = domain.quantities['elevation'].vertex_values
131XmomQ = domain.quantities['xmomentum'].vertex_values
132
133#StageQ = domain.quantities['stage'].centroid_values       
134StageQ = ElevationQ+HeightQ
135u,h,z,z_b = analytic_cannal(X.flat,0.0)
136linePlot(plot1,X,z,X,z_b,X,StageQ)
137raw_input()
138print "integral height %f"%domain.quantities['height'].get_integral()
139print "integral xmom %f"%domain.quantities['xmomentum'].get_integral()
140C = domain.centroids
141import time
142t0 = time.time()
143for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
144    domain.write_time()
145    print "integral height %f"%domain.quantities['height'].get_integral()
146    print "integral xmom %f"%domain.quantities['xmomentum'].get_integral()
147    u,h,z,z_b = analytic_cannal(X.flat,t)
148    StageQ = HeightQ+ElevationQ
149    #print "yield",HeightQ
150    linePlot(plot1,X,z,X,z_b,X,StageQ)
151    linePlot(plot2,X,u*h,X,u*h,X,XmomQ)
152    u,hc,z,z_b = analytic_cannal(C,t)
153    HeightQC = domain.quantities['height'].centroid_values
154    #print HeightQC
155    for k in range(N):
156        if hc[k] < 0.0:
157            hc[k] = 0.0
158    error = 1.0/(N)*sum(abs(hc-HeightQC))
159    print 'Error measured at centroids', error
160   # raw_input()
161   
162print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.