source: inundation/pyvolution/smf.py @ 2689

Last change on this file since 2689 was 2393, checked in by sexton, 18 years ago

Allowing domain to be passed in slump/slide function to allow georeferencing.

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