Changeset 3505

Timestamp:

Aug 17, 2006, 5:56:35 PM (19 years ago)

Author:

ole

Message:

Runup work

Files:

• documentation/user_manual/examples/runup.py (modified) (1 diff)
• documentation/user_manual/examples/runup_convergence.py (modified) (8 diffs)
• inundation/pyvolution/shallow_water.py (modified) (1 diff)

documentation/user_manual/examples/runup.py

@@ -49 +49 @@
 #------------------------------------------------------------------------------
 
-from math import sin, pi
+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: [(.1*sin(t*2*pi)-0.3), 0.0, 0.0])
+                   f=lambda t: [(.1*sin(t*2*pi)-0.3) * exp(-t), 0.0, 0.0])
 
 # Associate boundary tags with boundary objects

documentation/user_manual/examples/runup_convergence.py

@@ -2 +2 @@
 
 Water driven up a linear slope and time varying boundary,
-similar to a beach environment
+similar to a beach environment.
+
+The study area is discretised as a regular triangular grid 100m x 100m
 """
 
@@ -10 +12 @@
 #------------------------------------------------------------------------------
 
+import sys
+
+from pmesh.mesh_interface import create_mesh_from_regions
 from pyvolution.mesh_factory import rectangular_cross
 from pyvolution.shallow_water import Domain
@@ -15 +20 @@
 from pyvolution.shallow_water import Dirichlet_boundary
 from pyvolution.shallow_water import Time_boundary
-from pyvolution.shallow_water import Transmissive_boundary
+from pyvolution.shallow_water import Transmissive_Momentum_Set_Stage_boundary
+from pyvolution.util import file_function
+from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure, axis
+ion()
+
+
+#------------------------------------------------------------------------------
+# Model constants
+
+slope = -0.02     # 1:50 Slope
+highest_point = 3 # Highest elevation (m)
+sea_level = -2    # Mean sea level
+wave_height = 3   # Solitary wave height H
+simulation_name = 'runup_convergence'   
+
+
+# Basin dimensions (m)
+west = 0          # left boundary
+east = 500        # right boundary 
+south = 0         # lower boundary
+north = 100       # upper bourdary
 
 
@@ -22 +47 @@
 #------------------------------------------------------------------------------
 
-N = 10
-waveheight=5
-points, vertices, boundary = rectangular_cross(N, N, len1=100, len2=100) 
+# Structured mesh
+dx = 10           # Resolution: Lenght of subdivisions on x axis (length)
+dy = 10           # Resolution: Lenght of subdivisions on y axis (width)
+
+length = east-west
+width = north-south
+points, vertices, boundary = rectangular_cross(length/dx, width/dy, len1=length, len2=width,
+                                               origin = (west, south)) 
 
 domain = Domain(points, vertices, boundary) # Create domain
-domain.set_name('runup_convergence')        # Output to bedslope.sww
+
+
+
+# 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)
+
+
+
+domain.set_name(simulation_name)
 
 
@@ -35 +80 @@
 
 def topography(x,y):
-    slope = 0.02
-    elevation  = 10
-    return slope*x+elevation                      # linear bed slope
+    return slope*x+highest_point  # Return linear bed slope bathymetry as vector
 
 domain.set_quantity('elevation', topography) # Use function for elevation
-domain.set_quantity('friction', 0. )         # Constant friction 
-domain.set_quantity('stage', 0.0)            # Constant initial stage
+domain.set_quantity('friction', 0.0 )        # Constant friction 
+domain.set_quantity('stage', sea_level)      # Constant initial stage
+
 
 #------------------------------------------------------------------------------
@@ -47 +91 @@
 #------------------------------------------------------------------------------
 
-from math import sin, pi
+from math import sin, pi, cosh, sqrt
 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
+
+def waveform(t): 
+    return wave_height/cosh(t-5)**2
+
 Bw = Time_boundary(domain=domain,     # Time dependent boundary  
-                   f=lambda t: [((10<t<20)*waveheight), 0.0, 0.0])
+                   f=lambda t: [waveform(t), 0.0, 0.0])
+
+# Time dependent boundary for stage, where momentum is set automatically
+Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
+
 
 # Associate boundary tags with boundary objects
-domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Br, 'right': Bts, 'top': Br, 'bottom': Br})
 
 
@@ -63 +114 @@
 
 stagestep = []
-for t in domain.evolve(yieldstep = 10, finaltime = 100):
+for t in domain.evolve(yieldstep = 1, finaltime = 100):
     domain.write_time()
+
 
 #-----------------------------------------------------------------------------
@@ -70 +122 @@
 #-----------------------------------------------------------------------------
 
+
+# Define line of gauges through center of domain
+def gauge_line(west,east,north,south):
+    from Numeric import arange
+    x_vector = arange(west, (east-west)/2, 10) # Gauges every 1 meter from west to midpt
+    y = (north+south)/2.
+
+    gauges = []
+    for x in x_vector:
+        gauges.append([x,y])
+        
+    return gauges, x_vector
+
+gauges, x_vector = gauge_line(west,east,north,south)
+
+# Obtain interpolated timeseries at gauges
+f = file_function(domain.get_name()+'.sww',
+                  quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
+                  interpolation_points = gauges,
+                  verbose = True,
+                  use_cache = True)
+
+
+# Find runup distance from western boundary through a linear search
+max_stage = []
+min_stage = []
+runup_point = west
+coastline = east        
+for k, g in enumerate(gauges):
+    z = f(0, point_id=k)[1] # Elevation
+
+    min_w = sys.maxint
+    max_w = -min_w
+    for i, t in enumerate(f.get_time()):
+        w = f(t, point_id = k)[0]
+        if w > max_w: max_w = w
+        if w < min_w: min_w = w        
+
+    if max_w-z <= 0.01:  # Find first gauge where max depth > eps (runup)
+        runup_point = g[0]
+
+    if min_w-z <= 0.01:  # Find first gauge where min depth > eps (coastline)
+        coastline = g[0]        
+        
+    max_stage.append(max_w)
+    min_stage.append(min_w)    
+
+
+# Print
+
+print 'wave height [m]:                    ', wave_height
+runup_height = topography(runup_point, (north+south)/2.)
+print 'run up height [m]:                  ', runup_height 
+
+runup_distance = runup_point-coastline
+print 'run up distance from coastline [m]: ', runup_distance
+
+print 'Coastline (meters form west):       ', coastline
+
+
+
+# Store
+filename = 'maxrunup'+str(wave_height)+'.csv'
+fid = open(filename,'w')    
+s = 'Depth at eastern boundary,Waveheight,Runup distance,Runup height\n'
+fid.write(s)
+
+s = '%.2f,%.2f,%.2f\n' %(wave_height, runup_distance, runup_height)
+fid.write(s)
+
+fid.close()
+
+# Plot
+figure(1)
+plot(x_vector, max_stage, 'g+',
+     x_vector, min_stage, 'b+',     
+     x_vector, topography(x_vector,(north+south)/2.), 'r-')
+xlabel('x')
+ylabel('stage')
+title('Maximum stage for gauge line')
+#axis([33000, 47000, -1000, 3000])
+savefig('max_stage')
+
+close('all')
     
 

inundation/pyvolution/shallow_water.py

@@ -1059 +1059 @@
         q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
         value = self.function(self.domain.time)
-        q[0] = value[0]
+
+   
+        try:
+            # In case function returns a list of values
+            val = value[0]
+        except TypeError:
+            val = value
+
+        # Assign value to first component of q (stage)
+        try:
+            q[0] = float(val)
+        except:
+            msg = 'Supplied value could not be converted to float: val=', val     
+            raise Exception, msg    
+            
+            
+            
+        
         return q