source: development/pyvolution-1d/pt.py @ 3295

from math import sqrt, pi
from shallow_water_1d import *
from Numeric import allclose, array, zeros, ones, Float, take
from config import g, epsilon

def test_sqrt():
    for i in range (80000):
        a = sqrt(4356)
        b = sqrt(2031)

def test_flux1():
    #Use data from previous version of pyvolution
    #normal = array([1.,0])
    #ql = array([-0.2, 2, 3])
    #qr = array([-0.2, 2, 3])
    ql = array([-0.2, 2])
    qr = array([-0.2, 2])
    zl = zr = -0.5
    #flux, max_speed = flux_function(normal, ql, qr, zl, zr)
    flux, max_speed = flux_function(1.0, ql, qr, zl, zr)
    print 'flux', flux

def test_compute_fluxes0():
    #Do a full triangle and check that fluxes cancel out for
    #the constant stage case
    
    print 'check min time step in compute fluxes is ok, John'
    
    #a = [0.0, 0.0]
    #b = [0.0, 2.0]
    #c = [2.0,0.0]
    #d = [0.0, 4.0]
    #e = [2.0, 2.0]
    #f = [4.0,0.0]
    a=0.0
    b=2.0
    c=4.0
    d=6.0
    e=8.0
    f=10.0
    
    points = [a, b, c, d, e, f]
    #bac, bce, ecf, dbe
    #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    
    domain = Domain(points)
    #domain.set_quantity('stage', [[2,2,2], [2,2,2],
    #                              [2,2,2], [2,2,2]])
    domain.set_quantity('stage', [[2,2], [2,2], [2,2], [2,2], [2,2]])
    domain.check_integrity()
    
    #assert allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
    assert allclose(domain.neighbours, [[-1,1], [0,2], [1,3],[2,4], [3,-1]])
    #assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
    #assert allclose(domain.neighbour_edges, [[-1,0], [1,0], [1,0], [1,0], [1,-1]])
    assert allclose(domain.neighbour_vertices, [[-1,0], [1,0], [1,0], [1,0], [1,-1]])
    
    zl=zr=0. #Assume flat bed
    
    #Flux across right edge of volume 1
    #normal = domain.get_normal(1,0)
    #ql = domain.get_conserved_quantities(vol_id=1, edge=0)
    #qr = domain.get_conserved_quantities(vol_id=2, edge=2)
    #flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
    
    #Flux across right edge of element 1
    ql = domain.get_conserved_quantities(vol_id=1, vertex=1)
    qr = domain.get_conserved_quantities(vol_id=2, vertex=0)
    #print 'qr and ql 1'
    #print qr
    #print ql
    flux0, max_speed = flux_function(1.0, ql, qr, zl, zr)
    
    #Check that flux seen from other triangles is inverse
    tmp = qr; qr=ql; ql=tmp
    #print 'qr and ql 2'
    #print qr
    #print ql
    #normal = domain.get_normal(2,2)
    #flux, max_speed = flux_function(normal, ql, qr, zl, zr)
    flux, max_speed = flux_function(-1.0, ql, qr, zl, zr)
    #print 'fluxes'
    #print flux0
    #print flux
    assert allclose(flux + flux0, 0.)
    
    #Flux across upper edge of volume 1
    #normal = domain.get_normal(1,1)
    #ql = domain.get_conserved_quantities(vol_id=1, edge=1)
    #qr = domain.get_conserved_quantities(vol_id=3, edge=0)
    #flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
    
    #Flux across left edge of element 1
    #ql = domain.get_conserved_quantities(vol_id=1, edge=0)
    #qr = domain.get_conserved_quantities(vol_id=0, edge=1)
    ql = domain.get_conserved_quantities(vol_id=1, vertex=0)
    qr = domain.get_conserved_quantities(vol_id=0, vertex=1)
    flux1, max_speed = flux_function(-1.0, ql, qr, zl, zr)
    
    #Check that flux seen from other triangles is inverse
    tmp = qr; qr=ql; ql=tmp
    #normal = domain.get_normal(3,0)
    #flux, max_speed = flux_function(normal, ql, qr, zl, zr)
    flux, max_speed = flux_function(1.0, ql, qr, zl, zr)
    assert allclose(flux + flux1, 0.)
    
    #Flux across lower left hypotenuse of volume 1
    #normal = domain.get_normal(1,2)
    #ql = domain.get_conserved_quantities(vol_id=1, edge=2)
    #qr = domain.get_conserved_quantities(vol_id=0, edge=1)
    #flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
    
    #Check that flux seen from other triangles is inverse
    #tmp = qr; qr=ql; ql=tmp
    #normal = domain.get_normal(0,1)
    #flux, max_speed = flux_function(normal, ql, qr, zl, zr)
    #assert allclose(flux + flux2, 0.)
    
    
    #Scale by edgelengths, add up anc check that total flux is zero
    #e0 = domain.edgelengths[1, 0]
    #e1 = domain.edgelengths[1, 1]
    #e2 = domain.edgelengths[1, 2]
    
    #assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.)
    print 'flux0',flux0
    print 'flux1',flux1
    assert allclose(flux0+flux1, 0.)
    
    #Now check that compute_flux yields zeros as well
    domain.compute_fluxes()
    
    #for name in ['stage', 'xmomentum', 'ymomentum']:
    for name in ['stage', 'xmomentum']:
        #print name, domain.quantities[name].explicit_update
        assert allclose(domain.quantities[name].explicit_update[1], 0)
        
def test_1d_solution_I():
    print "TEST 1D-SOLUTION I"
    
    L = 2000.0     # Length of channel (m)
    N = 100        # Number of compuational cells
    cell_len = L/N # Origin = 0.0
    
    points = zeros(N+1,Float)
    for i in range(N+1):
        points[i] = i*cell_len
        
    domain = Domain(points)
    domain.order = 2
    def stage(x):
        for i in range(len(x)):
            if x[i]<=1000.0:
                x[i] = 10.0
            else:
                x[i] = 5.0
        return x

    domain.set_quantity('stage', stage)
    #L = domain.quantities['stage'].vertex_values
    #print "Initial Stage"
    #print L

    domain.set_boundary({'exterior': Reflective_boundary(domain)})
    
    import time
    t0 = time.time()
    yieldstep = 1.0
    finaltime = 50.0
    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
        #xmom = domain.quantities['xmomentum']
        #xmom = xmom.centroid_values
        #stage = domain.quantities['stage']
        #stage = stage.centroid_values
        #print 'stage', stage
        #print 'xmom', xmom
        #for i in range N
        #u = xmom/stage
        pass
    
    print 'That took %.2f seconds' %(time.time()-t0)
    
    #L = domain.quantities['stage'].vertex_values
    #print "Final Stage"
    #print L
    
    C = domain.quantities['stage'].vertex_values
    #print C
    
    f = file('test_solution_Ix.out', 'w')
    for i in range(N):
        f.write(str(C[i,1]))
        f.write("\n")
    f.close