"""Simple water flow example using ANUGA Water driven up a linear slope and time varying boundary, similar to a beach environment """ #------------------------------------------------------------------------------ # 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.shallow_water import Transmissive_boundary from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files from time import strftime, gmtime from os import sep, environ, getenv, getcwd,umask from anuga.utilities.polygon import Polygon_function #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ from anuga.pmesh.mesh_interface import create_mesh_from_regions name ='steep_lgflat' wave = [1, -1] #1 returns leading depression N-wave #-1 returns leading crest N-wave crest =[50, 450] crestdepth = [2, -2, -6] N = len (crest) for i in range(N): M = len (crestdepth) for k in range (M): B = len(wave) for l in range(B): length = (crest[i]+2500) width = 20. A = 1 T = 2700 umask(002) time = strftime('%Y%m%d_%H%M%S',gmtime()) output_dir = sep+'d'+sep+'xrd'+sep+'gem'+sep+'5'+sep+'nhi'+sep+'inundation'+sep+'data'+sep+'idealised_bathymetry_study'+sep+'final_models'+sep+'general'+sep+'steep_lgflat'+sep+str(wave[l])+'_'+str(length)+'_'+str(A)+'_'+str(T)+'_'+str(crestdepth[k])+'_'+str(name)+sep sww_file = str(name) copy_code_files(output_dir,__file__,__file__) start_screen_catcher(output_dir) dx = dy = .5 # Resolution: Length of subdivisions on both axes boundary_polygon = [[0,0],[length,0],[length,width],[0,width]] interior_polygon = [[140,0],[2280+crest[i],0],[2280+crest[i],20],[140,20]] interior_polygon2 =[[2290+crest[i],0],[2499+crest[i],0],[2499+crest[i],20],[2290+crest[i],20]] interior_regions = [[interior_polygon, 8], [interior_polygon2, 50]] meshname = str(name)+'.msh' create_mesh_from_regions(boundary_polygon, boundary_tags={'bottom': [0], 'right': [1], 'top': [2], 'left': [3]}, maximum_triangle_area=50, filename=meshname, interior_regions=interior_regions, use_cache=False, verbose=False) domain = Domain(meshname, use_cache=True, verbose=True) print 'Number of triangles = ', len(domain) print 'The extent is ', domain.get_extent() print domain.statistics() domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) domain.set_minimum_storable_height(0.01) domain.set_default_order(2) # Second order spatial approximation domain.set_name(sww_file) # Output name domain.set_datadir(output_dir) #------------------------------------------------------------------------------ # Setup initial conditions #------------------------------------------------------------------------------ def topography(x,y): """Complex topography defined by a function of vectors x and y """ z =(-0.84*x)+(0.84*(2220+crest[i]))+crestdepth[k]-4 N = len (x) for j in range(N): if x[j] < 180: z[j] = -4+crestdepth[k] elif 179 < x[j] < 200: z[j] = 0.2*x[j]-40+crestdepth[k] #Reef Flat elif 199 < x[j] < 2200: z[j] = -0.5+crestdepth[k] ##Crest elif 2199 < x[j] < (2200+crest[i]): z[j] = +crestdepth[k] ## #Curve down elif (2199+crest[i]) < x[j] < (2220+crest[i]): z[j] = -0.01*(x[j]-(2199+crest[i]))*(x[j]-(2199+crest[i]))+crestdepth[k] return z domain.set_quantity('elevation', topography) # Use function for elevation ## domain.set_quantity('friction', 0) # Constant friction domain.set_quantity('friction', Polygon_function( [(boundary_polygon, 0.05),(interior_polygon ,0.2), (interior_polygon2 , 0.05)] ) )#changing friction over two polygons domain.set_quantity('stage', 0.) # Constant negative initial stage #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ from math import sin, pi, exp, cos, cosh, sqrt Br = Reflective_boundary(domain) # Solid reflective wall Bt = Transmissive_boundary(domain) # Continue all values on boundary Bd = Dirichlet_boundary([0.,0.,0.]) # Constant boundary values Bw = Time_boundary(domain=domain, # Time dependent boundary f=lambda t: [sin(2*pi*(t)/1010), -37, 0.0]) g = 9.81 offshore_depth = 250 H_d_ratio = 0.008 Xo = 67000 po = 0.03 def waveform(t): return wave[l]*offshore_depth*(sqrt(g/offshore_depth)*t-Xo/offshore_depth)*sqrt(H_d_ratio*po)*H_d_ratio/cosh(sqrt(3*H_d_ratio*po/4)*(sqrt(g/offshore_depth)*t-Xo/offshore_depth))/cosh(sqrt(3*H_d_ratio*po/4)*(sqrt(g/offshore_depth)*t-Xo/offshore_depth)) Bf = Transmissive_Momentum_Set_Stage_boundary(domain, waveform) # Associate boundary tags with boundary objects domain.set_boundary({'left': Bd, 'right': Bf, 'top': Br, 'bottom': Br}) #------------------------------------------------------------------------------ # Evolve system through time #------------------------------------------------------------------------------ for t in domain.evolve(yieldstep = 5 , finaltime = length*2): domain.write_time() """ Generate time series of nominated "gauges" Note, this script will only work if pylab is installed on the platform Inputs: production dirs: dictionary of production directories with a association to that simulation run, eg high tide, magnitude, etc. Outputs: * figures stored in same directory as sww file * time series data stored in csv files in same directory as sww file * elevation at nominated gauges (elev_output) """ from os import getcwd, sep, altsep, mkdir, access, F_OK, remove from anuga.abstract_2d_finite_volumes.util import sww2timeseries # nominate directory location of sww file with associated attribute production_dirs = {output_dir: str(name)} # Generate figures swwfiles = {} for label_id in production_dirs.keys(): file_loc = label_id swwfile = file_loc + str(name)+'.sww' swwfiles[swwfile] = label_id print 'hello', swwfile texname, elev_output = sww2timeseries(swwfiles, sep+'d'+sep+'cit'+sep+'1'+sep+'cit'+sep+'natural_hazard_impacts'+sep+'inundation'+sep+'sandpits'+sep+'jbrowning'+sep+'anuga'+sep+'anuga_work'+sep+'development'+sep+'idealised_bathymetry_study'+sep+'final_models'+sep+'general'+sep+'gauges_lgflat.csv', production_dirs, report = False, reportname = '', plot_quantity = ['stage', 'speed'], generate_fig = False, surface = False, time_min = None, time_max = None, #time_unit = 'secs', title_on = True, verbose = True) remove(output_dir+sep+str(name)+'.sww')