source: trunk/anuga_work/development/2010-projects/anuga_1d/sww-sudi/run-program-numeric-sudi.py @ 8427

from scipy import sin, cos, sqrt, linspace, pi, dot
from Numeric import zeros, Float, array
from gaussPivot import *
from analytical_prescription import *
from parameter import *
import os, time, csv, pprint
from domain_sudi import *
from config import g, epsilon
from rootsearch import *
from bisect_function import *

#Analytical computations#################################################################
def root_g(a,b,t):
    dx = 0.01
    def g(u):
        return u + 2.0*A*pi/T*sin(2.0*pi/T*(t+u))
    while 1:
        x1,x2 = rootsearch(g,a,b,dx)
        if x1 != None:
            a = x2
            root = bisect(g,x1,x2,1)
        else:
            break
    return root
def shore(t):
    a = -1.0
    b = 1.0
    #dx = 0.01
    u = root_g(a,b,t)
    xi = -0.5*u*u + A*cos(2.0*pi/T*(t+u))
    position = 1.0 + xi
    return position, u




#Numerical computations###################################################################
def newtonRaphson2(f,q,tol=1.0e-13): ##1.0e-9 may be too large
    for i in range(50):                                                                                                                                    
        h = 1.0e-4      ##1.0e-4 may be too large.
        n = len(q)
        jac = zeros((n,n),Float)
        if 1.0+q[0]-x<0.0:
            print "PROBLEM OCCURS.......... 1.0+q[0]-x=",1.0+q[0]-x
            q[0] = x-1.0 +0.001
        f0 = f(q)
        for i in range(n):
            temp = q[i]
            q[i] = temp + h
            f1 = f(q)
            q[i] = temp
            jac[:,i] = (f1 - f0)/h
        if sqrt(dot(f0,f0)/len(q)) < tol: return q
        dq = gaussPivot(jac,-f0)
        q = q + dq
        if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q
    print 'Too many iterations'

def elevation(X):
    N = len(X)
    z = zeros(N,Float)
    for i in range(N):
        z[i] = (h_0/L)*X[i] - h_0
    return z

def height(X):
    N = len(X)
    z = zeros(N,Float)
    for i in range(N):
        z[i] = h_0 - (h_0/L)*X[i]
    return z

def velocity(X):
    N = len(X)
    return zeros(N,Float)
    
boundary = { (0,0): 'left',(N-1,1): 'right'}

domain = Domain(points,boundary)
domain.order = 2
domain.set_timestepping_method('rk2')
domain.set_CFL(1.0)
domain.beta = 1.0
domain.set_limiter("minmod")

#def wu_at_O_by_formal_expansion(t):
#    W1 = j0(4*pi/T)*cos(2*pi*t/T)
#    U1 = -j1(4*pi/T)*sin(2*pi*t/T)
#
#    gu = -2*pi/T*j0(4*pi/T)*sin(2*pi*t/T)
#    gw = -2*pi/T*j1(4*pi/T)*cos(2*pi*t/T)
#    fu = gw
#    fw = gu + j1(4*pi/T)*sin(2*pi*t/T)
#    W2 = gu*U1 + gw*W1 - 0.5*U1**2
#    U2 = fu*U1 + fw*W1
#
#    w = W1*A + W2*A*A
#    u = U1*A + U2*A*A
#    return w,u

def wu_at_O_by_fixed_point_iteration(t):
    w1 = A*j0(4.0*pi/T)*cos(2.0*pi*t/T)
    u1 = -1.0*A*j1(4.0*pi/T)*sin(2.0*pi*t/T)

    #w2 = -0.5*u1**2.0 + A*j0(4.0*pi*(w1+1.0)**0.5/T)*cos(2.0*pi*(u1+t)/T)
    #u2 = -1.0*A*j1(4.0*pi*(w1+1.0)**0.5/T)*(w1+1.0)**-0.5*sin(2.0*pi*(u1+t)/T)

    #w3 = -0.5*u2**2.0 + A*j0(4.0*pi*(w2+1.0)**0.5/T)*cos(2.0*pi*(u2+t)/T)
    #u3 = -1.0*A*j1(4.0*pi*(w2+1.0)**0.5/T)*(w2+1.0)**-0.5*sin(2.0*pi*(u2+t)/T)

    #w4 = -0.5*u3**2.0 + A*j0(4.0*pi*(w3+1.0)**0.5/T)*cos(2.0*pi*(u3+t)/T)
    #u4 = -1.0*A*j1(4.0*pi*(w3+1.0)**0.5/T)*(w3+1.0)**-0.5*sin(2.0*pi*(u3+t)/T)

    #w5 = -0.5*u4**2.0 + A*j0(4.0*pi*(w4+1.0)**0.5/T)*cos(2.0*pi*(u4+t)/T)
    #u5 = -1.0*A*j1(4.0*pi*(w4+1.0)**0.5/T)*(w4+1.0)**-0.5*sin(2.0*pi*(u4+t)/T)    
    return w1,u1

def f_SUDI(t):
    timing = t*sqrt(g*h_0)/L
    w, u = wu_at_O_by_fixed_point_iteration(timing)
    wO = w*h_0
    uO = u*sqrt(g*h_0)
    zO = -h_0
    hO = wO - zO
    pO = uO * hO
    #[    'stage', 'xmomentum', 'elevation', 'height', 'velocity']
    return [wO,  pO,   zO,      hO,      uO]

