source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/cg/parameter.py @ 7933

Last change on this file since 7933 was 7933, checked in by mungkasi, 14 years ago

Renaming a module "bisect.py" to "bisect_function.py" since it crashes with "bisect" in pylab.

File size: 2.2 KB
RevLine 
[7922]1from scipy.special import jn
2from scipy import pi, sqrt, linspace, pi, sin, cos
3from config import g
4from Numeric import Float
5from numpy import zeros
6from rootsearch import *
[7933]7from bisect_function import *
[7922]8
9def j0(x):
10    return jn(0.0, x)
11
12def j1(x):
13    return jn(1.0, x)
14
15def j2(x):
16    return jn(2.0, x)
17
18def j3(x):
19    return jn(3.0, x)
20
21def jm1(x):
22    return jn(-1.0, x)
23
24def jm2(x):
25    return jn(-2.0, x)
26
27
28"""
29##==========================================================================##
30#DIMENSIONAL PARAMETERS
31L = 5e4                 # Length of channel (m)
32h_0 = 5e2               # Height at origin when the water is still
33N = 550                 # Number of computational cells????????????????
34cell_len = 1.1*L/N      # Origin = 0.0
35Tp = 15.0*60.0                           # Period of oscillation
36a = 1.0                                 # Amplitude at origin
37X = linspace(-0.5*cell_len, 1.1*L+0.5*cell_len, N+2)             # Discretized spatial domain
38##=========================================================================##
39"""
40#DIMENSIONAL PARAMETERS
41L = 5e4                 # Length of channel (m)
42h_0 = 5e2               # Height at origin when the water is still
43numb = 550                 # initial Number of computational cells
44cell_len = 100.0#100.0 is actually 1.1*L/numb
45Tp = 15.0*60.0                           # Period of oscillation
46a =1.0                                  # Amplitude at origin
47               
48points = zeros(numb+2,Float)
49for i in range(numb+2):
50        points[i] = i*cell_len - 0.5*cell_len
51N = len(points)-1       # Number of computational cells
52
53
54#DIMENSIONLESS PARAMETERS
55eps = a/h_0
56T = Tp*sqrt(g*h_0)/L
57A = eps/j0(4.0*pi/T)
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.