Changeset 5544


Ignore:
Timestamp:
Jul 22, 2008, 4:44:03 AM (11 years ago)
Author:
ole
Message:

More testing of urs2sts and sts file boundary.
One test under construction and some issues were found with the initial tide level

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5541 r5544  
    59775977        os.remove(sww_file)
    59785978       
    5979     def write_mux2(self,lat_long_points, time_step_count, time_step,
     5979    def write_mux2(self, lat_long_points, time_step_count, time_step,
    59805980                   first_tstep, last_tstep,
    59815981                   depth=None, ha=None, ua=None, va=None):
     
    65476547        """
    65486548       
    6549         tide=0
    6550         time_step_count = 3
     6549        tide = 0.35
     6550        time_step_count = 6 # I made this a little longer (Ole)
    65516551        time_step = 2
    65526552        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
     
    66826682
    66836683       
    6684         assert allclose(2.0*transpose(ha_permutation), stage)  # Meters
     6684        assert allclose(2.0*transpose(ha_permutation)+tide, stage)  # Meters
    66856685
    66866686        #Check the momentums - ua
     
    67946794
    67956795       
    6796        
    6797        
    6798     def test_urs2sts_file_boundary_stsI(self):
     6796
     6797       
     6798    def test_urs2sts_ordering_different_sources(self):
     6799        """Test multiple sources with ordering file, different source files and weights.
     6800           This test also has more variable values than the previous ones
     6801        """
     6802       
     6803        from Numeric import sin, cos
     6804               
     6805        tide = 1.5
     6806        time_step_count = 10
     6807        time_step = 0.2
     6808       
     6809        times_ref = arange(0, time_step_count*time_step, time_step)
     6810        #print 'time vector', times_ref
     6811       
     6812        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
     6813        n = len(lat_long_points)
     6814       
     6815        # Create non-trivial weights
     6816        weights = [0.8, 1.5]
     6817       
     6818        # Create different timeseries starting and ending at different times
     6819        first_tstep=ones(n,Int)
     6820        first_tstep[0]+=2   # Point 0 starts at 2
     6821        first_tstep[1]+=4   # Point 1 starts at 4       
     6822        first_tstep[2]+=3   # Point 2 starts at 3
     6823       
     6824        last_tstep=(time_step_count)*ones(n,Int)
     6825        last_tstep[0]-=1    # Point 0 ends 1 step early
     6826        last_tstep[1]-=2    # Point 1 ends 2 steps early               
     6827        last_tstep[4]-=3    # Point 4 ends 3 steps early       
     6828       
     6829        #print
     6830        #print 'time_step_count', time_step_count
     6831        #print 'time_step', time_step
     6832        #print 'first_tstep', first_tstep
     6833        #print 'last_tstep', last_tstep               
     6834       
     6835       
     6836        # Create varying elevation data (positive values for seafloor)
     6837        gauge_depth=20*ones(n,Float)
     6838        for i in range(n):
     6839            gauge_depth[i] += i**2
     6840           
     6841        #print 'gauge_depth', gauge_depth
     6842       
     6843        # Create data to be written to first mux file       
     6844        ha0=2*ones((n,time_step_count),Float)
     6845        ha0[0]=arange(0,time_step_count)
     6846        ha0[1]=arange(time_step_count,2*time_step_count)
     6847        ha0[2]=arange(2*time_step_count,3*time_step_count)
     6848        ha0[3]=arange(3*time_step_count,4*time_step_count)
     6849        ua0=5*ones((n,time_step_count),Float)
     6850        va0=-10*ones((n,time_step_count),Float)
     6851
     6852        # Ensure data used to write mux file to be zero when gauges are
     6853        # not recording
     6854        for i in range(n):
     6855             # For each point
     6856             
     6857             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
     6858                 # For timesteps before and after recording range
     6859                 ha0[i][j] = ua0[i][j] = va0[i][j] = 0.0                                 
     6860
     6861
     6862        #print
     6863        #print 'using varying start and end time'
     6864        #print 'ha0', ha0
     6865        #print 'ua0', ua0
     6866        #print 'va0', va0       
     6867       
     6868        # Write first mux file to be combined by urs2sts
     6869        base_nameI, filesI = self.write_mux2(lat_long_points,
     6870                                             time_step_count, time_step,
     6871                                             first_tstep, last_tstep,
     6872                                             depth=gauge_depth,
     6873                                             ha=ha0,
     6874                                             ua=ua0,
     6875                                             va=va0)
     6876
     6877                                             
     6878                                             
     6879        # Create data to be written to second mux file       
     6880        ha1=zeros((n,time_step_count),Float)
     6881        ha1[0]=sin(times_ref)
     6882        ha1[1]=2*sin(times_ref - 3)
     6883        ha1[2]=5*sin(4*times_ref)
     6884        ha1[3]=sin(times_ref)
     6885        ha1[4]=sin(2*times_ref-0.7)
     6886               
     6887        ua1=zeros((n,time_step_count),Float)
     6888        ua1[0]=3*cos(times_ref)       
     6889        ua1[1]=2*sin(times_ref-0.7)   
     6890        ua1[2]=arange(3*time_step_count,4*time_step_count)
     6891        ua1[4]=2*ones(time_step_count)
     6892       
     6893        va1=zeros((n,time_step_count),Float)
     6894        va1[0]=2*cos(times_ref-0.87)       
     6895        va1[1]=3*ones(time_step_count)
     6896        va1[3]=2*sin(times_ref-0.71)       
     6897       
     6898       
     6899        # Ensure data used to write mux file to be zero when gauges are
     6900        # not recording
     6901        for i in range(n):
     6902             # For each point
     6903             
     6904             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
     6905                 # For timesteps before and after recording range
     6906                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0                                 
     6907
     6908
     6909        #print
     6910        #print 'using varying start and end time'
     6911        #print 'ha1', ha1
     6912        #print 'ua1', ua1
     6913        #print 'va1', va1       
     6914                                             
     6915                                             
     6916        # Write second mux file to be combined by urs2sts                                             
     6917        base_nameII, filesII = self.write_mux2(lat_long_points,
     6918                                               time_step_count, time_step,
     6919                                               first_tstep, last_tstep,
     6920                                               depth=gauge_depth,
     6921                                               ha=ha1,
     6922                                               ua=ua1,
     6923                                               va=va1)
     6924
     6925                                               
     6926        # Create ordering file
     6927        permutation = [4,0,2]
     6928
     6929        _, ordering_filename = tempfile.mkstemp('')
     6930        order_fid = open(ordering_filename, 'w') 
     6931        order_fid.write('index, longitude, latitude\n')
     6932        for index in permutation:
     6933            order_fid.write('%d, %f, %f\n' %(index,
     6934                                             lat_long_points[index][1],
     6935                                             lat_long_points[index][0]))
     6936        order_fid.close()
     6937       
     6938           
     6939
     6940                                               
     6941        # Call urs2sts with multiple mux files
     6942        urs2sts([base_nameI, base_nameII],
     6943                basename_out=base_nameI,
     6944                ordering_filename=ordering_filename,
     6945                weights=weights,
     6946                mean_stage=tide,
     6947                verbose=False)
     6948
     6949        # Now read the sts file and check that values have been stored correctly.
     6950        sts_file = base_nameI + '.sts'
     6951        fid = NetCDFFile(sts_file)
     6952
     6953        # Make x and y absolute
     6954        x = fid.variables['x'][:]
     6955        y = fid.variables['y'][:]
     6956
     6957        geo_reference = Geo_reference(NetCDFObject=fid)
     6958        points = geo_reference.get_absolute(map(None, x, y))
     6959        points = ensure_numeric(points)
     6960
     6961        x = points[:,0]
     6962        y = points[:,1]
     6963
     6964        for i, index in enumerate(permutation):
     6965            # Check that STS points are stored in the correct order
     6966           
     6967            # Work out the UTM coordinates sts point i
     6968            zone, e, n = redfearn(lat_long_points[index][0],
     6969                                  lat_long_points[index][1])             
     6970
     6971            #print i, [x[i],y[i]], [e,n]
     6972            assert allclose([x[i],y[i]], [e,n])
     6973           
     6974                       
     6975        # Check the time vector
     6976        times = fid.variables['time'][:]
     6977        assert allclose(ensure_numeric(times),
     6978                        ensure_numeric(times_ref))
     6979                       
     6980
     6981        # Check sts values
     6982        stage = fid.variables['stage'][:]
     6983        xmomentum = fid.variables['xmomentum'][:]
     6984        ymomentum = fid.variables['ymomentum'][:]
     6985        elevation = fid.variables['elevation'][:]
     6986
     6987       
     6988        # The quantities stored in the .sts file should be the weighted sum of the
     6989        # quantities written to the mux2 files subject to the permutation vector.
     6990       
     6991        ha_ref = weights[0]*take(ha0, permutation) + weights[1]*take(ha1, permutation)
     6992        ua_ref = weights[0]*take(ua0, permutation) + weights[1]*take(ua1, permutation)       
     6993        va_ref = weights[0]*take(va0, permutation) + weights[1]*take(va1, permutation)               
     6994
     6995        gauge_depth_ref = take(gauge_depth, permutation)                         
     6996
     6997       
     6998        assert allclose(transpose(ha_ref)+tide, stage)  # Meters
     6999
     7000        #Check the momentums - ua
     7001        #momentum = velocity*(stage-elevation)
     7002        # elevation = - depth
     7003        #momentum = velocity_ua *(stage+depth)
     7004
     7005        depth_ref = zeros((len(permutation), time_step_count), Float)
     7006        for i in range(len(permutation)):
     7007            depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i]
     7008
     7009           
     7010       
     7011
     7012        # The xmomentum stored in the .sts file should be the sum of the ua
     7013        # in the two mux2 files multiplied by the depth.
     7014        assert allclose(transpose(ua_ref*depth_ref), xmomentum)
     7015
     7016        #Check the momentums - va
     7017        #momentum = velocity*(stage-elevation)
     7018        # elevation = - depth
     7019        #momentum = velocity_va *(stage+depth)
     7020
     7021        # The ymomentum stored in the .sts file should be the sum of the va
     7022        # in the two mux2 files multiplied by the depth.
     7023        assert allclose(transpose(va_ref*depth_ref), ymomentum)
     7024
     7025        # check the elevation values.
     7026        # -ve since urs measures depth, sww meshers height,
     7027        assert allclose(-gauge_depth_ref, elevation)  #Meters
     7028
     7029        fid.close()
     7030        self.delete_mux(filesI)
     7031        self.delete_mux(filesII)
     7032        os.remove(sts_file)
     7033       
     7034
     7035       
     7036       
     7037       
     7038       
     7039               
     7040       
     7041    def test_file_boundary_stsI(self):
     7042        """test_file_boundary_stsI(self):
     7043        """
     7044       
    67997045        from anuga.shallow_water import Domain
    68007046        from anuga.shallow_water import Reflective_boundary
     
    68027048        from anuga.shallow_water import File_boundary
    68037049        from anuga.pmesh.mesh_interface import create_mesh_from_regions
    6804         from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
     7050
    68057051        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
    6806         tide=0.
     7052        tide = 0.37
    68077053        time_step_count = 5
    68087054        time_step = 2
     
    68117057        first_tstep=ones(n,Int)
    68127058        last_tstep=(time_step_count)*ones(n,Int)
    6813         gauge_depth=20*ones(n,Float)
    6814         ha=2*ones((n,time_step_count),Float)
    6815         ua=10*ones((n,time_step_count),Float)
    6816         va=-10*ones((n,time_step_count),Float)
     7059
     7060        h = 20       
     7061        w = 2
     7062        u = 10
     7063        v = -10
     7064        gauge_depth=h*ones(n,Float)
     7065        ha=w*ones((n,time_step_count),Float)
     7066        ua=u*ones((n,time_step_count),Float)
     7067        va=v*ones((n,time_step_count),Float)
    68177068        base_name, files = self.write_mux2(lat_long_points,
    6818                                    time_step_count, time_step,
    6819                                    first_tstep, last_tstep,
    6820                                    depth=gauge_depth,
    6821                                    ha=ha,
    6822                                    ua=ua,
    6823                                    va=va)
     7069                                           time_step_count, time_step,
     7070                                           first_tstep, last_tstep,
     7071                                           depth=gauge_depth,
     7072                                           ha=ha,
     7073                                           ua=ua,
     7074                                           va=va)
    68247075
    68257076        sts_file=base_name
    68267077        urs2sts(base_name,
    68277078                sts_file,
    6828                 mean_stage=tide,verbose=False)
     7079                mean_stage=tide,
     7080                verbose=False)
    68297081        self.delete_mux(files)
    68307082
     
    68407092                         interior_regions=interior_regions,verbose=False)
    68417093       
    6842         domain_fbound = pmesh_to_domain_instance(meshname, Domain)
     7094        domain_fbound = Domain(meshname)
    68437095        domain_fbound.set_quantity('stage', tide)
    68447096        Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon)
    68457097        Br = Reflective_boundary(domain_fbound)
    6846         Bd=Dirichlet_boundary([2.0,220,-220])
     7098
    68477099        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
    68487100        finaltime=time_step*(time_step_count-1)
    68497101        yieldstep=time_step
    68507102        temp_fbound=zeros(int(finaltime/yieldstep)+1,Float)
    6851         i=0
    6852         for t in domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,
    6853                                       skip_initial_step = False):
     7103
     7104        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,
     7105                                                   skip_initial_step = False)):
    68547106            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
    6855             i+=1
    6856        
    6857         domain_drchlt = pmesh_to_domain_instance(meshname, Domain)
     7107
     7108       
     7109        domain_drchlt = Domain(meshname)
    68587110        domain_drchlt.set_quantity('stage', tide)
    68597111        Br = Reflective_boundary(domain_drchlt)
    6860         Bd=Dirichlet_boundary([2.0,220,-220])
     7112        Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)])
    68617113        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
    68627114        temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float)
    6863         i=0
    6864         for t in domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,
    6865                                       skip_initial_step = False):
     7115
     7116        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,
     7117                                                   skip_initial_step = False)):
    68667118            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
    6867             i+=1
    6868 
     7119
     7120        #print domain_fbound.quantities['stage'].vertex_values
     7121        #print domain_drchlt.quantities['stage'].vertex_values
     7122                   
    68697123        assert allclose(temp_fbound,temp_drchlt)
     7124       
     7125        assert allclose(domain_fbound.quantities['stage'].vertex_values,
     7126                        domain_drchlt.quantities['stage'].vertex_values)
     7127                       
     7128        assert allclose(domain_fbound.quantities['xmomentum'].vertex_values,
     7129                        domain_drchlt.quantities['xmomentum'].vertex_values)                       
     7130                       
     7131        assert allclose(domain_fbound.quantities['ymomentum'].vertex_values,
     7132                        domain_drchlt.quantities['ymomentum'].vertex_values)                                               
     7133       
     7134       
    68707135        os.remove(sts_file+'.sts')
    68717136        os.remove(meshname)
    68727137
    6873     def test_urs2sts_file_boundary_stsII(self):
    6874         """ mux2 file has points not included in boundary
     7138    def test_file_boundary_stsII(self):
     7139        """test_file_boundary_stsII(self):
     7140       
     7141         mux2 file has points not included in boundary
    68757142         mux2 gauges are not stored with the same order as they are
    68767143         found in bounding_polygon. This does not matter as long as bounding
     
    68787145         the correct order).
    68797146         """
     7147         
    68807148        from anuga.shallow_water import Domain
    68817149        from anuga.shallow_water import Reflective_boundary
     
    68837151        from anuga.shallow_water import File_boundary
    68847152        from anuga.pmesh.mesh_interface import create_mesh_from_regions
    6885         from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
     7153
    68867154        bounding_polygon=[[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02],[6.0,97.0]]
    6887         tide=0.
    6888         time_step_count = 5
     7155        tide = 0.0 # FIXME (Ole): For some reason, this one has to be zero       
     7156        time_step_count = 20
    68897157        time_step = 2
    68907158        lat_long_points=bounding_polygon[0:2]
     
    69177185        interior_regions=None
    69187186        boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]}
    6919         #have to change boundary tags from last example because now bounding
    6920         #polygon starts in different place.
     7187        # have to change boundary tags from last example because now bounding
     7188        # polygon starts in different place.
    69217189        create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags,
    69227190                         maximum_triangle_area=extent_res,filename=meshname,
    69237191                         interior_regions=interior_regions,verbose=False)
    69247192       
    6925         domain_fbound = pmesh_to_domain_instance(meshname, Domain)
     7193        domain_fbound = Domain(meshname)
    69267194        domain_fbound.set_quantity('stage', tide)
    69277195        Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon)
    69287196        Br = Reflective_boundary(domain_fbound)
    6929         Bd=Dirichlet_boundary([2.0,220,-220])
     7197
    69307198        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
    69317199        finaltime=time_step*(time_step_count-1)
    69327200        yieldstep=time_step
    69337201        temp_fbound=zeros(int(finaltime/yieldstep)+1,Float)
    6934         i=0
    6935         for t in domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,
    6936                                       skip_initial_step = False):
     7202       
     7203        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,
     7204                                                   skip_initial_step = False)):
    69377205            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
    6938             i+=1
    6939        
    6940         domain_drchlt = pmesh_to_domain_instance(meshname, Domain)
     7206       
     7207        domain_drchlt = Domain(meshname)
    69417208        domain_drchlt.set_quantity('stage', tide)
    69427209        Br = Reflective_boundary(domain_drchlt)
    6943         Bd=Dirichlet_boundary([2.0,220,-220])
     7210        Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide])
    69447211        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
    69457212        temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float)
    6946         i=0
    6947         for t in domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,
    6948                                       skip_initial_step = False):
     7213
     7214        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,
     7215                                                   skip_initial_step = False)):
    69497216            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
    6950             i+=1
    6951 
    6952         assert allclose(temp_fbound,temp_drchlt)
     7217
     7218
     7219        assert allclose(temp_fbound,temp_drchlt)           
     7220           
     7221        #print domain_fbound.quantities['stage'].vertex_values
     7222        #print domain_drchlt.quantities['stage'].vertex_values
     7223                   
     7224           
     7225        assert allclose(domain_fbound.quantities['stage'].vertex_values,
     7226                        domain_drchlt.quantities['stage'].vertex_values)
     7227                       
     7228        assert allclose(domain_fbound.quantities['xmomentum'].vertex_values,
     7229                        domain_drchlt.quantities['xmomentum'].vertex_values)                       
     7230                       
     7231        assert allclose(domain_fbound.quantities['ymomentum'].vertex_values,
     7232                        domain_drchlt.quantities['ymomentum'].vertex_values)                                               
     7233       
     7234       
     7235           
     7236           
     7237
    69537238        os.remove(sts_file+'.sts')
    69547239        os.remove(meshname)
    69557240
    6956     def test_urs2sts_file_boundary_stsIII_ordering(self):
    6957         """Read correct points from ordering file
     7241       
     7242       
     7243    def test_file_boundary_stsIII_ordering(self):
     7244        """test_file_boundary_stsIII_ordering(self):
     7245        Read correct points from ordering file and apply sts to boundary
    69587246        """
    69597247        from anuga.shallow_water import Domain
     
    69627250        from anuga.shallow_water import File_boundary
    69637251        from anuga.pmesh.mesh_interface import create_mesh_from_regions
    6964         from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
    69657252
    69667253        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
    69677254        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
    6968         tide=0.
    6969         time_step_count = 5
     7255        tide = 0.0 # FIXME (Ole): For some reason, this one has to be zero
     7256        time_step_count = 50
    69707257        time_step = 2
    69717258        n=len(lat_long_points)
     
    69847271                                   va=va)
    69857272
    6986         #Write order file
     7273        # Write order file
    69877274        file_handle, order_base_name = tempfile.mkstemp("")
    69887275        os.close(file_handle)
     
    69917278        order_file=order_base_name+'order.txt'
    69927279        fid=open(order_file,'w')
    6993         #Write Header
     7280       
     7281        # Write Header
    69947282        header='index, longitude, latitude\n'
    69957283        fid.write(header)
     
    70127300        os.remove(order_file)
    70137301
    7014         #Append the remaining part of the boundary polygon to be defined by
    7015         #the user
     7302        # Append the remaining part of the boundary polygon to be defined by
     7303        # the user
    70167304        bounding_polygon_utm=[]
    70177305        for point in bounding_polygon:
     
    70397327        interior_regions=None
    70407328        boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]}
    7041         #have to change boundary tags from last example because now bounding
    7042         #polygon starts in different place.
     7329       
     7330        # have to change boundary tags from last example because now bounding
     7331        # polygon starts in different place.
    70437332        create_mesh_from_regions(boundary_polygon,boundary_tags=boundary_tags,
    70447333                         maximum_triangle_area=extent_res,filename=meshname,
    70457334                         interior_regions=interior_regions,verbose=False)
    70467335       
    7047         domain_fbound = pmesh_to_domain_instance(meshname, Domain)
     7336        domain_fbound = Domain(meshname)
    70487337        domain_fbound.set_quantity('stage', tide)
    70497338        Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=boundary_polygon)
    70507339        Br = Reflective_boundary(domain_fbound)
    7051         Bd=Dirichlet_boundary([2.0,220,-220])
     7340
    70527341        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
    70537342        finaltime=time_step*(time_step_count-1)
    70547343        yieldstep=time_step
    70557344        temp_fbound=zeros(int(finaltime/yieldstep)+1,Float)
    7056         i=0
    7057         for t in domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,
    7058                                       skip_initial_step = False):
     7345   
     7346        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,
     7347                                                   skip_initial_step = False)):
    70597348            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
    7060             i+=1
    7061        
    7062         domain_drchlt = pmesh_to_domain_instance(meshname, Domain)
     7349   
     7350       
     7351        domain_drchlt = Domain(meshname)
    70637352        domain_drchlt.set_quantity('stage', tide)
    70647353        Br = Reflective_boundary(domain_drchlt)
    7065         Bd=Dirichlet_boundary([2.0,220,-220])
     7354        Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide])
    70667355        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
    70677356        temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float)
    7068         i=0
    7069         for t in domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,
    7070                                       skip_initial_step = False):
     7357       
     7358        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,
     7359                                                   skip_initial_step = False)):
    70717360            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
    7072             i+=1
    7073 
     7361
     7362       
     7363        #print domain_fbound.quantities['stage'].vertex_values
     7364        #print domain_drchlt.quantities['stage'].vertex_values
     7365                   
    70747366        assert allclose(temp_fbound,temp_drchlt)
     7367
     7368       
     7369        assert allclose(domain_fbound.quantities['stage'].vertex_values,
     7370                        domain_drchlt.quantities['stage'].vertex_values)
     7371                       
     7372        assert allclose(domain_fbound.quantities['xmomentum'].vertex_values,
     7373                        domain_drchlt.quantities['xmomentum'].vertex_values)                       
     7374                       
     7375        assert allclose(domain_fbound.quantities['ymomentum'].vertex_values,
     7376                        domain_drchlt.quantities['ymomentum'].vertex_values)
     7377       
     7378        # Use known Dirichlet condition (if sufficient timesteps have been taken)
     7379        assert allclose(domain_drchlt.quantities['stage'].vertex_values[6], 2)       
     7380        assert allclose(domain_fbound.quantities['stage'].vertex_values[6], 2)
     7381       
     7382       
    70757383
    70767384        try:
     
    70817389       
    70827390        os.remove(meshname)
     7391       
     7392
     7393       
     7394    def Xtest_file_boundary_stsIV_sinewave_ordering(self):
     7395        """test_file_boundary_stsIV_sinewave_ordering(self):
     7396        Read correct points from ordering file and apply sts to boundary
     7397        This one uses a sine wave and compares to time boundary
     7398        """
     7399        #FIXME (Ole): Under construction
     7400       
     7401        from anuga.shallow_water import Domain
     7402        from anuga.shallow_water import Reflective_boundary
     7403        from anuga.shallow_water import Dirichlet_boundary
     7404        from anuga.shallow_water import File_boundary
     7405        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     7406
     7407        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
     7408        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
     7409        tide = 0.35
     7410        time_step_count = 50
     7411        time_step = 2
     7412       
     7413        n=len(lat_long_points)
     7414        first_tstep=ones(n,Int)
     7415        last_tstep=(time_step_count)*ones(n,Int)
     7416       
     7417        gauge_depth=20*ones(n,Float)
     7418        ha=2*ones((n,time_step_count),Float)
     7419        ua=10*ones((n,time_step_count),Float)
     7420        va=-10*ones((n,time_step_count),Float)
     7421       
     7422       
     7423        base_name, files = self.write_mux2(lat_long_points,
     7424                                           time_step_count, time_step,
     7425                                           first_tstep, last_tstep,
     7426                                           depth=gauge_depth,
     7427                                           ha=ha,
     7428                                           ua=ua,
     7429                                           va=va)
     7430
     7431        # Write order file
     7432        file_handle, order_base_name = tempfile.mkstemp("")
     7433        os.close(file_handle)
     7434        os.remove(order_base_name)
     7435        d=","
     7436        order_file=order_base_name+'order.txt'
     7437        fid=open(order_file,'w')
     7438       
     7439        # Write Header
     7440        header='index, longitude, latitude\n'
     7441        fid.write(header)
     7442        indices=[3,0,1]
     7443        for i in indices:
     7444            line=str(i)+d+str(lat_long_points[i][1])+d+\
     7445                str(lat_long_points[i][0])+"\n"
     7446            fid.write(line)
     7447        fid.close()
     7448
     7449        sts_file=base_name
     7450        urs2sts(base_name, basename_out=sts_file,
     7451                ordering_filename=order_file,
     7452                mean_stage=tide,
     7453                verbose=False)
     7454        self.delete_mux(files)
     7455       
     7456       
     7457       
     7458        # Now read the sts file and check that values have been stored correctly.
     7459        fid = NetCDFFile(sts_file + '.sts')
     7460
     7461        # Check the time vector
     7462        times = fid.variables['time'][:]
     7463       
     7464        #print times
     7465
     7466        # Check sts quantities
     7467        stage = fid.variables['stage'][:]
     7468        xmomentum = fid.variables['xmomentum'][:]
     7469        ymomentum = fid.variables['ymomentum'][:]
     7470        elevation = fid.variables['elevation'][:]
     7471
     7472        #print stage
     7473        #print xmomentum
     7474        #print ymomentum
     7475        #print elevation
     7476       
     7477       
     7478
     7479        # Create beginnings of boundary polygon based on sts_boundary
     7480        boundary_polygon = create_sts_boundary(base_name)
     7481       
     7482        os.remove(order_file)
     7483
     7484        # Append the remaining part of the boundary polygon to be defined by
     7485        # the user
     7486        bounding_polygon_utm=[]
     7487        for point in bounding_polygon:
     7488            zone,easting,northing=redfearn(point[0],point[1])
     7489            bounding_polygon_utm.append([easting,northing])
     7490
     7491        boundary_polygon.append(bounding_polygon_utm[3])
     7492        boundary_polygon.append(bounding_polygon_utm[4])
     7493
     7494        #print 'boundary_polygon', boundary_polygon
     7495       
     7496        plot=False
     7497        if plot:
     7498            from pylab import plot,show,axis
     7499            boundary_polygon=ensure_numeric(boundary_polygon)
     7500            bounding_polygon_utm=ensure_numeric(bounding_polygon_utm)
     7501            #plot(lat_long_points[:,0],lat_long_points[:,1],'o')
     7502            plot(boundary_polygon[:,0], boundary_polygon[:,1])
     7503            plot(bounding_polygon_utm[:,0],bounding_polygon_utm[:,1])
     7504            show()
     7505
     7506        assert allclose(bounding_polygon_utm,boundary_polygon)
     7507
     7508
     7509        extent_res=1000000
     7510        meshname = 'urs_test_mesh' + '.tsh'
     7511        interior_regions=None
     7512        boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]}
     7513       
     7514        # have to change boundary tags from last example because now bounding
     7515        # polygon starts in different place.
     7516        create_mesh_from_regions(boundary_polygon,
     7517                                 boundary_tags=boundary_tags,
     7518                                 maximum_triangle_area=extent_res,
     7519                                 filename=meshname,
     7520                                 interior_regions=interior_regions,
     7521                                 verbose=False)
     7522       
     7523        domain_fbound = Domain(meshname)
     7524        domain_fbound.set_quantity('stage', tide)
     7525        Bf = File_boundary(sts_file+'.sts',
     7526                           domain_fbound,
     7527                           boundary_polygon=boundary_polygon)
     7528        Br = Reflective_boundary(domain_fbound)
     7529
     7530        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
     7531        finaltime=time_step*(time_step_count-1)
     7532        yieldstep=time_step
     7533        temp_fbound=zeros(int(finaltime/yieldstep)+1,Float)
     7534   
     7535        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
     7536                                                   finaltime=finaltime,
     7537                                                   skip_initial_step=False)):
     7538            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
     7539   
     7540       
     7541        domain_drchlt = Domain(meshname)
     7542        domain_drchlt.set_quantity('stage', tide)
     7543        Br = Reflective_boundary(domain_drchlt)
     7544        w = 2.0+tide
     7545        h = 20+w
     7546        Bd = Dirichlet_boundary([w, 10*h,-10*h])
     7547        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
     7548       
     7549        temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float)
     7550       
     7551        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,
     7552                                                   skip_initial_step=False)):
     7553            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
     7554
     7555
     7556
     7557        print temp_fbound
     7558        print temp_drchlt
     7559
     7560        print domain_fbound.quantities['stage'].vertex_values
     7561        print domain_drchlt.quantities['stage'].vertex_values
     7562       
     7563        assert allclose(temp_fbound, temp_drchlt)               
     7564        assert allclose(domain_fbound.quantities['stage'].vertex_values,
     7565                        domain_drchlt.quantities['stage'].vertex_values)
     7566                       
     7567        assert allclose(domain_fbound.quantities['xmomentum'].vertex_values,
     7568                        domain_drchlt.quantities['xmomentum'].vertex_values)                       
     7569                       
     7570        assert allclose(domain_fbound.quantities['ymomentum'].vertex_values,
     7571                        domain_drchlt.quantities['ymomentum'].vertex_values)                                               
     7572       
     7573
     7574        try:
     7575            os.remove(sts_file+'.sts')
     7576        except:
     7577            # Windoze can't remove this file for some reason
     7578            pass
     7579       
     7580        os.remove(meshname)
     7581       
     7582       
     7583       
    70837584
    70847585    def test_lon_lat2grid(self):
     
    89449445if __name__ == "__main__":
    89459446
    8946     #suite = unittest.makeSuite(Test_Data_Manager,'test')
    8947     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_read_mux2_pyI')
    89489447    suite = unittest.makeSuite(Test_Data_Manager,'test')
     9448    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_sts')
     9449    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_ordering_different_sources')
    89499450    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
    89509451    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
Note: See TracChangeset for help on using the changeset viewer.