Changeset 4588 for anuga_core

Timestamp:

Jul 4, 2007, 6:09:12 PM (18 years ago)

Author:

ole

Message:

Added exception and helpful diagnostics if a boundary condition tries
to access file_boundary outside its area defined by the underlying
sww file.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (2 diffs)
• abstract_2d_finite_volumes/test_util.py (modified) (1 diff)
• fit_interpolate/interpolate.py (modified) (2 diffs)
• fit_interpolate/test_interpolate.py (modified) (2 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (3 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -1 +1 @@
 """boundary.py - Classes for implementing boundary conditions
 """
+
+from anuga.utilities.numerical_tools import NAN    
 
 
@@ -271 +273 @@
         if vol_id is not None and edge_id is not None:
             i = self.boundary_indices[ vol_id, edge_id ]
-            return self.F(t, point_id = i)
+            res = self.F(t, point_id = i)
+
+            if res == NAN:
+                x,y=self.midpoint_coordinates[i,:]
+                msg = 'NAN value found in file_boundary at '
+                msg += 'point id #%d: (%.2f, %.2f)' %(i, x, y)
+                #print msg
+                raise Exception, msg
+            
+            return res 
         else:
             #raise 'Boundary call without point_id not implemented'

anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

@@ -689 +689 @@
 
 
-    def xtest_spatio_temporal_file_function_time(self):
+    def NOtest_spatio_temporal_file_function_time(self):
+        # FIXME: This passes but needs some TLC
         # Test that File function interpolates correctly
         # When some points are outside the mesh

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -548 +548 @@
                                                   point_coordinates=\
                                                   self.interpolation_points)
+
+                    #assert len(result), len(interpolation_points)
                     self.precomputed_values[name][i, :] = result
 
@@ -645 +647 @@
                     if ratio > 0:
                         Q1 = Q[self.index+1, point_id]
-            
+
             #Linear temporal interpolation    
             if ratio > 0:

anuga_core/source/anuga/fit_interpolate/test_interpolate.py

@@ -14 +14 @@
 
 from Scientific.IO.NetCDF import NetCDFFile
-from Numeric import allclose, array, transpose, zeros, Float
+from Numeric import allclose, array, transpose, zeros, Float, sometrue, alltrue, take, where
 
 
@@ -1336 +1336 @@
                                    verbose = False)
 
+        
+        assert alltrue(I.precomputed_values['Attribute'][:,4] != NAN)
+        assert sometrue(I.precomputed_values['Attribute'][:,5] == NAN)
+
+        #X = I.precomputed_values['Attribute'][1,:]
+        #print X
+        #print take(X, X == NAN)
+        #print where(X == NAN, range(len(X)), 0)        
+        
         answer = linear_function(interpolation_points)
           

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -4274 +4274 @@
         #from domain 1. Call it domain2
 
-        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
+        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]]
-
+                   [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],
@@ -4386 +4387 @@
         R0 = Bf.evaluate(8,0)
         assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3 + mean_stage)
+
+
+        #Cleanup
+        os.remove(domain1.get_name() + '.' + domain1.format)
+
+
+    def test_spatio_temporal_boundary_outside(self):
+        """Test that field_boundary catches if a point is outside the sww that defines it
+        """
+
+        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)
+
+        domain1.reduction = mean
+        domain1.smooth = True #To mimic MOST output
+
+        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 (coinciding with some
+        #coordinates from domain 1, but one edge outside!),
+        #formed from the lower and right hand  boundaries and
+        #the sw-ne diagonal as in the previous test but scaled
+        #in the x direction by a factor of 2
+
+        points = [ [0,0],
+                   [2.0/3,0], [2.0/3,1.0/3],
+                   [4.0/3,0], [4.0/3,1.0/3], [4.0/3,2.0/3],
+                   [2,0],     [2,1.0/3],     [2,2.0/3],     [2,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)
+
+
+        #Read results for specific timesteps t=1 and t=2
+        from Scientific.IO.NetCDF import NetCDFFile
+        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        s1 = fid.variables['stage'][1,:]
+        s2 = fid.variables['stage'][2,:]
+        fid.close()
+
+        from Numeric import take, reshape, concatenate
+        shp = (len(x), 1)
+        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[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]])
+
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain2)
+        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
+        #                              domain2)
+        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
+                            domain2, mean_stage=1)
+        
+        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
+        domain2.check_integrity()
+
+        try:
+            for t in domain2.evolve(yieldstep = 1, finaltime = finaltime):
+                pass
+        except:
+            pass
+        else:
+            msg = 'This should have caught NAN at boundary'
+            raise Exception, msg
 
 
@@ -4619 +4744 @@
         
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Shallow_Water,'test')
+    suite = unittest.makeSuite(Test_Shallow_Water,'test')    
+    #suite = unittest.makeSuite(Test_Shallow_Water,'test_spatio_temporal_boundary_outside')
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_temp')