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

import os
from math import sqrt, pi
from shallow_water_domain_suggestion3 import *
from Numeric import allclose, array, zeros, ones, Float, take, sqrt
from config import g, epsilon

h_0 = 0.5
u_0 = 0.6
h_1 = 0.5
u_1 = 0.6


q = 0.3
g = 9.81

b = 4.0
xmax = 25.0
z_bmax = 0.2
Fr_0 = u_0/sqrt(g*h_0)

def stage(x):
    y,u,h = analytical_sol(x)
    return y
    

def elevation(x):
    b = 4.0
    xmax = 25.0
    z_bmax = 0.2
    z_b = zeros(len(x),Float)
    for i in range(len(x)):
        if (x[i] >= xmax-b) & (x[i] <= xmax+b):
            z_b[i] = z_bmax*(1.0-(x[i]-xmax)**2/b**2)
        else:
            z_b[i] = 0.0 
    return z_b
    #return 0.0

def xmomentum(x):
    return u_0*h_0


def bisection(func,xR,xL,H):
    while ((xR - xL) > epsilon):
  
        #Calculate midpoint of domain
        xM = xL + (xR - xL) / 2.0
  
        #Find f(xM)
        if ((f(xL,H) * f(xM,H)) > 0):
            #Throw away left half
            xL = xM
        else:
            #Throw away right half
            xR = xM
    return xR

def f(D,H):
    return D**3+D**2*(H-1.0-Fr_0**2/2.0)+Fr_0**2/2.0

def fprime(D,H):
    Fr_0 = u_0/sqrt(g*h_0)
    return 3*D**2+2*D*(H-1-Fr_0**2/2.0)

def analytical_sol(x):
    y = zeros(len(x),Float)
    u = zeros(len(x),Float)
    height = zeros(len(x),Float)
    for i in range(len(x)):
        if (x[i] >= xmax-b) & (x[i] <= xmax+b):
            zb = z_bmax*(1.0-(x[i]-xmax)**2/b**2)
        else:
            zb = 0.0 

        H = zb/h_0
        y1 = bisection(f,1.0,0.2/0.5,H)
        height[i] = y1*h_0
        y[i] = height[i]+zb
        u[i] = q/(height[i])
        

    return y,u,height

Fr_0 = u_0/sqrt(g*h_0)
zb=0.2
H = zb/h_0
D = bisection(f,1,0.2/0.5,H)

L = 50.0
N = 400
cell_len = L/N
points = zeros(N+1,Float)
for j in range(N+1):
    points[j] = j*cell_len
boundary = { (0,0): 'left',(N-1,1): 'right'}
domain = Domain(points,boundary)
D1 = Dirichlet_boundary([h_0,u_0*h_0,0.0,h_0,u_0])
D2 = Dirichlet_boundary([h_1,u_1*h_1,0.0,h_1,u_1])

domain.set_boundary({'left':D1,'right':D2})
domain.set_quantity('elevation',elevation)
domain.set_quantity('stage',stage)
domain.set_quantity('xmomentum',xmomentum)
X = domain.vertices
C = domain.centroids

Stage = domain.quantities['stage']
Xmom = domain.quantities['xmomentum']
Velocity = domain.quantities['velocity']

StageV = Stage.vertex_values
StageC = Stage.centroid_values
XmomV = Xmom.vertex_values
XmomC = Xmom.centroid_values
VelC = Velocity.centroid_values
VelV = Velocity.vertex_values
ElevationV = domain.quantities['elevation'].vertex_values

stage_bdry = Stage.boundary_values
xmom_bdry =  Xmom.boundary_values

w,u,h = analytical_sol(X.flat)
wc, uc, hc = analytical_sol(C)

h_error = zeros(1,Float)
uh_error = zeros(1,Float)
u_error = zeros(1,Float)
k = 0

domain.limiter = "minmod"
domain.order = 2
domain.set_timestepping_method('rk2')
domain.cfl = 1.0
print 'THE DOMAIN LIMITER is', domain.limiter
import time
yieldstep = finaltime = 0.2/71
t0=time.time()
while finaltime < 0.20001:
    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
        domain.write_time()
        print "integral", domain.quantities['stage'].get_integral()
        if t>= 0.0:      
            h_error[k] = cell_len*sum(abs(wc-StageC))
            uh_error[k] = cell_len*sum(abs(uc*hc-XmomC))
            u_error[k] = cell_len*sum(abs(uc-VelC))
            print "h_error %.10f" %(h_error[k])
            print "uh_error %.10f"% (uh_error[k])
            print "u_error %.10f"% (u_error[k])
            print "h_max %.10f"%max(abs(wc-StageC))
            print "uh max  %.10f"%max(abs(uc*hc-XmomC))
            print "u_max %.10f"%max(abs(uc-VelC))

            from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
            hold(False)
            clf()
        
            plot1 = subplot(311)
            plot(X,w,X,StageV,X,ElevationV)
            #xlabel('Position')
            ylabel('Stage')
            plot1.set_ylim([-0.1,1.2])
            legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'),
                   'upper right', shadow=False)
        
            plot2 = subplot(312)
            plot(X,u*h,X,XmomV)
            #xlabel('x (m)')
            ylabel('Momentum')
            #plot2.set_ylim([0.299,0.301])
        

            plot3 = subplot(313)
            plot(X,u, X,VelV)
            xlabel('Position')
            ylabel('Velocity')

            #filename = "steady_"
            #filename += str(t)
            #filename += ".png"
            #savefig(filename)
            #show()
    print 'That took %.2f seconds'%(time.time()-t0)
    # #raw_input("Press return key!")
    finaltime = finaltime + 0.2/71