"""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'