Changeset 5675


Ignore:
Timestamp:
Aug 22, 2008, 1:42:41 PM (16 years ago)
Author:
sexton
Message:

updates to smf to incorporate a scaling parameter

Location:
anuga_core/source/anuga/shallow_water
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/shallow_water/smf.py

    r5648 r5675  
    4646Geoscience Australia, June 2005
    4747"""
     48def find_min(x0, wa, kappad, dx):
     49    """Determine eta_min to scale eta(x,y)
     50
     51       eta(x,y) = n03d/nmin*[-f(y)]*g(x)
     52
     53       nmin = min (-f(y)*g(x) )
     54            = -f(ystar)*g(xstar)
     55
     56       ystar = min (-f(y) ), i.e. diff(-f(y))=0
     57       xstar = min ( g(x) ), i.e. diff(g(x))=0
     58
     59       ystar = y0 and -f(ystar)=1.0
     60    """
     61    from math import exp, cosh
     62
     63    step = 0.05
     64    x = x0+50.
     65    deriv = 10.0
     66    count_max = 1000000
     67    c = 0
     68    deriv = 10.
     69    f_ystar = 1.
     70   
     71    while c < count_max and deriv > 0:
     72        deriv = (x-x0)*exp(-((x-x0)/wa)**2.0) - \
     73                kappad*(x-dx-x0)*exp(-((x-dx-x0)/wa)**2.0)
     74       
     75        if deriv <= 0: xstar = x
     76        x -= step
     77        c += 1
     78   
     79    g_xstar = exp(-((xstar-x0)/wa)**2)-kappad*exp(-((xstar-dx-x0)/wa)**2)
     80
     81    etastar = g_xstar*f_ystar
     82
     83    return etastar
    4884
    4985def slide_tsunami(length, depth, slope, width=None, thickness=None, \
    5086                  x0=0.0, y0=0.0, alpha=0.0, \
    51                   gravity=9.8, gamma=None, \
    52                   massco=None, dragco=None, frictionco=0, psi=0, \
    53                   dx=None, kappa=None, kappad=None, zsmall=0.01, \
    54                   scale=None, domain=None, verbose=False):
     87                  gravity=9.8, gamma=1.85, \
     88                  massco=1, dragco=1, frictionco=0, psi=0, \
     89                  dx=None, kappa=3.0, kappad=0.8, zsmall=0.01, \
     90                  scale=None,
     91                  domain=None, verbose=False):
    5592
    5693    from math import sin, tan, radians, pi, sqrt, exp
    … …  
    70107        thickness = 0.01 * length
    71108
    72     if gamma is None:
    73         gamma = 1.85
    74 
    75     if massco is None:
    76         massco = 1.
    77 
    78     if dragco is None:
    79         dragco = 1.
    80 
    81     if kappa is None:
    82         kappa = 3.
    83 
    84     if kappad is None:
    85         kappa = 0.8
    86 
    87     if dx is None:
    88         dx = 0.
    89 
    90109    #calculate some parameters of the slide
    91110
    … …  
    108127             * (1 - exp(-2.2*(gamma-1)))
    109128    a3D = a2D / (1 + (15.5*sqrt(depth/(length*sint))))
     129
     130    from math import sqrt, log, e
     131    dx = 2.0 * (w * sqrt(-log((zsmall/a3D),e))) / 5.0
     132       
     133    # determine nmin for scaling of eta(x,y)
     134    nmin = find_min(x0,w,kappad,dx) 
     135   
     136    if scale is None:
     137        scale = a3D/nmin
    110138       
    111139    #a few temporary print statements
    … …  
    128156        print '\t2D amplitude: ', a2D
    129157        print '\t3D amplitude: ', a3D
     158        print '\tscale for eta(x,y):', scale
    130159
    131160    #keep an eye on some of the assumptions built into the maths
    … …  
    144173            print 'WARNING: gamma out of range (1.46 - 2.93) ', gamma
    145174
    146     return Double_gaussian(a3D=a3D, wavelength=w, width=width, \
    147                            x0=x0, y0=y0, alpha=alpha, \
    148                            dx=dx, kappa=kappa, kappad=kappad, zsmall=zsmall, scale=scale)
     175    return Double_gaussian(a3D, w, width, x0, y0, alpha, kappa, kappad, zsmall, dx, scale)
    149176
    150177#
    … …  
    199226def slump_tsunami(length, depth, slope, width=None, thickness=None, \
    200227                  radius=None, dphi=0.48, x0=0.0, y0=0.0, alpha=0.0, \
    201                   gravity=9.8, gamma=None, \
     228                  gravity=9.8, gamma=1.85, \
    202229                  massco=1, dragco=1, frictionco=0, \
    203                   dx=None, kappa=3.0, kappad=1.0, zsmall=0.01, \
     230                  dx=None, kappa=3.0, kappad=1.0, zsmall=0.01, scale=None, \
    204231                  domain=None,
    205232                  verbose=False):
    … …  
    244271    a3D = a2D / (1 + (2.06*sqrt(depth/length)))
    245272
     273    from math import sqrt, log, e
     274    dx = 2.0 * (w * sqrt(-log((zsmall/a3D),e))) / 5.0
     275       
     276    # determine nmin for scaling of eta(x,y)
     277    nmin = find_min(x0,w,kappad,dx) 
     278   
     279    if scale is None:
     280        scale = a3D/nmin
     281       
    246282    #a few temporary print statements
    247283    if verbose is True:
    … …  
    268304        print '\tsmall ', zsmall
    269305        print '\tKappa d ', kappad
     306        print '\tscale for eta(x,y):', scale
    270307
    271308    #keep an eye on some of the assumptions built into the maths
    … …  
    290327            print 'WARNING: gamma out of range (1.46 - 2.93) ', gamma
    291328
    292     return Double_gaussian(a3D=a3D, wavelength=w, width=width, \
    293                            x0=x0, y0=y0, alpha=alpha, \
    294                            dx=dx, kappa=kappa, kappad=kappad, zsmall=zsmall)
     329    return Double_gaussian(a3D, w, width, x0, y0, alpha, kappa, kappad, zsmall, dx, scale)
    295330
    296331#
    … …  
    324359
    325360    def __init__(self, a3D, wavelength, width, x0, y0, alpha, \
    326                  dx, kappa, kappad, zsmall, scale):
     361                 kappa, kappad, zsmall, dx, scale):
    327362        self.a3D = a3D
    328363        self.wavelength = wavelength
    … …  
    336371
    337372        if dx is None:
    338             self.determineDX(zsmall=zsmall)
    339         else:
    340             self.dx = dx
     373            from math import sqrt, log, e
     374            dx = 2.0 * (self.wavelength * sqrt(-log((zsmall/self.a3D),e))) / 5.0
     375        self.dx = dx
    341376
    342377    def __call__(self, x, y):
    … …  
    357392        am = self.a3D
    358393        am2 = 1.0
    359         #am2 = self.a2D
    360394        wa = self.wavelength
    361395        wi = self.width
    … …  
    367401        kappad = self.kappad
    368402        scale = self.scale
    369         #amin = self.find_min(x0,wa,kappad,kappa,dx,am)
    370         amin = 1.0
    371403
    372404        #double Gaussian calculation assumes water displacement is oriented
    … …  
    383415        maxz = 0.0
    384416        minz = 0.0
    385        
    386417        for i in range(N):
    387418            try:
    388                 #if i == 10: print 'hello', x[i], x0, xr[i]
    389419                z[i] =  -scale / ((cosh(kappa*(yr[i]-y0)/(wi+wa)))**2) \
    390420                            * (exp(-((xr[i]-x0)/wa)**2) - \
    … …  
    392422                if z[i] > maxz: maxz = z[i]
    393423                if z[i] < minz: minz = z[i]
    394                                
     424               
    395425            except OverflowError:
    396426                pass
    397 
    398         print 'max z', maxz
    399         print 'min z', minz
    400                        
     427               
    401428        return z
    402429
    … …  
    416443
    417444        from math import sqrt, log, e
    418 
     445       
    419446        a = self.a3D
    420447        c = self.wavelength
    … …  
    422449        self.dx = 2.0 * (c * sqrt(-log((zsmall/a),e))) / 5.0
    423450
    424     def find_min(self, x0, wa, kappad, kappa, dx, am):
    425         """Determine eta_min to scale eta(x,y)
    426         """
    427 
    428         from math import exp, cosh
    429 
    430         step = 10.0
    431         x = x0
    432         deriv = 10.0
    433         tol = 0.001
    434         count_max = 1000000
    435         c = 0
    436         am = 1.0
    437 
    438         while c < count_max and deriv > 0:
    439             deriv = (x-x0)*exp(-((x-x0)/wa)**2.0) - \
    440                     kappad*(x-dx-x0)*exp(-((x-dx-x0)/wa)**2.0)
    441             if deriv < 0: xstar = x
    442             x -= step
    443             c += 1
    444        
    445         etastar = -am * (exp(-((xstar-x0)/wa)**2) - \
    446                    kappad*exp(-((xstar-dx-x0)/wa)**2))
    447 
    448         return etastar
     451        return self.dx
     452

  • anuga_core/source/anuga/shallow_water/test_smf.py

    r5648 r5675  
    2424        scale = 1.0
    2525
    26         dg = Double_gaussian(a3D=a3D, wavelength=wavelength, width=width, \
    27                              x0=x0, y0=y0, alpha=alpha, dx = dx, \
    28                              kappa=kappa, kappad = kappad, zsmall = zsmall,
    29                              scale=scale)
     26        dg = Double_gaussian(a3D, wavelength, width, \
     27                             x0, y0, alpha, \
     28                             kappa, kappad, zsmall, dx, scale)
    3029
    3130        assert allclose(dg.a3D, a3D)
    … …  
    3938        assert allclose(dg.dx, dx)
    4039
    41 
     40   
    4241    def test_slide_tsunami(self):
    4342
    … …  
    4847        wid = 340.0
    4948        kappa = 3.0
     49        kappad = 0.8
     50        x0 = 100000.
    5051
    51         slide = slide_tsunami(length=len, depth=dep, slope=th, \
    52                               width = wid, thickness=thk, kappa=kappa)
     52        slide = slide_tsunami(length=len, depth=dep, slope=th, x0=x0, \
     53                              width = wid, thickness=thk, kappa=kappa, kappad=kappad, \
     54                              verbose=False)
    5355
    5456        assert allclose(slide.a3D, 0.07775819)
    5557        assert allclose(slide.wavelength, 2938.66695708)
    5658        assert allclose(slide.width, 340.0)
    57         assert allclose(slide.x0, 0.0)
    5859        assert allclose(slide.y0, 0.0)
    5960        assert allclose(slide.alpha, 0.0)
    60         assert allclose(slide.kappa, 3.0)
    61         assert allclose(slide.kappad, 0.8)
    6261
    6362
    64 ##    def test_slump_tsunami(self):
    65 ##
    66 ##        len = 4500.0
    67 ##        thk = 760.0
    68 ##        wid = 4500.0
    69 ##        dep = 1200.0
    70 ##        rad = 3330
    71 ##        dp = 0.23
    72 ##        th = 12
    73 ##        alpha = 0.0
    74 ##        x0 = 0
    75 ##        y0 = 0
    76 ##
    77 ##        slump = slump_tsunami(length=len, depth=dep, slope=th, thickness=thk,\
    78 ##                  radius=rad, dphi=dp, x0=x0, y0=y0, alpha=alpha)
    79 ##
    80 ##        assert allclose(slump.a3D, 9.82538623)
    81 ##        assert allclose(slump.wavelength, 3660.37606554)
    82 ##        assert allclose(slump.width, 4500.0)
    83 ##        assert allclose(slump.x0, 0.0)
    84 ##        assert allclose(slump.y0, 0.0)
    85 ##        assert allclose(slump.alpha, 0.0)
    86 ##        assert allclose(slump.kappa, 3.0)
    87 ##        assert allclose(slump.kappad, 1.0)
    88 ##
     63    def test_slump_tsunami(self):
     64
     65        length = 4500.0
     66        thk = 760.0
     67        wid = 4500.0
     68        dep = 1200.0
     69        rad = 3330
     70        dp = 0.23
     71        th = 12
     72        alpha = 0.0
     73        x0 = 0
     74        y0 = 0
     75        gamma = 1.85
     76
     77        slump = slump_tsunami(length, dep, th, wid, thk, rad, dp, x0, y0, alpha, gamma, scale=1.0)
     78
     79        assert allclose(slump.a3D, 9.82538623)
     80        assert allclose(slump.wavelength, 3660.37606554)
     81        assert allclose(slump.width, 4500.0)
     82        assert allclose(slump.x0, 0.0)
     83        assert allclose(slump.y0, 0.0)
     84        assert allclose(slump.alpha, 0.0)
     85        assert allclose(slump.kappa, 3.0)
     86        assert allclose(slump.kappad, 1.0)
     87
     88    def test_slide_tsunami_domain(self):
     89
     90        length = 600.0
     91        dep = 150.0
     92        th = 9.0
     93        thk = 15.0
     94        wid = 340.0
     95        kappa = 3.0
     96        kappad = 0.8
     97        x0 = 100000.
     98        y0 = x0
     99       
     100        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     101        polygon = [[0,0],[200000,0],[200000,200000],[0,200000]]
     102        create_mesh_from_regions(polygon,
     103                                 {'e0': [0], 'e1': [1], 'e2': [2], 'e3': [3]},
     104                                 maximum_triangle_area=5000000000,
     105                                 filename='test.msh',
     106                                 verbose = False)
     107
     108        from anuga.shallow_water import Domain
     109        domain = Domain('test.msh', use_cache = True, verbose = False)
     110
     111        slide = slide_tsunami(length, dep, th, x0, y0, \
     112                              wid, thk, kappa, kappad, \
     113                              domain=domain,verbose=False)
     114
     115        domain.set_quantity('stage', slide)
     116        stage = domain.get_quantity('stage')
     117        w = stage.get_values()
     118
     119##      check = [[-0.0 -0.0 -0.0],
     120##                 [-.189709745 -517.877716 -0.0],
     121##                 [-0.0 -0.0 -2.7695931e-08],
     122##                 [-0.0 -2.7695931e-08 -1.897097e-01]
     123##                 [-0.0 -517.877716 -0.0],
     124##                 [-0.0 -0.0 -0.0],
     125##                 [-0.0 -0.0 -0.0],
     126##                 [-0.0 -0.0 -0.0]]
     127
     128        assert allclose(min(min(w)), -517.877771593)
     129        assert allclose(max(max(w)), 0.0)
     130        assert allclose(slide.a3D, 518.38797486)
     131
    89132
    90133#-------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.