"""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, mkdir from anuga.utilities.polygon import Polygon_function from anuga.pmesh.mesh_interface import create_mesh_from_regions #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ crest =[50] crestdepth = [0.5, 0, -0.5, -4, -6] N = len (crest) for i in range(N): M = len (crestdepth) for k in range (M): length = (crest[i]+500) width = 20. A = 1 T = 1800 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+str(length)+'_'+str(A)+'_'+str(T)+'_'+str(crestdepth[k])+'_big_mesh'+sep sww_file = 'reef' 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 = [[170,0],[230+crest[i],0],[230+crest[i],20],[170,20]] interior_polygon2 =[[231+crest[i],0],[499+crest[i],0],[499+crest[i],20],[231+crest[i],20]] ## interior_regions = [[interior_polygon, 1250], [interior_polygon2, 1250]] meshname = 'reef'+'.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=False, 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.026*x)+(0.026*(200+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] ##Crest elif 199 < x[j] < (200+crest[i]): z[j] = +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.),(interior_polygon ,0.05), (interior_polygon2 , 0.)] ) )#changing friction over two polygons domain.set_quantity('stage', 0.) # Constant negative initial stage #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ from math import sin, pi, exp, cos 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]) A = 1 T = 1800 def waveform(t): return A*sin(2*pi*(t)/T) 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 from anuga.abstract_2d_finite_volumes.util import sww2timeseries # nominate directory location of sww file with associated attribute production_dirs = {output_dir: 'reef'} # Generate figures swwfiles = {} for label_id in production_dirs.keys(): file_loc = label_id swwfile = file_loc + 'reef.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+'mesh_test'+sep+'gauges_mesh.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)