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

Last change on this file since 5589 was 4504, checked in by sexton, 17 years ago

minor change to inundation_damage + updates to damage scripts for all scenarios

File size: 13.5 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=1.0, 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        print '\tDelta x ', dx
251        print '\tsmall ', zsmall
252        print '\tKappa d ', kappad
253
254    #keep an eye on some of the assumptions built into the maths
255
256    if ((slope < 10) or (slope > 30)):       
257        if verbose is True:
258            print 'WARNING: slope out of range (10 - 30 degrees) ', slope
259    if ((depth/length < 0.34) or (depth/length > 0.5)):     
260        if verbose is True:
261            print  'WARNING: d/b out of range (0.34 - 0.5) ', depth/length
262    if ((thickness/length < 0.10) or (thickness/length > 0.15)):     
263        if verbose is True:
264            print 'WARNING: T/b out of range (0.10 - 0.15) ', thickness/length
265    if ((radius/length < 1.0) or (radius/length > 2.0)):     
266        if verbose is True:
267            print 'WARNING: R/b out of range (1 - 2) ', radius/length
268    if ((dphi < 0.10) or (dphi > 0.52)):     
269        if verbose is True:
270            print 'WARNING: del_phi out of range (0.10 - 0.52) ', dphi
271    if ((gamma < 1.46) or (gamma > 2.93)):     
272        if verbose is True:
273            print 'WARNING: gamma out of range (1.46 - 2.93) ', gamma
274
275    return Double_gaussian(a3D=a3D, wavelength=w, width=width, \
276                           x0=x0, y0=y0, alpha=alpha, \
277                           dx=dx, kappa=kappa, kappad=kappad, zsmall=zsmall)
278
279#
280# Double_gaussian class
281#
282
283"""This is a callable class representing the initial water displacment
284   generated by a sediment slide or slump.
285
286Using input parameters:
287
288Required
289 w       initial wavelength of tsunami
290 a3D     3D initial amplitude of tsunami
291 width   width of smf
292
293Optional
294 x0      x origin of smf
295 y0      y origin of smf
296 alpha   angular orientation of smf in xy plane (0)
297 dx      offset of second Gaussian (0.2*width of first Gaussian)
298 kappa   multiplier for sech^2 function (3.0)
299 kappad  multiplier for second Gaussian function (0.8)
300 zsmall  an amount near to zero (0.01)
301
302Adrian Hitchman
303Geoscience Australia, June 2005
304"""
305
306class Double_gaussian:
307
308    def __init__(self, a3D, wavelength, width, x0, y0, alpha, \
309                 dx, kappa, kappad, zsmall):
310        self.a3D = a3D
311        self.wavelength = wavelength
312        self.width = width
313        self.x0 = x0
314        self.y0 = y0
315        self.alpha = alpha
316        self.kappa = kappa
317        self.kappad = kappad
318
319        if dx is None:
320            self.determineDX(zsmall=zsmall)
321        else:
322            self.dx = dx
323
324    def __call__(self, x, y):
325        """Make Double_gaussian a callable object.
326
327        If called as a function, this object returns z values representing
328        the initial 3D distribution of water heights at the points (x,y)
329        produced by a submarine mass failure.
330        """
331
332        from math import sin, cos, radians, exp, cosh
333        from Numeric import zeros, Float
334
335        #ensure vectors x and y have the same length
336        N = len(x)
337        assert N == len(y)
338
339        am = self.a3D
340        am2 = 1.0
341        #am2 = self.a2D
342        wa = self.wavelength
343        wi = self.width
344        x0 = self.x0
345        y0 = self.y0
346        alpha = self.alpha
347        dx = self.dx
348        kappa = self.kappa
349        kappad = self.kappad
350        #amin = self.find_min(x0,wa,kappad,kappa,dx,am)
351        amin = 1.0
352
353        #double Gaussian calculation assumes water displacement is oriented
354        #E-W, so, for displacement at some angle alpha clockwise from the E-W
355        #direction, rotate (x,y) coordinates anti-clockwise by alpha
356
357        cosa = cos(radians(alpha))
358        sina = sin(radians(alpha))
359
360        xr = ((x-x0) * cosa - (y-y0) * sina) + x0
361        yr = ((x-x0) * sina + (y-y0) * cosa) + y0
362
363        z = zeros(N, Float)
364        maxz = 0.0
365        minz = 0.0
366        for i in range(N):
367            try:
368                if i == 10: print 'hello', x[i], x0, xr[i]
369                z[i] =  -(am/am2) / (amin*(cosh(kappa*(yr[i]-y0)/(wi+wa)))**2) \
370                            * (exp(-((xr[i]-x0)/wa)**2) - \
371                                kappad*exp(-((xr[i]-dx-x0)/wa)**2))
372                if z[i] > maxz: maxz = z[i]
373                if z[i] < minz: minz = z[i]
374               
375            except OverflowError:
376                pass
377
378        print 'max z', maxz
379        print 'min z', minz
380               
381        return z
382
383    def determineDX(self, zsmall):
384        """Determine a suitable offset for the second Gaussian function.
385
386        A suitable offset for the second Gaussian function is taken to
387        be some fraction of the 'width' of the first Gaussian function.
388
389        The 'width' of the first Gaussian is obtained from the range of
390        the x coordinates over which the function takes values from
391        'near zero', through 1, and back to 'near zero'.
392
393        The parameter zsmall passed to this function specifies how much
394        'near zero' is.
395        """
396
397        from math import sqrt, log, e
398
399        a = self.a3D
400        c = self.wavelength
401
402        self.dx = 2.0 * (c * sqrt(-log((zsmall/a),e))) / 5.0
403
404    def find_min(self, x0, wa, kappad, kappa, dx, am):
405        """Determine eta_min to scale eta(x,y)
406        """
407
408        from math import exp, cosh
409
410        step = 10.0
411        x = x0
412        deriv = 10.0
413        tol = 0.001
414        count_max = 1000000
415        c = 0
416        am = 1.0
417
418        while c < count_max and deriv > 0:
419            deriv = (x-x0)*exp(-((x-x0)/wa)**2.0) - \
420                    kappad*(x-dx-x0)*exp(-((x-dx-x0)/wa)**2.0)
421            if deriv < 0: xstar = x
422            x -= step
423            c += 1
424       
425        etastar = -am * (exp(-((xstar-x0)/wa)**2) - \
426                   kappad*exp(-((xstar-dx-x0)/wa)**2))
427
428        return etastar
Note: See TracBrowser for help on using the repository browser.