Ignore:
Timestamp:
Sep 3, 2008, 10:28:15 PM (14 years ago)
Author:
steve
Message:

Setting up some more unit tests for shallow_water_domain

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/anuga_1d/test_shallow_water.py

    r5587 r5727  
    66
    77from shallow_water_domain import *
     8from shallow_water_domain import flux_function as domain_flux_function
     9
    810from Numeric import allclose, array, ones, Float
    911
     
    3234        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
    3335
    34         stage_ud, xmom_ud = compute_fluxes_python(domain)
     36        stage_ud, xmom_ud = local_compute_fluxes(domain)
    3537       
    3638        domain.compute_fluxes()
    3739
    38         #print doamin.quantities['xmomentum'].explicit_update
    39         #print compute_fluxes_python(domain)
     40        #print domain.quantities['stage'].explicit_update
     41        #print domain.quantities['xmomentum'].explicit_update
     42        print stage_ud
    4043
    4144        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
     
    5154
    5255        #This assumes h0 = 1.0e-3!!
    53         edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
     56        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
    5457        #print maxspeed
    5558        #print edgeflux
    5659       
    57         assert allclose(array([1.998002, 8.89201198],Float), edgeflux)
    58         assert allclose(5.1284971665, maxspeed)
     60        assert allclose(array([2.0, 8.9],Float), edgeflux, rtol=1.0e-005)
     61        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
    5962
    6063        normal = -1.0
     
    6467        zr = 0.0
    6568
    66         edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
     69        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
    6770
    6871
     
    7073        #print edgeflux       
    7174       
    72         assert allclose(array([-1.998002, -8.89201198],Float), edgeflux)
    73         assert allclose(5.1284971665, maxspeed)
     75        assert allclose(array([-2.0, -8.9],Float), edgeflux, rtol=1.0e-005)
     76        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
    7477
    7578    def test_domain_flux_function(self):
     
    8083        zr = 0.0
    8184
    82         edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
     85        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
    8386
    8487        #print edgeflux
    8588
    86         from shallow_water_domain import flux_function as domain_flux_function
     89
    8790
    8891        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
     
    9396        assert allclose(domainmaxspeed, maxspeed)
    9497
     98    def test_gravity(self):
     99        """
     100        Compare shallow_water_domain gravity calculation
     101        """
     102               
     103        def slope_one(x):
     104            return x
     105                       
     106        domain = Domain(self.points)
     107        domain.set_quantity('stage',4.0)
     108        domain.set_quantity('elevation',slope_one)
     109        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
     110
     111        gravity(domain)
     112       
     113        #print domain.quantities['stage'].vertex_values
     114        #print domain.quantities['elevation'].vertex_values
     115        #print domain.quantities['xmomentum'].explicit_update       
     116
     117        assert allclose( array([-34.3, -24.5, -14.7], Float), domain.quantities['xmomentum'].explicit_update )
     118
     119
     120    def test_evolve(self):
     121        """
     122        Compare shallow_water_domain gravity calculation
     123        """
     124               
     125        def slope_one(x):
     126            return x
     127                       
     128        domain = Domain(self.points)
     129        domain.set_quantity('stage',4.0)
     130        domain.set_quantity('elevation',slope_one)
     131        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
     132
     133        yieldstep=0.01
     134        finaltime=0.01
     135
     136        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
     137            domain.write_time()
     138       
     139        print domain.quantities['stage'].vertex_values
     140        print domain.quantities['elevation'].vertex_values
     141        print domain.quantities['xmomentum'].vertex_values
     142
     143
     144        print domain.quantities['stage'].centroid_values
     145        print domain.quantities['elevation'].centroid_values
     146        print domain.quantities['xmomentum'].centroid_values   
     147
     148        assert allclose( array([-34.3, -24.5, -14.7], Float), domain.quantities['xmomentum'].explicit_update )
    95149
    96150
    97151#==============================================================================
    98152
    99 def compute_fluxes_python(domain):
     153def local_compute_fluxes(domain):
    100154    """Compute all fluxes and the timestep suitable for all volumes
    101155    in domain.
     
    122176    N = domain.number_of_elements
    123177
    124     tmp0 = zeros((N,),Float)
    125     tmp1 = zeros((N,),Float)
     178    tmp0 = zeros(N,Float)
     179    tmp1 = zeros(N,Float)
    126180
    127181    #Shortcuts
     
    222276
    223277
    224 def flux_function(normal, ql, qr, zl, zr):
     278def local_flux_function(normal, ql, qr, zl, zr):
    225279    """Compute fluxes between volumes for the shallow water wave equation
    226280    cast in terms of w = h+z using the 'central scheme' as described in
Note: See TracChangeset for help on using the changeset viewer.