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

Last change on this file since 7632 was 7632, checked in by habili, 14 years ago

Fixed bug with logger

File size: 15.0 KB
RevLine 
[5675]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"""
[6157]48
[7276]49import numpy as num
[6157]50
[7632]51import anuga.utilities.log as lg
[6157]52
[7317]53
[5675]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:
[7632]147        lg.critical('\nThe slide ...')
148        lg.critical('\tLength: %s' % str(length))
149        lg.critical('\tDepth: %s' % str(depth))
150        lg.critical('\tSlope: %s' % str(slope))
151        lg.critical('\tWidth: %s' % str(width))
152        lg.critical('\tThickness: %s' % str(thickness))
153        lg.critical('\tx0: %s' % str(x0))
154        lg.critical('\ty0: %s' % str(y0))
155        lg.critical('\tAlpha: %s' % str(alpha))
156        lg.critical('\tAcceleration: %s' % str(a0))
157        lg.critical('\tTerminal velocity: %s' % str(ut))
158        lg.critical('\tChar time: %s' % str(t0))
159        lg.critical('\tChar distance: %s' % str(s0))
160        lg.critical('\nThe tsunami ...')
161        lg.critical('\tWavelength: %s' % str(w))
162        lg.critical('\t2D amplitude: %s' % str(a2D))
163        lg.critical('\t3D amplitude: %s' % str(a3D))
164        lg.critical('\tscale for eta(x,y): %s' % str(scale))
[5675]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:
[7632]170            lg.critical('WARNING: slope out of range (5 - 30 degrees) %s'
[7317]171                         % str(slope))
[5675]172    if ((depth/length < 0.06) or (depth/length > 1.5)):
173        if verbose is True:
[7632]174            lg.critical('WARNING: d/b out of range (0.06 - 1.5) %s'
[7317]175                         % str(depth/length))
[5675]176    if ((thickness/length < 0.008) or (thickness/length > 0.2)):
177        if verbose is True:
[7632]178            lg.critical('WARNING: T/b out of range (0.008 - 0.2) %s'
[7317]179                         % str(thickness/length))
[5675]180    if ((gamma < 1.46) or (gamma > 2.93)):
181        if verbose is True:
[7632]182            lg.critical('WARNING: gamma out of range (1.46 - 2.93) %s'
[7317]183                         % str(gamma))
[5675]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:
[7632]294        lg.critical('\nThe slump ...')
295        lg.critical('\tLength: %s' % str(length))
296        lg.critical('\tDepth: %s' % str(depth))
297        lg.critical('\tSlope: %s' % str(slope))
298        lg.critical('\tWidth: %s' % str(width))
299        lg.critical('\tThickness: %s' % str(thickness))
300        lg.critical('\tRadius: %s' % str(radius))
301        lg.critical('\tDphi: %s' % str(dphi))
302        lg.critical('\tx0: %s' % str(x0))
303        lg.critical('\ty0: %s' % str(y0))
304        lg.critical('\tAlpha: %s' % str(alpha))
305        lg.critical('\tAcceleration: %s' % str(a0))
306        lg.critical('\tMaximum velocity: %s' % str(um))
307        lg.critical('\tChar time: %s' % str(t0))
308        lg.critical('\tChar distance: %s' % str(s0))
309        lg.critical('\nThe tsunami ...')
310        lg.critical('\tWavelength: %s' % str(w))
311        lg.critical('\t2D amplitude: %s' % str(a2D))
312        lg.critical('\t3D amplitude: %s' % str(a3D))
313        lg.critical('\tDelta x %s' % str(dx))
314        lg.critical('\tsmall %s' % str(zsmall))
315        lg.critical('\tKappa d  %s' % str(kappad))
316        lg.critical('\tscale for eta(x,y): %s' % str(scale))
[5675]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:
[7632]322            lg.critical('WARNING: slope out of range (10 - 30 degrees) %s'
[7317]323                         % str(slope))
[5675]324    if ((depth/length < 0.34) or (depth/length > 0.5)):     
325        if verbose is True:
[7632]326            lg.critical('WARNING: d/b out of range (0.34 - 0.5) %s'
[7317]327                         % str(depth/length))
[5675]328    if ((thickness/length < 0.10) or (thickness/length > 0.15)):     
329        if verbose is True:
[7632]330            lg.critical('WARNING: T/b out of range (0.10 - 0.15) %s'
[7317]331                         % str(thickness/length))
[5675]332    if ((radius/length < 1.0) or (radius/length > 2.0)):     
333        if verbose is True:
[7632]334            lg.critical('WARNING: R/b out of range (1 - 2) %s'
[7317]335                         % str(radius/length))
[5675]336    if ((dphi < 0.10) or (dphi > 0.52)):     
337        if verbose is True:
[7632]338            lg.critical('WARNING: del_phi out of range (0.10 - 0.52) %s'
[7317]339                         % str(dphi))
[5675]340    if ((gamma < 1.46) or (gamma > 2.93)):     
341        if verbose is True:
[7632]342            lg.critical('WARNING: gamma out of range (1.46 - 2.93) %s'
[7317]343                         % str(gamma))
[5675]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
[7276]429        z = num.zeros(N, num.float)
[5675]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.