Changeset 6173


Ignore:
Timestamp:
Jan 14, 2009, 5:47:17 PM (10 years ago)
Author:
ole
Message:

Introduced time_limit in Field_boundary, File_boundary and file_function

Location:
anuga_core/source/anuga
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6166 r6173  
    158158    def __init__(self, filename, domain,
    159159                 time_thinning=1,
     160                 time_limit=None,
    160161                 boundary_polygon=None,   
    161162                 default_boundary=None,
     
    219220                               interpolation_points=self.midpoint_coordinates,
    220221                               time_thinning=time_thinning,
     222                               time_limit=time_limit,
    221223                               use_cache=use_cache,
    222224                               verbose=verbose,
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

    r6171 r6173  
    10521052        assert num.allclose(domain.starttime, start+delta)
    10531053
    1054 
     1054        assert num.allclose(F.get_time(), [-23., 37., 97., 157., 217.,
     1055                                            277., 337., 397., 457., 517.,
     1056                                            577., 637., 697., 757., 817.,
     1057                                            877., 937., 997., 1057., 1117.,
     1058                                            1177.])
    10551059
    10561060
     
    10831087        os.remove(filename + '.txt')               
    10841088
     1089       
     1090
     1091    def test_file_function_time_with_domain_different_start_and_time_limit(self):
     1092        """Test that File function interpolates correctly
     1093        between given times. No x,y dependency here.
     1094        Use domain with a starttime later than that of file
     1095
     1096        ASCII version
     1097       
     1098        This test also tests that time can be truncated.
     1099        """
     1100
     1101        # Write file
     1102        import os, time, calendar
     1103        from anuga.config import time_format
     1104        from math import sin, pi
     1105        from domain import Domain
     1106
     1107        finaltime = 1200
     1108        filename = 'test_file_function'
     1109        fid = open(filename + '.txt', 'w')
     1110        start = time.mktime(time.strptime('2000', '%Y'))
     1111        dt = 60  #One minute intervals
     1112        t = 0.0
     1113        while t <= finaltime:
     1114            t_string = time.strftime(time_format, time.gmtime(t+start))
     1115            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
     1116            t += dt
     1117
     1118        fid.close()
     1119
     1120        # Convert ASCII file to NetCDF (Which is what we really like!)
     1121        timefile2netcdf(filename)       
     1122
     1123        a = [0.0, 0.0]
     1124        b = [4.0, 0.0]
     1125        c = [0.0, 3.0]
     1126
     1127        points = [a, b, c]
     1128        vertices = [[0,1,2]]
     1129        domain = Domain(points, vertices)
     1130
     1131        # Check that domain.starttime isn't updated if later than file starttime but earlier
     1132        # than file end time
     1133        delta = 23
     1134        domain.starttime = start + delta
     1135        F = file_function(filename + '.tms', domain,
     1136                          time_limit=600,
     1137                          quantities=['Attribute0', 'Attribute1', 'Attribute2'])       
     1138        assert num.allclose(domain.starttime, start+delta)
     1139
     1140        assert num.allclose(F.get_time(), [-23., 37., 97., 157., 217.,
     1141                                            277., 337., 397., 457., 517.,
     1142                                            577.])       
     1143
     1144
     1145
     1146        # Now try interpolation with delta offset
     1147        for i in range(20):
     1148            t = i*10
     1149            q = F(t-delta)
     1150
     1151            #Exact linear intpolation
     1152            assert num.allclose(q[0], 2*t)
     1153            if i%6 == 0:
     1154                assert num.allclose(q[1], t**2)
     1155                assert num.allclose(q[2], sin(t*pi/600))
     1156
     1157        # Check non-exact
     1158        t = 90 #Halfway between 60 and 120
     1159        q = F(t-delta)
     1160        assert num.allclose( (120**2 + 60**2)/2, q[1] )
     1161        assert num.allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
     1162
     1163
     1164        t = 100 # Two thirds of the way between between 60 and 120
     1165        q = F(t-delta)
     1166        assert num.allclose( 2*120**2/3 + 60**2/3, q[1] )
     1167        assert num.allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
     1168
     1169
     1170        os.remove(filename + '.tms')
     1171        os.remove(filename + '.txt')               
     1172
     1173       
     1174       
     1175       
    10851176
    10861177
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r6171 r6173  
    4848                  interpolation_points=None,
    4949                  time_thinning=1,
     50                  time_limit=None,
    5051                  verbose=False,
    5152                  use_cache=False,
     
    132133              'interpolation_points': interpolation_points,
    133134              'domain_starttime': domain_starttime,
    134               'time_thinning': time_thinning,                   
     135              'time_thinning': time_thinning,     
     136              'time_limit': time_limit,                                 
    135137              'verbose': verbose,
    136138              'boundary_polygon': boundary_polygon}
     
    195197                   domain_starttime=None,
    196198                   time_thinning=1,
     199                   time_limit=None,
    197200                   verbose=False,
    198201                   boundary_polygon=None):
     
    221224                                        domain_starttime,
    222225                                        time_thinning=time_thinning,
     226                                        time_limit=time_limit,
    223227                                        verbose=verbose,
    224228                                        boundary_polygon=boundary_polygon)
     
    246250                             interpolation_points=None,
    247251                             domain_starttime=None,                           
    248                              time_thinning=1,                             
     252                             time_thinning=1,                 
     253                             time_limit=None,           
    249254                             verbose=False,
    250255                             boundary_polygon=None):
     
    329334        starttime = fid.starttime[0]
    330335    except ValueError:
    331         msg = 'Could not read starttime from file %s' %filename
     336        msg = 'Could not read starttime from file %s' % filename
    332337        raise msg
    333338
     
    336341    time = fid.variables['time'][:]   
    337342
     343    # FIXME(Ole): Is time monotoneous?
     344   
     345    # Apply time limit if requested
     346    upper_time_index = len(time)   
     347    msg = 'Time vector obtained from file %s has length 0' % filename
     348    assert upper_time_index > 0, msg
     349   
     350    if time_limit is not None:
     351        for i, t in enumerate(time):
     352            if t > time_limit:
     353                upper_time_index = i
     354                break
     355               
     356        msg = 'Time vector is zero. Requested time limit is %f' % time_limit
     357        assert upper_time_index > 0, msg               
     358
     359    time = time[:upper_time_index]
     360
     361
     362   
     363   
    338364    # Get time independent stuff
    339365    if spatial:
     
    345371        x = fid.variables['x'][:]
    346372        y = fid.variables['y'][:]
    347         if filename[-3:] == 'sww':
     373        if filename.endswith('sww'):
    348374            triangles = fid.variables['volumes'][:]
    349375
     
    363389                for j in range(len(x)):
    364390                    if num.allclose(vertex_coordinates[j],boundary_polygon[i],1e-4):
    365                         #FIX ME:
     391                        #FIXME:
    366392                        #currently gauges lat and long is stored as float and
    367393                        #then cast to double. This cuases slight repositioning
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r6152 r6173  
    10281028                                                      point_coordinates=\
    10291029                                                      self.interpolation_points,
    1030                                                       verbose=False) #No clutter
     1030                                                      verbose=False) # No clutter
    10311031                    elif triangles is None and vertex_coordinates is not None:
    10321032                        result = \
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r6166 r6173  
    12971297                 mean_stage=0.0,
    12981298                 time_thinning=1,
     1299                 time_limit=None,
    12991300                 boundary_polygon=None,   
    13001301                 default_boundary=None,                 
     
    13301331        self.file_boundary = File_boundary(filename, domain,
    13311332                                           time_thinning=time_thinning,
     1333                                           time_limit=time_limit,
    13321334                                           boundary_polygon=boundary_polygon,
    13331335                                           default_boundary=default_boundary,
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r6157 r6173  
    82448244
    82458245
    8246 
    8247        
    8248        
    8249        
    8250            
    8251 
    82528246       
    82538247        domain_drchlt = Domain(meshname)
     
    88458839        This one uses a sine wave and compares to time boundary
    88468840        """
    8847         #FIXME (Ole): Under construction
    88488841       
    88498842        from anuga.shallow_water import Domain
     
    89948987        domain_time.set_quantity('stage', tide)
    89958988        Br = Reflective_boundary(domain_time)
    8996         Bw=Time_boundary(domain=domain_time,
     8989        Bw = Time_boundary(domain=domain_time,
    89978990                         f=lambda t: [num.sin(t)+tide,3.*(20.+num.sin(t)+tide),2.*(20.+num.sin(t)+tide)])
    89988991        domain_time.set_boundary({'ocean': Bw,'otherocean': Br})
     
    90329025       
    90339026       
    9034        
     9027
     9028       
     9029       
     9030    def test_file_boundary_sts_time_limit(self):
     9031        """test_file_boundary_stsIV_sinewave_ordering(self):
     9032        Read correct points from ordering file and apply sts to boundary
     9033        This one uses a sine wave and compares to time boundary
     9034       
     9035        This one tests that times used can be limited by upper limit
     9036        """
     9037       
     9038        from anuga.shallow_water import Domain
     9039        from anuga.shallow_water import Reflective_boundary
     9040        from anuga.shallow_water import Dirichlet_boundary
     9041        from anuga.shallow_water import File_boundary
     9042        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     9043
     9044        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
     9045        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
     9046        tide = 0.35
     9047        time_step_count = 50
     9048        time_step = 0.1
     9049        times_ref = num.arange(0, time_step_count*time_step, time_step)
     9050       
     9051        n=len(lat_long_points)
     9052        first_tstep=num.ones(n,num.Int)
     9053        last_tstep=(time_step_count)*num.ones(n,num.Int)
     9054       
     9055        gauge_depth=20*num.ones(n,num.Float)
     9056       
     9057        ha1=num.ones((n,time_step_count),num.Float)
     9058        ua1=3.*num.ones((n,time_step_count),num.Float)
     9059        va1=2.*num.ones((n,time_step_count),num.Float)
     9060        for i in range(n):
     9061            ha1[i]=num.sin(times_ref)
     9062       
     9063       
     9064        base_name, files = self.write_mux2(lat_long_points,
     9065                                           time_step_count, time_step,
     9066                                           first_tstep, last_tstep,
     9067                                           depth=gauge_depth,
     9068                                           ha=ha1,
     9069                                           ua=ua1,
     9070                                           va=va1)
     9071
     9072        # Write order file
     9073        file_handle, order_base_name = tempfile.mkstemp("")
     9074        os.close(file_handle)
     9075        os.remove(order_base_name)
     9076        d=","
     9077        order_file=order_base_name+'order.txt'
     9078        fid=open(order_file,'w')
     9079       
     9080        # Write Header
     9081        header='index, longitude, latitude\n'
     9082        fid.write(header)
     9083        indices=[3,0,1]
     9084        for i in indices:
     9085            line=str(i)+d+str(lat_long_points[i][1])+d+\
     9086                str(lat_long_points[i][0])+"\n"
     9087            fid.write(line)
     9088        fid.close()
     9089
     9090        sts_file=base_name
     9091        urs2sts(base_name, basename_out=sts_file,
     9092                ordering_filename=order_file,
     9093                mean_stage=tide,
     9094                verbose=False)
     9095        self.delete_mux(files)
     9096       
     9097       
     9098       
     9099        # Now read the sts file and check that values have been stored correctly.
     9100        fid = NetCDFFile(sts_file + '.sts')
     9101
     9102        # Check the time vector
     9103        times = fid.variables['time'][:]
     9104       
     9105        #print times
     9106
     9107        # Check sts quantities
     9108        stage = fid.variables['stage'][:]
     9109        xmomentum = fid.variables['xmomentum'][:]
     9110        ymomentum = fid.variables['ymomentum'][:]
     9111        elevation = fid.variables['elevation'][:]
     9112
     9113       
     9114
     9115        # Create beginnings of boundary polygon based on sts_boundary
     9116        boundary_polygon = create_sts_boundary(base_name)
     9117       
     9118        os.remove(order_file)
     9119
     9120        # Append the remaining part of the boundary polygon to be defined by
     9121        # the user
     9122        bounding_polygon_utm=[]
     9123        for point in bounding_polygon:
     9124            zone,easting,northing=redfearn(point[0],point[1])
     9125            bounding_polygon_utm.append([easting,northing])
     9126
     9127        boundary_polygon.append(bounding_polygon_utm[3])
     9128        boundary_polygon.append(bounding_polygon_utm[4])
     9129
     9130        #print 'boundary_polygon', boundary_polygon
     9131       
     9132
     9133        assert num.allclose(bounding_polygon_utm,boundary_polygon)
     9134
     9135
     9136        extent_res=1000000
     9137        meshname = 'urs_test_mesh' + '.tsh'
     9138        interior_regions=None
     9139        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
     9140       
     9141        # have to change boundary tags from last example because now bounding
     9142        # polygon starts in different place.
     9143        create_mesh_from_regions(boundary_polygon,
     9144                                 boundary_tags=boundary_tags,
     9145                                 maximum_triangle_area=extent_res,
     9146                                 filename=meshname,
     9147                                 interior_regions=interior_regions,
     9148                                 verbose=False)
     9149       
     9150        domain_fbound = Domain(meshname)
     9151        domain_fbound.set_quantity('stage', tide)
     9152       
     9153       
     9154        Bf = File_boundary(sts_file+'.sts',
     9155                           domain_fbound,
     9156                           boundary_polygon=boundary_polygon)
     9157        time_vec = Bf.F.get_time()
     9158        assert num.allclose(time_vec, times_ref)                                   
     9159       
     9160        for time_limit in [0.1, 0.2, 0.5, 1.0, 2.2, 3.0, 4.3, 6.0, 10.0]:
     9161            Bf = File_boundary(sts_file+'.sts',
     9162                               domain_fbound,
     9163                               time_limit=time_limit,
     9164                               boundary_polygon=boundary_polygon)
     9165       
     9166            time_vec = Bf.F.get_time()
     9167            assert num.alltrue(time_vec < time_limit)
     9168           
     9169           
     9170        try:   
     9171            Bf = File_boundary(sts_file+'.sts',
     9172                               domain_fbound,
     9173                               time_limit=-1,
     9174                               boundary_polygon=boundary_polygon)           
     9175            time_vec = Bf.F.get_time()   
     9176            print time_vec   
     9177        except AssertionError:
     9178            pass
     9179        else:
     9180            raise Exception, 'Should have raised Exception here'
    90359181
    90369182    def test_lon_lat2grid(self):
     
    1110311249
    1110411250    suite = unittest.makeSuite(Test_Data_Manager,'test')
    11105     #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV_sinewave_ordering')
     11251    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_sts_time_limit')
    1110611252    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
    1110711253    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
Note: See TracChangeset for help on using the changeset viewer.