source: anuga_work/development/sudi/sw_1d/periodic_waves/johns/parameter.py @ 7837

Last change on this file since 7837 was 7837, checked in by mungkasi, 12 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

File size: 2.2 KB
Line 
1from scipy.special import jn
2from scipy import pi, sqrt, linspace, pi, sin, cos
3from config import g
4from Numeric import zeros, Float
5from rootsearch import *
6from bisect import *
7
8def j0(x):
9    return jn(0.0, x)
10
11def j1(x):
12    return jn(1.0, x)
13
14def j2(x):
15    return jn(2.0, x)
16
17def j3(x):
18    return jn(3.0, x)
19
20def jm1(x):
21    return jn(-1.0, x)
22
23def jm2(x):
24    return jn(-2.0, x)
25
26
27"""
28##==========================================================================##
29#DIMENSIONAL PARAMETERS
30L = 5e4                 # Length of channel (m)
31h_0 = 5e2               # Height at origin when the water is still
32N = 550                 # Number of computational cells????????????????
33cell_len = 1.1*L/N      # Origin = 0.0
34Tp = 15.0*60.0                           # Period of oscillation
35a = 1.0                                 # Amplitude at origin
36X = linspace(-0.5*cell_len, 1.1*L+0.5*cell_len, N+2)             # Discretized spatial domain
37##=========================================================================##
38"""
39#DIMENSIONAL PARAMETERS
40L = 5e4                 # Length of channel (m)
41h_0 = 5e2               # Height at origin when the water is still
42numb = 550                 # initial Number of computational cells
43cell_len = 100.0#100.0 is actually 1.1*L/numb
44Tp = 15.0*60.0                           # Period of oscillation
45a = 1.0                                  # Amplitude at origin
46               
47points = zeros(numb+2,Float)
48for i in range(numb+2):
49        points[i] = i*cell_len - 0.5*cell_len
50N = len(points)-1       # Number of computational cells
51
52
53#DIMENSIONLESS PARAMETERS
54eps = a/h_0
55T = Tp*sqrt(g*h_0)/L
56A = eps/j0(4.0*pi/T)
57
58
59def w_at_O_JOHNS(t):
60    return eps*cos(2.0*pi*t/T)
61
62def u_at_O_JOHNS(t):
63    a = -1.01#-1.0
64    b = 1.01#1.0
65    dx = 0.01
66    w = w_at_O_JOHNS(t)
67    def fun(u):
68        return u + A*j1(4.0*pi/T*(1.0+w)**0.5)*sin(2.0*pi/T*(t+u))/(1+w)**0.5
69    while 1:
70        x1,x2 = rootsearch(fun,a,b,dx)
71        if x1 != None:
72            a = x2
73            root = bisect(fun,x1,x2,1)
74        else:
75            break
76    return root
77
78
79def prescribe_at_O_JOHNS(t):
80    w = w_at_O_JOHNS(t)
81    u = u_at_O_JOHNS(t)
82    return w, u
Note: See TracBrowser for help on using the repository browser.