Changeset 896

Timestamp:

Feb 15, 2005, 5:25:14 PM (20 years ago)

Author:

ole

Message:

Fiddle

Location:

inundation/ga/storm_surge/pyvolution

Files:

• data_manager.py (modified) (3 diffs)
• test_shallow_water.py (modified) (4 diffs)
• test_util.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -1010 +1010 @@
                minlat = None, maxlat =None,
                minlon = None, maxlon =None,
-               mint = None, maxt = None, mean_stage = 0):
+               mint = None, maxt = None, mean_stage = 0,
+               origin = (0,0,0)):
     """Convert 'Ferret' NetCDF format for wave propagation to
     sww format native to pyvolution.
@@ -1021 +1022 @@
     assumed to be given in the GDA94 datum.
 
-
     min's and max's: If omitted - full extend is used.
     To include a value min may equal it, while max must exceed it.
+
+    origin is a 3-tuple with geo referenced
+    UTM coordinates (zone, easting, northing)  
     """ 
 
@@ -1195 +1198 @@
     volumes = array(volumes)        
 
-    origin = (min(x), min(y))
-    outfile.minx = origin[0]
-    outfile.miny = origin[1]    
-
-    #print x - origin[0]
-    #print y - origin[1]
-
-    #print len(x)
-    #print number_of_points
-    
-    outfile.variables['x'][:] = x - origin[0]
-    outfile.variables['y'][:] = y - origin[1]
+    if origin == None:
+        zone = refzone        
+        xllcorner = min(x)
+        yllcorner = min(y)
+    else:
+        zone = origin[0]
+        xllcorner = origin[1]
+        yllcorner = origin[2]                
+
+    
+    outfile.xllcorner = xllcorner
+    outfile.yllcorner = yllcorner
+    outfile.zone = zone
+
+    outfile.variables['x'][:] = x - xllcorner
+    outfile.variables['y'][:] = y - yllcorner
     outfile.variables['z'][:] = 0.0
     outfile.variables['elevation'][:] = 0.0    

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -2399 +2399 @@
         from util import mean
         domain1.reduction = mean 
-        domain1.smooth = True #FIXME: This may be a requirement for the 
-                              #interpolation process 
+        domain1.smooth = False  #Exact result
+
+        domain1.default_order = 2
+        domain1.store = True    
+        domain1.set_datadir('.')
+        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
+        
+        #FIXME: This is extremely important!
+        #How can we test if they weren't stored?
+        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+                
+        #Bed-slope and friction at vertices (and interpolated elsewhere)
+        domain1.set_quantity('elevation', 0)
+        domain1.set_quantity('friction', 0)     
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain1)
+        Bd = Dirichlet_boundary([0.3,0,0])
+        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
+        #Initial condition
+        domain1.set_quantity('stage', 0)
+        domain1.check_integrity()
+
+        finaltime = 5
+        #Evolution  (full domain - large steps)
+        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
+            pass
+            #domain1.write_time()
+
+        cv1 = domain1.quantities['stage'].centroid_values
+        
+        
+        #Create an triangle shaped domain (reusing coordinates from domain 1), 
+        #formed from the lower and right hand  boundaries and 
+        #the sw-ne diagonal 
+        #from domain 1. Call it domain2
+        
+        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
+                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
+                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
+        
+        vertices = [ [1,2,0], 
+                     [3,4,1], [2,1,4], [4,5,2], 
+                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]       
+
+        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
+                     (4,2):'right', (6,2):'right', (8,2):'right',
+                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}    
+                     
+        domain2 = Domain(points, vertices, boundary)
+        
+        domain2.reduction = domain1.reduction
+        domain2.smooth = False
+        domain2.default_order = 2
+        
+        #Bed-slope and friction at vertices (and interpolated elsewhere)
+        domain2.set_quantity('elevation', 0)
+        domain2.set_quantity('friction', 0)
+        domain2.set_quantity('stage', 0)        
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain2)        
+        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
+                                      domain2)
+        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
+        domain2.check_integrity()        
+        
+
+
+        #Evolution (small steps)
+        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
+            pass
+
+        
+        #Use output from domain1 as spatio-temporal boundary for domain2
+        #and verify that results at right hand side are close.  
+
+        cv2 = domain2.quantities['stage'].centroid_values
+
+        #print take(cv1, (12,14,16))  #Right
+        #print take(cv2, (4,6,8))
+        #print take(cv1, (0,6,12))  #Bottom     
+        #print take(cv2, (0,1,4))
+        #print take(cv1, (0,8,16))  #Diag       
+        #print take(cv2, (0,3,8))
+
+        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
+        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
+        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
+
+        #Cleanup
+        os.remove(domain1.filename + '.' + domain1.format)        
+        
+        
+
+    def test_spatio_temporal_boundary_2(self):
+        """Test that boundary values can be read from file and interpolated
+        in both time and space.
+        This is a more basic test, verifying that boundary object 
+        produces the expected results
+
+
+        """
+        import time
+       
+        #Create sww file of simple propagation from left to right 
+        #through rectangular domain
+        
+        from mesh_factory import rectangular
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(3, 3)
+
+        #Create shallow water domain
+        domain1 = Domain(points, vertices, boundary)
+        
+        from util import mean
+        domain1.reduction = mean 
+        domain1.smooth = True #To mimic MOST output
 
         domain1.default_order = 2
