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

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

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

File size: 1.9 KB
RevLine 
[7922]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 = -0.2#-1.0
19    b = 0.2#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))
23    #position = 1.0 + xi
24    return [xi, u]#[position, u]
25
26
27
28def prescribe(x,t):
29    from Numeric import Float
30    from numpy import array, zeros
31    if x<1.0:
32        q = zeros(2,Float)
33    else:
34        q = shore(t)
35    def fun(q): #Here q=(w, u)
36        f = zeros(2,Float)
37        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]))
38        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
39        return f
40    def newtonRaphson2(f,q,tol=1.0e-15):
41        for i in range(30):                                                                                                                                   
42            h = 1.0e-4      ##1.0e-4 may be too large.
43            n = len(q)
44            jac = zeros((n,n),Float)
45            if 1.0+q[0]-x<0.0:
46                #temp1 = 1.0+q[0]-x
47                #q[0] = q[0]-temp1
48                q[1] = SomethingWRONG
49                #return q
50                #q[0] =0.25
51            f0 = f(q)
52            for i in range(n):
53                temp = q[i]
54                q[i] = temp + h
55                f1 = f(q)
56                q[i] = temp
57                jac[:,i] = (f1 - f0)/h
58            if sqrt(dot(f0,f0)/len(q)) < tol: return q
59            dq = gaussPivot(jac,-f0)
60            q = q + dq
61            if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q
62        print 'Too many iterations'
63    q = newtonRaphson2(fun,q)   
64    return q[0], q[1]
65
66#w,u=prescribe(0.0,0.0)
Note: See TracBrowser for help on using the repository browser.