"""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.pmesh.mesh_interface import create_mesh_from_regions from anuga.abstract_2d_finite_volumes.util 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 # Application specific imports import aticket160_project as project_urs def run_model(**kwargs): tide = kwargs['tide'] alpha = kwargs['alpha'] friction = kwargs['friction'] time_thinning = kwargs['time_thinning'] scenario_name = kwargs['aa_scenario_name'] #------------------------------------------------------------------------------ # Copy scripts to time stamped output directory and capture screen # output to file #------------------------------------------------------------------------------ #copy script must be before screen_catcher print 'tide',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 + '.pts' kwargs['boundary_file']=project_urs.boundaries_dir_name + '.sww' # kwargs['Completed']='' print 'output_dir',kwargs['output_dir'] print 'myid',myid 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() # barrier() # filenames #boundaries_name = project_urs.boundaries_name meshes_dir_name = project_urs.meshes_dir_name+'.msh' #boundaries_dir_name = project_urs.boundaries_dir_name # tide = project_urs.tide # 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 #-------------------------------------------------------------------------- if myid == 0: print 'start create mesh from regions' create_mesh_from_regions(project_urs.poly_all, boundary_tags={'back': [4, 5], 'side': [3, 6], 'ocean': [0, 1, 2]}, maximum_triangle_area=project_urs.res_poly_all, interior_regions=project_urs.interior_regions, filename=meshes_dir_name, use_cache=False, verbose=True) barrier() #------------------------------------------------------------------------- # Setup computational domain #------------------------------------------------------------------------- print 'Setup computational domain' #from caching import cache #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True) #above don't work domain = Domain(meshes_dir_name, use_cache=False, verbose=True) # print 'domain id', id(domain) # print 'myhash', myhash(domain) print domain.statistics() # domain. print 'triangles',len(domain) kwargs['act_num_trigs']=len(domain) # print len(domain.mesh) #------------------------------------------------------------------------- # 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 = tide, geo_reference = domain.geo_reference) domain.set_quantity('stage', IC) domain.set_quantity('friction', friction) print 'Start Set quantity' domain.set_quantity('elevation', filename = kwargs['bathy_file'], use_cache = True, verbose = True, alpha = 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(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) domain.beta_h = 0 domain.tight_slope_limiters = 1 domain.H0=0.01 #------------------------------------------------------------------------- # 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=time_thinning, mean_stage=tide, # use_cache=True, verbose=True) kwargs['input_start_time']=domain.starttime print 'finished reading boundary file' Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([tide,0,0]) Bf = Dirichlet_boundary([5.,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() print "in time" for t in domain.evolve(yieldstep = 10, finaltime = 20): print "yeah", domain.write_time() if myid == 0: store_parameters(**kwargs) barrier() #------------------------------------------------------------- if __name__ == "__main__": run_model(file_name=project_urs.home+'detail.csv', aa_scenario_name=project_urs.scenario_name, ab_time=project_urs.time, res_factor= project_urs.res_factor, tide=project_urs.tide, user=project_urs.user, alpha = project_urs.alpha, friction=project_urs.friction, time_thinning = project_urs.time_thinning, dir_comment=project_urs.dir_comment)