Changeset 3247

Timestamp:

Jun 27, 2006, 6:10:14 PM (18 years ago)

Author:

sexton

Message:

updates for Year 10 student problem

Location:

documentation/user_manual/examples

Files:

• harbour.py (modified) (5 diffs)
• runuptest.py (modified) (5 diffs)

documentation/user_manual/examples/harbour.py

@@ -51 +51 @@
 
 domain = Domain(meshname,use_cache=True,verbose = True)
-domain.set_name('test')                     # Output to test.sww
+domain.set_name('harbour')                  # Output to harbour.sww
 domain.set_datadir('.')                     # Use current directory for output
 domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities
@@ -76 +76 @@
                 if yi < ycentredown: z[i,j] = 0.0
                 if yi > ycentreup: z[i,j] = 0.0
-                if yi > ycentredown and yi < ycentreup: z[i,j] = 
-harbour_depth
-           if xi > centre-harbour_length and xi < centre: z[i,j] = 
-harbour_depth
+                if yi > ycentredown and yi < ycentreup: z[i,j] = harbour_depth
+           if xi > centre-harbour_length and xi < centre: z[i,j] = harbour_depth
     return z
 
@@ -108 +106 @@
     domain.write_time()
 
+#-----------------------------------------------------------------------------
+# Interrogate solution
+#-----------------------------------------------------------------------------
+
 # Gauge line
 def gauge_line(west,east,north,south):
@@ -118 +120 @@
     return gauges
 
-#-----------------------------------------------------------------------------
-# Interrogate solution
-#-----------------------------------------------------------------------------
-
 gauges = gauge_line(west,east,north,south)
 
-f = file_function(gauges)
+from pyvolution.util import file_function
+
+f = file_function('harbour.sww',
+                  quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
+                  interpolation_points = gauges,
+                  verbose = True,
+                  use_cache = True)
+
 maxw = 0.0
 themaxw = maxw
@@ -130 +135 @@
 for k, g in enumerate(gauges):
     for i, t in enumerate(f.get_time()):
-         if w > maxw: maxw = w
+        w = f(t, point_id = k)[0]
+        elevation = f(t, point_id = k)[1]
+        if w > maxw: maxw = w
     maxws.append(maxw) 
-    if themaxw < maxw: 
-       themaxw = maxw
-       gaugeloc = k
+    if themaxw < maxw:
+        themaxw = maxw
+        gaugeloc = k
 
-from pylab import plot, xlabel, ylabel, title
+from pylab import plot, xlabel, ylabel, title, ion, close, savefig
+ion()
 plot(gaugex,maxws,'g-',gaugestar[0],gaugestar[1],'r*')
 xlabel('x')
 ylabel('stage')
 title('Maximum stage for gauge line')
+savefig('max_stage_harbour')
+close('all')

documentation/user_manual/examples/runuptest.py

@@ -17 +17 @@
 from pyvolution.shallow_water import Transmissive_boundary
 
-waveheight = 1.0
+
 #------------------------------------------------------------------------------
 # Setup computational domain
 #------------------------------------------------------------------------------
 
-west = 340000.
-east = 380000.
-south = 6265000.
-north = 6250000.
+waveheight = 1000.0
+depth_east_edge = -4000.
+west = 0.
+east = 80000.
+south = 0.
+north = 5000.
 polygon = [[east,north],[west,north],[west,south],[east,south]]
 meshname = 'test.msh'
@@ -34 +36 @@
                                         'bottom': [2],
                                         'right': [3]},
-                         maximum_triangle_area=50000,
+                         maximum_triangle_area=100000,
                          filename=meshname)
                          #interior_regions=interior_regions)
@@ -51 +53 @@
 
 def topography(x,y):
-    slope = -1./400.
-    c = 850.0
-    return slope*x+c                              # linear bed slope
+    slope = depth_east_edge/((east-west)/2.)
+    c = -(west+east)/2.*slope
+    return slope*x+c                         # linear bed slope
 
 domain.set_quantity('elevation', topography) # Use function for elevation
