Changeset 3530

Timestamp:

Aug 25, 2006, 1:08:54 PM (18 years ago)

Author:

ole

Message:

Merged in version developed with Rosh

File:

• anuga_core/source/anuga/examples/runup_convergence.py (modified) (11 diffs)

anuga_core/source/anuga/examples/runup_convergence.py

@@ -29 +29 @@
 # Model constants
 
-slope = -0.02     # 1:50 Slope
-highest_point = 3 # Highest elevation (m)
-sea_level = -1    # Mean sea level
-amplitude = 2     # Solitary wave height H
+slope = -0.05       # 1:20 Slope
+highest_point = 6   # Highest elevation (m)
+sea_level = 0       # Mean sea level
+min_elevation = -20 # Lowest elevation (elevation of offshore flat part)
+offshore_depth=sea_level-min_elevation # offshore water depth
+amplitude = 2       # Solitary wave height H
+normalized_amplitude = amplitude/offshore_depth 
 simulation_name = 'runup_convergence'   
 
@@ -61 +64 @@
 
 # Unstructured mesh
-#polygon = [[east,north],[west,north],[west,south],[east,south]]
-#meshname = simulation_name + '.msh'
-#create_mesh_from_regions(polygon,
-#                         boundary_tags={'top': [0], 'left': [1], 'bottom': [2], 'right': [3]},
-#                         maximum_triangle_area=dx*dy/2,
-#                         filename=meshname)
-#                         #interior_regions=interior_regions)
-#domain = Domain(meshname, use_cache=True, verbose = True)
+polygon = [[east,north],[west,north],[west,south],[east,south]]
+interior_polygon = [[200,north-10],[west+10,north-10],[west+10,south+10],[200,south+10]]
+meshname = simulation_name + '.msh'
+create_mesh_from_regions(polygon,
+                         boundary_tags={'top': [0], 'left': [1], 'bottom': [2], 'right': [3]},
+                         maximum_triangle_area=dx*dy/4, # Triangle area commensurate with structured mesh
+                         filename=meshname,
+                         interior_regions=[[interior_polygon,dx*dy/32]])
+domain = Domain(meshname, use_cache=True, verbose = True)
 
 
@@ -79 +83 @@
 #------------------------------------------------------------------------------
 
+#def topography(x,y):
+#    return slope*x+highest_point  # Return linear bed slope bathymetry as vector
+
 def topography(x,y):
-    return slope*x+highest_point  # Return linear bed slope bathymetry as vector
+    """Two part topography - slope and flat part
+    """
+
+    from Numeric import zeros, Float
+
+    z = zeros(len(x), Float) # Allocate space for return vector
+    for i in range(len(x)):
+
+        z[i] = slope*x[i]+highest_point # Linear bed slope bathymetry
+
+        if z[i] < min_elevation:        # Limit depth
+            z[i] = min_elevation
+
+    return z       
+        
+
+        
 
 domain.set_quantity('elevation', topography) # Use function for elevation
@@ -96 +119 @@
 
 def waveform(t): 
-    return sea_level + amplitude/cosh(t-10)**2
+    return sea_level + normalized_amplitude/cosh(t-25)**2
 
 Bw = Time_boundary(domain=domain,     # Time dependent boundary  
@@ -114 +137 @@
 
 stagestep = []
-for t in domain.evolve(yieldstep = 1, finaltime = 150):
+for t in domain.evolve(yieldstep = 1, finaltime = 100):
     domain.write_time()
 
@@ -126 +149 @@
 def gauge_line(west,east,north,south):
     from Numeric import arange
-    x_vector = arange(west, east, dx/4) # Gauges every dx/2 meter from west to east
+    x_vector = arange(west, (east-west)/2, 10) # Gauges every 10 meter from west to midpt
     y = (north+south)/2.
 
@@ -142 +165 @@
                   interpolation_points = gauges,
                   verbose = True,
-                  use_cache = False)
+                  use_cache = True)
 
 
@@ -172 +195 @@
 # Print
 print 'wave height [m]:                    ', amplitude
-runup_height = topography(runup_point, (north+south)/2.)
+runup_height = topography([runup_point], [(north+south)/2.])[0]
 print 'run up height [m]:                  ', runup_height 
 
@@ -179 +202 @@
 
 print 'Coastline (meters form west):       ', coastline
+
 
 
@@ -203 +227 @@
          
 savefig('snapshots')
-        
 
 
@@ -218 +241 @@
 fid.close()
 
-
-
-
-
 # Plot max runup etc
 ion()
-hold(False)
-figure(2)
+figure(1)
 plot(x_vector, max_stage, 'g+',
      x_vector, min_stage, 'b+',