source: anuga_work/production/sydney_2006/smf_approx.py @ 4282

"""Script for running a tsunami inundation scenario for Sydney, NSW, Australia.

Source data such as elevation and boundary data is assumed to be available in
directories specified by project.py
The output sww file is stored in project.outputdir 

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data and a simulated submarine landslide.

Ole Nielsen and Duncan Gray, GA - 2005 and Adrian Hitchman and Jane Sexton, GA - 2006
"""

#-------------------------------------------------------------------------------# Import necessary modules
#-------------------------------------------------------------------------------

# Standard modules
import os
import time

from pylab import *
ion()
hold(True)

# Application specific imports
import project                 # Definition of file names and polygons
from anuga.shallow_water.smf import slide_tsunami, slump_tsunami  # Function for submarine mudslide 

def stats_slide(depth,length,width,thickness,slope,gamma):
    from math import sin, tan, radians, pi, sqrt, exp

    gravity = 9.8
    massco = 1.0
    dragco = 1.0
    psi = 0.0
    sint = sin(radians(slope))
    tant = tan(radians(slope))
    tanp = tan(radians(psi)) 
    
    a0 = gravity * sint * ((gamma-1)/(gamma+massco)) * (1-(tanp/tant))
    ut = sqrt((gravity*depth) * (length*sint/depth) \
                    * (pi*(gamma-1)/(2*dragco)) * (1-(tanp/tant)))
    s0 = ut**2 / a0
    t0 = ut / a0

    #calculate some parameters of the water displacement produced by the slide

    w = t0 * sqrt(gravity*depth)
    a2D = s0 * (0.0574 - (0.0431*sint)) \
             * thickness/length \
             * ((length*sint/depth)**1.25) \
             * (1 - exp(-2.2*(gamma-1)))
    a3D = a2D / (1 + (15.5*sqrt(depth/(length*sint))))
    
    #scaled_amp = a3D/a2D
    
    return a3D, w

def stats_slump(depth,length,width,thickness,slope,gamma):
    from math import sin, radians, sqrt

    radius = length**2 / (8.0 * thickness)    
    massco = 1.0
    dphi = 0.48
    gravity = 9.8
    sint = sin(radians(slope))

    s0 = radius * dphi / 2
    t0 = sqrt((radius*(gamma+massco)) / (gravity*(gamma-1)))
    a0 = s0 / t0**2
    um = s0 / t0

    #calculate some parameters of the water displacement produced by the slump

    w = t0 * sqrt(gravity*depth)
    a2D = s0 * (0.131/sint) \
             * thickness/length \
             * (length*sint/depth)**1.25 \
             * (length/radius)**0.63 * dphi**0.39 \
             * (1.47 - (0.35*(gamma-1))) * (gamma-1)
    a3D = a2D / (1. + (2.06*sqrt(depth/length)))

    #scaled_amp = a3D/a2D
    
    return a3D, w

which_one = 'slide'
#which_one = 'slump'
#par = 'depth'
#par = 'density'
par = 'slope'

depths = range(750,2500,50)
allgammas = arange(1.45,1.85,.05)
all_lengths=[7000,17000,1000]
all_slopes = arange(2,10,.2)
lengths = [16840,13500,7050]
widths = [8860,4350,3080]
thicknesses = [424,165,144]
slopes = [4,4,3]
gammas = [1.46,1.49,1.48]
smftype = ['Bulli','Shovel','Yacaaba']
smfdepths = [2087,968,1119]

