Changeset 7090


Ignore:
Timestamp:
May 26, 2009, 5:41:49 PM (11 years ago)
Author:
ole
Message:

Added tests for volumetric conservation in the presence of inflow or rainfall.

Location:
anuga_core/source/anuga/shallow_water
Files:
2 edited

Legend:

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

    r7043 r7090  
    880880        """Create volumetric balance report suitable for printing or logging.
    881881        """
    882        
    883         boundary_flows, total_boundary_inflow, total_boundary_outflow = self.compute_boundary_flows()
     882
     883        X = self.compute_boundary_flows()
     884        boundary_flows, total_boundary_inflow, total_boundary_outflow = X
    884885       
    885886        s = '---------------------------\n'       
    886887        s += 'Volumetric balance report:\n'
    887888        s += '--------------------------\n'
    888         s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
    889         s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow       
     889        s += 'Total boundary inflow [m^3/s]: %.2e\n' % total_boundary_inflow
     890        s += 'Total boundary outflow [m^3/s]: %.2e\n' % total_boundary_outflow       
    890891        s += 'Net boundary flow by tags [m^3/s]\n'
    891892        for tag in boundary_flows:
    892             s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
    893        
    894         s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow)
    895         s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
     893            s += '    %s [m^3/s]: %.2e\n' % (tag, boundary_flows[tag])
     894       
     895        s += 'Total net boundary flow [m^3/s]: %.2e\n' % (total_boundary_inflow + total_boundary_outflow)
     896        s += 'Total volume in domain [m^3]: %.2e\n' % self.compute_total_volume()
    896897       
    897898        # The go through explicit forcing update and record the rate of change for stage and
     
    22782279                  the inflow region and then applied to update the
    22792280                  stage quantity.
     2281                 
     2282                  Note also that any reference to the internal attribute 'rate'
     2283                  later will use the one converted to m/s.
    22802284
    22812285    polygon: Specifies a polygon to restrict the rainfall.
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r7043 r7090  
    67456745       
    67466746
    6747         #------------------------------------------------------------------------------
     6747        #----------------------------------------------------------------------
    67486748        # Import necessary modules
    6749         #------------------------------------------------------------------------------
     6749        #----------------------------------------------------------------------
    67506750        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    67516751        from anuga.shallow_water import Domain
    67526752
    67536753
    6754         #------------------------------------------------------------------------------
     6754        #----------------------------------------------------------------------
    67556755        # Setup computational domain
    6756         #------------------------------------------------------------------------------
     6756        #----------------------------------------------------------------------
    67576757
    67586758        length = 100.
     
    67656765               
    67666766
    6767         #------------------------------------------------------------------------------
     6767        #----------------------------------------------------------------------
    67686768        # Simple flat bottom bathtub
    6769         #------------------------------------------------------------------------------
     6769        #----------------------------------------------------------------------
    67706770
    67716771        d = 1.0
     
    67766776
    67776777
    6778         #------------------------------------------------------------------------------
     6778        #----------------------------------------------------------------------
    67796779        # Slope
    6780         #------------------------------------------------------------------------------
     6780        #----------------------------------------------------------------------
    67816781               
    67826782        slope = 1.0/10  # RHS drops to -10m
     
    68126812       
    68136813
    6814         #------------------------------------------------------------------------------
     6814        #---------------------------------------------------------------------
    68156815        # Import necessary modules
    6816         #------------------------------------------------------------------------------
     6816        #---------------------------------------------------------------------
    68176817        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    68186818        from anuga.shallow_water import Domain
     
    68226822        from anuga.shallow_water.data_manager import get_flow_through_cross_section
    68236823
    6824         #------------------------------------------------------------------------------
     6824        #----------------------------------------------------------------------
    68256825        # Setup computational domain
    6826         #------------------------------------------------------------------------------
     6826        #----------------------------------------------------------------------
    68276827        finaltime = 300.0
    68286828
     
    68466846               
    68476847
    6848         #------------------------------------------------------------------------------
     6848        #----------------------------------------------------------------------
    68496849        # Setup initial conditions
    6850         #------------------------------------------------------------------------------
     6850        #----------------------------------------------------------------------
    68516851        slope = 0.0
    68526852        def topography(x, y):
     
    68606860                            expression='elevation')
    68616861
    6862         #------------------------------------------------------------------------------
     6862        #----------------------------------------------------------------------
    68636863        # Setup boundary conditions
    6864         #------------------------------------------------------------------------------
     6864        #----------------------------------------------------------------------
    68656865
    68666866        Br = Reflective_boundary(domain)      # Solid reflective wall
     
    68756875
    68766876       
    6877         #------------------------------------------------------------------------------
     6877        #----------------------------------------------------------------------
    68786878        # Evolve system through time
    6879         #------------------------------------------------------------------------------
     6879        #----------------------------------------------------------------------
    68806880        for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
    68816881       
     
    68876887           
    68886888       
    6889        
     6889           
     6890    def test_volume_conservation_inflow(self):
     6891        """test_volume_conservation
     6892       
     6893        Test that total volume in domain is as expected, based on questions
     6894        raised by Petar Milevski in May 2009.
     6895       
     6896        This test adds inflow at a known rate and verifies that the total
     6897        terminal volume is as expected.
     6898       
     6899        """
     6900
     6901        verbose = False
     6902       
     6903
     6904        #---------------------------------------------------------------------
     6905        # Import necessary modules
     6906        #---------------------------------------------------------------------
     6907        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     6908        from anuga.shallow_water import Domain
     6909        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     6910        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     6911        from anuga.shallow_water.shallow_water_domain import Inflow
     6912        from anuga.shallow_water.data_manager import get_flow_through_cross_section
     6913
     6914        #----------------------------------------------------------------------
     6915        # Setup computational domain
     6916        #----------------------------------------------------------------------
     6917        finaltime = 200.0
     6918
     6919        length = 300.
     6920        width  = 20.
     6921        dx = dy = 5       # Resolution: of grid on both axes
     6922       
     6923
     6924        points, vertices, boundary = rectangular_cross(int(length/dx),
     6925                                                       int(width/dy),
     6926                                                       len1=length, len2=width)
     6927
     6928
     6929        domain = Domain(points, vertices, boundary)   
     6930        domain.set_name('Inflow_volume_test')              # Output name
     6931               
     6932
     6933        #----------------------------------------------------------------------
     6934        # Setup initial conditions
     6935        #----------------------------------------------------------------------
     6936        slope = 0.0
     6937        def topography(x, y):
     6938            z=-x * slope
     6939            return z
     6940
     6941        domain.set_quantity('elevation', topography) # Use function for elevation
     6942        domain.set_quantity('friction', 0.0)         # Constant friction
     6943               
     6944        domain.set_quantity('stage',
     6945                            expression='elevation')  # Dry initially
     6946                           
     6947
     6948        #--------------------------------------------------------------
     6949        # Setup Inflow
     6950        #--------------------------------------------------------------
     6951
     6952        # Fixed Flowrate onto Area
     6953        fixed_inflow = Inflow(domain,
     6954                              center=(10.0, 10.0),
     6955                              radius=5.00,
     6956                              rate=10.00)                               
     6957                           
     6958        domain.forcing_terms.append(fixed_inflow)                           
     6959       
     6960        #----------------------------------------------------------------------
     6961        # Setup boundary conditions
     6962        #----------------------------------------------------------------------
     6963
     6964        Br = Reflective_boundary(domain) # Solid reflective wall
     6965        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     6966
     6967       
     6968        #----------------------------------------------------------------------
     6969        # Evolve system through time
     6970        #----------------------------------------------------------------------
     6971        ref_volume = 0.0
     6972        ys = 10.0  # Yieldstep
     6973        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
     6974       
     6975            # Check volume
     6976            assert num.allclose(domain.compute_total_volume(), ref_volume)
     6977       
     6978            if verbose :
     6979                print domain.timestepping_statistics()
     6980                print domain.volumetric_balance_statistics()
     6981                print 'reference volume', ref_volume
     6982           
     6983           
     6984            # Update reference volume
     6985            ref_volume += ys * fixed_inflow.rate
     6986       
     6987           
     6988           
     6989
     6990           
     6991    def test_volume_conservation_rain(self):
     6992        """test_volume_conservation
     6993       
     6994        Test that total volume in domain is as expected, based on questions
     6995        raised by Petar Milevski in May 2009.
     6996       
     6997        This test adds rain at a known rate and verifies that the total
     6998        terminal volume is as expected.
     6999       
     7000        """
     7001
     7002        verbose = False
     7003       
     7004
     7005        #---------------------------------------------------------------------
     7006        # Import necessary modules
     7007        #---------------------------------------------------------------------
     7008        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     7009        from anuga.shallow_water import Domain
     7010        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     7011        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     7012        from anuga.shallow_water.shallow_water_domain import Rainfall
     7013        from anuga.shallow_water.data_manager import get_flow_through_cross_section
     7014
     7015        #----------------------------------------------------------------------
     7016        # Setup computational domain
     7017        #----------------------------------------------------------------------
     7018        finaltime = 200.0
     7019
     7020        length = 300.
     7021        width  = 20.
     7022        dx = dy = 5       # Resolution: of grid on both axes
     7023       
     7024
     7025        points, vertices, boundary = rectangular_cross(int(length/dx),
     7026                                                       int(width/dy),
     7027                                                       len1=length, len2=width)
     7028
     7029
     7030        domain = Domain(points, vertices, boundary)   
     7031        domain.set_name('Rain_volume_test')              # Output name
     7032               
     7033
     7034        #----------------------------------------------------------------------
     7035        # Setup initial conditions
     7036        #----------------------------------------------------------------------
     7037        slope = 0.0
     7038        def topography(x, y):
     7039            z=-x * slope
     7040            return z
     7041
     7042        domain.set_quantity('elevation', topography) # Use function for elevation
     7043        domain.set_quantity('friction', 0.0)         # Constant friction
     7044               
     7045        domain.set_quantity('stage',
     7046                            expression='elevation')  # Dry initially
     7047                           
     7048
     7049        #--------------------------------------------------------------
     7050        # Setup rain
     7051        #--------------------------------------------------------------
     7052
     7053        # Fixed rain onto small circular area
     7054        fixed_rain = Rainfall(domain,
     7055                              center=(10.0, 10.0),
     7056                              radius=5.00,
     7057                              rate=10.00)   # 10 mm/s                           
     7058                           
     7059        domain.forcing_terms.append(fixed_rain)                           
     7060       
     7061        #----------------------------------------------------------------------
     7062        # Setup boundary conditions
     7063        #----------------------------------------------------------------------
     7064
     7065        Br = Reflective_boundary(domain) # Solid reflective wall
     7066        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     7067
     7068       
     7069        #----------------------------------------------------------------------
     7070        # Evolve system through time
     7071        #----------------------------------------------------------------------
     7072        ref_volume = 0.0
     7073        ys = 10.0  # Yieldstep
     7074        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
     7075       
     7076            # Check volume
     7077            V = domain.compute_total_volume()
     7078            msg = 'V = %e, Ref = %e' % (V, ref_volume)
     7079            assert num.allclose(V, ref_volume), msg
     7080       
     7081            if verbose :
     7082                print domain.timestepping_statistics()
     7083                print domain.volumetric_balance_statistics()
     7084                print 'reference volume', ref_volume
     7085                print V
     7086           
     7087           
     7088            # Update reference volume.
     7089            # FIXME: Note that rate has now been redefined
     7090            # as m/s internally. This is a little confusing
     7091            # when it was specfied as mm/s.
     7092           
     7093            delta_V = fixed_rain.rate*fixed_rain.exchange_area
     7094            ref_volume += ys * delta_V
     7095       
     7096           
     7097           
     7098           
     7099    def Xtest_rain_conservation_and_runoff(self):
     7100        """test_rain_conservation_and_runoff
     7101       
     7102        Test that total volume in domain is as expected, based on questions
     7103        raised by Petar Milevski in May 2009.
     7104       
     7105        This test adds rain at a known rate and verifies that the total
     7106        volume and outflows are as expected.
     7107       
     7108        """
     7109
     7110        # FIXME (Ole): Does not work yet. Investigate boundary flows
     7111       
     7112        verbose = True #False
     7113       
     7114
     7115        #---------------------------------------------------------------------
     7116        # Import necessary modules
     7117        #---------------------------------------------------------------------
     7118        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     7119        from anuga.shallow_water import Domain
     7120        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     7121        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     7122        from anuga.shallow_water.shallow_water_domain import Rainfall
     7123        from anuga.shallow_water.data_manager import get_flow_through_cross_section
     7124
     7125        #----------------------------------------------------------------------
     7126        # Setup computational domain
     7127        #----------------------------------------------------------------------
     7128        finaltime = 500.0
     7129
     7130        length = 300.
     7131        width  = 20.
     7132        dx = dy = 5       # Resolution: of grid on both axes
     7133       
     7134
     7135        points, vertices, boundary = rectangular_cross(int(length/dx),
     7136                                                       int(width/dy),
     7137                                                       len1=length, len2=width)
     7138
     7139
     7140        domain = Domain(points, vertices, boundary)   
     7141        domain.set_name('Rain_volume_runoff_test')         # Output name
     7142               
     7143
     7144        #----------------------------------------------------------------------
     7145        # Setup initial conditions
     7146        #----------------------------------------------------------------------
     7147        slope = 0.0
     7148        def topography(x, y):
     7149            z=-x * slope
     7150            return z
     7151
     7152        domain.set_quantity('elevation', topography) # Use function for elevation
     7153        domain.set_quantity('friction', 0.0)         # Constant friction
     7154               
     7155        domain.set_quantity('stage',
     7156                            expression='elevation')  # Dry initially
     7157                           
     7158
     7159        #--------------------------------------------------------------
     7160        # Setup rain
     7161        #--------------------------------------------------------------
     7162
     7163        # Fixed rain onto small circular area
     7164        fixed_rain = Rainfall(domain,
     7165                              center=(10.0, 10.0),
     7166                              radius=5.00,
     7167                              rate=10.00)   # 10 mm/s                           
     7168                           
     7169        domain.forcing_terms.append(fixed_rain)                           
     7170       
     7171        #----------------------------------------------------------------------
     7172        # Setup boundary conditions
     7173        #----------------------------------------------------------------------
     7174
     7175        Br = Reflective_boundary(domain) # Solid reflective wall
     7176        Bt = Transmissive_stage_zero_momentum_boundary(domain)
     7177        Bd = Dirichlet_boundary([-10, 0, 0])
     7178        domain.set_boundary({'left': Bt, 'right': Bd, 'top': Bt, 'bottom': Bt})
     7179
     7180       
     7181        #----------------------------------------------------------------------
     7182        # Evolve system through time
     7183        #----------------------------------------------------------------------
     7184        ref_volume = 0.0
     7185        ys = 10.0  # Yieldstep
     7186        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
     7187       
     7188            # Check volume
     7189            V = domain.compute_total_volume()
     7190            msg = 'V = %e, Ref = %e' % (V, ref_volume)
     7191            #assert num.allclose(V, ref_volume) or V < ref_volume, msg
     7192       
     7193            if verbose:
     7194                print domain.timestepping_statistics()
     7195                print domain.volumetric_balance_statistics()
     7196                print 'reference volume', ref_volume
     7197                print V
     7198           
     7199           
     7200            # Update reference volume.
     7201            # FIXME: Note that rate has now been redefined
     7202            # as m/s internally. This is a little confusing
     7203            # when it was specfied as mm/s.
     7204           
     7205            delta_V = fixed_rain.rate*fixed_rain.exchange_area
     7206            ref_volume += ys * delta_V
     7207       
     7208            # Compute outflow at right hand downstream boundary
     7209            boundary_flows, inflow , outflow = domain.compute_boundary_flows()
     7210            net_outflow = outflow - inflow
     7211       
     7212            outflow = boundary_flows['right']
     7213            if verbose:
     7214                print 'Outflow', outflow
     7215                print 'Net outflow', net_outflow
     7216       
     7217            # Update reference volume
     7218            ref_volume += ys * outflow           
     7219
     7220       
     7221               
    68907222       
    68917223    def test_inflow_using_flowline(self):
     
    69297261        dx = dy = 5          # Resolution: of grid on both axes
    69307262
    6931         points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
     7263        points, vertices, boundary = rectangular_cross(int(length/dx),
     7264                                                       int(width/dy),
    69327265                                                       len1=length, len2=width)
    69337266
     
    70137346                    #if verbose :
    70147347                    #    print domain.timestepping_statistics()
    7015                     #    print domain.volumetric_balance_statistics()                                   
    7016 
     7348                    #    print domain.volumetric_balance_statistics()
    70177349
    70187350                #--------------------------------------------------------------
     
    70217353                   
    70227354                # Square on flowline at 200m
    7023                 q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]])
    7024                 msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
     7355                q=domain.get_flow_through_cross_section([[200.0,0.0],
     7356                                                         [200.0,20.0]])
     7357                msg = 'Predicted flow was %f, should have been %f' %\
     7358                    (q, ref_flow)
    70257359                if verbose:
    7026                     print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
     7360                    print '90 degree flowline: ANUGA = %f, Ref = %f' %\
     7361                        (q, ref_flow)
    70277362                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
    70287363
    70297364                           
    70307365                # 45 degree flowline at 200m
    7031                 q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]])
    7032                 msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
     7366                q=domain.get_flow_through_cross_section([[200.0,0.0],
     7367                                                         [220.0,20.0]])
     7368                msg = 'Predicted flow was %f, should have been %f' %\
     7369                    (q, ref_flow)
    70337370                if verbose:
    7034                     print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
     7371                    print '45 degree flowline: ANUGA = %f, Ref = %f' %\
     7372                        (q, ref_flow)
    70357373                   
    70367374                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     
    71557493                    #if verbose :
    71567494                    #    print domain.timestepping_statistics()
    7157                     #    print domain.volumetric_balance_statistics()                                   
     7495                    #    print domain.volumetric_balance_statistics()
     7496                                                       
    71587497
    71597498
Note: See TracChangeset for help on using the changeset viewer.