"""Simple water flow example using ANUGA Water flowing down a channel with more complex topography """ #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross from anuga.shallow_water import Domain from anuga.shallow_water import Reflective_boundary from anuga.shallow_water import Dirichlet_boundary from anuga.shallow_water import Time_boundary from anuga.pmesh.mesh_interface import create_mesh_from_regions from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary from math import tan, sqrt, sin, pi #------------------------------------------------------------------------------ # Project file #------------------------------------------------------------------------------ #Set up file structure anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep meshes_dir = anuga_dir+'meshes'+sep meshname = 'circular_mesh.msh' output_dir = anuga_dir+'output'+sep ## #------------------------------------------------------------------------------ ## # Copy scripts to time stamped output directory and capture screen ## # output to file ## #------------------------------------------------------------------------------ ## print "Processor Name:",get_processor_name() ## ## #copy script must be before screen_catcher ## #print kwargs ## ## print 'output_dir', output_dir ## if myid == 0: ## copy_code_files(kwargs['output_dir'],__file__, ## dirname(project.__file__)+sep+ project.__name__+'.py' ) ## ## store_parameters(**kwargs) ## ## barrier() ## ## start_screen_catcher(kwargs['output_dir'], myid, numprocs) ## ## print "Processor Name:",get_processor_name() ## ## # filenames ### meshes_dir_name = project.meshes_dir_name+'.msh' ## ## # creates copy of code in output dir ## print 'min triangles', project.trigs_min, ## print 'Note: This is generally about 20% less than the final amount' #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ length = 30. width = 25. Cx = 12.96 # centre of island on the x axis Cy = 13.8 # centre of island on the y axis dx = dy = 1 # Resolution: Length of subdivisions on both axes water_depth = 0.32 # Can be 0.32 or 0.42 #boundary poly_domain = [[0,0],[length,0],[length,width],[0,width]] # exporting asc grid xmin = 0 xmax = length ymin = 0 ymax = width #Create interior region Dome = [[(Cx)-4,(Cy)-4],[(Cx)+4,(Cy)-4,], [(Cx)+4,(Cy)+4],[(Cx)-4,(Cy)+4]] remainder_res=1 Dome_res = .01 interior_dome = [[Dome, Dome_res]] #Create mesh print 'start create mesh from regions' create_mesh_from_regions(poly_domain, boundary_tags={'wavemaker': [0], 'right': [1], 'top': [2], 'left': [3]}, maximum_triangle_area = remainder_res, filename=meshname, interior_regions = interior_dome, use_cache=False, verbose=False) # Setup computational domain domain = Domain(meshes_dir+sep+meshname, use_cache=False, verbose = True) domain.set_name('circular') # Output name print 'memory usage before del domain',mem_usage() print domain.statistics() print 'triangles',len(domain) #------------------------------------------------------------------------------ # Setup initial conditions #------------------------------------------------------------------------------ def topography(x,y): """Complex topography defined by a function of vectors x and y """ z= 0*x # defining z for all values other than the if statements r= 3.6 # radius, provided in document angle = 14 # angle, provided in document h= r*tan(angle/57.2957795) # finding height of cone if not truncated N = len(x) for i in range(N): #truncated top if (x[i]-Cx)**2 + (y[i]-Cy)**2 <1.1**2: z[i] += 0.625 # cone if (x[i]-Cx)**2 + (y[i]-Cy)**2 1.1**2: z[i] = -(sqrt(((x[i]-Cx)**2+(y[i]-Cy)**2)/((r/h)**2))-h) return z domain.set_quantity('elevation', topography, verbose=True) # Use function for elevation domain.set_quantity('friction', 0.01) # Constant friction domain.set_quantity('stage',water_depth) # Dry initial condition #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ # Create boundary function from timeseries provided in file ##boundary_filename="ts2cnew1_input_20_80sec_new" ##prepare_timeboundary(boundary_filename+'.txt') ## ##function = file_function(boundary_filename+'.tms', ## domain, verbose=True) def wave_form(t): return 0.1*sin(2*pi*t/50.) # Create and assign boundary objects Bw = Dirichlet_boundary([water_depth, 0, 0]) #wall Bt = Transmissive_Momentum_Set_Stage_boundary(domain, wave_form) #wavemaker domain.set_boundary({'left': Bw, 'right': Bw, 'top': Bw, 'wavemaker': Bt}) #------------------------------------------------------------------------------ # Evolve system through time #------------------------------------------------------------------------------ for t in domain.evolve(yieldstep = 0.2, finaltime = 100.0): print domain.timestepping_statistics()