Changeset 7579

Timestamp:

Dec 7, 2009, 6:18:23 PM (15 years ago)

Author:

ole

Message:

Cleanup of school exercise. Optimisation to follow

File:

• anuga_work/development/classroom/WORKING-RipCurrent.py (modified) (8 diffs)

anuga_work/development/classroom/WORKING-RipCurrent.py

@@ -5 +5 @@
 # Import necessary modules
 #------------------------------------------------------------------------------
-from anuga.abstract_2d_finite_volumes.mesh_factory import *
-from anuga.shallow_water import *
-from anuga.shallow_water.shallow_water_domain import *
-from math import *
-from numpy import *
+from anuga.interface import create_domain_from_regions
+from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+from anuga.shallow_water.shallow_water_domain import Time_boundary
+from pylab import figure, quiver, show, cos, sin, pi
 import numpy
-from matplotlib import *
-from pylab import *
 import csv
-
-#Parameters
-
+import time
+
+
+#------------------------------------------------------------------------------
+# Parameters
+#------------------------------------------------------------------------------
 filename = "WORKING-RIP-LAB_Expt-Geometry_Triangular_Mesh"
-location_of_shore = 140    #the position along the y axis of the shorefront
-sandbar = 1.2 #height of sandbar
-sealevel = 0 #height of coast above sea level
-steepness = 8000 #period of sandbar - larger number gives smoother slope - longer period 
+location_of_shore = 140 # The position along the y axis of the shorefront
+sandbar = 1.2           # Height of sandbar
+sealevel = 0            # Height of coast above sea level
+steepness = 8000        # Period of sandbar - 
+                        # larger number gives smoother slope - longer period 
 halfchannelwidth = 5
 bank_slope = 0.1
@@ -27 +29 @@
 timestep = 1
 
+
 #------------------------------------------------------------------------------
 # Setup computational domain
@@ -32 +35 @@
 length = 120
 width = 170
-seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh for seafloor
+seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh
 feature_resolution = 1.0
 beach_resolution = 10.0
@@ -42 +45 @@
 meshname = str(filename)+'.msh'
 
-#interior regions
+# Interior regions
 feature_regions = [[feature_boundary_polygon, feature_resolution],
                    [beach_interior_polygon, beach_resolution]]
 
-create_mesh_from_regions(sea_boundary_polygon,
-                         boundary_tags={'bottom': [0],
-                                        'right' : [1],
-                                        'top' :   [2],
-                                        'left':   [3]},
-                         maximum_triangle_area = seafloor_resolution,
-                         filename = meshname,
-                         interior_regions = feature_regions,
-                         use_cache = False,
-                         verbose = False)
-domain = Domain(meshname, use_cache=True, verbose=True)
+domain = create_domain_from_regions(sea_boundary_polygon,
+                                    boundary_tags={'bottom': [0],
+                                                   'right' : [1],
+                                                   'top' :   [2],
+                                                   'left':   [3]},
+                                    maximum_triangle_area = seafloor_resolution,
+                                    mesh_filename = meshname,
+                                    interior_regions = feature_regions,
+                                    use_cache = True,
+                                    verbose = True)
+                                    
 domain.set_name(filename)                  # Output name
 print domain.statistics()
 
+
 #------------------------------------------------------------------------------
 # Setup initial conditions
 #------------------------------------------------------------------------------
-
 def topography(x,y):
     """Complex topography defined by a function of vectors x and y."""
 
-    #general slope and buildings
+    # General slope and buildings
     z=0.05*(y-(location_of_shore))
     
@@ -76 +79 @@
     for i in range(N):
         if y[i]>150:
-            z[i] = (0.1*(y[i]-150)) + 0.05*(y[i] - (location_of_shore))
+            z[i] = (0.1*(y[i]-150)) + 0.05*(y[i]-(location_of_shore))
 
     return z
@@ -85 +88 @@
     
     N = len(x)
-    for i in range(N):
-        if -1*(bank_slope)*x[i] + 112 < y[i] < -1*(bank_slope)*x[i] + 124 and 0<x[i]<((length/2)-halfchannelwidth):
+    
+    # It would be great with a comment about what this does
+    for i in range(N):
+        ymin = -1*(bank_slope)*x[i] + 112 
+        ymax = -1*(bank_slope)*x[i] + 124
+        xmin = 0
+        xmax = (length/2)-halfchannelwidth
+        if ymin < y[i] < ymax and xmin < x[i]< xmax:
             z[i] += sandbar*((cos((y[i]-118)/steepness)))
-    for i in range(N):
-        if -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+112) < y[i] < -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+124) and ((length/2)+halfchannelwidth)<x[i]<183:
+    
+    # It would be great with a comment about what this does        
+    for i in range(N):
+        ymin = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+112)
+        ymax = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+124)
+        xmin = (length/2)+halfchannelwidth
+        xmax = 183
+        if ymin < y[i] < ymax and xmin < x[i] < xmax:
             z[i] += sandbar*(cos((y[i]-118)/steepness))
+            
     return z
 
-domain.set_quantity('elevation', topography)           # elevation is a function
-domain.set_quantity('friction', 0.01)                  # Constant friction
-
-domain.add_quantity('elevation', topography3)
-
-domain.set_quantity('stage', 0)   # Dry initial condition
+domain.set_quantity('elevation', topography)  # Apply base elevation function
+domain.add_quantity('elevation', topography3) # Add elevation modification
+domain.set_quantity('friction', 0.01)         # Constant friction
+domain.set_quantity('stage', 0)               # Constant initial condition at mean sea level
+
 
 #------------------------------------------------------------------------------
