source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py

import os
from math import sqrt, pi
import numpy
import time
#from Numeric import allclose, array, zeros, ones, Float, take, sqrt



from anuga_1d.sww.sww_domain import *
from anuga_1d.config import g, epsilon
from anuga_1d.base.generic_mesh import uniform_mesh

h1 = 10.0
h0 = 0.0

def analytical_sol(C,t):
    
    #t  = 0.0     # time (s)
    # gravity (m/s^2)
    #h1 = 10.0    # depth upstream (m)
    #h0 = 0.0     # depth downstream (m)
    L = 2000.0   # length of stream/domain (m)
    n = len(C)    # number of cells

    u = zeros(n,Float)
    h = zeros(n,Float)
    x = C-3*L/4.0
    

    for i in range(n):
        # Calculate Analytical Solution at time t > 0
        u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t) 
        h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
        u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1))
        h3_ = 1.0/(9.0*g)*((x[i]+L/2.0)/t+2*sqrt(g*h1))*((x[i]+L/2.0)/t+2*sqrt(g*h1))

        if ( x[i] <= -1*L/2.0+2*(-sqrt(g*h1)*t)):
            u[i] = 0.0
            h[i] = h0
        elif ( x[i] <= -1*L/2.0-(-sqrt(g*h1)*t)):
            u[i] = u3_
            h[i] = h3_

        elif ( x[i] <= -t*sqrt(g*h1) ):
            u[i] = 0.0 
            h[i] = h1 
        elif ( x[i] <= 2.0*t*sqrt(g*h1) ):
            u[i] = u3 
            h[i] = h3 
        else:
            u[i] = 0.0 
            h[i] = h0 

    return h , u*h, u



def stage(x):
    import numpy

    y = numpy.where( (x>=L/4.0) & (x<=3*L/4.0), h1 , h0)

#    for i in range(len(x)):
#        if x[i]<=L/4.0:
#            y[i] = h0
#        elif x[i]<=3*L/4.0:
#            y[i] = h1
#        else:
#            y[i] = h0
    return y



print "TEST 1D-SOLUTION III -- DRY BED"


finaltime = 20.0
yieldstep = 1.0
L = 2000.0     # Length of channel (m)

k = 0

N = 2000
print "Evaluating domain with %d cells" %N
domain = Domain(*uniform_mesh(N,x_1=L))
    
domain.set_quantity('stage', stage)

Br = Reflective_boundary(domain)

domain.set_boundary({'left': Br, 'right' : Br})
domain.set_spatial_order(2)
domain.set_timestepping_method('rk2')
domain.set_CFL(1.0)
domain.set_limiter("minmod_kurganov")
domain.set_beta(1.5)
#domain.h0=0.0001

t0 = time.time()

for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    domain.write_time()

print 'That took %.2f seconds' %(time.time()-t0)


N = float(N)
HeightC    = domain.quantities['height'].centroid_values
StageC     = domain.quantities['stage'].centroid_values
BedC       = domain.quantities['elevation'].centroid_values
C          = domain.centroids

HeightV    = domain.quantities['height'].vertex_values
StageV     = domain.quantities['stage'].vertex_values
BedV       = domain.quantities['elevation'].vertex_values
VelocityV  = domain.quantities['velocity'].vertex_values
X          = domain.vertices


import pylab

pylab.hold(False)

plot1 = pylab.subplot(211)

pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat)

plot1.set_ylim([-1,11])

pylab.xlabel('Position')
pylab.ylabel('Stage')
pylab.legend(('Analytical Solution', 'Numerical Solution'),
           'upper right', shadow=True)


plot2 = pylab.subplot(212)

pylab.plot(X.flat,VelocityV.flat)
plot2.set_ylim([-20,20])

pylab.xlabel('Position')
pylab.ylabel('Velocity')

pylab.show()