Changeset 6122


Ignore:
Timestamp:
Jan 8, 2009, 12:18:48 PM (10 years ago)
Author:
ole
Message:

Added tests and a debug example for outflow conservation issue.

Files:
1 added
1 edited

Legend:

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

    r6086 r6122  
    1313from anuga.abstract_2d_finite_volumes.quantity import Quantity
    1414from anuga.geospatial_data.geospatial_data import Geospatial_data
    15 
     15from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    1616from shallow_water_domain import *
    1717
     
    30073007
    30083008
    3009        
    3010        
     3009
     3010    def Xtest_inflow_outflow_conservation(self):
     3011        """test_inflow_outflow_conservation
     3012       
     3013        Test what happens if water is abstracted from one area and
     3014        injected into another - especially if there is not enough
     3015        water to match the abstraction.
     3016        This tests that the total volume is kept constant under a range of
     3017        scenarios.
     3018        """
     3019       
     3020        from math import pi, cos, sin
     3021       
     3022        length = 20.
     3023        width = 10.
     3024
     3025        dx = dy = 2  # 1 or 2 OK
     3026        points, vertices, boundary = rectangular_cross(int(length/dx),
     3027                                                       int(width/dy),
     3028                                                       len1=length,
     3029                                                       len2=width)
     3030        domain = Domain(points, vertices, boundary)   
     3031        domain.set_name('test_inflow_conservation')  # Output name
     3032        domain.set_default_order(2)
     3033       
     3034
     3035        # Flat surface with 1m of water
     3036        stage = 1.0
     3037        domain.set_quantity('elevation', 0)
     3038        domain.set_quantity('stage', stage)
     3039        domain.set_quantity('friction', 0)
     3040
     3041        Br = Reflective_boundary(domain)
     3042        domain.set_boundary({'left': Br, 'right': Br, 'bottom': Br, 'top': Br})
     3043
     3044        # Setup one forcing term, constant inflow of 2 m^3/s on a circle
     3045        domain.forcing_terms = []
     3046        domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))
     3047
     3048        domain.compute_forcing_terms()
     3049        #print domain.quantities['stage'].explicit_update
     3050       
     3051        # Check that update values are correct
     3052        for x in domain.quantities['stage'].explicit_update:
     3053            assert allclose(x, 2.0/pi) or allclose(x, 0.0)
     3054
     3055       
     3056        # Check volumes without inflow
     3057        domain.forcing_terms = []       
     3058        initial_volume = domain.quantities['stage'].get_integral()
     3059       
     3060        assert allclose(initial_volume, width*length*stage)
     3061       
     3062        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
     3063            volume =  domain.quantities['stage'].get_integral()
     3064            assert allclose (volume, initial_volume)
     3065           
     3066           
     3067        # Now apply the inflow and check volumes for a range of stage values
     3068        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
     3069            domain.time = 0.0
     3070            domain.set_quantity('stage', stage)       
     3071                   
     3072            domain.forcing_terms = []       
     3073            domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))       
     3074            initial_volume = domain.quantities['stage'].get_integral()
     3075            predicted_volume = initial_volume
     3076            dt = 0.05
     3077            for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
     3078                volume = domain.quantities['stage'].get_integral()
     3079               
     3080                assert allclose (volume, predicted_volume)           
     3081                predicted_volume = predicted_volume + 2.0/pi/100/dt # Why 100?
     3082           
     3083           
     3084        # Apply equivalent outflow only and check volumes for a range of stage values
     3085        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
     3086            print stage
     3087           
     3088            domain.time = 0.0
     3089            domain.set_quantity('stage', stage)       
     3090            domain.forcing_terms = []       
     3091            domain.forcing_terms.append(Inflow(domain, rate=-2.0, center=(15,5), radius=1))       
     3092            initial_volume = domain.quantities['stage'].get_integral()
     3093            predicted_volume = initial_volume
     3094            dt = 0.05
     3095            for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
     3096                volume = domain.quantities['stage'].get_integral()
     3097               
     3098                print t, volume, predicted_volume
     3099                assert allclose (volume, predicted_volume)           
     3100                predicted_volume = predicted_volume - 2.0/pi/100/dt # Why 100?           
     3101           
     3102           
     3103        # Apply both inflow and outflow and check volumes being constant for a
     3104        # range of stage values
     3105        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:       
     3106            print stage
     3107           
     3108            domain.time = 0.0
     3109            domain.set_quantity('stage', stage)       
     3110            domain.forcing_terms = []       
     3111            domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))       
     3112            domain.forcing_terms.append(Inflow(domain, rate=-2.0, center=(15,5), radius=1))               
     3113            initial_volume = domain.quantities['stage'].get_integral()
     3114
     3115            dt = 0.05
     3116            for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
     3117                volume = domain.quantities['stage'].get_integral()
     3118               
     3119                print t, volume
     3120                assert allclose (volume, initial_volume)           
     3121
     3122           
     3123           
    30113124
    30123125    #####################################################
     
    64486561if __name__ == "__main__":
    64496562
    6450     suite = unittest.makeSuite(Test_Shallow_Water, 'test')
     6563    #suite = unittest.makeSuite(Test_Shallow_Water, 'test')
    64516564    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
    64526565    #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_energy_through_cross_section_with_g')   
    64536566    #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')   
    64546567    #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
    6455     #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')
    6456     #suite = unittest.makeSuite(Test_Shallow_Water,'test_time_dependent_rainfall')   
     6568    suite = unittest.makeSuite(Test_Shallow_Water,'test_inflow_outflow_conservation')
     6569    #suite = unittest.makeSuite(Test_Shallow_Water,'test_outflow_conservation_problem_temp')   
    64576570   
    64586571
Note: See TracChangeset for help on using the changeset viewer.