Changeset 1657

Timestamp:

Jul 29, 2005, 3:33:32 PM (20 years ago)

Author:

ole

Message:

Comments and work on benchmark 2 (lwru2.py)
New unit test of least squares populating domain

Location:

inundation/ga/storm_surge

Files:

• pmesh/load_mesh/loadASCII.py (modified) (1 diff)
• pyvolution/data_manager.py (modified) (2 diffs)
• pyvolution/general_mesh.py (modified) (1 diff)
• pyvolution/generic_boundary_conditions.py (modified) (2 diffs)
• pyvolution/quantity.py (modified) (1 diff)
• pyvolution/shallow_water.py (modified) (1 diff)
• pyvolution/test_least_squares.py (modified) (2 diffs)
• pyvolution/test_quantity.py (modified) (3 diffs)
• validation/LWRU2/lwru2.py (modified) (3 diffs)

inundation/ga/storm_surge/pmesh/load_mesh/loadASCII.py

@@ -997 +997 @@
             msg = 'Could not open file %s ' %ofile
             raise IOError, msg
-
-
-
-
-
-
-
-
         except IOError:
             # Catch this to add an error message

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -1 @@
-"""Functions to store and retrieve data for the Shallow Water Wave equation.
-There are two kinds of data
-
-  1: Constant data: Vertex coordinates and field values. Stored once
-  2: Variable data: Conserved quantities. Stored once per timestep.
-
-All data is assumed to reside at vertex locations.
-
-
-Formats used within ANUGA:
+"""datamanager.py - input output for AnuGA
+
+
+This module takes care of reading and writing datafiles such as topograhies,
+model output, etc
+
+
+Formats used within AnuGA:
 
 .sww: Netcdf format for storing model output
@@ -209 +207 @@
 #Class for storing output to e.g. visualisation
 class Data_format_sww(Data_format):
-    """Interface to native NetCDF format (.sww)
+    """Interface to native NetCDF format (.sww) for storing model output
+
+    There are two kinds of data
+
+    1: Constant data: Vertex coordinates and field values. Stored once
+    2: Variable data: Conserved quantities. Stored once per timestep.
+
+    All data is assumed to reside at vertex locations.
     """
 

inundation/ga/storm_surge/pyvolution/general_mesh.py

@@ -180 +180 @@
         if obj is True, the x/y pairs are returned in a 3*N x 2 array.
         FIXME, we might make that the default.
+        FIXME Maybe use keyword: continuous?
 
         

inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py

@@ -125 +125 @@
     """
 
+    #FIXME: Is this still necessary
 
     def __init__(self, filename, domain):
@@ -170 +171 @@
 
     The full solution history is not exactly the same as if
-    the models were couple:
+    the models were coupled:
     File boundary must read and interpolate from *smoothed* version
     as stored in sww and cannot work with the discontinuos triangles.

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -410 +410 @@
             assert A.shape[0] == len(indexes)
             vertex_list = indexes
+            
         #Go through list of unique vertices
-        for i_index,unique_vert_id in enumerate(vertex_list):
+        for i_index, unique_vert_id in enumerate(vertex_list):
             triangles = self.domain.vertexlist[unique_vert_id]
 

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -94 +94 @@
 
         #Reduction operation for get_vertex_values
-        #from util import mean
-        #self.reduction = mean
-        self.reduction = min  #Looks better near steep slopes
+        from util import mean
+        self.reduction = mean
+        #self.reduction = min  #Looks better near steep slopes
 
         self.quantities_to_be_stored = ['stage']

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -831 +831 @@
 
         f = interp.fit(z)
-        #f will be different from answerr due to smoothing
+        #f will be different from answer due to smoothing
         assert allclose(f, answer,atol=5)
 
@@ -923 +923 @@
         answer = linear_function(data_points2)
         assert allclose(z1, answer)
+
+
+
+
 
 

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -8 +8 @@
 from config import epsilon
 from Numeric import allclose, array, ones, Float
