"""Script for running a dam break simulation of UQ's dam break tank. Ole Nielsen and Duncan Gray, GA - 2006 """ #---------------------------------------------------------------------------- # Import necessary modules #---------------------------------------------------------------------------- # Standard modules import time import sys from shutil import copy from os import path # Related major packages from anuga.shallow_water import Domain, Reflective_boundary, \ Dirichlet_boundary, Time_boundary, File_boundary from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \ copy_code_files from anuga.abstract_2d_finite_volumes.region import Set_region from anuga.fit_interpolate.interpolate import interpolate_sww2csv import create_mesh #from anuga.pyvolution.data_manager import timefile2netcdf from math import cos from math import cosh from math import sin from math import sqrt from math import fabs # Application specific imports import project # Definition of file names and polygons def main(): #solitary wave solution parameters slopes = [0.1] #plane beach slope (0 if bathymetry is read from file) H_d_ratios = [0.00092] #ratio wave height to initial depth # pos = [0.002, 0.02, 0.08, 0.32] for slope in slopes: for H_d_ratio in H_d_ratios: # for po in pos: scenario(slope, H_d_ratio) def scenario(slope, H_d_ratio): longshore_length = 150 offshore_depth = 98 offshore_length = offshore_depth/slope bay_length = offshore_depth/4/slope period = 20 x1 = 4000 #------------------------------------------------------------------------- # Setup archiving of simulations #------------------------------------------------------------------------- id = 'Sine_wave_plane_beach'+'_s'+ str(slope)+ '_wr'+ str(H_d_ratio) copy (project.codedirname, project.outputtimedir + 'project.py') run_name = 'run_dam.py' run_name_out = 'run_dam'+id+'.py' copy (project.codedir + 'run_dam.py', project.outputtimedir + run_name_out) copy (project.codedir + 'create_mesh.py', project.outputtimedir + 'create_mesh.py') print'output dir', project.outputtimedir #FIXME this isn't working #normal screen output is stored in screen_output_name = project.outputtimedir + "screen_output.txt" screen_error_name = project.outputtimedir + "screen_error.txt" #------------------------------------------------------------------------- # Create the triangular mesh #------------------------------------------------------------------------- create_mesh.generate(project.mesh_filename, offshore_length, bay_length, longshore_length, is_course=True) # this creates the mesh #is_course=False) # this creates the mesh head,tail = path.split(project.mesh_filename) copy (project.mesh_filename, project.outputtimedir + tail ) #------------------------------------------------------------------------- # Setup computational domain #------------------------------------------------------------------------- domain = Domain(project.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(project.basename + id) domain.set_datadir(project.outputtimedir) domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) domain.set_minimum_storable_height(0.01) domain.set_store_vertices_uniquely(True) # for writing to sww #------------------------------------------------------------------------- # Setup initial conditions #------------------------------------------------------------------------- g = 9.81 def elevation_tilt(x, y): return x*slope domain.set_quantity('stage', filename = 'stage.xya', alpha = 0.001,verbose = True, use_cache = True) domain.set_quantity('friction', 0) domain.set_quantity('elevation', elevation_tilt) #bathymetry read from file #domain.set_quantity('elevation', filename = 'elevation.xya', # alpha = 10.0,verbose = True, use_cache = True) print 'Available boundary tags', domain.get_boundary_tags() Br = Reflective_boundary(domain) #bore described by function (solitary wave solution) W = Time_boundary(domain = domain, f=lambda t: [-offshore_depth*H_d_ratio*cos(2*3.14*t/period), 0, 0]) #sine wave # f=lambda t: [offshore_depth*H_d_ratio/cosh(sqrt(3*H_d_ratio*po/4)* # (sqrt(g/offshore_depth)*t-x1/offshore_depth))/cosh(sqrt(3*H_d_ratio*po/4)* # (sqrt(g/offshore_depth)*t-x1/offshore_depth)), 0, 0]) domain.set_boundary( {'wall': Br, 'edge': Br, 'back':W} ) #------------------------------------------------------------------------- # Evolve system through time #------------------------------------------------------------------------- t0 = time.time() for t in domain.evolve(yieldstep = 1, finaltime =200): #enter timestep and final time domain.write_time() print 'That took %.2f seconds' %(time.time()-t0) print 'finished' points = [[-offshore_length, longshore_length / 2], [(-offshore_length)*0.95, longshore_length / 2], [(-offshore_length)*0.9, longshore_length / 2], [(-offshore_length)*0.85, longshore_length / 2], [(-offshore_length)*0.8, longshore_length / 2], [(-offshore_length)*0.75, longshore_length / 2], [(-offshore_length)*0.7, longshore_length / 2], [(-offshore_length)*0.65, longshore_length / 2], [(-offshore_length)*0.6, longshore_length / 2], [(-offshore_length)*0.55, longshore_length / 2], [(-offshore_length)*0.5, longshore_length / 2], [(-offshore_length)*0.45, longshore_length / 2], [(-offshore_length)*0.4, longshore_length / 2], [(-offshore_length)*0.35, longshore_length / 2], [(-offshore_length)*0.3, longshore_length / 2], [(-offshore_length)*0.25, longshore_length / 2], [(-offshore_length)*0.2, longshore_length / 2], [(-offshore_length)*0.15, longshore_length / 2], [(-offshore_length)*0.1, longshore_length / 2], [(-offshore_length)*0.09, longshore_length / 2], [(-offshore_length)*0.08, longshore_length / 2], [(-offshore_length)*0.07, longshore_length / 2], [(-offshore_length)*0.06, longshore_length / 2], [(-offshore_length)*0.05, longshore_length / 2], [(-offshore_length)*0.04, longshore_length / 2], [(-offshore_length)*0.03, longshore_length / 2], [(-offshore_length)*0.02, longshore_length / 2], [(-offshore_length)*0.01, longshore_length / 2], [(-offshore_length)*0, longshore_length / 2], [(offshore_length)*0.01, longshore_length / 2], [(offshore_length)*0.02, longshore_length / 2], [(offshore_length)*0.02, longshore_length / 2], [(offshore_length)*0.03, longshore_length / 2], [(offshore_length)*0.04, longshore_length / 2]] #------------------------------------------------------------------------- # Calculate gauge info #------------------------------------------------------------------------- interpolate_sww2csv(project.outputtimedir + project.basename \ + id+".sww", points, project.depth_filename + id + '.csv', project.velocity_x_filename + id + '.csv', project.velocity_y_filename + id + '.csv') #------------------------------------------------------------- if __name__ == "__main__": main()