source: anuga_work/development/sudi/sw_1d/periodic_waves/cg/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.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 array, zeros, Float
48    q_shore = shore(t)
49    q = zeros(2,Float)
50    if x<1.0-abs(q_shore[0]):
51        q[0] = w_at_O(t)
52        q[1] = u_at_O(t)
53    else:
54        q = q_shore
55    def fun(q): #Here q=(w, u)
56        f = zeros(2,Float)
57        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]))
58        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
59        return f
60    def newtonRaphson2(f,q,tol=1.0e-15):
61        for i in range(30):                                                                                                                                   
62            h = 1.0e-4      ##1.0e-4 may be too large.
63            n = len(q)
64            jac = zeros((n,n),Float)
65            if 1.0+q[0]-x<0.0:
66                #temp1 = 1.0+q[0]-x
67                #q[0] = x-1.0+0.01
68                #print "1.0+q[0]-x=",1.0+q[0]-x
69                print "SomethingWRONG"
70                #return q
71                #q[0] =0.25
72            f0 = f(q)
73            for i in range(n):
74                temp = q[i]
75                q[i] = temp + h
76                f1 = f(q)
77                q[i] = temp
78                jac[:,i] = (f1 - f0)/h
79            if sqrt(dot(f0,f0)/len(q)) < tol: return q
80            dq = gaussPivot(jac,-f0)
81            q = q + dq
82            if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q
83        print 'Too many iterations'
84    q = newtonRaphson2(fun,q)   
85    return q[0], q[1]
86
87#w,u=prescribe(0.0,0.0)
Note: See TracBrowser for help on using the repository browser.