@@ -107 +122 @@
 Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
 
-
-amplitude = 0.4 #amplitude of wave (wave height m)
-period =  5     #wave period (sec)
-
 def wave(t):
-
-    A = amplitude # Amplitude [m] (Wave height)
-    T = period  # Wave period [s]
+    """Define wave driving the system
+    """
+    
+    A = 0.4  # Amplitude of wave [m] (wave height)
+    T = 5    # Wave period [s]
 
     if t < 30000000000:
@@ -126 +139 @@
 domain.set_boundary({'left': Br, 'right': Br, 'top': Bo, 'bottom': Bt})
 
+
 #------------------------------------------------------------------------------
 # Evolve system through time
 #------------------------------------------------------------------------------
 
-list1 = [x for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ]
-
-print list1
-
-u = numpy.zeros(len(list1), 'd')
-v = numpy.zeros(len(list1), 'd')
-
-
-for t in domain.evolve(yieldstep = (timestep), finaltime = (simulation_length)):
+# Read in gauge locations for interpolation and convert them to floats
+gauge_file = open('New_gauges.csv', 'r')
+G = [(float(x[0]), float(x[1])) for x in csv.reader(gauge_file,
+                                                    dialect='excel', 
+                                                    delimiter=',')]
+gauges = numpy.array(G)
+number_of_gauges = len(gauges)
+print gauges                               
+
+# Allocate space for velocity values
+u = numpy.zeros(number_of_gauges)
+v = numpy.zeros(number_of_gauges)
+
+t0 = time.time()
+for t in domain.evolve(yieldstep = timestep, finaltime = simulation_length):
     print domain.timestepping_statistics()
-    S = domain.get_quantity('stage').get_values(interpolation_points=list1)
-    E = domain.get_quantity('elevation').get_values(interpolation_points=list1)
+    S = domain.get_quantity('stage').get_values(interpolation_points=gauges)
+    E = domain.get_quantity('elevation').get_values(interpolation_points=gauges)
     depth = S-E
 
-    uh = domain.get_quantity('xmomentum').get_values(interpolation_points=list1)
-    vh = domain.get_quantity('ymomentum').get_values(interpolation_points=list1)
+    uh = domain.get_quantity('xmomentum').get_values(interpolation_points=gauges)
+    vh = domain.get_quantity('ymomentum').get_values(interpolation_points=gauges)
     u += uh/depth
     v += vh/depth
-
-n_time_intervals = (simulation_length)/(timestep)
-u_average = u / (n_time_intervals)
-v_average = v / (n_time_intervals)
-
-print "there were", n_time_intervals, "time steps"
-
-print "sum y velocity", v
-print "average y velocity", v_average
-print "sum x velocity", u
-print "average x velocity", u_average
-
-x_output = file("x_velocity.txt", 'w')
-y_output = file("y_velocity.txt", 'w')
-
-print >> x_output, " "
-print >> y_output, " "
-
-print >> x_output, u_average
-print >> y_output, v_average
-
-X = [x[0] for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ]
-Y = [x[1] for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ]
+print 'Evolution took %.2f seconds' % (time.time()-t0)
+
+    
+#------------------------------------------------------------------------------
+# Post processing
+#------------------------------------------------------------------------------    
+
+# Use this
+
+#     from anuga.fit_interpolate.interpolate import Interpolation_function
+
+#     # Get mesh and quantities from sww file
+#     X = get_mesh_and_quantities_from_file(filename,
+#                                           quantities=quantity_names,
+#                                           verbose=verbose)
+#     mesh, quantities, time = X
+
+#     # Find all intersections and associated triangles.
+#     segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
+
+#     # Get midpoints
+#     interpolation_points = segment_midpoints(segments)
+
+#     # Interpolate
+#     if verbose:
+#         log.critical('Interpolating - total number of interpolation points = %d'
+#                      % len(interpolation_points))
+
+#     I = Interpolation_function(time,
+#                                quantities,
+#                                quantity_names=quantity_names,
+#                                vertex_coordinates=mesh.nodes,
+#                                triangles=mesh.triangles,
+#                                interpolation_points=interpolation_points,
+#                                verbose=verbose)
+
+#     return segments, I
+
+
+
+
+
+
+n_time_intervals = simulation_length/timestep
+u_average = u/n_time_intervals
+v_average = v/n_time_intervals
+
+#print "there were", n_time_intervals, "time steps"
+
+#print "sum y velocity", v
+#print "average y velocity", v_average
+#print "sum x velocity", u
+#print "average x velocity", u_average
+
+x_output = file('x_velocity.txt', 'w')
+y_output = file('y_velocity.txt', 'w')
+
+#print >> x_output, " "
+#print >> y_output, " "
+
+#print >> x_output, u_average
+#print >> y_output, v_average
+
+
+X = gauges[:,0] 
+Y = gauges[:,1] 
+                              
 U = u_average.tolist()
 V = v_average.tolist()
 
-print "U = ", U
-print "U has type", type(U)
-
-
-from matplotlib import *
-from pylab import *
+#print "U = ", U
+#print "U has type", type(U)
+
 
 figure()