Changeset 5657

Timestamp:

Aug 14, 2008, 10:26:06 AM (16 years ago)

Author:

ole

Message:

Implemented default_boundary option in File_boundary and Field_boundary as
per ticket:293 and added a note in the user manual.

Location:

anuga_core

Files:

• documentation/user_manual/anuga_user_manual.tex (modified) (1 diff)
• source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (5 diffs)
• source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py (modified) (2 diffs)
• source/anuga/abstract_2d_finite_volumes/test_util.py (modified) (4 diffs)
• source/anuga/fit_interpolate/interpolate.py (modified) (2 diffs)
• source/anuga/shallow_water/shallow_water_domain.py (modified) (3 diffs)
• source/anuga/shallow_water/test_data_manager.py (modified) (3 diffs)

anuga_core/documentation/user_manual/anuga_user_manual.tex

@@ -2253 +2253 @@
 The boundary values are obtained from a file and interpolated to the
 appropriate segments for each conserved quantity.
+
+Optional argument \code{default_boundary} can be used to specify another boundary object to be used in case model time exceeds the time availabel in the file used by \code{File\_boundary}.
+The default_boundary could be a simple Dirichlet condition or even another File\_boundary 
+typically using data pertaining to another time interval.  
 \end{classdesc}
 

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -2 +2 @@
 """
 
+from warnings import warn
+
 from anuga.utilities.numerical_tools import NAN    
+from anuga.fit_interpolate.interpolate import Modeltime_too_late
+from anuga.fit_interpolate.interpolate import Modeltime_too_early
 
 
@@ -183 +187 @@
     Note that the resulting solution history is not exactly the same as if
     the models were coupled as there is no feedback into the source model.
+    
+    Optional keyword argument default_boundary must be either None or 
+    an instance of class descending from class Boundary.
+    This will be used in case model time exceeds that available in the 
+    underlying data.
        
     """
@@ -189 +198 @@
     # rather than File_boundary
 
-    def __init__(self, filename, domain, time_thinning=1, 
-                 use_cache=False, verbose=False, boundary_polygon=None):
+    def __init__(self, filename, domain, 
+                 time_thinning=1, 
+                 boundary_polygon=None,    
+                 default_boundary=None,
+                 use_cache=False, 
+                 verbose=False): 
+
         import time
         from Numeric import array, zeros, Float
@@ -251 +265 @@
                                verbose=verbose,
                                boundary_polygon=boundary_polygon)
-
+                             
+        # Check and store default_boundary
+        msg = 'Keyword argument default_boundary must be either None '
+        msg += 'or a boundary object.\n I got %s' %(str(default_boundary))
+        assert default_boundary is None or isinstance(default_boundary, Boundary), msg
+        self.default_boundary = default_boundary
+        self.default_boundary_invoked = False    # Flag
+
+        # Store pointer to domain
         self.domain = domain
+        
+        self.verbose = verbose
 
         # Test
@@ -293 +317 @@
 
         t = self.domain.time
+        
         if vol_id is not None and edge_id is not None:
             i = self.boundary_indices[ vol_id, edge_id ]
-            res = self.F(t, point_id = i)
+            
+            
+            #res = self.F(t, point_id = i)            
+            try:
+                res = self.F(t, point_id = i)
+            except (Modeltime_too_late, Modeltime_too_early), e:
+                if self.default_boundary is None:
+                    raise Exception, e # Reraise exception
+                else:
+                    if self.default_boundary_invoked is False:
+                        # Issue warning the first time
+                        msg = '%s' %str(e)
+                        msg += 'Instead I will using the default boundary: %s\n'\
+                            %str(self.default_boundary) 
+                        msg += 'Note: Further warnings will be supressed'
+                        warn(msg)
+                   
+                        # FIXME (Ole): Replace this crude flag with
+                        # Python's ability to print warnings only once.
+                        # See http://docs.python.org/lib/warning-filter.html
+                        self.default_boundary_invoked = True
+                    
+                    # Pass control to default boundary
+                    res = self.default_boundary.evaluate(vol_id, edge_id)
+                    
+                
             if res == NAN:
                 x,y=self.midpoint_coordinates[i,:]

anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

@@ -104 +104 @@
         """Test that boundary values can be read from file and interpolated
         This is using the .tms file format
+        
+        See also test_util for comprenhensive testing of the underlying 
+        file_function and also tests in test_datamanager which tests 
+        file_function using the sts format
         """
 
@@ -269 +273 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Generic_Boundary_Conditions,'test')
+    suite = unittest.makeSuite(Test_Generic_Boundary_Conditions, 'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

@@ -521 +521 @@
 
 
-    def qtest_spatio_temporal_file_function_time(self):
+    def test_spatio_temporal_file_function_time(self):
         """Test that File function interpolates correctly
         between given times.
@@ -554 +554 @@
         points, vertices, boundary =\
                 rectangular(4, 4, 15, 30, origin = (0, -20))
-        print "points", points
+        #print "points", points
 
         #print 'Number of elements', len(vertices)
@@ -639 +639 @@
                 #print ' ', q0
                 #print ' ', q1
-                print "q",q
-                print "actual", actual
+                #print "q",q
+                #print "actual", actual
                 #print
                 if q0 == NAN:
@@ -696 +696 @@
 
 
-    def NOtest_spatio_temporal_file_function_time(self):
+    def Xtest_spatio_temporal_file_function_time(self):
         # FIXME: This passes but needs some TLC
         # Test that File function interpolates correctly

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -39 +39 @@
 from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
 from anuga.abstract_2d_finite_volumes.util import file_function
+
+# Interpolation specific exceptions
+
+class Modeltime_too_late(Exception): pass
+class Modeltime_too_early(Exception): pass
+
 
 
@@ -856 +862 @@
         msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
         msg += ' does not match model time: %.16f\n' %t
-        if t < self.time[0]: raise Exception(msg)
-        if t > self.time[-1]: raise Exception(msg)
+        if t < self.time[0]: raise Modeltime_too_early(msg)
+        if t > self.time[-1]: raise Modeltime_too_late(msg)
 
         oldindex = self.index #Time index

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -1088 +1088 @@
     def __init__(self, filename, domain,
                  mean_stage=0.0,
-                 time_thinning=1, 
+                 time_thinning=1,
+                 boundary_polygon=None,    
+                 default_boundary=None,                  
                  use_cache=False,
-                 verbose=False,
-                 boundary_polygon=None):
+                 verbose=False):
         """Constructor
 
@@ -1107 +1108 @@
                        the speed of a model run that you are setting up 
                        and testing.
+                       
+        default_boundary: Must be either None or an instance of a 
+                          class descending from class Boundary.
+                          This will be used in case model time exceeds 
+                          that available in the underlying data.
+                                                
         use_cache:
         verbose:
@@ -1115 +1122 @@
         self.file_boundary = File_boundary(filename, domain,
                                            time_thinning=time_thinning,
+                                           boundary_polygon=boundary_polygon,
+                                           default_boundary=default_boundary,
                                            use_cache=use_cache,
-                                           verbose=verbose,
-                                           boundary_polygon=boundary_polygon) 
+                                           verbose=verbose)
+
         
         # Record information from File_boundary

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -7643 +7643 @@
         """
         
+        # FIXME (Ole): These tests should really move to 
+        # test_generic_boundaries.py
+        
         from anuga.shallow_water import Domain
         from anuga.shallow_water import Reflective_boundary
@@ -7774 +7777 @@
         os.remove(sts_file+'.sts')
         os.remove(meshname)
+        
+        
+        
+                
+        
+    def test_file_boundary_stsI_beyond_model_time(self):
+        """test_file_boundary_stsI(self):
+        
+        Test that file_boundary can work when model time
+        exceeds available data.
+        This is optional and requires the user to specify a default 
+        boundary object.
+        """
+        
+        # Don't do warnings in unit test
+        import warnings
+        warnings.simplefilter('ignore')
+        
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water import Reflective_boundary
+        from anuga.shallow_water import Dirichlet_boundary
+        from anuga.shallow_water import File_boundary
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
+        tide = 0.37
+        time_step_count = 5
+        time_step = 2
+        lat_long_points =bounding_polygon[0:3]
+        n=len(lat_long_points)
+        first_tstep=ones(n,Int)
+        last_tstep=(time_step_count)*ones(n,Int)
+
+        h = 20        
+        w = 2
+        u = 10
+        v = -10
+        gauge_depth=h*ones(n,Float)
+        ha=w*ones((n,time_step_count),Float)
+        ua=u*ones((n,time_step_count),Float)
+        va=v*ones((n,time_step_count),Float)
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha,
+                                           ua=ua,
+                                           va=va)
+
+        sts_file=base_name
+        urs2sts(base_name,
+                sts_file,
+                mean_stage=tide,
+                verbose=False)
+        self.delete_mux(files)
+
+        #print 'start create mesh from regions'
+        for i in range(len(bounding_polygon)):
+            zone,\
+            bounding_polygon[i][0],\
+            bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],
+                                            bounding_polygon[i][1])
+                                            
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags=boundary_tags,
+                                 maximum_triangle_area=extent_res,
+                                 filename=meshname,
+                                 interior_regions=interior_regions,
+                                 verbose=False)
+        
+        domain_fbound = Domain(meshname)
+        domain_fbound.set_quantity('stage', tide)
+        
+        Br = Reflective_boundary(domain_fbound)
+        Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)])        
+        Bf = File_boundary(sts_file+'.sts', 
+                           domain_fbound, 
+                           boundary_polygon=bounding_polygon,
+                           default_boundary=Bd) # Condition to be used 
+                                                # if model time exceeds 
+                                                # available data
+
+        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
+        
+        # Check boundary object evaluates as it should
+        for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects):
+            if B is Bf:
+            
+                qf = B.evaluate(vol_id, edge_id)  # File boundary
+                qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary
+
+                assert allclose(qf, qd) 
+                
+        
+        # Evolve
+        data_finaltime = time_step*(time_step_count-1)
+        finaltime = data_finaltime + 10 # Let model time exceed available data
+        yieldstep = time_step
+        temp_fbound=zeros(int(finaltime/yieldstep)+1, Float)
+
+        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step=False)):
+                                                   
+            D = domain_fbound
+            temp_fbound[i]=D.quantities['stage'].centroid_values[2]
+
+            # Check that file boundary object has populated 
+            # boundary array correctly  
+            # FIXME (Ole): Do this for the other tests too!
+            for j, val in enumerate(D.get_quantity('stage').boundary_values):
+            
+                (vol_id, edge_id), B = D.boundary_objects[j]
+                if isinstance(B, File_boundary):
+                    #print j, val
+                    assert allclose(val, w + tide)
+
+
+
+        
+        
+        
+            
+
+        
+        domain_drchlt = Domain(meshname)
+        domain_drchlt.set_quantity('stage', tide)
+        Br = Reflective_boundary(domain_drchlt)
+
+        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
+        temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float)
+
+        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step=False)):
+            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
+
+        #print domain_fbound.quantities['stage'].vertex_values
+        #print domain_drchlt.quantities['stage'].vertex_values
+                    
+        assert allclose(temp_fbound,temp_drchlt)
+        
+        assert allclose(domain_fbound.quantities['stage'].vertex_values,
+                        domain_drchlt.quantities['stage'].vertex_values)
+                        
+        assert allclose(domain_fbound.quantities['xmomentum'].vertex_values,
+                        domain_drchlt.quantities['xmomentum'].vertex_values)                        
+                        
+        assert allclose(domain_fbound.quantities['ymomentum'].vertex_values,
+                        domain_drchlt.quantities['ymomentum'].vertex_values)                                                
+        
+        
+        os.remove(sts_file+'.sts')
+        os.remove(meshname)
+
+        
 
     def test_file_boundary_stsII(self):
@@ -10105 +10268 @@
 
     suite = unittest.makeSuite(Test_Data_Manager,'test')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsI')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsI_beyond_model_time')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV_sinewave_ordering')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')