if par == 'depth':
    
    for i in range(len(lengths)):
        amp = []
        wavelength = []
        for depth in depths:
            
            if which_one == 'slide':
                a, w = stats_slide(depth,lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
            else:
                a, w = stats_slump(depth,lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])

            amp.append(a)
            wavelength.append(w/1000.)

        figure(1)
        plot(depths,amp)
        xlabel('Depth (m)')
        ylabel('amplitude (m)')

        figure(2)
        plot(depths,wavelength)
        xlabel('Depth (m)')
        ylabel('wavelength (km)')

    figure(1)
    legend(smftype)
    savefig('depth_vs_amp')
    
    figure(2)
    legend(smftype)
    savefig('depth_vs_wavelength')

    adata = []
    wdata = []
    for i in range(len(smfdepths)):
        if which_one == 'slide':
            a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
        else:
            a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
        adata.append(a)
        wdata.append(w/1000.)

    figure(1)
    plot(smfdepths,adata,'^')
    savefig('depth_vs_amp_data%s'%which_one)
    figure(2)
    plot(smfdepths,wdata,'^')
    savefig('depth_vs_wavelength_data%s'%which_one)

if par == 'density':
    
    for i in range(len(lengths)):
        amp = []
        wavelength = []
        for gamma in allgammas:
            
            if which_one == 'slide':
                a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gamma)
            else:
                a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gamma)

            amp.append(a)
            wavelength.append(w/1000.)

        figure(1)
        plot(allgammas,amp)
        xlabel('Density (g/cm3)')
        ylabel('amplitude (m)')

        figure(2)
        plot(allgammas,wavelength)
        xlabel('Density (g/cm3)')
        ylabel('wavelength (km)')

    figure(1)
    legend(smftype)
    savefig('density_vs_amp')
    figure(2)
    legend(smftype)
    savefig('density_vs_wavelength')

    adata = []
    wdata = []
    for i in range(len(gammas)):
        if which_one == 'slide':
            a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
        else:
            a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
        adata.append(a)
        wdata.append(w/1000.)

    figure(1)
    plot(gammas,adata,'^')
    savefig('density_vs_amp_data%s'%which_one)
    figure(2)
    plot(gammas,wdata,'^')
    savefig('density_vs wavelength_data%s'%which_one)

if par == 'slope':
    
    for i in range(len(lengths)):
        amp = []
        wavelength = []
        for slope in all_slopes:
            
            if which_one == 'slide':
                a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slope,gammas[i])
            else:
                a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slope,gammas[i])

            amp.append(a)
            wavelength.append(w/1000.)

        figure(1)
        plot(all_slopes,amp)
        xlabel('Slope (deg)')
        ylabel('amplitude (m)')

        figure(2)
        plot(all_slopes,wavelength)
        xlabel('Slope (deg)')
        ylabel('wavelength (km)')

    figure(1)
    legend(smftype)
    savefig('slope_vs_amp')
    figure(2)
    legend(smftype)
    savefig('slope_vs_wavelength')

    adata = []
    wdata = []
    for i in range(len(lengths)):
        if which_one == 'slide':
            a, w = stats_slide(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
        else:
            a, w = stats_slump(smfdepths[i],lengths[i],widths[i],thicknesses[i],slopes[i],gammas[i])
        adata.append(a)
        wdata.append(w/1000.)

    figure(1)
    plot(slopes,adata,'^')
    savefig('slope_vs_amp_data%s'%which_one)
    figure(2)
    plot(slopes,wdata,'^')
    savefig('slope_vs_wavelength_data%s'%which_one)


close('all')

##import sys; sys.exit()
##tsunami_source = slide_tsunami(length=16840.0,
##                               width=8860.0,
##                               depth=2087.0,
##                               slope=4.0,
##                               thickness=424.0,
##                               gamma=1.46,
##                               #dphi=0.23,
##                               x0=project.slump_origin2[0], 
##                               y0=project.slump_origin2[1], 
##                               #alpha=0.0, 
##                               #domain=domain,
##                               verbose=True)
##        tsunami_source = slump_tsunami(length=16840.0,
##                                       width=8860.0,
##                                       #depth=2087.0,
##                                       slope=4.0,
##                                       thickness=424.0,
##                                       gamma=1.45,
##                                       #dphi=0.23,
##                                       #radius = 3330,
##                                       x0=project.slump_origin2[0], 
##                                       y0=project.slump_origin2[1], 
##                                       #alpha=0.0, 
##                                       #domain=domain,
##                                       verbose=True)