source: inundation/pyvolution/smf.py @ 2380

Last change on this file since 2380 was 2219, checked in by duncan, 18 years ago

add verbose setting

File size: 11.9 KB
Line 
1#
2# slide_tsunami function
3#
4
5"""This function returns a callable object representing an initial water
6   displacement generated by a submarine sediment slide.
7
8Using input parameters:
9
10Required
11 length  downslope slide length
12 depth   water depth to slide centre of mass
13 slope   bathymetric slope
14
15Optional
16 x0      x origin (0)
17 y0      y origin (0)
18 alpha   angular orientation of slide in xy plane (0)
19 w       slide width (0.25*length)
20 T       slide thickness (0.01*length)
21 g       acceleration due to gravity (9.8)
22 gamma   specific density of sediments (1.85)
23 Cm      added mass coefficient (1)
24 Cd      drag coefficient (1)
25 Cn      friction coefficient (0)
26 psi     (0)
27 dx      offset of second Gaussian (0.2*width of first Gaussian)
28 kappa   multiplier for sech^2 function (3.0)
29 kappad  multiplier for second Gaussian function (0.8)
30 zsmall  an amount near to zero (0.01)
31
32The following parameters are calculated:
33
34 a0      initial acceleration
35 ut      theoretical terminal velocity
36 s0      charactistic distance of motion
37 t0      characteristic time of motion
38 w       initial wavelength of tsunami
39 a2D     2D initial amplitude of tsunami
40 a3D     3D initial amplitude of tsunami
41
42The returned object is a callable double Gaussian function that represents
43the initial water displacement generated by a submarine sediment slide.
44
45Adrian Hitchman
46Geoscience Australia, June 2005
47"""
48
49def slide_tsunami(length, depth, slope, width=None, thickness=None, \
50                  x0=0.0, y0=0.0, alpha=0.0, \
51                  gravity=9.8, gamma=1.85, \
52                  massco=1, dragco=1, frictionco=0, psi=0, \
53                  dx=None, kappa=3.0, kappad=0.8, zsmall=0.01,
54                  verbose=False):
55
56    from math import sin, tan, radians, pi, sqrt, exp
57
58    #if width not provided, set to typical value
59    if width is None:
60        width = 0.25 * length
61
62    #if thickness not provided, set to typical value
63    if thickness is None:
64        thickness = 0.01 * length
65
66    #calculate some parameters of the slide
67
68    sint = sin(radians(slope))
69    tant = tan(radians(slope))
70    tanp = tan(radians(psi))
71
72    a0 = gravity * sint * ((gamma-1)/(gamma+massco)) * (1-(tanp/tant))
73    ut = sqrt((gravity*depth) * (length*sint/depth) \
74                    * (pi*(gamma-1)/(2*dragco)) * (1-(tanp/tant)))
75    s0 = ut**2 / a0
76    t0 = ut / a0
77
78    #calculate some parameters of the water displacement produced by the slide
79
80    w = t0 * sqrt(gravity*depth)
81    a2D = s0 * (0.0574 - (0.0431*sint)) \
82             * (thickness/length) \
83             * ((length*sint/depth)**1.25) \
84             * (1 - exp(-2.2*(gamma-1)))
85    a3D = a2D / (1 + (15.5*sqrt(depth/(length*sint))))
86
87    #a few temporary print statements
88    if verbose is True:
89        print '\nThe slide ...'
90        print '\tLength: ', length
91        print '\tDepth: ', depth
92        print '\tSlope: ', slope
93        print '\tWidth: ', width
94        print '\tThickness: ', thickness
95        print '\tx0: ', x0
96        print '\ty0: ', y0
97        print '\tAlpha: ', alpha
98        print '\tAcceleration: ', a0
99        print '\tTerminal velocity: ', ut
100        print '\tChar time: ', t0
101        print '\tChar distance: ', s0
102        print '\nThe tsunami ...'
103        print '\tWavelength: ', w
104        print '\t2D amplitude: ', a2D
105        print '\t3D amplitude: ', a3D
106
107    #keep an eye on some of the assumptions built into the maths
108
109    if ((slope < 5) or (slope > 30)):
110        if verbose is True:
111            print 'WARNING: slope out of range (5 - 30 degrees) ', slope
112    if ((depth/length < 0.06) or (depth/length > 1.5)):
113        if verbose is True:
114            print  'WARNING: d/b out of range (0.06 - 1.5) ', depth/length
115    if ((thickness/length < 0.008) or (thickness/length > 0.2)):
116        if verbose is True:
117            print 'WARNING: T/b out of range (0.008 - 0.2) ', thickness/length
118    if ((gamma < 1.46) or (gamma > 2.93)):
119        if verbose is True:
120            print 'WARNING: gamma out of range (1.46 - 2.93) ', gamma
121
122    return Double_gaussian(a3D=a3D, wavelength=w, width=width, \
123                           x0=x0, y0=y0, alpha=alpha, \
124                           dx=dx, kappa=kappa, kappad=kappad, zsmall=zsmall)
125
126#
127# slump_tsunami function
128#
129
130"""This function returns a callable object representing an initial water
131   displacement generated by a submarine sediment slump.
132
133Using input parameters:
134
135Required
136 length  downslope slump length
137 depth   water depth to slump centre of mass
138 slope   bathymetric slope
139
140Optional
141 x0      x origin (0)
142 y0      y origin (0)
143 alpha   angular orientation of slide in xy plane (0)
144 w       slump width (1.0*length)
145 T       slump thickness (0.1*length)
146 R       slump radius of curvature (b^2/(8*T))
147 del_phi slump angular displacement (0.48)
148 g       acceleration due to gravity (9.8)
149 gamma   specific density of sediments (1.85)
150 Cm      added mass coefficient (1)
151 Cd      drag coefficient (1)
152 Cn      friction coefficient (0)
153 dx      offset of second Gaussian (0.2*width of first Gaussian)
154 kappa   multiplier for sech^2 function (3.0)
155 kappad  multiplier for second Gaussian function (0.8)
156 zsmall  an amount near to zero (0.01)
157
158The following parameters are calculated:
159
160 a0      initial acceleration
161 um      maximum velocity
162 s0      charactistic distance of motion
163 t0      characteristic time of motion
164 w       initial wavelength of tsunami
165 a2D     2D initial amplitude of tsunami
166 a3D     3D initial amplitude of tsunami
167
168The returned object is a callable double Gaussian function that represents
169the initial water displacement generated by a submarine sediment slump.
170
171Adrian Hitchman
172Geoscience Australia, June 2005
173"""
174
175def slump_tsunami(length, depth, slope, width=None, thickness=None, \
176                  radius=None, dphi=0.48, x0=0.0, y0=0.0, alpha=0.0, \
177                  gravity=9.8, gamma=1.85, \
178                  massco=1, dragco=1, frictionco=0, \
179                  dx=None, kappa=3.0, kappad=0.8, zsmall=0.01,
180                  verbose=False):
181
182    from math import sin, radians, sqrt
183
184    #if width not provided, set to typical value
185    if width is None:
186        width = length
187
188    #if thickness not provided, set to typical value
189    if thickness is None:
190        thickness = 0.1 * length
191
192    #if radius not provided, set to typical value
193    if radius is None:
194        radius = length**2 / (8.0 * thickness)
195
196    #calculate some parameters of the slump
197
198    sint = sin(radians(slope))
199
200    s0 = radius * dphi / 2
201    t0 = sqrt((radius*(gamma+massco)) / (gravity*(gamma-1)))
202    a0 = s0 / t0**2
203    um = s0 / t0
204
205    #calculate some parameters of the water displacement produced by the slump
206
207    w = t0 * sqrt(gravity*depth)
208    a2D = s0 * (0.131/sint) \
209             * (thickness/length) \
210             * (length*sint/depth)**1.25 \
211             * (length/radius)**0.63 * dphi**0.39 \
212             * (1.47 - (0.35*(gamma-1))) * (gamma-1)
213    a3D = a2D / (1 + (2.06*sqrt(depth/length)))
214
215    #a few temporary print statements
216    if verbose is True:
217        print '\nThe slump ...'
218        print '\tLength: ', length
219        print '\tDepth: ', depth
220        print '\tSlope: ', slope
221        print '\tWidth: ', width
222        print '\tThickness: ', thickness
223        print '\tRadius: ', radius
224        print '\tDphi: ', dphi
225        print '\tx0: ', x0
226        print '\ty0: ', y0
227        print '\tAlpha: ', alpha
228        print '\tAcceleration: ', a0
229        print '\tMaximum velocity: ', um
230        print '\tChar time: ', t0
231        print '\tChar distance: ', s0
232        print '\nThe tsunami ...'
233        print '\tWavelength: ', w
234        print '\t2D amplitude: ', a2D
235        print '\t3D amplitude: ', a3D
236
237    #keep an eye on some of the assumptions built into the maths
238
239    if ((slope < 10) or (slope > 30)):       
240        if verbose is True:
241            print 'WARNING: slope out of range (10 - 30 degrees) ', slope
242    if ((depth/length < 0.34) or (depth/length > 0.5)):     
243        if verbose is True:
244            print  'WARNING: d/b out of range (0.34 - 0.5) ', depth/length
245    if ((thickness/length < 0.10) or (thickness/length > 0.15)):     
246        if verbose is True:
247            print 'WARNING: T/b out of range (0.10 - 0.15) ', thickness/length
248    if ((radius/length < 1.0) or (radius/length > 2.0)):     
249        if verbose is True:
250            print 'WARNING: R/b out of range (1 - 2) ', radius/length
251    if ((dphi < 0.10) or (dphi > 0.52)):     
252        if verbose is True:
253            print 'WARNING: del_phi out of range (0.10 - 0.52) ', dphi
254    if ((gamma < 1.46) or (gamma > 2.93)):     
255        if verbose is True:
256            print 'WARNING: gamma out of range (1.46 - 2.93) ', gamma
257
258    return Double_gaussian(a3D=a3D, wavelength=w, width=width, \
259                           x0=x0, y0=y0, alpha=alpha, \
260                           dx=dx, kappa=kappa, kappad=kappad, zsmall=zsmall)
261
262#
263# Double_gaussian class
264#
265
266"""This is a callable class representing the initial water displacment
267   generated by a sediment slide or slump.
268
269Using input parameters:
270
271Required
272 w       initial wavelength of tsunami
273 a3D     3D initial amplitude of tsunami
274 width   width of smf
275
276Optional
277 x0      x origin of smf
278 y0      y origin of smf
279 alpha   angular orientation of smf in xy plane (0)
280 dx      offset of second Gaussian (0.2*width of first Gaussian)
281 kappa   multiplier for sech^2 function (3.0)
282 kappad  multiplier for second Gaussian function (0.8)
283 zsmall  an amount near to zero (0.01)
284
285Adrian Hitchman
286Geoscience Australia, June 2005
287"""
288
289class Double_gaussian:
290
291    def __init__(self, a3D, wavelength, width, x0, y0, alpha, \
292                 dx, kappa, kappad, zsmall):
293        self.a3D = a3D
294        self.wavelength = wavelength
295        self.width = width
296        self.x0 = x0
297        self.y0 = y0
298        self.alpha = alpha
299        self.kappa = kappa
300        self.kappad = kappad
301
302        if dx is None:
303            self.determineDX(zsmall=zsmall)
304        else:
305            self.dx = dx
306
307    def __call__(self, x, y):
308        """Make Double_gaussian a callable object.
309
310        If called as a function, this object returns z values representing
311        the initial 3D distribution of water heights at the points (x,y)
312        produced by a submarine mass failure.
313        """
314
315        from math import sin, cos, radians, exp, cosh
316        from Numeric import zeros, Float
317
318        #ensure vectors x and y have the same length
319        N = len(x)
320        assert N == len(y)
321
322        am = self.a3D
323        wa = self.wavelength
324        wi = self.width
325        x0 = self.x0
326        y0 = self.y0
327        alpha = self.alpha
328        dx = self.dx
329        kappa = self.kappa
330        kappad = self.kappad
331
332        #double Gaussian calculation assumes water displacement is oriented
333        #E-W, so, for displacement at some angle alpha clockwise from the E-W
334        #direction, rotate (x,y) coordinates anti-clockwise by alpha
335
336        cosa = cos(radians(alpha))
337        sina = sin(radians(alpha))
338
339        xr = ((x-x0) * cosa - (y-y0) * sina) + x0
340        yr = ((x-x0) * sina + (y-y0) * cosa) + y0
341
342        z = zeros(N, Float)
343
344        for i in range(N):
345            try:
346                z[i] =  -am / ((cosh(kappa*(yr[i]-y0)/(wi+wa)))**2) \
347                            * (exp(-((xr[i]-x0)/wa)**2) - \
348                                kappad*exp(-((xr[i]-dx-x0)/wa)**2))
349            except OverflowError:
350                pass
351
352        return z
353
354    def determineDX(self, zsmall):
355        """Determine a suitable offset for the second Gaussian function.
356
357        A suitable offset for the second Gaussian function is taken to
358        be some fraction of the 'width' of the first Gaussian function.
359
360        The 'width' of the first Gaussian is obtained from the range of
361        the x coordinates over which the function takes values from
362        'near zero', through 1, and back to 'near zero'.
363
364        The parameter zsmall passed to this function specifies how much
365        'near zero' is.
366        """
367
368        from math import sqrt, log, e
369
370        a = self.a3D
371        c = self.wavelength
372
373        self.dx = 2.0 * (c * sqrt(-log((zsmall/a),e))) / 5.0
Note: See TracBrowser for help on using the repository browser.