source: anuga_work/development/sudi/sw_1d/avalanche/A_velocity/debris_avalanche_RUN.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: 5.4 KB
Line 
1import os
2import time
3from shallow_water_domain_avalanche import *
4from Numeric import array, zeros, Float, sqrt
5from config import g, epsilon
6from numpy import sin, cos, tan
7from scipy import linspace
8from parameters import *
9
10def analytical_sol(X,t):
11    N = len(X)
12    u = zeros(N,Float)
13    h = zeros(N,Float)
14    w = zeros(N,Float)
15    z = zeros(N,Float)
16    mom = zeros(N,Float)
17    for i in range(N):
18        # Calculate Analytical Solution at time t > 0
19        if X[i] <= -2.0*c0*t + 0.5*m*t**2.0:
20            u[i] = 0.0
21            h[i] = 0.0
22        elif X[i] <= c0*t + 0.5*m*t**2.0:
23            u[i] = 2.0/3.0 * (X[i]/t - c0 + m*t)
24            h[i] = 1.0/(9.0*g) * (X[i]/t + 2.0*c0 - 0.5*m*t)**2.0
25        else:
26            u[i] = m*t
27            h[i] = h_0
28        z[i] = bed_slope*X[i]
29        w[i] = h[i] + z[i]
30        mom[i] = u[i]*h[i]
31    return u,h,w,z,mom
32
33def height(X):
34    N = len(X)
35    y = zeros(N,Float)
36    for i in range(N):
37        if X[i]<=0.0:
38            y[i] = 0.0
39        else:
40            y[i] = h_0
41    return y
42
43def stage(X):
44    N = len(X)
45    w = zeros(N,Float)
46    for i in range(N):
47        if X[i]<=0.0:
48            w[i] = bed_slope*X[i]
49        else:
50            w[i] = bed_slope*X[i] + h_0
51    return w       
52
53def elevation(X):
54    N = len(X)
55    y=zeros(N, Float)
56    for i in range(N):
57        y[i] = bed_slope*X[i]
58    return y
59"""
60#=========================================================================#
61#The following values are set in parameters.py
62h_0 = 10.0         # depth upstream. Note that the depth downstream is 0.0
63L = 2000.0         # length of stream/domain
64n = 100             # number of cells
65cell_len = L/n     # length of each cell
66points = zeros(n+1, Float)
67for i in range (n+1):
68    points[i] = i*cell_len - 0.5*L
69
70bed_slope = 0.005      # bottom slope, positive if it is increasing bottom.
71c0 = sqrt(g*h_0)   # sound speed
72m = -1.0*g*bed_slope   # auxiliary variable
73#==========================================================================#
74"""
75boundary = { (0,0): 'left',(n-1,1): 'right'}
76domain = Domain(points,boundary)
77domain.order = 2
78domain.set_timestepping_method('rk2')
79domain.set_CFL(1.0)
80domain.beta = 1.0
81domain.set_limiter("minmod")
82
83def f_right(t):
84    z_r = bed_slope*(0.5*L)
85    h_r = h_0 #+ bed_slope*cell_len
86    w_r = z_r + h_r
87    u_r = m*t   
88    #['stage', 'xmomentum', 'elevation', 'height', 'velocity']
89    return [w_r,  u_r*h_r,  z_r,  h_r,  u_r]
90
91T_right = Time_boundary(domain,f_right)
92#T_right = Transmissive_boundary(domain)
93#D_right = Dirichlet_boundary([bed_slope*(0.5*L)+h_0,  (m*domain.time)*h_0,  bed_slope*(0.5*L),  h_0,  m*domain.time])
94D_left = Dirichlet_boundary([-1.0*bed_slope*(0.5*L),  0.0,  -1.0*bed_slope*(0.5*L),  0.0,  0.0])
95domain.set_boundary({'left':D_left,'right':T_right})
96domain.set_quantity('stage',stage)
97domain.set_quantity('elevation',elevation)
98
99X = domain.vertices
100C = domain.centroids
101
102yieldstep = finaltime = 15.0
103t0 = time.time()
104
105while finaltime <= 15.1:
106    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
107        domain.write_time()
108       
109    #The following is for computing the error and plotting the result
110    Mom = domain.quantities['xmomentum']
111    Height = domain.quantities['height']
112    Stage = domain.quantities['stage']
113    Velocity = domain.quantities['velocity']
114    Elevation = domain.quantities['elevation']
115
116    #The following is for computing the error
117    Uc,Hc,Wc,Zc,Mc = analytical_sol(C,domain.time)
118    HeightC = Height.centroid_values
119    MomC = Mom.centroid_values
120    StageC = Stage.centroid_values
121    VelC = Velocity.centroid_values
122    ElevationC = Elevation.centroid_values
123
124    print "number of cells=",n
125    W_error = (1.0/n)*sum(abs(Wc-StageC))
126    M_error = (1.0/n)*sum(abs(Mc-MomC))
127    U_error = (1.0/n)*sum(abs(Uc-VelC))
128    print "stage_error %.10f" %(W_error)
129    print "momentum_error %.10f"%(M_error)
130    print "velocity_error %.10f" %(U_error)
131
132
133    #The following is for plotting the result.
134    Uv,Hv,Wv,Zv,Mv = analytical_sol(X.flat,domain.time)
135    HeightV = Height.vertex_values
136    MomV = Mom.vertex_values
137    StageV = Stage.vertex_values
138    VelV = Velocity.vertex_values
139    ElevationV = Elevation.vertex_values
140    from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
141    hold(False)
142    clf()
143
144    plot1 = subplot(311)
145    plot(X.flat,Wv,'b-', X.flat,StageV.flat,'k--', X.flat,ElevationV.flat,'k:')
146    #plot(X.flat,Wv, X.flat,StageV.flat, X.flat,ElevationV.flat)   
147    #xlabel('Position')
148    ylabel('Stage')
149    #plot1.set_ylim([-1.0,21.0])
150    #plot1.set_xlim([-480.0,-420.0])#([-9.0e-3,9.0e-3])
151    #legend(('analytical solution', 'numerical solution', 'discretized bed'), 'upper left', shadow=False)
152   
153    plot2 = subplot(312)
154    plot(X.flat,Mv,'b-', X.flat,MomV.flat,'k--')
155    #xlabel('Position')
156    ylabel('Momentum')
157    #plot2.set_xlim([-300.0,300.0])
158    plot2.set_ylim([-310.0,10.0])#([-90.0,10.0])
159    #legend(('analytical solution', 'numerical solution'), 'lower right', shadow=False)   
160
161    plot3 = subplot(313)
162    plot(X.flat,Uv,'b-', X.flat,VelV.flat,'k--')
163    xlabel('Position')
164    ylabel('Velocity')
165    #plot3.set_xlim([-300.0,300.0])
166    plot3.set_ylim([-45.0,5.0])#([-30.0,5.0])
167    #legend(('analytical solution', 'numerical solution'), 'lower right', shadow=False)
168   
169    finaltime = finaltime + 10.0
170
171    file = "A-case3-"
172    file += str(n)
173    file += ".eps"
174    savefig(file)     
175   
176    show()
177
178
Note: See TracBrowser for help on using the repository browser.