Changeset 6173

Timestamp:

Jan 14, 2009, 5:47:17 PM (16 years ago)

Author:

ole

Message:

Introduced time_limit in Field_boundary, File_boundary and file_function

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) (2 diffs)
• abstract_2d_finite_volumes/util.py (modified) (9 diffs)
• fit_interpolate/interpolate.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (2 diffs)
• shallow_water/test_data_manager.py (modified) (5 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -158 +158 @@
     def __init__(self, filename, domain, 
                  time_thinning=1, 
+                 time_limit=None,
                  boundary_polygon=None,    
                  default_boundary=None,
@@ -219 +220 @@
                                interpolation_points=self.midpoint_coordinates,
                                time_thinning=time_thinning,
+                               time_limit=time_limit,
                                use_cache=use_cache, 
                                verbose=verbose,

anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

@@ -1052 +1052 @@
         assert num.allclose(domain.starttime, start+delta)
 
-
+        assert num.allclose(F.get_time(), [-23., 37., 97., 157., 217.,
+                                            277., 337., 397., 457., 517.,
+                                            577., 637., 697., 757., 817.,
+                                            877., 937., 997., 1057., 1117.,
+                                            1177.])
 
 
@@ -1083 +1087 @@
         os.remove(filename + '.txt')                
 
+        
+
+    def test_file_function_time_with_domain_different_start_and_time_limit(self):
+        """Test that File function interpolates correctly
+        between given times. No x,y dependency here.
+        Use domain with a starttime later than that of file
+
+        ASCII version
+        
+        This test also tests that time can be truncated.
+        """
+
+        # Write file
+        import os, time, calendar
+        from anuga.config import time_format
+        from math import sin, pi
+        from domain import Domain
+
+        finaltime = 1200
+        filename = 'test_file_function'
+        fid = open(filename + '.txt', 'w')
+        start = time.mktime(time.strptime('2000', '%Y'))
+        dt = 60  #One minute intervals
+        t = 0.0
+        while t <= finaltime:
+            t_string = time.strftime(time_format, time.gmtime(t+start))
+            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
+            t += dt
+
+        fid.close()
+
+        # Convert ASCII file to NetCDF (Which is what we really like!)
+        timefile2netcdf(filename)        
+
+        a = [0.0, 0.0]
+        b = [4.0, 0.0]
+        c = [0.0, 3.0]
+
+        points = [a, b, c]
+        vertices = [[0,1,2]]
+        domain = Domain(points, vertices)
+
+        # Check that domain.starttime isn't updated if later than file starttime but earlier
+        # than file end time
+        delta = 23
+        domain.starttime = start + delta
+        F = file_function(filename + '.tms', domain,
+                          time_limit=600,
+                          quantities=['Attribute0', 'Attribute1', 'Attribute2'])        
+        assert num.allclose(domain.starttime, start+delta)
+
+        assert num.allclose(F.get_time(), [-23., 37., 97., 157., 217.,
+                                            277., 337., 397., 457., 517.,
+                                            577.])        
+
+
+
+        # Now try interpolation with delta offset
+        for i in range(20):
+            t = i*10
+            q = F(t-delta)
+
+            #Exact linear intpolation
+            assert num.allclose(q[0], 2*t)
+            if i%6 == 0:
+                assert num.allclose(q[1], t**2)
+                assert num.allclose(q[2], sin(t*pi/600))
+
+        # Check non-exact
+        t = 90 #Halfway between 60 and 120
+        q = F(t-delta)
+        assert num.allclose( (120**2 + 60**2)/2, q[1] )
+        assert num.allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
+
+
+        t = 100 # Two thirds of the way between between 60 and 120
+        q = F(t-delta)
+        assert num.allclose( 2*120**2/3 + 60**2/3, q[1] )
+        assert num.allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
+
+
+        os.remove(filename + '.tms')
+        os.remove(filename + '.txt')                
+
+        
+        
+        
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -48 +48 @@
                   interpolation_points=None,
                   time_thinning=1,
+                  time_limit=None,
                   verbose=False,
                   use_cache=False,
@@ -132 +133 @@
               'interpolation_points': interpolation_points,
               'domain_starttime': domain_starttime,
-              'time_thinning': time_thinning,                   
+              'time_thinning': time_thinning,      
+              'time_limit': time_limit,                                 
               'verbose': verbose,
               'boundary_polygon': boundary_polygon}
@@ -195 +197 @@
                    domain_starttime=None,
                    time_thinning=1,
+                   time_limit=None,
                    verbose=False,
                    boundary_polygon=None):
