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

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

Replaced 'print' statements with log.critical() calls.

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