source: trunk/anuga_core/source/anuga/culvert_flows/test_new_culvert_class.py @ 7858

#!/usr/bin/env python


import unittest
import os.path
import sys

import anuga

from anuga.utilities.system_tools import get_pathname_from_package
from anuga.geometry.polygon_function import Polygon_function
        
from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
from anuga.abstract_2d_finite_volumes.quantity import Quantity

from anuga.culvert_flows.new_culvert_class import Culvert_flow
from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
     
from math import pi,pow,sqrt

import numpy as num


# Helper functions
def run_culvert_flow_problem(depth):
    """Run flow with culvert given depth
    """

    length = 40.
    width = 5.

    dx = dy = 1           # Resolution: Length of subdivisions on both axes

    points, vertices, boundary = rectangular_cross(int(length/dx),
                                                   int(width/dy),
                                                   len1=length, 
                                                   len2=width)
    domain = anuga.Domain(points, vertices, boundary)   
    domain.set_name('Test_culvert_shallow') # Output name
    domain.set_default_order(2)


    #----------------------------------------------------------------------
    # Setup initial conditions
    #----------------------------------------------------------------------

    def topography(x, y):
        """Set up a weir
        
        A culvert will connect either side
        """
        # General Slope of Topography
        z=-x/1000
        
        N = len(x)
        for i in range(N):

           # Sloping Embankment Across Channel
            if 5.0 < x[i] < 10.1:
                # Cut Out Segment for Culvert face                
                if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
                   z[i]=z[i]
                else:
                   z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
            if 10.0 < x[i] < 12.1:
               z[i] +=  2.5                    # Flat Crest of Embankment
            if 12.0 < x[i] < 14.5:
                # Cut Out Segment for Culvert face                
                if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
                   z[i]=z[i]
                else:
                   z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
                           
                
        return z


    domain.set_quantity('elevation', topography) 
    domain.set_quantity('friction', 0.01)         # Constant friction 
    domain.set_quantity('stage',
                        expression='elevation + %f' % depth) # Shallow initial condition

    # Boyd culvert
    culvert = Culvert_flow(domain,
                           label='Culvert No. 1',
                           description='This culvert is a test unit 1.2m Wide by 0.75m High',   
                           end_point0=[9.0, 2.5], 
                           end_point1=[13.0, 2.5],
                           width=1.20, height=0.75,
                           culvert_routine=boyd_generalised_culvert_model,
                           number_of_barrels=1,
                           update_interval=2,
                           verbose=False)
    

    domain.forcing_terms.append(culvert)
    

    #-----------------------------------------------------------------------
    # Setup boundary conditions
    #-----------------------------------------------------------------------

    # Inflow based on Flow Depth and Approaching Momentum

    Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
    domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})


    
    #-----------------------------------------------------------------------
    # Evolve system through time
    #-----------------------------------------------------------------------

    #print 'depth', depth
    ref_volume = domain.get_quantity('stage').get_integral()
    for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
        new_volume = domain.get_quantity('stage').get_integral()
        
        msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
               % (new_volume, ref_volume))
        assert num.allclose(new_volume, ref_volume), msg        
        

    os.remove('Test_culvert_shallow.sww')

