""" 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.012, # planed wood. http://www.lmnoeng.com/manningn.htm outputdir_name=None, run_type=0, width=1.0, use_limits=True, end_tag = '_limiterD'): basename = 'zz_' + metadata_dic['scenario_id'] if run_type == 1: yieldstep = 1.0 finaltime = 15. maximum_triangle_area=0.1 outputdir_name += '_test' elif run_type == 2: yieldstep = 0.5 finaltime = None maximum_triangle_area=0.01 outputdir_name += '_test_long_time' elif run_type == 3: yieldstep = 0.1 finaltime = None maximum_triangle_area=0.01 #outputdir_name += '_yieldstep_0.1' elif run_type == 4: # 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 #outputdir_name += '_good' elif run_type == 5: # 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 #outputdir_name += '_good' elif run_type == 6: # 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 #outputdir_name += '_good' if use_limits is True: outputdir_name += '_lmts' else: outputdir_name += '_nolmts' outputdir_name += '_wdth_' + str(width) outputdir_name += '_z_' + str(friction) outputdir_name += '_ys_' + str(yieldstep) outputdir_name += '_mta_' + str(maximum_triangle_area) outputdir_name += end_tag 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) #return pro_instance if finaltime is None: finaltime = boundary_final_time - 0.1 # Edge boundary problems # 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, width=width, 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.0001) if use_limits is True: domain.set_default_order(2) # Use second order spatial scheme domain.set_timestepping_method('rk2') domain.use_edge_limiter = True domain.tight_slope_limiters = True domain.beta_w = 0.6 domain.beta_uh = 0.6 domain.beta_vh = 0.6 #------------------------------------------------------------------------- # 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 #!!! what was this doing? for t in domain.evolve(yieldstep, finaltime): domain.write_time() print 'That took %.2f seconds' %(time.time()-t0) print 'finished' flume_y_middle = 0.0 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'] + ".csv" interpolate_sww2csv(pro_instance.outputdir + basename +".sww", points, pro_instance.outputdir + "depth_" + id, pro_instance.outputdir + "velocity_x_" + id, pro_instance.outputdir + "velocity_y_" + id, pro_instance.outputdir + "stage_" + id, pro_instance.outputdir + "froude_" + id) 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 # 1 is fast and dirty # 4 is 0.01 # 5 is 0.001 # 6 is 0.0001 #run_type = 1 run_type = 5 #for run_data in [scenarios[5]]: #scenarios = scenarios[2:] #scenarios = [scenarios[0]] width = 1.0 width = 0.1 for run_data in scenarios: pro_instance = main( run_data['scenario_id'] + '_boundary.tsm' , run_data, width=width, run_type=run_type, outputdir_name=run_data['scenario_id'], use_limits=False, friction=0.012, end_tag='_G') #gauges_for_slope(pro_instance.outputdir,[run_data])