"""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 = 'curved_down_slope_2' shelf = [300000] slope = [250000] wave = [0.5, -0.5] #1 returns leading depression N-wave #-1 returns leading crest N-wave N = len (shelf) for i in range(N): M = len (slope) for k in range (M): B = len(wave) for l in range(B): length = (shelf[i]+slope[k]) width = 400. A = 1 T = 2700 umask(002) time = strftime('%Y%m%d_%H%M%S',gmtime()) ## output_dir = 'C:'+sep+'anuga_data'+sep+'topography'+sep+str(name)+sep+str(name)+'_'+str(wave[l])+'_'+str(shelf[i])+'_'+str(slope[k])+sep output_dir = sep+'d'+sep+'sim'+sep+'1'+sep+'mpittard'+sep+'idealised_bathymetry_study'+sep+'topography'+sep+str(name)+sep+str(name)+'_'+str(wave[l])+'_'+str(shelf[i])+'_'+str(slope[k])+sep sww_file = str(name) copy_code_files(output_dir,__file__,__file__) start_screen_catcher(output_dir) boundary_polygon = [[0,0],[length,0],[length,width],[0,width]] meshname = str(name)+'.msh' create_mesh_from_regions(boundary_polygon, boundary_tags={'bottom': [0], 'right': [1], 'top': [2], 'left': [3]}, maximum_triangle_area=20000, filename=meshname, 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) 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 """ o = 2500/(slope[k]*slope[k]/4) print str(2500/(slope[k]*slope[k]/4)) z = o*(x-(shelf[i]+slope[k]))*(x-(shelf[i]+slope[k]))-5125-10 S = len (x) for j in range(S): if x[j] < shelf[i]: z[j] = -125/(shelf[i]*shelf[i])*x[j]*x[j]-10 elif shelf[i] <= x[j] < (shelf[i]+slope[k]*0.5) : z[j] = (-o)*(x[j]-shelf[i])*(x[j]-shelf[i])-125-10 return z domain.set_quantity('elevation', topography) # Use function for elevation domain.set_quantity('friction', 0) # Constant friction domain.set_quantity('stage', 0) # Constant negative initial stage domain.tight_slope_limiters = 1 #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ from math import sin, pi, exp, cos, sqrt, cosh 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 g = 9.81 offshore_depth = 5145 H_d_ratio = 0.0004 Xo = 303000 po = 12 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 = 45, finaltime = -2000+((length/50000)+1)*600+((shelf[i]/25000+1)*1000)+1500): domain.write_time() for t in domain.evolve(yieldstep = 45, finaltime = 2700+((length/50000)+1)*600+((shelf[i]/25000+1)*1000)+1500, skip_initial_step = True): domain.write_time() for t in domain.evolve(yieldstep = 120, finaltime = -10000+(length/25)+2700+((length/50000)+1)*600+((shelf[i]/25000+1)*1000)+1500, skip_initial_step = True): 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+'sim'+sep+'1'+sep+'mpittard'+sep+'anuga'+sep+'anuga_work'+sep+'development'+sep+'idealised_bathymetry_study'+sep+'continental_shelves'+sep+'gauges.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) ## print (output_dir+sep+str(name)+'.sww') ## remove(output_dir+sep+str(name)+'.sww')