@@ -221 +224 @@
                                         domain_starttime,
                                         time_thinning=time_thinning,
+                                        time_limit=time_limit,
                                         verbose=verbose,
                                         boundary_polygon=boundary_polygon)
@@ -246 +250 @@
                              interpolation_points=None,
                              domain_starttime=None,                            
-                             time_thinning=1,                             
+                             time_thinning=1,                 
+                             time_limit=None,            
                              verbose=False,
                              boundary_polygon=None):
@@ -329 +334 @@
         starttime = fid.starttime[0]
     except ValueError:
-        msg = 'Could not read starttime from file %s' %filename
+        msg = 'Could not read starttime from file %s' % filename
         raise msg
 
@@ -336 +341 @@
     time = fid.variables['time'][:]    
 
+    # FIXME(Ole): Is time monotoneous?
+    
+    # Apply time limit if requested
+    upper_time_index = len(time)    
+    msg = 'Time vector obtained from file %s has length 0' % filename
+    assert upper_time_index > 0, msg
+    
+    if time_limit is not None:
+        for i, t in enumerate(time):
+            if t > time_limit:
+                upper_time_index = i
+                break
+                
+        msg = 'Time vector is zero. Requested time limit is %f' % time_limit
+        assert upper_time_index > 0, msg                
+
+    time = time[:upper_time_index]
+
+
+    
+    
     # Get time independent stuff
     if spatial:
@@ -345 +371 @@
         x = fid.variables['x'][:]
         y = fid.variables['y'][:]
-        if filename[-3:] == 'sww':
+        if filename.endswith('sww'):
             triangles = fid.variables['volumes'][:]
 
@@ -363 +389 @@
                 for j in range(len(x)):
                     if num.allclose(vertex_coordinates[j],boundary_polygon[i],1e-4):
-                        #FIX ME:
+                        #FIXME:
                         #currently gauges lat and long is stored as float and
                         #then cast to double. This cuases slight repositioning

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -1028 +1028 @@
                                                       point_coordinates=\
                                                       self.interpolation_points,
-                                                      verbose=False) #No clutter
+                                                      verbose=False) # No clutter
                     elif triangles is None and vertex_coordinates is not None:
                         result = \

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -1297 +1297 @@
                  mean_stage=0.0,
                  time_thinning=1,
+                 time_limit=None,
                  boundary_polygon=None,    
                  default_boundary=None,                  
