source: trunk/anuga_work/development/sudi/sw_1d/avalanche/B_momentum/debris_avalanche_RUN.py @ 7916

Last change on this file since 7916 was 7916, checked in by mungkasi, 13 years ago
File size: 5.4 KB
Line 
1import os
2import time
3from shallow_water_domain_avalanche import *
4from Numeric import Float, sqrt
5from config import g, epsilon
6from numpy import sin, cos, tan, array, zero
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 = 1.0
103t0 = time.time()
104
105while finaltime <= 1.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([-500.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 = "B-case2-"
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.