source: anuga_core/source/anuga/shallow_water/smf.py @ 6489

Last change on this file since 6489 was 6157, checked in by rwilson, 15 years ago

Change Numeric imports to general form - ready to change to NumPy?.

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