"""Simple water flow example using ANUGA Water driven up a linear slope and time varying boundary, similar to a beach environment Test case for MISG 2008, University of Wollongong """ #------------------------------------------------------------------------------ # 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.pmesh.mesh_interface import create_mesh_from_regions from Numeric import zeros, Float from RandomArray import normal, seed seed(13, 17) # Ensure random number sequence is reproducible #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ y_extent=2000.0 bounding_polygon = [[0,0],[2000,0],[2000,2000],[0,2000]] beach_polygon = [[5,900],[225,900],[225,1100],[5,1100]] interior_regions = [[beach_polygon, 50]] create_mesh_from_regions(bounding_polygon, boundary_tags={'bottom':[0], 'right':[1], 'top':[2], 'left':[3]}, maximum_triangle_area=2000, interior_regions= interior_regions, filename='misg.msh', use_cache=False, verbose=True) domain = Domain('misg.msh',use_cache=False, verbose=True) domain.set_name('sensitivity') # Output to file runup.sww domain.set_datadir('.') # Use current directory for output domain.tight_slope_limiters = 1 domain.beta_h = 0 domain.set_minimum_storable_height(0.01) print 'number of elements: ', len(domain) #------------------------------------------------------------------------------ # Setup initial conditions #------------------------------------------------------------------------------ def topography(x,y): N = len(x) z = zeros(N, Float) for i in range(N): if x[i] <= (depth-c)/m: z[i] = m*x[i]+c else: # linear bed slope #z[i] = slope*x[i] + depth-(slope)*(depth-c)/m z[i] = slope*x[i] + depth-(slope)*(depth-c)/m + 1.0*sin((x[i]-(depth-c)/m)/T) + 1.0*sin((y[i]-(depth-c)/m)/T) # IID noise #z[i] += normal(-100.0, 5.0) return z #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ from math import sin, pi, exp Br = Reflective_boundary(domain) # Solid reflective wall Bt = Transmissive_boundary(domain) # Continue all values on boundary Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values Bw = Time_boundary(domain=domain, # Time dependent boundary f=lambda t: [0.0, 0.0, 0.0]) from Numeric import zeros, Float, arange start=-1.0/40.0 finish=-1.8/40.0 n_vec=arange(start,finish,-2.0/400.0) N=len(n_vec) #N=1 # use for 2.5% and 5% change tests max0=zeros(N, Float) max1=zeros(N, Float) max2=zeros(N, Float) max3=zeros(N, Float) mmax0=zeros(N, Float) mmax1=zeros(N, Float) mmax2=zeros(N, Float) mmax3=zeros(N, Float) output_file ='misg_influence_data.csv' domain.store=True fid_out = open(output_file, 'w') s = 'variable, h0, h1, h2, h3, v0, v1, v2, v3\n' fid_out.write(s) for i in range(N): manning = 0.01*1.0 A=1.0*1.0 w=1.0/100.0 m=-1.0/20.0 depth=-10.0 c=10.0 T=100.0*1.0 #slope=-1.0/40.0*1.0 slope = n_vec[i] var = slope domain.set_quantity('elevation', topography) # Use function for elevation domain.set_quantity('friction', manning) domain.set_quantity('stage', 0.0) # Constant negative initial stage domain.set_quantity('xmomentum', 0.0) # Use function for elevation domain.set_quantity('ymomentum', 0.0) # Use function for elevation swwname = 'sensitivity_2d_wavey_bottom'#+str(i) domain.set_name(swwname) Btd = Transmissive_Momentum_Set_Stage_boundary(domain=domain, function=lambda t: [(0max0[i]: max0[i]=stage_centroid[tri_id0]-elev0 if stage_centroid[tri_id1]-elev1>max1[i]: max1[i]=stage_centroid[tri_id1]-elev1 if stage_centroid[tri_id2]-elev2>max2[i]: max2[i]=stage_centroid[tri_id2]-elev2 if stage_centroid[tri_id3]-elev3>max3[i]: max3[i]=stage_centroid[tri_id3]-elev3 if abs(xmom_centroid[tri_id0])>mmax0[i]: mmax0[i]=abs(xmom_centroid[tri_id0]) if abs(xmom_centroid[tri_id1])>mmax1[i]: mmax1[i]=abs(xmom_centroid[tri_id1]) if abs(xmom_centroid[tri_id2])>mmax2[i]: mmax2[i]=abs(xmom_centroid[tri_id2]) if abs(xmom_centroid[tri_id3])>mmax3[i]: mmax3[i]=abs(xmom_centroid[tri_id3]) print 'finished cycle: ', i print 'That took %.2f seconds' %(time.time()-t0) s = '%.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f\n' %(var, max0[i], max1[i], max2[i], max3[i], mmax0[i], mmax1[i], mmax2[i], mmax3[i]) fid_out.write(s) from pylab import plot, ion, savefig, xlabel, ylabel, title, close, legend,hist,subplot ion() plot(n_vec,max2,'bo-',n_vec,max1,'go-',n_vec,max0,'ro-',n_vec,max3,'co-') title('MISG testing: varying slope') xlabel('slope') ylabel('depth (m)') #legend(['loc2','loc1','loc0','loc3']) name='bed_slope' savefig(name) ''' #title('MISG testing: varying slope PDF') subplot(2,2,1) [h0,bins,patches]=hist(max0,10) subplot(2,2,2) [h1,bins,pa3ches]=hist(max1,10) subplot(2,2,3) [h2,bins,patches]=hist(max2,10) subplot(2,2,4) [h3,bins,patches]=hist(max3,10) #plot(bins,h0,'bo-',bins,h1,'ro-',bins,h2,'go-',bins,h3,'co-') xlabel('Runup') ylabel('Frequency') name='bed_slope_PDF' savefig(name) ''' close('all')