source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/cg/analytical_prescription.py @ 7922

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

Modifying codes so that arrays are compatible with numpy instead of Numeric.

File size: 2.5 KB
Line 
1from scipy import sin, cos, sqrt, pi, dot
2from gaussPivot import *
3from parameter import *
4
5def root_g(a,b,t):
6    dx = 0.01
7    def g(u):
8        return u + 2.0*A*pi/T*sin(2.0*pi/T*(t+u))
9    while 1:
10        x1,x2 = rootsearch(g,a,b,dx)
11        if x1 != None:
12            a = x2
13            root = bisect(g,x1,x2,1)
14        else:
15            break
16    return root
17def shore(t):
18    a = -1.0
19    b = 1.0
20    #dx = 0.01
21    u = root_g(a,b,t)
22    xi = -0.5*u*u + A*cos(2.0*pi/T*(t+u)) #This is the displacement
23    #position = 1.0 + xi
24    return [xi, u]#[position, u]
25
26def w_at_O(t):
27    return eps*cos(2.0*pi*t/T)
28
29def u_at_O(t):
30    a = -1.01
31    b = 1.01
32    dx = 0.01
33    w = w_at_O(t)
34    def fun(u):
35        return u + A*j1(4.0*pi/T*(1.0+w)**0.5)*sin(2.0*pi/T*(t+u))/(1.0+w)**0.5
36    while 1:
37        x1,x2 = rootsearch(fun,a,b,dx)
38        if x1 != None:
39            a = x2
40            root = bisect(fun,x1,x2,1)
41        else:
42            break
43    return root
44
45
46def prescribe(x,t):
47    from Numeric import Float
48    from numpy import array, zeros
49    q_shore = shore(t)
50    q = zeros(2,Float)
51    if x<1.0-abs(q_shore[0]):
52        q[0] = w_at_O(t)
53        q[1] = u_at_O(t)
54    else:
55        q = q_shore
56    def fun(q): #Here q=(w, u)
57        f = zeros(2,Float)
58        f[0] = q[0] + 0.5*q[1]**2.0 - A*j0(4.0*pi/T*(1.0+q[0]-x)**0.5)*cos(2.0*pi/T*(t+q[1]))
59        f[1] = q[1] + A*j1(4.0*pi/T*(1.0+q[0]-x)**0.5)*sin(2.0*pi/T*(t+q[1]))/(1+q[0]-x)**0.5
60        return f
61    def newtonRaphson2(f,q,tol=1.0e-15):
62        for i in range(30):                                                                                                                                   
63            h = 1.0e-4      ##1.0e-4 may be too large.
64            n = len(q)
65            jac = zeros((n,n),Float)
66            if 1.0+q[0]-x<0.0:
67                #temp1 = 1.0+q[0]-x
68                #q[0] = x-1.0+0.01
69                #print "1.0+q[0]-x=",1.0+q[0]-x
70                print "SomethingWRONG"
71                #return q
72                #q[0] =0.25
73            f0 = f(q)
74            for i in range(n):
75                temp = q[i]
76                q[i] = temp + h
77                f1 = f(q)
78                q[i] = temp
79                jac[:,i] = (f1 - f0)/h
80            if sqrt(dot(f0,f0)/len(q)) < tol: return q
81            dq = gaussPivot(jac,-f0)
82            q = q + dq
83            if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q
84        print 'Too many iterations'
85    q = newtonRaphson2(fun,q)   
86    return q[0], q[1]
87
88#w,u=prescribe(0.0,0.0)
Note: See TracBrowser for help on using the repository browser.