"""Script for running a tsunami inundation scenario for Boca do Rio. Adopted from cairns files Joaquim Luis, 2007 """ #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ # Standard modules import os import time import sys # Related major packages from anuga.shallow_water import Domain from anuga.shallow_water import Transmissive_boundary from anuga.shallow_water import Reflective_boundary from anuga.shallow_water import Dirichlet_boundary from anuga.shallow_water import Time_boundary from anuga.shallow_water import File_boundary from anuga.shallow_water import Field_boundary from anuga.pmesh.mesh_interface import create_mesh_from_regions from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf from anuga.shallow_water.data_manager import dem2pts from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary from anuga.abstract_2d_finite_volumes.util import file_function #------------------------------------------------------------------------------ # Define scenario as either ... #------------------------------------------------------------------------------ scenario = 'br' basename = scenario + 'source' #------------------------------------------------------------------------------ # Preparation of topographic data # Convert ASC 2 DEM 2 PTS using source data and store result in source data #------------------------------------------------------------------------------ # Filenames dem_name = 'sub_region_bat' meshname = 'sub_region.msh' # Create DEM from asc data convert_dem_from_ascii2netcdf(dem_name, use_cache=False, verbose=True) # Create pts file for onshore DEM dem2pts(dem_name, use_cache=False, verbose=True) from anuga.utilities.polygon import read_polygon, plot_polygons, \ polygon_area, is_inside_polygon ############################### # Domain definitions ############################### # bounding polygon for study area bounding_polygon = read_polygon('out_rect.csv') print 'Area of bounding polygon in km?', polygon_area(bounding_polygon)/1000000.0 ############################### # Interior region definitions ############################### # interior polygons poly_shallow = read_polygon('shallow.csv') #------------------------------------------------------------------------------ # Create the triangular mesh based on overall clipping polygon with a tagged # boundary and interior regions defined in project.py along with # resolutions (maximal area of per triangle) for each polygon #------------------------------------------------------------------------------ remainder_res = 2500 shallow_res = 250 interior_regions = [[bounding_polygon, remainder_res], [poly_shallow, shallow_res]] create_mesh_from_regions(bounding_polygon, boundary_tags={'west': [0], 'south': [1], 'east': [2], 'north': [3]}, maximum_triangle_area=remainder_res, filename=meshname, interior_regions=interior_regions, use_cache=False, verbose=True) #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ domain = Domain(meshname, use_cache=False, verbose=True) print 'Number of triangles = ', len(domain) print 'The extent is ', domain.get_extent() print domain.statistics() domain.set_name(basename) domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) domain.set_minimum_storable_height(0.01) #------------------------------------------------------------------------------ # Setup initial conditions #------------------------------------------------------------------------------ tide = 0.0 domain.set_quantity('stage', tide) domain.set_quantity('friction', 0.025) domain.set_quantity('elevation', filename=dem_name + '.pts', use_cache=False, verbose=True, alpha=0.1) #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ # Create boundary function from timeseries provided in file #function1 = file_function('mareg.tms', domain, verbose=True) # Create and assign boundary objects #Bts1 = Transmissive_Momentum_Set_Stage_boundary(domain, function1) print 'Available boundary tags', domain.get_boundary_tags() #Bt = Transmissive_boundary(domain) # Continue all values on boundary Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([tide,0,0]) Bf = Field_boundary('br_swan.sww', domain, time_thinning=1, mean_stage=tide, use_cache=False, verbose=True) #Boundary condition for sww feed at the east boundary domain.set_boundary({'west': Br,'south':Bf,'east': Br,'north': Br}) #------------------------------------------------------------------------------ # Evolve system through time #------------------------------------------------------------------------------ import time t0 = time.time() from Numeric import allclose from anuga.abstract_2d_finite_volumes.quantity import Quantity # save every 1sec leading up to wave approaching land for t in domain.evolve(yieldstep = 1, finaltime = 960): domain.write_time() domain.write_boundary_statistics(tags = 'east') print 'That took %.2f seconds' %(time.time()-t0)