source: anuga_work/development/sudi/sw_1d/periodic_waves/johns/analytical_prescription.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.0 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 = -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 array, zeros, Float
30    if x<1.0:
31        q = zeros(2,Float)
32    else:
33        q = shore(t)
34    def fun(q): #Here q=(w, u)
35        f = zeros(2,Float)
36        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]))
37        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
38        return f
39    def newtonRaphson2(f,q,tol=1.0e-15):
40        for i in range(30):                                                                                                                                   
41            h = 1.0e-4      ##1.0e-4 may be too large.
42            n = len(q)
43            jac = zeros((n,n),Float)
44            if 1.0+q[0]-x<0.0:
45                #temp1 = 1.0+q[0]-x
46                #q[0] = q[0]-temp1
47                q[1] = SomethingWRONG
48                #return q
49                #q[0] =0.25
50            f0 = f(q)
51            for i in range(n):
52                temp = q[i]
53                q[i] = temp + h
54                f1 = f(q)
55                q[i] = temp
56                jac[:,i] = (f1 - f0)/h
57            if sqrt(dot(f0,f0)/len(q)) < tol: return q
58            dq = gaussPivot(jac,-f0)
59            q = q + dq
60            if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q
61        print 'Too many iterations'
62    q = newtonRaphson2(fun,q)   
63    return q[0], q[1]
64
65#w,u=prescribe(0.0,0.0)
Note: See TracBrowser for help on using the repository browser.