source: anuga_work/production/sydney_2006/smf_approx.py @ 4347

Last change on this file since 4347 was 4316, checked in by sexton, 18 years ago

few more updates for slide report

File size: 8.7 KB
Line 
1"""Script for running a tsunami inundation scenario for Sydney, NSW, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project.py
5The output sww file is stored in project.outputdir
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and a simulated submarine landslide.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Adrian Hitchman and Jane Sexton, GA - 2006
11"""
12
13#-------------------------------------------------------------------------------# Import necessary modules
14#-------------------------------------------------------------------------------
15
16# Standard modules
17import os
18import time
19
20from pylab import *
21ion()
22hold(True)
23
24# Application specific imports
25import project                 # Definition of file names and polygons
26from anuga.shallow_water.smf import slide_tsunami, slump_tsunami  # Function for submarine mudslide
27
28def stats_slide(depth,length,width,thickness,slope,gamma):
29    from math import sin, tan, radians, pi, sqrt, exp
30
31    gravity = 9.8
32    massco = 1.0
33    dragco = 1.0
34    psi = 0.0
35    sint = sin(radians(slope))
36    tant = tan(radians(slope))
37    tanp = tan(radians(psi)) 
38   
39    a0 = gravity * sint * ((gamma-1)/(gamma+massco)) * (1-(tanp/tant))
40    ut = sqrt((gravity*depth) * (length*sint/depth) \
41                    * (pi*(gamma-1)/(2*dragco)) * (1-(tanp/tant)))
42    s0 = ut**2 / a0
43    t0 = ut / a0
44
45    #calculate some parameters of the water displacement produced by the slide
46
47    w = t0 * sqrt(gravity*depth)
48    a2D = s0 * (0.0574 - (0.0431*sint)) \
49             * thickness/length \
50             * ((length*sint/depth)**1.25) \
51             * (1 - exp(-2.2*(gamma-1)))
52    a3D = a2D / (1 + (15.5*sqrt(depth/(length*sint))))
53   
54    #scaled_amp = a3D/a2D
55   
56    return a3D, w
57
58def stats_slump(depth,length,width,thickness,slope,gamma):
59    from math import sin, radians, sqrt
60
61    radius = length**2 / (8.0 * thickness)   
62    massco = 1.0
63    dphi = 0.48
64    gravity = 9.8
65    sint = sin(radians(slope))
66
67    s0 = radius * dphi / 2
68    t0 = sqrt((radius*(gamma+massco)) / (gravity*(gamma-1)))
69    a0 = s0 / t0**2
70    um = s0 / t0
71
72    #calculate some parameters of the water displacement produced by the slump
73
74    w = t0 * sqrt(gravity*depth)
75    a2D = s0 * (0.131/sint) \
76             * thickness/length \
77             * (length*sint/depth)**1.25 \
78             * (length/radius)**0.63 * dphi**0.39 \
79             * (1.47 - (0.35*(gamma-1))) * (gamma-1)
80    a3D = a2D / (1. + (2.06*sqrt(depth/length)))
81
82    #scaled_amp = a3D/a2D
83   
84    return a3D, w
85
86which_one = 'slide'
87#which_one = 'slump'
88#par = 'depth'
89#par = 'density'
90par = 'depth'
91
92depths = range(750,2500,50)
93allgammas = arange(1.45,1.85,.05)
94all_lengths=[7000,17000,1000]
95all_slopes = arange(2,10,.2)
96lengths = [16840,13500,4189,9903]
97widths = [8860,4350,2898,4150]
98thicknesses = [424,165,53,140]
99slopes = [4,4,2.3,3.7]
100gammas = [1.46,1.49,1.48,1.48]
101smftype = ['Bulli','Shovel','Yacaaba','Birubi']
102#smfdepths = [2087,968,1119]
103smfdepths = [1470,877,1320,938]
104
105if par == 'depth':
106   
107    for i in range(len(lengths)):
108        amp = []
109        wavelength = []
110        for depth in depths:
111           
112            if which_one == 'slide':
113                a, w = stats_slide(depth,lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
114            else:
115                a, w = stats_slump(depth,lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
116
117            amp.append(a)
118            wavelength.append(w/1000.)
119
120        figure(1)
121        plot(depths,amp)
122        xlabel('Depth (m)')
123        ylabel('amplitude (m)')
124
125        figure(2)
126        plot(depths,wavelength)
127        xlabel('Depth (m)')
128        ylabel('wavelength (km)')
129
130    figure(1)
131    legend(smftype)
132    savefig('depth_vs_amp')
133   
134    figure(2)
135    legend(smftype)
136    savefig('depth_vs_wavelength')
137
138    adata = []
139    wdata = []
140    for i in range(len(smfdepths)):
141        if which_one == 'slide':
142            a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
143        else:
144            a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
145        adata.append(a)
146        wdata.append(w/1000.)
147
148    figure(1)
149    plot(smfdepths,adata,'^')
150    axis([750,2500,-0.5,6])
151    savefig('depth_vs_amp_data%s'%which_one)
152    figure(2)
153    plot(smfdepths,wdata,'^')
154    savefig('depth_vs_wavelength_data%s'%which_one)
155
156if par == 'density':
157   
158    for i in range(len(lengths)):
159        amp = []
160        wavelength = []
161        for gamma in allgammas:
162           
163            if which_one == 'slide':
164                a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gamma)
165            else:
166                a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gamma)
167
168            amp.append(a)
169            wavelength.append(w/1000.)
170
171        figure(1)
172        plot(allgammas,amp)
173        xlabel('Density (g/cm3)')
174        ylabel('amplitude (m)')
175
176        figure(2)
177        plot(allgammas,wavelength)
178        xlabel('Density (g/cm3)')
179        ylabel('wavelength (km)')
180
181    figure(1)
182    legend(smftype)
183    savefig('density_vs_amp')
184    figure(2)
185    legend(smftype)
186    savefig('density_vs_wavelength')
187
188    adata = []
189    wdata = []
190    for i in range(len(gammas)):
191        if which_one == 'slide':
192            a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
193        else:
194            a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
195        adata.append(a)
196        wdata.append(w/1000.)
197
198    figure(1)
199    plot(gammas,adata,'^')
200    savefig('density_vs_amp_data%s'%which_one)
201    figure(2)
202    plot(gammas,wdata,'^')
203    savefig('density_vs wavelength_data%s'%which_one)
204
205if par == 'slope':
206   
207    for i in range(len(lengths)):
208        amp = []
209        wavelength = []
210        for slope in all_slopes:
211           
212            if which_one == 'slide':
213                a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slope,gammas[i])
214            else:
215                a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slope,gammas[i])
216
217            amp.append(a)
218            wavelength.append(w/1000.)
219
220        figure(1)
221        plot(all_slopes,amp)
222        xlabel('Slope (deg)')
223        ylabel('amplitude (m)')
224
225        figure(2)
226        plot(all_slopes,wavelength)
227        xlabel('Slope (deg)')
228        ylabel('wavelength (km)')
229
230    figure(1)
231    legend(smftype)
232    savefig('slope_vs_amp')
233    figure(2)
234    legend(smftype)
235    savefig('slope_vs_wavelength')
236
237    adata = []
238    wdata = []
239    for i in range(len(lengths)):
240        if which_one == 'slide':
241            a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
242        else:
243            a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
244        adata.append(a)
245        wdata.append(w/1000.)
246
247    figure(1)
248    plot(slopes,adata,'^')
249    savefig('slope_vs_amp_data%s'%which_one)
250    figure(2)
251    plot(slopes,wdata,'^')
252    savefig('slope_vs_wavelength_data%s'%which_one)
253
254
255close('all')
256
257##import sys; sys.exit()
258##tsunami_source = slide_tsunami(length=16840.0,
259##                               width=8860.0,
260##                               depth=2087.0,
261##                               slope=4.0,
262##                               thickness=424.0,
263##                               gamma=1.46,
264##                               #dphi=0.23,
265##                               x0=project.slump_origin2[0],
266##                               y0=project.slump_origin2[1],
267##                               #alpha=0.0,
268##                               #domain=domain,
269##                               verbose=True)
270##        tsunami_source = slump_tsunami(length=16840.0,
271##                                       width=8860.0,
272##                                       #depth=2087.0,
273##                                       slope=4.0,
274##                                       thickness=424.0,
275##                                       gamma=1.45,
276##                                       #dphi=0.23,
277##                                       #radius = 3330,
278##                                       x0=project.slump_origin2[0],
279##                                       y0=project.slump_origin2[1],
280##                                       #alpha=0.0,
281##                                       #domain=domain,
282##                                       verbose=True)
Note: See TracBrowser for help on using the repository browser.