Changeset 5864

Timestamp:

Oct 24, 2008, 5:03:31 PM (15 years ago)

Author:

ole

Message:

Re-use of the interpolation matrix when interpolation points in get_values
are the same as in previous calls. This gives a very significant speed up
in get_energy_through_cross_sections.

It is a little rough around the edges, but the idea is very sound and
could be considered as a part of the interpolation code.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (3 diffs)
• culvert_flows/Test_Culvert_Flat_Water_Lev.py (modified) (1 diff)
• fit_interpolate/interpolate.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -16 +16 @@
 
 from Numeric import array, zeros, Float, less, concatenate, NewAxis,\
-     argmax, argmin, allclose, take, reshape
+     argmax, argmin, allclose, take, reshape, alltrue
 
 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
@@ -1029 +1029 @@
 
 
-
-    def get_interpolated_values(self, interpolation_points,
+    def obsolete_get_interpolated_values(self, interpolation_points,
                                 verbose=False):
         """ Get values at interpolation points
@@ -1075 +1074 @@
         # Use caching to reuse interpolation information
         from anuga.fit_interpolate.interpolate import interpolate
+        
+        return interpolate(vertex_coordinates, 
+                           triangles, 
+                           vertex_values, 
+                           interpolation_points,
+                           use_cache=True, 
+                           verbose=verbose)
+        
+        
+
+    def get_interpolated_values(self, interpolation_points,
+                                verbose=False):
+        """ Get values at interpolation points
+        
+        If interpolation points have been given previously, the 
+        associated matrices will be reused to save time.
+        
+        The argument interpolation points must be given as either a 
+        list of absolute UTM coordinates or a geospatial data object.
+        """
+
+        # FIXME (Ole): Could do with an input check (should be generalised 
+        # and used widely)
+        # That will check that interpolation points is either a list of 
+        # points, Nx2 array, or geospatial
+        
+        # Ensure points are converted to coordinates relative to mesh origin
+        # FIXME (Ole): This could all be refactored using the 
+        # 'change_points_geo_ref' method of Class geo_reference. 
+        # The purpose is to make interpolation points relative
+        # to the mesh origin.
+        #
+        # Speed is also a consideration here.
+        
+        # FIXME (Ole): The reuse of I is a little rough, but it sure does speed up
+        # frequent interpolations such as those found in culvert_class.
+        
+        use_cache = True
+        
+        # Get internal (discontinuous) triangles for use with the 
+        # interpolation object.
+        x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
+                                                                smooth=False)
+        # FIXME (Ole): This concat should roll into get_vertex_values
+        vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
+                                         axis=1)                                                                
+            
+            
+        if isinstance(interpolation_points, Geospatial_data):        
+            # Ensure interpolation points are in absolute UTM coordinates
+            interpolation_points = interpolation_points.get_data_points(absolute=True)
+                
+        # Reconcile interpolation points with georeference of domain
+        interpolation_points = self.domain.geo_reference.get_relative(interpolation_points) 
+        interpolation_points = ensure_numeric(interpolation_points)
+        
+        # Check if same interpolation points have been used before
+        if hasattr(self, 'stored_I') and\
+                alltrue(interpolation_points == self.stored_interpolation_points):
+            I = self.stored_I
+        else:    
+        
+            # Use caching to reuse interpolation information
+            from anuga.fit_interpolate.interpolate import Interpolate
+            from anuga.caching import cache
+    
+
+            # Create interpolation object with matrix
+            args = (ensure_numeric(vertex_coordinates, Float), 
+                    ensure_numeric(triangles))
+                      
+            if use_cache is True:
+                I = cache(Interpolate, args,
+                          verbose=verbose)
+            else:
+                I = apply(Interpolate, args)
+                          
+            self.stored_I = I    
+            self.stored_interpolation_points = interpolation_points
+            
+            
+    
+        # Call interpolate method with interpolation points
+        result = I.interpolate_block(vertex_values, interpolation_points,
+                                     use_cache=use_cache,
+                                     verbose=verbose)
+                               
+        return result
+        
         
         return interpolate(vertex_coordinates, 

anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

@@ -158 +158 @@
 #------------------------------------------------------------------------------
 
-#for t in domain.evolve(yieldstep = 0.01, finaltime = 45):
-#     if int(domain.time*100) % 100 == 0:
-#            domain.write_time()
+for t in domain.evolve(yieldstep = 0.01, finaltime = 45):
+     if int(domain.time*100) % 100 == 0:
+             domain.write_time()
     
 
+    
+    
+import sys; sys.exit() 
 # Profiling code
 import time

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -373 +373 @@
         See interpolate for doc info.
         """
-        if isinstance(point_coordinates,Geospatial_data):
+        if isinstance(point_coordinates, Geospatial_data):
             point_coordinates = point_coordinates.get_data_points(\
                 absolute=True)