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

import os
from math import sqrt
from shallow_water_domain_h_wellbalanced import *
from Numeric import zeros, Float
from scipy import linspace
#from analytic_dam_sudi import AnalyticDam

h0=5.0
h1=10.0

#analytical_sol=AnalyticDam(h0,h1)


print "TEST 1D-SOLUTION I"

L=2000.0
N=16000

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 = "minmod"



def stage(x):
    y=zeros(len(x), Float)
    for i in range (len(x)):
        if x[i]<=L/4.0:
            y[i]=0.0 #h0
        elif x[i]<=3*L/4.0:
            y[i]=h1
        else:
            y[i]=h0
    return y

def elevation(x):
    y=zeros(len(x), Float)
    for i in range (len(x)):
        if x[i]<=500:
            y[i]=4*x[i]*x[i]/250000
        elif x[i]<=1500:
            y[i]=(500-x[i])/1000 + 4
        else:
            y[i]=-3*(x[i]-1000)*(x[i]-2000)/250000
    return y


P=linspace(0,2000,2000)
Elev_bed=elevation(P)

domain.set_quantity('stage',stage)
domain.set_quantity('elevation', elevation)
print "domain order", domain.order

domain.set_boundary({'exterior':Reflective_boundary(domain)})

X=domain.vertices
C=domain.centroids
#plot1x=newLinePlot("Height")
#plot2x=newLinePlot("Momentum")

import time
yieldstep=finaltime=0.0001
t0=time.time()

i=1
while finaltime < 0.00011:
    for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
        domain.write_time()
    #if t>0.0:
        #N = float(N)
        #StageC = domain.quantities['stage'].centroid_values
        #XmomC = domain.quantities['xmomentum'].centroid_values
        #VelC = domain.quantities['velocity'].centroid_values
        #C = domain.centroids
        #hC, uhC, uC = analytical_sol(C,domain.time)


        #h_error = sum(abs(hC-StageC))/N
        #uh_error = sum(abs(uhC-XmomC))/N
        #u_error = sum(abs(uC-VelC))/N
        #print "h_error %.10f" %(h_error)
        #print "uh_error %.10f"%(uh_error)
        #print "u_error %.10f" %(u_error)
        #print 'That took %.2f seconds' %(time.time()-t0)
        X = domain.vertices
        StageQ = domain.quantities['stage'].vertex_values
        XmomQ = domain.quantities['xmomentum'].vertex_values
        VelQ = domain.quantities['velocity'].vertex_values

        #h, uh, u = analytical_sol(X.flat, domain.time)
        #linePlot(plot1x, X, HeightQ, X, h)
        #linePlot(plot2x, X, MomentumQ, X, uh)
        #print "press return"
        #pass
   
        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
        hold(False)
           
        plot1 = subplot(311)
        plot(X,StageQ, P,Elev_bed) #plot(X,StageQ,'k-', P,Elev_bed,'k:')
        plot1.set_ylim([-1,11])
        plot1.set_xlim([0.0,2000.0])
        #xlabel('Position')
        ylabel('Stage')
        #legend(('Analytical Solution', 'Numerical Solution'),
        #   'lower right', shadow=False)
        
        plot2 = subplot(312)
        plot(X,XmomQ) #plot(X,XmomQ,'k-')
        #plot2.set_ylim([-5,35])
        plot2.set_xlim([0.0,2000.0])
        #legend(('Analytical Solution', 'Numerical Solution'),
        #        'lower right', shadow=False)
        #xlabel('Position')
        ylabel('Momentum')
        plot3 = subplot(313)
        plot(X,VelQ) #plot(X,VelQ,'k-')
        #plot2.set_ylim([-5,35])
        plot3.set_xlim([0.0,2000.0])
        legend(('Penyelesaian Numeris', ' '),
            'lower right', shadow=False)
        xlabel('Posisi')
        ylabel('Kecepatan')

        #file = "dam"
        #file += str(finaltime)
        #file += ".png"
        #savefig(file)
        #show()
        
        filename = "%s%02i%s" %("dam_init2_", i, ".eps")
        savefig(filename)
        finaltime = finaltime + 0.25
        i = i + 1