source: trunk/anuga_work/development/sudi/sw_1d/shock_detector/shv/dam_h_2sides_RUN.py @ 7910

import os
from math import sqrt
from sww_domain_shm import *
from Numeric import Float
from numpy import zeros
#from analytic_dam_sudi import AnalyticDam
#Note:apply analytical_sol given in debris avalanche solution
from parameters import *

N = int(N) # number of cells
print "number of cells=",N
#analytical_sol=AnalyticDam(h_0,h_1)
boundary = {(0,0):'left', (N-1,1): 'right'}
domain = Domain(points,boundary)
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]=bed_slope*x[i]#0.0
        elif x[i]<=3*L/4.0:
            y[i]=h_1
        else:
            y[i]=h_0
    return y
domain.set_quantity('stage',stage)

def elevation(x):
    y=zeros(len(x), Float)
    for i in range (len(x)):
        y[i] = bed_slope*x[i]
    return y
domain.set_quantity('elevation',elevation)

### ================ Define the boundary function =========================
#def f_right(t):
#    z_r = bed_slope*(0.5*L)
#    h_r = h_0 #+ bed_slope*cell_len
#    w_r = z_r + h_r
#    u_r = m*t   
#    #['stage', 'xmomentum', 'elevation', 'height', 'velocity']
#    return [w_r,  u_r*h_r,  z_r,  h_r,  u_r]

#T_right = Time_boundary(domain,f_right)
#D_right = Dirichlet_boundary([bed_slope*(0.5*L)+h_0,(m*domain.time)*h_0,bed_slope*(0.5*L),h_0,m*domain.time])
#D_left = Dirichlet_boundary([-1.0*bed_slope*(0.5*L),  0.0,  -1.0*bed_slope*(0.5*L),  0.0,  0.0])
D_right = Dirichlet_boundary([h_0,  0.0,  0.0,  h_0,  0.0])
D_left = Dirichlet_boundary([0.0,  0.0,  0.0,  0.0,  0.0])
domain.set_boundary({'left':D_left,'right':D_right})
### ================  End of the definition of boundary function ===========

X=domain.vertices
C=domain.centroids
import time
yieldstep=finaltime=1.0
t0=time.time()
i=1

while finaltime < 1.01:
    for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
        domain.write_time()
    finaltime = finaltime + 10.0 
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
    #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.flat
    StageQ = domain.quantities['stage'].vertex_values.flat
    XmomQ = domain.quantities['xmomentum'].vertex_values.flat
    VelQ = domain.quantities['velocity'].vertex_values.flat
    BedQ = domain.quantities['elevation'].vertex_values.flat
    #h, uh, u = analytical_sol(X.flat, domain.time)

    SD = domain.shock_detector

    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
    hold(False)
    plot1 = subplot(211)
    #plot(X,h, X,StageQ, X,BedQ) #plot(X,h,'k-',X,StageQ,'k--')
    plot(X,StageQ, X,BedQ)
    plot1.set_ylim([-1,11])
    plot1.set_xlim([0.0,2000.0])
    legend(('Numerical solution', 'Bed elevation'),
            'upper left', shadow=False)    
    #xlabel('Position')
    ylabel('Stage')
    

    plot2 = subplot(212)
    plot(points,SD)
    plot2.set_xlim([0.0,2000.0])
    xlabel('Position')
    ylabel('Smoothness indicator')

    show()
    
    #filename = "%s%04i%s" %("dam_", i, ".eps")
    #savefig(filename)
    #finaltime = finaltime + 0.25
    #print "finaltime=", finaltime
    #i = i + 1
    #print "The domain.limiter is", domain.limiter
    #print 'That took %.2f seconds'%(time.time()-t0)
    #print '============================================================================='