"""Script for running a tsunami inundation scenario for carnarvon, WA, Australia. The scenario is defined by a triangular mesh created from project.polygon, the elevation data is compiled into a pts file through build_carnarvon.py and a simulated tsunami is generated through an sts file from build_boundary.py. Input: sts file (build_boundary.py for respective event) pts file (build_carnarvon.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 Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006 Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008 """ #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ # Standard modules from os import sep import os from os.path import dirname, basename from os import mkdir, access, F_OK from shutil import copy import time import sys # Related major packages from anuga.shallow_water import Domain from anuga.shallow_water import Dirichlet_boundary from anuga.shallow_water import File_boundary from anuga.shallow_water import Reflective_boundary from anuga.shallow_water import Field_boundary from Numeric import allclose from anuga.shallow_water.data_manager import export_grid, create_sts_boundary from anuga.pmesh.mesh_interface import create_mesh_from_regions from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters from anuga_parallel.parallel_abstraction import get_processor_name from anuga.caching import myhash from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage from anuga.fit_interpolate.benchmark_least_squares import mem_usage from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter from polygon import Polygon_function # Application specific imports import project # Definition of file names and polygons numprocs = 1 myid = 0 def run_model(**kwargs): #------------------------------------------------------------------------------ # Copy scripts to time stamped output directory and capture screen # output to file #------------------------------------------------------------------------------ print "Processor Name:",get_processor_name() #copy script must be before screen_catcher print 'output_dir',kwargs['output_dir'] copy_code_files(kwargs['output_dir'],__file__, dirname(project.__file__)+sep+ project.__name__+'.py' ) store_parameters(**kwargs) start_screen_catcher(kwargs['output_dir'], myid, numprocs) print "Processor Name:",get_processor_name() #----------------------------------------------------------------------- # Domain definitions #----------------------------------------------------------------------- # Read in boundary from ordered sts file urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir_event,project.scenario_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 # boundary tags refer to project.landward 4 points equals 5 segments start at N boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)} #-------------------------------------------------------------------------- # 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 #-------------------------------------------------------------------------- # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it # causes problems with the ability to cache set quantity which takes alot of times print 'start create mesh from regions' create_mesh_from_regions(bounding_polygon, boundary_tags=boundary_tags, maximum_triangle_area=project.res_poly_all, interior_regions=project.interior_regions, filename=project.meshes_dir_name, use_cache=True, verbose=True) #------------------------------------------------------------------------- # Setup computational domain #------------------------------------------------------------------------- print 'Setup computational domain' domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True) print 'memory usage before del domain',mem_usage() print domain.statistics() print 'triangles',len(domain) kwargs['act_num_trigs']=len(domain) #------------------------------------------------------------------------- # Setup initial conditions #------------------------------------------------------------------------- print 'Setup initial conditions' # sets the initial stage in the offcoast region only IC = Polygon_function( [(project.poly_mainland, 0),(project.poly_island1, 0) ,(project.poly_island2, 0)], default = kwargs['tide'], geo_reference = domain.geo_reference) domain.set_quantity('stage', IC) #domain.set_quantity('stage',kwargs['tide'] ) domain.set_quantity('friction', kwargs['friction']) print 'Start Set quantity',kwargs['elevation_file'] domain.set_quantity('elevation', filename = kwargs['elevation_file'], use_cache = False, verbose = True, alpha = kwargs['alpha']) print 'Finished Set quantity' ## #------------------------------------------------------ ## # Distribute domain to implement parallelism !!! ## #------------------------------------------------------ ## ## if numprocs > 1: ## domain=distribute(domain) #------------------------------------------------------ # Set domain parameters #------------------------------------------------------ print 'domain id', id(domain) domain.set_name(kwargs['scenario_name']) domain.set_datadir(kwargs['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 domain.set_store_vertices_uniquely(False) domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) domain.tight_slope_limiters = 1 print 'domain id', id(domain) #------------------------------------------------------------------------- # Setup boundary conditions #------------------------------------------------------------------------- print 'Available boundary tags', domain.get_boundary_tags() print 'domain id', id(domain) boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name Br = Reflective_boundary(domain) Bs = Dirichlet_boundary([kwargs['tide'],0,0]) print 'Available boundary tags', domain.get_boundary_tags() Bf = Field_boundary(boundary_urs_out+'.sts', # Change from file_boundary domain, mean_stage=project.tide, time_thinning=1, default_boundary=Dirichlet_boundary(0,0,0]), use_cache=True, verbose=True, boundary_polygon=bounding_polygon) domain.set_boundary({'back': Br, 'side': Bs, 'ocean': Bf}) kwargs['input_start_time']=domain.starttime print'finish set boundary' #---------------------------------------------------------------------------- # Evolve system through time #-------------------------------------------------------------------- t0 = time.time() for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime'] ,skip_initial_step = False): domain.write_time() domain.write_boundary_statistics(tags = 'ocean') # these outputs should be checked with the resultant inundation map x, y = domain.get_maximum_inundation_location() q = domain.get_maximum_inundation_elevation() print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q) print 'Simulation took %.2f seconds' %(time.time()-t0) #kwargs 'completed' must be added to write the final parameters to file kwargs['completed']=str(time.time()-t0) store_parameters(**kwargs) print 'memory usage before del domain1',mem_usage() #------------------------------------------------------------- if __name__ == "__main__": kwargs={} kwargs['file_name']=project.dir_comment kwargs['finaltime']=project.finaltime kwargs['output_dir']=project.output_run_time_dir kwargs['elevation_file']=project.combined_dir_name+'.pts' kwargs['scenario_name']=project.scenario_name kwargs['tide']=project.tide kwargs['alpha'] = project.alpha kwargs['friction']=project.friction #kwargs['num_cpu']=numprocs #kwargs['host']=project.host #kwargs['starttime']=project.starttime #kwargs['yieldstep']=project.yieldstep #kwargs['user']=project.user #kwargs['time_thinning'] = project.time_thinning run_model(**kwargs)