def func(q): #Here q=(w, u)
    f = zeros(2,Float)
    f[0] = q[0] + 0.5*q[1]**2.0 - A*j0(4.0*pi/T*(1.0+q[0]-x)**0.5)*cos(2.0*pi/T*(tim+q[1]))
    f[1] = q[1] + A*j1(4.0*pi/T*(1.0+q[0]-x)**0.5)*sin(2.0*pi/T*(tim+q[1]))/(1+q[0]-x)**0.5
    return f


Ver = domain.vertices
n_V = len(Ver)
AnalitW_V = zeros((n_V,2), Float)
AnalitP_V = zeros((n_V,2), Float)
AnalitZ_V = zeros((n_V,2), Float)
AnalitH_V = zeros((n_V,2), Float)
AnalitU_V = zeros((n_V,2), Float)

T1 = Time_boundary(domain,f_SUDI)
D2 = Dirichlet_boundary([(h_0/L)*Ver.flat[len(Ver.flat)-1]-h_0,  0.0,  (h_0/L)*Ver.flat[len(Ver.flat)-1]-h_0,  0.0,  0.0])
domain.set_boundary({'left':T1,'right':D2})

domain.set_quantity('height',height)
domain.set_quantity('elevation',elevation)
domain.set_quantity('velocity',velocity)

Cen = domain.centroids
n_C = len(Cen)
AnalitW_C = zeros(n_C, Float)
AnalitP_C = zeros(n_C, Float)
AnalitZ_C = zeros(n_C, Float)
AnalitH_C = zeros(n_C, Float)
AnalitU_C = zeros(n_C, Float)

waktu = Tp/6.0
WAKTU = 14.0*Tp
yieldstep = finaltime = waktu
t0 = time.time()
counter=1

shorelines_numerical_sudi = zeros(int(WAKTU/waktu), Float)
time_instants = zeros(int(WAKTU/waktu), Float)


while finaltime < WAKTU+0.001:    
    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
        domain.write_time()
    time_instants[counter-1] = domain.time
    
    Stage = domain.quantities['stage']
    Momentum = domain.quantities['xmomentum']
    Elevation = domain.quantities['elevation']
    Height = domain.quantities['height']
    Velocity = domain.quantities['velocity']

    StageV = Stage.vertex_values
    MomV = Momentum.vertex_values
    ElevationV = Elevation.vertex_values
    HeightV = Height.vertex_values
    VelV = Velocity.vertex_values
    
    StageC = Stage.centroid_values
    MomC = Momentum.centroid_values
    ElevationC = Elevation.centroid_values
    HeightC = Height.centroid_values
    VelC = Velocity.centroid_values
    

    #######Numerical C-G at centroids########
    table = zeros((n_C,6),Float)
    for r in range(n_C):
        for c in range(6):
            if c==0:
                table[r][c] = Cen[r]
            elif c==1:
                table[r][c] = StageC[r]
            elif c==2:
                table[r][c] = MomC[r]
            elif c==3:
                table[r][c] = ElevationC[r]
            elif c==4:
                table[r][c] = HeightC[r]                
            else:
                table[r][c] = VelC[r]

    outname = "%s%04i%s%f%s" %("900_w1u1_Cen_numerical_", counter, "_", domain.time, ".csv")                
    outfile = open(outname, 'w')
    writer = csv.writer(outfile)
    for row in table:
        writer.writerow(row)        
    outfile.close()


    #######Numerical C-G at vertices########
    table = zeros((len(Ver.flat),6),Float)
    for r in range(len(Ver.flat)):
        for c in range(6):
            if c==0:
                table[r][c] = Ver.flat[r]
            elif c==1:
                table[r][c] = StageV.flat[r]
            elif c==2:
                table[r][c] = MomV.flat[r]
            elif c==3:
                table[r][c] = ElevationV.flat[r]
            elif c==4:
                table[r][c] = HeightV.flat[r]                
            else:
                table[r][c] = VelV.flat[r]

    outname = "%s%04i%s%f%s" %("900_w1u1_Ver_numerical_", counter, "_", domain.time, ".csv")                
    outfile = open(outname, 'w')
    writer = csv.writer(outfile)
    for row in table:
        writer.writerow(row)        
    outfile.close()
    
    
    #######Numerical shoreline########
    for s in range(2*n_V):
        heiL = HeightV.flat[s]
        momR = MomV.flat[s+1]
        if heiL >= 1e-6:
            if abs(momR)==0.0: #<1e-15:
                break
    shorelines_numerical_sudi[counter-1] = Ver.flat[s+1]


    #######Setting for the next loop########
    counter = counter+1
    finaltime = finaltime + waktu
    

########Saving the numerical shoreline#########
table = zeros((int(WAKTU/waktu), 2),Float)
for r in range(int(WAKTU/waktu)):
    for c in range(2):
        if c==0:
            table[r][c] = time_instants[r]               
        else:
            table[r][c] = shorelines_numerical_sudi[r]
outname = "%s" %("900_w1u1_shore_numerical.csv")                
outfile = open(outname, 'w')
writer = csv.writer(outfile)
for row in table:
    writer.writerow(row)        
outfile.close()