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

import os
from math import sqrt, pi, sin, cos
import numpy
import time


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

#===============================================================================
# setup problem
#===============================================================================


z_infty = 10.0       ## max equilibrium water depth at lowest point.
L_x = 2500.0         ## width of channel
A0 = 0.5*L_x                  ## determines amplitudes of oscillations
omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation

def analytic_canal(x,t):
    u = numpy.zeros_like(x)    ## water velocity
    h = numpy.zeros_like(x)    ## water depth

    ## Define Basin Bathymetry
    z = numpy.zeros_like(x) ## elevation of basin
    w = numpy.zeros_like(x)   ## elevation of water surface

    z[:] = z_infty*(x**2/L_x**2)
    u[:] = -A0*omega*sin(omega*t)
    w[:] = numpy.maximum(z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x-0.5*A0/(L_x)*cos(omega*t)),z)
    h[:] = numpy.maximum(w-z, 0.0)

    T = 2.0*pi/omega

    return u,h,w,z, T


def stage(x):
    t=0.0
    u,h,w,z,T = analytic_canal(x,t)

    return w

def elevation(x):
    t=0.0
    u,h,w,z,T = analytic_canal(x,t)

    return z


def height(x):
    t=0.0
    u,h,w,z,T = analytic_canal(x,t)

    return h


#===============================================================================
finaltime = 500.0
yieldstep = 10.0



N = 100
print "Evaluating domain with %d cells" %N
domain = Domain(*uniform_mesh(N, x_0 = -2.0*L_x, x_1 = 2.0*L_x))

domain.set_spatial_order(2)
domain.set_timestepping_method('rk2')
domain.set_CFL(1.0)
#domain.set_limiter("minmod_kurganov")
domain.set_limiter("vanleer")

domain.set_beta(1.0)

domain.set_quantity('stage', stage)
domain.set_quantity('elevation',elevation)

domain.quantities['elevation'].extrapolate_second_order()
domain.quantities['stage'].extrapolate_second_order()

Br = Reflective_boundary(domain)

domain.set_boundary({'left': Br, 'right' : Br})

#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(311)

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

plot1.set_ylim([-1,40])

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


plot2 = pylab.subplot(312)

pylab.plot(X.flat,HeightV.flat)

plot2.set_ylim([-1,10])

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

plot3 = pylab.subplot(313)

pylab.plot(X.flat,VelocityV.flat)
plot3.set_ylim([-15,15])

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

pylab.show()