class Test_Culvert(unittest.TestCase):
    def setUp(self):
        pass

    def tearDown(self):
        pass


    def test_that_culvert_runs_rating(self):
        """test_that_culvert_runs_rating
        
        This test exercises the culvert and checks values outside rating curve
        are dealt with       
        """

        path = get_pathname_from_package('anuga.culvert_flows')    
        
        length = 40.
        width = 5.

        dx = dy = 1           # Resolution: Length of subdivisions on both axes

        points, vertices, boundary = rectangular_cross(int(length/dx),
                                                       int(width/dy),
                                                       len1=length, 
                                                       len2=width)
        domain = anuga.Domain(points, vertices, boundary)   
        domain.set_name('Test_culvert')                 # Output name
        domain.set_default_order(2)


        #----------------------------------------------------------------------
        # Setup initial conditions
        #----------------------------------------------------------------------

        def topography(x, y):
            """Set up a weir
            
            A culvert will connect either side
            """
            # General Slope of Topography
            z=-x/1000
            
            N = len(x)
            for i in range(N):

               # Sloping Embankment Across Channel
                if 5.0 < x[i] < 10.1:
                    # Cut Out Segment for Culvert face                
                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
                       z[i]=z[i]
                    else:
                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
                if 10.0 < x[i] < 12.1:
                   z[i] +=  2.5                    # Flat Crest of Embankment
                if 12.0 < x[i] < 14.5:
                    # Cut Out Segment for Culvert face                
                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
                       z[i]=z[i]
                    else:
                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
                                   
                        
            return z


        domain.set_quantity('elevation', topography) 
        domain.set_quantity('friction', 0.01)         # Constant friction 
        domain.set_quantity('stage',
                            expression='elevation')   # Dry initial condition

        filename=os.path.join(path, 'example_rating_curve.csv')
        culvert = Culvert_flow(domain,
                               culvert_description_filename=filename,        
                               end_point0=[9.0, 2.5], 
                               end_point1=[13.0, 2.5],
                               width=1.00,
                               use_velocity_head=True,
                               verbose=False)

        domain.forcing_terms.append(culvert)
        

        #-----------------------------------------------------------------------
        # Setup boundary conditions
        #-----------------------------------------------------------------------

        # Inflow based on Flow Depth and Approaching Momentum
        Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
        Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
        
        # Upstream and downstream conditions that will exceed the rating curve
        # I.e produce delta_h outside the range [0, 10] specified in the the 
        # file example_rating_curve.csv
        Btus = anuga.Time_boundary(domain, lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
        Btds = anuga.Time_boundary(domain, lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
        domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})


        #-----------------------------------------------------------------------
        # Evolve system through time
        #-----------------------------------------------------------------------

        min_delta_w = sys.maxint 
        max_delta_w = -min_delta_w
        for t in domain.evolve(yieldstep = 1, finaltime = 25):
            delta_w = culvert.inlet.stage - culvert.outlet.stage
            
            if delta_w > max_delta_w: max_delta_w = delta_w
            if delta_w < min_delta_w: min_delta_w = delta_w            
            
            pass

        # Check that extreme values in rating curve have been exceeded
        # so that we know that condition has been exercised
        assert min_delta_w < 0
        assert max_delta_w > 10        


        os.remove('Test_culvert.sww')

    def test_that_culvert_dry_bed_rating_does_not_produce_flow(self):
        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
        
        Test that culvert on a sloping dry bed doesn't produce flows
        although there will be a 'pressure' head due to delta_w > 0

        This one is using the rating curve variant
        """

        path = get_pathname_from_package('anuga.culvert_flows')    
        
        length = 40.
        width = 5.

        dx = dy = 1           # Resolution: Length of subdivisions on both axes

        points, vertices, boundary = rectangular_cross(int(length/dx),
                                                       int(width/dy),
                                                       len1=length, 
                                                       len2=width)
        domain = anuga.Domain(points, vertices, boundary)   
        domain.set_name('Test_culvert_dry') # Output name
        domain.set_default_order(2)


        #----------------------------------------------------------------------
        # Setup initial conditions
        #----------------------------------------------------------------------

        def topography(x, y):
            """Set up a weir
            
            A culvert will connect either side
            """
            # General Slope of Topography
            z=-x/1000
            
            N = len(x)
            for i in range(N):

               # Sloping Embankment Across Channel
                if 5.0 < x[i] < 10.1:
                    # Cut Out Segment for Culvert face                
                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
                       z[i]=z[i]
                    else:
                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
                if 10.0 < x[i] < 12.1:
                   z[i] +=  2.5                    # Flat Crest of Embankment
                if 12.0 < x[i] < 14.5:
                    # Cut Out Segment for Culvert face                
                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
                       z[i]=z[i]
                    else:
                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
                                   
                        
            return z


        domain.set_quantity('elevation', topography) 
        domain.set_quantity('friction', 0.01)         # Constant friction 
        domain.set_quantity('stage',
                            expression='elevation')   # Dry initial condition


        filename = os.path.join(path, 'example_rating_curve.csv')
        culvert = Culvert_flow(domain,
                               culvert_description_filename=filename,
                               end_point0=[9.0, 2.5], 
                               end_point1=[13.0, 2.5],
                               height=0.75,
                               verbose=False)

        domain.forcing_terms.append(culvert)
        

        #-----------------------------------------------------------------------
        # Setup boundary conditions
        #-----------------------------------------------------------------------

        # Inflow based on Flow Depth and Approaching Momentum

        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})


        #-----------------------------------------------------------------------
        # Evolve system through time
        #-----------------------------------------------------------------------

        ref_volume = domain.get_quantity('stage').get_integral()
        for t in domain.evolve(yieldstep = 1, finaltime = 25):
            new_volume = domain.get_quantity('stage').get_integral()
            
            msg = 'Total volume has changed'
            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
            pass
    

    
        os.remove('Test_culvert_dry.sww')
        
    def test_that_culvert_flows_conserves_volume(self):
        """test_that_culvert_flows_conserves_volume

        Test that culvert on a sloping dry bed limits flows when very little water
        is present at inlet.

        Uses helper function: run_culvert_flow_problem(depth):

        """

        # Try this for a range of depths
        for depth in [0.1, 0.2, 0.5, 1.0]:
            run_culvert_flow_problem(depth)


    def OBSOLETE_XXXtest_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self):
        """test_that_culvert_rating_limits_flow_in_shallow_inlet_condition
        
        Test that culvert on a sloping dry bed limits flows when very little water
        is present at inlet

        This one is using the rating curve variant
        """

        

        path = get_pathname_from_package('anuga.culvert_flows')    
        
        length = 40.
        width = 5.

        dx = dy = 1           # Resolution: Length of subdivisions on both axes

        points, vertices, boundary = rectangular_cross(int(length/dx),
                                                       int(width/dy),
                                                       len1=length, 
                                                       len2=width)
        domain = anuga.Domain(points, vertices, boundary)   
        domain.set_name('Test_culvert_shallow') # Output name
        domain.set_default_order(2)


        #----------------------------------------------------------------------
        # Setup initial conditions
        #----------------------------------------------------------------------

        def topography(x, y):
            """Set up a weir
            
            A culvert will connect either side
            """
            # General Slope of Topography
            z=-x/1000
            
            N = len(x)
            for i in range(N):

               # Sloping Embankment Across Channel
                if 5.0 < x[i] < 10.1:
                    # Cut Out Segment for Culvert face                
                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
                       z[i]=z[i]
                    else:
                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
                if 10.0 < x[i] < 12.1:
                   z[i] +=  2.5                    # Flat Crest of Embankment
                if 12.0 < x[i] < 14.5:
                    # Cut Out Segment for Culvert face                
                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
                       z[i]=z[i]
                    else:
                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
                                   
                        
            return z


        domain.set_quantity('elevation', topography) 
        domain.set_quantity('friction', 0.01)         # Constant friction 
        domain.set_quantity('stage',
                            expression='elevation + 0.1') # Shallow initial condition

        # Boyd culvert
        culvert = Culvert_flow(domain,
                               label='Culvert No. 1',
                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
                               end_point0=[9.0, 2.5], 
                               end_point1=[13.0, 2.5],
                               width=1.20, height=0.75,
                               culvert_routine=boyd_generalised_culvert_model,
                               number_of_barrels=1,
                               update_interval=2,
                               verbose=False)
        
        # Rating curve
        #filename = os.path.join(path, 'example_rating_curve.csv')
        #culvert = Culvert_flow(domain,
        #                       culvert_description_filename=filename,
        #                       end_point0=[9.0, 2.5], 
        #                       end_point1=[13.0, 2.5],
        #                       trigger_depth=0.01,
        #                       verbose=False)

        domain.forcing_terms.append(culvert)
        

        #-----------------------------------------------------------------------
        # Setup boundary conditions
        #-----------------------------------------------------------------------

        # Inflow based on Flow Depth and Approaching Momentum

        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})


        
        #-----------------------------------------------------------------------
        # Evolve system through time
        #-----------------------------------------------------------------------

        print 'depth', 0.1
        ref_volume = domain.get_quantity('stage').get_integral()
        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
            new_volume = domain.get_quantity('stage').get_integral()

            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
                   % (new_volume, ref_volume))
            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg        
        
        
        return
        # Now try this again for a depth of 10 cm and for a range of other depths
        for depth in [0.1, 0.2, 0.5, 1.0]:
            print 'depth', depth
            domain.set_time(0.0)
            
            domain.set_quantity('elevation', topography) 
            domain.set_quantity('friction', 0.01)         # Constant friction 
            domain.set_quantity('stage',
                                expression='elevation + %f' % depth)
            
        
            ref_volume = domain.get_quantity('stage').get_integral()
            for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
                new_volume = domain.get_quantity('stage').get_integral()
            
                msg = 'Total volume has changed: Is %.8f m^3 should have been %.8f m^3'\
                    % (new_volume, ref_volume)

                assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
    
    
    
    def test_that_culvert_dry_bed_boyd_does_not_produce_flow(self):
        """test_that_culvert_in_dry_bed_boyd_does_not_produce_flow(self):
        
        Test that culvert on a sloping dry bed doesn't produce flows
        although there will be a 'pressure' head due to delta_w > 0

        This one is using the 'Boyd' variant        
        """

        path = get_pathname_from_package('anuga.culvert_flows')    
        
        length = 40.
        width = 5.

        dx = dy = 1           # Resolution: Length of subdivisions on both axes

        points, vertices, boundary = rectangular_cross(int(length/dx),
                                                       int(width/dy),
                                                       len1=length, 
                                                       len2=width)
        domain = anuga.Domain(points, vertices, boundary)   
        domain.set_name('Test_culvert_dry') # Output name
        domain.set_default_order(2)


        #----------------------------------------------------------------------
        # Setup initial conditions
        #----------------------------------------------------------------------

        def topography(x, y):
            """Set up a weir
            
            A culvert will connect either side
            """
            # General Slope of Topography
            z=-x/1000
            
            N = len(x)
            for i in range(N):

               # Sloping Embankment Across Channel
                if 5.0 < x[i] < 10.1:
                    # Cut Out Segment for Culvert face                
                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
                       z[i]=z[i]
                    else:
                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
                if 10.0 < x[i] < 12.1:
                   z[i] +=  2.5                    # Flat Crest of Embankment
                if 12.0 < x[i] < 14.5:
                    # Cut Out Segment for Culvert face                
                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
                       z[i]=z[i]
                    else:
                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
                                   
                        
            return z


        domain.set_quantity('elevation', topography) 
        domain.set_quantity('friction', 0.01)         # Constant friction 
        domain.set_quantity('stage',
                            expression='elevation')   # Dry initial condition


        filename = os.path.join(path, 'example_rating_curve.csv')


        culvert = Culvert_flow(domain,
                               label='Culvert No. 1',
                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
                               end_point0=[9.0, 2.5], 
                               end_point1=[13.0, 2.5],
                               width=1.20, height=0.75,
                               culvert_routine=boyd_generalised_culvert_model,
                               number_of_barrels=1,
                               update_interval=2,
                               verbose=False)
        
        domain.forcing_terms.append(culvert)
        

        #-----------------------------------------------------------------------
        # Setup boundary conditions
        #-----------------------------------------------------------------------

        # Inflow based on Flow Depth and Approaching Momentum

        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})


        #-----------------------------------------------------------------------
        # Evolve system through time
        #-----------------------------------------------------------------------

        ref_volume = domain.get_quantity('stage').get_integral()
        for t in domain.evolve(yieldstep = 1, finaltime = 25):
            
            new_volume = domain.get_quantity('stage').get_integral()

            msg = 'Total volume has changed'
            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
            pass
    

    
    

    def test_predicted_boyd_flow(self):
        """test_predicted_boyd_flow
        
        Test that flows predicted by the boyd method are consistent with what what
        is calculated in engineering codes.
        The data was supplied by Petar Milevski
        """

        # FIXME(Ole) this is nowhere near finished
        path = get_pathname_from_package('anuga.culvert_flows')    
        
        length = 12.
        width = 5.

        dx = dy = 0.5           # Resolution: Length of subdivisions on both axes

        points, vertices, boundary = rectangular_cross(int(length/dx),
                                                       int(width/dy),
                                                       len1=length, 
                                                       len2=width)
        domain = anuga.Domain(points, vertices, boundary)   
        
        domain.set_name('test_culvert')                 # Output name
        domain.set_default_order(2)


        #----------------------------------------------------------------------
        # Setup initial conditions
        #----------------------------------------------------------------------

        def topography(x, y):
            # General Slope of Topography
            z=-x/10
            
            return z


        domain.set_quantity('elevation', topography) 
        domain.set_quantity('friction', 0.01)         # Constant friction 
        domain.set_quantity('stage', expression='elevation')


        Q0 = domain.get_quantity('stage')
        Q1 = Quantity(domain)
        
        # Add depths to stage        
        head_water_depth = 0.169
        tail_water_depth = 0.089
        
        inlet_poly = [[0,0], [6,0], [6,5], [0,5]]
        outlet_poly = [[6,0], [12,0], [12,5], [6,5]]        
        
        Q1.set_values(Polygon_function([(inlet_poly, head_water_depth),
                                        (outlet_poly, tail_water_depth)]))
        
        domain.set_quantity('stage', Q0 + Q1)



        culvert = Culvert_flow(domain,
                               label='Test culvert',
                               description='4 m test culvert',   
                               end_point0=[4.0, 2.5], 
                               end_point1=[8.0, 2.5],
                               width=1.20, 
                               height=0.75,
                               culvert_routine=boyd_generalised_culvert_model,        
                               number_of_barrels=1,
                               verbose=False)
                               
        
        domain.forcing_terms.append(culvert)
        
        # Call
        culvert(domain)
    

    

    def test_momentum_jet(self):
        """test_momentum_jet
        
        Test that culvert_class can accept keyword use_momentum_jet
        This does not yet imply that the values have been tested. FIXME
        """


        length = 40.
        width = 5.

        dx = dy = 1           # Resolution: Length of subdivisions on both axes

        points, vertices, boundary = rectangular_cross(int(length/dx),
                                                       int(width/dy),
                                                       len1=length, 
                                                       len2=width)
        domain = anuga.Domain(points, vertices, boundary)   
        domain.set_name('Test_culvert_shallow') # Output name
        domain.set_default_order(2)


        #----------------------------------------------------------------------
        # Setup initial conditions
        #----------------------------------------------------------------------

        def topography(x, y):
            """Set up a weir
            
            A culvert will connect either side
            """
            # General Slope of Topography
            z=-x/1000
            
            N = len(x)
            for i in range(N):

               # Sloping Embankment Across Channel
                if 5.0 < x[i] < 10.1:
                    # Cut Out Segment for Culvert face                
                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
                       z[i]=z[i]
                    else:
                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
                if 10.0 < x[i] < 12.1:
                   z[i] +=  2.5                    # Flat Crest of Embankment
                if 12.0 < x[i] < 14.5:
                    # Cut Out Segment for Culvert face                
                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
                       z[i]=z[i]
                    else:
                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
                                   
                        
            return z


        domain.set_quantity('elevation', topography) 
        domain.set_quantity('friction', 0.01)         # Constant friction 
        domain.set_quantity('stage',
                            expression='elevation + 1.0') # Shallow initial condition

        # Boyd culvert
        culvert = Culvert_flow(domain,
                               label='Culvert No. 1',
                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
                               end_point0=[9.0, 2.5], 
                               end_point1=[13.0, 2.5],
                               width=1.20, height=0.75,
                               culvert_routine=boyd_generalised_culvert_model,
                               number_of_barrels=1,
                               use_momentum_jet=True,
                               update_interval=2,
                               verbose=False)
    

        domain.forcing_terms.append(culvert)

        
        # Call
        culvert(domain)
    

        #-----------------------------------------------------------------------
        # Setup boundary conditions
        #-----------------------------------------------------------------------


        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})

        #-----------------------------------------------------------------------
        # Evolve system through time
        #-----------------------------------------------------------------------

        ref_volume = domain.get_quantity('stage').get_integral()
        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
            new_volume = domain.get_quantity('stage').get_integral()
            
            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
                   % (new_volume, ref_volume))
            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg        
        

        

               
#-------------------------------------------------------------

if __name__ == "__main__":
    suite = unittest.makeSuite(Test_Culvert, 'test')
    runner = unittest.TextTestRunner() #verbosity=2)
    runner.run(suite)