+
+
+#Aux for least_squares example
+def linear_function(point):
+    point = array(point)
+    return point[:,0]+point[:,1]
 
 
@@ -248 +254 @@
                         [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
 
+
     def test_set_vertex_values_using_general_interface(self):
         quantity = Quantity(self.mesh4)
@@ -266 +273 @@
                                                [2.5, 3.5, 2]])
 
+
+
+
+    def test_set_vertex_values_using_least_squares(self):
+        from least_squares import Interpolation, fit_to_mesh
+        
+        quantity = Quantity(self.mesh4)
+
+        #Get (enough) datapoints
+        data_points = [[ 0.66666667, 0.66666667],
+                       [ 1.33333333, 1.33333333],
+                       [ 2.66666667, 0.66666667],
+                       [ 0.66666667, 2.66666667],
+                       [ 0.0, 1.0],
+                       [ 0.0, 3.0],
+                       [ 1.0, 0.0],
+                       [ 1.0, 1.0],
+                       [ 1.0, 2.0],
+                       [ 1.0, 3.0],
+                       [ 2.0, 1.0],
+                       [ 3.0, 0.0],
+                       [ 3.0, 1.0]]
+
+
+        interp = Interpolation(quantity.domain.coordinates, quantity.domain.triangles,
+                               data_points, alpha=0.0)
+
+        z = linear_function(data_points)
+        answer = linear_function(quantity.domain.coordinates)
+
+        f = interp.fit(z)
+
+        #print "f",f
+        #print "answer",answer
+        assert allclose(f, answer)
+
+
+        quantity.set_values(f)
+
+
+        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
+        #print quantity.vertex_values, answer
+        assert allclose(quantity.vertex_values.flat, answer)
+
+        #Now try using the general interface
+
+        vertex_attributes = fit_to_mesh(quantity.domain.coordinates,
+                                        quantity.domain.triangles,
+                                        data_points,
+                                        z,
+                                        alpha = 0,
+                                        verbose=False)
+
+        #print vertex_attributes
+        quantity.set_values(vertex_attributes)
+        assert allclose(quantity.vertex_values.flat, answer)        
 
 

inundation/ga/storm_surge/validation/LWRU2/lwru2.py

@@ -7 +7 @@
 http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
 
