"""Script for running a tsunami inundation scenario for Patong Beach, Thailand. The scenario is defined by a triangular mesh created from project.polygon, the elevation data is compiled into a pts file through build_patong.py and a simulated tsunami for the 2004 event is generated through an sts file from build_boundary.py. Input: sts file (build_boundary.py) pts file (build_patong.py) information from project file Outputs: sww file stored in project.output_run_time_dir The export_results_all.py and get_timeseries.py is reliant on the outputs of this script """ #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ # Standard modules from os import sep import os from os.path import dirname import time # Related major packages from anuga.interface import create_domain_from_regions from anuga.interface import Dirichlet_boundary from anuga.interface import Transmissive_stage_zero_momentum_boundary from anuga.interface import Reflective_boundary from anuga.interface import Field_boundary from anuga.interface import create_sts_boundary from anuga.interface import csv2building_polygons from file_length import file_length from anuga.caching import cache from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters from anuga.utilities.polygon import read_polygon from polygon import Polygon_function # Application specific imports import project # Definition of file names and polygons numprocs = 1 myid = 0 #------------------------------------------------------------------------------ # Copy scripts to time stamped output directory and capture screen # output to file #------------------------------------------------------------------------------ #copy script must be before screen_catcher output_dir = project.output_run_time_dir print 'output_dir', output_dir copy_code_files(output_dir,__file__, dirname(project.__file__)+sep+ project.__name__+'.py' ) start_screen_catcher(output_dir, myid, numprocs) #----------------------------------------------------------------------- # Domain definitions #----------------------------------------------------------------------- # Read in boundary from ordered sts file urs_boundary_name = os.path.join(project.boundaries_dir, project.scenario_name) urs_bounding_polygon = create_sts_boundary(urs_boundary_name) # Reading the landward defined points, this incorporates the original clipping # polygon minus the 100m contour landward_bounding_polygon = read_polygon(project.landward_dir) # Combine sts polyline with landward points bounding_polygon = urs_bounding_polygon + landward_bounding_polygon # Counting segments N = len(urs_bounding_polygon)-1 # Number of landward_boundary points M = file_length(project.landward_dir) # Boundary tags refer to project.landward_boundary # 4 points equals 5 segments start at N boundary_tags={'back': range(N+1,N+M), 'side': [N,N+M], 'ocean': range(N)} #-------------------------------------------------------------------------- # Create the computational domain 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 #-------------------------------------------------------------------------- domain = create_domain_from_regions(bounding_polygon, boundary_tags=boundary_tags, maximum_triangle_area=project.res_poly_all, interior_regions=project.interior_regions, mesh_filename=project.meshes_dir_name, use_cache=True, verbose=True) print domain.statistics() #------------------------------------------------------------------------- # Setup initial conditions #------------------------------------------------------------------------- print 'Setup initial conditions' # sets the initial stage in the offcoast region only IC = Polygon_function([(project.poly_mainland, 0)], default=project.tide, geo_reference=domain.geo_reference) print 'stage' domain.set_quantity('stage', IC, use_cache=True, verbose=True) print 'friction' domain.set_quantity('friction', project.friction) print 'elevation' domain.set_quantity('elevation', filename=project.combined_dir_name+'.pts', alpha=project.alpha, use_cache=True, verbose=True) if project.use_buildings: # Add buildings from file print 'Reading building polygons' building_polygons, building_heights = csv2building_polygons(project.building_polygon_file) #clipping_polygons=project.building_area_polygons) print 'Creating %d building polygons' % len(building_polygons) def create_polygon_function(building_polygons, geo_reference=None): L = [] for i, key in enumerate(building_polygons): if i%100==0: print i poly = building_polygons[key] elev = building_heights[key] L.append((poly, elev)) buildings = Polygon_function(L, default=0.0, geo_reference=geo_reference) return buildings buildings = cache(create_polygon_function, building_polygons, {'geo_reference': domain.geo_reference}, verbose=True) print 'Adding buildings' domain.add_quantity('elevation', buildings, use_cache=True, verbose=True) #------------------------------------------------------ # Distribute domain to implement parallelism !!! #------------------------------------------------------ # if numprocs > 1: # domain=distribute(domain) #------------------------------------------------------ # Set domain parameters #------------------------------------------------------ domain.set_name(project.scenario_name) domain.set_datadir(output_dir) domain.set_default_order(2) # Apply second order scheme domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm #------------------------------------------------------------------------- # Setup boundary conditions #------------------------------------------------------------------------- print 'Set boundary conditions' Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([project.tide,0,0]) Bt = Transmissive_stage_zero_momentum_boundary(domain) Bf = Field_boundary(urs_boundary_name+'.sts', domain, mean_stage= project.tide, time_thinning=project.time_thinning, default_boundary=Bd, use_cache=True, verbose=True, boundary_polygon=bounding_polygon) domain.set_boundary({'back': Br, 'side': Bt, 'ocean': Bf}) #---------------------------------------------------------------------------- # Evolve system through time #-------------------------------------------------------------------- t0 = time.time() # Skip over the first 6000 seconds for t in domain.evolve(yieldstep=2000, finaltime=6000): print domain.timestepping_statistics() print domain.boundary_statistics(tags='ocean') # Start detailed model for t in domain.evolve(yieldstep=project.yieldstep, finaltime=project.finaltime, skip_initial_step=True): print domain.timestepping_statistics() print domain.boundary_statistics(tags='ocean') print 'Simulation took %.2f seconds' %(time.time()-t0)