Changeset 5657


Ignore:
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:
7 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/documentation/user_manual/anuga_user_manual.tex

    r5650 r5657  
    22532253The boundary values are obtained from a file and interpolated to the
    22542254appropriate segments for each conserved quantity.
     2255
     2256Optional 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}.
     2257The default_boundary could be a simple Dirichlet condition or even another File\_boundary
     2258typically using data pertaining to another time interval. 
    22552259\end{classdesc}
    22562260
  • anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r5373 r5657  
    22"""
    33
     4from warnings import warn
     5
    46from anuga.utilities.numerical_tools import NAN   
     7from anuga.fit_interpolate.interpolate import Modeltime_too_late
     8from anuga.fit_interpolate.interpolate import Modeltime_too_early
    59
    610
     
    183187    Note that the resulting solution history is not exactly the same as if
    184188    the models were coupled as there is no feedback into the source model.
     189   
     190    Optional keyword argument default_boundary must be either None or
     191    an instance of class descending from class Boundary.
     192    This will be used in case model time exceeds that available in the
     193    underlying data.
    185194       
    186195    """
     
    189198    # rather than File_boundary
    190199
    191     def __init__(self, filename, domain, time_thinning=1,
    192                  use_cache=False, verbose=False, boundary_polygon=None):
     200    def __init__(self, filename, domain,
     201                 time_thinning=1,
     202                 boundary_polygon=None,   
     203                 default_boundary=None,
     204                 use_cache=False,
     205                 verbose=False):
     206
    193207        import time
    194208        from Numeric import array, zeros, Float
     
    251265                               verbose=verbose,
    252266                               boundary_polygon=boundary_polygon)
    253 
     267                             
     268        # Check and store default_boundary
     269        msg = 'Keyword argument default_boundary must be either None '
     270        msg += 'or a boundary object.\n I got %s' %(str(default_boundary))
     271        assert default_boundary is None or isinstance(default_boundary, Boundary), msg
     272        self.default_boundary = default_boundary
     273        self.default_boundary_invoked = False    # Flag
     274
     275        # Store pointer to domain
    254276        self.domain = domain
     277       
     278        self.verbose = verbose
    255279
    256280        # Test
     
    293317
    294318        t = self.domain.time
     319       
    295320        if vol_id is not None and edge_id is not None:
    296321            i = self.boundary_indices[ vol_id, edge_id ]
    297             res = self.F(t, point_id = i)
     322           
     323           
     324            #res = self.F(t, point_id = i)           
     325            try:
     326                res = self.F(t, point_id = i)
     327            except (Modeltime_too_late, Modeltime_too_early), e:
     328                if self.default_boundary is None:
     329                    raise Exception, e # Reraise exception
     330                else:
     331                    if self.default_boundary_invoked is False:
     332                        # Issue warning the first time
     333                        msg = '%s' %str(e)
     334                        msg += 'Instead I will using the default boundary: %s\n'\
     335                            %str(self.default_boundary)
     336                        msg += 'Note: Further warnings will be supressed'
     337                        warn(msg)
     338                   
     339                        # FIXME (Ole): Replace this crude flag with
     340                        # Python's ability to print warnings only once.
     341                        # See http://docs.python.org/lib/warning-filter.html
     342                        self.default_boundary_invoked = True
     343                   
     344                    # Pass control to default boundary
     345                    res = self.default_boundary.evaluate(vol_id, edge_id)
     346                   
     347               
    298348            if res == NAN:
    299349                x,y=self.midpoint_coordinates[i,:]
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

    r5221 r5657  
    104104        """Test that boundary values can be read from file and interpolated
    105105        This is using the .tms file format
     106       
     107        See also test_util for comprenhensive testing of the underlying
     108        file_function and also tests in test_datamanager which tests
     109        file_function using the sts format
    106110        """
    107111
     
    269273#-------------------------------------------------------------
    270274if __name__ == "__main__":
    271     suite = unittest.makeSuite(Test_Generic_Boundary_Conditions,'test')
     275    suite = unittest.makeSuite(Test_Generic_Boundary_Conditions, 'test')
    272276    runner = unittest.TextTestRunner()
    273277    runner.run(suite)
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

    r5442 r5657  
    521521
    522522
    523     def qtest_spatio_temporal_file_function_time(self):
     523    def test_spatio_temporal_file_function_time(self):
    524524        """Test that File function interpolates correctly
    525525        between given times.
     
    554554        points, vertices, boundary =\
    555555                rectangular(4, 4, 15, 30, origin = (0, -20))
    556         print "points", points
     556        #print "points", points
    557557
    558558        #print 'Number of elements', len(vertices)
     
    639639                #print ' ', q0
    640640                #print ' ', q1
    641                 print "q",q
    642                 print "actual", actual
     641                #print "q",q
     642                #print "actual", actual
    643643                #print
    644644                if q0 == NAN:
     
    696696
    697697
    698     def NOtest_spatio_temporal_file_function_time(self):
     698    def Xtest_spatio_temporal_file_function_time(self):
    699699        # FIXME: This passes but needs some TLC
    700700        # Test that File function interpolates correctly
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r5566 r5657  
    3939from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
    4040from anuga.abstract_2d_finite_volumes.util import file_function
     41
     42# Interpolation specific exceptions
     43
     44class Modeltime_too_late(Exception): pass
     45class Modeltime_too_early(Exception): pass
     46
    4147
    4248
     
    856862        msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1])
    857863        msg += ' does not match model time: %.16f\n' %t
    858         if t < self.time[0]: raise Exception(msg)
    859         if t > self.time[-1]: raise Exception(msg)
     864        if t < self.time[0]: raise Modeltime_too_early(msg)
     865        if t > self.time[-1]: raise Modeltime_too_late(msg)
    860866
    861867        oldindex = self.index #Time index
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r5644 r5657  
    10881088    def __init__(self, filename, domain,
    10891089                 mean_stage=0.0,
    1090                  time_thinning=1,
     1090                 time_thinning=1,
     1091                 boundary_polygon=None,   
     1092                 default_boundary=None,                 
    10911093                 use_cache=False,
    1092                  verbose=False,
    1093                  boundary_polygon=None):
     1094                 verbose=False):
    10941095        """Constructor
    10951096
     
    11071108                       the speed of a model run that you are setting up
    11081109                       and testing.
     1110                       
     1111        default_boundary: Must be either None or an instance of a
     1112                          class descending from class Boundary.
     1113                          This will be used in case model time exceeds
     1114                          that available in the underlying data.
     1115                                               
    11091116        use_cache:
    11101117        verbose:
     
    11151122        self.file_boundary = File_boundary(filename, domain,
    11161123                                           time_thinning=time_thinning,
     1124                                           boundary_polygon=boundary_polygon,
     1125                                           default_boundary=default_boundary,
    11171126                                           use_cache=use_cache,
    1118                                            verbose=verbose,
    1119                                            boundary_polygon=boundary_polygon)
     1127                                           verbose=verbose)
     1128
    11201129       
    11211130        # Record information from File_boundary
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5627 r5657  
    76437643        """
    76447644       
     7645        # FIXME (Ole): These tests should really move to
     7646        # test_generic_boundaries.py
     7647       
    76457648        from anuga.shallow_water import Domain
    76467649        from anuga.shallow_water import Reflective_boundary
     
    77747777        os.remove(sts_file+'.sts')
    77757778        os.remove(meshname)
     7779       
     7780       
     7781       
     7782               
     7783       
     7784    def test_file_boundary_stsI_beyond_model_time(self):
     7785        """test_file_boundary_stsI(self):
     7786       
     7787        Test that file_boundary can work when model time
     7788        exceeds available data.
     7789        This is optional and requires the user to specify a default
     7790        boundary object.
     7791        """
     7792       
     7793        # Don't do warnings in unit test
     7794        import warnings
     7795        warnings.simplefilter('ignore')
     7796       
     7797        from anuga.shallow_water import Domain
     7798        from anuga.shallow_water import Reflective_boundary
     7799        from anuga.shallow_water import Dirichlet_boundary
     7800        from anuga.shallow_water import File_boundary
     7801        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     7802
     7803        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
     7804        tide = 0.37
     7805        time_step_count = 5
     7806        time_step = 2
     7807        lat_long_points =bounding_polygon[0:3]
     7808        n=len(lat_long_points)
     7809        first_tstep=ones(n,Int)
     7810        last_tstep=(time_step_count)*ones(n,Int)
     7811
     7812        h = 20       
     7813        w = 2
     7814        u = 10
     7815        v = -10
     7816        gauge_depth=h*ones(n,Float)
     7817        ha=w*ones((n,time_step_count),Float)
     7818        ua=u*ones((n,time_step_count),Float)
     7819        va=v*ones((n,time_step_count),Float)
     7820        base_name, files = self.write_mux2(lat_long_points,
     7821                                           time_step_count, time_step,
     7822                                           first_tstep, last_tstep,
     7823                                           depth=gauge_depth,
     7824                                           ha=ha,
     7825                                           ua=ua,
     7826                                           va=va)
     7827
     7828        sts_file=base_name
     7829        urs2sts(base_name,
     7830                sts_file,
     7831                mean_stage=tide,
     7832                verbose=False)
     7833        self.delete_mux(files)
     7834
     7835        #print 'start create mesh from regions'
     7836        for i in range(len(bounding_polygon)):
     7837            zone,\
     7838            bounding_polygon[i][0],\
     7839            bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],
     7840                                            bounding_polygon[i][1])
     7841                                           
     7842        extent_res=1000000
     7843        meshname = 'urs_test_mesh' + '.tsh'
     7844        interior_regions=None
     7845        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
     7846        create_mesh_from_regions(bounding_polygon,
     7847                                 boundary_tags=boundary_tags,
     7848                                 maximum_triangle_area=extent_res,
     7849                                 filename=meshname,
     7850                                 interior_regions=interior_regions,
     7851                                 verbose=False)
     7852       
     7853        domain_fbound = Domain(meshname)
     7854        domain_fbound.set_quantity('stage', tide)
     7855       
     7856        Br = Reflective_boundary(domain_fbound)
     7857        Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)])       
     7858        Bf = File_boundary(sts_file+'.sts',
     7859                           domain_fbound,
     7860                           boundary_polygon=bounding_polygon,
     7861                           default_boundary=Bd) # Condition to be used
     7862                                                # if model time exceeds
     7863                                                # available data
     7864
     7865        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
     7866       
     7867        # Check boundary object evaluates as it should
     7868        for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects):
     7869            if B is Bf:
     7870           
     7871                qf = B.evaluate(vol_id, edge_id)  # File boundary
     7872                qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary
     7873
     7874                assert allclose(qf, qd)
     7875               
     7876       
     7877        # Evolve
     7878        data_finaltime = time_step*(time_step_count-1)
     7879        finaltime = data_finaltime + 10 # Let model time exceed available data
     7880        yieldstep = time_step
     7881        temp_fbound=zeros(int(finaltime/yieldstep)+1, Float)
     7882
     7883        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
     7884                                                   finaltime=finaltime,
     7885                                                   skip_initial_step=False)):
     7886                                                   
     7887            D = domain_fbound
     7888            temp_fbound[i]=D.quantities['stage'].centroid_values[2]
     7889
     7890            # Check that file boundary object has populated
     7891            # boundary array correctly 
     7892            # FIXME (Ole): Do this for the other tests too!
     7893            for j, val in enumerate(D.get_quantity('stage').boundary_values):
     7894           
     7895                (vol_id, edge_id), B = D.boundary_objects[j]
     7896                if isinstance(B, File_boundary):
     7897                    #print j, val
     7898                    assert allclose(val, w + tide)
     7899
     7900
     7901
     7902       
     7903       
     7904       
     7905           
     7906
     7907       
     7908        domain_drchlt = Domain(meshname)
     7909        domain_drchlt.set_quantity('stage', tide)
     7910        Br = Reflective_boundary(domain_drchlt)
     7911
     7912        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
     7913        temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float)
     7914
     7915        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
     7916                                                   finaltime=finaltime,
     7917                                                   skip_initial_step=False)):
     7918            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
     7919
     7920        #print domain_fbound.quantities['stage'].vertex_values
     7921        #print domain_drchlt.quantities['stage'].vertex_values
     7922                   
     7923        assert allclose(temp_fbound,temp_drchlt)
     7924       
     7925        assert allclose(domain_fbound.quantities['stage'].vertex_values,
     7926                        domain_drchlt.quantities['stage'].vertex_values)
     7927                       
     7928        assert allclose(domain_fbound.quantities['xmomentum'].vertex_values,
     7929                        domain_drchlt.quantities['xmomentum'].vertex_values)                       
     7930                       
     7931        assert allclose(domain_fbound.quantities['ymomentum'].vertex_values,
     7932                        domain_drchlt.quantities['ymomentum'].vertex_values)                                               
     7933       
     7934       
     7935        os.remove(sts_file+'.sts')
     7936        os.remove(meshname)
     7937
     7938       
    77767939
    77777940    def test_file_boundary_stsII(self):
     
    1010510268
    1010610269    suite = unittest.makeSuite(Test_Data_Manager,'test')
    10107     #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsI')
     10270    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsI_beyond_model_time')
    1010810271    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV_sinewave_ordering')
    1010910272    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
Note: See TracChangeset for help on using the changeset viewer.