@@ -2426 +2544 @@
         finaltime = 5
         #Evolution  (full domain - large steps)
-        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
+        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
             pass
             #domain1.write_time()
 
-        cv1 = domain1.quantities['stage'].centroid_values
-        
-        
-        #Create an triangle shaped domain (reusing coordinates from domain 1), 
+        
+        #Create an triangle shaped domain (coinciding with some
+        #coordinates from domain 1), 
         #formed from the lower and right hand  boundaries and 
         #the sw-ne diagonal 
@@ -2453 +2570 @@
         
         domain2.reduction = domain1.reduction
-        domain2.smooth = True
-        domain2.default_order = 2
-        
-        #Bed-slope and friction at vertices (and interpolated elsewhere)
-        domain2.set_quantity('elevation', 0)
-        domain2.set_quantity('friction', 0)
-        domain2.set_quantity('stage', 0)        
-
-        # Boundary conditions
-        Br = Reflective_boundary(domain2)        
-        Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
-                                      domain2)
-        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
-        domain2.check_integrity()        
-        
-
-
-        #Evolution (small steps)
-        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
-            pass
-
-        
-        #Use output from domain1 as spatio-temporal boundary for domain2
-        #and verify that results at right hand side are close.  
-
-        cv2 = domain2.quantities['stage'].centroid_values
-
-        #print take(cv1, (12,14,16))  #Right
-        #print take(cv2, (4,6,8))
-        #print take(cv1, (0,6,12))  #Bottom     
-        #print take(cv2, (0,1,4))
-        #print take(cv1, (0,8,16))  #Diag       
-        #print take(cv2, (0,3,8))
-
-        assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
-        assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
-        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
-
-        #Cleanup
-        os.remove(domain1.filename + '.' + domain1.format)        
-        
-        
-
-    def test_spatio_temporal_boundary_2(self):
-        """Test that boundary values can be read from file and interpolated
-        in both time and space.
-        This is a more basic test, verifying that boundary object 
-        produces the expected results
-
-
-        """
-        import time
-       
-        #Create sww file of simple propagation from left to right 
-        #through rectangular domain
-        
-        from mesh_factory import rectangular
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(3, 3)
-
-        #Create shallow water domain
-        domain1 = Domain(points, vertices, boundary)
-        
-        from util import mean
-        domain1.reduction = mean 
-        domain1.smooth = True #FIXME: This may be a requirement for the 
-                              #interpolation process 
-
-        domain1.default_order = 2
-        domain1.store = True    
-        domain1.set_datadir('.')
-        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
-        
-        #FIXME: This is extremely important!
-        #How can we test if they weren't stored?
-        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
-
-                
-        #Bed-slope and friction at vertices (and interpolated elsewhere)
-        domain1.set_quantity('elevation', 0)
-        domain1.set_quantity('friction', 0)     
-
-        # Boundary conditions
-        Br = Reflective_boundary(domain1)
-        Bd = Dirichlet_boundary([0.3,0,0])      
-        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
-        #Initial condition
-        domain1.set_quantity('stage', 0)
-        domain1.check_integrity()
-
-        finaltime = 5
-        #Evolution  (full domain - large steps)
-        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
-            pass
-            #domain1.write_time()
-
-        
-        #Create an triangle shaped domain (reusing coordinates from domain 1), 
-        #formed from the lower and right hand  boundaries and 
-        #the sw-ne diagonal 
-        #from domain 1. Call it domain2
-        
-        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
-                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
-                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
-        
-        vertices = [ [1,2,0], 
-                     [3,4,1], [2,1,4], [4,5,2], 
-                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]       
-
-        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
-                     (4,2):'right', (6,2):'right', (8,2):'right',
-                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}    
-                     
-        domain2 = Domain(points, vertices, boundary)
-        
-        domain2.reduction = domain1.reduction
-        domain2.smooth = True
+        domain2.smooth = False
         domain2.default_order = 2
         
@@ -2594 +2593 @@
         points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)       
         #The diagonal points of domain 1 are 0, 5, 10, 15
-        #print points[0], points[5], points[15]
+        
+        #print points[0], points[5], points[10], points[15]
         assert allclose( take(points, [0,5,10,15]), 
                          [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])  

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -358 +358 @@
         from util import mean
         domain1.reduction = mean 
-        domain1.smooth = True #NOTE: This is a strict requirement for the 
-                              #interpolation process 
+        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
+                              # only one value.
 
         domain1.default_order = 2