source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/cg/run-discrepancy.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.7 KB
Line 
1import os
2from scipy.special import jn
3from scipy import sin, cos, sqrt, linspace, pi, zeros
4from rootsearch import *
5from bisect_function import *
6from Numeric import Float
7from numpy import zeros,dot
8from gaussPivot import *
9from config import g
10from analytical_prescription import *
11
12
13def j0(x):
14    return jn(0.0, x)
15
16def j1(x):
17    return jn(1.0, x)
18
19def j2(x):
20    return jn(2.0, x)
21
22def j3(x):
23    return jn(3.0, x)
24
25def jm1(x):
26    return jn(-1.0, x)
27
28def jm2(x):
29    return jn(-2.0, x)
30
31def bed(x):
32    return x-1.0
33
34
35def w_at_O(t):
36    return eps*cos(2.0*pi*t/T)
37
38def u_at_O(t):
39    a = -1.01#-1.0
40    b = 1.01#1.0
41    dx = 0.01
42    w = w_at_O(t)
43    def fun(u):
44        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
45    while 1:
46        x1,x2 = rootsearch(fun,a,b,dx)
47        if x1 != None:
48            a = x2
49            root = bisect(fun,x1,x2,1)
50        else:
51            break
52    return root
53
54"""
55##==========================================================================##
56#DIMENSIONAL PARAMETERS
57L = 5e4                 # Length of channel (m)
58h_0 = 5e2               # Height at origin when the water is still
59Tp = 15.0*60.0                           # Period of oscillation
60a = 1.0                                 # Amplitude at origin
61##=========================================================================##
62#DIMENSIONLESS PARAMETERS
63eps = a/h_0
64T = Tp*sqrt(g*h_0)/L
65A = eps/j0(4.0*pi/T)
66"""
67
68
69Time = linspace(0.0,T,100)
70N_T = len(Time)
71
72
73Stage = zeros(N_T, Float)   
74Veloc = zeros(N_T, Float)
75for i in range(N_T):
76    t=Time[i]
77    zet, vel = prescribe(0.0,t)
78    Stage[i] = zet
79    Veloc[i] = vel
80
81
82Stage_johns = zeros(N_T, Float)   
83Veloc_johns = zeros(N_T, Float)
84for i in range(N_T):
85    t=Time[i]
86    Stage_johns[i] = w_at_O(t)
87    Veloc_johns[i] = u_at_O(t)
88
89"""
90num=len(Stage)
91error_w=(1.0/num)*sum(abs(Stage-Stage_johns))*h_0
92error_u=(1.0/num)*sum(abs(Veloc-Veloc_johns))*sqrt(g*h_0)
93print "error_w=", error_w
94print "error_u=", error_u
95
96"""
97
98from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
99
100hold(False)
101clf()
102plot1 = subplot(211)
103plot(Time/T,Stage*h_0,'b-', Time/T,Stage_johns*h_0,'k--')
104xlabel('t/T')
105ylabel('Stage')
106#plot1.set_xlim([0.000,0.030])
107#plot1.set_ylim([0.980,1.005]) #([-9.0e-3,9.0e-3])
108legend(('C-G', 'Johns'),
109       'lower left', shadow=False)
110
111plot2 = subplot(212)
112plot(Time/T,Veloc*sqrt(g*h_0),'b-', Time/T,Veloc_johns*sqrt(g*h_0),'k--')
113xlabel('t/T')
114ylabel('Velocity')
115#plot1.set_xlim([0.0,1.1])
116#plot2.set_ylim([-0.05,0.05]) #([-1.0e-12,1.0e-12])
117legend(('C-G', 'Johns'),
118       'upper right', shadow=False)
119
120
121#filename = "discrepancy-closer"
122#filename += str(i)
123#filename += ".eps"
124#savefig(filename)
125#show()
126
127#plot(Time,Vel_at_O)
128#show()
Note: See TracBrowser for help on using the repository browser.