""" 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 = 0.1 finaltime = 15. maximum_triangle_area=0.01 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 >= 2: 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__": #slopes = [[-4.5,0.0],[0.0,0.0],[1.285,0.090],[16.1,.960]] # Note, gauge A has been removed, since it is used as the # boundary. run_type = 6 # T1R5 run_data = {'xleft':[-3.106,0.0], # Av' of ADV and Gauge A 'xtoe':[0.0,0.0], 'xbeach':[1.285,0.090], 'xright':[16.1,.960], 'offshore_water_depth':.4, 'scenario_id':'T1R5', 'gauge_names':['B','1','2','3','4','5','6','7','8', '9','10','11','12','13','14'], 'gauge_x':[-0.68, 1.572, 2.572, 3.572, 4.572, 5.572, 6.572, 7.572, 8.572, 9.572, 10.572, 11.572, 12.572, 13.572, 14.572], 'gauge_bed_elevation':[-0.400000, -0.293158, -0.234473, -0.175788, -0.117104, -0.058419, 0.000266, 0.058950, 0.117635, 0.176320, 0.235004, 0.293689, 0.352374, 0.411058, 0.469743] } main( run_data['scenario_id'] + '_boundary.tsm' , run_data, run_type = run_type, outputdir_name=run_data['scenario_id']) # T1R3 run_data = {'xleft':[-3.106,0.0], # Av' of ADV and Gauge A 'xtoe':[0.0,0.0], 'xbeach':[1.285,0.090], 'xright':[16.1,.960], 'offshore_water_depth':.4, 'scenario_id':'T1R3', 'gauge_names':['B','1','2','3','4','5','6','7','8', '9','10','11','12','13','14'], 'gauge_x':[-0.68, 1.572, 2.572, 3.572, 4.572, 5.572, 6.572, 7.572, 8.572, 9.572, 10.572, 11.572, 12.572, 13.572, 14.572], 'gauge_bed_elevation':[-0.400000, -0.293158, -0.234473, -0.175788, -0.117104, -0.058419, 0.000266, 0.058950, 0.117635, 0.176320, 0.235004, 0.293689, 0.352374, 0.411058, 0.469743] } main( run_data['scenario_id'] + '_boundary.tsm' , run_data, run_type = run_type, outputdir_name=run_data['scenario_id']) # #T2R7 # xleft is different run_data = {'xleft':[-4.586,0.0], # Av' of ADV and Gauge A 'xtoe':[0.0,0.0], 'xbeach':[1.285,0.090], 'xright':[16.1,.960], 'offshore_water_depth':.4, 'scenario_id':'T2R7', 'gauge_names':['B','1','2','3','4','5','6','7','8', '9','10','11','12','13','14'], 'gauge_x':[-0.68, 1.572, 2.572, 3.572, 4.572, 5.572, 6.572, 7.572, 8.572, 9.572, 10.572, 11.572, 12.572, 13.572, 14.572], 'gauge_bed_elevation':[-0.400000, -0.293158, -0.234473, -0.175788, -0.117104, -0.058419, 0.000266, 0.058950, 0.117635, 0.176320, 0.235004, 0.293689, 0.352374, 0.411058, 0.469743] } main( run_data['scenario_id'] + '_boundary.tsm' , run_data, run_type = run_type, outputdir_name=run_data['scenario_id']) # #T2R8 # xleft is different run_data = {'xleft':[-4.586,0.0], # Av' of ADV and Gauge A 'xtoe':[0.0,0.0], 'xbeach':[1.285,0.090], 'xright':[16.1,.960], 'offshore_water_depth':.4, 'scenario_id':'T2R8', 'gauge_names':['B','1','2','3','4','5','6','7','8', '9','10','11','12','13','14'], 'gauge_x':[-0.68, 1.572, 2.572, 3.572, 4.572, 5.572, 6.572, 7.572, 8.572, 9.572, 10.572, 11.572, 12.572, 13.572, 14.572], 'gauge_bed_elevation':[-0.400000, -0.293158, -0.234473, -0.175788, -0.117104, -0.058419, 0.000266, 0.058950, 0.117635, 0.176320, 0.235004, 0.293689, 0.352374, 0.411058, 0.469743] } main( run_data['scenario_id'] + '_boundary.tsm' , run_data, run_type = run_type, outputdir_name=run_data['scenario_id']) # #T3R29 # xleft is different run_data = {'xleft':[-3.875,0.0], # Av' of ADV and Gauge A 'xtoe':[0.0,0.0], 'xbeach':[1.285,0.090], 'xright':[16.1,.440], 'offshore_water_depth':.336, 'scenario_id':'T3R29', 'gauge_names':['1','2','3','4','5','6','7','8', '9','10','11','12','13','14','B'], 'gauge_x':[1.572, 2.572, 3.572, 4.572, 5.572, 6.572, 7.572, 8.572, 9.572, 10.572, 11.572, 12.572, 13.572, 14.572, -0.325], 'gauge_bed_elevation':[-0.237263, -0.213789, -0.190315, -0.166841, -0.143368, -0.119894, -0.096420, -0.072946, -0.049472, -0.025998, -0.002524, 0.020949, 0.044423, 0.067897, -0.336000] } main( run_data['scenario_id'] + '_boundary.tsm' , run_data, run_type = run_type, outputdir_name=run_data['scenario_id']) # #T3R28 run_data = {'xleft':[-3.875,0.0], # Av' of ADV and Gauge A 'xtoe':[0.0,0.0], 'xbeach':[1.285,0.090], 'xright':[16.1,.440], 'offshore_water_depth':.336, 'scenario_id':'T3R28', 'gauge_names':['1','2','3','4','5','6','7','8', '9','10','11','12','13','14','B'], 'gauge_x':[1.572, 2.572, 3.572, 4.572, 5.572, 6.572, 7.572, 8.572, 9.572, 10.572, 11.572, 12.572, 13.572, 14.572, -0.325], 'gauge_bed_elevation':[-0.237263, -0.213789, -0.190315, -0.166841, -0.143368, -0.119894, -0.096420, -0.072946, -0.049472, -0.025998, -0.002524, 0.020949, 0.044423, 0.067897, -0.336000] } main( run_data['scenario_id'] + '_boundary.tsm' , run_data, run_type = run_type, outputdir_name=run_data['scenario_id']) # #T4R31 # xleft is different run_data = {'xleft':[-2.43,0.0], # Av' of ADV and Gauge A 'xtoe':[0.0,0.0], 'xbeach':[1.285,0.090], 'xright':[16.1,.440], 'offshore_water_depth':.336, 'scenario_id':'T4R31', 'gauge_names':['1','2','3','4','5','6','7','8', '9','10','11','12','13','14','B'], 'gauge_x':[1.572, 2.572, 3.572, 4.572, 5.572, 6.572, 7.572, 8.572, 9.572, 10.572, 11.572, 12.572, 13.572, 14.572, -0.325], 'gauge_bed_elevation':[-0.237263, -0.213789, -0.190315, -0.166841, -0.143368, -0.119894, -0.096420, -0.072946, -0.049472, -0.025998, -0.002524, 0.020949, 0.044423, 0.067897, -0.336000] } main( run_data['scenario_id'] + '_boundary.tsm' , run_data, run_type = run_type, outputdir_name=run_data['scenario_id']) # #T4R32 run_data = {'xleft':[-2.43,0.0], # Av' of ADV and Gauge A 'xtoe':[0.0,0.0], 'xbeach':[1.285,0.090], 'xright':[16.1,.440], 'offshore_water_depth':.336, 'scenario_id':'T4R32', 'gauge_names':['1','2','3','4','5','6','7','8', '9','10','11','12','13','14','B'], 'gauge_x':[1.572, 2.572, 3.572, 4.572, 5.572, 6.572, 7.572, 8.572, 9.572, 10.572, 11.572, 12.572, 13.572, 14.572, -0.325], 'gauge_bed_elevation':[-0.237263, -0.213789, -0.190315, -0.166841, -0.143368, -0.119894, -0.096420, -0.072946, -0.049472, -0.025998, -0.002524, 0.020949, 0.044423, 0.067897, -0.336000] } main( run_data['scenario_id'] + '_boundary.tsm' , run_data, run_type = run_type, outputdir_name=run_data['scenario_id'])