"""Run a tsunami inundation scenario for Busselton, 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_elevation.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_elevation.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 import os import os.path import time from time import localtime, strftime, gmtime # Related major packages from Scientific.IO.NetCDF import NetCDFFile #import numpy as num from anuga.interface import create_domain_from_regions from anuga.interface import Transmissive_stage_zero_momentum_boundary from anuga.interface import Dirichlet_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 anuga.utilities.system_tools import file_length from anuga.shallow_water.data_manager import start_screen_catcher from anuga.shallow_water.data_manager import copy_code_files from anuga.shallow_water.data_manager import urs2sts from anuga.utilities.polygon import read_polygon, Polygon_function from anuga.caching import cache import anuga.utilities.log as log # Application specific imports from setup_model import project import build_urs_boundary as bub #------------------------------------------------------------------------------- # Copy scripts to time stamped output directory and capture screen # output to file. Copy script must be before screen_catcher #------------------------------------------------------------------------------- copy_code_files(project.output_run, [__file__, os.path.join(os.path.dirname(project.__file__), project.__name__+'.py'), os.path.join(os.path.dirname(project.__file__), 'setup_model.py')], verbose=True ) #start_screen_catcher(project.output_run, 0, 1) #------------------------------------------------------------------------------- # 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 #------------------------------------------------------------------------------- log.critical('Create computational domain') # Create the STS file # FIXME (Ole): This is deadly dangerous if buildcode changes (as was the case 24th March 2009) # We need to use caching instead! log.critical( 'project.mux_data_folder=%s' % project.mux_data_folder) if not os.path.exists(project.event_sts + '.sts'): bub.build_urs_boundary(project.mux_input_filename, project.event_sts) # Read in boundary from ordered sts file event_sts = create_sts_boundary(project.event_sts) # Reading the landward defined points, this incorporates the original clipping # polygon minus the 100m contour landward_boundary = read_polygon(project.landward_boundary) # Combine sts polyline with landward points bounding_polygon_sts = event_sts + landward_boundary # Number of boundary segments num_ocean_segments = len(event_sts) - 1 # Number of landward_boundary points num_land_points = file_length(project.landward_boundary) # Boundary tags refer to project.landward_boundary # 4 points equals 5 segments start at N boundary_tags={'back': range(num_ocean_segments+1, num_ocean_segments+num_land_points), 'side': [num_ocean_segments, num_ocean_segments+num_land_points], 'ocean': range(num_ocean_segments)} # Build mesh and domain domain = create_domain_from_regions(bounding_polygon_sts, boundary_tags=boundary_tags, maximum_triangle_area=project.bounding_maxarea, interior_regions=project.interior_regions, mesh_filename=project.meshes, use_cache=True, verbose=True) log.critical(domain.statistics()) # FIXME(Ole): How can we make this more automatic? domain.geo_reference.zone = project.zone domain.set_name(project.scenario_name) domain.set_datadir(project.output_run) domain.set_minimum_storable_height(0.01) # Don't store depth less than 1cm #------------------------------------------------------------------------------- # Setup initial conditions #------------------------------------------------------------------------------- log.critical('Setup initial conditions') # Set the initial stage in the offcoast region only if project.land_initial_conditions: IC = Polygon_function(project.land_initial_conditions, default=project.tide, geo_reference=domain.geo_reference) else: IC = 0 domain.set_quantity('stage', IC, use_cache=True, verbose=True) domain.set_quantity('friction', project.friction) domain.set_quantity('elevation', filename=project.combined_elevation+'.pts', use_cache=True, verbose=True, alpha=project.alpha) if project.use_buildings: # Add buildings from file log.critical('Reading building polygons') building_polygons, building_heights = csv2building_polygons(project.building_polygon) #clipping_polygons=project.building_area_polygons) log.critical('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: log.critical(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 log.critical('Creating %d building polygons' % len(building_polygons)) buildings = cache(create_polygon_function, building_polygons, {'geo_reference': domain.geo_reference}, verbose=True) log.critical('Adding buildings') domain.add_quantity('elevation', buildings, use_cache=True, verbose=True) #------------------------------------------------------------------------------- # Setup boundary conditions #------------------------------------------------------------------------------- log.critical('Set boundary - available tags:' % domain.get_boundary_tags()) Br = Reflective_boundary(domain) Bs = Transmissive_stage_zero_momentum_boundary(domain) Bf = Field_boundary(project.event_sts+'.sts', domain, mean_stage=project.tide, time_thinning=1, default_boundary=Dirichlet_boundary([0, 0, 0]), boundary_polygon=bounding_polygon_sts, use_cache=True, verbose=True) domain.set_boundary({'back': Br, 'side': Bs, 'ocean': Bf}) #------------------------------------------------------------------------------- # Evolve system through time #------------------------------------------------------------------------------- t0 = time.time() # Skip over the first 6000 seconds for t in domain.evolve(yieldstep=2000, finaltime=6000): log.critical(domain.timestepping_statistics()) log.critical(domain.boundary_statistics(tags='ocean')) # Start detailed model for t in domain.evolve(yieldstep=project.yieldstep, finaltime=project.finaltime, skip_initial_step=True): log.critical(domain.timestepping_statistics()) log.critical(domain.boundary_statistics(tags='ocean')) log.critical('Simulation took %.2f seconds' %(time.time()-t0))