""" Script for running a breaking wave simulation of Jon Hinwoods wave tank. Note: this is based on the frinction_ua_flume_2006 structure. Duncan Gray, GA - 2007 """ #---------------------------------------------------------------------------- # Import necessary modules #---------------------------------------------------------------------------- # Standard modules import time from time import localtime, strftime import sys from shutil import copy from os import path, sep from os.path import dirname, join #, basename from Numeric import zeros, size, Float # Related major packages from anuga.shallow_water import Domain, Reflective_boundary, \ Dirichlet_boundary, Time_boundary, \ File_boundary, \ Transmissive_Momentum_Set_Stage_boundary from anuga.fit_interpolate.interpolate import interpolate_sww2csv from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \ file_function from anuga.shallow_water.data_manager import copy_code_files from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ import File_boundary_time # Scenario specific imports import project # Definition of file names and polygons import create_mesh from prepare_time_boundary import prepare_time_boundary from interp import interp class Elevation_function: def __init__(self, slope): self.xslope_position = [slope['xleft'][0],slope['xtoe'][0], slope['xbeach'][0],slope['xright'][0]] self.yslope_height = [slope['xleft'][1],slope['xtoe'][1], slope['xbeach'][1],slope['xright'][1]] def __call__(self, x,y): z = interp(self.yslope_height, self.xslope_position, x) return z def main(boundary_file, metadata_dic, boundary_path=None, friction=0.01, outputdir_name=None, run_type=0): basename = 'zz_' + metadata_dic['scenario_id'] if run_type == 1: outputdir_name += '_test' yieldstep = 1.0 finaltime = 15. maximum_triangle_area=0.1 elif run_type == 2: outputdir_name += '_test_long_time' yieldstep = 0.5 finaltime = None maximum_triangle_area=0.01 elif run_type == 3: outputdir_name += '_test_good_time_mesh' yieldstep = 0.1 finaltime = None maximum_triangle_area=0.001 elif run_type == 4: outputdir_name += '_good_tri_area_0.01_A' # this is not a test # Output will go to a file # The sww file will be interpolated yieldstep = 0.01 finaltime = None maximum_triangle_area=0.01 elif run_type == 5: outputdir_name += '_good_tri_area_0.001_A' # this is not a test # Output will go to a file # The sww file will be interpolated yieldstep = 0.01 finaltime = None maximum_triangle_area=0.001 elif run_type == 6: outputdir_name += '_good_tri_area_0.0001_A' # this is not a test # Output will go to a file # The sww file will be interpolated yieldstep = 0.01 finaltime = None maximum_triangle_area=0.0001 metadata_dic = set_z_origin_to_water_depth(metadata_dic) pro_instance = project.Project(['data','flumes','Hinwood_2008'], outputdir_name=outputdir_name) print "The output dir is", pro_instance.outputdir copy_code_files(pro_instance.outputdir,__file__, dirname(project.__file__) \ + sep + project.__name__+'.py') copy (pro_instance.codedir + 'run_dam.py', pro_instance.outputdir + 'run_dam.py') copy (pro_instance.codedir + 'create_mesh.py', pro_instance.outputdir + 'create_mesh.py') boundary_final_time = prepare_time_boundary(metadata_dic, pro_instance.raw_data_dir, pro_instance.boundarydir) if finaltime is None: finaltime = boundary_final_time # Boundary file manipulation if boundary_path is None: boundary_path = pro_instance.boundarydir boundary_file_path = join(boundary_path, boundary_file) # # Convert the boundary file, .csv to .tsm # try: # temp = open(boundary_file_path) # temp.close() # except IOError: # prepare_time_boundary(boundary_file_path) mesh_filename = pro_instance.meshdir + basename + '.msh' #-------------------------------------------------------------------------- # Copy scripts to output directory and capture screen # output to file #-------------------------------------------------------------------------- # creates copy of code in output dir if run_type >= 2: #start_screen_catcher(pro_instance.outputdir, rank, pypar.size()) start_screen_catcher(pro_instance.outputdir) print 'USER: ', pro_instance.user #------------------------------------------------------------------------- # Create the triangular mesh #------------------------------------------------------------------------- # this creates the mesh #gate_position = 12.0 create_mesh.generate(mesh_filename, metadata_dic, maximum_triangle_area=maximum_triangle_area) head,tail = path.split(mesh_filename) copy (mesh_filename, pro_instance.outputdir + tail ) #------------------------------------------------------------------------- # Setup computational domain #------------------------------------------------------------------------- domain = Domain(mesh_filename, use_cache = False, verbose = True) print 'Number of triangles = ', len(domain) print 'The extent is ', domain.get_extent() print domain.statistics() domain.set_name(basename) domain.set_datadir(pro_instance.outputdir) domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) domain.set_minimum_storable_height(0.001) #domain.set_store_vertices_uniquely(True) # for writting to sww #------------------------------------------------------------------------- # Setup initial conditions #------------------------------------------------------------------------- domain.set_quantity('stage', 0.) #the origin is the still water level domain.set_quantity('friction', friction) elevation_function = Elevation_function(metadata_dic) domain.set_quantity('elevation', elevation_function) print 'Available boundary tags', domain.get_boundary_tags() # Create boundary function from timeseries provided in file #function = file_function(project.boundary_file, domain, verbose=True) #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) try: function = file_function(boundary_file_path, domain, verbose=True) except IOError: msg = 'Run prepare_time_boundary.py. File "%s" could not be opened.'\ %(pro_instance.boundary_file) raise msg Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([0.3,0,0]) Bts = Time_boundary(domain, function) domain.set_boundary( {'wall': Br, 'wave': Bts} ) #domain.set_boundary( {'wall': Br, 'wave': Bd} ) #------------------------------------------------------------------------- # Evolve system through time #------------------------------------------------------------------------- t0 = time.time() # It seems that ANUGA can't handle a starttime that is >0. domain.starttime = 1.0 for t in domain.evolve(yieldstep, finaltime): domain.write_time() print 'That took %.2f seconds' %(time.time()-t0) print 'finished' flume_y_middle = 0.5 points = [] for gauge_x in metadata_dic['gauge_x']: points.append([gauge_x, flume_y_middle]) print "points",points #------------------------------------------------------------------------- # Calculate gauge info #------------------------------------------------------------------------- if run_type >= 1: id = metadata_dic['scenario_id'] interpolate_sww2csv(pro_instance.outputdir + basename +".sww", points, pro_instance.outputdir + "depth_" + id + ".csv", pro_instance.outputdir + "velocity_x_" + id + ".csv", pro_instance.outputdir + "velocity_y_" + id + ".csv") return pro_instance def set_z_origin_to_water_depth(seabed_coords): offset = seabed_coords['offshore_water_depth'] keys = ['xleft', 'xtoe', 'xbeach', 'xright'] for x in keys: seabed_coords[x][1] -= offset return seabed_coords #------------------------------------------------------------- if __name__ == "__main__": from scenarios import scenarios from slope import gauges_for_slope #from plot import plot #4 is 0.01 5 is 0.001 run_type = 5 #for run_data in [scenarios[5]]: for run_data in scenarios: main( run_data['scenario_id'] + '_boundary.tsm' , run_data, run_type = run_type, outputdir_name=run_data['scenario_id']) gauges_for_slope() #plot()