"""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 from __future__ import division #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ from anuga.pmesh.mesh_interface import create_mesh_from_regions name ='steep_smflat_gap' crest =[50, 450] A = [0.5, 2, -0.5, -2] gap = [10000] E = len (crest) for i in range(E): M = len (gap) for k in range (M): B = len(A) for l in range(B): length = (crest[i]+2500) width = (gap[k]*3) crestdepth = (-3) umask(002) time = strftime('%Y%m%d_%H%M%S',gmtime()) ## output_dir = 'c:'+sep+'anuga_data'+sep output_dir = sep+'d'+sep+'xrd'+sep+'gem'+sep+'5'+sep+'nhi'+sep+'inundation'+sep+'data'+sep+'idealised_bathymetry_study'+sep+'final_models'+sep+'gap'+sep+str(length)+'_'+str(A[l])+'_'+str(gap[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],[780+crest[i],0],[780+crest[i],width],[140,width]] interior_polygon2 = [[810+crest[i],0],[2499+crest[i],0],[2499+crest[i],width],[810+crest[i],width]] interior_polygon3 = [[700,0], [700,gap[k]], [700+crest[i],gap[k]], [700+crest[i],0]] interior_polygon4 = [[700,2*gap[k]], [700,width], [700+crest[i],width], [700+crest[i],2*gap[k]]] ## interior_regions = [[interior_polygon, 400], [interior_polygon2, 2500], [interior_polygon3, 400], [interior_polygon4, 400]] meshname = str(name)+'.msh' create_mesh_from_regions(boundary_polygon, boundary_tags={'bottom': [0], 'right': [1], 'top': [2], 'left': [3]}, maximum_triangle_area=1250, 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 = (-81/(820+crest[i]-180)*x)+crestdepth-7+(180*(81/(820+crest[i]-180))) N = len (x) for j in range(N): ##Lagoon if x[j] < 180: z[j] = -7+crestdepth ##gap elif 179 < x[j] < (821+crest[i]) and (gap[k]) < y[j] < (2*gap[k]): z[j] = (-81/(820+crest[i]-180)*x[j])+crestdepth-7+(180*(81/(820+crest[i]-180))) ## ##back of gap ## elif 180 < x[j] < 200 and (gap[k]) < y[j] < (2*gap[k]): ## z[j] = -x[j]+173+crestdepth ## ##Side of Gap Right ## elif 179< x[j] < (821+crest[i]) and (gap[k]*2)-4 < y[j] < (2*gap[k])+1: ## z[j] = x[j]*20+crestdepth ## #### ##Side of Gap left ## elif 179 < x[j] < (821+crest[i]) and ((gap[k])-1) < y[j] < (gap[k]+4): ## z[j] = -x[j]*20+crestdepth #crest to lagoon elif 179 < x[j] < 200: z[j] = 0.2*x[j]-40+crestdepth #Flat elif 199 < x[j] < 700: z[j] = -0.5+crestdepth ##Crest elif 699 < x[j] < (700+crest[i]): z[j] = +crestdepth #Curve down elif (699+crest[i]) < x[j] < (720+crest[i]): z[j] = -0.01*(x[j]-(699+crest[i]))*(x[j]-(699+crest[i]))+crestdepth ##steep slope elif (719+crest[i]) < x[j] < (820+crest[i]): z[j] =(-0.84*x[j])+(0.84*(720+crest[i]))+crestdepth-4 return z domain.set_quantity('elevation', topography, alpha=0.1) # Use function for elevation ## domain.set_quantity('friction', 0) # Constant friction domain.set_quantity('friction', Polygon_function( [(boundary_polygon, 0.05),(interior_polygon, 0.05), (interior_polygon2, 0.05), (interior_polygon3, 0.2), (interior_polygon4, 0.2)] ) )#changing friction over two polygons domain.set_quantity('stage', 0.) # Constant negative initial stage domain.tight_slope_limiters = 1 #------------------------------------------------------------------------------ # 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 = 288 H_d_ratio = 0.008 Xo = 71700 po = 0.036 def waveform(t): return A[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+'gap'+sep+'gauges_smflat_gap.csv', production_dirs, report = False, reportname = '', plot_quantity = ['stage','speed','bearing'], 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')