source: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py @ 8254

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

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

def width(x):
    return numpy.ones_like(x)

def top(x):
    return 4.0 * numpy.ones_like(x)

def area(x):
    return height(x)*width(x)


#===============================================================================

def get_domain(dom):
    N = 100
    print "Evaluating domain with %d cells" %N
    domain = dom.Domain(*uniform_mesh(N, x_0 = -2.0*L_x, x_1 = 2.0*L_x), bulk_modulus = 75.0)

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

    domain.set_beta(1.0)

    domain.set_quantity('area', area)
    domain.set_quantity('stage', stage)
    domain.set_quantity('elevation',elevation)
    domain.set_quantity('width',width)
    domain.set_quantity('top',top)

    Br = dom.Reflective_boundary(domain)

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

    return domain

def animate_domain(domain, yieldstep, finaltime):
    import pylab
    pylab.ion()

    x = domain.get_centroids()
    z = domain.get_quantity('elevation', 'centroids')
    w = domain.get_quantity('stage', 'centroids')
    h = domain.get_quantity('height', 'centroids')
    v = domain.get_quantity('velocity', 'centroids')
    t = domain.get_quantity('top', 'centroids')

    plot1 = pylab.subplot(311)

    zplot, = pylab.plot(x, z)
    wplot, = pylab.plot(x, w)
    ztplot, = pylab.plot(x, z+t)

    plot1.set_ylim([-1,50])
    pylab.xlabel('Position')
    pylab.ylabel('Stage')

    plot2 = pylab.subplot(312)

    hplot, = pylab.plot(x, h)
    tplot, = pylab.plot(x, t)

    plot2.set_ylim([-1,10])
    pylab.xlabel('Position')
    pylab.ylabel('Height')

    plot3 = pylab.subplot(313)

    vplot, = pylab.plot(x, v)

    plot3.set_ylim([-15,15])
    pylab.xlabel('Position')
    pylab.ylabel('Velocity')

    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
        z = domain.get_quantity('elevation', 'centroids')
        w = domain.get_quantity('stage', 'centroids')
        h = domain.get_quantity('height', 'centroids')
        t = domain.get_quantity('top', 'centroids')
        v = domain.get_quantity('velocity', 'centroids')

        zplot.set_ydata(z)
        ztplot.set_ydata(z+t)
        wplot.set_ydata(w)
        hplot.set_ydata(h)
        tplot.set_ydata(t)
        vplot.set_ydata(v)

        pylab.draw()



def plot_domain(domain):
    import pylab

    x = domain.get_centroids()
    z = domain.get_quantity('elevation', 'centroids')
    w = domain.get_quantity('stage', 'centroids')
    h = domain.get_quantity('height', 'centroids')
    t = domain.get_quantity('top', 'centroids')
    u = domain.get_quantity('velocity', 'centroids')

    pylab.ioff()
    pylab.hold(False)

    plot1 = pylab.subplot(311)
    pylab.plot(x, z, x, w, x, z+t)
    plot1.set_ylim([-1,40])
    pylab.xlabel('Position')
    pylab.ylabel('Stage')

    plot2 = pylab.subplot(312)
    pylab.plot(x, h, x, t)
    plot2.set_ylim([-1,10])
    pylab.xlabel('Position')
    pylab.ylabel('Height')

    plot3 = pylab.subplot(313)
    pylab.plot(x, u)
    plot3.set_ylim([-15,15])
    pylab.xlabel('Position')
    pylab.ylabel('Velocity')

    pylab.show()

def write_domain(domain, outfile):
    x = domain.get_centroids()
    z = domain.get_quantity('elevation', 'centroids')
    w = domain.get_quantity('stage', 'centroids')

    f = open(outfile, 'w')
    for i in range(len(x)):
        f.write("%s %s %s\n" % (x[i], z[i], w[i]))
    f.close()