Changeset 5066 for anuga_work/development/MISG_2008

Timestamp:

Feb 20, 2008, 4:01:00 PM (17 years ago)

Author:

sexton

Message:

 

File:

• anuga_work/development/MISG_2008/sensitivity_study.py (modified) (5 diffs)

anuga_work/development/MISG_2008/sensitivity_study.py

@@ -3 +3 @@
 Water driven up a linear slope and time varying boundary,
 similar to a beach environment
+
+Test case for MISG 2008, University of Wollongong
 """
 
@@ -16 +18 @@
 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
@@ -25 +30 @@
 # Setup computational domain
 #------------------------------------------------------------------------------
-
-points, vertices, boundary = rectangular_cross(90, 30,
-                                               len1=3.0,
-                                               len2=1.0) # Basic mesh
-
-domain = Domain(points, vertices, boundary) # Create 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)
 
 #------------------------------------------------------------------------------
@@ -44 +63 @@
     N = len(x)
     z = zeros(N, Float)
-
-
+    
     for i in range(N):
-        # linear bed slope            
-        z[i] = -x[i]/15
-
-        # IID noise
-        z[i] += normal(0.0, 0.01)
+        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
-
-domain.set_quantity('elevation', topography) # Use function for elevation
-domain.set_quantity('friction', 0.01)        # Constant friction 
-domain.set_quantity('stage', -.4)            # Constant negative initial stage
-
 
 #------------------------------------------------------------------------------
@@ -69 +86 @@
 Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
 Bw = Time_boundary(domain=domain,     # Time dependent boundary  
-                   f=lambda t: [(.1*sin(t*2*pi)-0.1) , 0.0, 0.0])
-
-# Associate boundary tags with boundary objects
-domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
-
-
-import time
-t0 = time.time()
-#------------------------------------------------------------------------------
-# Evolve system through time
-#------------------------------------------------------------------------------
-
-for t in domain.evolve(yieldstep = 0.1, finaltime = 50.0):
-    domain.write_time()
-
-print 'That took %.2f seconds' %(time.time()-t0)
-
+                   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: [(0<t<pi/w)*A*sin(w*t), 0.0, 0.0])
+    domain.set_boundary({'left': Br, 'right': Btd, 'top': Br, 'bottom': Br})
+
+    # locations to record wave height and velocity
+    tri_id0 = domain.get_triangle_containing_point([-c/m,y_extent*0.5])
+    tri_id1 = domain.get_triangle_containing_point([-c/m-10,y_extent*0.5])
+    tri_id2 = domain.get_triangle_containing_point([-c/m-20,y_extent*0.5])
+    tri_id3 = domain.get_triangle_containing_point([-c/m+50,y_extent*0.5])
+      
+    stage_centroid     = domain.quantities['stage'].centroid_values
+    xmom_centroid      = domain.quantities['xmomentum'].centroid_values
+    elevation_centroid = domain.quantities['elevation'].centroid_values
+    
+    elev0 = elevation_centroid[tri_id0]
+    elev1 = elevation_centroid[tri_id1]
+    elev2 = elevation_centroid[tri_id2]
+    elev3 = elevation_centroid[tri_id3]
+
+    import time
+    t0 = time.time()
+    #------------------------------------------------------------------------------
+    # Evolve system through time
+    #------------------------------------------------------------------------------
+    domain.time=0.0
+    finaltime=4.0 #pi/w+150.0
+    max=0
+    for t in domain.evolve(yieldstep = 1.0, finaltime = finaltime):
+
+        if  stage_centroid[tri_id0]-elev0>max0[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')