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