-
 """
 
-######################
+
 # Module imports
-
-from pyvolution.shallow_water import Domain, Reflective_boundary,\
-     Dirichlet_boundary,Transmissive_boundary, Constant_height, Constant_stage
-
+from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary
 from pyvolution.mesh_factory import rectangular
 from Numeric import array, zeros, Float, allclose
 
+picklefile = 'domain.pck'
+from cPickle import dump, load
+
+try:
+    fid = open(picklefile)
+    domain = load(fid)
+    print 'Read pickled domain'
+except:
+
+    #Preparing points
+    from pyvolution.data_manager import xya2pts
+    xya2pts('bathymetry.txt', verbose = True)
 
 
-#Preparing points
-from pyvolution.data_manager import xya2pts
-xya2pts('Benchmark_2_Bathymetry.txt', verbose = True)
+    #######################
+    # Domain
+    print 'Creating domain'
+    points, vertices, boundary = rectangular(200, 200/5*3,
+                                             len1=5.448, len2=3.402)
+    domain = Domain(points, vertices, boundary)
+
+    domain.check_integrity()
+    domain.default_order = 2
+
+    print "Number of triangles = ", len(domain)
+    domain.store = True    #Store for visualisation purposes
+
+    import sys, os
+    base = os.path.basename(sys.argv[0])
+    domain.filename, _ = os.path.splitext(base)
+
+    #LS code to be included in set_quantity
+    from pyvolution import util, least_squares
+    points, attributes = util.read_xya('bathymetry.pts')
+    
+    #Fit attributes to mesh
+    vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
+                                                  domain.triangles,
+                                                  points,
+                                                  attributes['elevation'],
+                                                  verbose=True)
+
+    print 'Initial values'
+    domain.set_quantity('elevation', vertex_attributes)
+    domain.set_quantity('friction', 0.0)
+    domain.set_quantity('stage', 0.0)
+
+
+    fid = open('domain.pck', 'w')
+    dump(domain, fid)
+    fid.close()
 
 
 
-#######################
-# Domain
-#
-
-
-print 'Creating domain'
-#Create basic mesh
-#
-#The initial condition extends 50km off shore
-#and 5,000m is allowed on shore for wetting
-#(only about 200m is expected, though)
-
-#points, vertices, boundary = rectangular(400, 40,
-#                                         len1=55000, len2=5000,
-#                                         origin = (-5000, 0.0))
-
-points, vertices, boundary = rectangular(100, 100,
-                                         len1=5.448, len2=3.402)
-
-
-#Create shallow water domain
-domain = Domain(points, vertices, boundary)
-
-domain.check_integrity()
-domain.default_order = 2
-
-#Output params
-domain.smooth = True
-domain.reduction = min  #Looks a lot better on top of steep slopes
-print "Number of triangles = ", len(domain)
-
-domain.visualise = False
-domain.store = True    #Store for visualisation purposes
-domain.format = 'sww'   #Native netcdf visualisation format
-
-import sys, os
-base = os.path.basename(sys.argv[0])
-domain.filename, _ = os.path.splitext(base)
-
-
-#LS code to be included in set_quantity
-from pyvolution import util, least_squares
-
-points, attributes = util.read_xya('Benchmark_2_Bathymetry.pts')
-    
-#Fit attributes to mesh
-vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
-                                              domain.triangles,
-                                              points,
-                                              attributes['elevation'],
-                                              verbose=True)
-
-
-print 'Initial values'
-domain.set_quantity('elevation', vertex_attributes)
-#domain.set_quantity('elevation', 'Benchmark_2_Bathymetry.pts')
-domain.set_quantity('friction', 0.0)
-domain.set_quantity('stage', 0.0)
-
-#import sys; sys.exit()
-
-#print domain.quantities['stage'].centroid_values
 
 ######################
@@ -95 +75 @@
 print 'Boundaries'
 Br = Reflective_boundary(domain)
-Bt = Transmissive_boundary(domain)
-
-#Constant inflow
-Bd = Dirichlet_boundary([0.0, 0.0, 0.0])
+Bf = File_boundary('input_wave.txt', domain)
 
 #Set boundary conditions
-domain.set_boundary({'left': Br, 'right': Bt, 'bottom': Bt, 'top': Bt})
+domain.set_boundary({'left': Bf, 'right': Br, 'bottom': Br, 'top': Br})
 
 
@@ -108 +85 @@
 t0 = time.time()
 
-
-
-
-pt = []
-xes = []
-y = 2500 
-x0 = -100
-step = 50
-for i in range(1000):
-    x = x0+i*step
-    xes.append(x)
-    pt.append( [x,y] )
-
-from pylab import *
-from pyvolution.least_squares import Interpolation
-
- 
-V = domain.get_vertex_coordinates(obj=True)
-T = domain.get_triangles(obj=True)
-
-I = Interpolation(V,
-                  T,
-                  point_coordinates = pt,
-                  verbose = True)
-
-
-f = domain.quantities['elevation'].vertex_values.flat
-z = I.interpolate( f )
-
-print 'xxxxx'
-
-ion()
-hold(False)
-f = domain.quantities['stage'].vertex_values.flat
-y = I.interpolate( f )
-plot(xes, y, '-b', xes, z, '-k')
-set( gca(), Ylim=(-10,10) )
-
-raw_input('go')
-for t in domain.evolve(yieldstep = 1, finaltime = 300.0):
-
-    f = domain.quantities['stage'].vertex_values.flat
-    y = I.interpolate( f )
-    ioff()
-    plot(xes, y, '-b', xes, z, '-k')
-    ion()
-    set( gca(), Ylim=(-10,10) )
-    
-    
-    #print y[:], y.shape
-
-
+for t in domain.evolve(yieldstep = 0.1, finaltime = 22.5):
     domain.write_time()
 
 
 print 'That took %.2f seconds' %(time.time()-t0)
-show()
+