source: anuga_work/development/2010-projects/anuga_1d/sww/parabolic_cannal_vel.py @ 7839

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

Changing name of 1d projects so that it will be easy to moveto the trunk

File size: 4.4 KB
Line 
1import os
2from math import sqrt, pi, sin, cos
3from shallow_water_vel_domain import *
4from numpy import allclose, array, zeros, ones, numpy.float, take, sqrt
5from config import g, epsilon
6
7def analytic_cannal(C,t):
8    N = len(C)
9    u = zeros(N,numpy.float)    ## water velocity
10    h = zeros(N,numpy.float)    ## water depth
11    x = C
12    g = 9.81
13    ## Define Basin Bathymetry
14    z_b = zeros(N,numpy.float) ## elevation of basin
15    z = zeros(N,numpy.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,numpy.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),numpy.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,numpy.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),numpy.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', 'numpyal 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.