source: anuga_work/development/anuga_1d/parabolic_cannal_vel.py @ 7825

Last change on this file since 7825 was 5845, checked in by steve, 15 years ago

Adding shallow_water_vel code

File size: 4.4 KB
Line 
1import os
2from math import sqrt, pi, sin, cos
3from shallow_water_vel_domain import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6
7def analytic_cannal(C,t):
8    N = len(C)
9    u = zeros(N,Float)    ## water velocity
10    h = zeros(N,Float)    ## water depth
11    x = C
12    g = 9.81
13    ## Define Basin Bathymetry
14    z_b = zeros(N,Float) ## elevation of basin
15    z = zeros(N,Float)   ## elevation of water surface
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    for i in range(N):
21        z_b[i] = z_infty*(x[i]**2/L_x**2)
22        u[i] = -A0*omega*sin(omega*t)
23        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))
24    h = z-z_b
25    T = 2.0*pi/omega
26    return u,h,z,z_b, T
27
28
29L_x = 2500.0     # Length of channel (m)
30N = 400         # Number of compuational cells
31cell_len = 4*L_x/N   # Origin = 0.0
32
33points = zeros(N+1,Float)
34for i in range(N+1):
35        points[i] = -2*L_x +i*cell_len
36
37def stage(x):
38    z_infty = 10.0       ## max equilibrium water depth at lowest point.
39    L_x = 2500.0         ## width of channel
40    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
41    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
42    t=0.0
43    y = zeros(len(x),Float)
44    for i in range(len(x)):
45        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))
46        #y[i] = 12.0
47    return y
48
49def elevation(x):
50    N = len(x)
51    z_infty = 10.0
52    z = zeros(N,Float)
53    L_x = 2500.0
54    A0 = 0.5*L_x
55    omega = sqrt(2*g*z_infty)/L_x
56    for i in range(N):
57        z[i] = z_infty*(x[i]**2/L_x**2)
58    return z
59
60def height(x):
61    z_infty = 10.0       ## max equilibrium water depth at lowest point.
62    L_x = 2500.0         ## width of channel
63    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
64    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
65    t=0.0
66    y = zeros(len(x),Float)
67    for i in range(len(x)):
68        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)
69    return y
70
71domain = Domain(points)
72domain.order = 2
73domain.set_timestepping_method('euler')
74domain.set_CFL(1.0)
75domain.beta = 1.0
76domain.set_limiter("vanleer")
77
78   
79domain.set_quantity('stage', stage)
80domain.set_quantity('elevation',elevation)
81domain.set_boundary({'exterior': Reflective_boundary(domain)})
82 
83X = domain.vertices
84u,h,z,z_b,T = analytic_cannal(X.flat,domain.time)
85print 'T = ',T
86
87yieldstep = finaltime = 10.0 #T/2.0
88StageQ = domain.quantities['stage'].vertex_values
89XmomQ = domain.quantities['xmomentum'].vertex_values
90
91import time
92t0 = time.time()
93for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
94    domain.write_time()
95    print "integral", domain.quantities['stage'].get_integral()
96    if t>0.0:
97       
98        x = X.flat
99        print 'domain.time=',domain.time
100       
101        w_V = domain.quantities['stage'].vertex_values.flat
102        uh_V = domain.quantities['xmomentum'].vertex_values.flat
103        u_V = domain.quantities['velocity'].vertex_values.flat
104        h_V = domain.quantities['height'].vertex_values.flat
105       
106        u,h,z,z_b,T = analytic_cannal(x,domain.time)
107        w = z
108        for k in range(len(h)):
109            if h[k] < 0.0:
110                h[k] = 0.0
111                w[k] = z_b[k]
112       
113        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
114        hold(False)
115
116        #print 'size X',X.shape
117        #print 'size w',w.shape
118        signh = (h_V>0)
119
120        plot1 = subplot(311)
121        #plot(x,w,  x,w_V,  x,z_b)
122        plot(x,signh)
123        plot1.set_xlim([-6000,6000])
124        xlabel('Position')
125        ylabel('Stage')
126        legend(('Analytic Solution', 'Numerical Solution', 'Bed'),
127               'upper center', shadow=True)
128
129       
130        plot2 = subplot(312)
131        plot(x,u*h,x,uh_V)
132        #plot2.set_ylim([-1,25])
133        xlabel('Position')
134        ylabel('Xmomentum')
135
136        print u_V
137       
138        plot3 = subplot(313)
139        plot(x,u, x,u_V)
140        #plot2.set_ylim([-1,25])
141        xlabel('Position')
142        ylabel('Velocity')
143        show()
144raw_input("Press the return key!")
Note: See TracBrowser for help on using the repository browser.