"""Script for running tsunami inundation scenario for Dampier, WA, Australia. Source data such as elevation and boundary data is assumed to be available in directories specified by project_urs.py The output sww file is stored in project_urs.output_run_time_dir The scenario is defined by a triangular mesh created from project_urs.polygon, the elevation data and a simulated tsunami generated with URS code. Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006 """ #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ # Standard modules from os import sep 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 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_api import distribute, numprocs, myid, barrier 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 # Application specific imports import project_urs # Definition of file names and polygons 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 kwargs ''' print 'tide',kwargs['tide'] kwargs['est_num_trigs']=project_urs.trigs_min kwargs['num_cpu']=numprocs kwargs['host']=project_urs.host kwargs['res_factor']=project_urs.res_factor kwargs['starttime']=project_urs.starttime kwargs['yieldstep']=project_urs.yieldstep kwargs['finaltime']=project_urs.finaltime kwargs['output_dir']=project_urs.output_run_time_dir kwargs['bathy_file']=project_urs.combined_dir_name+'.txt' # kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts' kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww' ''' print 'output_dir',kwargs['output_dir'] if myid == 0: copy_code_files(kwargs['output_dir'],__file__, dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' ) store_parameters(**kwargs) barrier() start_screen_catcher(kwargs['output_dir'], myid, numprocs) print "Processor Name:",get_processor_name() # filenames # meshes_dir_name = project_urs.meshes_dir_name+'.msh' # creates copy of code in output dir print 'min triangles', project_urs.trigs_min, print 'Note: This is generally about 20% less than the final amount' #-------------------------------------------------------------------------- # Create the triangular mesh based on overall clipping polygon with a # tagged # boundary and interior regions defined in project_urs.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 if myid == 0: print 'start create mesh from regions' create_mesh_from_regions(project_urs.poly_all, boundary_tags=project_urs.boundary_tags, maximum_triangle_area=project_urs.res_poly_all, interior_regions=project_urs.interior_regions, filename=project_urs.meshes_dir_name+'.msh', use_cache=False, verbose=True) barrier() #------------------------------------------------------------------------- # Setup computational domain #------------------------------------------------------------------------- print 'Setup computational domain' #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True) #above don't work domain = Domain(project_urs.meshes_dir_name+'.msh', 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 #------------------------------------------------------------------------- if myid == 0: print 'Setup initial conditions' from polygon import Polygon_function #following sets the stage/water to be offcoast only # IC = Polygon_function( [(project_urs.poly_mainland, -1.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['bathy_file'] domain.set_quantity('elevation', filename = kwargs['bathy_file'], use_cache = True, verbose = True, alpha = kwargs['alpha']) print 'Finished Set quantity' barrier() #------------------------------------------------------ # Distribute domain to implement parallelism !!! #------------------------------------------------------ if numprocs > 1: domain=distribute(domain) #------------------------------------------------------ # Set domain parameters #------------------------------------------------------ print 'domain id', id(domain) domain.set_name(kwargs['aa_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.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK) print 'domain id', id(domain) domain.beta_h = 0 #domain.tight_slope_limiters = 1 #------------------------------------------------------------------------- # Setup boundary conditions #------------------------------------------------------------------------- print 'Available boundary tags', domain.get_boundary_tags() print 'domain id', id(domain) #print 'Reading Boundary file',project_urs.boundaries_dir_namea + '.sww' Bf = Field_boundary(kwargs['boundary_file'], domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], use_cache=True, verbose=True) kwargs['input_start_time']=domain.starttime print 'finished reading boundary file' Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([kwargs['tide'],0,0]) print'set_boundary' domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bf}) print'finish set boundary' #---------------------------------------------------------------------------- # Evolve system through time #-------------------------------------------------------------------- t0 = time.time() for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): domain.write_time() domain.write_boundary_statistics(tags = 'ocean') for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 21600 ,skip_initial_step = True): domain.write_time() domain.write_boundary_statistics(tags = 'ocean') for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime'] ,skip_initial_step = True): domain.write_time() domain.write_boundary_statistics(tags = 'ocean') 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 'That took %.2f seconds' %(time.time()-t0) #kwargs 'completed' must be added to write the final parameters to file kwargs['completed']=str(time.time()-t0) if myid==0: store_parameters(**kwargs) barrier print 'memory usage before del domain1',mem_usage() def export_model(**kwargs): #store_parameters(**kwargs) print 'memory usage before del domain',mem_usage() #del domain print 'memory usage after del domain',mem_usage() swwfile = kwargs['output_dir']+kwargs['aa_scenario_name'] print'swwfile',swwfile export_grid(swwfile, extra_name_out = 'town', quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation timestep = None, reduction = max, cellsize = 25, NODATA_value = -9999, easting_min = project_urs.eastingmin, easting_max = project_urs.eastingmax, northing_min = project_urs.northingmin, northing_max = project_urs.northingmax, verbose = False, origin = None, datum = 'GDA94', format = 'asc') inundation_damage(swwfile+'.sww', project_urs.buildings_filename, project_urs.buildings_filename_out) #------------------------------------------------------------- if __name__ == "__main__": kwargs={} kwargs['est_num_trigs']=project_urs.trigs_min kwargs['num_cpu']=numprocs kwargs['host']=project_urs.host kwargs['res_factor']=project_urs.res_factor kwargs['starttime']=project_urs.starttime kwargs['yieldstep']=project_urs.yieldstep kwargs['finaltime']=project_urs.finaltime kwargs['output_dir']=project_urs.output_run_time_dir kwargs['bathy_file']=project_urs.combined_dir_name+'.txt' # kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts' kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww' kwargs['file_name']=project_urs.home+'detail.csv' kwargs['aa_scenario_name']=project_urs.scenario_name kwargs['ab_time']=project_urs.time kwargs['res_factor']= project_urs.res_factor kwargs['tide']=project_urs.tide kwargs['user']=project_urs.user kwargs['alpha'] = project_urs.alpha kwargs['friction']=project_urs.friction kwargs['time_thinning'] = project_urs.time_thinning kwargs['dir_comment']=project_urs.dir_comment run_model(**kwargs) if myid==0: export_model(**kwargs) barrier