source: anuga_work/development/anuga_1d/anuga_1d_suggestion1_2_3/horizontal_dry_dam_h_1side.py @ 7576

import os
from math import sqrt, sin, cos, pi, exp
from shallow_water_domain_h_nonwellbalanced import *
from Numeric import zeros, Float

def analytic_dry_dam(C,t):
    #t  = 0.0     # time (s)
    g  = 9.81     # 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)
    z = zeros(n,Float)
    x = C-L/2.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)) 

        if ( 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+z , u*(h+z)


def stage_dam(x):
    h1 = 10.0
    h0 = 0.0
    y=zeros(len(x), Float)
    for i in range(len(x)):
        if x[i] <= 1000.0:
            y[i] = h1
        else:
            y[i] = h0
    return y



'''
def newLinePlot(title='Simple Plot'):
    import Gnuplot
    gg=Gnuplot.Gnuplot(persist=0)
    gg.title(title)
    gg('set data style linespoints')
    gg.xlabel('x')
    gg.ylabel('y')
    return gg

def linePlot(gg, x1, y1, x2, y2):
    import Gnuplot
    plot1=Gnuplot.PlotItems.Data(x1.flat, y1.flat, with="linespoints")
    plot2=Gnuplot.PlotItems.Data(x2.flat, y2.flat, with="lines 3")
    gg.plot(plot1, plot2)
'''


print "TEST 1D-SOLUTION I"

L=2000.0
N=400
cell_len=L/N
points=zeros(N+1, Float)
for i in range(N+1):
    points[i]=i*cell_len
domain=Domain(points)
domain.order = 2
domain.set_timestepping_method('rk2')
domain.cfl = 1.0
domain.limiter = "vanleer"
domain.set_quantity('stage',stage_dam)
domain.set_boundary({'exterior':Reflective_boundary(domain)})

X=domain.vertices
C=domain.centroids

#plot1x=newLinePlot("Stage")
#plot2x=newLinePlot("Momentum")

print 'THE DOMAIN ORDER is', domain.order
print 'The Time Stepping Method is', domain.set_timestepping_method
print 'THE DOMAIN LIMITER is', domain.limiter
import time
t0=time.time()
yieldstep=20.0
finaltime=20.0
print "integral", domain.quantities['stage'].get_integral()
for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
    domain.write_time()
    print "integral", domain.quantities['stage'].get_integral()
    if t>0.0:
        StageQ=domain.quantities['stage'].vertex_values
        MomentumQ=domain.quantities['xmomentum'].vertex_values
        ElevationQ=domain.quantities['elevation'].vertex_values
        
        #h, uh=analytical_sol(X.flat, domain.time)
        #linePlot(plot1x, X, StageQ, X, ElevationQ)
        #linePlot(plot2x, X, MomentumQ, X, MomentumQ)#, X, uh)   #(plot2x, X, MomentumQ, X, uh)
        #print "press return"
        #pass


        #This is only for dry dam
        N = float(N)
        StageC = domain.quantities['stage'].centroid_values
        XmomC = domain.quantities['xmomentum'].centroid_values
        ElevationC = domain.quantities['elevation'].centroid_values
        C = domain.centroids
        h, uh = analytic_dry_dam(C,domain.time)
        h_error = 1.0/(N)*sum(abs(h-StageC))
        uh_error = 1.0/(N)*sum(abs(uh-XmomC))
        print "h_error %.10f" %(h_error)
        print "uh_error %.10f"% (uh_error)
        print 'That took %.2f seconds' %(time.time()-t0)
        X = domain.vertices
        StageQ = domain.quantities['stage'].vertex_values
        XmomQ = domain.quantities['xmomentum'].vertex_values
        h, uh = analytic_dry_dam(X.flat,domain.time)
        #End for dry dam
    
        
        from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
        #ion()
        hold(False)
        clf()

        plot1 = subplot(211)
        plot(X,StageQ, X,h, X,ElevationQ)   #(X,ElevationQ,X,StageQ, X,h)
        plot1.set_ylim([-1.0,11.0]) #([9.999,10.006]) #([9.9999,10.0006]) #([-1.0,11.0])
        xlabel('Position')
        ylabel('Stage')
        legend( ('Numerical Solution', 'Analytical Solution', 'Bed Elevation'), 'upper right', shadow=False)      

        plot2 = subplot(212)
        plot(X,MomentumQ, X,uh)
        plot2.set_ylim([-5.0,30.0])
        legend( ('Numerical Solution', 'for momentum'), 'upper right', shadow=False)    
        xlabel('Position')
        ylabel('Xmomentum')
        
        #file = "dam_h_"
        #file += str(number_of_cells[i])
        #file += ".eps"
        #savefig(file)
        show()
        
print 'That took %.2f seconds'%(time.time()-t0)
raw_input("Press return key")