source: anuga_work/development/landslide_4/benchmark_problem4.py @ 5567

"""Script for running the 3rd Long Wave Workshop Benchmark Problem 4.

Bridgette Lewis and Jane Sexton, July 2008
"""

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

# Standard modules
import os
import time
from shutil import copy
from os.path import dirname, basename
from os import mkdir, access, F_OK, sep
import sys

# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary 
from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files

#----------------------------------------------------------------------------
# Create the triangular mesh based on the problem statement
#-------------------------------------------------------------------------------

from anuga.pmesh.mesh_interface import create_mesh_from_regions
remainder_res = 250000.
#local_res = 50000.
#coast_res = 500.
#interior_regions = [[project_slide.poly_syd1, local_res],
#                    [project_slide.poly_coast, coast_res]]
length = 10.
width = 10.
depth = 10.
bounding_polygon = [[0,width],[0,0],[length,0],[length,width]]
meshname = 'slide.msh'
from caching import cache
_ = cache(create_mesh_from_regions,
          bounding_polygon,
           {'boundary_tags': {'top': [0], 'left': [1], 
                              'bottom': [2], 'right': [3]},
           'maximum_triangle_area': remainder_res,
           'filename': meshname},
           #'interior_regions': interior_regions},
          verbose = True, evaluate=False)
print 'created mesh'

#-------------------------------------------------------------------------------                                 
# Setup computational domain
#-------------------------------------------------------------------------------                                 
domain = Domain(meshname, use_cache = True, verbose = True)

print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()
print domain.statistics()

domain.set_name('slide')
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
domain.set_minimum_storable_height(0.01)

#-------------------------------------------------------------------------------                                 
# Setup initial conditions
#-------------------------------------------------------------------------------
from math import tan
def topography(x,y):
    z = -depth + tan(0.5)*x
    return z
domain.set_quantity('stage', 0.0)
domain.set_quantity('friction', 0.0) 
domain.set_quantity('elevation', topography,
                    use_cache = True,
                    verbose = True,
                    alpha = 0.1
                    )

#-------------------------------------------------------------------------------
# Set up scenario (tsunami_source is a callable object used with set_quantity)
#-------------------------------------------------------------------------------
from smf import slide_tsunami

tsunami_source = slide_tsunami(length=2.,
                               width=1.,
                               depth=1.,
                               slope=45,
                               thickness=0.1, 
                               x0=length/2, 
                               y0=width/2, 
                               alpha=0.,
                                         #dx=0.001,
                               domain=domain)
                               #verbose=True)

#-------------------------------------------------------------------------------                                 
# Setup boundary conditions
#-------------------------------------------------------------------------------
print 'Available boundary tags', domain.get_boundary_tags()

Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([tide,0,0])

domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd} )


#-------------------------------------------------------------------------------                                 
# Evolve system through time
#-------------------------------------------------------------------------------
import time
from Numeric import allclose
from anuga.abstract_2d_finite_volumes.quantity import Quantity

t0 = time.time()

for t in domain.evolve(yieldstep = 1, finaltime = 10): 
    domain.write_time()
    stagestep = domain.get_quantity('stage') 

    if allclose(t, 1):
        slide = Quantity(domain)
        slide.set_values(tsunami_source)
        domain.set_quantity('stage', slide + stagestep)
    
print 'That took %.2f seconds' %(time.time()-t0)

print 'finished'