@@ -79 +81 @@
 
 stagestep = []
-for t in domain.evolve(yieldstep = 10, finaltime = 100.0):
+for t in domain.evolve(yieldstep = 10, finaltime = 1000.0):
     domain.write_time()
-
-# Gauge line
-def gauge_line(west,east,north,south):
-    from Numeric import arange
-    gaugex = arange(west,east,1000.)
-    gauges = []
-    for x in gaugex:
-        y = (north+south)/2.
-        gauges.append([x,y])
-    return gauges
 
 #-----------------------------------------------------------------------------
@@ -96 +88 @@
 #-----------------------------------------------------------------------------
 
-gauges = gauge_line(west,east,north,south)
+# Gauge line
+def gauge_line(west,east,north,south):
+    from Numeric import arange
+    gaugex = arange(east,west,-1000.)
+    gauges = []
+    gaugey = []
+    for x in gaugex:
+        y = (north+south)/2.
+        gaugey.append(y)
+        gauges.append([x,y])
+    return gauges, gaugex, gaugey
 
-f = file_function(gauges)
-maxw = 0.0
-themaxw = maxw
-maxws = []
+gauges, gaugex, gaugey = gauge_line(west,east,north,south)
+
+from pyvolution.util import file_function
+
+f = file_function('test.sww',
+                  quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
+                  interpolation_points = gauges,
+                  verbose = True,
+                  use_cache = True)
+
+from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure
+ion()
+
+maxw = []
+minw = []
+maxd = []
+count = 0
 for k, g in enumerate(gauges):
+    stage = []
+    elevation = []
+    depths = []
     for i, t in enumerate(f.get_time()):
-        if elevation > 0:
-            if w > maxw: maxw = w
-    maxws.append(maxw) 
-    if elevation > 0:
-       if themaxw < maxw: 
-           themaxw = maxw
-           gaugeloc = k
+        w = f(t, point_id = k)[0]
+        elev = f(t, point_id = k)[1]
+        depth = w-elev
+        stage.append(w)
+        elevation.append(elev)
+        depths.append(elev)
+        if g[0] < (west+east)/2. and depth <= 0 and count == 0:
+            count += 1
+            posx = g[0]
+            loc = k
+    maxw.append(max(stage))
+    minw.append(min(stage))
+    maxd.append(max(depths))
 
-filename = 'maxrunup'+waveheight+'.csv'
+filename = 'maxrunup'+str(waveheight)+str(depth_east_edge)+'.csv'
 fid = open(filename,'w')    
-s = 'Waveheight,Runup distance,Runup height\n'
+s = 'Depth at eastern boundary,Waveheight,Runup distance,Runup height\n'
 fid.write(s)
-gaugestar = gauges[gaugeloc]
-s = '%.2f, %.2f/n' %(waveheight,gaugestar[0],gaugestar[1])
+runup_distance = posx 
+runup_height = topography(runup_distance,(north+south)/2.)
+print 'run up height:   ', runup_height 
+if runup_distance < (east+west)/2.:
+    runup_distance_from_coast = (east+west)/2.-runup_distance
+else:
+    runup_distance_from_coast = runup_distance-(east+west)/2.
+print 'run up distance from coastline: ', runup_distance_from_coast
+s = '%.2f,%.2f,%.2f,%.2f\n' %(depth_east_edge,waveheight,runup_distance_from_coast,runup_height)
 fid.write(s)
 
-from pylab import plot, xlabel, ylabel, title
-plot(gaugex,maxws,'g-',gaugestar[0],gaugestar[1],'r*')
+figure(1)
+plot(gaugex[:loc],maxw[:loc],'g-')
 xlabel('x')
 ylabel('stage')
 title('Maximum stage for gauge line')
+savefig('max_stage')
+
+close('all')