"""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, sep from os.path import dirname #, basename # Related major packages from anuga.shallow_water import Domain, Reflective_boundary, \ Dirichlet_boundary, Time_boundary, File_boundary from anuga.abstract_2d_finite_volumes.region import Set_region from anuga.fit_interpolate.interpolate import interpolate_sww2csv from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \ copy_code_files # Scenario specific imports import project # Definition of file names and polygons import create_mesh def main(friction): #-------------------------------------------------------------------------- # Copy scripts to time stamped output directory and capture screen # output to file #-------------------------------------------------------------------------- # creates copy of code in output dir print "The output dir is", project.outputtimedir copy_code_files(project.outputtimedir,__file__, dirname(project.__file__)+sep+ project.__name__+'.py' ) copy (project.codedir + 'create_mesh.py', project.outputtimedir + 'create_mesh.py') myid = 0 numprocs = 1 #start_screen_catcher(project.outputtimedir, myid, numprocs) print 'USER: ', project.user #------------------------------------------------------------------------- # Create the triangular mesh #------------------------------------------------------------------------- gate_position = 0.85 create_mesh.generate(project.mesh_filename, gate_position, #is_coarse=True) # this creates the mesh is_coarse=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) 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 writting to sww #------------------------------------------------------------------------- # Setup initial conditions #------------------------------------------------------------------------- slope = 0.05 def elevation_tilt(x, y): return x*slope domain.set_quantity('stage', elevation_tilt) domain.set_quantity('friction', friction) domain.set_quantity('elevation',elevation_tilt) print 'Available boundary tags', domain.get_boundary_tags() domain.set_region('dam','stage',0.20, location = 'unique vertices') #domain.set_region(Set_region('dam','stage',0.2, # location = 'unique vertices')) Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([0,0,0]) # to drain the water out. domain.set_boundary( {'wall': Br, 'edge': Bd} ) #------------------------------------------------------------------------- # Evolve system through time #------------------------------------------------------------------------- t0 = time.time() for t in domain.evolve(yieldstep = 0.01, finaltime = 30): domain.write_time() print 'That took %.2f seconds' %(time.time()-t0) print 'finished' points = [[gate_position - 0.65,0.2], [gate_position - 0.55,0.2], [gate_position - 0.45,0.2], [gate_position - 0.35,0.2], [gate_position - 0.25,0.2] ] #------------------------------------------------------------------------- # Calculate gauge info #------------------------------------------------------------------------- interpolate_sww2csv(project.outputtimedir + project.basename +".sww", points, project.depth_filename, project.velocity_x_filename, project.velocity_y_filename) #------------------------------------------------------------- if __name__ == "__main__": main(0.000)