@@ -1330 +1331 @@
         self.file_boundary = File_boundary(filename, domain,
                                            time_thinning=time_thinning,
+                                           time_limit=time_limit,
                                            boundary_polygon=boundary_polygon,
                                            default_boundary=default_boundary,

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -8244 +8244 @@
 
 
-
-        
-        
-        
-            
-
         
         domain_drchlt = Domain(meshname)
@@ -8845 +8839 @@
         This one uses a sine wave and compares to time boundary
         """
-        #FIXME (Ole): Under construction
         
         from anuga.shallow_water import Domain
@@ -8994 +8987 @@
         domain_time.set_quantity('stage', tide)
         Br = Reflective_boundary(domain_time)
-        Bw=Time_boundary(domain=domain_time,
+        Bw = Time_boundary(domain=domain_time,
                          f=lambda t: [num.sin(t)+tide,3.*(20.+num.sin(t)+tide),2.*(20.+num.sin(t)+tide)])
         domain_time.set_boundary({'ocean': Bw,'otherocean': Br})
@@ -9032 +9025 @@
         
         
-        
+
+        
+        
+    def test_file_boundary_sts_time_limit(self):
+        """test_file_boundary_stsIV_sinewave_ordering(self):
+        Read correct points from ordering file and apply sts to boundary
+        This one uses a sine wave and compares to time boundary
+        
+        This one tests that times used can be limited by upper limit
+        """
+        
+        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
+
+        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
+        tide = 0.35
+        time_step_count = 50
+        time_step = 0.1
+        times_ref = num.arange(0, time_step_count*time_step, time_step)
+        
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.Int)
+        last_tstep=(time_step_count)*num.ones(n,num.Int)
+        
+        gauge_depth=20*num.ones(n,num.Float)
+        
+        ha1=num.ones((n,time_step_count),num.Float)
+        ua1=3.*num.ones((n,time_step_count),num.Float)
+        va1=2.*num.ones((n,time_step_count),num.Float)
+        for i in range(n):
+            ha1[i]=num.sin(times_ref)
+        
+        
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha1,
+                                           ua=ua1,
+                                           va=va1)
+
+        # Write order file
+        file_handle, order_base_name = tempfile.mkstemp("")
+        os.close(file_handle)
+        os.remove(order_base_name)
+        d=","
+        order_file=order_base_name+'order.txt'
+        fid=open(order_file,'w')
+        
+        # Write Header
+        header='index, longitude, latitude\n'
+        fid.write(header)
+        indices=[3,0,1]
+        for i in indices:
+            line=str(i)+d+str(lat_long_points[i][1])+d+\
+                str(lat_long_points[i][0])+"\n"
+            fid.write(line)
+        fid.close()
+
+        sts_file=base_name
+        urs2sts(base_name, basename_out=sts_file,
+                ordering_filename=order_file,
+                mean_stage=tide,
+                verbose=False)
+        self.delete_mux(files)
+        
+        
+        
+        # Now read the sts file and check that values have been stored correctly.
+        fid = NetCDFFile(sts_file + '.sts')
+
+        # Check the time vector
+        times = fid.variables['time'][:]
+        
+        #print times
+
+        # Check sts quantities
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+
+        
+
+        # Create beginnings of boundary polygon based on sts_boundary
+        boundary_polygon = create_sts_boundary(base_name)
+        
+        os.remove(order_file)
+
+        # Append the remaining part of the boundary polygon to be defined by
+        # the user
+        bounding_polygon_utm=[]
+        for point in bounding_polygon:
+            zone,easting,northing=redfearn(point[0],point[1])
+            bounding_polygon_utm.append([easting,northing])
+
+        boundary_polygon.append(bounding_polygon_utm[3])
+        boundary_polygon.append(bounding_polygon_utm[4])
+
+        #print 'boundary_polygon', boundary_polygon
+        
+
+        assert num.allclose(bounding_polygon_utm,boundary_polygon)
+
+
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        
+        # have to change boundary tags from last example because now bounding
+        # polygon starts in different place.
+        create_mesh_from_regions(boundary_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)
+        
+        
+        Bf = File_boundary(sts_file+'.sts', 
+                           domain_fbound, 
+                           boundary_polygon=boundary_polygon)
+        time_vec = Bf.F.get_time()
+        assert num.allclose(time_vec, times_ref)                                   
+        
+        for time_limit in [0.1, 0.2, 0.5, 1.0, 2.2, 3.0, 4.3, 6.0, 10.0]:
+            Bf = File_boundary(sts_file+'.sts', 
+                               domain_fbound, 
+                               time_limit=time_limit,
+                               boundary_polygon=boundary_polygon)
+        
+            time_vec = Bf.F.get_time()
+            assert num.alltrue(time_vec < time_limit)
+            
+            
+        try:    
+            Bf = File_boundary(sts_file+'.sts', 
+                               domain_fbound, 
+                               time_limit=-1,
+                               boundary_polygon=boundary_polygon)            
+            time_vec = Bf.F.get_time()    
+            print time_vec    
+        except AssertionError:
+            pass
+        else:
+            raise Exception, 'Should have raised Exception here'
 
     def test_lon_lat2grid(self):
@@ -11103 +11249 @@
 
     suite = unittest.makeSuite(Test_Data_Manager,'test')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV_sinewave_ordering')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_sts_time_limit')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
     #suite = unittest.makeSuite(Test_Data_Manager,'covered_')