source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/profile_dry_dam.py @ 7884

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


def run_evolve():

    h1 = 10.0
    h0 = 0.01

    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 = 4.0
    yieldstep = 0.1
    L = 2000.0     # Length of channel (m)

    k = 0

    N = 800
    print "Evaluating domain with %d cells" %N
    domain = Domain(*uniform_mesh(N))

    domain.set_quantity('stage', stage)

    Br = Reflective_boundary(domain)

    domain.set_boundary({'left': Br, 'right' : Br})
    domain.order = 2
    domain.set_timestepping_method('rk2')
    domain.set_CFL(1.0)
    domain.set_limiter("minmod")
    #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)


import cProfile
cProfile.run('run_evolve()', 'dry_dam_prof')


import pstats
p = pstats.Stats('dry_dam_prof')

#p.strip_dirs().sort_stats(-1).print_stats(20)

p.sort_stats('cumulative').print_stats(30)